Numerical Analysis of Curing Residual Stress and Deformation in Thermosetting Composite Laminates with Comparison between Different Constitutive Models

A multi-physics coupling numerical model of the curing process is proposed for the thermosetting resin composites in this paper, and the modified “cure hardening instantaneously linear elastic (CHILE)” model and viscoelastic model are adopted to forecast residual stress and deformation during the curing process. The thermophysical properties of both models are evolved in line with temperature and degree of cure (DOC). Accordingly, the numerical simulation results are improved to be more accurate. Additionally, the elastic modulus of the materials is calibrated to be equal to the modulus of viscoelastic relaxation by a defined function of time in the CHILE model. Subsequently, this work effectuates the two proposed models in a three-dimensional composite laminate structure. Through comparing the two numerical outcomes, it is customary that the residual stress and deformation acquired by the modified model of CHILE conform to those ones assessed through adopting the viscoelastic model.


Introduction
By virtue of advantaged properties in mechanical aspects, the thermosetting resin composites have been extensively adopted in numerous fields, including civil infrastructure, ship industries, and the aerospace industry. However, during the manufacturing process of the composites, several phenomena will consequently generate the residual stress and deformation of the composite structures [1][2][3], and therefore decrease the mechanical properties taken on by the composites. These phenomena include material relaxation or degradation, chemical shrinkage, thermal expansion, and inherent anisotropy. In this regard, evaluating the residual stress and deformation in the design and manufacturing of the composite materials is imperative. The procedure of trial-and-error method is the industrial method for estimating the partial deformation. This method is costly, time-consuming, and inaccurate. Comparatively, the simulation of finite element analysis (FEA) is a method capable of accurately, simply, and effectively estimating the residual stress and deformation.

Heat Transfer Equation
In the course of cure period of the thermosetting resin composites, the conduction and generation of heat originates from heat being internally nonlinear in the curing of resin. In this regard, sources of internal heat and external heat flux are involved in the model of transient heat transfer. The model of transient heat transfer is based on Fourier's heat conduction equation, and can be described as [1,2]:  where Q is the internal heat source; T is the instant temperature taken on by composite materials at time t, and t refers to the absolute time; k x , k y and k z denote the anisotropic thermal conductivities taken on by composite materials in the directions of x, y, and z; C pc and ρ c refer to the specific heat and density of the composites. The internal heat source Q is the instantaneous heat generated by the cross-linking polymerization of the resin, which is governed by: where H R is the total release heat during curing; dα/dt and α represent the curing rate of the resin matrix and DOC, respectively; v f represents the fiber volume fraction taken on by composite materials; ρ r is the density of the resin matrix.

Chemical Reaction
The standard differential scanning calorimetry (DSC) counts as the common experiment-based method adopted to study the curing-related kinetics taken on by the thermoset resin composites [31,32]. In the course of the curing period, the resin matrix appears as a bulk of complicated chemical reactions with the development of complicated cross-linked structures and a massive release of heat. In this regard, to shed light on the curing course via a rate function, this work employs a phenomenological model relative to the curing-related kinetics. For 3501-6 epoxy resin, the rate of conversion dα/dt is given by the following equations [1,2,33]: where K i (i = 1, 2, 3) indicates the curing rate constants, which can be given by the Arrhenius functions as follows: where A i are the Arrhenius constant; ∆E i are the activation energies, and R is the gas constant. Table 1 lists the parameters of curing-related kinetics of 3501-6 epoxy resin.

Thermo-Physical Properties
The instantaneous heat conduction model and the solidification kinetics model can effectively simulate the thermochemical process during solidification. This means that the thermo-physical properties of composites and their composition can accurately predict the solidification behaviors of composite materials which mainly contain thermal conductivity, specific heat, and density. Evidently, this temperature and the change in the thermal physical properties of DOC are considered. The thermal properties of AS4 carbon fiber and 3501-6 epoxy resin are listed in Table 2. The thermal physical properties of composites are able to be obtained through adopting the mixture rule in Equations (A1)-(A4) proposed in Appendix A. Viscoelastic model. Currently, the viscoelastic model, even if there are some disadvantages (such as running time consumption and complex formulas) [26][27][28]31,33], remains the most accurate method to forecast the viscoelastic behavior of resins and thermosetting materials. This model factors in the time, temperature, and the change of DOC, a Boltzmann superposition principle. Under the method of time-temperature superposition, the stress of viscoelastic materials is represented by the constitutive function below, in the form of genetic integral [26,27,31]: where the free thermo-chemical strain ε tc is defined as where ϕ and φ refer to the coefficient of shrinkage expansion (CSE) and thermal expansion (CTE) in each orientation. Thermal expansion strain, ε th , and cure shrinkage strain, ε ch , are contained as free thermo-chemical strain. ε tot , τ and t refer to total strain, past time, and current time. C denotes the relaxation stiffness matrix as functions of DOC, temperature, and time. This work adopts single underlines and double underlines to indicate the vectors and matrices. This work can denote Equation (5), as below, as being in line with the features taken on by thermorheologically simple materials at a constant DOC: in which: where α T is the shift factor which allows for time-temperature superposition, α ref references DOC, and ξ(t) and ξ(τ) indicate the past reduced time and current reduced time, respectively. The relaxation modulus of a viscoelastic material is denoted via a series of Prony, as shown below, in light of the generalized Maxwell model with n Maxwell elements in parallel [26]: Here, W i is weigh factor, and E u and E ∞ indicate the adequately unrelaxed and relaxed modulus. The non-relaxation materials properties and proposed relaxation modulus, inclusive of CSE, CTE, Poisson's ratio, shear modulus, and modulus related to elasticity, can be adopted to acquire the performance relative to viscoelasticity taken on by the anisotropic composite laminates precisely in the light of the micro-mechanics theory [1]. Table 3 lists the unrelaxed mechanical properties (namely, mechanical properties relative to elasticity) of AS4/3501-6 prepreg constituents. Furthermore, this work is able to acquire the mechanical properties taken on byAS4/3501-6 prepreg through adopting micromechanical Equations (A.8)-(A. 16), as exhibited in Appendix A.

Modified CHILE Model
The CHILE model continues to be the essential, critical elasticity-related model. A modified CHILE model is able to be acquired to forecast residual stress and deformation of the viscoelastic materials as the work adopts an equivalent modulus corresponding to the relaxation modulus, related to viscoelasticity at certain frequencies or a specific time, though the CHILE model fails to shed light on the state of the polymer being in the viscoelastic regime in the course of the cure period. The acquired outcomes turn out to be nearly equal to those forecasted through adopting the viscoelasticity-related method. The work can denote the integral form of the CHILE constitutive function as: where E indicates the instant elasticity-related modulus that is estimated. And more notably, the temperature and DOC change along with time in the course of the cure cycle. Thus, the constitutive equation has the form of the following equation: Evidently, Equations (7) and (11) shall be identical ones as the function below is effectuated [30]: Here, we can see that Equation (12) seems odd, because the part to the right of the equal sign is a relaxation modulus related to both t and τ, whereas the left-hand side is a modulus related to elasticity, which is merely an equation of τ. This function shall be expounded in the paragraphs of this section below. To assess Equation (11), modulus related to elasticity E, at any integration time τ, is required to be assessed to meet Equation (12), i.e., the modulus related to elasticity, and E can be defined as the value of the relaxation modulus at time τ, C, at an appropriate undetermined time, t e , under the DOC and temperature at time τ. Accordingly, this work can denote Equation (12) as: where t e indicates an undetermined and appropriate time. It is noteworthy that the DOC and temperature are functions of τ for the relaxation modulus, C, and the modulus related to elasticity. The work can denote the relaxation modulus at time t e as below, in light of the principle of time-temperature superposition: Evidently, the function below is required for Equations (12)- (14): A necessary and sufficient condition for Equation (15) is given for Equation (16): Substituting Equation (8) into Equation (16), the function can be reorganized below: Meanwhile, the determination of t e is counted as a problem that is critically difficult. To reckon with this, the integral in the Equation (17) falls into two parts: the cool-down stage and the hold stage.
where t f denotes the time at the preliminary cool-down stage and the end of the hold stage. The second integral belongs to the cool-down regime, whereas the first one belongs to the hold regime. Accordingly, this work can assess the time t e as below through deducing the two integrals [30]: Eventually, the state of stress for the composites is able to be accurately acquired through adopting Equations (10), (13), and (19).

Numerical Formulation
Residual stress strongly effects the properties of composite parts manufactured with thermosetting polymers. The residual stress and deformation for the composites was simulated through adopting the viscoelastic approach and the modified CHILE method. This work draws a comparison between the two sets of outcomes, and also references White and Kim's research [34,35]. That it is conducive to anatomize the simulation of optimization for the evolution of residual stress and deformation.
Two models of simulation were encompassed by three coupled modules; mechanical modules, instant heat conduction, and cure kinetics. Given that the output data of some modules count as the input data of others, the order of loading these modules is of critical significance in the model. Figure 1 presents flow of the viscoelastic model for the simulation of residual stress and deformation, in light of the loading order of these modules. Instant heat conduction, in light of Equation (1), is loaded, and it is the first loading module that provides curing temperature for other modules and counts as the critical basis in curing process. The cure kinetics is loaded on the second module, adopting Equation (3), and offering DOC for the mechanical module. As DOC and temperature are correlated, the modified CHILE model is in light of Equation (10), and the viscoelastic model adopts Equation (5). At last, the mechanical module is loaded. It is noteworthy that, as several thermo-physical properties transmit to heat conduction module require DOC captured from the second module, the first and second modules call data with each other in each time increment.
In COMSOL Multiphysics (Version 4.3b), such a process should be a repeated iteration of each increment until the desired precision is reached (for this article we selected 0.001). In this respect, when the thermodynamic analysis of coupling is performed, that is, the local error, the error is within one time. The global error is the estimate of the sum of local errors, in excess of, or under, the sum of local errors. Therefore, after the proper accuracy test, the calculation precision can be improved effectively.
This work adopts some normal COMSOL Multiphysics modules to realize the two models of process. This work effectuates the instant heat conduction module through introducing the thermo-physical property parameters presented in Table 2, and through adopting the "Heat Transfer" Application Module. Additionally, in this Application Module, this work establishes the "Heat Source" term through factoring in chemical heat release, bound by the DOC (Equation (2)). In this work we implemented the module of Cure kinetics through adopting the "Coefficient format PDF" Application Module, taking on user-defined function. The general function of this module is defined below: where f is the source term; c, a, e a , d a , α, β and γ are all coefficients; u and t are dependent variable and instant time, respectively. The cure rate functions (Equations (3) and (4)) are established in Equation (20) in this paper through adopting the appropriate coefficients. The DOC can be acquired by following the cure kinetic parameters listed in Table 1. Eventually, this work adopts the "Viscoelastic Material" item to effectuate the mechanical module in viscoelastic model in line with Equation (5). Such item is selected from "Structural Mechanics" Application Module. In the "Long-Term Elastic Properties" of the "Domain Selection", "Young's modulus and shear modulus" item is specified as well as the defined relaxation modulus, which is set in the "Generalized Maxwell Model" tab, presented from Equation (9), and the material parameters listed in Table 3 are accessed to the related fields. In COMSOL Multiphysics (Version 4.3b), such a process should be a repeated iteration of each increment until the desired precision is reached (for this article we selected 0.001). In this respect, when the thermodynamic analysis of coupling is performed, that is, the local error, the error is within one time. The global error is the estimate of the sum of local errors, in excess of, or under, the sum of local errors. Therefore, after the proper accuracy test, the calculation precision can be improved effectively.
This work adopts some normal COMSOL Multiphysics modules to realize the two models of process. This work effectuates the instant heat conduction module through introducing the thermophysical property parameters presented in Table 2, and through adopting the "Heat Transfer" Application Module. Additionally, in this Application Module, this work establishes the "Heat Source" term through factoring in chemical heat release, bound by the DOC (Equation (2)). In this work we implemented the module of Cure kinetics through adopting the "Coefficient format PDF" Application Module, taking on user-defined function. The general function of this module is defined below:  (3) and (4)) are established in Equation (20) in this paper through adopting the appropriate coefficients. The DOC can be acquired by following the cure kinetic parameters listed in Table 1. Eventually, this work adopts the "Viscoelastic Material" item to effectuate the mechanical module in viscoelastic model in line with Equation (5). Such item is selected from "Structural Mechanics" Application Module. In the "Long-Term Elastic Properties" of the "Domain Selection", "Young's modulus and shear modulus" item is specified as well as the defined relaxation modulus，which is set in the "Generalized Maxwell The mechanical module in the modified CHILE model is implemented in the "Elastic Mechanics" item chosen from "Structural Mechanics" Application Module in the light of Equation (11). In line with Equations (13) and (19), time variable t e is set to participate in the definition of modulus related to elasticity for obtaining an equivalent modulus that corresponds to the viscoelastic relaxation modulus.

Material Model Verification
To judge the accuracy and efficiency of the two proposed models for simulation, this paper adopts a same-thickness-plies orthotropic composite laminate and a stacking sequence of [0 • /90 • /90 • /0 • ] (as presented in Figure 2) as an instance, with the size ascertained as 10.16 cm × 10.16 cm × 2.54 cm.
An ideal contact interface is set between layers, that is, the natural boundary (Newman boundary) condition. The work adopts the AS4/3501-6 prepreg encompassing AS4 fibers pre-impregnated with 3501-6 epoxy resin as the major material, where the fiber volume fraction taken on by composite materials reaches 50%. The cross-ply laminate of AS4/3501-6 is acquired through introducing the related parameters of the material to the proposed simulation models.
The period of curing recommended from the manufacturer is used in these simulations which contains five stages: first heat-up period with a ramp rate of 2.5 • C/min from 25 • C to 116 • C, first 1-h hold period at 116 • C, second heat-up period with the same ramp rate to 177 • C, second 2-h hold period at 177 • C, eventually, cool-down period with a ramp rate of −2.5 • C /min from 177 • C to 25 • C. Two boundary conditions, heat transfer and mechanical pressure, are used to simulate the tool-part interfaces, which can be set in the instant heat conduction module and the mechanical module, respectively. The mold heating temperature is defined on the outside surface of the composite laminates, while the pressure loaded on the three planes is kept during the whole cure cycle (as presented from Figure 3). relaxation modulus.

Material Model Verification
To judge the accuracy and efficiency of the two proposed models for simulation, this paper adopts a same-thickness-plies orthotropic composite laminate and a stacking sequence of (as presented in Figure 2) as an instance, with the size ascertained as 10.16 cm × 10.16 cm × 2.54 cm. An ideal contact interface is set between layers, that is, the natural boundary (Newman boundary) condition. The work adopts the AS4/3501-6 prepreg encompassing AS4 fibers pre-impregnated with 3501-6 epoxy resin as the major material, where the fiber volume fraction taken on by composite materials reaches 50%. The cross-ply laminate of AS4/3501-6 is acquired through introducing the related parameters of the material to the proposed simulation models.
The period of curing recommended from the manufacturer is used in these simulations which contains five stages: first heat-up period with a ramp rate of 2.5 °C/min from 25 °C to 116 °C, first 1h hold period at 116 °C, second heat-up period with the same ramp rate to 177 °C, second 2-h hold period at 177 °C, eventually, cool-down period with a ramp rate of −2.5 °C /min from 177 °C to 25 °C. Two boundary conditions, heat transfer and mechanical pressure, are used to simulate the tool-part interfaces, which can be set in the instant heat conduction module and the mechanical module, respectively. The mold heating temperature is defined on the outside surface of the composite laminates, while the pressure loaded on the three planes is kept during the whole cure cycle (as presented from Figure 3).   In addition, it is necessary to verify the mesh to obtain an appropriate mesh, the shape and refinement of which impact on the convergence and accuracy of simulation. In this model, a hexahedral mesh is employed as refined mesh (as presented in Figure 3). The number of elements per layer is 2304, and the total number is 9216.

Thermo-Chemical Analysis
The temperature and DOC in central point (5.08, 5.08, 1.27) for the laminate between the modified CHILE model and the viscoelastic model are compared in Figures 4 and 5 in the case of recommended period of curing for manufacture. The work also exhibits the other outcome provided In addition, it is necessary to verify the mesh to obtain an appropriate mesh, the shape and refinement of which impact on the convergence and accuracy of simulation. In this model, a hexahedral mesh is employed as refined mesh (as presented in Figure 3). The number of elements per layer is 2304, and the total number is 9216.

Thermo-Chemical Analysis
The temperature and DOC in central point (5.08, 5.08, 1.27) for the laminate between the modified CHILE model and the viscoelastic model are compared in Figures 4 and 5 in the case of recommended period of curing for manufacture. The work also exhibits the other outcome provided by Kim and White [34]. It is noteworthy that the laminate is centered (5.08, 5.08, 1.27) between layer 3 and 2. The center temperature goes under the room temperature in the first heat-up period, (period of curing) as the surface is heated first (as presented from Figure 6a), mean time, the cross-linking rate of resin is slow. The temperature at the center rises sharply with the first peak at 46 min in the vicinity of 121 • C (as presented in Figure 6b) in the following hold period, while the cross-linking rate accelerates as the outcome of heat conduction and reaction heat accumulation. In the second heat-up period, the center temperature remains lower than the room temperature, arising from the low coefficient of heat conduction and heat transfer lag (as illustrated in Figure 6c), but the cross linking is more intensive, and the reaction rate increases significantly. In the second holding stage, arising from the internal chemical heat release and the accumulation of heat conduction of the model, high temperature area transfers to internal and the center temperature reaches the maximum temperature peak of about 121 • C at 133 min (as presented in Figure 6d). Eventually, the center temperature of the model is gradually close to room temperature (as presented in Figure 6e,f), and the DOC nearly tends to 1. The simulation results of the two models are identical and coincide with those attained in Kim and White's study. It indicates that the heat conduction and the cure kinetics modules are quite qualified to simulate thermo-chemical response for the composite laminate during the cure cycle.

Residual Stress Analysis
The major factor leading to residual stress generation is the cure shrinkage, i.e., the chemical shrinkage and thermal strains arising from CTE (Equation (6)). Figure 7 shows the altered state of

Residual Stress Analysis
The major factor leading to residual stress generation is the cure shrinkage, i.e., the chemical shrinkage and thermal strains arising from CTE (Equation (6)). Figure 7 shows the altered state of thermal strain in different directions at the points 5.08, 0, and 1.27 for the composite laminate in the course of cure cycle. Thermal strain in the longitudinal direction is negligible, while those in through-thickness and transverse directions change noticeably, particularly during two heat-up periods. Additionally, the modified and viscoelastic CHILE models are acquired to well conform to each other. In Figure 8, the two simulation models perform well in estimating cure shrinkage from the simulated DOC in the course of cure cycle. Arising from the lateral pressure, the cure shrinkage in the through-thickness direction is evident, whereas the cure shrinkage in x and y directions can be neglected. As indicated, the cure shrinkage in the through-thickness direction computed by the modified CHILE model is in excess of that acquired by the viscoelastic model. Additionally, the modified and viscoelastic CHILE models are acquired to well conform to each other. In Figure 8, the two simulation models perform well in estimating cure shrinkage from the simulated DOC in the course of cure cycle. Arising from the lateral pressure, the cure shrinkage in the through-thickness direction is evident, whereas the cure shrinkage in x and y directions can be neglected. As indicated, the cure shrinkage in the through-thickness direction computed by the modified CHILE model is in excess of that acquired by the viscoelastic model.     Figure 9 shows the simulated outcomes of modulus related to elasticity development, acquired by the linear elastic model [36] and the modified CHILE model in through-thickness direction in the course of cure cycle. Obviously, the modulus related to elasticity, acquired by the modified CHILE approach, increases slowly during curing and continues growing during the cool-down period. This outcome is contrary to the one estimated by the linear elastic model.  Figure 9 shows the simulated outcomes of modulus related to elasticity development, acquired by the linear elastic model [36] and the modified CHILE model in through-thickness direction in the course of cure cycle. Obviously, the modulus related to elasticity, acquired by the modified CHILE approach, increases slowly during curing and continues growing during the cool-down period. This outcome is contrary to the one estimated by the linear elastic model. To prove whether the stress constitutive functions (Equations (5) and (10)) are feasible and accurate in the two proposed models, this work acquired transverse stress 2  in the  0 ply at x = 5.08 and y = 5.08 of the laminate and the development of interlaminar normal stress 3  at the points 5.08, 0, and 1.27, presented from Figures 10 and 11. As indicated from the foregoing figures, the outcomes of the viscoelastic model and the modified CHILE model conform to the outcomes presented by White and Kim [36]. Particularly, the excellent consistency between the viscoelastic model and the results presented by White and Kim was found. As illustrated in Figure 10, the final stress for the normal stress from the modified CHILE model is 24.5 MPa, which is 6.1% in excess of the viscoelastic solution, while the final stress for the transverse stress from the modified CHILE model in Figure 11 is 33 MPa, which is 9.1% in excess of the viscoelastic solution. The von Minses stress distribution of the composite laminate after cool-down, forecasted respectively by the modified and viscoelastic CHILE models, is illustrated in Figure 12. To prove whether the stress constitutive functions (Equations (5) and (10)) are feasible and accurate in the two proposed models, this work acquired transverse stress σ 2 in the 0 • ply at x = 5.08 and y = 5.08 of the laminate and the development of interlaminar normal stress σ 3 at the points 5.08, 0, and 1.27, presented from Figures 10 and 11. As indicated from the foregoing figures, the outcomes of the viscoelastic model and the modified CHILE model conform to the outcomes presented by White and Kim [36]. Particularly, the excellent consistency between the viscoelastic model and the results presented by White and Kim was found. As illustrated in Figure 10, the final stress for the normal stress from the modified CHILE model is 24.5 MPa, which is 6.1% in excess of the viscoelastic solution, while the final stress for the transverse stress from the modified CHILE model in Figure 11 is 33 MPa, which is 9.1% in excess of the viscoelastic solution. The von Minses stress distribution of the composite laminate after cool-down, forecasted respectively by the modified and viscoelastic CHILE models, is illustrated in Figure 12.
model and the results presented by White and Kim was found. As illustrated in Figure 10, the final stress for the normal stress from the modified CHILE model is 24.5 MPa, which is 6.1% in excess of the viscoelastic solution, while the final stress for the transverse stress from the modified CHILE model in Figure 11 is 33 MPa, which is 9.1% in excess of the viscoelastic solution. The von Minses stress distribution of the composite laminate after cool-down, forecasted respectively by the modified and viscoelastic CHILE models, is illustrated in Figure 12.

Curing Deformation Analysis
This work introduces the Maximum deformation, defined as the maximum absolute displacement in the through-thickness direction, to present the gap of the deformation between the

Curing Deformation Analysis
This work introduces the Maximum deformation, defined as the maximum absolute displacement in the through-thickness direction, to present the gap of the deformation between the modified CHILE model and the viscoelastic CHILE model. The curing deformation contour forecasted by the two models after cool-down is indicated in Figure 13. Evidently, the laminate curves identically from four corners arising from the symmetrical and uniform laying of the fibers, as Figure  2 shows. The degree of deformation forecasted by the modified CHILE model is slightly higher than

Curing Deformation Analysis
This work introduces the Maximum deformation, defined as the maximum absolute displacement in the through-thickness direction, to present the gap of the deformation between the modified CHILE model and the viscoelastic CHILE model. The curing deformation contour forecasted by the two models after cool-down is indicated in Figure 13. Evidently, the laminate curves identically from four corners arising from the symmetrical and uniform laying of the fibers, as Figure 2 shows. The degree of deformation forecasted by the modified CHILE model is slightly higher than that of the viscoelastic model. The maximum deformation estimated by the modified CHILE model is 0.79 cm, which is 0.06 cm more than the viscoelastic model (0.73 cm).

Discussion
The previous sections discuss the differences between the viscoelastic model and the modified CHILE model. Besides forecast accuracy, a significant difference of computation time is also observed between the two models. In heat-up and hold periods, i.e., during curing, the modified CHILE model runs approximately 10 times faster than the viscoelastic model, while the multiple approximately reaches 20 times in the cool-down period. Clearly, the more complex the composite structure is, the more prominent the time benefit of the modified CHILE approach is.

Conclusions
In this paper, both the gold standard viscoelastic model and the modified CHILE model have been presented to forecast process-induced residual stress and deformation for the thermosetting resin composite laminates. The outcomes from these two kinds of fully 3D coupled simulation models are compared with each other to show their trade-offs. To improve the prediction accuracy, the evolution of the thermo-physical properties, with the temperature and DOC, is also considered.
This work establishes a four-layer composite laminate to validate two simulation models. A variable time parameter is introduced in the modified CHILE model to factor in relaxation features taken on by the materials. The outcomes show that the forecasted curing temperature, DOC, and remaining stress by the two proposed models are in good agreement with the results obtained by Kim and White. However, there is a significant difference in modulus related to elasticity between the modified and elastic CHILE models. Furthermore, the maximum deformation attained from the modified CHILE model is evidently in excess of that acquired by the viscoelastic model within acceptable limits.
Additionally, although the modified CHILE model has a lower accuracy than the viscoelastic model, the viscoelastic model, within the margin of error, can be replaced by the modified model in the numerical modeling of process-induced residual stress and deformation for a shorter computation time, and a far more efficient numerical implementation. As previously mentioned, in

Discussion
The previous sections discuss the differences between the viscoelastic model and the modified CHILE model. Besides forecast accuracy, a significant difference of computation time is also observed between the two models. In heat-up and hold periods, i.e., during curing, the modified CHILE model runs approximately 10 times faster than the viscoelastic model, while the multiple approximately reaches 20 times in the cool-down period. Clearly, the more complex the composite structure is, the more prominent the time benefit of the modified CHILE approach is.

Conclusions
In this paper, both the gold standard viscoelastic model and the modified CHILE model have been presented to forecast process-induced residual stress and deformation for the thermosetting resin composite laminates. The outcomes from these two kinds of fully 3D coupled simulation models are compared with each other to show their trade-offs. To improve the prediction accuracy, the evolution of the thermo-physical properties, with the temperature and DOC, is also considered.
This work establishes a four-layer composite laminate to validate two simulation models. A variable time parameter is introduced in the modified CHILE model to factor in relaxation features taken on by the materials. The outcomes show that the forecasted curing temperature, DOC, and remaining stress by the two proposed models are in good agreement with the results obtained by Kim and White. However, there is a significant difference in modulus related to elasticity between the modified and elastic CHILE models. Furthermore, the maximum deformation attained from the modified CHILE model is evidently in excess of that acquired by the viscoelastic model within acceptable limits.
Additionally, although the modified CHILE model has a lower accuracy than the viscoelastic model, the viscoelastic model, within the margin of error, can be replaced by the modified model in the numerical modeling of process-induced residual stress and deformation for a shorter computation time, and a far more efficient numerical implementation. As previously mentioned, in other generalized processing cases, the modified CHILE method was able to be adopted to delve into the impacts exerted by the thickness of the composite laminates, the stacking sequence, and the fiber volume fraction on curing-induced residual stress and deformation development.