On the Alternate Direction Implicit (ADI) Method for Solving Heat Transfer in Composite Stamping

Thermostamping of thermoplastic matrix composites is a process where a preheated blank is rapidly shaped in a cold matching mould. Predictive modelling of the main physical phenomena occurring in this process requires an accurate prediction of the temperature field. In this paper, a numerical method is proposed to simulate this heat transfer. The initial three-dimensional heat equation is handled using an additive decomposition, a thin shell assumption, and an operator splitting strategy. An adapted resolution algorithm is then presented. It results in an alternate direction implicit decomposition: the problem is solved successively as a 2D surface problem and several onedimensional through thickness problems. The strategy was fully validated versus a 3D calculation on a simple test case and the proposed strategy is shown to enable a tremendous calculation speed up. The limits of applicability of this method are investigated with two parametric studies, one on the thickness to width ratio and the other one on the effect of curvature. These conditions are usually fulfilled in industrial cases. Finally, even though the method was developed under linear assumption (constant material properties), the strategy validity is extended to multiply, temperature dependant (nonlinear) case using an industrial test case. Because of the standard methods involved, the proposed ADI method can readily be implemented in existing software.


Introduction 1.Context
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 [1] [2].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 [3]), 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.

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. [4] or Cao et al. [5] proposed to take the blank temperature into account in the mechanical predictions of thermostamping process.Cao et al. [5] 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. [6], de Luca et al. [4] propose a modelling of the through thickness heat transfer using finite volume but are only able to predict the average temperature 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 approximations 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. [7] using a finite difference method.Nonetheless they neglected the in-plane effects and thus considered only unidirectional throughthickness heat transfer.On the contrary, in real industrial processes, in-plane diffusion and 3D effects cannot be neglected, especially when boundary condition sharply evolve (in the vicinity of cavity edges) or in case of curved geometries.In the present paper, a fine description of the through thickness temperature 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 [8], Aniform [9] or PAM-Form [4]).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 [10], it writes With this decomposition, the accuracy of through thickness description depends 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 [11] [12] [13] [14].Nonetheless, using one single shell element in the thickness highly restricts the possible throughthickness temperature profile description.Even with the parabolic shape presupposed by Alves Do Carmo and Rocha De Faria [15] or the higher order interpolation proposed by Surana and Abusaleh [13], 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. [16] 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 i 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 decomposition (PGD) was proposed for the elastic problem on a shell like domain.
More recently, the [17] 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.

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 composite 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 incremental iterative time integration strategy.Over 50 years ago, Douglas [18] and Douglas and Rachford [19] suggested to treat separately, within one time step, the different spatial directions.This led to the so called locally one-dimensional methods [20] 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 [21] [22] [23] [24].
Following these ideas, the present paper proposes an operator splitting strategy adapted to the composite shell problems to solve the reduced heat transfer 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 application.

Domain
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 ( ) e e e can be attached to each point.Coordinate z enables the location of points along the thickness direction, that is to say along z e , the normal vector to the shell mid-plane (see Figure 1).In this domain the composite material is considered to be a continuous medium with effective homogeneous properties.

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 q is written as: where K is the thermal conductivity tensor, T the temperature field and ∇ ⋅ the spatial derivative operator.In the present work, it is assumed that the through thickness direction z is a principal direction of the thermal conduc- tivity.This is a classical assumption in the case of standard composite laminates [10] [25].Thus, in the ( ) e e e basis, it writes  2) can be separated into a through thickness and an in-plane fluxes: .
In the case of a flat shell Ω , the coordinate system ( ) , , p q z is the natural cartesian coordinate system ( ) , , X Y Z , 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 s ∇ ⋅.This demons- tration shows that in the case of a thin shell with small curvature, the operator s ∇ ⋅ does not depend on the through thickness position z .Using this separation, without internal heat source in the domain Ω , the energy balance typically writes ( ) ρ being the density of the composite material and p c its specific heat.Once again, for a flat shell, the surface divergence s curved shell, it is defined in Appendix and it is constant through thickness.

Boundary and Initial Conditions
The domain Ω is bounded by the boundaries lat Γ , sup Γ and inf Γ , as defined in Figure 1.For the sake of simplicity, the lateral boundaries are considered insulated: n being the outward normal to each surface.Conversely, in order to accura- tely model temperature history imposed on the upper and lower boundaries sup Γ and inf Γ , a mixed boundary condition is assumed: where The initial temperature field, assumed given, is defined as:

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 decomposition of the temperature field and a thin shell assumption.An implementation strategy is then proposed to numerically solve this problem.Here, for the sake of clarity, the heat transfer problem is assumed linear (the material properties ρ , p c and K do not depend on the temperature T ).The extension to nonlinear case will be discussed with a test case in Section 3.3.

Additive Decomposition
The first step in the proposed model reduction is to seek the solution T of the system of Equations ( 5) to (8) as a sum of a through thickness averaged field and of a fluctuation field:

x y z t T x y t T x y z t
where the operator is the through thickness average, h being the local shell thickness.It is obvious that using this additive decomposition, the average field z T does not depend on the z -coordinate whereas the fluctuation field T  has a zero thickness average.This decomposition is intuitive and does not necessitate any assumption.Substituting this decomposition (9) in the heat Equation ( 5), considering constant material properties, and noting that z T and the operator s ∇ do not depend on the z -coordinate gives: Applying the average operator z ⋅ 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 temperature 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 h L  , 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 ( ) 14) and (18), and adding the term This equation, along with boundary and initial conditions ( 6), ( 7) and (8) defines the reduced boundary value problem ( P ).In the bulk Equation ( 19), the first spatial differential operator of the right hand side acts on the average parts of the temperature field z T 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.

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 n t , the solution ( ) Any conventional time integration scheme, such as for example explicit or implicit schemes, can be used to determined , , n T x y z , so that the developments detailed hereunder will easily be implemented in such software environment.
Operator splitting.To solve Equation ( 19), an operator splitting method is used.This numerical method enables to solve evolution equations that involve a sum of differential operators (see for example [20]).Adopting the splitting initially suggested by Douglas [18] and later called locally one-dimensional (LOD) method (see for instance [26] and references therein), the two differential operators in the right hand side of Equation ( 19) are considered separately.Note that in this linear case, the proposed splitting does not introduce additional numerical error beside the time integration error [20].As illustrated in Figure 2, a so-called fractional time step method is adopted, where two problems are solved successively, each one containing one of the operators:  gives the intermediate result ( ) where the initial condition 1 2 n T + is the value of the field computed in step 1.
The solution of this second step at the end of the time step ( 1 d Whereas the system ( z P ) defined in step .0 on at .
Therefore, the fluctuating part T  of this second step is simply the fluctuating part of 1 2 n T + computed in the first step.In other terms, this second step does not introduce any additional out-of-plane fluctuation to the solution.

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 ( ) The dimension of the 3D discretized field is then s z N N × .
Resolution scheme.
n n n n z z T p q z n T p q z T p q T p q z T p q In this sum, • 1 2 n T + is obtained by solving the fluctuation 1D boundary value problem ( z P ) (Equation ( 20)).This problem is parametrized by the in-plane position ( ) , p q through the dependency of the thickness h and boundary conditions Nonetheless each resolution has the complexity of a 1D boundary value problem.Furthermore, each resolution is independent, and can be solved in a parallel manner as illustrated in Figure 3.
T + is obtained as a post-processing by averaging the above is obtained by solving one single in plane 2D boundary value problem ( m P ) (Equation ( 22)) using the 2D field 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 50 z N ∼ .In this case, the number of degree of freedom reaches 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 s N in plane nodes and z N 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, s N independent 1D boundary value problem ( z P ) with z N degrees of freedom can be solved in parallel, followed by one single 2D boundary value problem ( m P ) with s N degrees of freedom.This strategy should result in a greatly reduced computational cost with a preserved accuracy, which opens the way for integrating 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 h L  , the ratio between characteristic in-plane diffusion time p t and charac- teristic through thickness diffusion time z t writes , this ratio drops below 3 10 − .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 fluctuating problem ( z P ) and the in-plane problem ( m P ).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 dt , the through-thickness problems are calculated by a sub-integration with smaller time steps dt′ of the order ( )

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).

Validation
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 ® ).

Test Conditions
A square flat plate of dimensions 2 2 0.1 0.1 m L = × and thickness 5 mm h = is considered.The origin of the ( ) , , x y z 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: p x = and q y = .
Material properties.In this test case, a PA66/glass fibre composite material is considered.The homogenized material properties are adapted from the literature [27].The in plane conductivity s K 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 temperature init 20 C T =  .A different heating condition is imposed on the upper and lower surfaces with . It is representative of the temperature imposed by a hot mould contact.In order to add in-plane variability to the problem, the exchange coefficients sup h and inf h artificially depend on space position ( ) The problem is solved on the time interval [ ] 0, 20 s t = .

Numerical Parameters
Numerical methods.The 1D transient boundary value problems ( z P ) and the 2D transient boundary value problem ( m P ) 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 ( z P ) problems.
Table 1.The material properties used in the test case are adapted from Faraj et al. [27].For the proposed separated method, the mesh consists of the same 31 nodes through the thickness for the z P problems and of a triangular regular mesh with the same 3721 nodes for the m P problem.
The interpolations used in every finite element methods (3D in COMSOL, 2D in m P and 1D in z P ) 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 0.01 s .The time integration scheme is a standard backward difference scheme.On the contrary, a constant time step d 0.01 s t = is used in the presented method.In this first test case, the time steps for both z P and m P 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 [28].

Comparison
Figure 4 shows the in-plane temperature profiles at three different heights, at final time 20 s t = . Figure 5 represents the through thickness temperature

s t =
).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 [29] regarding this issue).
The maximum residual relative error is defined, where 3D T is the field computed with the 3D model in COMSOL and T is the field computed with the presented method.At final time 20 s t = , the error err does not exceed 2.5% which represents around 6.5 C  .

Speed up
The reference finite element simulation was computed in 10000 s on a desktop computer (see Table 3).The solving time per time step was about 5 s .The separated form solution was computed on the same computer in no more than 356 s , with about 0.178 s per time step.This represents a speed up of over 28

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 ( ) ( ) 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 h .This ensures that the metrics g given in the Appendix do not depend on the z hr = coordinate.
Thick part.In the test case presented above, the aspect ratio for conduction ( ) ( ) ) which explains the good applicability 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 A .With this aim, the plate dimension L was decreased.The plate is still flat and square.As shown in Figure 6 if A stays smaller that 0.01, 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 5% .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 is kept (see Figure 7).
The boundary conditions on the upper and lower surfaces are now such that  .The through-thickness profile along the first diagonal schematized in Figure 7 is plotted in Figure 8.
As the blank thickness to radius of curvature h R ratio gets larger, the metric tensor g given in Appendix by Equation ( 30) depends on the through thickness position z .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 1 h R = .1 h R = .
To identify the limit of applicability, several simulations with varying radius of curvature R were performed.As shown in Figure 9 if h R stays below 0.2, which is usually the case in industrial geometries, the error err between the 3D finite element reference solution and the separated form is below 0.3% .

Problem Definition
The proposed ADI resolution method was applied to an industrial case representative of the thermostamping process.A 2 mm thick laminate comprised of 16 anisotropic plies stacked on a 0,90 The temperature dependant thermal properties are adapted from carbon fibre reinforced PEEK and are given in Table 4.The initially hot laminate (at 400 C  ) comes in contact with a cold matrix and punch set, as illustrated in Figure 10, at time 0 s t = .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 m P problem consists of a 1D homogenized through thickness problem.Because of the 0,90 stacking sequence, the in-plane thermal conductivity tensor s K 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 5 s with 150 time steps.In COMSOL, the exact multiply description is implemented.Using symmetry, only half of the geometry is considered and presented hereunder.

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 h R .The error is computed at final time 20 s t = . As h R increases, the metric tensor g 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 3 s 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.
Transverse thermal conductivity  punch and the matrix).Still the proposed ADI method is able to partly discribe this tridimensional effect thanks to the m P problem that considers in plane diffusion.
Temperature profiles at three different positions at time • 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 onedimensional approach under predicts the temperature field.On the contrary, the ADI proposed method, enables a partial description of the three-dimensional effects ( 5 C ±  ).Nonetheless, the method results in overpredicting temperature at height 1 mm y = . At this worst position, three-dimensional effects are enhanced, 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 industrial 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 0 ,90 anisotropic stacking.The ADI, which considers an homogenized in-plane conductivity for the in-plane m P 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.

Proposed Integration in a Global Procedure
Several thermostamping simulation tools exist which handle the mechanics.This is the case, for instance, of Plasfib [30], Aniform [31] or PAMForm [4].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 m P 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 z P 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 z T is readily available as the solution of the m P 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 description, the full three-dimensional field is also provided at each time step (Equation ( 25)).

Conclusions
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 ( z P ) with a fine time step and a good accounting of thermal shocks problems.
-solving of one 2D problem ( m P ).
The strong potential of this numerical strategy for computational costs reduction 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 A that includes both geometrical aspect ratio h L and anisotropy of the conductivity tensor.
-a small thickness to radius of curvature ratio h R .
These two conditions are fulfilled in standard thermo-stamping industrial cases.
• The proposed formulation is such that the problems m P and z P 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.

Appendix. Arbitrary Curvilinear Shell Description
In this Appendix, the surface operator s ∇ is defined in the arbitrary curved shell domain illustrated in Figure 1.This definition follows the framework adopted by Benveniste [32] in the case of a thin interphase.A similar approach is fully detailed in three dimensions by Bognet et al. [33] in their appendix.

A.1. Mid-Surface Description
Mapping.The reference global cartesian system is denoted as ( ) , , X Y Z with its origin O .First, the mid-surface Γ of the shell like domain is considered.A position ( ) , , G X Y Z = on this surface Γ is parametrized in a reference dimensionless coordinate system ( ) , p q using the mapping function [ ] ( ) ( ) This mapping φ is such that p and q are dimensionless.
Basis.The natural basis at point G consists of the two tangent vectors

e e e g e e e e
The unit normal to the tangent surface at point G is also defined as 0 0 0 0 .
p q p q ∧ = ∧ e e n e e (28) The second order tensor b , representing the second fundamental form, which components are defined as Because r is the distance to the mid-surface, the ( ) , , p q r coordinate sys- tem is parallel curvilinear as defined by Benveniste [32].
Basis.At point M , the natural basis associated to this curvilinear coordinate system consists of the three vectors , , p q r using: (i) the mid-surface Γ mapping (the function p p r p p q q q q r q q r r r h r r h r h where the second order tensor c represents the extrinsic third fundamental form.This equation shows that the local metric tensor g is obtained as a correction of the metric tensor 0 g at the mid surface Γ .This correction depends on the distance rh from Γ and gets larger as the curvature b increases.But it also depends on the shell thickness variation ( )

2
,ij h 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 h , the second term is negligible compared to the first one.
Because ij c is a product including the curvature ij b (see for instance Equation (59) by Bognet et al. [33]), the third term also vanishes besides 0 ij g when the curvature of Γ tends to 0. If, in addition, the shell thickness variations ,i h are small, the last term can also be neglected.Then, the fundamental metric tensor  (31) and is thus independent of the through thickness position r in the shell.Furthermore, the inverse of this metric tensor is also block-diagonal and writes

A.3. Surface Differential Operators
Gradient.Following [34], the gradient of a scalar β is a first order tensor. 1 The expression (30) for ij g differs from that of [33] because, in our case, r h = e n, where h de- pends on the coordinates p and q .
In the contravariant basis ( ) e e e , it writes ( ) In the case where the out-of plane coordinate r is redimensionalized, as z hr = , the normal vector z r e h = = e n, and Divergence.First, the following scalar magnitude is defined: The determinant of g is thus ( ) 2 0 det .h g = g Following [34], the divergence of a vector ( ) where the Einstein summation notation is used on the index k .Since 0 h g does not depend on r , this sum can be decomposed into in-plane and an out-of-plane terms: where the surface divergence s ∇ ⋅v writes, in the basis ( p e , q e ): ( ) ( ) s p q p q v h g v h g h g h g In the case where the out-of plane coordinate r is redimensionalized, as z hr = Submit or recommend next manuscript to SCIRP and we will provide best service for you: Accepting pre-submission inquiries through Email, Facebook, LinkedIn, Twitter, etc.A wide selection of journals (inclusive of 9 subjects, more than 200 journals) Providing 24-hour high-quality service User-friendly online submission system Fair and swift peer-review system Efficient typesetting and proofreading procedure Display of the result of downloads and visits, as well as the number of cited articles Maximum dissemination of your research work Submit your manuscript at: http://papersubmission.scirp.org/Or contact msa@scirp.org

Figure 1 .
Figure 1.Shell like domain Ω on which the heat transfer problem is solved.z e de- notes the out of plane direction and h the thickness of the laminate.A typical in-plane dimension is 1 m L .
Step 1: solve the following 1D boundary value problem called ( z P ) over one full time step dt

Figure 2 .
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.

T
. Thus, the problem ( z P ) has to be solved s N times.

Figure 3 .
Figure 3. Resolution strategy.At each time step, the solution is obtained with two successive steps: solving a set of s N fluctuation problems ( z P ) and solving one single 3D heat transfer problem defined in Section 2.1 using standard methods would result in solving a transient problem with

A
being a dimensionless parameter characteristic of the so-called conduction aspect ratio.In a typical industrial case, where

Table 2 .
Initial and boundary conditions used in the test case.Mesh.For the reference simulation, a 3D regular mesh made of 3600 hexahedron is obtained by extruding a regular in-plane 2D mesh that consists of 60 60 × quadrangular elements.There are thus 30 elements in the thickness, and in terms of nodes,

Figure 4 .
Figure 4. Temperature profile at 0 y = versus x for three different heights z in the plate.The plot is at final time 20 s t = .The reference 3D finite element solution (plain lines) is accurately recovered with the proposed methodology (markers).

Figure 5 ..
Figure 5. Temperature profile at 0 y = versus z for two different in plane positionsx and two different instants t .Once again, the 3D finite element solution (plain lines) is accurately recovered with the proposed methodology (markers).
second condition discussed above, a curved shell was considered.The domain Ω is now an L-shape blank of length 0.05 m L = , with two flanges of length 0.1 m f L = and a radius of curvature of the mid-plane surface 0.005 m R = .The blank thickness 0.005 m h =

.Figure 6 .
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 A .The error is computed at final time 20 s t = .As A increases, the thin plate assumption fails, and the separated form resolution cannot predict 3D effects.

Figure 7 .
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 R and the flange thickness h .The arrow represents the section along which the of Figure 9 is plotted.

Figure 8 .
Figure 8. Through-thickness temperature profile at time 20 s t = obtained with the full 3D finite element solution and the proposed strategy.Case of a strong curvature: 1 h R = .

Figure 12 .
Figure 12.Industrial test.Temperature profiles through thickness at three different positions.Plots are at time 3 s for the full 2D resolution, the ADI method and the one-dimensional method.
of the Society for Industrial and Applied Mathematics, 3, 42-65.[19] Douglas, J. and Rachford, H.H. (1956) On the Numerical Solution of Heat Conduc- comma notation denotes derivation.Metric tensor.The first fundamental metric tensor of this 2D surface writes, in the local basis,( the surface Γ .A.2.Shell Domain ParametrizationMapping.A position M in the thin shell domain Ω is parametrized as described in FigureA1.The projection G of M on the mid-surface Γ is first defined.Therefore, OM OG GM = +   .G is parametrized using the map- ping(27) and the third dimensionless coordinate where h is the local thin shell thickness.It quantifies the distance between M and the mid-surface Γ .Thus the coordinate dimensionless.In summary, the shell domain mapping writes [ ]

Figure A1 .
Figure A1.Illustration of the mapping used to parametrize the shell domain Ω .The position ( ) , , M X Y Z = in the physical Euclidean space is obtained from the dimensionless coordinates ( ) ϕ → ) and (ii) the distance from the mid-surface hr .

(
decomposed, using Equation (32) into an in-plane and an out-of- As described in section for a conductivity tensor which has a principal direction in the out of plane direction (Equation (3)), the flux in-plane/out-ofplane decomposition (4) is recovered.
As given in Section 2.1.2,the heat equation decomposition (5) is recovered.
1 is a well posed unidimensional partial differential equation, it is somewhat disturbing that both T and z T appear in the problem (21) defined in step 2. ADI model.To ensure the well-posedness of this step 2, the additive decomposition (9) is again substituted in system (21).Applying the average operator z ⋅ gives the in-plane boundary value problem ( m P ) ( ) lat 1 2

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 z P problems, this speed up even reaches 33 times.