Fully integrated multi-scale modelling of damage and time-dependency in thermoplastic-based woven composites

In this work, a multi-scale model established from the concept of periodic homogenization is utilized to predict the cyclic and time-dependent response of thermoplastic-based woven composites. The macroscopic behaviour of the composite is determined from finite element simulations of the representative unit cell of the periodic microstructure, where the local non-linear constitutive laws of the components are directly integrated, namely, the matrix and the yarns. The thermoplastic matrix is described by a phenomenological multi-mechanisms constitutive model accounting for viscoelasticity, viscoplasticity and ductile damage. For the yarns, a hybrid micromechanical–phenomenological constitutive model accounting for anisotropic damage and anelasticity induced by the presence of a diffuse micro-crack network is utilized. The capabilities of the overall multi-scale model are validated by comparing the numerical predictions with experimental data. Further illustrative examples are also provided, where the composite undergoes time-dependent deformations under uni-axial and non-proportional multi-axial loading paths. The multi-scale model is also employed to analyze the influence of the local deformation processes on the macroscopic response of the composite.


Introduction
In the automotive industry, thermoplastic-based woven composites appear to be particularly attractive, due their lightweight and good mechanical properties associated with a relative ease of manufacturing. However, the wide-scale use of these composites in structural applications has been hampered by the difficulty to predict their complex behaviour, especially over long periods of time.
Although the overall behaviour of thermoplastic-based composites can be described in a purely phenomenological manner (Feld et al., 2018;Launay et al., 2011;Mandys et al., 2019;Vasiukov et al., 2015), such a modelling strategy generally leads to a large number of parameters that may be difficult to identify and their values are valid only for a unique microstructure configuration. This is a real limitation towards capturing accurately the local behaviour and deformation processes. To overcome this issue, the composites response can be described with the help of multi-scale modelling approaches, by considering a representative volume element (RVE) of the microstructure. The main advantage of this modelling strategy is that the local constitutive laws of the constituents are directly integrated along with the geometric definition of the microstructure, leading to a more physical description of the deformation processes.
The multi-scale modelling of heterogeneous materials, in which the components exhibit nonlinear mechanical response, requires advanced homogenization techniques. Among the existing modelling strategies, two of them appear to be particularly attractive for this purpose, as they can easily account for many types of constitutive laws and microstructures: • The first one considers the so-called mean-field approaches (Mura, 1987;Nemat-Nasser and Hori, 1999;Qu and Cherkaoui, 2006) based on the Eshelby's solution (Eshelby, 1957) such as the micromechanical scheme of Mori-Tanaka (Benveniste, 1987;Mori and Tanaka, 1973), the self-consistent method (Hill, 1965) or the method proposed by Ponte Castañeda and Willis (1995). These theories usually deal with randomly distributed ellipsoidal inhomogeneities embedded in a matrix material. Therefore, they are particularly adapted to the cases of short and long fibres reinforced composites. Woven composites may be eventually considered if the yarns waviness is neglected, or by discretizing the wavy yarns into smaller but straight segments and replacing them with an equivalent system of ellipsoids (Abdin et al., 2016). While initially developed for linear elasticity, these methods can be extended to the case of non-linear behaviour through an incremental formulation (Gavazzi and Lagoudas, 1990;Lagoudas et al., 1991). Although these methods are based on semi-analytical solutions with fast computations, they unfortunately keep a limited accuracy and must be used with caution, especially when the matrix phase exhibits a non-linear response (Lagoudas et al., 1991). This inaccuracy is mainly caused by the local fields in the matrix that are considered in an average way. This is likely not to give good enough estimates when the non-linearity is matrix-dominated. More accurate solutions are sometimes obtained through specific enhancements that are mostly dedicated to short or long fibres reinforced composites (Barral et al., 2020;Chaboche et al., 2005;Doghri and Ouaar, 2003) rather than the case of woven composites. • The second strategy deals with full-field approaches, whose most representative case is the periodic homogenization (Michel et al., 1999;Sanchez-Palencia, 1974;Suquet, 1987). By definition, this method is particularly suitable for composites with periodic microstructure, such as woven composites (Courtois et al., 2019;Hofer et al., 2020;Ullah et al., 2019;Xu et al., 2018). Moreover, periodic homogenization allows to properly define the concept of homogenized behaviour through a rigorous framework, which is fundamentally independent of the choice of the local constitutive laws of the different phases Chatzigeorgiou et al., 2015). Due to its periodicity, the RVE is represented by the smallest repeating element, commonly referred to as the unit cell, on which well-defined periodic boundary conditions are applied (Li, 1999(Li, , 2001Li and Wongsto, 2004). Except for particular cases, the resolution of such a problem requires the use of numerical techniques, mostly the finite element (FE) method.
The geometry of the unit cell is then represented through a FE mesh integrating the local constitutive equations of each sub-domain. Compared to mean-field approaches, the use of periodic homogenization may lead to an important computational cost, which is nevertheless necessary to accurately and properly describe the behaviour of certain types of composite materials such as woven composites.
In the present work, a fully integrated multi-scale model is utilized to describe and study the cyclic and time-dependent response of thermoplastic-based woven composites, under the small strain assumption and isothermal conditions. Based on the framework of periodic homogenization (Michel et al., 1999;Sanchez-Palencia, 1974;Suquet, 1987), a proper transition is established between the macroscopic and microscopic scales, enabling a direct integration of the local nonlinear constitutive laws of the components, namely, the thermoplastic polymer matrix and the yarns. The proposed multi-scale model is thus referred as 'fully integrated' as the macroscopic response of the composite is determined from FE simulations on its unit cell. The local behaviour of the thermoplastic matrix is described by a phenomenological multi-mechanisms constitutive model accounting for viscoelasticity, viscoplasticity and ductile damage (Praud, 2018;Praud et al., 2017a). For the yarns, a hybrid micromechanical-phenomenological constitutive model is utilized to represent the anisotropic damage and anelasticity induced by the presence of a diffuse micro-crack network (Praud, 2018;Praud et al., 2017b).
By considering the theory of periodic homogenization and the established constitutive models for the matrix and the yarn phases, the novelty of this work lies within the development of a comprehensive and robust fully integrated multi-scale model to describe the cyclic and time-dependent response of thermoplastic-based woven composites. The accuracy of the proposed multi-scale model is evaluated by comparing the numerical predictions with experiments carried out on a woven composite made of a polyamide 6-6 matrix with an E-glass balanced 2-2 twill weave reinforcement, referred to as the studied composite. The multi-scale model is furthermore employed to understand the local damage mechanisms and deformation processes, and to analyze their influence on the macroscopic response of the composite in relation with its microstructure. Up to the author's knowledge, such a comprehensive multi-scale analysis, integrating highly complicated constitutive laws for the matrix and the yarns, has never been performed before for this type of materials.
In this paper, the following notation is adopted. Bold lower-case Latin symbols are vectors (a, b. . .). Bold upper-case Latin and Greek symbols are second-order tensors (A, B. . . a; b. . .) and I denotes the second-order tensors identity tensor. Blackboard symbols are fourth-order tensors (A; B . . .), and I denotes the fourth-order tensors identity tensor. All the others standard symbols (neither in bold nor in blackboard) are scalars (a, b . . . A, B . . . a, b . . .). Furthermore, Á; : and represent the single contracted, twice contracted and dyadic products, respectively.

Theoretical background
A periodic medium is defined by a unit cell (the smallest repeating unit element) that is representative of the microstructure. If the dimensions of the unit cell are small enough compared to the dimensions of the macroscopic medium, then a scale separation can be assumed with two scales: a macroscopic one, defined by the macroscopic (or global) coordinates, denoted by x 2 V, and a microscopic one, defined by the microscopic (or local) coordinates, denoted by x 2 V. The assumption of scale separation ( V >> V) allows to consider that the unit cell is representative of a macroscopic material point. At the macroscopic scale, the heterogeneous medium is replaced by an equivalent homogenized medium ( Figure 1). Under quasi-static conditions (the inertia effects are considered negligible), the motion of any macroscopic and microscopic material point Mð xÞ and Mð x; xÞ, respectively, is governed by the macro-scale and micro-scale equations given in Table 1. The homogenization consists in defining the global constitutive behaviour, described by the operator F (see the left part of Table 1), which is implicitly defined by the local operator F (see the right part of Table 1) and the geometrical characteristics of the various constituents in the unit cell. The local models may be non-linear and involve internal state variables. They may also be timedependent (viscoelasticity, viscoplasticity, etc.).
In order to link the micro-scale with the macro-scale equations, a connection between scales is required. To do so, the macroscopic stress and strain fields are identified through volume average of their microscopic counterparts over the unit cell (Hill, 1967;Mura, 1987;Nemat-Nasser and Hori, 1999). Moreover, from the divergence theorem, it can be demonstrated that the stress and strain averages within a unit cell are also connected to the traction and displacement vectors applied on its boundaries rð x; x; tÞ Á nðxÞ x dS (1) Figure 1. Schematic representation of a heterogeneous material with a periodic microstructure by considering a scale separation.  where n stands for the outward normal defined on each point of the unit cell borders (8x 2 @V).
In Table 1, the macro-scale equations are supplemented by boundary conditions specific to the treated problem. For the micro-scale equations, the boundary conditions arise from the assumption of periodicity (Suquet, 1987) implying that, within the unit cell, the displacement vector u of any microscopic material point Mð x; xÞ can be written under the following form where the first term stands for the affine part of the local displacement field. The latter is directly related to the macroscopic strain e. The second term u 0 represents a periodic fluctuation, while the last term u 0 depicts a rigid body motion that comes out from the macroscopic problem and consequently does not depend on the microscopic problem. As u 0 is periodic, it takes the same value on each pair of opposite points x 1 and x À of the unit cell boundaries (x 1 ; x À 2 @V) By substituting (3) into (4), these periodicity conditions can be expressed in terms of u instead of u 0 , while involving the macroscopic strain tensor e uð x; x 1 ; tÞ À uð x; x À ; tÞ ¼ eð x; tÞ Á ðx þ À x À Þ Note that, due to its periodic aspect, the strain average (within the unit cell) produced by u 0 is null. Therefore, the total strain average is well equal to the macroscopic strain as expressed in (2).

Finite element resolution of the unit cell problem
The relationship between the macroscopic stress and strain (operator F) results from the solution of the microscopic problem (given in Table 1). The latter can be easily solved by the FE method, using a mesh of the unit cell, on which the periodicity conditions (5) are applied using the generalized methodology proposed by Li (1999Li ( , 2001, Li and Wongsto (2004) and Praud (2018), which is described in this section.
Considering a cubic FE unit cell, 1 the aim of this method is to apply a macroscopic strain e, by taking into account the periodic boundary conditions previously described in (5). A displacement gradient is then applied between each pair of opposite boundary nodes, respectively, denoted by the indices i and j. This gradient is directly related to the macroscopic strain tensor e sym: In the above expression, u 1 , u 2 , u 3 and x 1 , x 2 , x 3 are the components of the displacement and position vectors, u and x, respectively. The proposed method introduces the six components of the macroscopic strain tensor as 'additional degrees of freedom' that are directly involved in the boundary conditions (Li, 1999(Li, , 2001Li and Wongsto, 2004). These 'additional degrees of freedom' are linked to the mesh of the unit cell using kinematic constraint equations obtained from (6) and are thus used as 'constraint drivers'. From a practical point of view, the 'constraint drivers' can be introduced by adding six new nodes among which only the first degree of freedom is used. Each displacement of those 'constraint drivers', noted u cd 11 ; u cd 22 ; u cd 33 ; u cd 12 ; u cd 13 and u cd 23 , takes the values of each component of the macroscopic strain tensor: e 11 ; e 22 ; e 33 ; 2 e 12 ; 2 e 13 and 2 e 23 , respectively. By proceeding this way, after properly formulating the kinematic constraint equations (Praud, 2018), the dual forces associated with the 'constraint drivers', noted F cd 11 ; F cd 22 ; F cd 33 ; F cd 12 ; F cd 13 and F cd 23 , directly provide the product of the unit cell's volume with the corresponding components of the macroscopic stress tensor, r 11 ; r 22 ; r 33 ; r 12 ; r 13 and r 23 , respectively, as illustrated in Figure 2. This approach enables to apply any state of macroscopic strain or stress on the unit cell. Since each 'constraint driver' can be independently set either in terms of displacement (macroscopic strain) or in terms of force (macroscopic stress). Mixed stress-strain states can also be assigned. As previously mentioned, the solution of a periodic homogenization is determined up to a rigid body motion. For this reason, a node of the model needs to be clamped in order to guarantee the uniqueness of the solution and the solvability of the FE system.

Unit cells for woven composites
In order to apply the framework of periodic homogenization to the case of woven composites, it is necessary to provide a geometric representation of the unit cell along with its associated FE mesh. With this purpose, many works have been already undertaken (Cou egnat, 2008;Lin et al., 2011;Lomov et al., 2000Lomov et al., , 2001Lomov et al., , 2007Sherburn, 2007), leading to the development of several dedicated tools. Among these tools, the TexGen platform (Lin et al., 2011;Sherburn, 2007), which is utilized in this work, is an open source and free software developed at the University of Nottingham and shows very interesting capabilities in generating FE unit cell of any woven pattern.
Concerning the 2-2 twill weave pattern studied in the present case (see Figure 3(b)), the geometry of the woven microstructure is defined through the parameters l 1 , l 2 , l 3 , l 4 and l 5 , shown in Figure 3(a). For the examined composite, the characteristic dimensions l 1 , l 2 and l 3 are evaluated from experimental observations carried out by means of X-ray computed micro-tomography, as illustrated in Figure 4. The dimension l 4 is obtained by dividing the total thickness of the provided plates by the number of layers they contain. The parameter l 5 has been deliberately over-sized in the present work to ensure a proper enough mesh within the matrix domain. The average values considered for the multi-scale model are listed in Table 2. Note that, under this configuration, the unit cell yields to a yarn volume fraction of 59%. Furthermore, analyses from micrography have shown that the intrayarn fibre volume fraction is about 85%. Overall, the fibre volume fraction in the whole composite is about 50%, which is well in accordance with the supplier data.  Using the data in Table 2, the unit cell of the studied composite is generated and subsequently meshed through the TexGen utilities (Lin et al., 2011;Sherburn, 2007) (see Figure 5). The yarns, which are composed of numerous fibres unidirectionally oriented inside the matrix, are treated within the described modelling approach as an equivalent homogeneous medium with an anisotropic constitutive behaviour. The material orientation is defined for each point with respect to the yarns middle line whose waviness is automatically calculated according to the considered weaving pattern, as shown in Figure 5(d).
The mesh of the unit cell, shown in Figure 5, contains 8,522 nodes and 40,060 first-order tetrahedral elements (C3D4 in ABAQUS/Standard). Note that a spatial convergence analysis was carried out regarding the mesh refinement in order to obtain a good enough compromise between accuracy and reasonable computational cost.

Local constitutive models
The constitutive models (or sub-models) presented in this section are utilized to describe the behaviour of the matrix and yarn phases. Nevertheless, it should be pointed out that the homogenization framework previously introduced is completely independent of the local constitutive models, which Table 2. Characteristic dimensions of the microstructure and adopted values for the studied composite. means that the sub-models presented in this section work as 'bricks' that can be easily interchanged with other sub-models.

Local constitutive model for the matrix
Thermoplastics and more especially semi-crystalline polymers, such as polyamide 6-6, are known to exhibit a complex dissipative behaviour combining both fluid and solid properties. Moreover, experimental investigations (Detrez et al., 2011) showed that these materials, when loaded, present a reduction of their apparent stiffness induced by damage mechanisms. To account for this particular behaviour, the constitutive model developed for such type of materials by Praud (2018) and Praud et al. (2017a) is utilized to describe the response of the matrix phase. This model is written within the framework of thermodynamics Germain et al., 1983;Lemaitre and Chaboche, 1990) and integrates viscoelasticity, viscoplasticity and ductile damage through a multi-mechanisms formulation, as represented by the rheological scheme in Figure 6. For the matrix phase, the state laws are based on the following form of the Helmholtz free energy potential which describes the energetic state of any matrix material point as a function of the strain e, as observable state variable, and the internal state variables e v i ; e p , r and D that represent the strain related to the ith Kelvin-Voigt branch, the viscoplastic strain, the hardening variable and the damage variable, respectively. Through the thermodynamical derivation, the stress in the matrix phase is expressed as In the same manner, for the other associated variables related to the matrix phase, it yields and where, r v i , R(r) and Y denote the viscous stress in the ith Kelvin-Voigt branch, the hardening function and the energy density release rate, respectively. Note that, in equations (7) to (9), C e and C v i denote the initial fourth-order stiffness tensors of the single spring and the spring of the ith Kelvin-Voigt branch, respectively. These tensors are classically formulated for bulk isotropic materials and are defined by the Young modulus E e or E v i , respectively, as well as the Poisson ratio which is assumed to be the same in all stiffness tensors. Moreover, in equations (7) and (10), the hardening function R(r) is chosen as a power law, such as R r ð Þ ¼ Kr n , where K and n are material parameters.
The evolution of each viscoelastic strain e v i is governed by a dual dissipation potential written as whose derivative with respect to the associated variable r v i gives the following evolution law where V v i denotes the viscosity tensor of the linear dash-pot of the ith Kelvin-Voigt branch. As for the stiffness tensors, each viscosity tensor is classically formulated for bulk isotropic materials and is defined by a viscosity g v i and the Poisson ratio . The latter is assumed to be the same in each stiffness and viscosity tensor.
In the matrix phase, the damage mechanism is assumed to be ductile as it only evolves with the viscoplasticity, making those two mechanisms directly coupled and simultaneously activated. Such an evolution is represented through the normality rule of an indicative function (Krairi and Doghri, 2014;Lemaitre, 1985) given by In the above equation, the first term f r; R; D ð Þrepresents the yield function where eq r ð Þ designates the equivalent Von Mises stress 2 and R 0 the yield threshold. The second term f D Y; D ð Þdepicts the damage contribution in which S and b are material parameters. Thus, expressing the normality of the indicative function F r; R; Y; D ð Þleads to the following viscoplastic-damage flow rules where k depicts the viscoplastic-damage multiplier, which appears to be equal to the hardening variable (k ¼ r). Its evolution is obtained by considering that the positive part of the yield function f is linked to _ r by In this relationship, Q _ r ð Þ represents a viscous stress function, which is chosen under the form of a power law, such as Q _ r ð Þ ¼ H _ r m , where H and m are material parameters. It should be noticed that r is only allowed to increase, _ r > 0 if f > 0, or to remain constant, _ r ¼ 0 if f < 0. Note that the convexity of u Ã and F, given in equations (11) and (13), respectively, guarantees the thermodynamical admissibility of the evolution laws (12), (14) to (16).
The constitutive models of the matrix have been implemented into the FE solver ABAQUS/ Standard by means of User MATerial subroutines (UMAT). The numerical integration scheme is based on the 'convex cutting plane' form of the 'return mapping algorithm' (Ortiz and Simo, 1986;Simo and Hughes, 1998;Simo and Ortiz, 1985), which is described in Praud (2018) and Praud et al. (2017) for this model. Readers are referred to these references for more details regarding the numerical implementation.
Local constitutive model for the yarns As previously mentioned, at lower scale a yarn is composed of numerous unidirectionally oriented fibres embedded in the matrix. Consequently, their behaviour can be considered as equivalent to the one of a unidirectional composite. In media with continuous reinforcement (long fibres), the longitudinal behaviour (tension 11) appears to be linear elastic up to the material brittle failure due to fibres breakage. The transverse behaviour (tension 22 and/or in-plane shear 12) exhibits a more progressive degradation induced by the growth of a diffuse micro-crack network that initiates by debonding at the fibre/matrix interfaces and propagates by coalescence (Yang et al., 2020). To represent this behaviour, the constitutive model specifically developed for this type of materials (Praud, 2018;Praud et al., 2017b) is utilized to describe the response of the yarn phase. This model is also written within the framework of thermodynamics Germain et al., 1983;Lemaitre and Chaboche, 1990) while being based on a hybrid micromechanical-phenomenological formulation. Accordingly, the damage is introduced through a micromechanical description of a RVE containing a micro-crack density c c , as illustrated in Figure 7. Furthermore, in this type of material, damage is often accompanied by permanent strains caused by the non-closure effect of the micro-cracks. This aspect is phenomenologically described by an anelastic strain tensor e s that is added to the elastic part of the strain.
For the yarn phase, the state laws are based on the following form of the Helmholtz free energy potential qw e; e s ; c c ð Þ ¼ 1 2 e À e s ð Þ : ½C 0 À D c c ð Þ : e À e s ð Þ This potential depends on the strain e, as observable state variable, and the internal state variables e s and c c representing the anelastic strain and the micro-crack density, respectively.
Through the thermodynamical derivation, the stress in the yarn phase is expressed as The other associated variables related to the yarn phase are given by where, Y c depicts the energy density release rate.
In equations (18) to (20), the initial stiffness tensor C 0 , defined as transversely isotropic, is gradually lowered by the stiffness reduction tensor, denoted by D c c ð Þ . The latter is calculated as a function of the micro-crack density c c through the micromechanical scheme of Mori-Tanaka (Mori and Tanaka, 1973), where micro-cracks are idealized by quasi-flat ellipsoidal inclusions of void (with zero stiffness) whose normal is oriented along the second directionx 2 (see Figure 7). Indeed, due to the microstructure arrangement in the yarns, the micro-cracks are forced to propagate in a plane parallel to the fibre directionx 1 (Li et al., 2019). Moreover, by assuming that the yarns are mainly loaded in the plane (x 1 ;x 2 ), the propagation plane of the micro-cracks necessarily becomes perpendicular to the second directionx 2 and consequently always oriented in the same plane. Accordingly, it yields for the stiffness reduction tensor where T c is the interaction tensor of the void inclusion, itself defined from the Eshelby's tensor S E (Eshelby, 1957). The latter is numerically evaluated (Gavazzi and Lagoudas, 1990) from the stiffness tensor of the reference medium, namely C 0 , and the geometrical configuration of the void inclusion, which is the same as in Praud (2018) and Praud et al. (2017b). A 0 c c ð Þ designates the strain localization tensor, which allows to define the local strain and stress fields, and more specially the ones in the virgin part of the material e 0 e; c c ; e s ð Þ ¼ A 0 c c ð Þ : e À e s ð Þ ; r 0 e; c c ; e s ð Þ ¼ C 0 : A 0 c c ð Þ : e À e s ð Þ In the yarns, the damage is driven by Hill-like quadratic interaction criterion (Hill, 1948;Ottosen and Ristinmaa, 2005) expressed from the stress in the virgin part of the material, namely r 0 where H is a fourth-order tensor set in such a way that H c only depends on the components 22 and 12 of r 0 . H contains the parameters R 22 and R 12 , which represent the initial damage thresholds in pure transverse tension and in pure in-plane shear, respectively. The development of the micro-crack density c c initiates only when H c exceeds 1. After this stage, it follows a Weibull-like function g of the maximal value of H c reached in the whole loading history, noted hereafter as sup H c ð Þ. This gives where S and b are material parameters. Additionally, c 1 c depicts a saturation level that cannot be exceeded (c c < c 1 c ). Note that the development relationship (25) accounts well for the irreversible aspect of damage. Indeed, c c can only grow, The proposed model for the yarn phase considers that the anelasticity occurs as permanent strains caused by the micro-cracks non-closure upon unloading. Thus, the anelastic strain e s evolves with the micro-crack density c c , making the two mechanisms directly coupled and simultaneously activated. This evolution is governed by the normality rule of an indicative function given by with In the above two equations, H s r ð Þ is another quadratic interaction criterion, where the fourthorder tensor F is set in such a way that H s only depends on the components 22 and 12 of r. F contains the parameters a 22 and a 12 , which are transverse tension and in-plane shear anelasticity parameters, respectively. Therefore, applying the normality of the indicative function F r; Y c ð Þgives the following evolution laws where k is the damage-anelasticity multiplier, which appears to be equal to the micro-crack density c c (k ¼ c c ) whose evolution is governed by the damage development relationship (25). Note that the form of H s r ð Þ implies that only the components 22 and 12 of e s are active. Note that the convexity of F given in equation (27), as well as the irreversibility aspect of c c through the damage development relationship (25), ensures the thermodynamical admissibility of the evolution laws (28) and (29).
Like for the matrix, the constitutive model of the yarns has also been implemented into the FE solver ABAQUS/Standard by means of UMAT. The numerical integration scheme is based on the 'convex cutting plane' form of the 'return mapping algorithm' (Ortiz and Simo, 1986;Simo and Hughes, 1998;Simo and Ortiz, 1985), which is described in Praud (2018) and Praud et al. (2017b) for this model. Readers are referred to these references for more details regarding the numerical implementation.

Remarks and discussion
In the matrix phase, damage is assumed to be isotropic through the scalar internal state variable D. However, it is worth mentioning that, even if in reality the matrix is assumed to be initially isotropic, its response may exhibit anisotropic effects induced by damage depending on the directions of the principal stresses. Although some modelling approaches can take such behaviour into account (Ayadi et al., 2018;Chaboche, 1992;Ekh et al., 2003;Lemaitre and Desmorat, 2005), this requires unconventional characterization highlighting the directional effects of damage. In the absence of such data, the isotropic idealization of damage remains the most appropriate assumption.
Unlike the matrix, the definition of damage in the yarn phase is anisotropic through the stiffness reduction tensor D c c ð Þ . Based on micromechanical relationships, this tensor depends on the internal state variable c c , which represents a density of micro-cracks whose normal is oriented alongx 2 (see Figure 7). It is recalled that this assumes that the yarns are mainly loaded in the plane (x 1 ;x 2 ).
The constitutive model of the matrix phase is time-dependent whereas the one of the yarn phase is not. Since the yarns are composed of fibres and matrix, they should normally exhibit a timedependent response inherited from the matrix they contain. However, in the studied woven composite, the amount of matrix within the yarns is relatively low (about 15%). Therefore, it is assumed that the behaviour of the yarns is idealized as time-independent as it is mostly driven by the fibres. This assumption also allows to reduce the number of parameters for the identification, which is the topic of the next section.
It is important to mention that the constitutive models utilized for the matrix and the yarn phases make use of local continuum damage theories for modelling the degradation mechanisms caused by micro-defects within the unit cell of the woven composite. Thus, the damage mechanisms are described according to continuum quantities (internal state variables). This allows to conveniently formulate the problem according to the theory of continuum mechanics. This works well as long as the materials do not reach the strain softening regime, which was well verified in all the examples considered in this study. In the contrary case, the use of non-local modelling theories (Kiefer et al., 2018;Miehe et al., 2010;Peerlings et al., 1996;Pijaudier-Cabot and Bazant, 1987;Seupel and Kuna, 2019) would be necessary to objectively represent the localization phenomena beyond the strain softening onset (Bazant, 1984;De Borst et al., 1993).
Moreover, the macro-scale and micro-scale equations are both expressed in the small strains formalism. In this paper, all the studied cases and the related results are performed and presented at moderate strain levels. However, it is worth noticing that the proposed multi-scale model cannot capture geometrical non-linearities such as micropolar or chiral effects (Kadic et al., 2019) that are likely to occur at high strain levels in a woven microstructure. To account for these phenomena, the use of higher-order homogenization strategies (Forest, 2002;Goda et al., 2013;Goda and Ganghoffer, 2016;Larsson and Diebels, 2007) would be necessary. This would also require extending the local constitutive models of the matrix and the yarns to a large deformation formalism.

Experimental procedure and identification strategy
The matrix is identified from experiments directly performed on the unfilled polyamide 6-6, from specific experiments and a dedicated identification procedure, which are detailed in Praud (2018) and Praud et al. (2017b). The obtained parameters for the matrix are listed in Table 3. Note that polyamide and polyamide-based composites are well known to be highly sensitive to the environmental conditions (Arif et al., 2014;Launay et al., 2013;Malpot et al., 2015), especially the relative humidity (RH) and the temperature (T). In this work, the following environmental conditions were considered: RH ¼ 50% and T ¼ 23 C (room temperature) for both the studied composite and its unfilled polyamide 6-6 matrix.
Regarding the yarns, it is impossible to be separated from the matrix in order to perform mechanical tests on them. Consequently, the yarns parameters are identified from the macroscopic response of the entire woven composite through reverse techniques (Levenberg, 1944;Marquardt, 1963;Meraghni et al., 2011Meraghni et al., , 2014. Moreover, it is recalled that the behaviour of the yarns is assumed to be time-independent. Thus, for the multi-scale model, the overall time-dependency is driven exclusively by the rheology of the matrix phase. With this in mind, the parameters related to the yarns are identified only from strain-controlled off-axis tests. Those tests are performed on ½AEh s tensile laminated specimens 3 (see Figure 8) repeatedly loaded and unloaded at progressively increasing stress levels and at a relatively slow strain rate (about 3:5 Â 10 À3 s À 1 ) until failure. These tests, referred to as 'quasi-static', are performed with a servo-hydraulic tensile machine where the axial macroscopic strain e xx is measured by means of an extensometer, while the axial macroscopic stress r xx is monitored by a load cell. From these quasi-static tests, an optimization algorithm based on the Levenberg-Marquardt technique (Levenberg, 1944;Marquardt, 1963;Meraghni et al., 2011Meraghni et al., , 2014 is utilized to identify the parameters of the yarn phase. This is achieved by minimizing the least squares between the numerical and experimental macroscopic stress responses ( r xx ). Beforehand, the initially transversely Figure 8. Tensile laminated specimens ½AEh s axially loaded in thex direction. During the tests, the axial macroscopic strain e xx is measured by means of an extensometer, while the axial macroscopic stress r xx is monitored by a load cell. (a) Laminate: multi-layered composite; the thickness of a single layer is about 0.5 mm, making a total thickness of 2 mm for the whole laminated specimens. (b) Oriented layer: (x;ỹ;z) is the coordinate system of the laminate, while (x 1 ;x 2 ;x 3 ) is the material one within a single layer. (c) Geometry and adopted dimensions of the laminated specimens (dimensions in mm), rectangular specimens with large enough dimensions were considered to obtain experimental data that are reasonably representative of the macroscopic response of the composite.
isotropic stiffness of the yarns (C 0 ) is estimated by means of linear periodic homogenization (Praud, 2018). Therefore, only the parameters related to the damage and the anelasticity in the yarns are identified through the reverse engineering procedure. However, during the identification, the microcrack saturation parameter c 1 c is kept fixed at a sufficiently high level so that it is never reached and only the ascending part of the Weibull function (25) is acting (Praud, 2018;Praud et al., 2017b). In this manner, the identification remains well conditioned and robust. The obtained values of the yarns parameters are listed in Table 4.
Once all parameters are identified, the predictions obtained with the multi-scale model are compared with additional experimental data, where the laminated specimens are subjected to a sinusoidal stress signal oscillating at a frequency of 1 Hz during 100 s, i.e. 100 cycles. These tests are referred to as cyclic. Under cyclic loading, thermoplastic-based composites may exhibit a significant temperature increase due to the self-heating phenomenon, caused by the dissipation (Benaarbia et al., 2014(Benaarbia et al., , 2015Praud, 2018). For this reason, during the cyclic tests, thermal measurements were carried out by means of an IR thermal camera. Figure 9(a) to (d) shows the numerical simulations provided by the multi-scale model and the experimental data for the quasi-static tests performed with the ½0 4 ; ½AE15 s ; ½AE30 s and ½AE45 s specimens, respectively. It is recalled that these tests were used for the identification. First of all, it is observed that the response of the ½0 4 specimen is quasi-elastic, as only a small amount of apparent damage and anelasticity are generated before the composite brutally fails at a relatively high stress level (about 430 MPa with the experiments (see Figure 9(a) to (c)). For the ½AE45 s specimen, the multi-scale model overestimates the quasi-static response (see Figure 9(d)), retaining though the same order of stress levels and a similar tendency. The deviation of the present model may be caused by not accounting for other damage mechanisms occurring under significant in-plane shear stress levels, like for instance the debonding at the yarns/matrix interface. Figures 10 and 11 show the predictions provided by the multi-scale model and the experimental data obtained for the stress-controlled cyclic tests performed with the ½0 4 and ½AE45 s laminated specimens, respectively. These tests were not used for the identification. It is noted that the ½0 s specimen was loaded at a maximum stress level which corresponds to about 80% of the quasi-static failure point. For the ½AE45 s specimens, as the response is much more ductile, a maximum stress level of 95 MPa is taken instead, in order to keep a moderate strain amplitude. The failure occurred after 80 cycles for the ½0 4 specimen, while for the ½AE45 s no failure appeared before the completion of the 100 cycles. When the composite is subjected to a cyclic loading, phenomena directly inherited from the viscous and damage mechanisms of the matrix phase can be clearly observed, namely, the accumulation of strain accompanied by a progressive stiffness reduction. According to the layers angles, the amplitudes of these phenomena are quite different and, once more, bring out the anisotropy induced by the microstructure. Overall, similar tendencies are predicted by the multi-scale model, compared to the experimental results. The macroscopic strain responses are in relatively good agreement for the ½0 4 specimen, while the overall accumulation of strain and damage is underestimated for the ½AE45 s specimen.

Results and comparisons with experimental data
Moreover, it is worth pointing out that the hysteresis loops upon the loading/unloading stages of each single cycle are not well reproduced by the multi-scale model, especially for the ½AE45 s strain ( e xx ) for the ½AE45 s laminated specimen.
laminated specimen (see Figure 11). This is likely due to the sub-model of the yarn phase, which, unlike the one of the matrix, does not integrate any viscoelastic mechanism and time-dependent features in general. Moreover, from the thermal measurements (see Figure 12(a) and (b)), it can be observed that, while the self-heating is relatively limited for the ½0 4 specimen (about 3 C of temperature increase), this phenomenon becomes important for the ½AE45 s specimens (about 10 C of temperature increase). This can also explain the deviation between the experimental data and the model (see Figure 11), which assumes isothermal conditions in the present case.

Macroscopic response of the woven composite
In order to provide a better understanding of the multi-scale model, additional examples of simulations carried out on a single macroscopic material point, so-called virtual tests, are performed. In these examples, the parameters previously identified for the matrix and the yarns in Tables 3 and 4 are utilized. The composite is subjected to several loading configurations: warp tension, in-plane   shear and combined warp tension and in-plane shear stress states, as illustrated in Figure 13(a) to (c), respectively. These loading configurations are applied according to different paths including, creep and strain recovery in warp tension and in in-plane shear, and non-proportional combined warp tension and in-plane shear.

Creep and strain recovery
In the creep and strain recovery in warp tension 'virtual test', a normal stress r 11 of 350 MPa is first applied on the composite in 5 s. This stress is then held for 300 s before returning back to zero in 5 s.   In the final stage, the composite is kept at zero stress for another 300 s (see Figure 14(a)). The results of this 'virtual test' are presented in Figures 14 to 16. In the case of the creep and strain recovery in-plane shear, a shear stress r 12 of 40 MPa is applied on the composite with the same loading path as for the warp tension (see Figure 17(a)). The results of this 'virtual test' are presented in Figures 17 to 19. When the composite is loaded in the warp direction (0 s < t < 5 s), most of the load is carried by the warp yarns in their longitudinal direction. Meanwhile, a small part of this load is also carried by the weft yarns in their transverse direction. This brings about the occurrence of micro-cracks only in these yarns (see Figure 16(a)). In the matrix, most of the stresses are concentrated in the yarnscrossing areas, leading to the development of some damage in these locations (see Figure 15(a)). The degradations occurring in the weft yarns and in the matrix seem to have limited consequences at the macroscopic scale, as most of the load is carried by the warp yarns in their longitudinal direction, which behave elastically with an important stiffness. For this reason, the response of the composite at the macroscopic scale remains quasi-linear during the first loading stage (see Figure 14(c)).
When the composite is loaded in in-plane shear (0 s < t < 5 s), most of the load is carried by the yarns in in-plane shear. This brings about a fast growth of the micro-cracking in both warp and weft yarns (see Figure 19(a)). Then, the overall load is progressively transferred to the matrix, in the inter-yarns areas, where the stresses are mainly concentrated and where damage develops (see Figure 18(a)). Note that those degradation mechanisms have a significant influence on the macroscopic in-plane shear response of the composite, which appears to be more matrix-dominated and strongly non-linear (see Figure 17(c)).
Note that identical degradation mechanisms were experimentally observed by means of X-ray computed micro-tomography for a similar thermoplastic-based woven composite (Pomar ede et al., 2018).
During the creep stages, when the stress is held (5 s < t < 305 s), the macroscopic strain increases under the action of a constant macroscopic stress (see Figures 14(b) and (c), 17(b) and (c)). The macroscopic creep response of the composite is caused by the time-dependent behaviour of the matrix phase, as well as the microstructure interactions between the matrix and the yarn phases.
Although the sub-model of the yarns does not account for any time-dependent feature, a microcrack density growth is observed within the yarn phase (see Figures 16 and 19). Indeed, when the matrix locally creeps and gets damaged, a part of its sustained load is gradually transferred to the yarns, leading to an increase of the micro-crack density in this phase. Moreover, in warp tension, the macroscopic creep is quite low, in contrast to the in-plane shear. Indeed, in the latter case, the behaviour of the composite is mainly matrix-dominated.

Non-proportional combined warp tension and in-plane shear
In these examples, the composite is subjected to a combined warp tension and in-plane shear stress state, which is applied through two different loading paths having the same amplitudes. In the first loading path, referred to as path 1, a normal stress in the warp direction r 11 of 200 MPa is first applied in 5 s. In the next stage, r 11 is held constant, while an in-plane shear stress r 12 of 30 MPa is applied in 5 s. Afterwards, both r 11 and r 12 are held constant for 5 s before being released in another 5 s (see Figure 20(a)). In the second loading path, referred to as path 2, the normal and shear stresses are interchanged with respect to the path 1 (see Figure 20(b)).
The results of these simulations are presented in Figure 20. These 'virtual tests' highlight the importance of the loading path when the composite is subjected to combined stress states. Although the same stress amplitudes are applied for both paths 1 and 2, it can be noticed that the amplitudes of the strain responses are quite different (see Figure 20(c)). The maximum value of the normal strain e 11 resulting from the path 1 is slightly greater than the one obtained from the path 2. Indeed, the normal stress r 11 is held longer for the path 1 than for the path 2, which means that the composite is exposed for more time to creep in the warp direction. Similarly, the maximum value of the in-plane shear strain 2 e 12 resulting from the path 2 is far greater than the one obtained from the path 1. This time, the in-plane shear stress r 12 is held longer for the path 2 than for the path 1 and, consequently, the composite has more time to creep in in-plane shear. The difference in the strain amplitudes is much more important in in-plane shear, as in this case the behaviour of the composite is more matrix-dominated and therefore exhibits more creep.