A Multi-Scale Method for Predict Residual Stress Appearance in the Process of On-Line Consolidation of Thermoplastic Composites

. On-line consolidation process is a promising non-autoclave technique to design large size parts in aeronautical industry. Residual stress growth during the process is locking signif-icantly the process validation. In this paper, we present a multi-scale modeling based on the constrained natural element method (C-NEM) and on the proper orthogonal decomposition method which allows to evaluate residual stress growth. The method allows to take into account the thermomechanical properties evolution of the material during the process, which is a dominating factor in residual stress development.


Introduction
Processing of thermoplastic composites involves residual stress development. Residual stresses can lead to distortion of finished components, matrix cracking, and inter-ply delamination. Tensile stresses may even reach to a significant fraction of the tensile strength. Therefore, to ensure reliable and satisfactory performance, residual stresses within a composite part should be known before the part is put into use. The tape placement process is one of the few techniques that have the potential to continuously process thermoplastic composites in large scale industrial production. In the process, an incoming composite tape is bonded to the previously laid and consolidated laminate under heat and pressure locally applied to the interface (figure 1). By laying additional layers in different directions, a part with desired thickness and properties can be fabricated. Since heat is mainly supplied to the region where the tape and the laminate meet, the process is highly non-isothermal and, residual stresses may become critical. The process poses a very challenging problem of determining residual stresses for the following reason the tape placement is a continuous process. As such, when a new layer is being placed, the previously laid and consolidated layers are again subjected to heating. If the temperature of these layers exceeds the glass transition temperature, annealing will ensue and the residual stresses developed during previous tape placement could relax. Most of the previous studies on the prediction of residual stresses developed during thermoplastic composite processing concentrated on press molding (Chapman et al., 1990;Jeronimidis et al., 1988;Lawrence et al., 1990;Wang et al., 1992;Unger et al., 1993;Lee et al., 1994;Li et al., 1997;Wang et al., 1997;Domb et al., 1998;Ersoy et al., 2000;Ersoy, 1998). On the other hand, Nejhad et al. (1992) proposed a model to predict process induced stresses in filament winding. They modeled filament winding as a transient process. Since this led to excessive computational time, residual stresses could not be obtained within a reasonable time. Moreover, Residual stress analysis using a thermoviscoelastic finite element model was carried out by Sonmez et al. (2002). Sonmez et al., 2002]

Description of the occurring phenomena
During the welding, the changes of state of the matrix lead to the coexistence of three components (fiber, amorphous matrix and crystalline matrix) with different thermomechanical properties in a zone with high thermal gradients (see figure 2).

Figure 2. Schematic diagram of the on-line consolidation process and a modeling domain diagram representative of the phenomena occurring during the process
Unacceptable high levels of stress may arise because of two factors. First of all, thermal expansion behavior of fiber-reinforced thermoplastic composites is highly anisotropic. This is due to the large difference in the thermal expansion coefficients of the matrix and fiber materials. In a typical carbon-reinforced unidirectional thermoplastic composite, a temperature increase induces considerable expansion of the laminate in the transverse direction but almost none along the fiber direction. Besides, fibers, being themselves anisotropic, contribute to this effect. The other factor is the temperature gradients induced during the cooling process. This is because of the fact that, when different regions of the composite experience different temperature histories and stiffness properties of the material vary with the temperature, deformation behavior will be incompatible. Although elastic properties of fibers are almost unaffected by temperature, influence of temperature on the properties of the polymer matrix is drastic. Moreover, high processing temperatures typical for high performance thermoplastic composites exacerbate both the effect of temperature gradients and the effect of anisotropic thermal expansion behavior. Crystallization and melting phenomena occurring in the composite during process due to welding cause important evolution of the material properties. Indeed material thermomechanical properties depend on temperature and rate of crystallinity. These evolutions must be taken into account in the calculation, by using a multi-scale modeling.

Introduction
We consider that the viscous relaxation times are greater than the characteristic process time. Therefore, a thermoelastic behavior can be considered for the simulations. Moreover, a plane strain 2D analysis will be performed. The process simulation at the macroscopic scale required homogenized thermal parameters at each point and time. On the other hand, a microscopic modeling will be required to define the homogenized thermomechanical parameters (to be employed in the macroscopic analysis) as well as to computed the microscopic thermomechanical fields evolution allowing to check the criteria of default appearance. Thus, both scale analysis must be performed simultaneously, the macroscopic one in the whole macroscopic domain, and the microscopic one in some representative domains, from which the solutions could be interpolated everywhere in the whole domain. Thus, the proposed strategy is not a purely homogenization technique, but a multi-scale analysis where the microscopic analysis also serves to computed the homogenized (and then macroscopic) thermomechanical parameters. This modeling results coarser that a fully microscopic analysis as the one performed in (Ladeveze et al., 2001;Ladeveze et al., 2003;Farhat et al., 1991) by using high performance computing (domain decomposition running in parallel computing platforms). On the other hand, it differs also of standard multi-scale homogenization techniques in which the microscopic domain only serves to compute homogenized parameters but not for defining the whole time evolution of microscopic fields. Moreover it is well known that these last techniques have serious difficulties to represent real boundary conditions. In the technique here proposed the microscopic domains are strongly connected with the macroscopic mesh, in fact they result by zooming some finite elements of the macroscopic mesh. Thus, the scales connection is very simple. Of course in the case of an eventual incompatibility between the micro and macro meshes, the use of the Arlequin approach (Ben Dhia, 1998) or any superposition method (S-FEM) (Takano et al., 2004) could be envisaged. It can be noticed that the particular choice of the microscopic domains allows for (i) the imposition of real boundary conditions on their boundaries ; (ii) compute the whole time evolution of the microscopic thermomechanical fields ; and (iii) compute the homogenized thermomechanical parameters related to the elements allocating the microscopic domains, from which these parameters can be interpolated everywhere in the macroscopic domain. The microscopic zones must be representative, and then they must contain at least the tape which will be laid and the previous one (figure 3). Firstly, we define the locations of the microscopic domains which will allow us to describe the microscopic state of the composite during the process. A mesh is assigned to each microscopic domain. The large nodal density differences between the microscopic and macroscopic are accurately accounted thanks to the meshless character of the natural element approximation considered for the interpolation of the different fields (Sukumar et al., 1998;Yvonnet et al., 2005;Yvonnet et al., 2004). Microscopic domains are coupled to the macroscopic one through their boundaries. As described later, the natural element interpolation is based on the use of Voronoi cells instead of Delaunay triangles (used in the finite element framework) for defining the fields approximation. It has been proved in some of our former works (Martinez et al., 2001) that the quality of the Delaunay triangulation (dual of the Voronoi diagram) does not affect the quality of the computed solution in the NEM framework, even when one proceeds in very distorted meshes.

General algorithm
The global solution procedure is presented on figure 4. After the welding of two composite layers we assume known the initial temperature (which results from the numerical modelling of the welding operation) as well as the boundary conditions related to the particular industrial process considered. We suppose that we want to describe the state of the composite in the time interval {t 0 , t f }. The different steps of the calculation are presented below. We realize a partition of the time interval {t 0 , t f } in m time intervals by means of a time step Δt Macro

Macroscopic thermal calculation
We consider a regular mesh composed of bilinear quadrilateral elements. We perform the thermal calculation in the macroscopic time intervals {t k , t k + Δt Macro } , k ∈ {0, . . . , m − 1} by taking into account the crystallization kinetics as well as the homogenized thermomechanical parameters. This thermal calculation allows us to compute the macroscopic thermal field that is then transfered to the microscopic calculation.

Microscopic calculation
We consider a time step Δt micro , in general much lower than Δt Macro , i.e. Δt micro << Δt Macro . We perform p microscopic calculations (time steps) in a time   figure 5. On the boundary of each microscopic cell considered for the microscopic calculation we prescribe the temperature from the temperatures known at the four nodes defining that cell at t n = nΔt Macro and t n+1 = (n + 1) Δt Macro , assuming a linear evolution in both space -on the cell boundary -and time. Thus, at time t = t n + iΔt micro , for a point x located on the segment defined by nodes k and l we can write where s is the distance between the considered point and the node k and L kl is the distance between both nodes. We assume that a macroscopic time step Δt Macro corresponds with p microscopic time steps, i.e. Δt Macro = pΔt micro . In this manner the solution of the thermal problem offers the field T (X, , where X denotes the macroscopic position related to the microscopic domain (center of gravity of the microscopic cell) an x the point in the microscopic domain. Moreover we can evaluate the microscopic stress fields in any point of the microscopic domains and for any time, in order to evaluate the failure criteria.

Homogenization
In the considered microscopic time intervals, thermomechanical properties practically don't evolve. Therefore we can realize the homogenization of the thermomechanical properties in the time interval {t k , t k + Δt Macro } , k ∈ {0, . . . , m − 1} by means of the defined representative volume elements. The homogenized thermomechanical properties are then used in the macroscopic scale at next time increment in order to evaluate the macroscopic fields. The principle of the homogenization method and the transfer of the microscopic information towards the macroscopic scale is presented later.

Macroscopic stress calculation
The macroscopic stresses are calculated by means of the thermomechanical homogenized properties. This calculation will be described later.

Thermal modeling
In this section, we present the thermal modeling used for the simulation.

Energy balance
The equation to be solved, the transient 2D heat transfer equation including the heat generated by crystallisation, may be written as where ρ is the material density (kg.m −3 ), C p is the material specific heat (KJ.g −1 .K −1 ), k is the material thermal conductivity (W.m −1 .K −1 ), χ the relative degree of crystallinity, m m the matrix mass fraction, and H c the crystallisation enthalpy (KJ.g −1 ).

Crystallisation kinetics
An Avrami crystallisation kinetics law has been used to describe time crystallisation evolution where E a is the activation energy (KJ.mol −1 ), n the Avrami coefficient and t c the time elapsed from the beginning of the polymer crystallisation.

Matrix formulation
The thermal equation [3] can be written as N is the vector containing shape function, B contains the shape function derivatives, h is the thermal transfer coefficient (W.m −2 .K), and T amb is the ambient temperature (K).

Model reduction
Computing times can be significantly alleviated by using model reduction, like the proper orthogonal decomposition. The aim of the model reduction is to define a transformation matrix to get few basis functions to describe the spatially distributed state, with few degree of freedom. We present in this section a method which has been used to reduce computing times based on the proper orthogonal decomposition.

The proper orthogonal decomposition
We assume that the evolution of a certain field T (x, t) is known. In practical applications, this field is defined at the nodes of a spatial mesh x i and for some times . T m defines the vector containing the nodal degrees of freedom (temperatures) at time t m . The main idea of the proper orthogonal decomposition is to obtain the most typical or characteristic structure φ (x) among these T m (x) , ∀m. The maximization of which can be rewritten in the form Defining the vector φ such that its i-component is φ (x i ), [12] results in the eigenvalue problem [13], whose eigenvectors related to the highest eigenvalues define the characteristic structure of T m (x), where the two points correlation matrix is given by which is symmetric and positive definite. If we define the matrix Q containing the discrete field history then it is easy to verify that the matrix c in [13] results

A posteriori reduced modelling
We solve the eigenvalue problem defined by [13] selecting the eigenfunctions φ k associated with the eigenvalues belonging to the interval defined by the highest eigenvalue and that value divided by a large enough value. In practice n is much lower than N . Now, we could try to use these n eigenfunctions φ k for approximating the solution Now, if we consider the linear system of equations resulting from an implicit or semiimplicit time discretization of the thermal model defined by [6] and by multiplying both terms by B T we obtain where the size of B T G m B is n × n, instead N × N . When n << N as is the case in numerous physical models, the solution of [21] is preferred because its reduced size.

Enriching the approximation basis
Obviously, accurate simulations require an error evaluation as well as the possibility of adapting the approximation basis by introducing new functions able to describe the solution features. Ryckelynck proposed in (Ryckelynck, 2005) an adaptive procedure, able to construct or enrich the reduced approximation basis. For this purpose, he proposed to enrich the reduced approximation basis by adding some Krylov subspaces generated by the equation residual. Despite the fact that this enrichment tends to increases the number of approximation function, when it is combined with a POD decomposition that continuously reduces this number, the size of the problems is quickly stabilized. This strategy, which was successfully used in  as well as in , is summarized in the present section. We assume the solution accurately described in the interval ] 0, t s ] by using the reduced basis B. Now, the solution evolution is performed in the time interval ] t s , t s+1 ], with t s+1 = t s + SΔt solving [21] making use of the reduced approximation basis defined by matrix B When time t s+1 is reached, a control step is performed in order to evaluate the accuracy of the solution computed using the reduced basis. The control step is performed at t s+1 by computing the residual If the norm of the residual is small enough R < (being a small enough parameter) the computed solution can be assumed as good, and the time integration goes on in ] t s+1 , t s+2 ], with t s+2 = t s+1 + SΔt, using [22] without changing the reduced approximation basis.
On the contrary, if R > Bζ M , the approximation must be improved. For this purpose, we propose to enrich the reduced approximation basis by introducing the just computed residual (or some Krylov's subspaces generated by the residual, as proposed in (Martinez et al., 2001))

B ← [BR]
[24] and now, the evolution is recomputed in the ] t s , t s+1 ] using [22] with the just updated reduced basis B. When the convergence criterion is satisfied (i.e. R < ) we apply a POD on the entire past history ] 0, t s+1 ] for extracting the reduced basis able to represent all the past evolution of the unknown field, which results in an updated B that can be considered as the optimal reduced basis to describe the evolution of the unknown field in ] 0, t s+1 ]. Using the just updated reduced basis B the time evolution given by [22] is performed in ] t s+1 , t s+2 ] and the solution accuracy is checked at t = t s+2 .
The enrichment tends to increase the size of the reduced basis, but the POD reduces its size. The combination of both procedures allows to stabilize the size of the model as we noticed in some of some former works .

Microscopic description : dealing with highly heterogenous nodal distributions
In the present method, we define in the macroscopic mesh some microscopic domains which will allow us to determine the microscopic thermomechanical fields evolution as well as the homogenized thermomechanical properties during the simulation. In order to use the information coming from the macroscopic scale we must ensure the transition between the macroscopic and the microscopic scales. A microscopic domain is localized in a macroscopic element. The size of the representative volume element must be smaller then the macroscopic element. So the microscopic mesh nodal density evolves suddenly when ones moves from a microscopic domain to the transition one (which corresponds to the macroscopic scale). The nodal density inside the microscopic domain must be enough to represent accurately the microstructure, whereas in the transition domain it must be enough to represent the structure response. Thus, usually, the nodal distribution becomes highly heterogenous leading to very distorted meshes inadmissible if one proceeds in the finite element framework. For this reason, the use of a meshless discretization technique seems appropriate, because it is much less sensible to the aforementioned mesh distortions (see figure 6). In the next paragraphs, we briefly touch upon the basis of the natural element method that will be used for discretizing the variational formulation given by [3].
The natural element method is an appealing choice among the different meshless methods, because essential boundary conditions can be imposed directly, without detriment to the other properties (linear consistency, smoothness,...). Thus the finite element and natural element interpolations coupling becomes direct. It was originally proposed by Traversoni (Traversoni, 1994), Sambridge et al. (Sambridge et al., 1995) and was widely investigated in elastostatic problems by Sukumar (Sukumar et al., 1998). Moreover, it was successfully applied in fluid dynamic simulations involving polymer injection and metal forming simulations using updated Lagrangian formulations (Martinez et al., 2001;Martinez et al., 2004;Alfaro et al., 2006), in problems involving cracks (Yvonnet et al., 2004) and in phase change problems involving moving discontinuities (Yvonnet et al., 2005).

The constrained natural element interpolation revisited
We briefly touch upon the foundation of Sibson's natural neighbor coordinates (shape functions) that are used in the natural element method. For a more in-depth discussion on the Sibson interpolant and its application for solving second-order partial differential equations, the interested reader can refer to Braun and Sambridge (Sambridge et al., 1995) and Sukumar et al. (Sukumar et al., 1998). The NEM in-

Figure 6. Mesh of a microscopic cell and of the intermediate mesh realized with the method of the constrained natural elements (C-NEM)
terpolant is constructed on the underlying Voronoi diagram, the Delaunay tessellation being the topological dual of the Voronoi diagram. For the sake of simplicity we only consider the 2D case, the 3D case being a direct extension. Let S = {n 1 , n 2 , . . . , n N } be a set of nodes in 2 . The Voronoi diagram is the subdivision of 2 into regions T i (Voronoi cells) defined by where x i are the coordinates of node n i and d (., .) denotes the Euclidean distance. The Sibson coordinates of x with respect to a natural neighbor n i (see figure 7) is defined as the ratio of the overlap area (volume in 3D) of their Voronoi cells to the total area (volume in 3D) of the Voronoi cell associated with the point x If the point x coincides with the node n i , i.e. x = x i then ψ i (x i ) = 1, and all the other shape functions vanish, i.e. ψ i (x j ) = δ ij (δ ij being the Kronecker's delta). The properties of positivity, interpolation and partition of unity are then verified (Li et al., 1997

Figure 7. Construction of the Sibson shape functions
where n is the number of neighbor nodes of point x.
The natural neighbor interpolation satisfies the local coordinate property (Sibson, 1980), namely which combined with [27], implies that the natural neighbor interpolant spans the space of linear polynomials (linear completeness). Natural neighbor shape functions are C ∞ at any point except at the nodes, where they are only C 0 , and on the boundaries of the Delaunay circles (spheres in 3D) where they are C 1 , because of the discontinuity in the neighbor nodes across these boundaries. Another important property of this interpolation is its ability to reproduce linear functions on the boundary of convex domains. The proof can be found in Sukumar et al. (Sukumar et al., 1998), that we illustrate in figure 7(right) due to the fact that the Voronoi cells areas associated to points on the boundary become infinite, the contribution of internal points vanishes in the limit when the point approaches the convex boundary, and the shape functions associated with nodes n 1 and n 2 become linear on the segment (n 1 − n 2 ). This is not the case when non convex boundaries are considered where appropriate treatments are required (Cueto et al., 2000;Yvonnet et al., 2004). This drawback will be considered later. Consider an interpolation scheme for a vectorvalued function U h (x) : Ω ⊂ 2 → 2 , in the form where U i are the vectors defining the nodal degrees of freedom at the n natural neighbors of point x and ψ i (x) are the Sibson coordinates defined in [26] associated with each node n i . It can be noticed that [29] defines a local interpolation scheme which will be used to define both the trial and the test functions considered in the discretization of different variational formulations. A recent development in the NEM, the constrained natural element method (C-NEM), was proposed in (Yvonnet et al., 2005;Yvonnet et al., 2004) in order to circumvent the problems induced by non convex domains. In this approach, a visibility criterion is introduced to restrict the natural neighbors (influent nodes). For this purpose the constrained Voronoi diagram is constructed, from which the shape functions can be easily computed. In this manner, linear interpolation is recovered along the boundary of non-convex domains, making possible the introduction of essential boundary conditions as well as the treatment of fixed or moving discontinuities.
In the C-NEM framework, the interpolation can be expressed by where V is the number of natural neighbors visible from point x and ψ C i is the constrained natural neighbor shape function, which is actually the Sibson interpolant computed using the constrained Voronoi diagram (Yvonnet et al., 2004).

Transferring microscopic information towards the macroscopic scale
In this section we present the homogenization strategy. One of the contributions of the present work is the calculation of the thermomechanical homogenized properties by the application of the model reduction method presented in section 4.3. This method allows to use real boundary conditions to calculate the homogenized properties. The thermomechanical properties are then extended to the transition domain (where the macroscopic analysis is performed) by using an interpolation method. The transition medium has a particular role in the multi-scale model, because it ensures the connection between all the microscopic patterns.

Homogenization strategy
Let's Ω be the domain of interest, i.e. the whole macroscopic domain. We define in this macroscopic domain the elements Z i=1...N for which we will place a microscopic pattern (i.e. a representative volume element) in order to determine the microscopic state and the thermomechanical homogenized properties (see figure 8 left). An independent mesh Γ i=1...N is built for each macroscopic elements Z i=1...N with the constrained natural element method (C-NEM) (see figure 8 right Each Ω i is a microscopic pattern. These patterns are either surrounded by a transition medium Ξ i (see figure 8 right). The transfer between macroscopic scale and microscopic scale is realized by the four nodes of the macroscopic element Z i (see figure 8 right) and the transition domain Ξ i=1...N which surrounds the representative volume element (Ω i ) i=1...N . The thermomechanical properties are assumed to be perfectly defined at the microscopic scale inside each one of the microscopic patterns (Ω i ). In the transition domain we assume that the spatial variations of the averaged (or macroscopic) thermomechanical properties is smooth enough. In the next section we present a new method to realize the

Thermal conductivity homogenization
We assume known the temperature field T (X, x, t = t n + iΔt micro ), i ∈ [1, . . . , p], where X denotes the macroscopic position related to the microscopic domain (center of gravity of the microscopic cell) an x the point in the microscopic domain. We construct the matrix Q as presented in [15]. From the application of the model reduction method presented in section 4.3 we select the two most representative eigenfunctions φ 1 and φ 2 (related to the highest eigenvalues). Now, in a fist approximation the temperature evolution at the microscopic level can be approximated as where the β 1 and β 2 coefficients can be computed by minimizing the integral Now, with the beta coefficients just computed we can define the first and second temperature contributions Φ 1 and Φ 2 respectively, with Φ 1 = β 1 φ 1 and Φ 2 = β 2 φ 2 . We define the corresponding local gradients as well as the average ones The scale transfer tensor M (X, x) is then defined by and applying this expression to both gradient contributions, the four components of tensor M (X, x) at each point x ∈ Ω i can be computed Where (X) serves to locate Ω i in Ω. Finally, the homogenized conductivity tensor corresponding to the microscopic domain Ω i results The thermal conductivity at a point M in the transition domain Ξ i of the microscopic pattern Γ i is The simplest manner for defining the behavior at any point in the macroscopic domain Ω consists of using C-NEM interpolation constructed on the basis of the Voronoi diagram related to the cloud of points consisting of the center of gravity of the macroscopic elements Z i . Let's ψ i be the shape function related to the center of gravity of Z i which is perfectly defined in the whole domain Ω. Thus, the macroscopic law results Therefore the thermal conductivity is perfectly defined at each point M of the macroscopic domain Ω.
Remark 1 : if the macroscopic thermal gradient G 2 is too small, we construct an artificial G 2 in order to calculate the scale transfer tensor. Indeed [36] leads to The tensor G (X) must be invertible, i.e. det [G (X)] = 0 i.e. G 1X (X) ·G 2Y (X) = G 1Y (X) ·G 2X (X). For this purpose, we define the microscopic thermal gradient g 2 (x) on each point x of the RVE such as g 2 (x) ⊥g 1 (x). Thus we can calculate the transfer scale tensor M (X, x) and the homogenized conductivity with [37].

Remark 2 :
We can apply the model reduction directly on the microscopic thermal gradients. For this purpose, we construct the matrix of the microscopic and build the matrix c = Q.Q T . The application of the model reduction gives the eigenvectors φ 1 and φ 2 and therefore the homogenized conductivity with [37].

Remark 3 :
The same procedure can be used to compute the homogenized constitutive law. Obviously, in the linear case one could perform a standard homogenization by solving three linear boundary problems. However, the non-linear case could be addressed within the framework just described.

Stress calculation
The residual stress calculation requires the calculation of the thermal stress field and the crystallization induced stress field within the composite during the process. These stresses are induced during the cooling phase.

Thermal induced stress calculation
The macroscopic thermal stress at a point M of the macroscopic domain is defined by the relationship with C |M the macroscopic constitutive behavior at the point M and th |M the macroscopic thermal strain at a point M . The macroscopic thermal strain at a point M of the macroscopic domain is calculated by the relationship with α th |M is the thermal expansion at the point M , T |M the temperature at the point M and T ref the reference temperature for which there is no thermal expansion. The thermal expansion α th |M at the point M is calculated by the relationship The average thermal expansion on each microscopic domain Ω i is calculated by evaluating the average thermal strain on the microscopic domain. The average thermal strain on the microscopic domain Ω i can be expressed as with the microscopic thermal strain Moreover the average thermal strain can be expressed by And the average temperature variation on the microscopic domain Ω i can be expressed by Homogenized thermal expansion α th i and the thermal stress field can thus be calculated. The thermal stress field induces nodal forces in the macroscopic mesh Ω which can be explained by with B is the matrix which contains the shape function derivatives.

Crystallization induced stress calculation
The volume change ΔV c|M at a point M of the macroscopic domain Ω can be expressed by By considering an isotropic material, the volume change at the point M is and therefore the crystallization induced strain at the point M is with 1 is the identity tensor. Finally, the crystallization induced strain at the point M is At the microscopic level, only the matrix density varies. As the density of amorphous and crystalline phases is known, the actual density of a point P in the matrix is where ρ a is the density of the amorphous phase and ρ c the density of the crystalline phase. The average material density of a microscopic domain Ω i is calculated by and the macroscopic material density at each point M of the macroscopic domain Ω is calculated by The macroscopic induced crystallization stress at each point M of the macroscopic domain Ω is then defined by The induced crystallization stress field induces nodal forces in the macroscopic mesh Ω which can be explained by with B is the matrix which contains the shape functions derivatives.

Macroscopic stress calculation
In order to predict the residual stress field in the composite at the end of the process, we have to evaluate the macroscopic stress field in the structure. The global stress is calculated by the equation where |M is the macroscopic global strain, crist |M is the macroscopic crystallization induced strain and th |M is the macroscopic thermal induced strain.
[59] can be rewritten where σ crist |M is the macroscopic crystallization induced stress ([57]) and σ th |M is the macroscopic thermal induced stress ([41]). The macroscopic global strain |M is calculated by the relationship with U |M is the global displacement field at the point M , obtained by solving the equation

Numerical examples
In this section we present the results of the simulation of a cooling of two prepreg layers laid one on the other. We consider a graphite fiber reinforced composite similar to APC2 (Cogswell, 1992) (see table 1). The geometrical and boundary conditions of the simulation are given on figure 9. The results are presented for a time of 10 −1 s. The temperature field is presented on figure 10. The degree of crystallinity field is presented on figure 11. The thermal stress field is presented on figure 12. According to 41, the thermal stress components are negative. The crystallization stress field is presented on figure 13. We note that the values of the crystallization stress field are low. Residual stress field is then mainly induced by thermal stress field. Finally the global stress field [59] (σ Y Y component) is presented on figure 14. We see that in agreement with theory the top prepreg layer is in compression (σ Y Y < 0).

Conclusion
In this paper, we have proposed a coupled two-scales modelling for thermomechanical problems involving microstructure, which allows to evaluate residual stress growth during the on-line consolidation process. The method allows to take into account the thermomechanical properties evolution of the material during the process, which is a dominating factor in residual stress development. We use a meshless constrained natural element interpolation for taking into account the large and very localized variations in the nodal density. This technique allows considering simultaneously both scales, the one related to the microstructure description and the one related to its evolution (the macroscopic one). For this purpose the domain of interest is divided in some microscopic zones representing the microstructure evolution in the whole macroscopic domain, which do not cover the entire domain, and a transition domain. The thermomechanical properties are defined only in the microscopic zones at the microscopic scale. The thermomechanical properties are homogenized with the real boundary conditions by addressing proper orthogonal decomposition. The homogenized laws accurately defined at the microscopic level are then extended, using an appropriate technique, to that transition domain. In contrast to the vast majority of homogenization techniques, our approach allows an accurate description of the boundary conditions, because the microscopic domains can be located on the domain boundary. The prospects for the study are to realize a global simulation of the industrial process and then to compare the results of the numerical simulation with the industrial experimental results. This constitutes a work in progress.