Multiscale Study of the Effect of Fiber Twist Angle and Interface on the Viscoelasticity of 2D Woven Composites

Time and temperature affect the viscoelasticity of woven composites, and thus affect their long-term mechanical properties. We develop a multiscale method considering fiber twist angle and interfaces to predict viscoelasticity. The multiscale approach is based on homogenization theory and the time–temperature superposition principle (TTSP). It is carried out in two steps. Firstly, the effective viscoelasticity properties of yarn are calculated using microscale homogenization; yarn comprises elastic fibers, interface, and a viscoelastic matrix. Subsequently, the effective viscoelasticity properties of woven composites are computed by mesoscale homogenization; it consists of homogenized viscoelastic yarns and matrix. Moreover, the multiscale method is verified using the Mechanics of Structure genome (MSG) consequence. Finally, the effect of temperature, fiber twist angle, fiber array, and coating on either the yarn’s effective relaxation stiffness or the relaxation moduli of the woven composite is investigated. The results show that increased temperature shortens the relaxation time of viscoelastic woven composites, and fiber twist angle affects tensors in the relaxation stiffness matrix of the yarn; the coating affects the overall mechanical properties of woven composites as well.


Introduction
Woven composites have the advantage of high strength, high stiffness, good oxidation resistance, and excellent thermal stability. They have been widely utilized in aerospace, national defense, biomedical, and other industrial fields. Woven composites are often subjected to long-term load in engineering practice, which seriously affects their viscoelastic behavior leading to the failure of the structure, especially at high temperature. Stress relaxation occurs when viscoelastic materials are subjected to static or variable strain continuously and this phenomenon is called creep [1]. Researchers have now developed a number of approaches to investigate the viscoelasticity of composites. Hashin [1,2] proposes the elastic viscoelastic correspondence principle. Via this principle, the parametric three-dimensional finite-volume direct averaging micromechanics (FVDAM) by Chen et al. [3] and the elastic-based locally exact homogenization theory (LEHT) by Wang and Pindera [4] are developed to accommodate linearly viscoelastic phase response. As an extension of the Eshelby-based Mori-Tanaka (MT) model, Weng et al. [5] investigated the overall viscoelastic behavior of composites with different shapes according to different aspect ratios. Katouzian et al. [6] employ the MT method to study the response of viscoelastic composites vs. the time compared with the experiment data. Moreover, Yang et al. [7][8][9] experimentally studied the long-term creep behavior of fiber-reinforced composite tubes subjected to flexural loading. Martynenko et al. [10,11] conduct numerical simulations and experiments to investigate the effect of the temperature and the lasting time on the effective viscoelasticity of fiber-reinforced composites. Kwok et al. [12] establish a viscoelastic model of single-layer plain woven carbon fiber-reinforced epoxy composites after Materials 2023, 16, 2689 2 of 23 being folded for a period of time compared with experimental measurements to study the deployment of tape springs at different times and temperatures. However, the classical models, elasticity-based homogenization approaches, or mechanics experiments among techniques are presently limited. For instance, the MT model may provide a reasonable estimate of elastic homogenized moduli but not the local stress distributions [3]; the LEHT is limited to composites reinforced by long cylindrical fibers, and experiments cost much money, time, and effort.
In recent years, homogenization theory (HT), the finite element method (FEM), and multiscale analysis have been widely used to predict woven composites' stiffness and strength. For example, Zhao et al. [13] used multiphysics locally exact homogenization theory to investigate the multiscale homogenized thermal conductivity and thermomechanical properties of Advanced European superconductors filament groups (EAS). Pathan et al. [14,15] studied the sensitivity of the viscoelastic response of fiber-reinforced polymers (FRPs) to fiber shapes, fiber volume fraction, interphase volume fraction, and interphase properties via FEM. Also, they compare the Monte Carlo simulations' results on random RVEs to those obtained by analyzing periodic square and hexagonal unit cells. Deviredy et al. [16] investigate the square array RVE (S-RVE) and hexagonal array RVE (H-RVE) of FPRs with the circular and square cross section of fiber and calculate the elastic modulus and thermal conductivity. Moreover, Liu et al. [17] used mechanics of structure genome (MSG) solid and plate model to capture the long-term viscoelastic behaviors of textile composites. Rique et al. [18] extend MSG to construct a linear thermo-viscoelastic model to analyze three-dimensional heterogeneous materials made of constituents with time-and temperature-dependent behavior. Seifert et al. [19] proposed a finite elementbased micromechanical model to obtain the viscoelastic properties of glass fiber composites at high temperatures. Cai et al. [20] propose a multiscale model of three-dimensional four-way woven composites to analyze the influence of braiding angles and fiber volume fractions on the viscoelastic properties by FEM and creep experiment.
Further, the microscopic parameters of woven composites relating to viscoelasticity should be studied more comprehensively. As shown in Figure 1a,b, a traditional twist yarn in a fabric consists of fibers at an angle with the direction of the yarn axis. Twisted fiber bundles characterized by the twist angel, distinguishing significantly from untwisted ones in the mechanical property, contribute to the mechanical properties of the composite through fiber-yarn-fabric sequence. In Figure 1d, a region exists at the boundary between fibers and matrix, which possesses mechanical properties different from those of the fibers and the matrix due to physical and chemical reactions between the two main phases; this region is usually modeled as a physical coating of the continuum [21], so that the thickness, the material property, and the bonding condition of the coating significantly affect the overall mechanical properties of the woven composite. Some relevant studies are as follows. Miao et al. [22,23] attempt to optimize the yarn structure of fiber-reinforced polymer composites and establish the relationship between the fiber twist angle and mechanical properties of unidirectional fiber composites. Xiong et al. [24] construct a multiscale mechanical model and study the influence of fiber twist angle on the mechanical properties of plain woven composites. Fisher and Brinson [25] analyze the mechanical response of fiber-reinforced polymer matrix composites with viscoelastic interface regions by Mori-Tanaka micromechanical model and study the physical aging of viscoelastic composites. To obtain the viscoelastic behavior of polymer-based heterogeneous materials, Huang et al. [26] developed 3D viscoelastic calculated grains (CGs) containing spherical inclusions, interfacial phases/coatings, and non-interfacial phases/coatings. Yang et al. [27] theoretically study the frequency-and temperature-dependent viscoelastic behavior of the short fiber-reinforced polymers (SFRPs) and consider interface/interface conditions. Although these applications provide various aspects of viscoelastic behavior of composites, they are limited to composites with infinitesimally small microstructures. consider interface/interface conditions. Although these applications provide various aspects of viscoelastic behavior of composites, they are limited to composites with infinitesimally small microstructures.  [24]; (b) The orientation of sisal fibers in the yarn [28]; (c) The fiber rod end surface [29]; (d) The fiber and coating in the yarn [30]. (Reproduced with permission ref. [24]. Copyright 2018 Elsevier; ref. [28]. Copyright 2016 Elsevier; ref. [29]. Copyright 2020 MDPI; ref. [30]. Copyright 2010 Elsevier).
Usually, multiscale methods are classified as hierarchical and concurrent methods by Belystchko and Song [31]. Based on the structure of woven composite having hierarchy and feature, this paper develops a hierarchical multiscale method to predict the viscoelastic of woven composite. The technique can study the influence of microscopic parameters on the long-term viscoelastic behavior of two-dimensional woven composites. Considering the non-isotropic viscoelastic materials cannot be calculated directly by commercial finite element software, the mesoscale estimate in the multiscale method introduces the discretization theory. Therefore, the long-term viscoelastic properties of yarns are divided into multiple transient properties with equal intervals. The long-term anisotropic viscoelastic simulation of RVE2 is equivalent to the instantaneous elastic simulation. This method is easy to calculate the viscoelastic properties of 2D woven composites. We were, moreover, using the multiscale approach to investigate the effect of more factors, such as twist yarn, coating interface, fiber array, and ambient temperature.
The remainder part of the paper is organized as follows: Section 2 establishes a multiscale RVE structure of a 2D woven composite and defines the fiber twist angle. Section 3 describes a viscoelastic multiscale model considering fiber twist angle and coating interface. Section 4 verifies the accuracy of the multiscale method via MSG methods. In Section 5 we investigate the combined effects of temperature, fiber twist angle, coating interface, and array type on the homogenized viscoelasticity moduli of woven composites, reporting new results. Conclusions are presented in Section 6.  [24]; (b) The orientation of sisal fibers in the yarn [28]; (c) The fiber rod end surface [29]; (d) The fiber and coating in the yarn [30]. (Reproduced with permission ref. [24]. Copyright 2018 Elsevier; ref. [28]. Copyright 2016 Elsevier; ref. [29]. Copyright 2020 MDPI; ref. [30]. Copyright 2010 Elsevier).
Usually, multiscale methods are classified as hierarchical and concurrent methods by Belystchko and Song [31]. Based on the structure of woven composite having hierarchy and feature, this paper develops a hierarchical multiscale method to predict the viscoelastic of woven composite. The technique can study the influence of microscopic parameters on the long-term viscoelastic behavior of two-dimensional woven composites. Considering the non-isotropic viscoelastic materials cannot be calculated directly by commercial finite element software, the mesoscale estimate in the multiscale method introduces the discretization theory. Therefore, the long-term viscoelastic properties of yarns are divided into multiple transient properties with equal intervals. The long-term anisotropic viscoelastic simulation of RVE2 is equivalent to the instantaneous elastic simulation. This method is easy to calculate the viscoelastic properties of 2D woven composites. We were, moreover, using the multiscale approach to investigate the effect of more factors, such as twist yarn, coating interface, fiber array, and ambient temperature.
The remainder part of the paper is organized as follows: Section 2 establishes a multiscale RVE structure of a 2D woven composite and defines the fiber twist angle. Section 3 describes a viscoelastic multiscale model considering fiber twist angle and coating interface. Section 4 verifies the accuracy of the multiscale method via MSG methods. In Section 5 we investigate the combined effects of temperature, fiber twist angle, coating interface, and array type on the homogenized viscoelasticity moduli of woven composites, reporting new results. Conclusions are presented in Section 6.

Multiscale Framework and RVEs
The 2D and 3D woven composites have complex hierarchical structures composed with matrix and interwoven yarns and the fill yarn and warp yarn have periodical characteristics, respectively. The mechanical properties of woven composites mainly depend on the intrinsic properties of the fiber and matrix, the yarn, and the woven composite structure. The yarn structure, such as its geometry, internal fiber orientation, and fiber volume fraction, determines the main properties of the woven composite, especially mechanical properties [32]. As shown in Figure 2, it can be seen that the hierarchical structures of 2D woven composites consist of three scales: microscale, mesoscale, and macroscale. We start from the modeling of the yarn as RVE1 at microscale, which consists of twisted or untwisted fibers distributed in the matrix and there is the interface between the fiber and matrix; in this paper, we account the fiber with hexagonal or square array in matrix. The mesoscale as RVE2 consists of warp and fill yarns embedding in the matrix. The macroscale is 2D woven composites which is composed with RVE2s. The numerical examples of the latter two scales used in this paper were generated and analyzed by Texgen software (v3.12.0, University of Nottingham, Nottingham, UK) [33] and finite element analysis software ABAQUS 2020. Based on the "bottom-up" analysis process, we predict the overall performance of macroscopic 2D woven composites from the inherent properties of microscopic fibers and matrices. Fibers are assumed to be linear elastic and transversely isotropic, while the matrix is viscoelastic and isotropic.

Multiscale Framework and RVEs
The 2D and 3D woven composites have complex hierarchical structures composed with matrix and interwoven yarns and the fill yarn and warp yarn have periodical characteristics, respectively. The mechanical properties of woven composites mainly depend on the intrinsic properties of the fiber and matrix, the yarn, and the woven composite structure. The yarn structure, such as its geometry, internal fiber orientation, and fiber volume fraction, determines the main properties of the woven composite, especially mechanical properties [32]. As shown in Figure 2, it can be seen that the hierarchical structures of 2D woven composites consist of three scales: microscale, mesoscale, and macroscale. We start from the modeling of the yarn as RVE1 at microscale, which consists of twisted or untwisted fibers distributed in the matrix and there is the interface between the fiber and matrix; in this paper, we account the fiber with hexagonal or square array in matrix. The mesoscale as RVE2 consists of warp and fill yarns embedding in the matrix. The macroscale is 2D woven composites which is composed with RVE2s. The numerical examples of the latter two scales used in this paper were generated and analyzed by Texgen software (v3.12.0, University of Nottingham, Nottingham, UK) [33] and finite element analysis software ABAQUS 2020. Based on the "bottom-up" analysis process, we predict the overall performance of macroscopic 2D woven composites from the inherent properties of microscopic fibers and matrices. Fibers are assumed to be linear elastic and transversely isotropic, while the matrix is viscoelastic and isotropic.

Twist Angle and Coating of Fiber in the Yarn
The viscoelastic properties of every scale model are subject to the components and structure. Yarns play an important role in the woven composites. Effective properties of yarn are computed based on fibers and matrix. The twisted yarn axis has an angle  with the fiber axis, and the angle is called a torsion angle or twist angle as shown in Figure 1a,b. During the injection molding and extrusion process [34], the twist angle of a fiber, related to the fiber orientation, depends on its radical location in the yarn's cross section and changes gradually and continuously from the interior of the yarn to its surface, as shown in Figure 3, so the fiber orientation distribution function can be defined as seen in the following. Xiong et al. [24] used a spatial coordinate system to define the orientation of

Twist Angle and Coating of Fiber in the Yarn
The viscoelastic properties of every scale model are subject to the components and structure. Yarns play an important role in the woven composites. Effective properties of yarn are computed based on fibers and matrix. The twisted yarn axis has an angle θ with the fiber axis, and the angle is called a torsion angle or twist angle as shown in Figure 1a,b. During the injection molding and extrusion process [34], the twist angle of a fiber, related to the fiber orientation, depends on its radical location in the yarn's cross section and changes gradually and continuously from the interior of the yarn to its surface, as shown in Figure 3, so the fiber orientation distribution function can be defined as seen in the following. Xiong et al. [24] used a spatial coordinate system to define the orientation of fibers in the yarn in order to describe the fiber orientation distribution (FOD). As shown in Figure 3b [34], 1-axis is the principal axis of the fiber, and the 2-axis and the 3-axisareradialaxes perpendicular to the 1-axis. x 1 , y 1 and z 1 are the coordinate system of the RVE1. θ is defined as the angle between 1-axis and the z 1 -axis of the spatial coordinate system; the value range is (−θ 0 , θ 0 ). θ 0 is the twisting of the fibers on the surface of the yarn. ϕ between the projection of the fiber on the x 1 − y 1 plane and the y 1 -axis; the value range is (0, 2π).
in Figure 3b [34], 1-axis is the principal axis of the fiber, and the 2-axis and the 3axisareradial-axes perpendicular to the 1-axis. 1 x , 1 y and 1 z are the coordinate system of the RVE1.  is defined as the angle between 1-axis and the 1 z -axis of the spatial coordinate system; the value range is  is the twisting of the fibers on the surface of the yarn.  between the projection of the fiber on the 1 1 x y  plane and the 1 y -axis; the value range is (0, 2 )  . The twist angle of the fiber in the yarn depends on the radial position of the yarn cross section, and the twist angle of the fiber near the surface of the yarn is greater than the twist angle inside the yarn, as shown in Figure 3b. The fiber orientation distribution function can be defined as [24]: Although fibers in the yarn are randomly arranged, the microstructure of the yarn is usually idealized in the literature as a square or hexagonal pack, and representative volume element is S-RVE or H-RVE, as shown in Figure 4. To be considered are the interfaces between the fiber and matrix (shown in Figure 1d) and it is assumed the interfaces or coating are perfect enough to connect the fiber and matrix together. a is the diameter of the fiber and b is the exterior diameter of the coating. 1 l , 2 l and 3 l are the lengths of RVEs. In Figure 5, the mesoscale RVE2 model was created by Texgen software.  The twist angle of the fiber in the yarn depends on the radial position of the yarn cross section, and the twist angle of the fiber near the surface of the yarn is greater than the twist angle inside the yarn, as shown in Figure 3b. The fiber orientation distribution function can be defined as [24]: Although fibers in the yarn are randomly arranged, the microstructure of the yarn is usually idealized in the literature as a square or hexagonal pack, and representative volume element is S-RVE or H-RVE, as shown in Figure 4. To be considered are the interfaces between the fiber and matrix (shown in Figure 1d) and it is assumed the interfaces or coating are perfect enough to connect the fiber and matrix together. a is the diameter of the fiber and b is the exterior diameter of the coating. l 1 , l 2 and l 3 are the lengths of RVEs. In Figure 5, the mesoscale RVE2 model was created by Texgen software. w and h are the width and thickness of the yarn embedding in the matrix of L z thickness and L x length in the periodical interval c.

Viscoelastic Constitutive Model
The constitutive relationship of viscoelastic materials depends on the time and the temperature, and the uniaxial stress-strain relationship of isotropic viscoelastic materials can be represented by the Boltzmann's heredity integral: D t is the creep modulus, and the elastic modulus ( ) E t can be defined by the Prony series over a wide range of timescale [35]: E  is the long-term modulus, i E are the Prony coefficients, and i  are the relaxation times, whose value is the larger, the slower the stress relaxation decays.

Viscoelastic Constitutive Model
The constitutive relationship of viscoelastic materials depends on the time and the temperature, and the uniaxial stress-strain relationship of isotropic viscoelastic materials can be represented by the Boltzmann's heredity integral: is the time-dependent strain, ( ) D t is the creep modulus, and the elastic modulus ( ) E t can be defined by the Prony series over a wide range of timescale [35]: E  is the long-term modulus, i E are the Prony coefficients, and i  are the relaxation times, whose value is the larger, the slower the stress relaxation decays.

Viscoelastic Constitutive Model
The constitutive relationship of viscoelastic materials depends on the time and the temperature, and the uniaxial stress-strain relationship of isotropic viscoelastic materials can be represented by the Boltzmann's heredity integral: t is a time variable, τ is a passed time, σ(t) is the time-dependent stress, ε(t) is the timedependent strain, D(t) is the creep modulus, and the elastic modulus E(t) can be defined by the Prony series over a wide range of timescale [35]: E ∞ is the long-term modulus, E i are the Prony coefficients, and ρ i are the relaxation times, whose value is the larger, the slower the stress relaxation decays. Temperature has a significant influence on viscoelastic behavior, and TTSP (the timetemperature superposition principle) provides a theoretical foundation to study the longterm behavior of viscoelastic materials by short-term characterization experiments. This principle relates the temperature-dependent relaxation modulus to the time-dependent one, just as E(t, T) = E(t , T 0 ), which indicates that the modulus at temperature T and time t is equal to that at the reference temperature T 0 and shortened time t . Therefore, we can cut down the observation time to study the long-term viscoelastic behavior of the polymer by increasing the temperature of the experimental environment, and this equivalence can be achieved by means of a conversion factor a T , which is the ratio of the relaxation time at two different temperatures: where T is the general temperature, and T 0 is the reference temperature that is often taken as the glass transition temperature. In order to illustrate the dependence of viscoelasticity on temperature, the viscoelastic material studied in this paper is a thermos-rheological simplicity material and the conversion factor a T of the material is applicable over the whole relaxation time. When the temperature changes with time, the shortened time is obtained by integration: According to TTSP, the temperature shift function can be expressed by the Williams-Landel-Ferry (WLF) equation [36]: here c 1 and c 2 are material-dependent constants and the logarithm base is 10. Thus, the viscoelastic constitutive relation Equation (5) at temperature T can be rewritten as: For the anisotropic material, the stress-strain relationship can be expressed in the time domain as: Then, the relaxation modulus tensor can be expressed as

Multiscale Homogenization
The homogenization theory (HT) is employed by [3,4,13]. The theory can harvest the equivalent material properties of the composite material having periodic feature. In the twisted yarn, the helical fibers in different radial positions have different twist angles (shown in Figure 3). Thus, the stiffness matrix of each fiber and coating at the twist angle can be calculated as [24]: Since the fiber is an elastic material, its stiffness matrix does not change with time. A series of helical fibers located at the point of the same radius in the yarn cross section can be assumed as a layer, whose orientation changes continuously from the center to the surface of the yarn. Based on the directional average method and the iso-strain hypothesis [34,37], the stiffness matrix of each helix layer is integrated by helix angle to obtain the average stiffness matrix after fiber and coating torsion [24]: The yarn is composed of fibers, coatings, and the matrix. The first two of these are elastic and the matrix is made from the polymer of viscoelastic elasticity, so the yarn is transversely isotropic. According to the model of the RVE1 established in the previous section, the overall average stress component of the twisted yarn can be expressed as: where f , int and m represent the fiber, interfacial phase and matrix, and the superscript "~" in the formulation indicates the homogenized sign of linear viscoelasticity. σ RVE1 i (t) is the average time-varying stress field in RVE1, C y i (t) is the time-varying relaxation stiffness tensor in RVE1, and V RVE1 is the volume.
The in-plane structure of the 2D woven composite material is symmetrical and periodical, and the mesoscale model is in-plane isotropic. The relaxation stiffness matrix of warp and fill yarns in the global coordinate system can be described as: Based on the RVE2 built in the previous section, the average stress component of the 2D woven composite can be expressed as: In the formulation, the subscript λ means that the variable belongs to the corresponding phase, such as the warp, the fill, and the matrix.
We can apply the corresponding normal strain or shear strain to the single-cell model of different scales, and then solve the viscoelastic partial differential equation according to the above analytical models to restore the stress distribution in the numerical microelement, and finally calculate the macroscopic overall average stress component of the 2D woven composite material in the form of integration: The constitutive equations of 2D woven composites are represented in the form of relaxation stiffness matrix engineering:

Periodic Boundary Conditions of RVEs
Periodic boundary conditions must be applied to RVE of different scales to ensure uniform deformation and continuity of RVE under various loading conditions. In this paper, the periodic boundary condition proposed by [38] is adopted, and the displacement equation for each pair of nodes on the parallel boundary surface of RVE is as follows: where the superscripts j+ and j− refer to the two corresponding nodes on the boundary of RVE and ε ik is the average strain of RVE. x k is the coordinate value of k-axis and u * i is the periodic displacement of i-axis. The relative displacement between two corresponding nodes on a parallel boundary surface can be expressed as: Equation (26) can be easily implemented by setting the displacement constraint on each pair of coupled nodes on the boundary of RVE.

Validation
To verify the result based on the multiscale model and homogenization viscoelasticity to every scale RVE, the results at the microscale predicted in this paper are compared with that of the literature [17]. In this case, the yarn has no coating and twist angel, and the fiber volume fraction v f is fixed at 0.64 in the yarn. The elastic properties of the fiber are shown in Table 1, the viscoelasticity of the matrix is expressed by the relaxation time and Prony coefficient, as shown in Table 2, and the Poisson's ratio is assumed to be constant at 0.33.   Figure 4a. Each analysis consists of two steps and the relaxation stiffness is defined with the option * VISCO in ABAQUS. The first step lasting a short time period (0.1 s) is to apply a unit strain on RVE1. In the second step, the strain remains constant and lasts 10 10 s.
In Figure 6, the results of microscale S-RVE1 and H-RVE1 in this paper are compared with that of Liu-MSG in [17], whose work only focuses on S-RVE1 and both results agree well. Additionally, results of two RVE1s have a small difference in axial relaxation stiffness tensor C RVE1 11 , the maximum difference of which is between 0.05~2.2%, while the maximum difference in radial relaxation stiffness tensor C RVE1 22 is about 2.1~8.4%. The effect of time for tensor C RVE1 11 is almost negligible, because the behavior in this direction is dominated by the fibers, whose elastic behavior is independent of time. Other tensors of the relaxation stiffness decrease with the time increased, exhibiting a trend of variation over time. More importantly, the microstructure has more significant effect on the radial relaxation tensor than that on the axial relaxation one, which can be explained by the dominated contribution of fibers to the axial direction of the yarn and the behavioral gap between the matrix and the fiber which is not obvious in the radial direction.
We employed the FEM to get the homogenization viscoelasticity prediction for the woven composite. The width and thickness L z of the warp or fill in Figure 5 are 0.9 mm and 0.06 mm, respectively, while the thickness of the fabric is 0.12 mm, and the interval of neighboring warps or fills is 1.75 mm. Progressive meshes with C3D8R and C3D10 solid elements are adopted in ABAQUS to guarantee the convergence, and co-node grids are set up at the interface between the yarn and the matrix. It is well established that yarn and RVE2 are not isotropic, and finite element software cannot directly define the viscoelasticity of the anisotropic material [17]. However, the long-term relaxation viscoelasticity of RVE2 has consisted of many transient responses. Based on the discretization theory, the long-term relaxation response of RVE2 can be decomposed into multiple transient reactions with the same interval, simplifying the finite element calculation. The yarn's transient properties are summarized in Tables A1 and A2 of Figure 7, the results of the paper are compared with that of Rique-MSG in [18], and both of results are in good agreement. Figure 7 shows the relaxation moduli of plain woven composites by viscoelastic multiscale homogenization. S-RVE2 and H-RVE2 represent the relaxation moduli of the mesoscale corresponding to S-RVE1 or H-RVE1, respectively. Rique-MSG represents the relaxation moduli of plain woven composites calculated by Rique et al. [18] using the Mechanics of Structure genome MSG under periodic boundary conditions. In Rique's research, the arrangement of fibers in the yarn is a square array. Therefore, in Figure 7, the results of S-RVE2 are closer to those of Rique-MSG, and the consequences of S-RVE2 and H-RVE2 have noticeable differences only in the in-plane tensile moduli. In contrast, the differences in other moduli are slight.

Numerical Investigations
Based on the multiscale homogenization, this section investigates the effect of the temperature and several other microstructural parameters, such as the twist angle of fibers, interfacial parameters, and fiber arrangements, on the relaxation modulus of twodimensional woven composites.

The Effect of Temperature
The relaxation modulus curve of PMT-F4 is obtained by creep experiment [20]. The material-dependent constants c 1 , c 2 are found by fitting data of the conversion factor vs. temperature data to the WLF equation (Equation (11)) and then we can infer that c 1 = 28.3816, c 2 = 93.291 at the reference temperature T 0 = 40 • C. This section assumes that the elastic properties of fibers are independent of the temperature and the fiber volume fraction v f and the twist angle θ in each yarn are chosen to be 0.64 and 0 • , taking no account of the coating's influence, in order to investigate the effect of the temperature on the viscoelastic properties of 2D woven composites. In Abaqus 6.14, the relationship between the viscoelasticity of the matrix and the temperature is defined in the form of the WLF equation and the temperature field is applied to the matrix. Figures 8 and 9 illustrate the relaxation response of yarns and 2D woven composites to different temperatures. With the temperature increasing, the relaxation time of yarn and woven composites gradually shortens, so the effect of temperature on the viscoelasticity of the yarn and the woven composite conforms to TTSP. Moreover, the RVE1's (yarn) relaxation tensor C RVE1 ij and the RVE2 s (woven composite) relaxation modulus E RVE2 ij (G RVE2 ij ) decrease sharply in the time domain with the temperature increasing, and the higher the temperature leads to the faster relaxation. The phenomenon is due to the molecular movement stimulated by the temperature. When the temperature rises, the viscous frictional resistance between molecules becomes smaller and the relaxation reaction of molecules is easier. In addition, there is no significant distinction between H-RVE1 and S-RVE1 in the trend of stiffness vs. time curve and the relaxation time. Figure 8a-c,f show that tensors C RVE1

22
, and C RVE1 55 in the H-RVE1 are significantly lower than that of S-RVE1, while the tensors C RVE1 23 and C RVE1 44 of H-RVE1 are significantly higher than that of S-RVE1 in Figure 8d,e. At the mesoscale, as illustrated in Figure 9, the impact of RVE1's structure on RVE2's homogenized viscoelastic moduli E RVE2   Figure 9. Effect of temperature on the effective relaxation modulus of 2D woven composites.

The Effect of the Fiber Twist Angle
After studying the effect of the temperature on the viscoelasticity of RVEs of different scales, the influence of the twist angle is investigated considering θ = 0 • , 30 • , and 60 • at the temperature fixed at 60 • C in this section. The effective stiffness matrix of the twist fiber calculated by Equation (15) in Section 3.2 is shown in Table 3. In spite of the twist angle, the transverse elastic moduli E f 2 equals to E f 3 as well as the shear moduli G f 12 equals to G f 13 , so the twist fiber can be considered transverse isotropic. Table 3 shows that the twist angle of the fiber is increasing, and the effective elastic properties of the fiber change dramatically where these radial properties are effectively improved.  Figure 10 shows the variation of the relaxation stiffness of the yarn affected by a wide range of the twist angle. The relaxation stiffness tensor C RVE1 11 decreases sharply with the increase of the twist angle, but the other relaxation stiffness components increase. In Figure 10a, the effective normal stiffness tensor C RVE1 11 of the yarn with the twist angle of 30 • is only about 84% of the twist less yarn, while the ratio decrease to 30% when the twist angle is 60 • . In Figure 10b,f, the twist angle promotes a large improvement on C RVE1 12 and C RVE1 55 . Surprisingly, the value of C RVE1 12 on the twist angle 60 • gradually decreases with the time increasing and finally is lower than that on the twist angle 30 • . The effect of the high twist angle is more prominent on C RVE1 22 , C RVE1 23 , and C RVE1 44 , but this effect will diminish after a period of time. Additionally, curve trends of the two microstructures H-RVE1 and S-RVE1 are almost the same, but their values are different. Figure 11 shows the effect of the fiber twist angle on the relaxation moduli of a twodimensional woven composite. It can be seen that with the increase of the fiber twist angle, the in-plane relaxation moduli E RVE2 xx and E RVE2 yy of the woven composite decrease sharply, but the variation of the out-plane relaxation modulus E RVE2 zz and the shear moduli G RVE2 xy , G RVE2 xz , and G RVE2 yz induced by the fiber twist angle is not obvious. Therefore, the fiber twist angle can be realized as a non-main parameter in the optimal design of 2D woven composites. In addition, there is no significant difference in E RVE2  Figure 11. Effect of the twist angle on the effective relaxation moduli of 2D woven composite.

The Effect of the Coating Thickness
The physical coating between the fiber and matrix is accepted as a kind of component existing in the yarn [21]. To demonstrate the effect of the coating thickness which is represented by the ratio b/a, b is the coating diameter and a is fiber diameter, respectively, shown in Figure 4, b/a = 1 indicates that there is no coating between the fiber and the matrix, while b/a > 1 means the coating is surrounding the fiber. The field temperature is fixed at 60 • C and the volume fraction v f of the untwist fiber in each yarn is 0.64. The mechanical properties of Pyrolytic carbon used for the coating are given in Table 4. The coating thickness significantly improves the viscoelastic properties of the yarn and affects the axial and radial relaxation components in varying degrees, just as shown in Figure 12. The normal stiffness tensor C RVE1 11 of b/a = 1.05 is greatly improved compared with that of b/a = 1,and the microstructure of the hexagonal array has greater impact on C RVE1 11 than that of the square array. However, the effect of the coating on C RVE1 12 is very insignificant. At the initial moment t = 1 s, C RVE1 12 at the ratio b/a = 1.05 is lower than that at b/a = 1; with the increase of the time, C RVE1 12 gradually decreases and becomes higher than the value of b/a = 1. The coating thickness also improves the radial relaxation stiffness tensor C RVE1 22 and C RVE1 23 to an extent as illustrated in Figure 12c,d. In Figure 12e,f, the coating has a greater effect on the axial shear tensor C RVE1 55 than the radial shear tensor C RVE1 44 , which depends on the property of the coating, whose axial mechanical stiffness is much higher than the radial one.
In Figure 13, the coating thickness significantly improves the in-plane relaxation moduli E RVE2 xx and E RVE2 yy of the 2D woven composite, but the out-plane relaxation modulus E RVE2 zz and the shear moduli G RVE2 xy , G RVE2 yz are slightly promoted. In general, excellent coating properties are necessary to improve the overall mechanical properties of the 2D woven composite. In Figure 13a Figure 13. Effect of the coating thinkness on effective relaxation moduli of 2D woven composites.

Conclusions
This paper investigates the viscoelastic properties of 2D woven composite employing elastic-based homogenization theory to establish a multiscale model considering the coating and the fiber twist angle. According to TTSP, the influence of temperature on the relaxation time of the 2D woven composites is studied. The effects of the fiber twist angle and the coating on the viscoelasticity are discussed in detail. The conclusions are summarized as below: (1) The results using the multiscale method show that the fiber array considerably affects stiffness relaxation of the yarn. C RVE1

22
of the square array are higher than H-RVE1, while the H-RVE1's tensors C RVE1 are more important, while the improvement on other modulus is minor, and gradually disappears with the time. In addition, the lower twist angle has a more significant effect on the axial stiffness of the yarn, while the radial stiffness is more sensitive to the higher angle.
(3) The coating, the material property, and the thickness can effectively improve the overall viscoelasticity of 2D woven composites, especially the in-plane relaxation moduli. When the stiffness of the coating is higher than that of the matrix, the coating will effectively improve the overall mechanical properties of the composite. Designing the coating is significant in exploiting the potentiality of 2D woven composites. (4) The multiscale method facilitates the calculation of the viscoelasticity of woven composites. Combining with the discrete theory and FEM, this paper provides an appropriate approach for analyzing non-isotropic composites. In addition, the effect of more microscopic parameters on the mesoscopic properties is considered and calculated.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Each fiber is assumed to be transverse isotropic and their stiffness matrix is given as follows [39]:  where m = cos α 0 , n = sin α 0 ,u = cos α 0 , v = sin α 0 and α 0 is undulation angle. The relaxation stiffness [ C RVE2 ] of RVE2 to calculate relaxation moduli is as follows: