Elasto‑viscoplastic model for rayon yarns

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.


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 Vol:.(1234567890) elements, in architectural membranes, and as preforms 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 important ones being cotton, rayon, polyamides, polyester (PET), aramids and polyethylene naphthalate (PEN).These materials are selected depending on the application envisaged.Rayon is a good choice for passenger 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 support the loads and contain the air pressure, providing the axial and lateral rigidity for acceleration, braking and cornering [5].The accuracy in predicting the performance 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 characterizing 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 fibers 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 behavior 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 textile 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 interactions 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 fiberreinforced rubber structures rely mainly on numerical methods.Structural reinforcements are generally modeled using two approaches.In the first one the entire structure is considered as a continuum and the modeling is based on the composite materials theory, see e.g.[13].In the second one the rubber matrix and the yarns/cords are considered separately and the reinforcements are modeled using the so-called "rebar" elements.The basic idea is to insert special bar elements with their own material properties into volume elements that should be strengthened [14][15][16].This latter technique is especially used for tires.
To feed the macroscopic models with a feasible 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 trajectories, 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 developed 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 anisotropic elastoplastic constitutive model for the flax fibers is considered; the reinforcement in this case 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 viscoplastic 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 threedimensional continuum of cylindrical shape, made of a material modeled by the non-isotropic constitutive law here developed.As cords are made up of singleply 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).

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 different 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 sample 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 conditions.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 −1 to maintain approximately the same Moisture Content (MC) during 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 −1 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 −1 , 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 testing 20 nominally identical specimens and the dashed curves represent the averages.Although the actual experimental results are obtained in terms of measured 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 2 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 −1 , 0.01 s −1 , 0.003 s −1 , 0.0017 s −1 and 0.0003 s −1 .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 −1 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.

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 1 , x 2 , x 3 } = Cartesian coordi- nates, with x 3 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: The unit tangent vector at each point a(x) hence reads . 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 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 2 ) for all the yarns considered.

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 ratedependent, 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 unloadingreloading 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 ( Ẽ,η ) and with a Bingham viscoplastic unit comprising three elements in parallel: a frictional device ( y ), a possibly nonlinear hardening spring (H) and a viscous dash-pot ( ).Based on this rheological model, the total strain is hence the sum of three contributions: elastic, viscoelastic and viscoplastic: The generalization of the 1D model to a three-dimensional continuum context is dealt with in Sect.

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.

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 e , viscoelastic ve and viscoplastic vp strains: As the response of the material depends on the fibers orientation, we assume that the Helmholtz freeenergy (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 e and i .Adopting the classic formalism of thermodynamics with internal variables, is thus a function of state and can be written as follows: where V k denotes a set of internal (scalar, vector or ten- sor) variables, which in the present model describe the material hardening.Specifically, in our case, we choose to define two internal variables and n , associated respectively to a linear and a nonlinear kinematic hardening.The observable variables are the total strain and the fiber direction a , while ve and 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 e is related to the spring E in the 1D rheological model, while the inelastic contribution i comes from the deformation of the viscolelastic Kelvin element, with the spring Ẽ in parallel to the dash-pot η , and of the viscoplastic Bingham element, composed of a slider with yield limit y in parallel with a spring H and a dash-pot .The inelastic term is thus further split into two contributions: where ve is the viscoelastic contribution, while the term vp represents the energy locked in the material due to microscopic rearrangements related to hardening.The current state of vp is assumed to be independent of the material orientation a. (5) As detailed in [29], since the sense of a is immate- rial, both the elastic and viscoelastic contributions of , respectively e and ve , must be even functions of a .Accordingly, they can be expressed as functions of the second order tensor a ⊗ a .Moreover, as the free energy must be unmodified if the material and its fibers undergo a rotation, then both e and ve must be scalar-valued isotropic functions of their arguments and can thus be written in terms of the five independent invariants respectively of e and a , and of ve and a .Considering a generic second order symmetric ten- sor S , representing either the elastic strain e or the viscoelastic strain ve , the five invariants are given by: where tr(•) is the trace of •.
To obtain linear stress-strain relations, both e and 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 e ( being the mass density), the most general quadratic form in e reads and depends on five material parameters: L , P are the shear moduli respectively along the fibers direction a (index L) and in the plane of isotropy (orthog- onal to a , index P), and , , are three material constants that can be derived from the more usual engineering elastic constants E L , E P , P , and PL , as shown in "Appendix".These 5 constants come from the generalization to 3D of the spring E. The isotropic contribution1 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: As relation (12) must hold for any process and the strain rate ̇ can be chosen arbitrarily, we have: The first condition defines the state law for the stress tensor , that is given by: with C being the elastic stiffness of the material: 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: The second condition in (13) expresses the requirement of non-negativeness of the dissipation and defines the thermodynamic forces associated to the viscoelastic strain ve and to the internal variables and n .Specifically, one has: with the latter two representing the current backstresses.From the first of relations (17), accounting for relation (9), one has: with Finally, from the last two relations in (17), using relations (10) one finds 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 nonnegative convex and zero at the origin dissipation potential Φ , defined as follows: ve ∶= e e − ve ve , X ∶= vp , ve In particular, assuming a linear viscoelastic material response, the dissipation potential Φ ve must be a posi- tive definite quadratic form of ̇ ve , while the depend- ence on the fiber direction is introduced again through the invariants defined by relations (7), such that: with L , P , , , being characteristic retardation times, and μL , μP , λ, α, β being the five material con- stants used for the viscoelastic free-energy (9).
The derivation of Φ ve gives the thermodynamic force ve defined in (17) and associated with ve , such that: with By combining relations (17), ( 18), ( 13) and ( 23), one has: The viscoplastic contribution to the dissipation process is more easily described in terms of the dual potential Φ vp * , obtained through the Legendre- Fenchel transform of the dissipation potential Φ vp .Accordingly, one has: To consider an elastic domain, the dissipation potential must depend on the positive part ⟨ ⟩ of the acti- vation (or yield) function .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 and the internal variables n in the expression of the potential (see [30][31][32]), considering it as a parameter. (26) Vol:. (1234567890) 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): where B is a positive scalar material parameter.Note that the term added to 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 where the equivalent stress eq represents the modulus of the projection of the over-stress ( − X − X n ) along the fibers, with 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: where sign(•) gives the sign of • .The term ⟨ ⟩ 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 Φ vp * defined by (27) and read ( 27) 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 con- sidered to be active only in the material direction.This is due to the particular choice for the yield function .

Parameters identification
The presented model has a total of 20 material parameters.These are 5 constants L , P , , , (or L , P , E L , E P , PL ) governing the instantaneous elas- tic behavior, 10 constants μL , μP , λ, α, β and L , P , , , describing the visco-elastic response of the material, 1 initial yield stress y , 1 visco-plastic parameter , and the 3 scalars H , H n and B describ- ing 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 μL , μP , λ, α, β be proportional to the corresponding L , P , , , through a single scalar parameter C which can be interpreted as the ratio of the spring stiffness Ẽ E of the 1D model.Namely: Similarly for the viscous effect, a single characteristic retardation time , which can be interpreted as the ratio η∕ Ẽ , can be postulated: With the above assumptions the independent viscoelastic 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 2 as shown in Fig. 3b.The initial elastic slope allows to fix a condition for E L and ẼL ; the elastic limit at low strain rate gives y ; the slope of the hardening phase depends on H , while H n modifies the initial nonlinear inelastic phase; the difference between the almost parallel hardening curves obtained at different strain rates gives the viscoplastic coefficient .
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 , P , LP and 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 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.

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 Δ , and the values , ve , vp , , n at time t n , the algorithm computes the updated quantities at time t n+1 .Rela- tions ( 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: 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 and ve are updated at time t n+1 , while vp , , n are kept constant from time t n .Otherwise, a return mapping step is applied and all the values , ve , vp , , n 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 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.

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 −1 .
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 ε = 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).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 eq , measured up to the total axial strain ε and defined as the integral in time of || ̇ 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 vp (as can be checked by using relations ( 20) and ( 30)).During a tensile 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 −1 , 0.01 s −1 , 0.003 s −1 , 0.0017 s −1 and 0.0003 s −1 ; b yarns Y0H and Y48H at strain rates of 0.017 s −1 : comparison between experimental data (shaded regions) and numerical results (continuous curves) Table 2 Identified material constants (humid conditions) 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 −1 , 0.003 s −1 , and 0.0003 s −1 .
The material model is able to predict the strain rate dependence that is particularly visible in the elasticviscoplastic 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 : 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 (34) p(S r ) = p d (1 − S r ) + p h S r a significant decrease of the initial modulus.This is captured fairly well by the geometrical model used for the numerical simulations.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 the twist of ± 3% .Figure 14 shows that the result- ing 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 ε = 0.13 .A more detailed exami- nation of such stresses at the Gauss points nearest to the yarn axis shows a maximum relative difference of 1.3% for three adopted discretizations.

Conclusions
In this paper, a new anisotropic viscoelastic-viscoplastic 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 Fig. 13 Tensile response of yarns with different twist levels at a strain rate of 0.017 s −1 in dry conditions.Comparison of the experimental (dashed) and numerical (continuous) responses Table 3 Identified material constants (dry conditions) 1.67 80.9 4500 73839 492.27 3000 Fig. 14 Numerical response of yarn Y48H under given monotonically increasing strain for ± 3% variation of twist level 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 elastoviscoplastic 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 nonuniform 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 , PL , P , and the coef- ficients , , , P of the invariant-based approach used to define the free-energy (5): with and In the above relations, E L is the Young's modulus along the direction of reinforcement, LP is the Poisson's coefficient governing the deformation in the plane of isotropy due to a stress along the direction of reinforcement, 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 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 , PL , P , L and then with the coef- ficients , , , L , P , fixing a certain direction of reinforcement a .Since both stiffnesses describe the same transversely isotropic material, their components must be equal.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/.

Fig. 1
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:

Fig. 2
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 4.2,(3) = e + ve + vp

Fig. 5 a
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

Fig. 7
Fig. 7 Finite element discretization of yarn Y48; the red arrows represent the directions of the fibers

Fig. 10 Fig. 11
Fig. 10 Radial contraction during a uniaxial tensile test of yarn Y48H: a microtomography images; b comparison of radius versus axial strain from experimental tests (markers) and numerical result (continuous)

Fig. 15
Fig. 15 Numerical response for three different FE discretizations for yarn Y48H (polar mesh, orthogonal mesh and unstructured mesh) at a strain of ε = 0.13.(Color figure online)

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 (35) = E P ( P + PL LP ) den = E P ( LP (1 + P − PL ) − P ) den = 1 − 2 P − 2 PL LP − 2 PL LP P (37) PL ∕E P = LP ∕E L .

Table 1
Mean values and standard deviation of the force at different strain levels for untwisted yarns at 0.017 s −1 strain rate 4ig.4Cyclic test on untwisted yarns at 0.003s −1 : a strain history; b stress versus strain response