Additive Manufacturing (AM) or 3D Printing is emerging as a prominent manufacturing technology  in many industrial sectors, such as the aerospace, defence, dental or biomedical. These sectors are on the lookout to find the way to exploit its many potentials, including vast geometrical freedom in design, access to new materials with enhanced properties or reduced time-to-market. However, this growth cannot be long-term sustained without the support from predictive computer simulation tools. Only them provide the appropriate means to jump the hurdle of slow and costly trial-and-error physical experimentation, in product design and qualification, and to improve the understanding of the process-structure-property-performance link.
This work concerns the numerical simulation at the part scale of metal AM processes by Powder-Bed Fusion (PBF) technologies, such as Direct Metal Laser Sintering (DMLS), described in [8, Fig. 1], Selective Laser Melting (SLM), or Electron Beam Melting (EBM). Compared to other metal technologies, powder-bed methods feature the thinnest layer thicknesses, from 60 to 20 microns (or even below), and the fastest solidification rates. On the one hand, building industrial components usually requires depositing thousands of layers; thus, from the modelling viewpoint, computational efficiency should be at the forefront. On the other hand, high heat fluxes concentrate in the last printed layers; away from that region, the thermal distribution is much smoother and less physically relevant.
Many researchers have used the Finite Element (FE) method to investigate metal AM processes, often aided by their knowledge of modelling other well-known processes, such as casting or welding [4, 17]. At the part scale, early FE models proved their applicability to many engineering problems, e.g. selection of process parameters , design of scanning path  or evaluation of distortions and residual stresses [10, 26, 20]. However, numerical tests were often limited to short single-part builds and small deposition volumes [28, 16].
More recently, the attention has turned to the design of strategies to accelerate simulations to tackle longer processes, multiple-part builds, higher deposition volumes and, eventually, deal with industrial-scale scenarios in reasonable simulation times. Some authors have attempted to exploit adaptive mesh refinement [12, 25, 15] and/or parallelization , while most of them consider surrogate models that are inevitably accompanied by some sacrifice of accuracy or physical representativeness.
Typical model simplifications , alone or combined, consist of time-averaging the history of the process, by lumping welds or layers [14, 19], or reduce the domain of analysis by, e.g. excluding the region of loose powder-bed surrounding the part and/or the base plate. In the case of thermal analyses, heat transfer from the part to the excluded region is then accounted for with an equivalent heat loss boundary condition at the solid-powder interface . However, determination of these boundary conditions has been limited to rather simple approximations and challenged by lack of experimental data, especially concerning heat conduction through the loose powder. In this sense, a path yet to be explored in metal AM is to develop physics-based alternatives, similar to those studied in the area of casting solidification to model the mould [11, 31, 13], to replace the heat flow model at regions of less physical interest with a much faster and less memory demanding model.
The purpose of this work is to establish a new physics-based rationale for domain reduction in the thermal FE analysis of metal AM processes by powder-bed methods. The new technique, referred to as the Virtual Domain Approximation (VDA) in Sect. 3, approximates the 3D transient heterogeneous heat flow problem (see Sect. 2) at low-relevance subregions (e.g. loose powder, plate) by a 1D heat conduction problem (see Fig. 1). This much simpler 1D problem can be analytically solved and reformulated as an equivalent boundary condition for the 3D reduced-domain problem.
Applying this formulation to a simple proof-of-concept example in Sect. 4.1, reduced Part-Plate (PP) and Part-only (P) models are derived and compared with respect to a full powder-bed-part-plate High Fidelity (HF) model. The results verify the ability of the spatially reduced models to approximate the reference response, but with a speed-up greater than ten. A second numerical test in Sect. 4.2 addresses experimental validation against physical tests carried out at the Monash Centre for Additive Manufacturing (MCAM) in Melbourne, Australia, using an EOSINT M280 machine and Ti-6Al-4V powder. In this case, a HF model is first calibrated and validated against the experimental data. This step is next repeated for two ancillary PP reduced models. In the first one, the Heat Transfer Coefficient (HTC) between the part and the powder is assumed constant, as in earlier works [8, 7]. In the second one, the PP adopts the VDA formulation.
Given the scale of the experiment, 17.5 of deposited volume in 647 layers and 3.5 of process, the HF model size amounts to 9.8 million unknowns, whereas the PP models to 0.7 million. Sensitivity analysis and parameter estimation at such level of complexity are only practical by means of an advanced computational framework. For this reason, the numerical tests are supported by a high-end parallel computing framework, unprecedented in the simulation of metal AM processes. The framework combines three tools: (1) FEMPAR-AM , a module of FEMPAR , a general-purpose object-oriented message-passing/multi-threaded scientific software for the fast solution of multiphysics problems governed by Partial Differential Equations; (2) Dakota , for the automatic parametric exploration of the models; and (3) TITANI 7, a High Performance Computing (HPC) cluster at the Universitat Politècnica de Catalunya (Barcelona, Spain) to support the calculations. Using this innovative methodology for the MCAM experiment, the physics-based VDA-PP model is shown to reproduce the response of the HF model with increased accuracy, with respect to the constant PP variant. However, the simulation times of both reduced models is practically the same. In other words, the extra cost devoted to evaluate the equivalent boundary condition with the VDA formulation is negligible, in front of the computational cost devoted for the solution of the problem at the reduced domain.
As a result, the VDA
formulation offers good compromise between accuracy and efficiency. Indeed, on the one hand, the computational benefit is clear, as the mesh covers a smaller region and the number of degrees of freedom with respect to the complete model is significantly reduced. On the other hand, the impact on accuracy is efficiently controlled, because the neglected physics are taken into account in the evaluation of the equivalent heat loss boundary condition, without affecting the computational cost of the simulation. Even though theVDA formulation does not totally get rid of cumbersome parameter estimations, it can be easily calibrated with respect to a HF model. This makes it an appealing alternative for optimization problems (e.g. design of minimum-distortion scanning path), or any other engineering problem involving a parametric exploration. More interestingly, by exploiting the locality of the PBF process, the VDA could eventually be useful to reduce the domain down to a few layers, as shown in Fig. 0(d), while keeping good relation with the thermal response of the HF model.
2. Heat transfer analysis in AM
This section describes the transient model for the thermal FE analysis of the printing process at the part scale. The contents centre upon (1) an overview of the model variants studied in this work, which depend on the domain of analysis and heat loss model (Sect. 2.1); (2) the governing equation and the determination of the material properties at the powder state (Sect.2.2); and (3) the treatment of boundary conditions, according to the heat loss model (Sect. 2.3). The reader is referred to [7, 8] for further details on the weak formulation and the FE modelling of the geometry growth during the metal deposition process.
2.1. Model variants
Let be an open bounded domain in , representing the system formed by the printed object , the building plate and the surrounding powder bed . grows in time during the build process. After the printing, it remains fixed, while cooling down to the room temperature.
Several variants of the thermal model, represented in Fig. 2, arise from considering different computational domains :
If , the model is referred to as High Fidelity (HF). The contour of is formed by a region in contact with the air in the chamber , the lateral wall and the bottom wall of the plate .
If the powder bed is excluded from the computational domain, i.e., , then the so-called PP model is obtained. In this case, the solid-powder interface is also in the boundary, i.e., .
Finally, taking only yields the P model, with .
From the computational viewpoint, HF is the most demanding model, because the part, the plate and the powder bed must be meshed, whereas P is the simplest, because it only needs to mesh the part. On the other hand, HF is the most physically representative model, among those studied here. Therefore, it is established as the reference model.
As a result of excluding the powder bed (and the building plate), heat conduction through the powder (and the plate) must be accounted for with an equivalent heat loss boundary condition at (or and ). For this purpose, two strategies are studied: (1) constant-valued HTC and (2) VDA (cf. Sect. 3). They lead to further submodels, such as HTC-PP, VDA-PP and VDA-P. Only these three are covered in this work. In particular, their computational cost and accuracy with respect to the reference HF is assessed in Sect. 4.
2.2. Governing equation
Heat transfer in at the part scale is governed by the balance of energy equation, expressed as
where is the heat capacity coefficient, given by the product of the density of the material and the specific heat , and is the thermal conductivity. Furthermore, is the rate of energy supplied to the system per unit volume by a very intense and concentrated laser or electron beam that moves in time according to a user-defined deposition sequence, referred to as the scanning path. It is computed as
with the laser power , the heat absorption coefficient, a measure of the laser efficiency, and is the volume swept by the laser during the time increment . Note that phase transformations occur much faster than the diffusion process and the amount of latent heat is much smaller than the energy input . That is why, given the scale of analysis, phase-change effects are neglected.
In case of the HF model, the powder-bed is included into the computational domain. As a result, there are two distinct material phases playing a role in (1), i.e. solid and powder, which are separated at the interface between the solid part-plate ensemble and the granular powder-bed. To model the material in powder state, the thermophysical properties are determined in terms of the bulk material data at solid state and the porosity of the powder-bed . In particular, density and specific heat are obtained as
while determination of thermal conductivity is more involved and often relies on empirical expressions. Among the models present in the literature, Sih and Barlow  establish, for a powder-bed composed of spherical particles , the relation
where is the thermal conductivity of the surrounding air or gas and is the contribution of radiation amongst the individual powder particles, given by Damköhler’s equation:
with the average diameter of the particles and the Stefan-Boltzmann constant, i.e. .
2.3. Boundary conditions
Eq. (1) is subject to the initial condition
with the pre-heating temperature of the build chamber, and the boundary conditions applied on the regions shown in Fig. 2, which depend on the model variant:
HF model: Heat convection and radiation through the free surface , heat conduction through the lateral wall on and heat conduction through the bottom wall of the plate on .
PP model: Heat convection and radiation only through , heat conduction through the powder bed along , heat conduction through the lateral wall only on and heat conduction through the bottom wall of the plate on .
P model: Heat convection and radiation through , heat conduction through the powder bed along and heat conduction through the building plate on .
After linearising the Stefan-Boltzmann law for heat radiation , all heat loss boundary conditions mentioned above can be expressed in terms of the Newton law of cooling:
where refers to the kind of heat loss mechanism and the boundary region where it applies. Alternatives to Eq. (8) can also be considered, e.g. the temperature on can be prescribed to , taking into account that the thermal inertia of the building plate is much larger than the printed part.
For the HTC model, and are taken as constant, due to the difficulty in describing experimentally the temperature dependency of these quantities. On the other hand, for the VDA case, and are both temperature dependent and computed as described in Sect. 3.
3. Virtual domain approximation
As mentioned in Sect. 2.1, if the powder-bed (and the building plate) are not included in the domain of analysis, then heat transfer through these regions must be accounted for with equivalent heat conduction boundary conditions. According to Eq. (8), this leads to the determination of and values for heat loss through the solid-powder interface and the part-plate interface.
However, lack of experimental data and physical modelling in the literature are a hurdle in the way to estimate these quantities. The authors are not aware of the existence of, for instance, well-established temperature-dependent empirical correlations for the HTC solid-powder or approximate time evolution laws for the temperature at the building plate.
To overcome these challenges, this work introduces a novel method, i.e. the Virtual Domain Approximation (VDA). The VDA enhances a state-of-the-art thermal contact model for metal casting analysis , such that it includes the effect of the thermal inertia of the powder bed (or the building plate). The result of the VDA procedure is a temperature-dependent physics-based way to evaluate and in Eq. (8), a clear improvement with respect to the HTC model.
Assuming the PP variant, the VDA method is outlined in Fig. 3 and explained in the following paragraphs. It is based on the assumption that the entire 3D heat transfer problem across the powder bed can be modelled with an equivalent 1D heat conduction problem across a contouring wall on . This problem is, in fact, equivalent to the classic 1D transient heat transfer through a plane wall .
According to this, the setting considers a 1D spring with material properties at powder state, subject to thermal contact at both ends. At one end of the spring, the powder-bed surface is in contact with the solid part/plate surface. As the powder-bed is a porous medium, it is reasonable to assume that the contacting surfaces do not match perfectly. In particular, the standard Fourier law does not hold. Hence, the heat flux is computed as the product of an HTC and the thermal gap between both surfaces, i.e.
The spring represents the region of the powder-bed subject to relevant thermal effects of the printing process (i.e. with presence of significant thermal gradients). It has average length , referred to as the effective thermal thickness. A point in the powder bed that distances itself more than with respect to the part is barely affected by the thermal gradient and undergoes little changes in temperature during the printing process, i.e. it mostly remains at .
At the other end of the spring, thermal powder-to-powder contact applies. The expression for the heat flux is analogous to Eq. (9); in this case, the HTC is and the thermal gap can be taken as , in agreement with the discussion at the previous paragraph. Hence,
Assuming now a VDA attached to each integration point located on a face in , is determined from the solution of Eq. (1), whereas the material properties of the powder, i.e. , and , and the parameters , and are assumed to be known data. This leaves the VDA model as a simple well-posed 1D transient heat transfer problem. The first step of the VDA method is to apply a suitable discretization of this 1D problem. The resulting linear system is then modified, by adding as an additional unknown with the extra equation given by Eq. (9). Following this, static condensation allows one to recover , such that it no longer depends on , only on , the known data and the type of discretization. In this way, an equivalent expression of Eq. (9) is obtained, that takes the same form as Eq. (8) and can be straightforwardly evaluated in the discrete form of Eq. (1).
Let us illustrate the procedure considering (1) a discretization in space of the spring with a single linear Lagrangian finite element along the thickness ; (2) a forward first order finite difference to approximate the time derivative; and (3) constant material properties. For this setting, the linear system of the VDA problem at the -th time step () is
where and are the coefficients of the mass and conductivity matrices for the 1D linear lagrangian FE of the VDA problem, with
Identifying now the blocks
Eq. (13) can be rewritten as
Following this, application of static condensation to Eq. (16) to remove the unknowns of , leads to the relation
The VDA method can be applied verbatim (with enhanced accuracy), if several and/or higher order FEs are considered, since the discretization leads to the same block structure of Eq. (16). Likewise, the part-plate thermal contact in the P model can be analogously constructed; in this case, the material properties are those of the plate.
But the main advantage of the VDA concerns computational efficiency. As the method ends up extracting a parametrized closed-form expression to evaluate the equivalent boundary condition, there is no need to solve a linear system at each integration point. The additional cost with respect to the HTC model merely consists in extra storage of the temperature values of the VDA at previous time steps (depending on the time integration scheme) and extra floating point operations to evaluate Eq. (18).
The main limitation of the method is that the contour of the part-powder interface is generally a smooth 2D shape and, as a result, heat transfer across the powder-bed is not unidimensional. For simple shapes (e.g. cylinder or sphere) the local curvature of can be taken into account by considering a modified 1D heat transfer problem with
with the geometrical correction factors and defined as
where, as shown in Fig. 4, is the volume of the fraction of powder bed modelled by the VDA model and and are the areas of the external surfaces of the VDA and the part.
Another drawback is that the VDA, in spite of adding physics into the computation of the boundary condition, still needs to estimate some unknown quantities, such as , or . Assuming that and , one can get rid of estimating the HTCs, because the boundary conditions at both ends of the VDA problem become of Dirichlet type. A relation like Eq. (8) is then obtained following exactly the same procedure outlined above. Additionally, the determination of is object of discussion in Sect. 4.1.
Remark: Assuming that , and is the temperature in the middle of the VDA at the previous time step, the expressions of and of the example above in terms of the VDA parameters are
for one linear FE along the VDA:
for two linear FEs along the VDA:
for one quadratic FE along the VDA:
for one linear FE along the VDA, and :
4. Numerical experiments
4.1. Verification of the VDA against the HF model
The purpose of the first numerical experiment is to verify the new VDA formulation and demonstrate its capabilities and benefits in the thermal simulation of a SLM or DMLS process. Three model variants (cf. Sect. 2.1) are considered: (1) HF, (2) VDA-PP and (3) VDA-P. The first one is taken as the numerical reference for the other two, due to its higher physical representativeness.
The example is designed to be simple, but suitable enough to compare the accuracy and efficiency of models (2) and (3) with respect to model (1). According to this, the object of simulation is the printing of a 10x10x10 cube on top of a 110x110x20 metal substrate. Different materials are selected for the cubic sample and the build plate: while the former is made of Maraging steel M300, the latter is composed of Stainless Steel (SS) 304L. Bulk temperature-dependent thermal properties of SS 304L  are represented in Fig. 5. On the other hand, Tab. 1 prescribes constant density and specific heat values for M300, due to lack of temperature-dependent data, and the M300 powder is assumed to have 54 % of relative density.
The printing process considers a layer thickness of 30 . Therefore, a total number of 333 layers are deposited to build the sample. The values of the remaining process parameters are detailed in Tab. 2. Note that there is no information regarding the scanning path, because the simulation considers a layer-by-layer deposition sequence. This means that the printing of a complete layer is simulated in a single time step. As a result, the simulation follows an alternating sequence in time consisting of:
Printing step: A new layer is activated, i.e. added into the computational domain, and the problem is solved applying the energy input necessary to fuse the powder of the whole layer. The time increment is given by the scanning time.
Cooling step: The problem is solved without heat application to account for the time lowering the plate, recoat time and laser relocation. Here, the time step is the recoating time.
Fig. 6 shows the FE discretizations of the three tested models. All three cases consider structured meshes of varying size to adequately capture the physics of the process. In particular, a finer layer-conforming mesh is prescribed for the fabricated part, whereas a coarser mesh is used away from the printing region. As observed in Tab. 3, the mesh size decreases one order of magnitude from the HF model to the VDA-PP and VDA-P ones.
|Model||Mesh size [DoFs]||Simulation time -[%]|
|HF||274.5k||17.6 - 100|
|VDA-PP||66.5k||1.3 - 7.4|
|VDA-P||49.9k||0.8 - 4.5|
Concerning the boundary conditions, the top surfaces of the three models are subject to heat transfer through the surrounding air with an HTC of 10 [W/m2C] and air temperature of 20 [C]. On the other hand, the same HTC and environment temperature are assigned at the lateral and bottom walls of the build chamber.
The last ingredient is to characterize the BCs at the solid-powder and part-plate contacting surfaces, by establishing the VDA parameters for the VDA solid-powder (virtual powder) and the VDA part-plate (virtual plate). Recall that only the first applies to the VDA-PP model, whereas both apply to the VDA-P one.
Both VDAs adopt the full Dirichlet variant described at the end of Sect. 3 with a single linear Lagrangian FE and without geometrical correction factors. This means that the only parameters that need estimation are the virtual domain thicknesses and and the virtual domain environment temperature .
The strategy for such estimation is based on simple calibration with respect to the HF model. Let us explain the procedure step-by-step to model the virtual powder (analogous for the virtual plate). First, the temperature distribution of the HF model is analysed to identify the region concentrating the strongest gradients. As shown in Fig. 7 for the virtual powder, at about 14 away from the part, the temperature drops to 90 [C]. Further away, the thermal gradients are apparently smoother. Therefore, is set to 90 [C] and is initially approximated by 14 . Calibration with respect to the thermal history of the HF model at a selected point follows to correct to the final value. VDA parameter values obtained with this approach are listed in Tab. 4.
The numerical experiments for this example are supported by the in-house research software COupled MEchanical and Thermal (COMET)  and GiD [9, 21] as a pre- and postprocessing software. Fig. 8 compares the temperature evolutions of the reduced models against the reference HF model, at a point located in the middle of the bottom surface of the part. As observed, both VDA-PP and VDA-P are capable of reproducing the thermal response of the HF model with errors bounded by 10 % and 20 %. The VDA-P is clearly a little less accurate than the former, as expected. However, as seen in Tab. 3, the computational running times of the VDA models are one order of magnitude less than the HF model. This showcases the ability of the new formulation to approximate the response of a high-fidelity model, but with significantly increased efficiency, and it opens the path for larger scale and experimentally-based simulations, such as the one in the next section.
4.2. Verification and validation against physical experiments
4.2.1. Experimental campaign
An experimental campaign took place at the Monash Centre for Additive Manufacturing (MCAM) in Melbourne, Australia, with the purpose of (1) calibrating experimentally the thermal analysis FE framework described in Sect. 2 and (2) assess the novel VDA model presented in Sect. 3.
The printing system employed for the experiments is the EOSINT M280 from Electro Optical Systems (EOS) GmbH. It uses an Yb-fibre laser with variable beam width and power up to 400 [W]. The printing process is carried out in a closed 250x250x325 chamber subject to a laminar flow of argon that prevents oxidation.
The printed specimen is an oblique square prism with the lower base located in the centre of the building plate and a 45-degrees slope, as shown in Fig. 9. Cross section dimensions are 30x30 and the height is 80 [mm].
Eight thermocouples for temperature measurements are inserted into 0.78 [mm] diameter holes at the upward and downward facing lateral surfaces of the prism (CH1-4 and CH7-8) or in the powder bed (CH5-6). The position of the channels is indicated in Tab. 5 and Fig. 10. The printing job was interrupted at 20.58 [mm] height to install the thermocouples on a set of supporting structures that were printed together with the prism, as shown in Fig. 11.
K-type thermocouples and a Graphtec GL-900 8 high-speed data-logger are used for data gathering. The sampling rate of the data logger is 1 and the time constant of the thermocouples is 7 . As thermocouples are not welded and can move inside the hole, their measurements can be perturbed.
The process parameters used for the printing job are described in Tab. 6. The layer thickness is set to 30
. This means that 647 layers are deposited in about 3.5 [h] to build the samples. As observed, the recoat and laser relocation time varies between odd and even layers.
|Recoat and relocation time (odd layers)|
|Recoat and relocation time (even layers)|
Regarding the scanning strategy, the laser travels along the y direction back and forth for each layer, as shown in Fig. 11. Note that the number of hatches drawn does not correspond to the actual number of hatches, which is a much higher value, according to the power beam size.
The printed samples are made of Ti6Al4V Titanium alloy. The temperature-dependent properties of the bulk material, covering the range from room temperature to fusion temperature, are available in . The base plate is made of CP Ti, a material with similar thermal properties as those of Ti64. Complementary experiments in  estimated the relative density of Ti64 powder at around 54 % of the density of the bulk material at room temperature.
Fig. 12 describes the experimental data gathered at the eight thermocouple channels. During the printing of the first layers, the thermocouples are close to the laser and a sharp and highly oscillatory temperature build-up takes place. This trend stabilizes a little before the hundredth layer, when the thermocouples are far enough from the laser spot. Then, the temperature evolution enters a quasi steady-state regime, until the process finishes and the temperature falls to room temperature.
The numerical experiments were supported by a framework consisting of: (1) FEMPAR , an advanced object-oriented parallel FE library for large scale computational science and engineering, (2) Dakota , a suite of iterative mathematical and statistical methods for parameter estimation, sensitivity analysis and optimization of computational models, and (3) TITANI 7, a High Performance Computer (HPC) cluster at the Universitat Politècnica de Catalunya (Barcelona, Spain).
Such advanced computational framework, integrating FEA tools with scientific software for parameter exploration on a HPC platform, has not been previously observed in the literature related to the numerical simulation of metal AM processes, but it is very convenient to carry out the verification and validation tasks at the scale of the experiment.
|Computing nodes||5 DELL PowerEdge R630|
|CPUs||2 x Intel Xeon E52650L v3 (1.8GHz)|
|Cores||12 cores per processor|
|Local disk||2 x 250 GB SATA|
An overview of the procedure followed during the numerical tests is:
Implementation of an interface to communicate Dakota with FEMPAR.
Design and implementation of a physically accurate reference (HF) heat transfer model.
Calibration and validation of the HF model against experimental data generated at the MCAM research centre.
Implementation, calibration and validation of two additional HTC-PP and VDA-PP reduced HPC models in TITANI.
According to this, three different numerical models were tested, the only difference among them being how heat loss through the powder bed is accounted for. In the HF model, the purpose is to maximize the accuracy of the model, by including the powder bed into the computational domain of analysis, to establish a reference simulation for the reduced models.
In the HTC-PP model, the powder-bed is excluded from the computational domain and heat loss through the powder bed is modelled with a constant heat conduction boundary condition at the solid-powder contact surfaces. The VDA-PP model adopts the same hypotheses of the HTC-PP model, except for the computation of the HTC at the solid-powder interfaces, which is derived from Eq. (18), considering a single quadratic element to solve the VDA 1D heat conduction problem.
The calibration is