An efficient reduced simulation of residual stresses in composite forming processes

An efficient simulation of thermal and mechanical models involved in thermosetting composites forming needs to overcome some numerical difficulties related to: (i) the multi-scale behaviour; (ii) the complex geometries involved needing too many degrees of freedom; (iii) the large time intervals where the solution has to be computed; (iv) the non-linearity of the involved evolution equations; (v) the numerous couplings… In this work, an efficient strategy based on a separated representation is proposed. This method enables to avoid the use of an incremental strategy and can lead to impressive computing time savings especially when the model involves fine meshes and very small time steps. The local non-linear chemical kinetics and its coupling with the global heat balance equation are naturally introduced in the separated representation algorithm. Moreover, the dependence of the thermal conductivity and the specific heat on the temperature and on the reaction advancement degree are also taken into account. Knowing the history of the temperature field, the separated representation is again used to solve the thermo-mechanical problem in order to determine the residual stresses.


Introduction
There are many recent issues concerning the modelling and simulation of composites manufacturing processes. Indeed, composite materials are more and more used for their mechanical properties and for their light weight. An efficient simulation of thermosetting curing process must be able to accurately predict the temperature and the degree of reticulation through the entire process in order to ensure that the curing is homogeneous and then to avoid some defaults in the finished parts. During the process some residual stresses also appear, coming in particular from local variations of the thermal expansion and the chemical shrinkage. Obviously, the minimisation of these residual stresses is of capital interest for the properties of composites structures.
Some relevant works using finite element approach have investigated the curing process, in which the thermal, the chemistry and the mechanics are coupled. Rabearison et al. [1] took into account these couplings in their numerical simulation of 3D structural parts during the curing. A satisfactory comparison is obtained between their numerical results and the experimental data of the local temperature history. The co-workers [2,3] analysed the autoclave curing process for thermoset matrix composites by a finite element method. In the first paper [2], they used an optimization technique to recurrently enhance the curing process. Furthermore, a fully 3D coupled thermo-chemo-viscoelastic model has been used to simulate the warpage and the springforward of thin cross-ply laminates in curved parts [3]. Cheung et al. [4] presented a coupled thermo-kinetic simulation to obtain the thermal equilibrium and chemical kinetics during the curing of Resin Transfer Molding (RTM) process. They developed a numerical strategy which allows solutions of large mesh problems with the resources of a personal computer.
From a numerical point of view, there are several difficulties related to the simulation of composites processing. A first one is related to the growing need of the size and complexity of composite parts. Thus, too many nodes are usually required to perform accurate simulations. The number of time steps needed in the time discretization also increases drastically due to the differences of the characteristic times of heat diffusion and chemical kinetics. Many other difficulties should be addressed, being the most important the multi-physics coupling and the strong non linearity of the resulting model.
In order to reduce the computation time related to numerical simulations, a new family of solvers, called the Proper Generalized Decomposition method, or in short the PGD method, has been developed [5][6][7][8]. It is an "a priori" resolution technique, that derives from a POD decomposition (Proper Orthogonal Decomposition) without requiring any basis or known solutions. The solution is progressed by resolving some spatial problems (time independent) and temporal problems (ordinary differential equation). A lot of recent works uses the PGD to solve some multidimensional problems [7][8][9][10].
The aim of this work is to propose a reduced modelling strategy to perform accurate and efficient simulations of composite forming processes. The considered model involves a weak coupling of the thermal and the mechanical models because the mechanic is strongly affected by the thermal effects. As a first approximation, the effects of mechanic on the thermal balance are ignored.

Themo-mechanical modelling
The model used in our study is the one proposed by Abou Msallem et al. [11] and Abou Msallem [12].

Thermal model
The heat balance equation in presence of curing kinetics writes: where T is the temperature, k is the homogenized thermal conductivity tensor, C p is the specific heat, ρ is the density of the material, H is the total heat generated during the reaction of the resin (J/kg), w is the mass fraction of resin in the composite and α is the degree of advancement of the chemical curing reaction. This equation is non linear because both the thermal conductivity and the specific heat depend on T and α.
In this paper, we are focusing on plate or shell geometries that only involve the heat transfer in the thickness direction. In this case, the thermal model reduces to the one dimensional heat equation, which involves that the conductivity is a scalar.
There are several models describing the evolution of the degree of curing. One of the most used models is the one proposed by Kamal and Sourour [13]. However, the present work deals with the Bailleul's model [14] that can be written in a simple polynomial form: where coefficients α i are material parameters and the prefactor reads: with K ref , A and T ref additional material parameters.

Mechanical model
Once the thermo-kinetic problem is determined, the deformations induced by the thermal expansion a T À T ref À Á Â Ã and the chemical shrinkage b a À a ref À Á Â Ã can be calculated. Therefore, the thermo-chemo-elastic constitutive equation involving the just referred ingredients writes: where H is the fourth order stiffness tensor and, α and β are the thermal expansion and the chemical shrinkage tensors, respectively.

The proper generalized decomposition revisited
For alleviating the computational cost related to numerical simulations, it has been proposed recently the utilization of reduced modelling based on the use of proper generalized decompositions (PGD) making use of spacetime separated representations, which leads to nonincremental strategies. This technique was successfully applied for treating some multi-dimensional models encountered in the kinetic theory description of complex fluids [7,8]. In what follows, and for the sake of simplicity, this technique is described in the context of the solution of the 1D heat equation (i.e. the part related to the specific heat has been voluntary omitted): Now, a separated representation of the temperature field is assumed involving N functional couples F i ; G i ð Þ f g i¼1;...;N such that: where F i f g i¼1;...;N and G i f g i¼1;...;N are defined in the time interval I ¼ 0; t max ½ and in the space domain Ω, respectively. The weak formulation of the heat equation can be expressed as follows: Find T verifying the essential boundary conditions, such that: for any test function T * defined in an appropriate functional space.
First we assume that, at iteration n, some of the functional couples are already computed: ...;n with 1 n < N . As these n functional couples are not enough to describe the solution, new functional products should be added. Thus, we must compute the couple RðtÞ; SðxÞ ð Þ . This calculation is performed by using an alternating directions fixed point algorithm [15], which is described in detail below.
First, it is expected that R(t) is known, therefore we can compute S(x), then R(t) again, and so on until reaching convergence. After reaching the convergence, we can identify the couple (R,S) to F nþ1 ; G nþ1 ð Þ . In what follows, we are illustrating the numerical procedure. We start by postulating that: which implies for the test function that: Introducing all these approximations into the weak formulation results in: The alternating directions fixed point schema proceeds as follows.

Computing the function S
By assuming that R is known, R * vanishes and Eq. 10 reduces to: As all the functions depending on time are assumed to be known, the time integration can be performed numerically. The remaining problem, defined in the space domain, can be solved by using an appropriate discretization technique.

Computing the function R
Once the function S is known, the next step consists in computing the function R. At this step, S * vanishes and then the weak formulation writes: At once the space integration can be performed leading to an ODE (Ordinary Differential Equation) related to the evolution of the function depending on time.

Checking the convergence
We need to define two convergence criteria to complete the algorithm. The first one is related to the alternating direction fixed point scheme used for computing the basis enrichment, that is, for each functional couple RðtÞ; SðxÞ ð Þ . For this purpose, we use: where ε 1 is a small enough parameter and ||·|| defines a norm on 0; t max ½ ÂΩ. The second convergence check concerns the enrichment stopping criterions. In our simulations, we employ a residual based criterion expressed as: where ε 2 is a small parameter.

Global simulation strategy
The text below summarizes the global reduced strategy that has been used to solve the thermal model previously described (Eq. 1). This strategy consists of 3 steps repeated until reaching the global convergence. The general algorithm is depicted in Fig. 1.

First step
We assume that the degree of curing and the temperature is known at each node position and at any time. Thus, k and C p can be determined everywhere and through the history of the curing process from: where the functions α k , b k , a C p et b C p depend linearly on α(x,t). In the 1D case (chosen for the sake of simplicity), the heat equation writes: This equation can be solved using the following separated approximation: In order to apply the strategy previously described, each function involved in the model must be written on a separated form in order to ensure the separated treatment of each space. For this purpose, we should describe the degree of curing in a separated form. The second difficulty concerns the non linearity that is addressed using an appropriate linearization. There are many possibilities. The simplest one consists in calculating all the linear terms at the previous enrichment iteration.

Second step
We assume now that T is known and we want to determine the prefactor K that appears in the reaction kinetics Eq. 3 in a space-time separated form. For this purpose, first, we consider the exponential term: which can be rewritten as: As T is known on a separated form, a proper generalized decomposition of e K can be performed using the same algorithm that is used to solve the PDE (Partial Differential Equation). Now, we come back to: which can be approximated by using power series. This is done by computing the powers of e K (until e K p =p! becomes small enough). Knowing a separated form of e K, we can easily find a separated form of e K p .

Third step
This step concerns the solution of Eq. 2 related to the curing reaction advancement itself. It is a local equation which does not involve any a priori particular difficulty. Nevertheless, when the number of nodes becomes too large, the integration of the associated ODE at each nodal position becomes expensive. Moreover, after computing all these solutions, one should perform a SVD (Singular Value Decomposition) of the degree of advancement of the curing reaction α in order to introduce this separated approximation in the heat equation. With the separated representation of the prefactor Initialisation of T and α n is set to 1 Step 1: Resolution of the heat equation with: Step 2: Computation of K(T) in a separated form

Results
We focus on a 1D thermal problem defined in the thickness of an epoxy resin plate. The plate thickness is 5 cm. Figure 2 depicts the curing cycle applied to the plate. It is assumed that the temperature on boundaries is perfectly controlled during the process and follows exactly the curing cycle. For the chemical reaction, α is set to 0 everywhere at the beginning of the curing process. The proposed algorithm (Fig. 1) is applied to solve this problem. Figure 3 depicts the temperature and the curing degree versus time and position during the process. Figure 4 focuses on the evolution at the center and at the edge of the plate. We can see that there is a peak in temperature originated by the exothermic chemical reaction. This is a critical point because if the peak is too high, the material can be damaged. We can also encounter from the simulation prediction that the plate will be well cured everywhere. Another interesting point is the cooling because it is in this phase that some residual stresses can appear. The curing cycle has to be optimized according to all these constraints.
From the numerical point of view, the convergence is reached quickly and only ten iterations are required to perform the proper generalized decomposition of the thermal field. The computation needs a couple of minutes on a standard personal laptop for a model involving 500 nodes in the thickness and 1,000 time steps. The degree of the interpolation functions used is one.
If we increase the number of nodes by a certain factor, the computation time of the proposed strategy will grow linearly, instead of the exponential growing of standard incremental strategies. The finer are the considered space and time discretizations, the more are the computing time savings. Thus, this strategy is expected to be especially suitable for simulating forming process involving very complex parts, when more experienced strategies (incremental ones) are almost inefficient.

On the proper generalized decomposition of the mechanical problem
This part deals with the mechanical problem that concerns the evaluation of residual stresses at the end of the process.

Model formulation
In order to solve the mechanical problem, the linear elasticity and the hypothesis of infinitesimal deformations are assumed. Firstly, let us introduce the equilibrium equation of motion as follows: where σ and f represents the Cauchy stress tensor and the body forces, respectively. Both quantities are assumed to be  . The deformation tensor ε and the displacement field u are related from the usual relationship: To completely formulate the mechanical problem, it remains to express a governing equation. For the sake of simplicity, the chemical shrinkage is neglected (by taking into account the chemical shrinkage, no further complexity is added to the problem). Therefore, the thermo-elasticity constitutive law may be written as: where H is a fourth-order elasticity tensor and α is a second-order tensor describing the thermal expansion. T ref represents the reference temperature, whereas T is the temperature at time t, which is computed according to the above simplified linear heat diffusion equation defined in a homogeneous isotropic medium, and where the eventual existence of chemical reactions is introduced by the presence of a simple space dependant source term q(x) (in this section we are only interested in analyzing the thermoelastic coupling): The initial and boundary conditions for both heat and mechanical problems write: For solving this thermo-mechanical model, the problem is decoupled: the heat equation is assumed solved (its treatment was addressed in the previous section) and the proper generalized decomposition strategy is used again to determine the mechanical properties.

Proper generalized decomposition
The main aim of this section is to present a discretization based on the separated representations in order to find the displacement field u and therefore the Cauchy stress tensor σ associated to the model formulation. To apply the separated representation method, the displacement solution u may be approximated using N couples of functions T m i ; X m i À Á È É i¼1;...;N (the index m refers to the mechanical model) such that:  The temperature is assumed given in the separated form: where x j ; ϕ j n o j¼1;...;p are the p couples of functions needed for approximating the temperature field solution of the thermal model. The weak formulation of the mechanical problem can be expressed as follows: Find u verifying the essential boundary conditions such that: for all the test functions u * in an appropriate functional space.
Then, the functions involved in Eq. 32 are computed by assuming that the set of functional couples T m i ; X m i À Á È É i¼1;...;n with 0 n < N are already known. To enrich the approximation basis with the couple (R,S), an alternating directions fixed point algorithm is used. When the convergence is reached, the resulting couple (R,S) is used to update the next functional couple T m nþ1 ; X m nþ1 À Á . Therefore the separated representation at the n iteration yields: which suggests that the test function can be written as: Including all these assumptions into the weak formulation gives after rearrangement: Then, the alternating directions fixed point algorithm is used to determine the unknown functions S and R.
Computing the spatial function By assuming that R is known, R * becomes null and Eq. 37 (supposing that the dilation tensor time independent) reduces to: where: Finally, as the weak formulation (37) is satisfied for all S * , the associated strong formulation becomes: that can be solved by using an appropriate discretization technique for computing the space function S. In what follows, a Galerkin finite element method is used for this purpose.

Computing the temporal function
Once the function S is known, the next step consists in computing function R. At this stage, S * vanishes in Eq. 36 that reduces to: Equation 45 is valid for all R * , therefore it is possible to come back to the following strong algebraic formulation: which gives directly the function R. This equation does not involve time derivatives because the inertia terms were neglected in the mechanical model. The computation of functions R and S continues until reaching convergence. We assume that Q n+1 iterations were needed. When the convergence is reached we assign: T m nþ1 ¼ R and X m nþ1 ¼ S. Then, the enrichment procedure continues until verifying the approximation convergence criterion u N À1 ð Þ À u ðN Þ ", ε being a small enough parameter (u (N) refers to the separated approximation of the displacement field using n functional couples).

Numerical results
For illustrating the approach described beforehand, a basic 2D square domain Ω ¼ 0; 1 ½Â 0; 1 ½ is chosen with a time interval 0; t max ¼ 0:15 ½ (in this part, the units are voluntary omitted since all of them are expressed in the metric system). Moreover, we considered that x=(x, y). The problem is investigated under the plane stress hypothesis. The material is assumed to be homogeneous and isotropic, which reduces the tensor H to two independent elastic constants, the Young modulus E=10 3 and the Poisson ratio υ=0.3. The thermal expansion becomes α=αI, where the coefficient α is set to 1 and I represents the identity tensor. A diffusion coefficient k of 1 is applied. Furthermore, for the sake of simplicity, the body forces f are neglected, and the source term has the following form: q x; y ð Þ ¼ sin 2px ð Þsin 2py ð Þ. The initial and boundary conditions have already been mentioned in Eqs. 26 to 31. The mesh of the domain consists of 40×40 nodes defining a regular mesh of linear triangular finite elements. Linear interpolation functions have been used. Time integration is performed with a fully implicit scheme with a time step of 0.001.  Figure 5 depicts the temperature field obtained at time t=t max (i.e. the quasi steady state solution). The results are obtained from the temperature fields associated to the thermal model given by Eq. 25. The temperature profile is intentionally chosen to be time and space dependant and it appears as an input field into the thermo-elasticity constitutive law Eq. 24.
For this test, the proper generalized decomposition algorithm needs three couples of functions T m i ; X m i À Á i¼1;...;3 . Figure 6 shows the three temporal modes T m 1 , T m 2 and T m 3 . As the spatial modes X m i À Á i¼1;...;3 depend on the x-and the y-coordinates, Figures 7 and 8  ...;3 , plus the ones associated to the initial conditions, the displacement field is directly computed from Eq. 32. Figure 9 presents the quasi steady state solutions of the displacement in the x-and y-directions, respectively.  The accuracy of these solutions has been compared to the ones (superscript ref) obtained with the COMSOL Multiphysics® (formerly called FEMLAB) software. For instance, the error of the displacement in the x-direction is defined as u ref x À u x 2 , where u ref x and u x are the reference and computed solutions, respectively (the error in the ydirection is similarly obtained by substituting x for y). The thermo-mechanical problem has been solved with the FEMLAB software by using 2808 degrees of freedom. The maximum error in the x-direction is found to be 8.12E-03, whereas it is 5.55E-03 in the y-direction. These errors seem satisfactory as only four couples of function are used to determine the solution.
Finally Fig. 10 shows the resultant stresses obtained from the computed displacements. Fig. 10 (top) and (middle) depict the x-and y-stress in the 2D square domain, respectively. The shear stress is plotted in the Fig. 10 (bottom). Again theses solutions have been compared to the COMSOL Multiphysics® software. In general, the results are in a good agreement but some slight discrepancies are observed. The stresses computed by using both methods (i.e. separated representations and FEMLAB) exhibit slight mesh dependence. Simulations using finer meshes are then needed.

Discussion
Proper generalized decompositions applied to the solution of thermo-elastic models seems to be an appealing strategy. PGD solves N×Q space and time problems, where N is the number of functional couples (3 in the case just addressed) and Q is the maximum of the iterations needed for computing the different functional couples Q 1 ; . . . ; Q N f g (3 in the example previously addressed). Thus, one should solve N×Q 2D problems (time independent) and N×Q 1D problems (time discretization) instead the 2D discretization at each time step needed when using standard implicit incremental strategies. We would like to mention that the solution of the transient 1D problems does not introduce any difficulty (even for extremely short time steps). Therefore, for solving fine 3D models involving too short time steps (as the ones encountered when the heat equation involves reaction kinetics) the PGD strategy described in this paper seems to be the best way as compared to the standard incremental solvers.

Conclusions
In this paper, an efficient strategy based on the use of separated representations within the context of a proper generalized decomposition is proposed to solve models involving thermo-kinetic and thermo-mechanical couplings, as the ones usually encountered in the advanced modeling of composite manufacturing processes. PGD based schemes could reduce computing times of some order of magnitude, making possible the efficient and accurate solutions of models never satisfactory solved until now.
In the thermal model, the local non-linearity of the chemical kinetics and its coupling with the global heat equation induced some issues successfully by-passed by using separated representations.
Again, PGD has been applied successfully to solve the thermo-mechanical model from which the residual stresses are easily determined.
Up to now, the thermo-kinetic and thermo-mechanical problems were addressed separately. The ongoing work will focus on the coupling between these two efficient PGD based solvers deeply described in the present work.