# Elasto-viscoplastic model for rayon yarns

Marco Moscatelli* · Lucas Pires da Costa · Paola Caracino · Simone Agresti · Giorgio Novati · Claudia Comi

Received: 5 October 2023 / Accepted: 11 March 2024 / Published online: 21 May 2024
© The Author(s) 2024

**Abstract** In this paper we develop a new model for the simulation of the mechanical behavior of rayon twisted yarns, at macroscopic level. A yarn with its continuous filaments is represented by an equivalent three-dimensional solid of cylindrical shape, discretized by finite elements, with properly defined local anisotropic material properties. The new constitutive model, inspired by experimental results on rayon untwisted yarns, is formulated in the framework of the thermodynamics of irreversible processes and includes visco-elastic and visco-plastic dissipation mechanisms. The effect of twist is taken into account by including the direction of the fibers in the free

energy definition. The overall model is validated comparing numerical and experimental results on twisted rayon yarns.

**Keywords** Yarns · Fibrous materials · Anisotropy · Fiber orientation · Generalized standard materials · Constitutive model · Visco-plasticity · Visco-elasticity · Rayon · Textile materials

## 1 Introduction

Polymeric yarns made of hundreds of filaments (or fibers) twisted together exhibit very high flexural flexibility while having high stiffness and strength in the longitudinal direction. Yarns (of the same or different filament materials) are often twisted together to form multi-ply yarns, also known as cords.

Polymer-based yarns and cords are extensively used as reinforcement in composite materials and, more specifically, in rubber composites. The latter materials are employed in various products, with tires, belting systems and hoses being the most notable examples [1]. The textile constituent combined with the elastomer matrix results in a structure that exhibits high stiffness in the direction of the reinforcement and high flexibility in the perpendicular plane [2].

Textile reinforcements (unidirectional or in the form of fabric) are very versatile and are employed in diverse contexts, for instance in concrete structural

M. Moscatelli (✉) · L. Pires da Costa · G. Novati ·

C. Comi

Department of Civil and Environmental Engineering,
Politecnico di Milano, Piazza Leonardo da Vinci 32,
20133 Milan, Italy

e-mail: marco.moscatelli@polimi.it

L. Pires da Costa

e-mail: lucas.piresdacosta@polimi.it

G. Novati

e-mail: giorgio.novati@polimi.it

C. Comi

e-mail: claudia.comi@polimi.it

P. Caracino · S. Agresti

R&D Materials, Reinforcements group, Pirelli Tyre Spa,
Via Piero e Alberto Pirelli 25, 20126 Milan, Italy

e-mail: paola.caracino@pirelli.com

S. Agresti

e-mail: simone.agresti@pirelli.com
elements, in architectural membranes, and as pre-
forms for composites. They have also been used for
artificial implants: in [3] a knit of yarns is combined
with a silicon-based matrix to create an aortic valve
with sufficient mechanical stiffness.

For rubber composites, different reinforcement
materials have been employed, the most impor-
tant ones being cotton, rayon, polyamides, polyester
(PET), aramids and polyethylene naphthalate (PEN).
These materials are selected depending on the appli-
cation envisaged. Rayon is a good choice for passen-
ger car tires, because of its low heat shrinkage, high
modulus and good adhesion properties (see e.g. [4]
for a review on the main properties of this material).
The reinforcements play a fundamental role to sup-
port the loads and contain the air pressure, providing
the axial and lateral rigidity for acceleration, braking
and cornering [5]. The accuracy in predicting the per-
formance of a tire is thus strongly dependent on the
validity of the material models chosen to describe the
mechanical behavior of the reinforcement component.

Experimental tests carried out in [6] on bundles of
straight rayon fibers (untwisted rayon yarns), under
monotonic and cyclic loading, were aimed at charac-
terizing the material nonlinearity of rayon and gave
evidence of an elastic-visco-plastic behavior.

For twisted yarns and cords, the knowledge of
the actual filament trajectories in a textile specimen
subjected to a tensile test represents an important
piece of information, as the orientation of the fib-
ers strongly influences the mechanical response of
the reinforcement and hence should be represented
with good accuracy in a numerical model to be used
for simulations. To gain such knowledge, different
experimental techniques can be employed. In [7], the
so-called “fiber tracer technique” was first developed
to determine the path of some fibers. This method
is not able to provide a full geometric description of
the reinforcement, as only a few fibers can be traced.
More recently, an X-ray microtomography and image
processing technique has been developed in [8] to
accurately reconstruct the geometry of single- and
multi-ply yarns.

Based on the knowledge of the material behav-
ior and the orientation of the fibers, the study and
modeling of polymeric yarns and cords have been
addressed at different scales in the literature.

Micromechanical approaches try to model the tex-
tile reinforcements at the level of each elementary

component, which is considered to be one filament
or a bundle of a few filaments together. Typically,
a finite element (FE) method is developed where
each elementary component is modeled either with
rod elements, as in [9], or beam elements, as in [10,
11]. Contact and nonlinear material behavior can
be accounted for with this approach. These models
are useful to gain information on the local interac-
tions between the different filaments, for instance to
determine the initial configuration of the fibers [12].
Nevertheless, these methods require a considerable
modeling effort. They are thus limited to small-size
problems, to check the validity of the assumptions
and hypotheses used for the upper scales.

At a macroscopic scale, the analyses of fiber-
reinforced rubber structures rely mainly on numeri-
cal methods. Structural reinforcements are gener-
ally modeled using two approaches. In the first one
the entire structure is considered as a continuum
and the modeling is based on the composite materi-
als theory, see e.g. [13]. In the second one the rub-
ber matrix and the yarns/cords are considered sepa-
rately and the reinforcements are modeled using the
so-called “rebar” elements. The basic idea is to insert
special bar elements with their own material proper-
ties into volume elements that should be strengthened
[14–16]. This latter technique is especially used for
tires.

To feed the macroscopic models with a feasi-
ble description of the textile (rebar) components,
an intermediate scale is necessary where the global
behavior of yarns and cords is considered. The
early approaches relied on an energy-based theory,
first developed in [17] for single-ply yarns and later
improved to account for more complex fibers trajecto-
ries, see e.g. [18]. These methods generally allow for
the prediction of the monotonic uniaxial stress–strain
response only. More efficient models based on the
continuum theory have thus been developed. In [19,
20], a one-dimensional thermodynamically consistent
viscoelastoplastic phenomenological model is devel-
oped for flax twisted yarns embedded in an epoxy
matrix. This model does not account explicitly for the
geometric effects caused by the torsion of the fibers
composing the yarns. The same composites are also
analyzed in [21], where matrix and fibers are modeled
separately and a thermodynamically consistent ani-
stropic elastoplastic constitutive model for the flax
fibers is considered; the reinforcement in this case
<figure><img src="image_1.png" /><figcaption>**Fig. 1** Microscopic image of a rayon yarn (480 twists per meter) with the definition of the fiber direction **a** and yarn equivalent representation through a homogeneous cylinder:</figcaption></figure>

the blue and red lines are helices representing an external and internal fiber, respectively

is composed of unidirectional non-twisted fibers. In [22], a hyperelastic anisotropic constitutive model is developed for multi-ply yarns. The inelastic effects are nevertheless disregarded and the phenomenological approach is built starting directly from the twisted yarns, i.e. without a fitting of the experimental curves valid for the fibers.

Yarns can also be arranged in woven fabrics, with further geometric complexity. Again with reference to the elastic behavior only, in [23], the effective properties of a fabric are obtained by numerical analyses of a unit cell at varying weaving patterns.

There is a need for a global three-dimensional model that, starting from the material behavior at the level of the fibers and including material inelastic effects together with geometric effects, can effectively predict the behavior of rayon yarns and cords in any stress condition. In the present work, we aim to contribute in this direction by developing a new approach integrating material inelastic effects and geometrical ones. The result is a novel phenomenological material model for twisted yarns and cords made up of polymeric fibers. We treat the textile reinforcements as a continuum, characterized by a viscoelastic and visco-plastic transversely-isotropic behavior, with material directions locally defined and depending on the filaments twist level. More precisely, we approximate the complex geometry of the twisted single-ply yarns as a bundle of filaments arranged as concentric helices, directed along the unit vector **a**(*x*, *t*), and we replace the real geometry of the yarns by an equivalent three-dimensional continuum of cylindrical shape, made of a material modeled by the non-isotropic constitutive law here developed. As cords are made up of single-ply yarns, the same material model could also be applied to them, with a proper definition of the material directions. Figure 1 schematically shows the definition of material directions on a rayon yarn with a

linear density of 1840 dtex and 480 twists per meter and the equivalent representation through a continuum cylinder with anisotropic properties. The blue and red lines are helices representing an external fiber (at the maximum radius) and an internal fiber (at an intermediate radius), respectively.

The paper is organized as follows. In Sect. 2, we experimentally characterize the rayon fibers and describe the main features of their mechanical behavior. In Sect. 3, we present the geometrical modeling of the twisted yarns, to explicitly take into account the fibers’ directions. In Sect. 4, we develop a novel anisotropic elasto-viscoplastic model in the thermodynamic framework. Section 5 traces a strategy for the identification of the material parameters. In Sect. 6 we describe the model implementation in a finite element (FE) analysis procedure and compare the experimental results on untwisted and twisted yarns with those obtained from three-dimensional (3D) FE simulations. Some conclusions and future prospects are finally given in Sect. 7.

*Notation:* Throughout this work, normal-face letters (a) stand for scalars, while boldface letters (a) can denote both vectors and tensors. In general, we drop everywhere the dependence from the space and time variables (*x*, *t*).

## 2 Experimental evidences

Polymeric twisted yarns are the basic constituents of the cords which are often used as reinforcing elements in tires, as a lighter alternative to steel cords. An extended experimental campaign has been recently carried out at the Pirelli laboratory (in Milan) on rayon yarns to obtain essential information for the development of a proper constitutive law at the macroscopic level.
Rayon yarns consisting of approximately 1000 filaments, with a linear density of 1840 dtex (1840 g per 10,000 m) are considered in this study. To distinguish between the material behavior of the filaments and the effects induced by twist, experimental tests were first carried out on untwisted yarns, here denoted as yarn Y0; successively, tests on twisted yarns were performed with two different twist levels (with Z-twist direction), 380 and 480 twists per meter, hereafter called yarn Y38 and Y48, respectively. In this section, we only comment on the results obtained on yarn Y0, which are the basis for the development of the constitutive model. Tests on yarns Y38 and Y48 will be presented in Sect. 6.2 together with the model predictions.

The aim of the experimental work is to obtain a
complete characterization for rayon, tested under dif-
ferent conditions.

Since rayon fibers are strongly affected by the level
of moisture, experimental tests were carried out both
in “Dry conditions” (denoted by “D”), with the sam-
ple kept into an oven at 100°C for at least two hours
before the test, and in “Humid conditions” (denoted
by “H”), with the sample conditioned for 24 h in
the same climatic room where the experiments take
place, at 20°C and 65% of relative humidity. These
tests were carried out according to the procedures
prescribed by the BISFA guidelines [24].

Uniaxial tensile tests, under monotonically
increasing displacement, were thus performed on
greige yarns specimens, both in dry and humid condi-
tions. Bollard grips were used in the tests to prevent
specimen breakage in the clamping area. The adopted
gauge length is 500 mm. The strain rate of the tests
in dry conditions is fixed at 0.017 s⁻¹ to maintain
approximately the same Moisture Content (MC) dur-
ing each test. This can be checked by considering
Fig. 2, where the humidity gain in time is reported
starting from dry conditions during the first 130 min.
Even though the humidity gain is rather fast in the
first phase, for a strain rate of 0.017 s⁻¹ and a maxi-
mum strain of about 9%, the total duration of the test
is approximately of 5.3 s; and the variation of moisture content can be considered negligible in such time interval.

Adopting a strain rate of 0.017 s⁻¹, the influence
of humidity on the material behavior is shown in
Fig. 3a, where the experimental response for yarn Y0
in dry conditions (in red) is compared to that obtained

in humid conditions (in blue). The shaded regions
contain the stress–strain curves obtained from test-
ing 20 nominally identical specimens and the dashed
curves represent the averages. Although the actual
experimental results are obtained in terms of meas-
ured force, axial stress appears on the ordinates of
Fig. 3a, b. Such stress has been computed assuming a
uniform distribution and a cross-section area equal to
0.17 mm² based on the microscopy images in Fig. 5.
The mean values and the standard deviation of the
experimental force at nominal strains of 1%, 3% 5%,
and 7% are reported in Table 1.

The material was then tested in humid conditions at various strain rates, namely of 0.017 s⁻¹, 0.01 s⁻¹, 0.003 s⁻¹, 0.0017 s⁻¹ and 0.0003 s⁻¹. Figure 3b shows the results obtained in terms of measured stress versus nominal strain. The material exhibits a rate-dependent behavior. From Fig. 3a, b, it appears that the material response is linear for small values of applied load, but becomes nonlinear for higher values.

Further insight into the material response can
be obtained from non-monotonic tensile tests. The
adopted strain rate value is ±0.003 s⁻¹ for increas-
ing and decreasing elongation. Five specimens were
tested. Figure 4 shows the imposed strain history and
the corresponding stress–strain response for one of
the 5 specimens.

Upon unloading permanent strains are recorded
and hysteresis loops during unloading and reloading
are visible.

From the described experimental results, we can
interpret the material response as elasto-visco-plastic
at the scale of the filaments constituting the yarn.

<figure><img src="image_2.png" /><figcaption>**Fig. 2** Percentage of moisture content (MC) versus time for a specimen of yarn Y0. The measurement starts from dry conditions and stops when saturation is almost reached</figcaption></figure>
<table><caption>**Table 1** Mean values and standard deviation of the force at different strain levels for untwisted yarns at 0.017 s⁻¹ strain rate</caption><thead><tr><th rowspan="2">Strain</th><th colspan="4">Humid yarn Y0</th><th colspan="4">Dry yarn Y0</th></tr><tr><th>1%</th><th>3%</th><th>5%</th><th>7%</th><th>1%</th><th>3%</th><th>5%</th><th>7%</th></tr></thead><tbody><tr><td>Mean force [N]</td><td>19.25</td><td>29.83</td><td>39.34</td><td>49.14</td><td>26.22</td><td>52.29</td><td>63.29</td><td>73.47</td></tr><tr><td>Std Dev [N]</td><td>0.64</td><td>0.87</td><td>0.83</td><td>3.11</td><td>0.62</td><td>1.42</td><td>1.47</td><td>1.59</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Fig. 4** Cyclic test on untwisted yarns at 0.003s⁻¹: **a** strain history; **b** stress versus strain response</figcaption></figure>

<figure><img src="image_4.png" /></figure>

### 3 Geometrical model of twisted single-ply yarns

The linear and nonlinear behavior of twisted yarns is strongly dependent on the local fiber orientation, in particular on the angle each fiber forms with the yarn axis.

In defining the geometry of a single yarn, the model usually adopted is that of coaxial-helixes (see e.g. [25]). The assumed idealized helical (Z-twist) yarn geometry is illustrated in Fig. 5a. Each fiber follows a helical path around one of the concentric cylinders forming the yarn. All the helical fibers (of equal or different radii) have the same pitch *h* (defined as the length along the yarn axis of one turn of twist). The following notation is adopted: *R* = yarn radius; *ρ* = radius of internal fiber; *a* = unit vector tangent

to filament trajectory; {x₁, x₂, x₃} = Cartesian coordinates, with x₃ taken along the yarn axis; α (or θ) = angle between **a**(*R*) (or **a**(*ρ*)) and the yarn axis; *N* = yarn twist = turns per unit length (with *N* = 1/*h*).

The aforementioned geometric model implies:

$$ \tan \theta = \frac{\rho}{R} \tan \alpha, \quad \tan \alpha = 2\pi RN. \qquad (1) $$

The unit tangent vector at each point **a**(*x*) hence reads

$$ \mathbf{a} = \frac{1}{\sqrt{\rho^2 + \left(\frac{1}{2\pi N}\right)^2}} \begin{cases} -x_2 \\ x_1 \\ 1/2\pi N \end{cases}. \qquad (2) $$
<figure><img src="image_5.png" /><figcaption>**Fig. 5** a Definition of yarn axis, trajectory of an external and internal fiber for a Z-twist, and geometric parameters; b cross section of rayon yarn embedded in rubber: definition of equivalent radius *R*</figcaption></figure>

As recently discussed in [8], the determination of the real trajectories would require experimental measurements, e.g. by ad-hoc micro tomography analyses. The experimental measures performed in that work on nylon yarns show that the inclination angle corresponding to the maximum density is a bit lower than the one predicted by the helix model (1) and that some fibers are more inclined than $\alpha$ in (1). However, the patterns of bivariate histograms of the orientation fiber combined with the radial position, indicate a preferential orientation which follows relation (1). For the above reason and since micro tomography data are not available for the rayon yarns here considered, the expression (2) will be used.

Furthermore, upon tensile loads, the fiber directions change and one should consider in general time dependent directions $a = a(x, t)$. In the present work, for the sake of simplicity, we will neglect this effect and we will refer to the initial directions, only related to the initial twist $N$ and to the initial radius $R$. This latter is determined as the radius of a circle of area equal to that obtained from the microscopic measurement of the yarn cross-section area shown in Fig. 5b. This information is fundamental for the determination of the stress from the force history obtained from the strain-driven experimental tests. We use a radius $R = 0.233$ mm (area $A = 0.17$ mm²) for all the yarns considered.

## 4 Anisotropic elasto-viscoplastic model

From the experimental tests described in Sect. 2 one can draw the following indication for a proper constitutive model: (i) there is an initial linear behavior that can be described by linear elasticity; (ii) the subsequent phase can be interpreted as elastoplastic, characterized by permanent strains; (iii) the yield limit initially grows nonlinearly and then almost linearly with the strain, suggesting nonlinear hardening, asymptotically tending to a linear one; (iv) the hardening response is rate-dependent, showing a visco-plastic behavior; (v) there are viscous effects also in the elastic regime as shown by the presence of hysteresis loops in the unloading-reloading cycles. To interpret all these phenomenological aspects, in [26] we proposed the one-dimensional six-parameter rheological model depicted in Fig. 6.

The model consists of a linear spring (E) in series with a Kelvin unit ($\tilde{E}, \tilde{\eta}$) and with a Bingham viscoplastic unit comprising three elements in parallel: a frictional device ($\sigma_y$), a possibly nonlinear hardening spring (H) and a viscous dash-pot ($\eta$). Based on this rheological model, the total strain is hence the sum of three contributions: elastic, viscoelastic and viscoplastic:

$$ \varepsilon = \varepsilon^e + \varepsilon^{ve} + \varepsilon^{vp} \quad (3) $$

The generalization of the 1D model to a three-dimensional continuum context is dealt with in Sect. 4.2,
<figure><img src="image_6.png" /><figcaption>Fig. 6 One-dimensional visco-elastic-visco-plastic rheological model (six-parameter model)</figcaption></figure>

while Sect. 4.1 summarizes the main assumptions and limitation of the model.

## 4.1 Basic assumptions

In order to obtain a 3D material model that can describe the yarn behavior in structural analyses of rubber composites, we make several simplifying assumptions that define the applicability domain. Some restrictions could be removed without particular difficulties, at the price of additional complexity.

1. The filament or fiber direction **a** is known from the yarn construction at each material point;

2. The yarn deforms in small strain regime; thus the fiber directions do not vary with time;

3. Isothermal conditions;

4. The activation of the plastic behavior only depends on the stress component acting on the fiber direction;

5. The fibers have a very low stiffness when compressed and to account for the differences of behavior in tension and compression, one should adopt piecewise elasticity in the formulation with additional material parameters (see e.g. [27]). This is not done here, as in this work we only consider numerical simulations of tensile tests.

## 4.2 Thermodynamic formulation

Let us consider a material composed of a single family of fibers, directed along the unit vector **a**(**x**) in the reference configuration, undergoing small strains, under isothermal conditions. Unlike in the case of a composite, in which the fibers are embedded in a matrix, here only the fibers are present, twisted together to form a single yarn and we use

the formulation to define an equivalent continuous anisotropic material. The model is developed in the framework of the thermodynamics of irreversible processes with internal variables, and one can refer to [28] for a general treatment.

As in the one dimensional case (see Fig. 6), the total strain **ε** is expressed through strain additivity as the sum of the elastic **ε**<sup>e</sup>, viscoelastic **ε**<sup>ve</sup> and viscoplastic **ε**<sup>vp</sup> strains:

$$ \varepsilon = \varepsilon^e + \varepsilon^{ve} + \varepsilon^{vp}. \qquad (4) $$

As the response of the material depends on the fibers orientation, we assume that the Helmholtz free-energy (per unit of mass) **ψ** depends on the unit vector **a**. We postulate a decoupled representation of the Helmholtz free-energy density **ψ** into “elastic” and “inelastic” contributions, respectively denoted as **ψ**<sub>e</sub> and **ψ**<sub>i</sub>. Adopting the classic formalism of thermodynamics with internal variables, **ψ** is thus a function of state and can be written as follows:

$$ \psi = \psi(\varepsilon, \varepsilon^{ve}, \varepsilon^{vp}, V_k; \mathbf{a}) = \psi_e(\varepsilon^e; \mathbf{a}) + \psi_i(\varepsilon^{ve}; V_k; \mathbf{a}), \qquad (5) $$

where $V_k$ denotes a set of internal (scalar, vector or tensor) variables, which in the present model describe the material hardening. Specifically, in our case, we choose to define two internal variables $\alpha_l$ and $\alpha_{nl}$, associated respectively to a linear and a nonlinear kinematic hardening. The observable variables are the total strain **ε** and the fiber direction **a**, while $\varepsilon^{ve}$ and $\varepsilon^{vp}$ play the role of internal variables. Note that, following the small strain hypothesis, **a** remains fixed in time and thus is a space dependent material parameter for our problem.

In relation (5), the elastic contribution **ψ**<sub>e</sub> is related to the spring **E** in the 1D rheological model, while the inelastic contribution **ψ**<sub>i</sub> comes from the deformation of the viscoelastic Kelvin element, with the spring $\tilde{E}$ in parallel to the dash-pot $\tilde{\eta}$, and of the viscoplastic Bingham element, composed of a slider with yield limit $\sigma_y$ in parallel with a spring **H** and a dash-pot $\eta$. The inelastic term is thus further split into two contributions:

$$ \psi_i(\varepsilon^{ve}, V_k; \mathbf{a}) = \psi_{ve}(\varepsilon^{ve}; \mathbf{a}) + \psi_{vp}(\alpha_l, \alpha_{nl}), \qquad (6) $$

where $\psi_{ve}$ is the viscoelastic contribution, while the term $\psi_{vp}$ represents the energy locked in the material due to microscopic rearrangements related to hardening. The current state of $\psi_{vp}$ is assumed to be independent of the material orientation **a**.
As detailed in [29], since the sense of **a** is immaterial, both the elastic and viscoelastic contributions of $\psi$, respectively $\psi_e$ and $\psi_{ve}$, must be even functions of **a**. Accordingly, they can be expressed as functions of the second order tensor **a** $\otimes$ **a**. Moreover, as the free energy must be unmodified if the material and its fibers undergo a rotation, then both $\psi_e$ and $\psi_{ve}$ must be scalar-valued isotropic functions of their arguments and can thus be written in terms of the five independent invariants respectively of $\epsilon^e$ and **a**, and of $\epsilon^{ve}$ and **a**. Considering a generic second order symmetric tensor **S**, representing either the elastic strain $\epsilon^e$ or the viscoelastic strain $\epsilon^{ve}$, the five invariants are given by:

$$ 
\begin{aligned}
I_1(\mathbf{S}) &= \mathrm{tr}(\mathbf{S}), & I_2(\mathbf{S}) &= \mathrm{tr}(\mathbf{S}^2), & I_3(\mathbf{S}) &= \mathrm{tr}(\mathbf{S}^3), \\
I_4(\mathbf{S}, \mathbf{a}) &= \mathrm{tr}\{\mathbf{S} \cdot \mathbf{a} \otimes \mathbf{a}\} &= \mathbf{a} \cdot S\mathbf{a}, & & \\
I_5(\mathbf{S}, \mathbf{a}) &= \mathrm{tr}\{\mathbf{S}^2 \cdot \mathbf{a} \otimes \mathbf{a}\} &= \mathbf{a} \cdot S^2\mathbf{a}, & &
\end{aligned}
\quad (7) $$

where tr(•) is the trace of •.

To obtain linear stress–strain relations, both $\psi_e$ and $\psi_{ve}$ must be positive definite quadratic forms of the strain. From relations (7), they will thus only depend on the four invariants ($I_1, I_2, I_4, I_5$), as $I_3$ is of third order in its argument **S**, representing a strain tensor.

For the elastic free energy $\rho\psi_e$ ($\rho$ being the mass density), the most general quadratic form in $\epsilon^e$ reads

$$ \rho\psi_{ve}(\epsilon^{ve}; \mathbf{a}) = \frac{1}{2}\tilde{\lambda}I_1^2(\epsilon^{ve}) + \tilde{\mu}_pI_2(\epsilon^{ve}) + \tilde{\alpha}I_4(\epsilon^{ve}, \mathbf{a})I_1(\epsilon^{ve}) \\ + 2(\tilde{\mu}_L - \tilde{\mu}_p)I_5(\epsilon^{ve}, \mathbf{a}) + \frac{1}{2}\tilde{\beta}I_4^2(\epsilon^{ve}, \mathbf{a}), \quad (9) $$

where $\tilde{\mu}_L, \tilde{\mu}_p, \tilde{\lambda}, \tilde{\alpha}, \tilde{\beta}$ are five constants representing the generalization to 3D of the spring $\tilde{E}$.

The viscoplastic contribution $\rho\psi_{vp}$ is assumed to be a quadratic form in $\alpha_\ell$ and $\alpha_{n\ell}$:

$$ \rho\psi_{vp} = \frac{1}{2}(H_\ell \alpha_\ell : \alpha_\ell + H_{n\ell} \alpha_{n\ell} : \alpha_{n\ell}), \quad (10) $$

with $H_\ell$ and $H_{n\ell}$ being positive scalars.

The Clausius-Duhem inequality under isothermal conditions can be written as:

$$ \sigma : \dot{\varepsilon} - \rho \left( \frac{\partial \psi_e}{\partial \varepsilon^e} : \dot{\varepsilon}^e + \frac{\partial \psi_{ve}}{\partial \varepsilon^{ve}} : \dot{\varepsilon}^{ve} + \frac{\partial \psi_{vp}}{\partial \alpha_\ell} : \dot{\alpha}_\ell + \frac{\partial \psi_{vp}}{\partial \alpha_{n\ell}} : \dot{\alpha}_{n\ell} \right) \geq 0. \quad (11) $$

or, equivalently, as

(8)

$$ \rho\psi_e(\epsilon^e; \mathbf{a}) = \underbrace{\frac{1}{2}\lambda I_1^2(\epsilon^e) + \mu_p I_2(\epsilon^e) + \alpha I_4(\epsilon^e, \mathbf{a}) I_1(\epsilon^e)}_{\text{isotropic contribution}} + 2(\mu_L - \mu_p) I_5(\epsilon^e, \mathbf{a}) + \underbrace{\frac{1}{2}\beta I_4^2(\epsilon^e, \mathbf{a})}_{\text{anisotropic contribution}}, $$

and depends on five material parameters: $\mu_L, \mu_p$ are the shear moduli respectively along the fibers direction **a** (index $L$) and in the plane of isotropy (orthogonal to **a**, index $P$), and $\lambda, \alpha, \beta$ are three material constants that can be derived from the more usual engineering elastic constants $E_L, E_P, \nu_P$, and $\nu_{PL}$, as shown in “Appendix”. These 5 constants come from the generalization to 3D of the spring $E$. The *isotropic contribution*<fnref n="1" /> indicated in relation (8) is added to the *anisotropic contribution* which is due to the presence of the fibers.

The viscoelastic free-energy is defined similarly to (8) as:

$$ \left( \sigma - \rho \frac{\partial \psi_e}{\partial \varepsilon^e} \right) : \dot{\varepsilon} + \rho \left( \frac{\partial \psi_e}{\partial \varepsilon^e} - \frac{\partial \psi_{ve}}{\partial \varepsilon^{ve}} \right) : \dot{\varepsilon}^{ve} + \rho \frac{\partial \psi_e}{\partial \varepsilon^e} : \dot{\varepsilon}^{vp} - \rho \frac{\partial \psi_{vp}}{\partial \alpha_\ell} : \dot{\alpha}_\ell - \rho \frac{\partial \psi_{vp}}{\partial \alpha_{n\ell}} : \dot{\alpha}_{n\ell} \geq 0 \quad (12) $$

As relation (12) must hold for any process and the strain rate $\dot{\varepsilon}$ can be chosen arbitrarily, we have:

$$ \left\{ 
\begin{array}{l}
\sigma = \rho \frac{\partial \psi_e}{\partial \varepsilon^e} \\
\rho \left( \frac{\partial \psi_e}{\partial \varepsilon^e} - \frac{\partial \psi_{ve}}{\partial \varepsilon^{ve}} \right) : \dot{\varepsilon}^{ve} + \rho \frac{\partial \psi_e}{\partial \varepsilon^e} : \dot{\varepsilon}^{vp} \\
- \rho \frac{\partial \psi_{vp}}{\partial \alpha_\ell} : \dot{\alpha}_\ell - \rho \frac{\partial \psi_{vp}}{\partial \alpha_{n\ell}} : \dot{\alpha}_{n\ell} \geq 0
\end{array}
\right. \quad (13) $$

<footnote n="1">¹ Note that, in an isotropic material, this would be the only term in the Helmholtz free-energy.</footnote>
The first condition defines the state law for the stress tensor $\sigma$, that is given by:

$$ \sigma = C : \varepsilon^e, \qquad (14) $$

with C being the elastic stiffness of the material:

$$ C = \lambda I \otimes I + 2\mu_P I + \alpha(a \otimes a \otimes I + I \otimes a \otimes a) + 2(\mu_L - \mu_P)A + \beta a \otimes a \otimes a \otimes a \qquad (15) $$

where I is the second order unit tensor, I is the fourth order unit tensor, and A is a fourth order tensor whose components are given by:

$$ A_{ijkl} = \frac{1}{2}(a_k a_i \delta_{jl} + a_i a_l \delta_{jk} + a_j a_k \delta_{il} + a_j a_l \delta_{ik}). \qquad (16) $$

The second condition in (13) expresses the requirement of non-negativeness of the dissipation and defines the thermodynamic forces associated to the viscoelastic strain $\varepsilon^{ve}$ and to the internal variables $\alpha_{\ell}$ and $\alpha_{n\ell}$. Specifically, one has:

$$ \begin{aligned} \sigma^{ve} &:= \rho \left( \frac{\partial \psi_e}{\partial \varepsilon^e} - \frac{\partial \psi_{ve}}{\partial \varepsilon^{ve}} \right), & X_{\ell} &:= \rho \frac{\partial \psi_{vp}}{\partial \alpha_{\ell}}, \\ X_{n\ell} &:= \rho \frac{\partial \psi_{vp}}{\partial \alpha_{n\ell}}, \end{aligned} \qquad (17) $$

with the latter two representing the current backstresses. From the first of relations (17), accounting for relation (9), one has:

$$ \rho \frac{\partial \psi_{ve}}{\partial \varepsilon^{ve}} = \tilde{C} : \varepsilon^{ve}, \qquad (18) $$

with

$$ \tilde{C} = \tilde{\lambda} I \otimes I + 2\tilde{\mu}_P I + \tilde{\alpha}(a \otimes a \otimes I + I \otimes a \otimes a) + 2(\tilde{\mu}_L - \tilde{\mu}_P)A + \tilde{\beta}a \otimes a \otimes a \otimes a. \qquad (19) $$

Finally, from the last two relations in (17), using relations (10) one finds

$$ X_{\ell} = H_{\ell} \alpha_{\ell}, \quad X_{n\ell} = H_{n\ell} \alpha_{n\ell}, \qquad (20) $$

corresponding to linear relations between the internal variables describing the kinematic hardening behavior, and their associated thermodynamic forces.

To automatically fulfill the Clausius-Duhem inequality (12), we adopt the hypothesis of generalized normality and we postulate the existence of a non-negative convex and zero at the origin dissipation potential $\Phi$, defined as follows:

$$ \begin{aligned} \Phi &= \Phi(\dot{\varepsilon}^{ve}, \dot{\varepsilon}^{vp}, \dot{\alpha}_{\ell}, \dot{\alpha}_{n\ell}:a, \alpha_{n\ell}) \\ &= \Phi^{ve}(\dot{\varepsilon}^{ve};a) + \Phi^{vp}(\dot{\varepsilon}^{vp}, \dot{\alpha}_{\ell}, \dot{\alpha}_{n\ell};a, \alpha_{n\ell}). \end{aligned} \qquad (21) $$

In particular, assuming a linear viscoelastic material response, the dissipation potential $\Phi^{ve}$ must be a positive definite quadratic form of $\dot{\varepsilon}^{ve}$, while the dependence on the fiber direction is introduced again through the invariants defined by relations (7), such that:

$$ \begin{aligned} \Phi^{ve}(\dot{\varepsilon}^{ve};a) = \frac{1}{2} &\theta_{\lambda} \tilde{\lambda} I_1^2(\dot{\varepsilon}^{ve}) + \theta_{\mu_P} \tilde{\mu}_P I_2(\dot{\varepsilon}^{ve}) + \theta_{\alpha} \tilde{\alpha} I_4(\dot{\varepsilon}^{ve}, a) I_1(\dot{\varepsilon}^{ve}) \\ &+ 2(\theta_{\mu_L} \tilde{\mu}_L - \theta_{\mu_P} \tilde{\mu}_P) I_5(\dot{\varepsilon}^{ve}, a) + \frac{1}{2} \theta_{\beta} \tilde{\beta} I_4^2(\dot{\varepsilon}^{ve}, a) \end{aligned} \qquad (22) $$

with $\theta_{\mu_L}, \theta_{\mu_P}, \theta_{\lambda}, \theta_{\alpha}, \theta_{\beta}$ being characteristic retardation times, and $\tilde{\mu}_L, \tilde{\mu}_P, \tilde{\lambda}, \tilde{\alpha}, \tilde{\beta}$ being the five material constants used for the viscoelastic free-energy (9).

The derivation of $\Phi^{ve}$ gives the thermodynamic force $\sigma^{ve}$ defined in (17) and associated with $\varepsilon^{ve}$, such that:

$$ \sigma^{ve} = \frac{\partial \Phi^{ve}}{\partial \dot{\varepsilon}^{ve}} = \Xi : \dot{\varepsilon}^{ve}, \qquad (23) $$

with

$$ \Xi = \theta_{\lambda} \tilde{\lambda} I \otimes I + 2\theta_{\mu_P} \tilde{\mu}_P I + \theta_{\alpha} \tilde{\alpha}(a \otimes a \otimes I + I \otimes a \otimes a) + 2(\theta_{\mu_L} \tilde{\mu}_L - \theta_{\mu_P} \tilde{\mu}_P)A + \theta_{\beta} \tilde{\beta}a \otimes a \otimes a \otimes a. \qquad (24) $$

By combining relations (17), (18), (13) and (23), one has:

$$ \dot{\varepsilon}^{ve} = \Xi^{-1} : (\sigma - \tilde{C} : \varepsilon^{ve}). \qquad (25) $$

The viscoplastic contribution to the dissipation process is more easily described in terms of the dual potential $\Phi^{vp}$, obtained through the Legendre-Fenchel transform of the dissipation potential $\Phi^{vp}$. Accordingly, one has:

$$ \dot{\varepsilon}^{vp} = \frac{\partial \Phi^{vp*}}{\partial \sigma}, \quad \dot{\alpha}_{\ell} = - \frac{\partial \Phi^{vp*}}{\partial X_{\ell}}, \quad \dot{\alpha}_{n\ell} = - \frac{\partial \Phi^{vp*}}{\partial X_{n\ell}}. \qquad (26) $$

To consider an elastic domain, the dissipation potential must depend on the positive part $\langle\varphi\rangle$ of the activation (or yield) function $\varphi$. Furthermore, to account for a nonlinear kinematic hardening in the framework of generalized standard materials, it is necessary to introduce an additional “recall term” depending on $X_{n\ell}$ and the internal variables $\alpha_{n\ell}$ in the expression of the potential (see [30–32]), considering it as a parameter.
Different expressions for the dual dissipation potential can be assumed. A common choice, here adopted to reduce the number of material parameters, is to define a quadratic potential, leading to a 3D extension of the Bingham model (with nonlinear kinematic hardening):

$$ \Phi^{vp*}(\sigma, X_{\ell}, X_{n\ell}; a, \alpha_{n\ell}) = \frac{1}{2\eta} \left\langle \varphi(\sigma, X_{\ell}, X_{n\ell}; a) + \frac{B}{2} \left( \frac{1}{H_{n\ell}} X_{n\ell} : X_{n\ell} - \alpha_{n\ell} : \alpha_{n\ell} \right) \right\rangle^2, \quad (27) $$

where B is a positive scalar material parameter. Note that the term added to $\varphi$ in (27) vanishes when the second of the hardening laws (20) is fulfilled.

As proposed in [21], it can be accepted that plastic deformations start when the stress in the fibers direction exceeds a certain value. Accordingly, we choose the following expression of the yield function

$$ \varphi = \varphi(\sigma, X_{\ell}, X_{n\ell}; a) = \underbrace{|(\sigma - X_{\ell} - X_{n\ell}) : (a \otimes a)|}_{\sigma_{eq}} - \sigma_y, \quad (28) $$

where the equivalent stress $\sigma_{eq}$ represents the modulus of the projection of the over-stress $(\sigma - X_{\ell} - X_{n\ell})$ along the fibers, with $\sigma_y$ being the initial yield stress in tension of the filament material. Note that a yield function similar to the one above (relation (28)) was also used in [33, 34] to describe the anisotropic elasto-plastic behavior of paper and paperboard.

The (associative) flow rule for the viscoplastic strain thus becomes:

$$ \dot{\varepsilon}^{vp} = \frac{\langle\varphi\rangle}{\eta} \operatorname{sign}((\sigma - X_{\ell} - X_{n\ell}) : (a \otimes a))(a \otimes a), \quad (29) $$

where $\operatorname{sign}(\cdot)$ gives the sign of $\cdot$. The term $\frac{\langle\varphi\rangle}{\eta}$ gives the magnitude of the visco-plastic strains and hence represents the visco-plastic multiplier. Similarly, the evolution laws for the kinematic internal variables are obtained from the second of (26) with $\Phi^{vp*}$ defined by (27) and read

$$ \left\{ \begin{array}{l} \dot{\alpha}_{\ell} = -\frac{\partial \Phi^{vp*}}{\partial X_{\ell}} = \dot{\varepsilon}^{vp}, \\ \dot{\alpha}_{n\ell} = -\frac{\partial \Phi^{vp*}}{\partial X_{n\ell}} = \dot{\varepsilon}^{vp} - \frac{\langle\varphi\rangle}{\eta} \frac{B}{H_{n\ell}} X_{n\ell}. \end{array} \right. \quad (30) $$

The nonlinear kinematic hardening here described was initially introduced by Armstrong and Frederick [35], and then further developed by other authors [36]. The non-negativeness of the dissipation can be proved, as shown e.g. in [30]. One can easily retrieve a linear kinematic hardening by setting the scalar B equal to 0.

**Remark 1** The associated flow rule (29) assumes that no plastic deformations are present in the plane of isotropy of the transversely isotropic material. In the absence of experimental data, this is the simplest choice. Accordingly, one can define a plastic Poisson’s coefficient giving the ratio between lateral and axial plastic deformations, that is null for the present associative model. If this was not the case, then one could consider a non-associative model to correctly account for the volumetric deformations during the plastic phase.

**Remark 2** With the above model, the current state of the microscopic rearrangements interpreted as kinematic hardening is considered to be independent of the material direction a, whereas their evolution is considered to be active only in the material direction. This is due to the particular choice for the yield function $\varphi$.

## 5 Parameters identification

The presented model has a total of 20 material parameters. These are 5 constants $\mu_L, \mu_p, \lambda, \alpha, \beta$ (or $\mu_L, \mu_p, E_L, E_p, v_{PL}$) governing the instantaneous elastic behavior, 10 constants $\tilde{\mu}_L, \tilde{\mu}_p, \tilde{\lambda}, \tilde{\alpha}, \tilde{\beta}$ and $\theta_{\mu_L}, \theta_{\mu_p}, \theta_{\lambda}, \theta_{\alpha}, \theta_{\beta}$ describing the visco-elastic response of the material, 1 initial yield stress $\sigma_y$, 1 visco-plastic parameter $\eta$, and the 3 scalars $H_{\ell}, H_{n\ell}$ and B describing the kinematic hardening.

Due to the anisotropic characteristic of the constitutive model, the identification of 15 independent visco-elastic material constants is unfeasible if based solely on the uniaxial tensile tests. Several assumptions can be made in order to reduce the number of these parameters. We here describe a possible simplifying strategy. The idea is that the non-isotropic effects are due to the fiber-inclination and hence are similar for the elastic and visco-elastic behavior. In particular one can assume that the 5 constants $\tilde{\mu}_L, \tilde{\mu}_p, \tilde{\lambda}, \tilde{\alpha}, \tilde{\beta}$ be proportional to the corresponding
$\mu_L$, $\mu_P$, $\lambda$, $\alpha$, $\beta$ through a single scalar parameter $\tilde{C}$ which can be interpreted as the ratio of the spring stiffness $\tilde{E}/\tilde{E}_E$ of the 1D model. Namely:

$$ \frac{\tilde{\mu}_L}{\mu_L} = \frac{\tilde{\mu}_P}{\mu_P} = \frac{\tilde{\lambda}}{\lambda} = \frac{\tilde{\alpha}}{\alpha} = \frac{\tilde{\beta}}{\beta} = \tilde{C} \quad (31) $$

Similarly for the viscous effect, a single characteristic retardation time $\theta$, which can be interpreted as the ratio $\tilde{\eta}/\tilde{E}$, can be postulated:

$$ \theta_{\mu_L} = \theta_{\mu_P} = \theta_{\lambda} = \theta_{\alpha} = \theta_{\beta} = \theta. \quad (32) $$

With the above assumptions the independent visco-elastic parameters are reduced to 7. Note that often two characteristic times are introduced even in the isotropic case to be identified with tension and shear creep tests. Of course this could also be done in the present case, but due to lack of experimental data we keep the previously described simplification.

The total number of free parameters of the material model for multi-filament yarns is hence reduced to 12. Their identification is performed in several steps.

First, the uniaxial tests performed at different strain rates on conditioned yarn Y0 specimens are considered. The stresses are obtained from the experimental data by dividing the force by an area of 0.17 mm² as shown in Fig. 3b. The initial elastic slope allows to fix a condition for $E_L$ and $\tilde{E}_L$; the elastic limit at low strain rate gives $\sigma_y$; the slope of the hardening phase depends on $H_\ell$, while $H_{n\ell}$ modifies the initial nonlinear inelastic phase; the difference between the almost parallel hardening curves obtained at different strain rates gives the viscoplastic coefficient $\eta$.

It should be noted that, even though the identification procedure on the basis of the uniaxial tests is not unique, the above described steps can be used as a guideline in the procedure.

Four independent elastic constants governing the transversal multiaxial instantaneous material behavior remain to be determined. These are, for instance, $E_P$, $\nu_P$, $\nu_{LP}$ and $\mu_L$. As the stress–strain curves for yarn Y0 are independent of these parameters, other tests, e.g. those on twisted yarns, have to be considered.

Actually, due to the fiber inclination, the uniaxial force applied to twisted yarns produces normal stresses in the fibers’ direction and in the transversal direction, together with shear stresses between the fibers. Accordingly, all elastic constants influence the reinforcement behavior. The Poisson’s coefficients, due to the

discontinuous nature of a fiber bundle and for the sake of simplicity, are assumed to be zero. Different pairs of $E_P$ and $\mu_L$ values can be identified to obtain a good agreement with the uniaxial experiments on twisted yarns.

**Remark 3** The available experimental data are not sufficient to ensure the uniqueness of the parameter identification. Other multiaxial tests would be required to carry out a more accurate parameter identification. A sensitivity analysis quantifying the influence of each material parameter on the final output, as those performed in elastic regime e.g. in [37], could give further guidance in the parameter identification process. This study is however outside the scope of this work and will be pursued elsewhere.

# 6 Numerical analyses

## 6.1 Implementation of the model

The non-linear elasto-viscoplastic model proposed in the present work has been implemented as a user routine UMAT in the finite element code ABAQUS.

The algorithm is based on the standard elastic predictor/plastic corrector format. The classic implicit backward difference scheme is assumed for time integration of the constitutive law. Let us consider a time interval $[t_n, t_{n+1}]$ at a generic Gauss point of the finite element mesh. Given the incremental strain $\Delta\epsilon$, and the values $\sigma, \epsilon^{ve}, \epsilon^{vp}, \alpha_\ell, \alpha_{n\ell}$ at time $t_n$, the algorithm computes the updated quantities at time $t_{n+1}$. Relations (14), (25), (29), (20), (30), and (28) govern the viscoelastic-viscoplastic rate constitutive law. The backward difference integration procedure gives the equations of the finite step:

$$ \left\{ \begin{aligned} \sigma_{n+1} &= \sigma_n + C : (\Delta\epsilon - \Delta\epsilon^{ve} - \Delta\epsilon^{vp}) \\ \Delta\epsilon^{ve} &= \Xi^{-1} : (\sigma_{n+1} - \tilde{C} : \epsilon_{n+1}^{ve})\Delta t \\ \Delta\epsilon^{vp} &= \frac{\langle\varphi_{n+1}\rangle}{\eta} \operatorname{sign}((\sigma_{n+1} - X_{\ell^{n+1}} - X_{n\ell^{n+1}}) : (a \otimes a))(a \otimes a)\Delta t \\ X_{\ell^{n+1}} &= H_\ell (\alpha_{\ell^n} + \Delta\alpha_\ell) \\ X_{n\ell^{n+1}} &= H_{n\ell} (\alpha_{n\ell^n} + \Delta\alpha_{n\ell}) \\ \Delta\alpha_\ell &= \Delta\epsilon^{vp} \\ \Delta\alpha_{n\ell} &= \Delta\epsilon^{vp} - \frac{\langle\varphi_{n+1}\rangle}{\eta} \frac{B}{H_{n\ell^{n+1}}} X_{n\ell} \Delta t \\ \varphi_{n+1} &= |(\sigma_{n+1} - X_{\ell^{n+1}} - X_{n\ell^{n+1}}) : (a \otimes a)| - \sigma_y \end{aligned} \right. \quad (33) $$
In the first phase (predictor) the numerical integration algorithm evaluates the elastic trial state, assuming the increment to be purely viscoelastic. The next phase (corrector) is a check for plastic consistency. If the assumption made during the trial step is correct, then $\sigma$ and $\varepsilon^{ve}$ are updated at time $t_{n+1}$, while $\varepsilon^{vp}$, $\alpha_{\ell}$, $\alpha_{n\ell}$ are kept constant from time $t_n$. Otherwise, a return mapping step is applied and all the values $\sigma$, $\varepsilon^{ve}$, $\varepsilon^{vp}$, $\alpha_{\ell}$, $\alpha_{n\ell}$ are updated at time $t_{n+1}$.

The single yarn is then modeled as a solid with a cylindrical shape, discretized by 3D finite elements.

The local fiber direction $\boldsymbol{a}$, depending on the yarn twist and on the distance between the considered point and the yarn axis according to relation (2), is assigned at each Gauss point, once for all, at the beginning of the analysis.

In practical engineering applications, as in the case of rubber-reinforced composites where multiple yarns are present, the association of the correct fiber direction to each Gauss point of a yarn FE discretization can still be performed. In fact, assuming the yarn axis geometry as known, the distance of each Gauss point in a yarn from the pertinent axis would still be computable, thus allowing the determination of the local fiber direction.

## 6.2 Finite element simulations

The available uniaxial experimental tests on untwisted and twisted yarns have been considered. First the material parameters are identified using the strategy explained in Sect. 5 and then used to simulate all tests.

For this, a 2 mm long fiber-reinforced cylindrical bar with circular cross-section was modeled in Abaqus (see Fig. 7). The nodes at the bottom face of

the bar were restrained in tangential and longitudinal directions, but not in radial direction. The nodes at the top face were subjected to an imposed longitudinal displacement (equal value for all nodes) while being restrained in the tangential direction and free in the radial direction. Twenty-noded solid elements with reduced integration (C3D20R in Abaqus) were used. The new material model, implemented in a UMAT, was assigned to the solid and fiber inclination was given at each Gauss point.

The results of the model are compared in Fig. 8 with some of the experimental results.

In Fig. 8a the experimental (dashed lines) and numerical (continuous lines) results for monotonic tests on yarn Y0 in humid conditions are plotted for different strain rates. Note that the experimental curves representing average values (relevant to 20 specimens) are the same already shown in Fig. 3b.

Figure 8b contains a comparison between the experimental data (shaded regions) and the numerical results (continuous curves) for the monotonic uniaxial tensile tests of yarns Y0 and Y48 in humid conditions (respectively Y0H and Y48H), at a strain rate of 0.017 s⁻¹.

The values of the 12 parameters of the model identified and used for the simulations of humid rayon are reported in Table 2.

The results for twisted yarns are given in terms of force versus strain as, in contrast with the result valid for the untwisted yarn, the cross-sectional axial stress for the twisted yarns is no more uniform. This can be seen in Fig. 9a, showing the axial stress in the yarn Y48H, for a total axial strain $\bar{\varepsilon} = 0.1$. In particular, the axial stress is higher when moving towards the core of the yarn, as the fibers get more and more aligned with the yarn axis (and thus with the direction of loading).

<figure><img src="image_7.png" /><figcaption>**Fig. 7 Finite element discretization of yarn Y48; the red arrows represent the directions of the fibers**</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 8** Uniaxial tensile tests of yarns in humid conditions a untwisted yarn Y0H: comparison between experimental data (dashed curves) and numerical results (continuous curves) at strain rates 0.017 s⁻¹, 0.01 s⁻¹, 0.003 s⁻¹, 0.0017 s⁻¹ and</figcaption></figure>

0.0003 s⁻¹; b yarns Y0H and Y48H at strain rates of 0.017 s⁻¹: comparison between experimental data (shaded regions) and numerical results (continuous curves)

<table><caption>**Table 2** Identified material constants (humid conditions)</caption><thead><tr><th>E<sub>p</sub> [MPa]</th><th>E<sub>L</sub> [MPa]</th><th>ν<sub>p</sub> [-]</th><th>ν<sub>LP</sub> [-]</th><th>μ<sub>L</sub> [MPa]</th><th>C̄ [-]</th></tr></thead><tbody><tr><td>450</td><td>13000</td><td>0</td><td>0</td><td>10</td><td>5.5</td></tr><tr><td>θ [s]</td><td>σ<sub>y</sub> [MPa]</td><td>H<sub>τ</sub> [MPa]</td><td>H<sub>ητ</sub> [MPa]</td><td>B [-]</td><td>η [s MPa]</td></tr><tr><td>2.80</td><td>50</td><td>3570</td><td>40000</td><td>1000</td><td>3000</td></tr></tbody></table>

This response results in higher visco-plastic deformations around the core of the yarn, as shown in Fig. 9a, containing a contour plot of the equivalent visco-plastic strain $\varepsilon^{eq}$, measured up to the total axial strain $\bar{\varepsilon}$ and defined as the integral in time of $\|\dot{\varepsilon}^{vp}\|$.

Tomographic measurements of the radius variation of a yarn Y48 during a uniaxial tensile test were performed in the Pirelli laboratory. The experimental data are reported in Fig. 10. More specifically, Fig. 10a shows the microscopic images of the yarn Y48H cross-section during the test, captured at 3 levels of axial strain. The tomography gives an initial value of the radius $R_0 = 0.24$ mm, in

reasonable agreement with, but slightly higher than the one obtained for the area measured through the microscopy reported in Fig. 5 ($R_0 = 0.233$ mm) and used in the simulations. The variation of the radius, normalized with the initial value, for varying imposed axial strain is reported in Fig. 10b (markers), together with the numerical results (continuous curve). The model is able to capture, at least qualitatively, the experimentally observed radius variation.

As apparent from Fig. 10b, the model predicts a nonlinear radial contraction during a uniaxial test. As pointed out also in [21], this is due to the chosen flow potential (27). In particular, the (associative) flow rule (29), derived from the flow potential, depends on the backstress. This causes a nonlinear relation between the strain along the direction of reinforcement and the strain in the plane perpendicular to it (plane of isotropy). The backstress tensor is indeed dependent on the viscoplastic strain tensor $\varepsilon^{vp}$ (as can be checked by using relations (20) and (30)). During a tensile

<figure><img src="image_9.png" /><figcaption>**Fig. 9** Contour plots from the numerical simulation of yarn Y48H, at the axial strain $\bar{\varepsilon} = 0.1$: a axial stress ($\sigma_3$) in MPa; b equivalent (visco-)plastic strain ($\varepsilon^{eq}$). (Color figure online)</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 10** Radial contraction during a uniaxial tensile test of yarn Y48H: a micro-tomography images; b comparison of radius versus axial strain from experimental tests (markers) and numerical result (continuous)</figcaption></figure>

test on a bundle of parallel fibers, the evolution of the plastic strain components in the plane of isotropy will thus depend on their instantaneous values, generating a nonlinear radial/axial strain relation.

To check the validity of the presented model, different loading conditions, twist and moisture levels have been considered.

In particular the non-monotonic test on yarn Y0 of Fig. 4 has been simulated: the numerical prediction (shown by a continuous line in Fig. 11) is in good agreement with the experimental results obtained for five specimens (dashed lines).

In Fig. 12, we show a comparison between experimental data (shaded regions) and numerical results for a tensile uniaxial test of a yarn Y48H, carried out at strain rates of 0.017 s⁻¹, 0.003 s⁻¹, and 0.0003 s⁻¹.

The material model is able to predict the strain rate dependence that is particularly visible in the elastic-viscoplastic range.

The yarns were also tested in dry conditions (c.f. Sect. 2). With reference to such conditions, Fig. 13 compares the experimental and numerical curves at different twist levels.

The numerical curves are obtained by assuming a linear dependence of a generic material constant p on the degree of saturation $S_r$:

$$p(S_r) = p_d(1 - S_r) + p_h S_r \quad (34)$$

where $p_d = p(0)$ and $p_h = p(1)$ are the values of the parameter in dry and fully saturated conditions, respectively (see e.g. [6]). The fixed material constants in dry conditions are collected in Table 3.

The experimental responses, see Fig. 13, show that the twisting makes the yarns more compliant, with

<figure><img src="image_11.png" /><figcaption>**Fig. 11** Uniaxial non-monotonic test on yarn Y0H at strain rate 0.003 s⁻¹: comparison between experimental data (dashed curves) and numerical results (continuous curves)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 12** Tensile response of yarn Y48H at strain rates 0.017 s⁻¹, 0.003 s⁻¹, and 0.0003 s⁻¹. Comparison of the experimental (shaded regions) and numerical (continuous) responses</figcaption></figure>

a significant decrease of the initial modulus. This is captured fairly well by the geometrical model used for the numerical simulations.
<figure><img src="image_13.png" /><figcaption>**Fig. 13** Tensile response of yarns with different twist levels at a strain rate of 0.017 s⁻¹ in dry conditions. Comparison of the experimental (dashed) and numerical (continuous) responses</figcaption></figure>

<table><caption>**Table 3** Identified material constants (dry conditions)</caption><thead><tr><th></th><th>E<sub>p</sub> [MPa]</th><th>E<sub>L</sub> [MPa]</th><th>ν<sub>p</sub> [-]</th><th>ν<sub>LP</sub> [-]</th><th>μ<sub>L</sub> [MPa]</th><th>ζ̃ [-]</th></tr></thead><tbody><tr><td>800</td><td>17617</td><td>0</td><td>0</td><td>10</td><td>6.7871</td></tr><tr><td>θ [s]</td><td>σ<sub>y</sub> [MPa]</td><td>H<sub>c</sub> [MPa]</td><td>H<sub>nε</sub> [MPa]</td><td>B [-]</td><td>η [s MPa]</td></tr><tr><td>1.67</td><td>80.9</td><td>4500</td><td>73839</td><td>492.27</td><td>3000</td></tr></tbody></table>

**Remark 4** In Fig. 12, the inelastic response of the yarn Y48 in humid conditions is characterized by a nonlinear behavior showing a self-stiffening response, which becomes more evident at large strains. We interpret the hardening effect as a reorientation of the fibers along the direction of the stretch, as it is typically observed in polymeric bulky materials due to the alignment of the polymer network [38]. This particular behavior seems to take place only in humid conditions, as it is absent in the experimental curves in dry conditions shown in Fig. 13 for twisted yarns.

**Remark 5** The knowledge of the local fiber inclination is essential for the proposed modeling. As discussed in Sect. 3 the coaxial helices model correctly reproduces the radial dependence for a given twist, but direct experimental measures show that there is a variation around this value. Furthermore, the actual twist of the yarn can be different from the nominal one. To see the sensitivity of the numerical results with respect to this parameter, we performed simulations of yarn Y48H (humid conditions) by changing

<figure><img src="image_14.png" /><figcaption>**Fig. 14** Numerical response of yarn Y48H under given monotonically increasing strain for ±3% variation of twist level</figcaption></figure>

the twist of ±3%. Figure 14 shows that the resulting variation of the force is around 3% as well, with higher axial force for smaller twist values.

**Remark 6** To check the overall consistency of the proposed procedure, three different meshes have been used for discretizing the cylindrical body representing yarn Y48H. These meshes, shown in Fig. 15, are characterized by a different structure in the cross-section and different element sizes. The three corresponding FE analyses produce global force-strain curves which are indistinguishable. The distributions of normal longitudinal stresses are in very good agreement as shown by the contour maps of Fig. 15 relevant to an imposed strain of $\bar{\epsilon} = 0.13$. A more detailed examination of such stresses at the Gauss points nearest to the yarn axis shows a maximum relative difference of 1.3% for three adopted discretizations.

## 7 Conclusions

In this paper, a new anisotropic viscoelastic-visco-plastic constitutive model for materials composed of a single family of continuous fibers with different orientations has been developed. In particular, the model is conceived to describe at the macro-scale the complex mechanical behavior of yarns and cords made of thousands of micrometric fibers of rayon. The peculiar geometry of twisted yarns is taken into account by prescribing at any point the fiber direction. The chosen continuum theory is based on the additive decomposition of the strain tensor into elastic, viscoelastic
<figure><img src="image_15.png" /><figcaption>**Fig. 15** Numerical response for three different FE discretizations for yarn Y48H (polar mesh, orthogonal mesh and unstructured mesh) at a strain of $\epsilon = 0.13$. (Color figure online)</figcaption></figure>

and viscoplastic parts. The approach here developed allows to include, in a unified manner, the effect of the peculiar geometry of the yarns and the non-linear, dissipative behavior of the single fibers.

The elastic part of the deformation is described by an anisotropic Helmholtz free-energy function, which is written using a fixed global system of reference. The main advantage of using this method is that no local reference systems have to be defined, as the stiffness can be written at a material point concerning the fixed global reference system.

The yield surface initiating the viscoplastic behavior is modeled by employing the stress measured along the fiber's direction. The material parameters needed by the model are obtained by fitting uniaxial force versus strain data on untwisted and twisted yarns. In this work, we focused on the development of a model to be used for yarns subjected to different loading conditions, twists and moisture levels. With the fitted parameters, the model was shown to be able to reproduce the overall behavior of the experimental responses.

The model allows us to interpret the effects of twist experimentally observed on the macroscopic uniaxial response, namely, by increasing the twist, (i) the overall elastic stiffness decreases; (ii) the limit force of the linear elastic behavior decreases; (iii) the change of slope between elastic and elasto-viscoplastic parts is more gradual. The experimental evidence (i) is interpreted as a combined effect of transverse isotropy (with higher stiffness in fiber

direction) and fiber inclination inside the yarn (fibers with higher inclination giving a lower contribution to the overall axial stiffness). The experimental evidences (ii) and (iii) can be explained by the non-uniform inclination of fibers and hence a non-uniform state of stress inside the yarn. The fibers close to the axis, which are less inclined, are subject to higher stress and hence enter into the plastic regime first, for lower values of the overall force, then gradually also more external fibers reach the yield limit progressively.

Even though in this work explicit reference is always made to polymeric yarns, the new transversally isotropic viscoelastic-viscoplastic model can also be employed for other materials characterized by a single family of continuous fibers.

The proposed model can be improved in at least two directions: by considering large deformations, and/or by introducing a reorientation of the fibers. This latter extension of the model will be pursued as a next step.

Albeit the aforementioned developments could improve the predictions, the current model can already form the basis of a three-dimensional viscoelastic-viscoplastic continuum mechanical model to be used in finite element procedures for the simulation of the mechanical behavior of yarn-reinforced rubber composites. A similar approach can also be used for multi-ply yarns, provided that the correct geometrical description of fiber trajectories is inserted into the model.
## Appendix: elastic constants

We report in the following the relations between the engineering constants $E_L$, $E_P$, $v_{PL}$, $v_P$, and the coefficients $\lambda$, $\alpha$, $\beta$, $\mu_P$ of the invariant-based approach used to define the free-energy (5):

$$ 
\begin{aligned}
\lambda &= \frac{E_P(v_P + v_{PL}v_{LP})}{\text{den}} \\
\alpha &= \frac{E_P(v_{LP}(1 + v_P - v_{PL}) - v_P)}{\text{den}} \\
\beta &= \frac{E_L(1 - v_P^2)}{\text{den}} - (2\alpha + \lambda) - 4\mu_L + 2\mu_P \\
\mu_P &= \frac{E_P}{2(1 + v_P)}
\end{aligned}
\quad (35) $$

with

$$ \text{den} = 1 - v_P^2 - 2v_{PL}v_{LP} - 2v_{PL}v_{LP}v_P \quad (36) $$

and

$$ v_{PL}/E_P = v_{LP}/E_L. \quad (37) $$

In the above relations, $E_L$ is the Young's modulus along the direction of reinforcement, $v_{LP}$ is the Poisson's coefficient governing the deformation in the plane of isotropy due to a stress along the direction of reinforcement, $v_{PL}$ is the Poisson's coefficient governing the deformation in the direction of reinforcement due to a stress in the plane of isotropy, $E_P$ and $v_P$ are respectively the Young's modulus and Poisson's coefficient in the plane of isotropy.

Relations (35) can be found by expressing the material stiffness tensor first with the engineering constants $E_L$, $E_P$, $v_{PL}$, $v_P$, $\mu_L$ and then with the coefficients $\lambda$, $\alpha$, $\beta$, $\mu_L$, $\mu_P$, fixing a certain direction of reinforcement $\alpha$. Since both stiffnesses describe the same transversely isotropic material, their components must be equal.

**Funding** Open access funding provided by Politecnico di Milano within the CRUI-CARE Agreement.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated

otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

1. Rodgers B (2015) Rubber compounding chemistry and applications, vol 624, 2nd edn. CRC Press. https://doi.org/10.1201/b18931

2. Horrocks AR (2016) Technical textiles in transport (land, sea, and air). In: Horrocks AR, Anand SC (2016). Handbook of Technical Textiles: Vol. 2. Duxford, UK: Woodhead Publishing. https://doi.org/10.1016/B978-1-78242-465-9.00011-2

3. Sodhani D, Reese S, Moreira R, Jockenhoevel S, Mela P, Stapleton SE (2016) Multi-scale modelling of textile reinforced artificial tubular aortic heart valves. Meccanica 52(3):677–693. https://doi.org/10.1007/s11012-016-0479-y

4. Mendes ISF, Prates A, Evtuguin DV (2021) Production of rayon fibres from cellulosic pulps: state of the art and current developments. Carbohydr Polym 273:118466. https://doi.org/10.1016/j.carbpol.2021.118466

5. Teodosio L, Alferi G, Genovese A, Farroni F, Mele B, Timpone F, Sakhnevych A (2021) A numerical methodology for thermo-fluid dynamic modelling of tyre inner chamber: towards real time applications. Meccanica 56(3):549–567. https://doi.org/10.1007/s11012-021-01310-w

6. Caracino P, Agresti S, Pires da Costa L, Novati G, Comi C (2022) The nonlinear behaviour of cords to be used in rayon-rubber composites. In: Marano C, Briatico Vangosa F, Andena L, Frassine R (eds) Constitutive models for rubber XII, CRC Press, pp 472–476. https://doi.org/10.1201/9781003310266-77

7. Morton WE, Yen KC (1952) The arrangement of fibres in fibro yarns. J Text Inst Trans 43(2):60–66. https://doi.org/10.1080/19447025208659646

8. Sibellas A, Adrien J, Durville D, Maire E (2020) Experimental study of the fiber orientations in single and multi-ply continuous filament yarns. J Text Inst 111(5):646–659. https://doi.org/10.1080/00405000.2019.1659471

9. Miao Y, Zhou E, Wang Y, Cheeseman BA (2008) Mechanics of textile composites: micro-geometry. Compos Sci Technol 68(7–8):1671–1678. https://doi.org/10.1016/j.compscitech.2008.02.018

10. Durville D (2010) Simulation of the mechanical behaviour of woven fabrics at the scale of fibers. Int J Mater Form 3(S2):1241–1251. https://doi.org/10.1007/s12289-009-0674-7

11. Vassiliadis S, Kallivretaki A, Provatidis C (2010) Mechanical modelling of multifilament twisted yarns. Fibers Polym 11(1):89–96. https://doi.org/10.1007/s12221-010-0089-6

12. Durville D, Baydoun I, Moustacas H, Périé G, Wielhorski Y (2018) Determining the initial configuration and characterizing the mechanical properties of 3d angle-interlock
fabrics using finite element simulation. Int J Solids Struct 154:97–103. https://doi.org/10.1016/j.ijsolstr.2017.06.026

13. Reese S (2003) Meso-macro modelling of fibre-reinforced rubber-like composites exhibiting large elastoplastic deformation. Int J Solids Struct 40(4):951–980. https://doi.org/10.1016/S0020-7683(02)00602-9

14. Meschke G, Helmwein P (1994) Large-strain 3d-analysis of fibre-reinforced composites using rebar elements: hyperelastic formulations for cords. Comput Mech 13(4):241–254. https://doi.org/10.1007/bf00350227

15. Korunović N, Fragassa C, Marinković D, Vitković N, Trajanović M (2019) Performance evaluation of cord material models applied to structural analysis of tires. Compos Struct 224:111006. https://doi.org/10.1016/j.compstruct.2019.111006

16. Zreid I, Kaliske M (2022) Modeling thermal shrinkage of tire cords and its application in FE analysis of post cure inflation. Finite Elem Anal Des 205:103757. https://doi.org/10.1016/j.finel.2022.103757

17. Treloar LRG, Riding G (1963) la theory of the stress-strain properties of continuous-filament yarns. J Text Inst Trans 54(4):156–170. https://doi.org/10.1080/194447026308660166

18. Sibellas A, Rusinowicz M, Adrien J, Durville D, Maire E (2021) The importance of a variable fibre packing density in modelling the tensile behaviour of single filament yarns. J Text Inst 112(5):733–741. https://doi.org/10.1080/00405000.2020.1781347

19. Poilâne C, Cherif ZE, Richard F, Vivet A, Doudou BB, Chen J (2014) Polymer reinforced by flax fibres as a viscoelastoplastic material. Compos Struct 112:100–112. https://doi.org/10.1016/j.compstruct.2014.01.043

20. Richard F, Poilâne C, Yang H, Gehring F, Renner E (2018) A viscoelastoplastic stiffening model for plant fibre unidirectional reinforced composite behaviour under monotonic and cyclic tensile loading. Compos Sci Technol 167:396–403. https://doi.org/10.1016/j.compscitech.2018.08.020

21. Bedzra R, Reese S, Simon J-W (2020) Hierarchical multi-scale modelling of flax fibre/epoxy composite by means of general anisotropic viscoelastic-viscoplastic constitutive models: part I—micromechanical model. Int J Solids Struct 202:58–74. https://doi.org/10.1016/j.ijsolstr.2020.05.020

22. Popa CM, Gebhardt C, Raje N, Steenwyk B, Kaliske M (2020) Investigation of cord-rubber composite durability by the material force method. Eng Fract Mech 229:106909. https://doi.org/10.1016/j.engfracmech.2020.106909

23. Abu Bakar IA, Kramer O, Bordas S, Rabczuk T (2013) Optimization of elastic properties and weaving patterns of woven composites. Compos Struct 100:575–591. https://doi.org/10.1016/j.compstruct.2012.12.043

24. BISFA Standard (2007) Testing methods for viscose, cupro, acetate, triacetate and lyocell filament yarns

25. Platt MM (1950) Mechanics of elastic performance of textile materials. Text Res J 20(1):1–15. https://doi.org/10.1177/004051755002000101

26. Pires da Costa L, Novati G, Caracino P, Comi C, Agresti S (2022) Tensile behaviour of rayon cords in different conditions. Materials Research Proceedings, Vol. 26, pp 133–138. https://doi.org/10.21741/9781644902431-22

27. Curnier A, He Q-C, Zysset P (1995) Conewise linear elastic materials. J Elast 37(1):1–38. https://doi.org/10.1007/bf00043417

28. Lemaitre J, Chaboche J-L (1990) Mechanics of solid materials. Cambridge University Press. https://doi.org/10.1017/CBO9781139167970

29. Spencer AJM (ed) (1984) Continuum theory of the mechanics of fibre-reinforced composites. Springer, Vienna. https://doi.org/10.1007/978-3-7091-4336-0

30. Germain P, Nguyen QS, Suquet P (1983) Continuum thermodynamics. J Appl Mech 50(4b):1010–1020. https://doi.org/10.1115/1.3167184

31. Francfort GA (2018) Recovering convexity in non-associated plasticity. C R Méc 346(3):198–205. https://doi.org/10.1016/j.crme.2017.12.005

32. Lucchetta A, Auslender F, Bornert M, Kondo D (2021) Incremental variational homogenization of elastoplastic composites with isotropic and Armstrong–Frederick type nonlinear kinematic hardening. Int J Solids Struct 222:223:111000. https://doi.org/10.1016/j.ijsolstr.2021.02.011

33. Xia QS, Boyce MC, Parks DM (2002) A constitutive model for the anisotropic elastic-plastic deformation of paper and paperboard. Int J Solids Struct 39(15):4053–4071. https://doi.org/10.1016/S0020-7683(02)00238-X

34. Li Y, Stapleton SE, Reese S, Simon J-W (2018) Anisotropic elastic-plastic deformation of paper: out-of-plane model. Int J Solids Struct 130–131:172–182. https://doi.org/10.1016/j.ijsolstr.2017.10.003

35. Armstrong PJ, Frederick CO (1966) A mathematical representation of the multiaxial Bauschinger effect. CEGB: Report RD/B/N 731

36. Chaboche J-L (1993) Cyclic viscoplastic constitutive equations, part I: a thermodynamically consistent formulation. J Appl Mech 60(4):813–821. https://doi.org/10.1115/1.2900988

37. Ilyani Akmar AB, Lahmer T, Bordas SPA, Beex LAA, Rabczuk T (2014) Uncertainty quantification of dry woven fabrics: a sensitivity analysis on material properties. Compos Struct 116:1–17. https://doi.org/10.1016/j.compstruct.2014.04.014

38. Johnsen J, Clausen AH, Grytten F, Benallal A, Hopperstad OS (2019) A thermo-elasto-viscoplastic constitutive model for polymers. J Mech Phys Solids 124:681–701. https://doi.org/10.1016/j.jmps.2018.11.018

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.