An Improved Analytical Solution for Process-Induced Residual Stresses and Deformations in Flat Composite Laminates Considering Thermo-Viscoelastic Effects

Dimensional control can be a major concern in the processing of composite structures. Compared to numerical models based on finite element methods, the analytical method can provide a faster prediction of process-induced residual stresses and deformations with a certain level of accuracy. It can explain the underlying mechanisms. In this paper, an improved analytical solution is proposed to consider thermo-viscoelastic effects on residual stresses and deformations of flat composite laminates during curing. First, an incremental differential equation is derived to describe the viscoelastic behavior of composite materials during curing. Afterward, the analytical solution is developed to solve the differential equation by assuming the solution at the current time, which is a linear combination of the corresponding Laplace equation solutions of all time. Moreover, the analytical solution is extended to investigate cure behavior of multilayer composite laminates during manufacturing. Good agreement between the analytical solution results and the experimental and finite element analysis (FEA) results validates the accuracy and effectiveness of the proposed method. Furthermore, the mechanism generating residual stresses and deformations for unsymmetrical composite laminates is investigated based on the proposed analytical solution.


Introduction
Due to the high specific strength and stiffness, the fiber-reinforced epoxy resin composite materials have been widely applied in many fields including aerospace, automobile, civil infrastructures, and ship industries. However, residual stresses are inevitable in the composite structures during the manufacturing process. These process-induced residual stresses may lead to undesirable shape distortions when the cured components are released from the mold and such distortions are often large enough to render the part unserviceable. Moreover, residual stress during fabrication greatly decreases the fatigue life and dimensional accuracy of composite parts [1][2][3]. It is important to accurately predict the development of residual stresses and deformations during curing [4][5][6].
Many factors may govern the evolution of the residual stresses and deformation of composite structures during the fabrication process. The main factors include tool-part interaction, the mismatch of thermal and chemical properties of constitutive materials, and material degradation or viscoelastic effects during curing [1,7,8]. Numerous investigations regarding the effects of these previously mentioned factors on the process-induced residual stresses and deformations have been implemented by using analytical models and numerical methods based on finite element analysis [9][10][11][12][13].
The analytical method can provide quick prediction of residual stresses and deformations with a certain level of accuracy and explain the underlying mechanisms. Furthermore, the analytical method can also be used to validate the FE results for the simple cases where analytical solution can provide accurate results. Consequently, the analytical method has been increasingly applied in the investigations of cure behavior of composite structures during manufacturing [14][15][16][17][18][19][20][21][22][23]. Nelson and Cairns proposed a simple equation to predict the spring-in of the curved plate during the cooling-down stage based on thermo-elastic theory [14]. Radford et al. [15] modified this analytical model by taking into account the effects of chemical shrinkage of the resin on the spring-in. Moreover, a modified analytical model has been proposed to consider the effects of part thickness on the spring-in of C-shaped curved composite parts based on path-dependent constitutive models by Wisnom et al. [16]. The tool-part interactions have important influence on the residual stresses and deformations of composite structures, which have been considerably investigated by the experimental and finite element methods [4,7,24,25]. However, the effects of tool-part interactions on the residual stresses and deformations of composite structures have not been taken into account in the previously mentioned analytical models [4][5][6][7][8][9][10][11][12][13][14][15][16]. Consequently, a simple analytical model has been proposed to regard effects of the tool-part interactions on the cure-induced residual stresses and deformation of flat composite laminates by Twigg et al. [17]. The deformations of flat composite laminates obtained by this presented analytical solution align well with the experimental results. However, the dependence of residual stresses on material properties has not been captured by this model. Consequently, Arafath A et al. [19][20][21] presented a set of more accurate analytical models to predict the effects of tool-part interactions on the distribution of residual stresses and deformations for flat and curved composite laminates. The predictions of residual stresses and deformations obtained by these previously mentioned analytical models agree well with the numerical results and measured results. However, the viscoelastic effects of resin on the residual stresses and deformations of composite laminates during curing have not been taken into consideration in any of these previously mentioned models.
In many recent studies, the viscoelastic effects of the resin regarding the cure-behavior of composite laminates have been considered through finite element analysis [26][27][28][29][30][31][32][33][34][35]. However, the finite element analysis is usually very time-consuming and it is difficult to gain any insight into the basic mechanism generating the residual stresses and distortions of composite parts during curing. In this paper, an improved analytical model is proposed by taking into account the viscoelastic effects of epoxy resin for the development of residual stresses and distortions of flat composite laminate during the manufacturing process. In the proposed model, the time-temperature superposition thermo-viscoelastic theory is introduced in this paper to include the viscoelastic effects of resin during the curing process. Moreover, the Timoshenko's beam theory is adopted to describe the mechanical behavior of the composite layers or the tool-part interactions. An incremental formulation of the governing equation for the process-induced residual stresses and distortions with the viscoelastic effects is then established. Furthermore, a novel iterative strategy is proposed to solve the governing equation. To validate the accuracy and efficiency, a set of simulations are implemented based on the proposed model and finite element analysis. The corresponding FEA code was developed and incorporated into ANSYS with a User MAT subroutine. The simulation results obtained by the proposed model are compared to both the finite element results and experimental results.

Theoretical Development of the Analytical Solution with Thermo-Viscoelastic Effects
In this section, the proposed analytical solution with thermo-viscoelastic effects on the development of process-induced residual stresses and deformations is formulated in detail. An incremental formulation of the governing differential equation is first derived based on the Timoshenko's beam theory and viscoelastic constitutive equation to describe the viscoelastic mechanical behavior of a unidirectional composite laminate during the curing process. Second, the analytical solution of the governing differential equation is developed. Third, the analytical solution of the residual stresses and deformations for multilayer composite laminates is derived based on the proposed analytical solution for unidirectional composite laminate. At last, the method to calculate the deformation of the composite laminates is given.

The Governing Differential Equation
The proposed analytical solution is based on the further improvement of the analytical solution proposed by Arafath A et al. [19][20][21]. Consequently, the same simplifications are adopted, which are described below: (1) The Poisson' effects are neglected.
(2) At the symmetrical y-z plane of the composite laminates shown in Figure 1a, the axial displacement equals zero. (3) At the end of the composite laminates in the axial direction, the axial stresses are assumed to be zero.
As shown in Figure 1a, a composite ply would be stretched due to the tool expansion or the interactions of adjacent composite plies, which is a classical beam undergoing tractions at the top and bottom surfaces shown in Figure 1b. According to the Timoshenko's beam theory [36], the equilibrium equation in the x direction in the absence of body forces at the time k∆t during the curing process can be given below: where ∆σ k xx and ∆τ k xy are the incremental axial stress and shear stress at time k∆t, respectively. analytical solution of the governing differential equation is developed. Third, the analytical solution of the residual stresses and deformations for multilayer composite laminates is derived based on the proposed analytical solution for unidirectional composite laminate. At last, the method to calculate the deformation of the composite laminates is given.

The Governing Differential Equation
The proposed analytical solution is based on the further improvement of the analytical solution proposed by Arafath A et al. [19][20][21]. Consequently, the same simplifications are adopted, which are described below: (1) The Poisson' effects are neglected.
(2) At the symmetrical y-z plane of the composite laminates shown in Figure 1a, the axial displacement equals zero. (3) At the end of the composite laminates in the axial direction, the axial stresses are assumed to be zero.
As shown in Figure 1a, a composite ply would be stretched due to the tool expansion or the interactions of adjacent composite plies, which is a classical beam undergoing tractions at the top and bottom surfaces shown in Figure 1b. According to the Timoshenko's beam theory [36], the equilibrium equation in the x direction in the absence of body forces at the time k t Δ during the curing process can be given below: Δσ and k xy Δτ are the incremental axial stress and shear stress at time kΔt , respectively. According to the integral form of polymer viscoelastic behavior and characteristics of thermo-rheologically simple materials, both the relaxed function and the constitutive relationship can be simplified to decouple the temperature effect using the time-temperature superposition (TTS) principle. A three dimensional (3D) incremental viscoelastic constitutive equation can be written as the equation below [27][28][29][30][31][32][33][34][35].  According to the integral form of polymer viscoelastic behavior and characteristics of thermo-rheologically simple materials, both the relaxed function and the constitutive relationship can be simplified to decouple the temperature effect using the time-temperature superposition (TTS) principle. A three dimensional (3D) incremental viscoelastic constitutive equation can be written as the equation below [27][28][29][30][31][32][33][34][35].
where ∆σ k and σ k−1 are the incremental stress vector and stress vector at time k∆t and (k − 1)∆t, respectively. w is the total number of the Maxwell elements. i is the serial number of each Maxwell element. C * is the time-dependent consistent tangent stiffness matrix at time k∆t, which can be defined by the equation below. where C u and C ∞ are the fully unrelaxed and fully relaxed stiffness matrix, respectively. ω i is the weigh factors for the ith Maxwell element. The stiffness A k 1i coefficient and stress relaxation coefficient A k 2i at time k∆t can be obtained as: where a T denotes the shift factor as the function of degree of cure and temperature. τ i is the discrete stress relaxation time for the ith Maxwell element. ∆ξ k is the incremental reduced time at current time k∆t and can be calculated by ∆ξ k = ∆t/a T . Without the Poisson's effects, the incremental viscoelastic constitutive relationship between stress and strain at time k∆t can be written by the formulas below.
where σ k−1 xyi and τ k−1 xyi are the stress and shear stress of the ith Maxwell element at time (k − 1)∆t, respectively. ∆u k x and ∆u k y are the first derivatives of axial mechanical displacement increment k∆t at time t with respect to the x and y, respectively. E k xx and G k xy denote the time-dependent consistent tangent axial and shear modulus at time k∆t, respectively. They can be obtained by using the equation below.
where E ∞ xx and G ∞ xy are the fully relaxed axial and shear modulus, respectively. E u xx and G u xy are the fully unrelaxed axial and shear modulus, respectively.
Substituting Equations (5) and (6) into Equation (1) and, after some rearrangements, the governing equation can be written by using the equation below.
where ∆u k xx and ∆u k yy are the second derivative of axial displacement increment ∆u k at time k∆t with respect to the x and y, respectively, and the parameters ∆1, ∆2 can be obtained by using the equations below.
where E xxi and G xyi are the unrelaxed axial and shear modulus of the ith Maxwell element, respectively.
They can be calculated by E xxi = ω i (E u xx − E ∞ xx ) and G xyi = ω i G u xy − G ∞ xy , respectively.

The Solution of the Governing Differential Equation
It can be seen that the governing Equation (7) is a second order non-homogeneous linear ordinary differential equation with constant coefficients. Usually, it is difficult to directly obtain the closed form solution for the governing Equation (7). However, there would be little viscoelastic effects when the material is in the fully-unrelaxed or fully relaxed state where both A k−1 1i and A k 2i equal to one or zero. Therefore, the previously mentioned parameters ∆1 and ∆2 would be equal to zero and Equation (7) would become the Laplace equation, which can be written by using the equation below.
E k xx ∆u k xx + G k xy ∆u k yy = 0 (9) The Laplace Equation (9) can be directly solved by the method of separation of variables [19][20][21]. The solution of Equation (8) can be given as the equation below.
sin(k n x) A k n exp β k n y + B k n exp −β k n y + ∆ε k the-chem x (10) where A k n and B k n are the unknown constants and can be found from the boundary conditions at the interfaces. ∆ε k the-chem is the axial free thermal-chemical strain increment at time k∆t and can be calculated by ∆ε kj the-chem = ∆Tα CTE + ∆Cγ CSE .α CTE and γ CSE are the coefficient of thermal expansion and the coefficient of shrinkage expansion, respectively. ∆T and ∆C are the temperature increment and cure degree increment. The characteristic constant β k n is time-dependent and can be calculated by the equation below.
where l is half of the length of the laminate in the x direction. It can be seen from Equation (7) that the axial displacement increment ∆u k at time k∆t is related to the axial displacement increment of all the previous (k − 1) time steps. In order to obtain the solution of Equation (7), the time j∆t is defined by assuming that the material transforms from the fully relaxed state into a viscoelastic solid state. Readers can refer to Reference [34] for the method to calculate the time j∆t. For the fully relaxed state, the mechanical behavior of material can be described by the Laplace Equation (9) so that the solution at the time (j − 1)∆t still conforms to the formulation as Equation (10). According to superposition principle of a differential equation with constant coefficients, it can be inferred that the solution of Equation (9) at time j∆t can be a linear combination of solutions of the corresponding Laplace equations at time (j − 1)∆t and j∆t, which can be written as: To keep the balance on both sides of the equation, the non-homogeneous solution of equation ∆u k must contain all the homogeneous solutions ∆u(q∆t) of all the previous time steps. Since the parameters ∆1 and ∆2 are a linear combination of ∆u q , it can be inferred that the non-homogeneous solution of equation ∆u k at time k∆t can be written as a linear combination of the homogeneous solution ∆u(q∆t) of all time.
Substituting Equation (13) into Equation (7) yields the coefficient a k,q below.
It can be seen that the coefficient a k,q is independent of the unknowns A k n and B k n and can be directly solved based on the calculation of material parameters of composite laminates during the curing process. The detailed derivation of the coefficient a k,q is shown in the Appendix A.

The Solution of Residual Stresses for Multi-Layers Composite Laminates
It is well known that the composite laminates are generally composed of laminate with different fiber orientations cured on a solid tool in autoclave, which is shown in Figure 2. To predict accurately the cure-induced residual stresses of multi-layers composite laminates, both the interactions of different composite layers and the tool-part interactions are considered in this scenario.
where q c is a time-dependent characteristic parameter and can be obtained by q q q x x x y c = E / G . The coefficient k-1,q i d can be iteratively solved using the equation below.
It can be seen that the coefficient k,q a is independent of the unknowns k n A and k n B and can be directly solved based on the calculation of material parameters of composite laminates during the curing process. The detailed derivation of the coefficient k,q a is shown in the Appendix.

The Solution of Residual Stresses for Multi-Layers Composite Laminates
It is well known that the composite laminates are generally composed of laminate with different fiber orientations cured on a solid tool in autoclave, which is shown in Figure 2. To predict accurately the cure-induced residual stresses of multi-layers composite laminates, both the interactions of different composite layers and the tool-part interactions are considered in this scenario.
where ∆u j (q∆t) is the corresponding Laplace solution of the jth layer at time q∆t, which can be calculated by using the equation below.
where A q,j n and B q,j n are unknowns associated with the jth layer. For the tool: As the tool is usually made of metallic materials, the relationship between the stress and strain can continually be described by the elastic theory and written as Equation (18) below [19][20][21].
where G 0 and t 0 are the shear modulus and thickness of tool material. Boundary conditions According to the layerwise approach [37], the layer axial displacement and shear stress should be equal on the interface between the two neighbor plies, which means that, at a generic jth interface (y = y j ), there are: It can be seen that Equation (18) constitutes 2(m − 1) equations for (m − 1) interfaces. For a stress free boundary condition at the top surface, the shear stress ∆τ k,m xy can be expressed as: Polymer film and a release agent from the mold usually separate the composite parts during autoclave processing. A compliant interface layer with shear modulus G s and thickness t s has been proposed to describe the tool-part interactions during curing by Arafath et al. [19][20][21]. The interface layer is very thin and the only stress that it transfers is a shear stress ∆τ s at time k∆t across its thickness, which is adopted in this scenario and can be formulated as: where ∆u k,1 and ∆u k,0 are the longitudinal displacements of the part and tool at the position y = 0. As Equation (19)-(21) contains 2m + 1 equations, the total 2m + 1 unknowns including the unknowns A q,j n and B q,j n of each composite layer and the unknowns D k 1n of the tool can be determined and then the stresses and displacement of each layer can be calculated.

Deformations Prediction
At the end of the curing process, the tool is removed and the deformations are calculated by the unbalanced moment due to axial stress distribution through thickness. The moment can be obtained by using the equations below.
where b is the beam width and h is the thickness of composite laminates. Then, the bending deformation can be calculated by the equation below.
where (EI) eff is the bending rigidity of the composite laminate at the end of the cure cycle.
In summary, the process-induced residual stresses and deformations of composite components during curing can be solved by the flow chart shown in Figure 3.
where ( ) eff EI is the bending rigidity of the composite laminate at the end of the cure cycle. In summary, the process-induced residual stresses and deformations of composite components during curing can be solved by the flow chart shown in Figure 3.

Material Properties and MODEL
Both numerical simulation and experimental verification for the symmetrical and unsymmetrical composite laminates [0 2 /90 2 ]s and [0 4 /90 4 ] would be carried out in this paper to verify the efficiency and accuracy of the proposed analytical model and compared with previous studies. The illustrations of the materials and model are described as follows.

Material Properties
The composite parts used in this study are made of graphite/epoxy (AS4/3501-06) unidirectional prepregs from Hexcel Composites, which has been used for thermo-viscoelastic analysis during the curing process in previous studies [25][26][27][28][29][30][31][32][33][34][35]. The viscoelastic material properties are listed in Tables 1 and 2, respectively. It is assumed that the initial resin modulus in the viscous state is 5-6 orders of magnitude less than that of the fully cured resin, which has been verified by experimental results [19][20][21][22][23] and is adopted in this case.
The composite specimen strips are generally manufactured by the autoclave assisting cure method, which is shown in Figure 4. It can be seen that a polymer release film and release agent separates the tool and the composite part during the curing process in the autoclave. During the curing process, the tool and composite part will interact with each other due to the mismatch of thermal-chemical properties between different materials. To improve the efficiency of modeling, the release film and the agent between the tool surface and the composite have been modeled as a shear layer, which is very thin, and the only stress that it transfers is a shear stress that is uniform across its thickness. The material properties of the tool and the shear layer are shown in Table 3 [19][20][21].    The simulations are implemented based on both the proposed new analytical solution and the viscoelastic finite element analysis. The linear element solid 186 with 20 nodes is used in the FEA analysis in which the three-dimensional thermo-viscoelastic constitutive Equation (2) is incorporated into finite element software ANSYS as a user subroutine USERMAT to describe viscoelastic behavior of materials. Considering the symmetrical geometry and boundary conditions of the composite laminate, only 1/4 part of the specimen is modeled. The symmetrical constraints are applied in the X-Z and Y-Z symmetrical planes. The longitude and transverse directions are defined as X and Y, respectively. The thickness direction is defined as Z. Each ply contains 2 elements and 40 elements along x directions for the composite laminate. The shear layer ply is meshed into 2 elements through the thickness direction. The total number of nodes and elements are 9608 and 1320, respectively. The chemical shrinkage of the resin and cure thermal loads are treated as equivalent thermal loads and applied in the FEA model. The finite element analysis model is shown as Figure 5. In the finite element analysis by ANSYS, the cycle solution strategy is applied to solve the cure-induced residual

Finite Element Model
The simulations are implemented based on both the proposed new analytical solution and the viscoelastic finite element analysis. The linear element solid 186 with 20 nodes is used in the FEA analysis in which the three-dimensional thermo-viscoelastic constitutive Equation (2) is incorporated into finite element software ANSYS as a user subroutine USERMAT to describe viscoelastic behavior of materials. Considering the symmetrical geometry and boundary conditions of the composite laminate, only 1/4 part of the specimen is modeled. The symmetrical constraints are applied in the X-Z and Y-Z symmetrical planes. The longitude and transverse directions are defined as X and Y, respectively. The thickness direction is defined as Z. Each ply contains 2 elements and 40 elements along x directions for the composite laminate. The shear layer ply is meshed into 2 elements through the thickness direction. The total number of nodes and elements are 9608 and 1320, respectively. The chemical shrinkage of the resin and cure thermal loads are treated as equivalent thermal loads and applied in the FEA model. The finite element analysis model is shown as Figure 5. In the finite element analysis by ANSYS, the cycle solution strategy is applied to solve the cure-induced residual stresses and deformations of the composite structure.
The thermo-chemical analyses should be developed to obtain the temperature and degree of cure distributions throughout the laminates [4,35]. However, the thickness of composite laminates used in this case is 1 mm, which is so thin that the cure degree and temperature field can be considered uniform. Thus, thermo-chemical analysis is not presented in this scenario.
viscoelastic behavior of materials. Considering the symmetrical geometry and boundary conditions of the composite laminate, only 1/4 part of the specimen is modeled. The symmetrical constraints are applied in the X-Z and Y-Z symmetrical planes. The longitude and transverse directions are defined as X and Y, respectively. The thickness direction is defined as Z. Each ply contains 2 elements and 40 elements along x directions for the composite laminate. The shear layer ply is meshed into 2 elements through the thickness direction. The total number of nodes and elements are 9608 and 1320, respectively. The chemical shrinkage of the resin and cure thermal loads are treated as equivalent thermal loads and applied in the FEA model. The finite element analysis model is shown as Figure 5. In the finite element analysis by ANSYS, the cycle solution strategy is applied to solve the cure-induced residual stresses and deformations of the composite structure.
The thermo-chemical analyses should be developed to obtain the temperature and degree of cure distributions throughout the laminates [4,35]. However, the thickness of composite laminates used in this case is 1 mm, which is so thin that the cure degree and temperature field can be considered uniform. Thus, thermo-chemical analysis is not presented in this scenario.

Case 1: Numerical Simulation for Symmetrical Composite Laminates
In this section, numerical simulation for the process-induced stresses and deformations of composite laminates with [02/902]s layups cured on a solid tool is implemented based on the analytical solutions and finite element analysis. The recommended manufacturing cure cycle presented in Reference [35] in the autoclave consists of a first ramp of 2.5 °C/min up to 116 °C hold

Case 1: Numerical Simulation for Symmetrical Composite Laminates
In this section, numerical simulation for the process-induced stresses and deformations of composite laminates with [0 2 /90 2 ]s layups cured on a solid tool is implemented based on the analytical solutions and finite element analysis. The recommended manufacturing cure cycle presented in Reference [35] in the autoclave consists of a first ramp of 2.5 • C/min up to 116 • C hold for 60 min, a second ramp of 2.5 • C/min up to 175 • C, which holds for 20 min, and a cool down of −2.5 • C /min up to 25 • C. Figure 6 shows the axial stress distribution obtained by the viscoelastic finite element analysis. Figure 7 shows the comparison of the numerical predictions and analytical predictions of axial stress distribution through-thickness direction. The analytical predictions are both implemented by the proposed analytical solution and the analytical solution presented by Arafath A et al. [19][20][21], respectively. It can be seen that the results obtained by the proposed model show good agreement with the results of viscoelastic finite element analysis. The maximum difference of axial stresses between the proposed analytical method and the viscoelastic finite element analysis is 1.3 Mpa, which corresponds to a relative error rate of 3.8%. On the contrary, it can be observed that there are large differences between the residual stresses obtained by the analytical solution presented by Arafath A et al. [19][20][21] and the proposed analytical solution. The maximum difference of axial stresses is 11.8 Mpa, which corresponds to a relative error rate of 40.3%.  Figure 6 shows the axial stress distribution obtained by the viscoelastic finite element analysis. Figure 7 shows the comparison of the numerical predictions and analytical predictions of axial stress distribution through-thickness direction. The analytical predictions are both implemented by the proposed analytical solution and the analytical solution presented by Arafath A et al. [19][20][21], respectively. It can be seen that the results obtained by the proposed model show good agreement with the results of viscoelastic finite element analysis. The maximum difference of axial stresses between the proposed analytical method and the viscoelastic finite element analysis is 1.3 Mpa, which corresponds to a relative error rate of 3.8%. On the contrary, it can be observed that there are large differences between the residual stresses obtained by the analytical solution presented by Arafath A et al. [19][20][21] and the proposed analytical solution. The maximum difference of axial stresses is 11.8 Mpa, which corresponds to a relative error rate of 40.3%.
These simulation results show that the proposed analytical solution can deal well with the viscoelastic effects of the resin on the residual stresses of symmetrical composite laminates and is more accurate than the analytical solution presented by Arafath A et al. [19][20][21].      Figure 9 shows the prediction of the deformations of the composite laminates at the end of the cure process, which are also calculated by the previously mentioned methods. It can be observed that the results obtained by the proposed analytical solution and the analytical solution by Arafath A et al. [19][20][21] are in good agreement with the results obtained by viscoelastic finite element analysis. Moreover, it can be seen that the residual stresses change sharply from compression to tension in the 0°-90° interface. This can be explained by the composite laminates contracting due to chemical shrinkage These simulation results show that the proposed analytical solution can deal well with the viscoelastic effects of the resin on the residual stresses of symmetrical composite laminates and is more accurate than the analytical solution presented by Arafath A et al. [19][20][21]. Figure 8 shows the prediction of the deformations obtained by the viscoelastic FEA. Figure 9 shows the prediction of the deformations of the composite laminates at the end of the cure process, which are also calculated by the previously mentioned methods. It can be observed that the results obtained by the proposed analytical solution and the analytical solution by Arafath A et al. [19][20][21] are in good agreement with the results obtained by viscoelastic finite element analysis. Moreover, it can be seen that the residual stresses change sharply from compression to tension in the 0 • -90 • interface. This can be explained by the composite laminates contracting due to chemical shrinkage of the resin and the cooling process when curing. As both the thermal and chemical shrinkage coefficients of 90 • composite layers are larger than those of 0 • composite layers, the 0 • layers would be compressed and, thus, the 90 • layers would be tensioned to keep balance. These results indicate that the proposed analytical solution can accurately predict the cure-induced deformations of symmetrical composite laminates compared with Arafath A et al. [19][20][21]. of the resin and the cooling process when curing. As both the thermal and chemical shrinkage coefficients of 90° composite layers are larger than those of 0° composite layers, the 0° layers would be compressed and, thus, the 90° layers would be tensioned to keep balance. These results indicate that the proposed analytical solution can accurately predict the cure-induced deformations of symmetrical composite laminates compared with Arafath A et al. [19][20][21].     Figure 10 shows the comparison of development of the resultant bending moment during the curing process between the proposed analytical solution and the analytical solution presented by Arafath A et al. [19][20][21]. It can be seen that the bending moment mainly occurs at the initial stages of the curing process when the epoxy resin is in the viscous state and shear modulus is very low. In a viscous state, the matrix is in the fully-relaxed state so that there are no obvious viscoelastic effects. The proposed analytical solution would degenerate into the elastic solution, which is the same as the analytical solution proposed by Arafath A et al. [19][20][21]. Therefore, the bending moment obtained by the proposed model is the same as that obtained by the analytical model proposed by Arafath A et al. [19][20][21] in the viscous state before the gelation point. Since the stresses produced in the viscoelastic state are basically symmetrical for the x-z plane, there is almost no bending moment. Therefore, the deformation obtained by the proposed analytical solution agrees with that of the analytical solution presented by Arafath A et al. in the viscoelastic state after gelation.
At length, the efficiency of the proposed analytical solution in calculating the residual stress and deformations is compared to that of the viscoelastic finite element analysis and the analytical solution presented by Arafath A et al. [19][20][21]. It can be seen from Table 4 that, although the run-time  Figure 10 shows the comparison of development of the resultant bending moment during the curing process between the proposed analytical solution and the analytical solution presented by Arafath A et al. [19][20][21]. It can be seen that the bending moment mainly occurs at the initial stages of the curing process when the epoxy resin is in the viscous state and shear modulus is very low. In a viscous state, the matrix is in the fully-relaxed state so that there are no obvious viscoelastic effects. The proposed analytical solution would degenerate into the elastic solution, which is the same as the analytical solution proposed by Arafath A et al. [19][20][21]. Therefore, the bending moment obtained by the proposed model is the same as that obtained by the analytical model proposed by Arafath A et al. [19][20][21] in the viscous state before the gelation point. Since the stresses produced in the viscoelastic state are basically symmetrical for the x-z plane, there is almost no bending moment. Therefore, the deformation obtained by the proposed analytical solution agrees with that of the analytical solution presented by Arafath A et al. in the viscoelastic state after gelation.
At length, the efficiency of the proposed analytical solution in calculating the residual stress and deformations is compared to that of the viscoelastic finite element analysis and the analytical solution presented by Arafath A et al. [19][20][21]. It can be seen from Table 4 that, although the run-time of the proposed analytical model is slightly greater than the model proposed by Arafath A et al. [19][20][21], the proposed analytical model is still very efficient when compared to the FEA model. The run-time of the proposed analytical solution is less than 20 s in this case. of the proposed analytical model is slightly greater than the model proposed by Arafath A et al. [19][20][21], the proposed analytical model is still very efficient when compared to the FEA model. The run-time of the proposed analytical solution is less than 20 s in this case.  The analytical solution by Arafath A et al [19][20][21] The proposed analytical solution Figure 10. Comparison of the development of the bending moment by different analytical models of symmetrical composite laminates. These previously mentioned results show that the proposed analytical solution is more accurate and efficient to predict the process-induced residual stresses and deformation of the symmetrical composite laminates compared to the solution presented by Arafath A et al. [19][20][21].

Illustration of the Model
Compared to the investigations for the cure behavior of symmetrical composite laminates, there is limited knowledge on the development of residual stresses and deformation of the unsymmetrical composite laminates during curing. However, it does not mean that coupled unsymmetrical laminates are useless. They can be used during unique conditions [38][39][40][41]. Considering that this paper focuses on the viscoelastic effects of the matrix on the process-induced residual stresses and distortions, a set of specimen strips with [0 4 /90 4 ] layup composite laminates, which has been experimentally investigated by Kim and Hahn in Reference [38] are taken as an example for simulations to further verify the accuracy and effectiveness of the proposed analytical model. In the study mentioned earlier [38], a set of unsymmetrical composite laminates made of graphite/epoxy (T300/3501-6) prepregs have been subjected to interrupted cure cycles shown in Figure 11. The heating and cooling rates are set to be approximately 3 • C/min. The unsymmetrical composite laminates are 152 mm long strips with a width of 25 mm. The assumption that materials properties of T300/3501-6 prepregs are the same as those of AS4/3501-6 prepregs in Reference [38] is adopted in the following investigations. of T300/3501-6 prepregs are the same as those of AS4/3501-6 prepregs in Reference [38] is adopted in the following investigations.  [19][20][21]. Figure 13 shows the deformations of unsymmetrical composite laminates obtained by the   Figure 13 shows the deformations of unsymmetrical composite laminates obtained by the viscoelastic FEA. The predicted dimensionless curvatures under the interrupted cure cycles A-E obtained by the proposed analytical solution are shown in Figure 14. The results are obtained by the analytical solution presented by Arafath A et al. [19][20][21] and the viscoelastic finite element analysis as well as the experimental data presented by Kim et al. [38], respectively. It can be seen that the predicted dimensionless curvatures under all cure cycles by the proposed analytical solution are in good agreement with both the experimental results and the results obtained by FEA with the viscoelastic constitutive law. However, all of the predicted values under all cure cycles obtained by the analytical solution presented by Arafath A et al. [19][20][21] are much larger than the predicted results of the proposed analytical solution and experimental values.
These previously mentioned results show that the proposed analytical solution can deal well with the viscoelastic effects of the polymer matrix and provide a more accurate prediction of the development of residual stresses and deformation of the unsymmetrical composite laminates during curing than the analytical solution presented by Arafath A et al. [19][20][21].

The Further Investigations on the Cure Mechanism
It is well known that the residual stresses and deformation of composite laminates are influenced by many factors including tool-part interactions, thermal effects, chemical shrinkage of the resin, and stress relaxation of matrix. However, the mechanism concerning the effects of these above factors on the development of residual stresses and deformation for the unsymmetrical composite laminates during curing has not been clearly clarified. Therefore, a further investigation is implemented to gain insight into the generation mechanism of the residual stresses and deformation based on the proposed analytical solution.

The Further Investigations on the Cure Mechanism
It is well known that the residual stresses and deformation of composite laminates are influenced by many factors including tool-part interactions, thermal effects, chemical shrinkage of the resin, and stress relaxation of matrix. However, the mechanism concerning the effects of these above factors on the development of residual stresses and deformation for the unsymmetrical composite laminates during curing has not been clearly clarified. Therefore, a further investigation is implemented to gain insight into the generation mechanism of the residual stresses and deformation based on the proposed analytical solution.

The Further Investigations on the Cure Mechanism
It is well known that the residual stresses and deformation of composite laminates are influenced by many factors including tool-part interactions, thermal effects, chemical shrinkage of the resin, and stress relaxation of matrix. However, the mechanism concerning the effects of these above factors on the development of residual stresses and deformation for the unsymmetrical composite laminates during curing has not been clearly clarified. Therefore, a further investigation is implemented to gain insight into the generation mechanism of the residual stresses and deformation based on the proposed analytical solution.

Effects of Chemical Shrinkage
The dimensionless curvatures of unsymmetrical composite laminates under the cure cycles A-E are calculated based on the proposed analytical solution and are shown in Figure 15 in which both the cases with thermal effects alone and chemical shrinkage of the resin are implemented. It can be seen that the effects of chemical shrinkage of the resin on the curvature are dependent on the degree of cure and cure cycles.
Moreover, it can be seen from Figure 15 that the chemical shrinkage of the resin has no significant influences on the deformation of unsymmetrical composite laminates under cure cycle A and B. Moreover, there is a better agreement between measured and predicted values obtained by the proposed model considering the chemical shrinkage of the resin when compared to the cases with thermal effects alone after gelation. For the results obtained by the proposed model under cure cycle C, it can be seen that the chemical shrinkage of the resin accounts for about 8% of the final curvature. For the results obtained by the proposed model under cure cycles D and E, it can be seen that the chemical shrinkage of the resin accounts for about 17% of the final curvature. It is very close to the experimental results 20% for a non-symmetrical [0/90] lay-up composite laminate made of AS4/8852 prepregs by R. Akkerman et al. [39].  [39].
The bending moment at the end of curing determines the deformation of the composite laminates. To further investigate the mechanism of cured-deformation, the developments of bending moments under the cure cycles numbered A-E are calculated based on the proposed analytical solution and the solution presented by Arafath A et al. [19][20][21] shown in Figure 16 in which both cases with thermal effects alone and chemical shrinkage of the resin are implemented.
It can be seen from Figure 16a and Figure 16b that there is no significant difference regarding the developments of the bending moment between the case with thermal effects alone and that with chemical shrinkage of the resin under the cure cycle A and B. This is similar to the results of dimensionless curvature shown in Figure 15. It can be explained by the fact that the modulus of the matrix is very low in the viscous state, which means no significant stress can be developed in the composite laminates. These negative bending moments can be caused by the shear interactions between tool and composite parts. In addition, it can be observed from Figure 16a and Figure 16b that there are some positive bending moments during the end stage of the cool-down process, which can be considered to be caused by thermal mismatch of different composite layers. Both the degree of cure is very low so that the relaxation time of the resin system is very low. Therefore, only some small residual stresses have been formed at the end of the cooling process due to interactions of different composite layers. Furthermore, the degree of cure stops increasing at this stage and, therefore, there is no effect of the chemical shrinkage of the resin on the residual stresses and the bending moment under cure cycle A and B.  Figure 15 are, thus, smaller than those under cure cycle D and E shown in Figure 15, respectively. In addition, it can be seen from Figures 16d and 16e that the effects of chemical shrinkage of the resin on the evolution of the bending moments begins from gelation to the end of the second dwelling under cure cycle D and E. The matrix has been almost fully cured at the end of cure cycles D and E so that the bending moments at the end of cure cycle D are nearly equal to those at the end of cure cycle E. Therefore, the effects of the chemical shrinkage of the final curvature under cure cycle D in Figure  15 are nearly the same as those under cure cycle E in Figure 15. Experimental data by Kim [38]  The bending moment at the end of curing determines the deformation of the composite laminates. To further investigate the mechanism of cured-deformation, the developments of bending moments under the cure cycles numbered A-E are calculated based on the proposed analytical solution and the solution presented by Arafath A et al. [19][20][21] shown in Figure 16 in which both cases with thermal effects alone and chemical shrinkage of the resin are implemented.
It can be seen from Figure 16a,b that there is no significant difference regarding the developments of the bending moment between the case with thermal effects alone and that with chemical shrinkage of the resin under the cure cycle A and B. This is similar to the results of dimensionless curvature shown in Figure 15. It can be explained by the fact that the modulus of the matrix is very low in the viscous state, which means no significant stress can be developed in the composite laminates. These negative bending moments can be caused by the shear interactions between tool and composite parts. In addition, it can be observed from Figure 16a,b that there are some positive bending moments during the end stage of the cool-down process, which can be considered to be caused by thermal mismatch of different composite layers. Both the degree of cure is very low so that the relaxation time of the resin system is very low. Therefore, only some small residual stresses have been formed at the end of the cooling process due to interactions of different composite layers. Furthermore, the degree of cure stops increasing at this stage and, therefore, there is no effect of the chemical shrinkage of the resin on the residual stresses and the bending moment under cure cycle A and B.
Moreover, it can be seen from Figure 16c-e that the cure shrinkage of the matrix has clear effects on the development of the bending moments after gelation. The degree of cure under the cure cycle is lower than those under the cure cycle D and E. Therefore, the effects of chemical shrinkage on the bending moment under the cure cycle C are smaller than those under cure cycles D and E. The effects of the chemical shrinkage of the final curvature under cure cycle C are lower than those under cure cycle E. The effects of chemical shrinkage on the final curvature under cure cycle C shown in Figure 15 are, thus, smaller than those under cure cycle D and E shown in Figure 15, respectively. In addition, it can be seen from Figure 16d,e that the effects of chemical shrinkage of the resin on the evolution of the bending moments begins from gelation to the end of the second dwelling under cure cycle D and E. The matrix has been almost fully cured at the end of cure cycles D and E so that the bending moments at the end of cure cycle D are nearly equal to those at the end of cure cycle E. Therefore, the effects of the chemical shrinkage of the final curvature under cure cycle D in Figure 15 are nearly the same as those under cure cycle E in Figure 15.

Effects of Tool-Part Interactions
It can also be seen from Figure 16 that the bending moment due to the tool-part interactions is mainly formed at a viscous state before gelation and is very small when compared to those due to different composite layer interactions that occurred at a viscoelastic state after gelation. The tool-part interactions contribute to the bending moment and deformation by about 4.8 % at the end of cure cycle E. However, it can be seen from Figure 12 that there is a significant stress gradient close to the tool-part interface due to tool-part interactions for all the cases under cure cycles A-E. It can be, thus, concluded that the effects of tool-part interactions on the residual stresses cannot be neglected for the unsymmetrical composite laminates.

Effects of Tool-Part Interactions
It can also be seen from Figure 16 that the bending moment due to the tool-part interactions is mainly formed at a viscous state before gelation and is very small when compared to those due to different composite layer interactions that occurred at a viscoelastic state after gelation. The tool-part interactions contribute to the bending moment and deformation by about 4.8 % at the end of cure cycle E. However, it can be seen from Figure 12 that there is a significant stress gradient close to the tool-part interface due to tool-part interactions for all the cases under cure cycles A-E. It can be, thus, concluded that the effects of tool-part interactions on the residual stresses cannot be neglected for the unsymmetrical composite laminates.

Conclusions
In this paper, an improved analytical solution considering the viscoelastic effects on the residual stresses and deformations of flat composite laminates has been established based on the modification of a previous analytical solution. Specifically speaking, the proposed analytical solution is derived based on the Timoshenko's beam theory and the time-temperature superposition thermo-viscoelastic theory. The accuracy and effectiveness of the proposed analytical solution is verified through a comparison of results with both experiment and viscoelastic FEA for symmetrical and unsymmetrical composite laminates. The numerical simulation and experimental investigation results show that many important factors inducing the residual stresses and deformations such as tool-part interaction, thermal mismatch, chemical shrinkage, and viscoelastic effects of the resin can be accurately analyzed by the proposed analytical solution. In the numerical simulations and experimental investigations, some important findings are summarized below.
(1) The proposed analytical solution can deal well with the viscoelastic effects of the polymer matrix on the residual stresses and deformations during the curing process. The proposed analytical solution is more efficient and has a similar accuracy to 3-D viscoelastic finite element analysis. Funding: Support for this study was provided by the National Natural Science Foundation of China (Grant Number: 51475377). In addition, I would like to thank all the colleagues and anonymous reviewers who helped improve the paper.

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

Appendix A
In this section, the detailed derivation of the solution of the coefficient a k,q and d k−1,q i are given. Substituting Equation (13) into Equation (8) and pursuing some rearrangements, the differential Equation (8)  It can be seen from Equation (A3) that the coefficient a k,q is related to the characteristic parameter β q n . To guarantee that the Equation (A3) can be established for any position y, it can be inferred that the corresponding items including β q n should be equal for both sides of Equation (A3). Therefore, the coefficient a k,q can be obtained as follows: It is known that the equivalent axial elastic modulus and shear modulus is kept unchanged in the fully relaxed state, which means G k xy − E k xx /c 2 q = 0. Therefore, it can be inferred that the coefficient a k,q is zero when G k xy − E k xx /c 2 q equals to zero. The coefficient a k,q can, thus, be rewritten as: It can be seen that the time complexity of d k−1,q i shown as Equation (A2) has been greatly reduced when compared to Equation (24). It is meaningful for improving the efficiency of the proposed analytical solution.