Thermoplastic composites offer new possibilities for the industry. Large struc- tures can be processed rapidly and more cost-effectively than when thermoset composites are used, since the latter need to undergo lengthy curing reactions. The ability to fuse thermoplastic resins gives new perspectives for forming processes.
The thermostamping process is derived from the metallic materials industry. Forming occurs in two steps. In a first step, a semi-finished thermoplastic flat laminate, called the blank, is heated above the processing temperature of the matrix, usually using infra-red lamps. In the second step, this hot blank is quickly transferred to a cooled mould where it is stamped and given its final shape   . The heating and cooling steps are therefore separated. This results in an high production rate that makes this process very attractive for the industry.
Even though metal stamping has been the subject of extensive research work in the past decades (see for instance the review by Karbasian and Tekkaya  ), thermostamping of composite materials adds a new level of complexity for two resons. Indeed, the mechanical deformation and heat transfer occurring in the blanks may result in a complex and unexpected behaviour, especially when dealing with textile composite laminates. Nonetheless, accurate modelling and prediction of the main physical phenomena involved are prerequisite for an efficient process optimization.
1.2. Heat Transfer in Composite Stamping
It is well established that the temperature evolution is of major importance in this forming process. Keeping this in mind, de Luca et al.  or Cao et al.  proposed to take the blank temperature into account in the mechanical predic- tions of thermostamping process. Cao et al.  considered only two possible state: a high temperature state before the blank comes in contact with the mould, and a low temperature state after contact occurs. Based on previous work by Pickett et al.  , de Luca et al.  propose a modelling of the through thickness heat transfer using finite volume but are only able to predict the average tem- perature per ply in the case of a composite laminate. In thermostamping process thought, and especially during the stamping step, because of the thermal shock between the cold mould and the hot blank, high through-thickness temperature gradients may arise. The models by these authors, based on rough approxima- tions of the through thickness temperature profiles, cannot accurately describe these high through-thickness variations.
A finer through thickness temperature distribution description was proposed by Thomann et al.  using a finite difference method. Nonetheless they ne- glected the in-plane effects and thus considered only unidirectional through- thickness heat transfer. On the contrary, in real industrial processes, in-plane diffusion and 3D effects cannot be neglected, especially when boundary con- dition sharply evolve (in the vicinity of cavity edges) or in case of curved geome- tries. In the present paper, a fine description of the through thickness tempera- ture profile, in conjunction with the in-plane transfer is proposed.
Furthermore, the proposed model is designed to be easily implemented in any existing industrial code (such as Plasfib  , Aniform  or PAM-Form  ). The heat transfer problem should then be solved within acceptable computational times. With this aim, the full three-dimensional heat transfer problem cannot be solved using standard methods. Instead, a model reduction is necessary.
Considering the composite blank as a thin shell, it is natural to decompose the 3D temperature solution into a shape function and an in-plane temperature. As suggested by Saetta and Rega  , it writes
With this decomposition, the accuracy of through thickness description de- pends on the type of shape functions chosen. Within this framework, some authors suggested to construct new 3D shell finite elements that integrate this through thickness heat transfer effects     . Nonetheless, using one single shell element in the thickness highly restricts the possible through- thickness temperature profile description. Even with the parabolic shape pre- supposed by Alves Do Carmo and Rocha De Faria  or the higher order interpolation proposed by Surana and Abusaleh  , sharp profiles that arise in case of the thermal shocks that occur in thermostamping, will not be accurately described.
Adopting a fine through-thickness discretization therefore seems a more flexible approach, though potentially time-consuming. In this idea, Bognet et al.  wrote the above decomposition as a sum of separated modes
where the shape functions, themselves, are described with a fine discretization involving hundreds of degrees of freedom. In this framework, Bognet et al. considered a series of multiplicative shape functions, where each mode is the product of an out-of-plane function by an in-plane function. The out-of-plane function is therefore identical for all the points of the shell. Using this in-plane/ out-of-plane separation, a solving strategy using the proper generalized decom- position (PGD) was proposed for the elastic problem on a shell like domain. More recently, the  the method has been extended to nonlinear thermal problems. Though possibly efficient in some cases, such a resolution strategy in the environments of existing codes might be challenging. In particular, dealing with space varying boundary conditions and material non-linearity requires complicated developments and a probably a high number of modes.
1.3. Alternate Direction Implicit (ADI) Decomposition
In this paper, starting from a very general approximation framework as given by Equation (1), we propose a reduced numerical scheme, adapted to thin compo- site shells, that preserves the three-dimensional nature of the heat transfer problem but takes advantage of the good physical separation between in-plane and out-of-plane phenomena, even in case of anisotropic thermal properties. The present method is based on an operator splitting technique that enables to simplify a time evolution problem implying several spatial dimensions. The general framework of operator splitting techniques always considers an incre- mental iterative time integration strategy. Over 50 years ago, Douglas  and Douglas and Rachford  suggested to treat separately, within one time step, the different spatial directions. This led to the so called locally one-dimensional methods  or alternate direction implicit (ADI) methods. Then, numerous extension were proposed to reduce the error of the splitting strategy, and to validate the convergence and stability of the schemes, in linear and nonlinear cases     .
Following these ideas, the present paper proposes an operator splitting strategy adapted to the composite shell problems to solve the reduced heat trans- fer model. In fine, this results in two separated problems. A solving algorithm and numerical implementation is then proposed. The approach is validated on a flat plate test case, and its limits are determined with parametric studies. The method validity is extended to nonlinear cases with an industrial appli- cation.
2.1. Initial Heat Transfer Problem
The heat transfer problem is solved in the domain representing a composite laminate blank. It is considered to be an arbitrary curved thin shell, where the local positions are located via a curvilinear parallel coordinate system. A local frame can be attached to each point. Coordinate enables the location of points along the thickness direction, that is to say along, the normal vector to the shell mid-plane (see Figure 1). In this domain the compo- site material is considered to be a continuous medium with effective homogene- ous properties.
2.1.2. Heat Equation
In the considered heat transfer problem, the conduction is assumed to be governed by an anisotropic Fourier law where the local heat flux is written as:
where is the thermal conductivity tensor, the temperature field and the spatial derivative operator. In the present work, it is assumed that the through thickness direction is a principal direction of the thermal conduc- tivity. This is a classical assumption in the case of standard composite laminates   . Thus, in the basis, it writes
Figure 1. Shell like domain on which the heat transfer problem is solved. de- notes the out of plane direction and the thickness of the laminate. A typical in-plane dimension is.
being the in-plane thermal conductivity tensor and the through thickness thermal conductivity. Note that this hypothesis fails in the case of complex 3D architectured composites. Defining the in-plane surface gradient, Equation (2) can be separated into a through thickness and an in-plane fluxes:
In the case of a flat shell, the coordinate system is the natural cartesian coordinate system, and. In the more complex case of an arbitrary curved shell, the reader should refer to Appendix for a proper definition of the surface gradient. This demons- tration shows that in the case of a thin shell with small curvature, the operator does not depend on the through thickness position.
Using this separation, without internal heat source in the domain, the energy balance typically writes
being the density of the composite material and its specific heat. Once again, for a flat shell, the surface divergence, but for curved shell, it is defined in Appendix and it is constant through thickness.
2.1.3. Boundary and Initial Conditions
The domain is bounded by the boundaries, and, as defined in Figure 1. For the sake of simplicity, the lateral boundaries are considered insulated:
being the outward normal to each surface. Conversely, in order to accura- tely model temperature history imposed on the upper and lower boundaries and, a mixed boundary condition is assumed:
where (respectively) is the temperature imposed on the upper (re- spectively lower) boundary and (respectively) is the heat exchange coefficient. This mixed boundary condition modelling can account for non ideal contact with the mould   . In its limit form, it is also suited to model both Dirichlet or Neumann boundary conditions. Note that the development pro- posed hereunder could seemlessly be conducted with any type of boundary conditions (temperature imposed, heat flux, radiating surface...).
The initial temperature field, assumed given, is defined as:
2.2. Alternate Direction Implicit (ADI) Model
This section presents a reduction of the heat transfer problem defined above. The reduced boundary value problem is obtained thanks to an intuitive decom- position of the temperature field and a thin shell assumption. An implementa- tion strategy is then proposed to numerically solve this problem. Here, for the sake of clarity, the heat transfer problem is assumed linear (the material pro- perties, and do not depend on the temperature). The extension to nonlinear case will be discussed with a test case in Section 3.3.
2.2.1. Additive Decomposition
The first step in the proposed model reduction is to seek the solution of the system of Equations (5) to (8) as a sum of a through thickness averaged field and of a fluctuation field:
where the operator
is the through thickness average, being the local shell thickness. It is obvious that using this additive decomposition, the average field does not depend on the -coordinate whereas the fluctuation field has a zero thickness average. This decomposition is intuitive and does not necessitate any assump- tion. Substituting this decomposition (9) in the heat Equation (5), considering constant material properties, and noting that and the operator do not depend on the -coordinate gives:
Applying the average operator on both hands of this equation leads to
By defining the upper and lower inward boundary fluxes
Equation (12) writes:
which is the average field heat equation. It rules the in-plane mean field tem- perature evolution. Subtracting this mean heat equation from Equation (11) results in the fluctuating heat equation:
which rules the through thickness temperature fluctuation.
Assuming a thin plate for which, the so called aspect ratio for conduc- tion:
and the dimensional analysis safely leads to
Equation (15) then reduces to the fluctuating field heat equation:
Equations (14) and (18) achieve a decomposition of the initial heat Equation (5) in the average and fluctuating contributions. Nonetheless, without further assumptions, these two equations are strongly coupled through the source terms.
Reduced model. Summing Equations (14) and (18), and adding the term, gives
This equation, along with boundary and initial conditions (6), (7) and (8) defines the reduced boundary value problem (). In the bulk Equation (19), the first spatial differential operator of the right hand side acts on the average parts of the temperature field only. The solving of this reduced boundary value problem is therefore not straightforward. In the next section, a numerical method is proposed to solve this original model. It will also confirm its well-posedness.
2.2.2. Operator Splitting
Time discretization. The time evolution problem given by Equation (19) is solved in the framework of a standard incremental iterative time integration scheme. At a given time, the solution is supposed to be known. Then, the solution of the reduced boundary value problem defined above is searched at next time step.
Any conventional time integration scheme, such as for example explicit or implicit schemes, can be used to determined in terms of, so that the developments detailed hereunder will easily be implemented in such software environment.
• Step 1: solve the following 1D boundary value problem called () over one full time step
Figure 2. Operator splitting strategy. Instead of solving the full evolution equation on one time step, each differential operator is addressed successively. The initial condition of the second step is the field obtained at the end of the first step.
gives the intermediate result at the end of the time step.
• Step 2: solve the 2D boundary value problem over one full time step
where the initial condition is the value of the field computed in step 1. The solution of this second step at the end of the time step () is identified to.
Whereas the system () defined in step 1 is a well posed unidimensional partial differential equation, it is somewhat disturbing that both and appear in the problem (21) defined in step 2.
ADI model. To ensure the well-posedness of this step 2, the additive decom- position (9) is again substituted in system (21). Applying the average operator gives the in-plane boundary value problem ()
Finally, subtracting (22) from (21) results in
which admits the trivial constant solution:
Therefore, the fluctuating part of this second step is simply the fluctuating part of computed in the first step. In other terms, this second step does not introduce any additional out-of-plane fluctuation to the solution.
2.3. Numerical Algorithm
To ensure spatial numerical integration of this problems, a spatial discretization has to be adopted. Within the defined shell like domain a natural extruded discretization is assumed. Thus, and without loss of generality, for each in-plane discrete position amongst nodes, there is out of plane nodes. The dimension of the 3D discretized field is then.
Following the above additive decomposition and operator splitting strategy,
In this sum,
• is obtained by solving the fluctuation 1D boundary value problem () (Equation (20)). This problem is parametrized by the in-plane position through the dependency of the thickness and boundary conditions, , and. Thus, the problem () has to be solved times. Nonetheless each resolution has the complexity of a 1D boundary value pro- blem. Furthermore, each resolution is independent, and can be solved in a parallel manner as illustrated in Figure 3.
• is obtained as a post-processing by averaging the above field through thickness.
• is obtained by solving one single in plane 2D boundary value problem () (Equation (22)) using the 2D field as an initial condi- tion. At the end of time step, it gives the field.
Expected computational speed up. A conventional in plane discretization of an industrial geometry would typically result in nodes.
Figure 3. Resolution strategy. At each time step, the solution is obtained with two successive steps: solving a set of fluctuation problems () and solving one single in-plane problem ().
Additionally, because of the high through thickness temperature gradients associated with thermal shocks that occur in thermo-stamping, a fine through thickness discretization is required, for instance. In this case, the number of degree of freedom reaches.
Solving the initial 3D heat transfer problem defined in Section 2.1 using standard methods would result in solving a transient problem with degrees of freedom and a three-dimensional complexity. It would quickly result in unrealistic computational costs. Moreover, in the case of a thin shell domain, the proposed mesh, involving in plane nodes and through thickness nodes, would result in anisotropic mesh that may lead to numerical errors.
On the contrary, in the proposed resolution strategy, at each time step, independent 1D boundary value problem () with degrees of freedom can be solved in parallel, followed by one single 2D boundary value problem () with degrees of freedom. This strategy should result in a greatly reduced computational cost with a preserved accuracy, which opens the way for integrat- ing such approach as sub-routine in industrial simulation tools. Moreover, the in-plane and out-of-plane mesh sizes appear in two different problems and thus saves from complicated anisotropic meshing techniques.
Asynchronous time integration. Because of the thin plate assumption where, the ratio between characteristic in-plane diffusion time and charac- teristic through thickness diffusion time writes
being a dimensionless parameter characteristic of the so-called conduction aspect ratio. In a typical industrial case, where, , and, this ratio drops below. Therefore, the characteristic through thickness diffusion time is way shorter than its in-plane counterpart. This context justifies the use of an asynchronous time integration scheme, where two different time steps are used respectively for the through thickness fluc- tuating problem () and the in-plane problem ().
In practice, the global resolution algorithm presented in Figure 3 is kept, and the global time stepping is based on the in-plane requirements (). During each time step, the through-thickness problems are calculated by a sub-integration with smaller time steps of the order.
3. Results and Discussion
In this section, first, the proposed separated model and resolution strategy is validated on a test case that largely fulfills the thin shell assumption. Then the speed up is discussed and the limits of the presented model are investigated with rougher cases (thick and curved shell).
In order to validate the proposed resolution strategy, the temperature fields obtained using the presented model are compared with the temperature fields obtained by solving the initial three-dimensional problem, using a commercial software (COMSOL Multiphysics 5.0®).
3.1.1. Test Conditions
A square flat plate of dimensions and thickness is considered. The origin of the cartesian coordinate system is taken in the centre of the plate. In such a flat plate case, the curvilinear coordinates are identified to the cartesian ones: and.
Material properties. In this test case, a PA66/glass fibre composite material is considered. The homogenized material properties are adapted from the litera- ture  . The in plane conductivity is considered isotropic and all the material properties are supposed constant and are listed in Table 1.
Boundary and initial conditions. The boundary and initial conditions are given in Table 2. The plate is supposed to be initially at uniform room tempera- ture.
A different heating condition is imposed on the upper and lower surfaces with. It is representative of the temperature imposed by a hot mould con- tact. In order to add in-plane variability to the problem, the exchange coeffi- cients and artificially depend on space position through the characteristic gaussian function
The problem is solved on the time interval.
3.1.2. Numerical Parameters
Numerical methods. The 1D transient boundary value problems () and the 2D transient boundary value problem () are solved using a finite element method with piecewise linear interpolation. An implicit (backward Euler) time integration scheme is used for all time integrations. The proposed algorithm was programmed in MATLAB®, which enables the parallel resolution of the () problems.
Table 1. The material properties used in the test case are adapted from Faraj et al.  .
Table 2. Initial and boundary conditions used in the test case.
Mesh. For the reference simulation, a 3D regular mesh made of 3600 hexa- hedron is obtained by extruding a regular in-plane 2D mesh that consists of quadrangular elements. There are thus 30 elements in the thickness, and in terms of nodes, and.
For the proposed separated method, the mesh consists of the same 31 nodes through the thickness for the problems and of a triangular regular mesh with the same 3721 nodes for the problem.
The interpolations used in every finite element methods (3D in COMSOL, 2D in and 1D in) are all linear, which enables to expect for the same precision.
Time step. Time stepping in the FEM reference simulation follows the COMSOL built-in algorithm and is forced not to exceed. The time integration scheme is a standard backward difference scheme. On the contrary, a constant time step is used in the presented method. In this first test case, the time steps for both and problems are the same.
The convergence of the numerical methods used was first validated on a standard one-dimensional test case by comparing the numerical solution with an analytical solution given by Jaeger  .
Figure 4. Temperature profile at versus for three different heights in the plate. The plot is at final time. The reference 3D finite element solution (plain lines) is accurately recovered with the proposed methodology (markers).
Figure 5. Temperature profile at versus for two different in plane positions and two different instants. Once again, the 3D finite element solution (plain lines) is accurately recovered with the proposed methodology (markers).
profiles in the centre and on the edge of the plate at and final time. The figures show a good superposition of the reference field obtained with the finite element simulation and the one obtained with the presented method. The same numerical artifact (a slight oscillation) is found with both methods in the through thickness profile at early time (). This is due to the finite element and time discretization that fail to accurately predict thermal shocks. this artifact does not influence the later time predictions (see for instance Fachinotti and Bellet  regarding this issue).
The maximum residual relative error
is defined, where is the field computed with the 3D model in COMSOL and is the field computed with the presented method. At final time, the error does not exceed which represents around.
3.2. Efficiency and Model Limits
3.2.1. Speed up
The reference finite element simulation was computed in on a desktop computer (see Table 3). The solving time per time step was about. The separated form solution was computed on the same computer in no more than, with about per time step. This represents a speed up of over 28
Table 3. Computational speeds, the proposed method results in a speed up of over 28 times. In case of asynchronous time stepping and parallel resolution of problems, this speed up even reaches 33 times.
times. Using asynchronous time steps for and results in an additional reduction in the total computational time. Moreover, the problems are solved in a sequential manner in this test case. Solving them in parallel results in additional speed-up.
3.2.2. Extreme Cases
The limits of the proposed resolution strategy are investigated in this section. It is reminded that two conditions were required in the model development:
1) A small aspect ratio for conduction such that Equation (18) stands. This corresponds to the thin-shell assumption in the case where the in-plane and through thickness conductivities are of the same order of magnitude.
2) In the case of a curved shell domain, the radii of curvature should be large compared to the shell thickness. This ensures that the metrics given in the Appendix do not depend on the coordinate.
Thick part. In the test case presented above, the aspect ratio for conduction is very small () which explains the good app- licability of the thin plate assumption and the presented reduced method. The limit imposed by the first condition above was investigated by performing additional simulations with larger values of. With this aim, the plate dimen- sion was decreased. The plate is still flat and square. As shown in Figure 6 if stays smaller that, the thin plate assumption stands and the error given by (26) between the 3D finite element reference solution and the separated form solution does not exceed. It would even fall to lower than 1% for typical part shape encountered in composites processing ().
Sharp curvature. In order to investigate the curvature limit imposed by the second condition discussed above, a curved shell was considered. The domain is now an L-shape blank of length, with two flanges of length and a radius of curvature of the mid-plane surface. The blank thickness is kept (see Figure 7).
The boundary conditions on the upper and lower surfaces are now such that and. The reference field computed with the full 3D formulation using COMSOL Multiphysics® and the
Figure 6. Maximum error between the temperature fields computed with COMSOL using a 3D model and with the presented approach vs. aspect ratio for conduction. The error is computed at final time. As increases, the thin plate assumption fails, and the separated form resolution cannot predict 3D effects.
Figure 7. Geometry of the L-shape domain considered in the sharp curvature study. The sharpness of the curvature is given by the ratio between the radius of curvature and the flange thickness. The arrow represents the section along which the profile of Figure 9 is plotted.
As the blank thickness to radius of curvature ratio gets larger, the metric tensor given in Appendix by Equation (30) depends on the through thickness position. Thus Equation (31) does not stand and the proposed decomposition strategy fails at predicting the initial 3D problem. This is the case for Figure 8 where the thickness to radius ratio.
Figure 8. Through-thickness temperature profile at time obtained with the full 3D finite element solution and the proposed strategy. Case of a strong curvature:.
To identify the limit of applicability, several simulations with varying radius of curvature were performed. As shown in Figure 9 if stays below 0.2, which is usually the case in industrial geometries, the error between the 3D finite element reference solution and the separated form is below.
3.3. Application to Industrial Nonlinear Case
3.3.1. Problem Definition
The proposed ADI resolution method was applied to an industrial case representative of the thermostamping process. A thick laminate comprised of 16 anisotropic plies stacked on a sequence is considered. The temperature dependant thermal properties are adapted from carbon fibre reinforced PEEK and are given in Table 4. The initially hot laminate (at) comes in contact with a cold matrix and punch set, as illustrated in Figure 10, at time.
The 2D heat transfer problem is solved using (i) a full 2D resolution in COMSOL (ii) the presented alternate direction implicit (ADI) method, and (iii) a series of independant one-dimensional through thickness problems. In the ADI method, the problem consists of a 1D homogenized through thickness problem. Because of the stacking sequence, the in-plane thermal conductivity tensor is isotropic and is an average of the longitudinal and transsverse properties given in Table 4. Nonlinear resolution is performed in MATLAB over a physical time of with 150 time steps. In COMSOL, the exact multiply description is implemented. Using symmetry, only half of the geometry is considered and presented hereunder.
3.3.2. Results and Discussion
Three-dimensional effect. The problem is nonlinear, and, as visible in Figure 11, highly three-dimensional at the vicinity of the shear edge (between the
Figure 9. Maximum error between the temperature fields computed with COMSOL using a 3D model and with the presented approach vs. thickness to radius of curvature ratio. The error is computed at final time. As increases, the metric tensor becomes not constant through thickness, and the separated form resolution fails at predicting 3D effects.
Figure 10. Industrial test case geometry and boundary conditions. The problem is solved on the multiply laminate domain.
Figure 11. Industrial test case. Close up on temperature fields at time computed with the full 2D resolution (up) and the ADI method (down). The three-dimensional effect is partly described with the ADI method.
Table 4. Material properties used in the industrial case, representative of carbon fibre/ PEEK composite.
punch and the matrix). Still the proposed ADI method is able to partly discribe this tridimensional effect thanks to the problem that considers in plane diffusion.
Temperature profiles at three different positions at time are plotted in Figure 12.
• Far from the shear edge (cut CC'), the temperature gradient is mostly through thickness and the three approaches prove efficient at describing the through thickness temperature field.
• In the centre of the shear edge zone (cut AA'), the ADI method enables an accurate recovery of the through thickness profile obtained with the full 2D method. On the contrary, at this position AA', the one-dimensional method highly overestimates the temperature since it does not account for the nearby cold moulds.
• Similarly in the intermediate region over the matrix (cut BB’), the one- dimensional approach under predicts the temperature field. On the contrary, the ADI proposed method, enables a partial description of the three-dimensional effects (). Nonetheless, the method results in overpredicting temperature at height. At this worst position, three-dimensional effects are en- hanced, the decomposition methods fails and this artifact (also visible in Figure 11) cannot be overcome.
Nonlinearity. In addition to this three-dimensional effect, the proposed in- dustrial case is nonlinear, since the properties are temperature dependant. In this nonlinear case, the ADI method still proved efficient at predicting the temperature field. The efficiency of the method is explained by the very smooth non-linearities of the thermal properties used in the test case (see Table 4). Given that this is the case for the majority of industrial thermoplastic composite, the decomposition ADI method will likely be efficient for other industrial materials.
Multiply. Finally, the industrial test case was performed with a 16 plies laminates, with a very harsh anisotropic stacking. The ADI, which considers an homogenized in-plane conductivity for the in-plane problem, still proves efficient at predicting the thermal fields. In conclusion, as far as the heat transfer is concerned, a multiply stacking representative of an industrial blank can be considered as homogeneous through thickness.
Figure 12. Industrial test. Temperature profiles through thickness at three different positions. Plots are at time for the full 2D resolution, the ADI method and the one-dimensional method.
3.4. Proposed Integration in a Global Procedure
Several thermostamping simulation tools exist which handle the mechanics. This is the case, for instance, of Plasfib  , Aniform  or PAMForm  . In order to improve the physical description of these tools, accurate prediction of heat transfer is required. Implementation using the presented method, is possible for several reasons:
1) In these tools, the global time integration scheme is incremental and therefore follows the same scheme as the one described in Section 2.2.2. The iterative time integration procedure is thus consistent between the existing mechanical algorithm, and the proposed heat transfer with operator splitting algorithm.
2) The two-dimensional problem is a standard partial differential equa- tion. The spatial integration can be integrated using standard numerical methods. The FEM tools developed for the other physics (in the above examples, mechanics) can easily be reused for this heat equation.
3) The problems are independent one-dimensional partial differential equations. Implementation is straightforward using standard numerical methods (finite difference, finite elements).
4) The through thickness average two-dimensional temperature field is readily available as the solution of the problem at each time step. Thus it can be used as an input for a rough evaluation of a through thickness equivalent mechanical behavior. Furthermore, should one want a finer mechanical descrip- tion, the full three-dimensional field is also provided at each time step (Equation (25)).
An alternate direction implict (ADI) solving strategy was proposed to predict the temperature field in thin shells. It is particularly adapted to simulate temperature effects in thermo-stamping processes. The main contributions of this work are the following:
• An in-plane/out-of-plane decomposition strategy was proposed. The initial 3D heat transfer problem can be solved in two successive steps:
-solving of a series of 1D problems () with a fine time step and a good accounting of thermal shocks problems.
-solving of one 2D problem ().
The strong potential of this numerical strategy for computational costs reduc- tion was clearly highlighted.
• The applicability of this solving strategy was investigated. Two conditions are to be fulfilled for the model to be predictive:
-a small aspect ratio for conduction dimensionless ratio that includes both geometrical aspect ratio and anisotropy of the conductivity tensor.
-a small thickness to radius of curvature ratio.
These two conditions are fulfilled in standard thermo-stamping industrial cases.
• The proposed formulation is such that the problems and are classical and can be solved within a standard incremental time integration scheme and FEM formulations. Thus, the ADI decomposition can readily be implemented in existing industrial simulation softwares.
This study is part of the COMMANDO-STAMP project managed by IRT Jules Verne (French Institute in Research and Technology in Advanced Manufactur- ing Technologies for Composite, Metallic and Hybrid Structures). The authors wish to associate the industrial and academic partners of this project; Respec- tively SAFRAN, Peugeot Citroën Automotive, SOLVAY, CEMCAT, LTN, GeM, LAMCOS and 3SR. Also, fruitful discussions with Philippe Boisse and Nahiene Hamila about the integration in a global procedure are to be acknowledged.
Appendix. Arbitrary Curvilinear Shell Description
In this Appendix, the surface operator is defined in the arbitrary curved shell domain illustrated in Figure 1. This definition follows the framework adopted by Benveniste  in the case of a thin interphase. A similar approach is fully detailed in three dimensions by Bognet et al.  in their appendix.
A.1. Mid-Surface Description
Mapping. The reference global cartesian system is denoted as with its origin. First, the mid-surface of the shell like domain is considered. A position on this surface is parametrized in a reference dimensionless coordinate system using the mapping function
This mapping is such that and are dimensionless.
Basis. The natural basis at point consists of the two tangent vectors
where the standard comma notation denotes derivation.
Metric tensor. The first fundamental metric tensor of this 2D surface writes, in the local basis,
The unit normal to the tangent surface at point is also defined as
The second order tensor, representing the second fundamental form, which components are defined as, gives the local mean curvature
and Gaussian curvature
of the surface.
A.2. Shell Domain Parametrization
Mapping. A position in the thin shell domain is parametrized as described in Figure A1. The projection of on the mid-surface is first defined. Therefore,. is parametrized using the map- ping (27) and the third dimensionless coordinate
is defined, where is the local thin shell thickness. It quantifies the distance between and the mid-surface. Thus the coordinate is also dimensionless. In summary, the shell domain mapping writes
Because is the distance to the mid-surface, the coordinate sys- tem is parallel curvilinear as defined by Benveniste  .
Basis. At point, the natural basis associated to this curvilinear coordinate system consists of the three vectors
Figure A1. Illustration of the mapping used to parametrize the shell domain. The position in the physical Euclidean space is obtained from the dimen- sionless coordinates using: (i) the mid-surface mapping (the function) and (ii) the distance from the mid-surface.
Metric tensor. The symmetric fundamental metric tensor is described in terms of its coordinates where and. By defini- tion
Because the system is parallel curvilinear,. Moreover, Equation (29) gives
because is a unit vector.
Following Equation (64) in  , the component, and write1
1The expression (30) for differs from that of  because, in our case, , where depends on the coordinates and.
where the second order tensor represents the extrinsic third fundamental form. This equation shows that the local metric tensor is obtained as a correction of the metric tensor at the mid surface. This correction depends on the distance from and gets larger as the curvature increases. But it also depends on the shell thickness variation that may occur in the case of blanks with ply drops.
In the case where the radii of curvature of the surface are large compared to the shell thickness, the second term is negligible compared to the first one. Because is a product including the curvature (see for instance Equation (59) by Bognet et al.  ), the third term also vanishes besides when the curvature of tends to 0. If, in addition, the shell thickness variations are small, the last term can also be neglected. Then, the fundamental metric tensor in the shell reduces to
and is thus independent of the through thickness position in the shell. Furthermore, the inverse of this metric tensor is also block-diagonal and writes
A.3. Surface Differential Operators
Gradient. Following  , the gradient of a scalar is a first order tensor. In the contravariant basis, it writes
which can be decomposed, using Equation (32) into an in-plane and an out-of- plane term:
where the surface gradient writes, in the basis (,):
In the case where the out-of plane coordinate is redimensionalized, as, the normal vector, and
As described in section 2.1.2, for a conductivity tensor which has a principal direction in the out of plane direction (Equation (3)), the flux in-plane/out-of- plane decomposition (4) is recovered.
Divergence. First, the following scalar magnitude is defined:
The determinant of is thus
Following  , the divergence of a vector writes
where the Einstein summation notation is used on the index. Since does not depend on, this sum can be decomposed into in-plane and an out-of-plane terms:
where the surface divergence writes, in the basis (,):
In the case where the out-of plane coordinate is redimensionalized, as, and. The divergence then writes
As given in Section 2.1.2, the heat equation decomposition (5) is recovered.