Yancheng Zhang
email: yancheng.zhang@mines-paristech.fr

Gildas Guillemot

Marc Bernacki

Michel Bellet

Macroscopic thermal finite element modelling of additive metal manufacturing by selective laser melting process

Keywords: Selective laser melting (SLM), Macroscopic heat transfer, Finite element, Level set

come

Macroscopic thermal finite element modelling of additive metal manufacturing by selective laser melting process

Introduction

Selective laser melting (SLM) is an advanced form of the selective laser sintering (SLS) process with full melting of the powder bed particles taking place by using one or more lasers. It was developed in the late 90's by Fockele and Schwarze in cooperation with the Fraunhofer Institute and then commercialized in 2004 by the company MCP HEK GmbH (now SLM Solutions GmbH) [START_REF] Bhavar | A review on powder bed fusion technology of metal additive manufacturing[END_REF].

SLM processes have received great interest for both research and industry aspects. However, the high level of residual stress induced by the large thermal gradients often leads to cracking phenomena in the part or the supporting structures during construction, or delamination from the build plate. Different studies aiming at avoiding such problems have been carried out for varied metallic alloys [START_REF] Thijs | A study of the microstructural evolution during selective laser melting of Ti-6Al-4V[END_REF][START_REF] Qiu | Microstructure and tensile properties of selectively laser-melted and of HIPed laser-melted Ti-6Al-4V[END_REF][START_REF] Zadpoor | Mechanics of additively manufactured biomaterials[END_REF][START_REF] Haigh | Additive manufacturing with polypropylene microfibers[END_REF].

In order to minimize costly trial and error approach by repeated experiments, numerical simulation has been introduced and continuously developed to model the thermal and mechanical responses during and after processing. Besides, end-use mechanical properties are also a point of interest and determined by the final microstructure and the porosities of pieces [START_REF] Martukanitz | Toward an integrated computational system for describing the additive manufacturing process for metallic materials[END_REF][START_REF] Denlinger | Thermal modelling of Inconel 718 processed with powder bed fusion and experimental validation using in situ measurements[END_REF]. Generally speaking, studies on numerical modelling work of SLM as reported in the literature are actually developed at various scales and usually considering several physical phenomena, leading to complex models. Recently, Megahed et al. [START_REF] Megahed | Metal additive-manufacturing process and residual stress modeling[END_REF] discussed the models required to span the scope of AM processes with the focus to predict built material characteristics and residual stresses of produced parts. The different models and numerical methods can be classified according to their characteristic scales, as detailed below.

For approaches at the powder particles scale, Steuben et al. [START_REF] Steuben | Discrete element modeling of particle-based additive manufacturing processes[END_REF] present the development of a methodology with the discrete element method extended to incorporate the physics of laser heating, leading to a computationally efficient approach to simulate SLS of metals at relatively large scale. A 3D mesoscopic powder model has also been developed using a hybrid finite element (FE) and finite volume formulation on unstructured meshes by Khairallah and Anderson [START_REF] Khairallah | Mesoscopic simulation model of selective laser melting of stainless steel powder[END_REF].

Concerning the material deposit scale, the powder bed is mostly assumed as a continuum medium enduring continuous evolution in order to predict its melting/solidification, the formation of the melt pool and the resulting shapes of the deposited tracks [START_REF] Chen | Three-dimensional finite element thermomechanical modeling of additive manufacturing by selective laser melting for ceramic materials[END_REF]. A thermal resolution is developed to predict the temperature field evolution. As a complementary part, the Navier-Stokes equations must be solved. Due to the process velocity and considering the rapid material evolution, the time step is in the level of µs. As a result, predictions on the scale of the constructed workpiece are totally unreachable by the two previous approaches due to prohibitive computational times.

Regarding the workpiece scale, macroscopic approaches are the most appropriate and the 3D finite element method (FEM) is usually adopted. In the work of Foroozmehr et al. [START_REF] Foroozmehr | Finite Element Simulation of Selective Laser Melting process considering Optical Penetration Depth of laser in powder bed[END_REF], the model of heat source takes into account the optical penetration depth of the laser beam into the powder bed and its dependency with the powder size. The melt pool depth prediction is validated by experiments.

Investigations concerning residual stress were also performed during the SLM simulation process by Gu et He [START_REF] Gu | Finite element simulation and experimental investigation of residual stresses in selective laser melted Ti-Ni shape memory alloy[END_REF]. They found that the maximum residual stresses exist at the end of the first track and last track, which was also confirmed by experimental results. In order to consider the microstructure effect, coupled constitutive and microstructure models are part of the simulation model of Lindgren et al. [START_REF] Lindgren | Simulation of additive manufacturing using coupled constitutive and microstructure models[END_REF], in which precipitation hardening and phase change are also included.

For such macroscale modelling, regular or simple shape of the workpiece is usually considered when the powder bed is taken into account [START_REF] Denlinger | Thermal modelling of Inconel 718 processed with powder bed fusion and experimental validation using in situ measurements[END_REF][START_REF] Riedlbauer | Macroscopic simulation and experimental measurement of melt pool characteristics in selective electron beam melting of Ti-6Al-4V[END_REF], while the powder bed is rarely considered when complex workpieces are addressed due to the difficulties to consider part/powder interfaces in such practical cases.

The objective of this paper is to present an innovative FEM thermal model. Its main features are: i) the deposition of energy and matter by layer or by fraction of layer, not following continuously and explicitly the laser scan path in time and space; ii) the numerical heat transfer modelling through the consolidated piece as well as non-melted powder bed; and iii) the procedures implemented for mesh metric calculation and time step adaptation. Numerical validation is presented and an application example with complex geometry is given in the last section. A second level set function, , is then defined to divide in two regions: during the process, the construction front, i.e., the interface between material (solid workpiece and remaining powder) and gas domains is identified thanks to the zero isovalue of the signed distance to this horizontal surface (red line in Fig. 1). The elaborated workpiece and the non-melted powder are found in the zone , while the gas occupies the zone . As a consequence, the workpiece under construction is split into a partially elaborated workpiece, , and a virtual workpiece (still to be built), . The complement to is the non-melted powder bed, . Therefore, in the material domain , the level set defines the interface between solid workpiece and powder, . Finally, the entire domain can be expressed as , where = presents the material domain containing the remaining powder, solid workpiece under construction and the substrate.

Modelling strategy

G-code interpreter module

Based on CAD information, the laser scan path is usually generated in the form of a G-code or other specific machine language file, by a software tool named as a "slicer". The slicer intersects the CAD model with a series of planes in order to define a set of cross sections. For each plane, both perimeter and infill laser beam trajectories are computed to produce a layer [START_REF] Steuben | Discrete element modeling of particle-based additive manufacturing processes[END_REF]. In the present work, the discretized laser beam trajectories are extracted from the G-code file. Then each trajectory is interpreted and transferred to the FE solver. Each layer can in turn be decomposed into several layer fractions, according to the effective scan path and additional user-defined parameters. For instance, such a parameter could be a critical value of the scan length which would serve as a maximum possible value for the scan length of each layer fraction. Inside a single layer, the different layer fractions can be juxtaposed, or separated, as illustrated in Fig. 2.

Some characteristic features of each layer fraction are straightforwardly provided by the G-code interpretation: the scan length inside each layer fraction, denoted , and the time interval associated with the motion of the beam to the next layer fraction, . Knowing the laser beam scanning velocity, , the scan time of a layer fraction, , is easily deduced:

(1)

The total treatment time for a layer fraction can be defined considering an addition with the time of motion of the beam, However, this added time, , is generally small compared to , so that it will be ignored afterwards, the treatment time, , being restricted to the scan time, . It should be mentioned that concept of layer fraction is quite important to capture the local effect of the scanning strategy for both thermal and mechanical analysis compared to an entire layer deposition.

Two extreme situations can be described. First, taking an arbitrary high value for the critical scan length will result in deposition full layer by full layer. In other words, each layer has only one layer fraction with . In this case, less accurate but fast calculations can be run. By contrary, the model degenerates into the "equivalent source" hypothesis if the critical scan length is taken equal to the nominal beam diameter, . In this case, we would have for each layer fraction , making the calculation more accurate but less computationally efficient. Therefore the critical scan length is necessarily fixed between these two critical values, depending on the desired precision/efficiency ratio. The 1D scheme of Fig. 3 a) describes the concept of layer fraction through a decomposition of the entire layer trajectory scan.

Determination of energy input

Applying the thermal input along the laser pathas proposed in small scale numerical modellingwould result in poor computational efficiency due to the very small time steps required as well as a dynamic local refinement of the mesh in the region exposed to the laser beam. In the present context of macroscopic numerical modelling, a different and innovative strategy is proposed. Each layer fraction is assumed to be entirely heated during a time interval equal to the effective time during which the powder bed is exposed to the laser beam. This time interval is estimated considering the time required to cover the length corresponding to the laser beam diameter:

(2)

Note that this time interval is not specific to a layer fraction, i.e., its value is the same, whatever the considered layer fractions. In the numerical simulation, the time step is thus considered for the heating interval. After being heated in the single time step , the layer fraction is supposed to cool, as well as the rest of the domain during the time difference between and :

(3)

Over this cooling interval, a small calculation time step is defined as , where is a number of time steps which is chosen a priori. Note that is generally much smaller than .

As a consequence, defines the minimum time step in the proposed strategy.

A uniform heat flux is developed during the heating time along the surface of the layer fraction.

Denoting ̇ the surface heat flux that should be applied to the layer fraction, of area , during the heating time , its expression can be easily calculated by the following energy balance:

( ) ̇ (4)
where P L is the nominal laser power, and R the reflection coefficient of the laser radiation at the surface of the powder bed. Considering the numerical implementation and the need to heat a specific volume affected by the laser, this heat input is uniformly applied to the whole thickness of the powder bed, using a volume heat source ̇ . This uniform volume source is then given by

̇ ̇ ( ) (5)
The heating/cooling simulation strategy for each layer fraction is summarized in Fig. 3 b). Considering an overlay ratio between contiguous scan tracks as shown in Fig. 3 c), the layer fraction area can be related to the laser scan length by the following approached expression, which is purely geometric:

( ) ( ) (6)
Injecting Eq. ( 2) and Eq. ( 6) in Eq. ( 5) yields a new expression for the volume heat source:

̇ ( ) ( ) (7)
At the end of the simulation of an entire layer treatment (in a series of layer fractions), an inter-layer dwell time is required to spread a new powder layer at the top surface of the construction. This time interval is known from the G-code file and is considered as a cooling stage in the modelling process. It is accordingly simulated before going on with the heating stage of the first fraction of the newly deposited powder layer.

Dynamic mesh control

To make the simulation accurate, efficient, and computationally sustainable, the local mesh size and the total number of elements must be controlled. The objective is here to refine the FE mesh in the deposited layer, and progressively make it coarser with the distance to the construction front, without affecting the shape of the constructed workpiece. In order to achieve this, and as this may require the definition of anisotropic elements, an aimed mesh metric field is continuously calculated during the simulation and communicated to the remeshing procedure in charge of updating the nodal positions and the mesh topology. At each mesh node, the calculation of the associated 3×3 matrix, hereafter denoted , containing the aimed mesh size in the different spatial directions, is essentially based on the two level set functions and defined above, as explained hereafter. First, a global isotropic background mesh metric is defined. Then a fine anisotropic mesh metric is defined in the vicinity of the construction front. Finally, the metric is obtained by the intersection of and .

Background mesh metric on the whole domain

Firstly, an isotropic mesh metric is defined, with a size increasing with the absolute distance to the upper limit zone , defined as (see Fig. 1). is diagonal and can be expressed by the simple relation:

( ) with | | (8)
where is the prescribed mesh size in the gas domain along the construction front ( ), and is a mesh coefficient.

A second isotropic mesh metric is then defined with a mesh size depending essentially on the level set function . Like , is a diagonal matrix, ( ⁄ ), with the mesh size defined by Eq. ( 9) and explained below.

{ ( ) | | ( ( ) ) | | | | (9)
As described in the thesis of Shakoor [START_REF] Shakoor | Three-dimensional numerical modeling of ductile fracture mechanisms at the microscale[END_REF], Eq. ( 9) defines a linear transition from a curvaturedependent mesh size along the interface , to a mesh size equal to at a distance from this interface. denotes the minimum curvature radiusthat is to say the smallest radius among the two principal curvature radiiassociated to the local iso-distance surface (i.e. corresponding to the value) of the interface. It is defined as the inverse of the maximum eigenvalue (in absolute value) of the Hessian matrix of the level set function , and is a scaling parameter. A precise P1 description of the Hessian matrix (linked to the second spatial derivatives of ) is obtained thanks to a Superconvergent Patch Recovery strategy detailed in [START_REF] Shakoor | A new body-fitted immersed volume method for the modeling of ductile fracture at the microscale: Analysis of void clusters and stress state effects on coalescence[END_REF]. In addition, is a lower bound limiting the mesh size when the geometry is singular (near sharp corners or edges of the constructed workpiece). Similarly, the upper bound limits the mesh size near planar regions of the constructed workpiece (where ). Outside this thin layer of thickness , a mesh size is defined.

Based on the node-wise isotropic metrics and , the final background mesh metric is defined for each node of by the following mixture law:

( ) ( ( )) (10)
where the stepwise Heaviside function is defined by:

( ) { (11)

Mesh metric in , in the vicinity of the construction front

To capture accurately the heat input in the part, a fine mesh zone is selected in the region defined by , where (see Fig. 1). Contrary to the previous strategy, an anisotropic mesh is defined in order to limit the number of mesh elements. The mesh metric is defined as a function of the level set function , and of another function ( ), which approaches the signed distance to the yellow rectangle of thickness in Fig. 1. In the considered region, an anisotropic metric is first defined by

[ ( ) ( ) ] (12)
where , and are the aimed mesh sizes along the horizontal and the vertical direction, respectively. These quantities are defined by the following expressions:

{ | | { ( ) | | ( ( ) ) | | | | (13)
Thus, is defined as a function of , through a power law, where is a scaling coefficient ( ). This coefficient is defined in order to increase rapidly the mesh size out of the heating zone of thickness . The mesh size is defined as a function of the level set function , similarly to Eq. ( 9), using a scaling parameter . The aim of Eq. ( 13) is to define a fine mesh in compared with the background mesh of the entire domain . Applying the constraint , the truncated metric is obtained for the fine mesh by

( ( )) ( ) (14)
As and can be diagonalized in the same orthonormal basis (main axes x, y, and z)

considering their local expression, a metric intersection is easily obtained by considering the maximum eigenvalue in each of these directions:

| = ( | | )| (15)
Then the final mesh metric can be defined by the corresponding maximum eigenvalues of metrics and in each direction, which actually takes the minimum mesh size in the intersection domain.

[

( ) ( ) ] [ ( ) ( ) ] (16)
The final mesh metric is transferred to the automatic meshing procedure named Fitz [START_REF] Shakoor | A new body-fitted immersed volume method for the modeling of ductile fracture at the microscale: Analysis of void clusters and stress state effects on coalescence[END_REF], in order to get a new mesh respecting the desired metric. In addition, Fitz allows respecting interfaces within the analysis domain , when identified by level set functions: the new mesh is conform along the interfaces while respecting the objective metric. Details on the remeshing strategy with Fitz are available in references [START_REF] Shakoor | Three-dimensional numerical modeling of ductile fracture mechanisms at the microscale[END_REF][START_REF] Shakoor | A new body-fitted immersed volume method for the modeling of ductile fracture at the microscale: Analysis of void clusters and stress state effects on coalescence[END_REF][START_REF] Shakoor | An adaptive level-set method with enhanced volume conservation for simulations in multiphase domains[END_REF].

As a result, the elements are kept in a refined state in the vicinity of the construction front (region ) and are coarsened in other zones. More precisely, an isotropic and relatively coarse mesh is obtained for the constructed workpiece while preserving a conform mesh and avoiding numerical smoothing or irregularity of this interface due to remeshing operations. This adaptive meshing/remeshing strategy significantly reduces the number of elements compared to a standard level set based front-capturing remeshing strategy [START_REF] Chen | Three-dimensional finite element thermomechanical modeling of additive manufacturing by selective laser melting for ceramic materials[END_REF].

Thermal solver

Local heat transfer equations

The standard heat transfer equation is used:

( ) ̇ (17)
where T denotes the temperature, the density, the specific heat, the thermal conductivity. Note that in the context of a macroscopic simulation, and as melting and solidification are localized in a very small region, it is not necessary to take into account the latent heat of fusion/solidification. The right hand side of Eq. ( 17) consists of the volume heat source defined in Eq. ( 5). Eq. ( 17) is subject to an initial condition, which consists of a uniform temperature in the whole domain at time s.

Different boundary conditions are defined along the boundary of the domain considering the local environment. Along the lateral boundary of the material domain, , the heat extraction is expressed by

( ) ( 18
)
where n is the unit outward normal vector, and is a coefficient characterizing the conductive heat exchange with the machine environment, which is supposed to be at a temperature . Along the bottom boundary of the model, , the heat extraction is expressed by

( ) (19)
where the coefficient characterizes the convective heat exchange coefficient with surrounding air at temperature . Finally, an adiabatic condition (null heat flux) is assumed along the boundaries of the gas domain . In the present model, the heat exchanges along interfaces and are purely diffusive, no convection is taken into account.

Finite element formulation

The Galerkin FE method applied to the problem defined by Eqs. ( 17) -( 19) yields the first-order differential equation:

(

) 20
where T is the vector of nodal unknown temperature values, C is the heat capacity matrix, K is the thermal diffusion matrix and F is the internal and external flux vector. The standard linear tetrahedral element is used for spatial discretization.

Thanks to the remeshing module Fitz, the mesh is maintained conform along interfaces. Accordingly, in the present work, the treatment of local and possibly abrupt changes of material properties is based on a P0 interpolation procedure, as detailed hereafter. The procedure begins by defining the value of level set functions and at the barycenter of each finite element e:

( ) ∑ ( ) ∑ (21)
where (respectively ) denotes the value of ( ) at node n, while ( ) (respectively ( ) ) denotes its averaged value on the tetrahedron. This being done, it is possible to define two Heaviside functions of P0 type which take the sole values 0 or 1. A first Heaviside function is attached to the gas domain, and a second one, , to the workpiece. They are defined by

( ) ( ( ) ) { ( ) ( ) ( ( ) ) { ( ) (22)
Then, in each element e, a P0 type estimation of any material property P (density, conductivity, or specific heat) can be expressed by the following expression, in which index (e) is omitted for the sake of clarity:

= ( ) ( ( ))[ ( ) ( ( )) ] (23)
where , and denote the value of variable P in gas, in the consolidated workpiece and in the non-melted powder, respectively. This methodology enables to avoid standard mixture laws for the material properties which are usual for front-capturing approaches. This is an advantage as they are always difficult to calibrate and generally require a significant mesh refinement around interfaces.

Values obtained by Eq. ( 23) are then used in each element e to form the element contributions to the different matrices and vectors of Eq. ( 20), according to the standard FE method. Detailed expressions are not given here; they can be found in any FEM reference book.

Methods to control thermal shock effects

It is well known that in heat transfer non steady state FEM, both spatial and temporal instabilities are encountered when the application of a boundary condition generates a steep gradient that cannot be captured by a too coarse mesh [START_REF] Putti | Finite element approximation of the diffusion operator on tetrahedra[END_REF]. In this work, two numerical methods are used in order to avoid such so-called thermal shock effects. They will be compared.

To avoid diffusion-type instabilities, the mesh size should satisfy the following condition [START_REF] Fachinotti | Linear tetrahedral finite elements for thermal shock problems[END_REF]:

√ ( 24
)
where is the chosen time step. In SLM simulations, two reasons may cause the occurrence of such shock. First, the application of the laser heat input to a region which possibly has a temperature significantly lower than the melting temperature, and second, the initialization of the temperature of the newly-spread powder layer. Therefore it is of utmost importance to guarantee that the criterion expressed in Eq. ( 24) is satisfied in region , which is the locus for heat and matter input. It can be seen that the small time step associated with heating of layer fractions, which is given by Eq. ( 2), plays a critical role.

The simplest method consists in decreasing the mesh size under the critical value defined above.

However, for a large scale model, this might be not sustainable regarding computation time, as too many elements should be introduced, especially in the powder zone because of its low thermal conductivity. To illustrate the point, in the present study the critical mesh size is about 0.027 mm. To minimize the computational effort, two methods have been investigated; they are briefly described in the two sub-sections hereafter.

 Scaling coefficient for the definition of the heating time

The first method consists in increasing arbitrarily the time interval (and accordingly, the time step) associated with the heating of a layer fraction, which is given by Eq. ( 2). For that, a scaling coefficient is introduced, expressing that this heating time is now times the time interval during which the powder bed is effectively exposed to the laser beam, yielding:

(25)

Proceeding this way, the critical mesh size is increased by a factor √ in , which decreases the number of elements in the model.

 Asynchronous method

An alternative or complementary method is the so-called asynchronous method initially proposed by Jaouen [START_REF] Jaouen | Modélisation tridimensionnelle par éléments finis pour l'analyse thermo-mécanique du refroidissement des pièces coulées[END_REF] in the context of modelling solidification processes. It consists in considering a larger time step , significantly higher than the limiting heating time step . At a given time t, the transient heat transfer problem is then solved using the time step , yielding a provisional solution . The incremental temperature variation is then linearly interpolated to express the solution at time :

( ) (26)
It should be noted that the asynchronous method is quite stable, as the thermal solver is fully implicit, and so, unconditionally stable whatever the amplitude of the time step. Using this method, the mesh size in region has now to respect Eq. ( 24), but using as the time step. Hence, the critical mesh size is increased by a factor √ , which may dramatically decrease the number of elements and the computational effort. The impact of the two methods presented, and their associated parameters ( , ), will be discussed further, in terms of accuracy and computational time.

Modelling of the construction process

After the calculation has been run on the time interval corresponding to the treatment of a whole layer, including the inter-layer dwell time, it is necessary to continue the simulation of the construction. This requires updating the mesh, transporting the level set function to the new mesh, modifying the value of the level set function on the mesh and defining the initial temperature in the newly created powder layer. Consequently, the level set function ( ) at the new time step is updated in the vertical (z) direction with respect to its value at the previous time step t:

(27)
In the newly deposited powder layer, i.e. for nodes such that , the temperature is reset to the value , as done in [START_REF] Michaleris | Modeling metal deposition in heat transfer analyses of additive manufacturing processes[END_REF], while the temperature of the previous time step is conserved in the material domain . In practice, the nodal temperature values are re-initialized using the previous value of the Heaviside function and :

( ) ( ( )) (28)

Implementation

The procedures described in the previous subsections were implemented in the C++ library CimLib developed in CEMEF. The execution flow chart is shown in Fig. 4. During the preparation stage, the CAD model of the analyzed geometry is created, and the G-code of the construction process (laser trajectory, dwell time…) is generated by the "Slicer". Thanks to the "Interpreter", relevant information for the numerical simulation is obtained and the determination of the fractions of each layer can be done. The developed numerical model is then launched to perform the simulation of the whole construction process, by chaining the treatment of each layer, fractions of each layer, inter-layer dwell time, deposition of a new layer and remeshing.

Model validation: energy conservation

The possible causes for a non-conservative balance for energy in the calculation can be listed as follows: space and time discretization; treatment of the interfacial zones in the model (mixing laws); discretization of the applied heat source Eq. ( 7); initiation of a newly deposited powder layer Eq. (28).

To validate the developed model regarding energy conservation, it is proposed to study the construction of a brick shape workpiece as given in Fig. 5 a).

In this demonstration, deposition of material and energy by entire layer is adopted to simplify the deposition process (i.e. division of layers into several fractions is not considered). Material properties, heat exchange and laser parameters are given in Table 1. It can be seen from density information that the initial porosity of the powder bed is taken as 61%. Heat exchange conditions are taken adiabatic for a straightforward validation of energy conservation.

In line with Fig. 3 a), the heating of each layer is calculated with a single time step. The cooling stage is also calculated with a unique time step ( ). A new layer of material is added at the end of the cooling stage during the dwell time. This deposition scheme is repeated all along the construction process. Associated characteristic times and lengths are given in Table 2.

Sensitivity to the mesh

The meshing/remeshing parameters used in this simulation are given in Table 3. In the zone with fine mesh, near the heating zone, it can be seen by comparison with Table 2 that they should yield meshes respecting the criterion on mesh size, Eq. ( 24). This will guarantee an accurate description of the input energy.

Fig. 6 shows the mesh evolution at interfaces and in the global domain. It can be seen that the mesh is maintained fine in the neighborhood of the heating zone. Once the construction front has passed, the fine mesh is replaced by a coarser mesh. The conform mesh strategy allows respecting precisely the constructed workpiece interface during de-refinement, while the fine mesh persists in the vicinity of the construction front. However, the mesh is relatively coarse at the interface between gas and material far from the heating zone, as shown in Fig. 6 b) and c). This could influence energy conservation when a new powder layer is deposited. This is why it is required to check energy conservation, and especially the influence of the factor which controls the mesh size variation near the construction front, in . Therefore a parametric sensitivity study is performed, the results being given in Fig. 7.

As it can be seen in Fig. 7 a), the element density decreases when increases. Regarding the energy of the system, note that the theoretical evaluation of its variation with time is quite easy, as the duration of laser exposition is known, and boundary conditions are taken as adiabatic. It is found that the relative difference between simulation and theory is in the range -3% to +4% (Fig. 7 b)). The stepwise shape of the energy evolution is due to the long interlayer dwell time (10 s), during which energy is unchanged. The simulated evolution superimposes the theoretical one for . From a practical point of view and regarding application to real process and materials, all the three results are acceptable, the relative error being small compared to other uncertainties on the values of material and process parameters. For the case with , the temperature distributions for the entire model, the workpiece with powder, and the sole workpiece are given in Fig. 8.

In addition, the evolution of the temperature profile along the vertical direction V1 (as defined in Fig. 5 b)) during the construction process is plotted in Fig. 9. Each curve is plotted after cooling and dwell time following the deposition of a certain number of layers. The first part of each curve is within the constructed workpiece, the second part lies in the gas domain. It can be seen that, under the defined process conditions, the vertical temperature profile in the workpiece tends to converge with time.

Study of the influence of

Keeping the same remeshing parameters, the method consisting of a time scaling of the heating time is now tested. The value of the scaling coefficient for heating time, (Eq. ( 25)), is changed from 1 (reference without time scaling) to 4 and 8. The corresponding heating and cooling times as well as the values of the applied volume heat source are given in Table 4.

To study the temperature evolution, two capture lines are selected in vertical and horizontal directions as described in Fig. 5 b). At the end of heating of the 40 th layer, looking at the central vertical profile, it can be seen that the peak of temperature is all the more marked as is small (Fig. 10 a)). A maximum temperature of 2697°C is obtained in the material at the gas interface for = 1. Using higher values of this scaling coefficient, this maximum temperature is reduced (1410°C for = 8) at the interface. At the end of the cooling stage just following deposition of layer 40, the difference between curves is greatly reduced, and after the dwell time, there is no more difference between the three temperature profiles, which means that the scaling coefficient method does not influence the simulation results. Similar remarks can be done for the horizontal temperature profile at mid-height of the built workpiece (Fig. 10 b)). As the studied zone is far from the heating zone, less difference is found between the three cases, even at the end of heating and cooling steps. After cooling down during the inter-layer dwell time, profiles are found perfectly superimposed for the three cases. In total, this confirms the relevance of the -method (time scaling applied to heating).

Coarse mesh with numerical methods for avoiding thermal shock

As discussed in Section 2.4.3 both the -method and the asynchronous method permit a coarser mesh definition. Therefore, the metric in the sensitive zone is now defined with parameters allowing for a coarser mesh than in the previous reference case. Two simulations will be done with these new meshing rules, one with the -method using the value , and one using the asynchronous method. The associated characteristic data are given in Table 5. The mesh parameters are given in Table 6; they permit a reduction of the number of elements to 0.6 million.

As shown in Fig. 11 a), the curves of and asynchronous methods are found superimposed for the vertical profiles at the end of heating. The same tendency keeps for cooling after the dwell time with a difference of 2.02% between the maximum values of the two proposed methods compared with the reference case. For the horizontal profiles, the two proposed methods still lead to similar results for both heating and cooling after dwell time, where a 2.8% difference is found for the two methods compared with the reference case for the maximum temperature.

In addition, the relative difference for the calculated energy is less than 3% compared with the reference value (Table 7), while the computational time is reduced by more than 2.5 for both methods.

As a conclusion, the -method and the asynchronous method can be recommended as they offer a good compromise between accuracy and computational cost.

Application

Geometrical parameters

Using the strategy presented above, the SLMadditive manufacturing of a part is simulated. As shown in Fig. 12, the size of the simulated system is 30×30×16 mm 3 , with a 2 mm thick substrate. The part to be constructed is shown immersed in the global mesh. It is an impeller with 4 spiral blades, the thickness of which is 1 mm. The radius of the lower plate and the height of the workpiece are 10 mm.

Moreover, a vertical central hole goes through the workpiece. The height of the argon gas domain is supposed to decrease continuously during construction to reach a final minimum value of 4 mm. The deposition of energy is considered by layer fractions, which are based on the laser trajectories from a G-code description as detailed previously.

Thermal modelling

The temperature history is calculated by the thermal solver presented in section 2.4. For the thermal exchange, both convection and radiation are considered for bottom and lateral surfaces, the latter evolving with the level set function during the construction process. The associated values can be found in Table 8. The part material is the nickel-based alloy IN718. The thermal properties of the different materials are considered as temperature independent and are the same as that of Table 1.

Mesh adaptation

Considering the heating time , which is equal to 3.6 ms, and using the criterion of Eq. ( 24), the critical mesh size is 0.017 mm, which is too small to mesh the present model. Therefore, themethod is adopted to increase critical lengths in the three domains. With = 40, the characteristic times and lengths, and the mesh parameters for constructing mesh metrics are given in Table 9 and Table 10, respectively.

As shown in Fig. 13 a), the mesh in the gas region far from the construction front is coarse, while the zone close to the construction front is kept fine enough to avoid thermal shock. For the constructed workpiece, Fig. 13 b), fine elements are present along edges to preserve the geometrical shape, while coarser ones can be found inside. Using remeshing after deposition of each layer, the number of elements in the system can be stabilized around 2 million, and the total computation time is 17 hours on a 60 processor cluster.

Using the above strategy, and defining specifically adiabatic boundary conditions along the whole boundary of the system, energy conservation can be tested for the prescribed model. A 3% error is obtained by comparison with the theoretical input value, which confirms that proposed methods are acceptable. This verification being made, the simulation is performed using the heat exchange conditions detailed in Table 8.

Simulation results

After deposition of 50 layers and cooling down with a final dwell time of 15 s, the calculated temperature field in the entire model is given in Fig. 14 a). By removing the upper gas zone, the material zone is presented in Fig. 14 b). This allows observing the temperature distribution in the powder bed and in the consolidated workpiece along the top and lateral surfaces. The highest temperature is around 600 °C, in and near the last heated zones, while the lowest temperature is around 50 °C, along edges of the lateral surface. . By removing the non-exposed powder, Fig. 14 c) reveals a vertical temperature gradient in the part, the temperature in the bottom plate and the substrate being around 400 °C.

To better evaluate the temperature distribution, vertical and horizontal profiles in the hot zone (V1: The capture lines in the cool zone go through one blade, while the ones in the hot zone cross through the center. The vertical line V1 crosses the substrate from 0 to 2 mm, the powder bed from 2 to 12 mm and argon gas from 12 to 16 mm in z direction, while the horizontal line H1 crosses several times through the powder and the workpiece when the central hole is encountered. The vertical line V2

crosses the substrate from 0 to 2 mm, the workpiece from 2 to 3 mm, the powder, one spiral blade, powder and gas, while H2 crosses the powder, one spiral blade and powder. All transitions between the different materials are marked with small red circles on the curves plotted in Fig. 15 c) and d).

As shown in Fig. 15 c) the temperature tends to be homogeneous in the condensed material (substrate and workpiece) for both V1 and V2. The temperature in V1 increases in the powder zone along z direction and reaches the maximum value of 600 °C near the interface between powder and gas, then decreases in the gas zone. As V2 is near the lateral surfaces, the temperature starts to decrease after a peak at 483 °C at mid-height of the formed workpiece. The temperature increases slightly in the gas, which is maybe affected by the hot gas in the center (red zone in Fig. 15 a)). The temperature profiles in the horizontal direction are presented in Fig. 15 d). Along H1, a maximum value of 580 °C is found in powder zone, in the vertical center hole. For the capture line H2, the maximum temperature locates at the interface of the powder and the spiral blade. Generally speaking, the temperature gradients in the powder zone are more marked than in the workpiece zone, because of the lower diffusivity of the granular material. For profiles (H1, V1) located in and crossing the center hole, the maximum temperature always locates in the powder zone. For the profiles not crossing the center (H2, V2), the maximum temperature is found at the interface between workpiece and powder. In total these results

show the relevance of the proposed approach in which heat transfer is resolved throughout the whole system, including powder. Contrary to simpler methods which consist in applying estimated heat exchange boundary conditions directly along the surface of the workpiece during its construction, the proposed method seems to offer extended perspectives regarding an in-depth understanding of heat transfer during SLM.

Conclusion

A 3D macroscopic finite element thermal model for SLM process has been developed on the scale of the constructed workpiece. The simulation encompasses the non-melted powder bed and the surrounding gas. A layer fraction strategy based on a description of the global laser scan plan has been implemented. Level set functions are used to track interfaces between workpiece and powder, and between material and gas. An adaptive remeshing technique, capable of local refinement and derefinement of the mesh, while preserving a conform mesh along interfaces, is used to control mesh definition in the whole model during the construction process. In addition, regarding heat transfer resolution, two methods have been tested and validated to reduce computational time: time scaling applied to heating of layer fractions, and the asynchronous resolution. It is demonstrated that both techniques offer a good compromise between accuracy and computation time.

The developed solution strategy has been applied to the simulation of the SLM of a part of a certain complexity. Results show that temperature gradients in the part zone and in the substrate tend to be lower than in the powder zone, due to different heat conductivities. The simulation results bring to the fore the influence of the powder zone on the global temperature distribution: higher gradients, hot spots due to heat accumulation.

The developed 3D macroscopic modelling of heat transfer in SLM of metal parts can be applied without a priori limitation regarding the geometrical complexity. However, the current computation times are still high. The optimization of the definition of mesh metrics in the conform remeshing strategy is certainly a way to decrease the computational effort. Future work will also focus on one hand on quantitative comparisons with experimental measurements, which are in progress, and on another hand on the implementation of a mechanical resolution module in order to predict the stress build-up in the part during its construction, and the residual stresses.

Fig. 1

 1 Fig.1presents a schematic of the FEM strategy. Firstly, the geometry of the workpiece,, is given by the information extracted from a computer-aided design (CAD) model. It includes the final part plus the supporting structures, which constitute the full consolidated geometry to be produced by SLM. It is immersed in a larger analysis domain , also embedding a bottom substrate, , and a gas domain, , to fully surround the CAD model.

[ 15 ,

 15 [START_REF] Riedlbauer | Macroscopic simulation and experimental measurement of melt pool characteristics in selective electron beam melting of Ti-6Al-4V[END_REF], 0] ⎼[START_REF] Riedlbauer | Macroscopic simulation and experimental measurement of melt pool characteristics in selective electron beam melting of Ti-6Al-4V[END_REF][START_REF] Riedlbauer | Macroscopic simulation and experimental measurement of melt pool characteristics in selective electron beam melting of Ti-6Al-4V[END_REF][START_REF] Steuben | Discrete element modeling of particle-based additive manufacturing processes[END_REF] mm, H1: [15, 0, 8] ⎼[START_REF] Riedlbauer | Macroscopic simulation and experimental measurement of melt pool characteristics in selective electron beam melting of Ti-6Al-4V[END_REF] 30,[START_REF] Megahed | Metal additive-manufacturing process and residual stress modeling[END_REF] mm) and cool zone (V2:[START_REF] Denlinger | Thermal modelling of Inconel 718 processed with powder bed fusion and experimental validation using in situ measurements[END_REF][START_REF] Chen | Three-dimensional finite element thermomechanical modeling of additive manufacturing by selective laser melting for ceramic materials[END_REF] 0] mm ⎼[START_REF] Denlinger | Thermal modelling of Inconel 718 processed with powder bed fusion and experimental validation using in situ measurements[END_REF][START_REF] Chen | Three-dimensional finite element thermomechanical modeling of additive manufacturing by selective laser melting for ceramic materials[END_REF][START_REF] Steuben | Discrete element modeling of particle-based additive manufacturing processes[END_REF] mm, H2: [7, 0, 5] mm ⎼ [7, 30, 5] mm) are selected in the entire model as shown in Fig.15 a).

(Fig. 1 .

 1 Fig. 1. Schematics of the domains defining the system for FEM analysis.

Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 Fig. 8 .

 2345678 Fig. 2. Decomposition of a single powder layer into several fractions to be heated consecutively.

Fig. 9 .Fig. 10 .Fig. 11 .Fig. 12 .Fig. 13 .Fig. 14 .Fig. 15 .

 9101112131415 Fig. 9. Evolution of temperature profile along V1 (see Fig. 5 b)) during construction.

Table 1

 1 Material properties and process parameters. The relative density of the powder bed (here 0.4) is slightly underestimated comparatively to usual values in SLM. However, this should have a limited impact on the calculated temperature distributions.

			Properties		Workpiece	Powder bed Argon gas
	Materials [24]	Density, ρ [kg m -3 ] Thermal conductivity, k [W m -1 K -1 ]	8185 10.8	3191.5* 0.13	1.3 0.024
		Specific heat, c p [J kg -1 K -1 ]		427	427	1000
		Machine environment temp.,	[°C]	20
		Initial temperature, [°C]			20
	Heat exchange	Ambient temperature, Conductive coefficient,	[°C] [W•m -2 •K -1 ]	20 0
		Convection coefficient,		[W•m -2 •K -1 ]	0
		Nominal power, [W]			115
		Scan speed, ν L [mm s -1 ]			55.03
		Reflection coefficient, R			0.55
	Laser and	Beam diameter,	[mm]			0.5
	process [7]	Thickness of powder layer,		[mm]	0.15
		Overlay coefficient,			0
		Interlayer dwell time [s]			10
		Scaling factor				1
	*

Table 2

 2 Characteristic times and lengths for the reference case.

				Critical mesh size
	Characteristic times [ms]			(Eq. (24)) [mm]	̇
							[W•	]
	Heating	Cooling	Inter layer	Powder	Workpiece	Gas
	9.1	136.3	10000	0.027	0.168	0.410	1380

Table 3

 3 Mesh parameters for the reference case.

	[mm]	Number of
		elements
		[million]

Table 4

 4 Characteristic times and lengths for the time scaling method.

		Characteristic times [ms]	̇
	Value of	Heating Cooling	Inter layer	[W mm -3 ]
	1 (ref.)	9.1	136.3	10000	1380
	4	36.3	109.0	10000	345
	8	72.7	72.7	10000	172.5

Table 5

 5 Characteristic times and lengths for methods to avoid thermal shock.

					Critical mesh size
		Characteristic times [ms]		(Eq. (24)) [mm]	̇
								[W mm -3 ]
	Test cases	Heating	Cooling Inter-layer Powder Workpiece	Gas
	Reference ( = 1)	9.1	136.3	10000	0.027	0.168	0.410	1380
	-method ( = 4)	36.3	109.0	10000	0.053	0.335	0.819	345
	Asynchr. ( = 1)	9.1	136.3	10000	-	-	-	1380

Table 6

 6 Mesh parameters for methods to avoid thermal shock.

				Mesh parameters [mm]
		(i = 2, Eq. (9))	0.05	0.3	1.0	0.6 0.005	-	10	-	0.5
		Reference ( = 1)	0.025 0.15 0.15 0.15 0.005 1.75	-	10	-
	(i = 3,
	Eq.	-method ( = 4) 0.05	0.3	0.3	0.3 0.005 2.0	-	10	-
	(13))
		Asynchr. ( = 1)	0.05	0.3	0.3	0.3 0.005 2.0	-	10	-

Table 7

 7 Simulation results.

	Test cases		Number of elements [million]	Calculated energy Value [J] Difference [%] Value [min] Accel. ratio Simulation time (20 cores)
	Reference (	)	0.97	300.93	-	675	-
	-method ( = 4)	0.61	308.41	+2.49	252	2.68
	Asynchr. (	)	0.61	308.28	+2.44	227	2.97

Table 8

 8 Heat exchange and process conditions for the application case.

	Properties	Values

Table 9

 9 Characteristic times and lengths for the application case.

		Characteristic times [ms]	Critical mesh size [mm]	̇
		Heating	Cooling	Inter-layer	Powder Workpiece	Gas	[W mm -3 ]
	= 1	3.6	varying	15000	0.017	0.105	0.250	6468.4
	= 40	145.4	varying	15000	0.106	0.670	1.638	161.71

Table 10

 10 Mesh parameters for the application case.

	[mm]	Number
		of
		elements
		[million]