Thermo-Viscoelastic Response of 3D Braided Composites Based on a Novel FsMsFE Method

A homogenization-based five-step multi-scale finite element (FsMsFE) simulation framework is developed to describe the time-temperature-dependent viscoelastic behavior of 3D braided four-directional composites. The current analysis was performed via three-scale finite element models, the fiber/matrix (microscopic) representative unit cell (RUC) model, the yarn/matrix (mesoscopic) representative unit cell model, and the macroscopic solid model with homogeneous property. Coupling the time-temperature equivalence principle, multi-phase finite element approach, Laplace transformation and Prony series fitting technology, the character of the stress relaxation behaviors at three scales subject to variation in temperature is investigated, and the equivalent time-dependent thermal expansion coefficients (TTEC), the equivalent time-dependent thermal relaxation modulus (TTRM) under micro-scale and meso-scale were predicted. Furthermore, the impacts of temperature, structural parameters and relaxation time on the time-dependent thermo-viscoelastic properties of 3D braided four-directional composites were studied.


Introduction
Three-dimensional (3D) braided composites are increasingly being used in the manufacturing of structural parts due to their incomparable advantages compared with traditional laminated composites in recent years. In order to ensure the safety and reliability in use of the 3D braided composite, the mesoscopic configuration modeling and mechanical performance prediction under multiple working conditions [1][2][3][4][5][6] of 3D braided composites have recently received considerable attention by many scholars. It is well known that the main body of the braided composites is matrix material, and polymer as matrix materials shows obvious viscoelasticity in composites. Hence, the braided composites will also exhibit considerable viscoelastic behavior. With the rapid development of science and technology, various approaches [7][8][9][10][11][12][13][14][15][16][17] have been devised to deal with the viscoelastic problems of 3D braided composites. Liu et al. [7] and Li et al. [8] studied the creep behavior of resin-based 3D braided composites by tensile creep test, and discussed the influence of braided structure and fiber volume fraction on the creep behavior of materials. Priyank et al. [9] set up a three-dimensional micromechanical analytical model based on the classical laminate theory, and studied the equivalent elastic and viscoelastic properties of a 3D braided composites. In addition, the multi-scale method was extensively to predict the viscoelastic behaviors of 3D braided composites. Yuan et al. [10] used the multi-scale method to study the viscoelastic properties of 3D braided composites, and gave the variation law of viscoelastic properties with process parameters. Cai et al. [11] studied the viscoelasticity of 3D braided composites based on the interior cell model and yarn model, and discussed the influence of braiding parameters on the viscoelasticity behavior. Mourid et al. [12] developed an analytical homogenization technique to predict the viscoelastic behavior of 3D braided composites by extending the elastic homogenization frame and coupling the Laplace-Carson transforms. Zhou et al. [13] studied the elastic constants and creep performance by a multiple scale viscoelastic model for the 3D braided composites. Liu et al. [14] studied the viscoelastic behaviors of 3D braided composites based on a yarn solid model and mechanics of structure genome plate model. Jia et al. [15] analyzed the nonlinear viscoelastic responses and damage mechanisms of 3D orthogonal composite under quasi-static tensile loading based on a micro/meso-scale repeated unit cells model. Tate et al. [16] presented the viscoelastic effects on fatigue behavior of 3D carbon/epoxy braided composites based on the experimental methods. Zhai et al. [17] proposed a multi-scale calculation method to study the stress relaxation response of 3D braided composites at three scales, and characterized the viscoelastic properties at different scales.
Thermoviscoelasticity is an important thermo mechanical property for failure and failure of resin matrix composites, which is often used in the long-term and variable temperature environment. Hwang et al. [18] predicted the stress and deformation histories of notchless and notched graphite/epoxy composites under mechanical and thermal loads based on finite-element formula. Muddasani et al. [19] presented some experimental works and finite element analyses of nonlinear thermo-viscoelastic behaviors for multilayered composites. Li et al. [20] established a 3D incremental viscoelastic constitutive model, the residual stresses and the curing deformation of T-shaped composite structures during curing studied. Similarly, it is particularly significant to study the thermo-viscoelastic behaviors before the 3D resin-based braided composite can be confidently used in the primary structures. Seifert et al. [21] studied the time-temperature dependent material properties of 3D braided composite material based on experimental and numerical methods. Hirsekorn et al. [22] obtained the viscoelastic behavior of 3D braided composites from the dependence of component on temperature and degree of cure by a homogenization strategy. Cai et al. [23] reported the thermo-viscoelastic properties of 3D braided composites in different temperature states based on the three-cell model. Sun et al. [24] predicted the dynamic thermo-viscoelasticity response of 3D braided composites to alternating stress load at room temperature, and studied the relationship between the complex compliance and angular frequency based on the mesoscopic representative unit cell (RUC).
These great results are of great significance to explore the thermo-viscoelastic performance of 3D braided composites, but the research on the thermo-viscoelasticity behavior of 3D braided composites are mostly based on the meso-scale RUC and the time-temperaturedependent stress distributions have not received special attention. The time-temperaturedependent mechanical behavior correlation at different scales is not clear, some alternative numerical methods are necessary. To address this drawback, this paper presented a homogenization-based five-step multi-scale finite element (FsMsFE) framework for handling the thermo-viscoelasticity behaviors of 3D braided four-directional composites based on the correspondence principle and time-temperature equivalence principle. In the micro-scale, each kind of fiber and matrix is modeled separately. In the meso-scale, the reinforcement structure is modeled and yarns composed of several thousand of fibers are regarded as homogeneous material. In the macro-scale, the structure is modeled as homogeneous material. Based upon the Prony series fitting, the characters of the thermal stress relaxation at three scales subject to variation in temperature are investigated, and the equivalent time-dependent thermal expansion coefficients (TTEC), the equivalent time-dependent thermal relaxation modulus (TTRM) in micro/meso-scale were predicted. Furthermore, the effects of braiding parameters, relaxation time and temperature on the time-dependent thermo-viscoelastic properties of 3D braided four-directional composites were analyzed. This research could also be extended to the study of thermo-viscoelasticity behaviors of 3D multidirectional polymer-based braided composites or the inhomogeneous materials with periodic structures.

Multi-Scale Coupling for 3D Braided Composites
In fact, the change of material properties at a lower scale plays an important role in the failure process of macro-structure, but it is often difficult to capture the mechanical response of materials from the perspective of calculation. According to the periodical characteristic of 3D four-directional braided architecture, the numerical models and thermo-viscoelasticity behavior analysis of 3D four-directional braided composites are carried out on three-scales, such as the fiber/matrix scale, the fiber tow/matrix scale, and the homogeneous 3D braided composites scale. As shown in Figure 1, the scheme for the multi-scale modeling framework of 3D braided four-directional composites is firstly presented. The coupling technique for the multi-scale thermo-viscoelasticity behaviors prediction of 3D braided four-directional composites is given as following.

Multi-Scale Coupling for 3D Braided Composites
In fact, the change of material properties at a lower scale plays an important role in the failure process of macro-structure, but it is often difficult to capture the mechanical response of materials from the perspective of calculation. According to the periodical characteristic of 3D four-directional braided architecture, the numerical models and thermo-viscoelasticity behavior analysis of 3D four-directional braided composites are carried out on three-scales, such as the fiber/matrix scale, the fiber tow/matrix scale, and the homogeneous 3D braided composites scale. As shown in Figure 1, the scheme for the multi-scale modeling framework of 3D braided four-directional composites is firstly presented. The coupling technique for the multi-scale thermo-viscoelasticity behaviors prediction of 3D braided four-directional composites is given as following.

Association of Multiple Scales
The homogenization method is usually regarded as an effective method to solve multi-scale problems of composite materials with periodic microstructures [25][26][27], the multi-scale thermo-viscoelastic problem of 3D braided composites is solved based on the homogenization method. Since characterization of the resin-based 3D braided composites requires scale transitions between more than two scales, different models were defined. On the macro scale, the meso-scale RUC is assumed periodic. Similar to the macro-scale model, meso-scale RUC is composed of yarns and matrix, yarns are assumed to be formed from micro scale fiber bundle RUC. Let x, y, z denote the coordinate system of the macro-

Association of Multiple Scales
The homogenization method is usually regarded as an effective method to solve multi-scale problems of composite materials with periodic microstructures [25][26][27], the multi-scale thermo-viscoelastic problem of 3D braided composites is solved based on the homogenization method. Since characterization of the resin-based 3D braided composites requires scale transitions between more than two scales, different models were defined. On the macro scale, the meso-scale RUC is assumed periodic. Similar to the macro-scale model, meso-scale RUC is composed of yarns and matrix, yarns are assumed to be formed from micro scale fiber bundle RUC. Let x, y, z denote the coordinate system of the macro-scale, the meso-scale and the micro-scale, respectively. They are related to one another scale by the positive real parameters ε and ω (0 < ε << 1, 0 < ω << 1) as follows. (1)

Equations for the Micro-Scale, Meso-Scale and Macro-Scale
When the yarn of meso-scale RUC is regarded as an in-homogeneous viscoelastic medium with periodical structure, the static thermo-viscoelastic problem in equilibrium condition can be effectively expressed in y and z coordinate systems as follows [28].
where u ε k is displacement components, α kl is thermal expansion coefficient, T ε is temperature increment, f i , t i are body force and traction force, v i is virtual displacement, Ω ε , Γ are the yarn solid part and the boundary of domain. G ijkl is the static relaxation tensor at T of the component materials, which can be calculated by translating along the time axis based on the time-temperature equivalence principle [23], as shown in Figure 2. α T is obtained by the Williams-Landel-Ferry equation and C 1 , C 2 are the material constants. y x y = , z = ω ε .

Equations for the Micro-Scale, Meso-Scale and Macro-Scale
When the yarn of meso-scale RUC is regarded as an in-homogeneous v medium with periodical structure, the static thermo-viscoelastic problem in eq condition can be effectively expressed in y and z coordinate systems as follows  According to the homogenization method [29], the yarn displacements in t are asymptotically expanded into Equation  The stress-strain relationship for a static thermo-viscoelastic material can as [28]. According to the homogenization method [29], the yarn displacements in the ε-space are asymptotically expanded into Equation (3). The terms u i n (n ∈ N 0 ) are z-periodic functions called correctors of order n of the displacement field, and u 0 i (y, t) is the effective displacement only depending on the meso-scopic coordinates y.
The stress-strain relationship for a static thermo-viscoelastic material can be written as [28].
where variables with the mark~show that they are Laplace transformed, s is the Laplace transformed parameter. Laplace transform is used to Equation (3) and take it into Equation (5), equating the terms with the same power of ε, the following expressions can be obtained.
As Equation (7) is linear with respect to u 0 k , its solution u 1 k can be expressed in terms of u 0 k as where χ kl i and Φ represent the z-periodic characteristic displacement for viscoelasticity and thermo-viscoelasticity problem.
Replacing Equation (9) in Equation (7), the following equations can be obtained Then, functions χ kl i and Φ i can be obtained according to Equations (10) and (11). According to Equation (9), Equation (4) can be written as where σ 0 ij (T, y, z, s) is the approximation of the micro-scale stress field in Laplace transformed domain.
Considering Equations (8) and (9) and taking the volume average by f H = 1 |Y| Y f dy, the following result is obtained Comparing with Equation (5) and according to Equations (10) and (11), the Laplace transformation of equivalent relaxation modulus of yarns G H1 ijkl and equivalent TTRM of yarns β H1 ij can be written as Thus, the equivalent stress of the yarn composites in Laplace domain can be written as Under small deformation assumption and according to Liu et al. [28], the linear viscoelastic behavior of the yarn composite can be described in terms of the strain-stress relation by the Laplace inverse transformation for Equation (16).
Hence, Equation (17) can be seen as the equivalent constitutive equation of thermal stress relaxation for viscoelastic composites.
Although the constitutive equation of viscoelastic composites is similar to that of single-phase material expressed by Equation (4), the thermal stress relaxation law of the two materials is different. The thermal stress relaxation law of single-phase materials is the same as that caused by the instantaneous thermal strain α kl T ε , but the thermal stress relaxation law of viscoelastic composites is different. To further illustrate, parameter α H1 kl is introduced and the equivalent TTRM β H1 ij can be written as Therefore, the following form of parameter α H1 kl can be obtained and it can be found that the parameter α H1 kl is time-dependent. At the same time, it can also be known that the strain caused by initial temperature increment has a development process, which cannot be completed instantaneously for viscoelastic composites.
Replacing Equation (17) by Equation (19), the following result is given and the parameter α H1 kl (t) can be defined as equivalent TTEC by comparing with Equation (4).
When the 3D braided composites of macro-scale is regarded as an in-homogeneous viscoelastic medium with periodical meso-scale RUC, the static thermo-viscoelastic problem in equilibrium condition can be effectively expressed in x and y coordinate systems as follows.
where f i is body force, p i , v i are traction force and virtual displacement, G ijkl is the static relaxation tensor of the component materials. u ω k is the displacement, α kl (t) is the timedependent thermal expansion coefficient, Ω ω , Γ are the macro-scale solid part and the boundary of the domain. Based on the homogenization method, a perturbation of the macro-scale displacement vector u 0 i (x, t) is performed with respect to the meso-scale coordinate system y.
Applying the correspondence principle to thermo-viscoelastic analysis of macro-scale and meso-scale, Laplace transform is made to Equations (21) and (22). By substituting Equation (22) into Equation (21) and sorting out the items about ω. Considering the coefficient of ω −1 and ω 0 , the following equations can be obtained according to the homogenization theory [28].
The separation variable method is used to solve Equation (23), where X kl m and Ψ k represent the y-periodic characteristic displacement for viscoelasticity and thermo-viscoelasticity problem, respectively.
Replacing Equation (25) by Equation (23), the characteristic displacement field vector X kl m and Ψ k are the solution of the following equations: Replacing Equation (24) in Equation (25), the Laplace transformation of equivalent relaxation modulus G H2 ijkl , equivalent TTRM β H2 ij and equivalent TTEC α H2 ij for the mesoscale RUC can be similarly expressed as: According to Equations (21) and (22), the approximation of the meso-scale stress field in Laplace transformed domain can be written as Based upon the above derivation, the macroscopic equivalent stress vector σ H2 ij can be expressed as

Inverse Laplace Transformation
Based on the previous derivation, the Laplace transform of homogenized viscoelastic behavior is solved by using the commonly used homogenization technique of elastic materials in the transformation space. In order to obtain this result, the inverse Laplace transform is needed to calculated the homogenized thermal-viscoelastic behavior in the time domain.
Assuming that the viscoelastic properties of 3D braided composites satisfy the following three-parameter model [31], the viscoelastic functions can be represented by the n-order Prony series expressed as follows where the G ∞ , G i , τ i are the long-term modulus, Prony coefficients and relaxation times, respectively.
In the Laplace transformed domain, the viscoelastic functions can be expressed as where the non-linear least square method [32] is used to calculate the coefficients G ∞ , G i , τ i . Then, the equivalent thermo-relaxation response in the time domain can be considered as Equation (33).

Computation of the Effective Properties for Micro-Scale and Meso-Scale
In order to analyze the thermo-viscoelastic behavior of 3D braided composites, the details of the FE method for computing the micro-scale characteristic functions χ kl i and meso-scale characteristic function X kl i are necessary. The FE discretized form of Equations (10) and (26) would be written as where where B is the shape function derivative matrix. G q (q = z, y) is the Laplace transformed relaxation tensor matrix of the micro or meso-scale [13]. φ and ϑ are the discretized form of the corrector χ and X. The right term of Equation (35) is a matrix F q (q = z, y).
To obtain the matrix φ and ϑ in Equation (35), the multi-phase finite element (MPFE) method [33] was employed. A micro-scale RUC and a meso-scale RUC are composed of six kinds of 20-node rectangular isoparametric elements in Figure 1. Based on the MPFE method, the Gauss quadrature formulae are employed to calculate all element integrals in Equations (36)-(39). In addition, 27 Gauss integral points are selected in each kind of element in this work.
where N is the number of 20-node rectangular isoparametric elements of micro-scale RUC or meso-scale RUC, respectively. G α is the material property matrix of components, W ii , W jj , W kk is the corresponding weights. J is the Jacobian matrix.
On the premise of fine mesh generation and convergence of finite element results, the matrix φ and the homogenized effective thermo-viscoelastic properties of the yarns can be obtained firstly. Similarly, once the matrix ϑ has been determined, the homogenized effective thermo-viscoelastic properties of the meso-scale RUC can be obtained.

FsMsFE Analysis Sequence
In the following, a detailed description of the proposed five-step multi-scale strategy for determining the thermo-viscoelastic behaviors of 3D braided composites will be given.
(1) Establish a micro-scale RUC model. Then, Equation (35) is solved to obtain the characteristic function φ. Thereafter, the equivalent thermo-viscoelastic properties α H1 kl and β H1 ij of the yarns are computed by Equations (15) and (19).

Numerical Models and Verification
In this section, the thermo-viscoelastic behaviors of 3D four-directional braided composites are studied based on the FsMsFE method presented in this paper. The reinforced fiber is T300 carbon and the matrix is epoxy resin ED-6, respectively. The material parameters of the carbon fiber and resin are shown in Table 1. In this work, the time-dependent thermo-viscoelastic response of the 3D braided composite clamped-free beam with an instantaneous temperature increment T ε and a certain forced displacement at the free edge are studied. As shown in Figure 3, the dimension of the beam is 90 × 20 × 6 mm 3 , the temperature increment is T ε and the forced displacement at the free edge is d = 0.1 mm.   When the fiber volume fraction of the yarn is 0.28, the microscale RUC consisted of 3024 20-node rectangular isoparametric elements is used to calculate the effective relaxation modulus of the yarn and verify the accuracy of the present FsMsFE method. To obtain the effective relaxation modulus of the yarn by Equation (14), it is necessary to calculate When the fiber volume fraction of the yarn is 0.28, the microscale RUC consisted of 3024 20-node rectangular isoparametric elements is used to calculate the effective relaxation modulus of the yarn and verify the accuracy of the present FsMsFE method. To obtain the effective relaxation modulus of the yarn by Equation (14), it is necessary to calculate the corrector χ via Equation (35). In this work, it is assumed that the reinforcement phase is considered as elastic transversely isotropic material. The volume deformation of the matrix is considered to be linear elastic and the shear deformation conforms to the three-parameter solid model. Therefore, the material property matrix G fi and G ma for reinforcement fiber and matrix can be written as where K, Y(t) are volume modulus and shear modulus, respectively. G 1 , G 2 , η 2 are the viscoelastic parameters of matrix. E 1 , E 3 , γ 21 , γ 13 , G 12 and G 13 are the elastic parameters of reinforcement phase.
In fof a beam under thermo-mec [34] are put forward to solve the micro-scale problem. When the corrector χ based on the microscale RUC is obtained, Figure 4 shows the comparison between equivalent viscoelastic properties of the yarns obtained by FsMsFE method and that evaluated by conventional FE method [35]. It is found that the present results are in good agreement with the existing data in the literature and the differences caused by different discretization methods is kept in a small range.
where K , () Yt are volume modulus and shear modulus, respectively.  In this section, periodicity boundary conditions [34] are put forward to solve the micro-scale problem. When the corrector  based on the microscale RUC is obtained, Figure 4 shows the comparison between equivalent viscoelastic properties of the yarns obtained by FsMsFE method and that evaluated by conventional FE method [35]. It is found that the present results are in good agreement with the existing data in the literature and the differences caused by different discretization methods is kept in a small range.

Effective Thermo-Viscoelastic Properties
As mentioned in Section 2, to obtain thermo-viscoelastic response of 3D four-directional braided composites, the effective thermo-viscoelastic properties of the yarn and meso-scale RUC should be firstly evaluated. Owing to the thermo-viscoelastic properties for composites are temperature dependent, the effective TTRM β H1 ij , β H2 ij and effective TTEC α H1 kl , α H2 11 of different scales at room temperature 25 • C (ground state) are presented in Figures 5 and 6.
In Figure 5a,b, the yarns are considering the unidirectional reinforced composite with a constant fiber volume fraction 0.722, 0.850 and 0.866, the transverse equivalent TTRM β H1 11 increases monotonically while the longitudinal equivalent TTRM β H1 33 decreases monotonically with the increasing of the time, the transverse equivalent TTRM β H1 11 and longitudinal equivalent TTRM β H1 33 show a downward trend with increasing fiber volume fraction. A similar scenario of effective thermal expansion coefficients α H1 kl in both transverse and longitudinal is shown in Figure 5c,d and they will be reached a steady state. As shown in Figure 6, the effective thermo-viscoelastic properties of the meso-scale RUC under different braiding angles (20 • , 30 • , 40 • ) and fiber volume fraction 0.37 are calculated. It can be seen that both the transverse and longitudinal equivalent TTRM β H2 ij are firstly decrease and then increase when the braiding angle increases and the transverse effective TTEC α H2 11 decrease monotonically with the increasing of the braiding angle when the longitudinal effective TTEC α H2 33 firstly increase and then decrease. Additionally, Figures 5 and 6 show that the different fiber fraction and braiding angle in the ground state will affect the thermo-viscoelastic properties of composite in time domain.

Effective Thermo-Viscoelastic Properties
As mentioned in Section 2, to obtain thermo-viscoelastic response of 3D four-directional braided composites, the effective thermo-viscoelastic properties of the yarn and meso-scale RUC should be firstly evaluated. Owing to the thermo-viscoelastic properties for composites are temperature dependent, the effective TTRM  In Figure 5a,b, the yarns are considering the unidirectional reinforced composite with a constant fiber volume fraction 0.722, 0.850 and 0.866, the transverse equivalent TTRM In order to observe the effects of temperature on the thermo-viscoelastic properties of 3D four directional braided composites, the analysis results of yarns (V f = 0.85) and meso-scale RUC 2 (V f = 0.37, ϕ = 30 • ) at different temperatures (25 • C, 35 • C, 45 • C) are shown in Figures 7 and 8. It can be found that the influence of temperature on thermoviscoelastic properties of 3D four directional braided composites is very obvious and higher the temperature the shorter the time required for the effective TTRM and effective TTEC to stabilize. This phenomenon is similar to the effect of temperature on the viscoelastic properties of 3D composites reported by Seifert et al. [21].

Time-Dependent Multiscale Stress Field Distributions
The initial damage of 3D braided composites is commonly attributed to the destruction of different scales, therefore a multiscale stress relaxation analysis is given in this section. Based on the different thermal and mechanical loading boundary conditions, the macro/meso/micro-scale time-dependent stresses can be obtained by using the FsMsFE analysis method presented in the Sections 2.2 and 3.1.
When the braiding parameters are determined (Vf = 0.37, φ = 30°) and the temperature is 45 °C , Figure 9 presents the macro-scale stress

Time-Dependent Multiscale Stress Field Distributions
The initial damage of 3D braided composites is commonly attributed to the destruction of different scales, therefore a multiscale stress relaxation analysis is given in this section. Based on the different thermal and mechanical loading boundary conditions, the macro/meso/micro-scale time-dependent stresses can be obtained by using the FsMsFE analysis method presented in the Sections 2.2 and 3.1.
When the braiding parameters are determined (V f = 0.37, ϕ = 30 • ) and the temperature is 45 • C, Figure 9 presents the macro-scale stress σ z of position A, the meso-scale stress σ z on the plane (Z = 0.5H) of a RUC and the micro-scale stress σ 3 (fiber direction) belong to the region B of the mesoscopic yarn. It can be noted that the stress relaxation response of 3D four-directional braided composites consists of a shorter relaxation stage and a longer steady stage. At the same moment, the smaller the selected scale the greater the stresses.

Time-Dependent Multiscale Stress Field Distributions
The initial damage of 3D braided composites is commonly attributed to the destruction of different scales, therefore a multiscale stress relaxation analysis is given in this section. Based on the different thermal and mechanical loading boundary conditions, the macro/meso/micro-scale time-dependent stresses can be obtained by using the FsMsFE analysis method presented in the Sections 2.2 and 3.1.
When the braiding parameters are determined (Vf = 0.37, φ = 30°) and the temperature is 45 °C , Figure 9 presents the macro-scale stress z  of position A, the meso-scale stress z  on the plane (Z = 0.5H) of a RUC and the micro-scale stress 3  (fiber direction) belong to the region B of the mesoscopic yarn. It can be noted that the stress relaxation response of 3D four-directional braided composites consists of a shorter relaxation stage and a longer steady stage. At the same moment, the smaller the selected scale the greater the stresses.  In order to determine the influence of temperature changes on the multi-scale stress of 3D braided composites, the initial and stable state stress distributions of different scales at different temperatures (25 • C, 35 • C, 45 • C) are presented in Figures 10-12. icantly before and after stabilization under the influence of temperature changes and the tensile load is mainly shared by the braided yarns. As shown in Figure 12, it can be observed that there is a more obvious micro-scale stress change for matrix than the fiber region under the influence of temperature changes. Based on the above phenomenon, it can be noted that the initial damage probability of 3D braided composites will increase and the viscoelastic effect in the micro-scale fiber/matrix model cannot be neglected when the temperature increment is considered.   tensile load is mainly shared by the braided yarns. As shown in Figure 12, it can be observed that there is a more obvious micro-scale stress change for matrix than the fiber region under the influence of temperature changes. Based on the above phenomenon, it can be noted that the initial damage probability of 3D braided composites will increase and the viscoelastic effect in the micro-scale fiber/matrix model cannot be neglected when the temperature increment is considered.

Conclusions
This article presents a novel FsMsFE calculation procedure for simulating the thermo-viscoelastic response of 3D braided four-directional composites. The integration of FE calculation programs is more flexible and the dependence on the microstructure type is less through the developed Fortran calculation procedure. The five-step homogenization is carried out at micro-scale, meso-scale and macro-scale and the equivalent thermo-viscoelastic properties and stress relaxation at three scales for 3D braided composite subject to temperature increment and forced displacement are investigated, respectively. The thermo-viscoelastic properties of the 3D braided composites are obviously anisotropic and time-temperature-dependent. The stress relaxation behavior of 3D braided composites is greatly affected by temperature and temperature increment will shorten the relaxation stage of 3D braided composites and aggravate the material failure. It shows obvious non-uniformity for the stress distribution of micro-scale RUC and mesoscale RUC, especially at the interfaces of the fiber/matrix or yarn/matrix. This method can be used for the structural design of braided composite materials when the viscoelasticity and temperature rise of the material are considered. In Figure 10, it is shown that the macroscopic stress peak value change is small, the temperature rise causes the macroscopic stress stabilization time to be significantly reduced. Simultaneously, it can be observed that the time-dependent stress under the fixed displacement boundary condition at ground state is consistent with the trend described by Tate et al. [16]. From Figure 11, it is found that both meso-scale stress increases significantly before and after stabilization under the influence of temperature changes and the tensile load is mainly shared by the braided yarns. As shown in Figure 12, it can be observed that there is a more obvious micro-scale stress change for matrix than the fiber region under the influence of temperature changes. Based on the above phenomenon, it can be noted that the initial damage probability of 3D braided composites will increase and the viscoelastic effect in the micro-scale fiber/matrix model cannot be neglected when the temperature increment is considered.

Conclusions
This article presents a novel FsMsFE calculation procedure for simulating the thermoviscoelastic response of 3D braided four-directional composites. The integration of FE calculation programs is more flexible and the dependence on the microstructure type is less through the developed Fortran calculation procedure. The five-step homogenization is carried out at micro-scale, meso-scale and macro-scale and the equivalent thermoviscoelastic properties and stress relaxation at three scales for 3D braided composite subject to temperature increment and forced displacement are investigated, respectively. The thermo-viscoelastic properties of the 3D braided composites are obviously anisotropic and time-temperature-dependent. The stress relaxation behavior of 3D braided composites is greatly affected by temperature and temperature increment will shorten the relaxation stage of 3D braided composites and aggravate the material failure. It shows obvious nonuniformity for the stress distribution of micro-scale RUC and mesoscale RUC, especially at the interfaces of the fiber/matrix or yarn/matrix. This method can be used for the structural design of braided composite materials when the viscoelasticity and temperature rise of the material are considered.