Efficient thermal simulation of large-scale metal additive manufacturing
using hot element addition
B. Nijhuis
⇑, H.J.M. Geijselaers, A.H. van den Boogaard
Chair of Nonlinear Solid Mechanics, Faculty of Engineering Technology, University of Twente, Drienerlolaan 5, 7522NB Enschede, the Netherlands
a r t i c l e i n f o
Article history:Received 13 September 2020 Accepted 11 December 2020 Available online 29 December 2020 Keywords:
Thermal analysis
Discontinuous Galerkin method Wire and arc additive manufacturing (WAAM)
Directed energy deposition
a b s t r a c t
Directed energy deposition processes can reduce material waste and manufacturing time of large metal parts through near net-shape production at high deposition rates. However, the localised high heat input gives rise to undesired heat accumulation, residual stresses and distortions. In this work, a fast thermal model is developed to aid in predicting and preventing these drawbacks by providing insight in the rela-tion between process settings, deposirela-tion strategy and thermal response. Material addirela-tion and heat input are efficiently combined by adding new elements at elevated temperature. Newly deposited ele-ments are assigned an artificially enhanced heat capacity to match the process heat input. The discontin-uous Galerkin finite element method is used for spatial discretisation. The resulting numerical scheme is fully explicit and can be solved element-wise. Unlike previous prescribed-temperature heat input mod-els, the proposed method correctly captures the process heat input, irrespective of substrate temperature and element size and without calibration. Comparison with experimental data shows that the thermal history of a large additively manufactured part can be accurately predicted.
Ó 2020 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http:// creativecommons.org/licenses/by/4.0/).
1. Introduction
Manufacturing large metal parts could result in a high amount of material waste when using subtractive manufacturing tech-niques such as milling or turning. Casting or forging processes require expensive part specific tooling. In additive manufacturing (AM), material is added in a layer-by-layer fashion[1]. This can alleviate these drawbacks by producing near-net-shape parts with-out any part specific tooling[2]. AM also offers an unprecedented design freedom and high levels of functional integration[3].
Within the range of AM techniques, Directed Energy Deposition (DED) processes are capable of producing large parts at high depo-sition rates[3]. These processes, however, face challenges inherent to metal AM. Localised heating and cooling causes high thermal gradients[4], which leads to residual stresses and distortions due to constrained thermal expansion[5]. Residual stresses added to service loads could negatively affect the life time of the component
[6]. Distortions may cause the final part to exceed its geometrical tolerances. Heat accumulation could lead to non-uniform beads
[7]. Differences in cooling rate at different locations or partial remelting upon depositing new material additionally gives rise to anisotropic and heterogeneous material properties[8].
A good understanding of the aforementioned phenomena is required to produce parts that meet their design requirements. Even though this can be done on a trial-and-error experimental basis, a comprehensive and fast computational model could iden-tify the final part properties for various process settings more easily. This allows process optimisation and provides an insight into the evolution of relevant physical quantities that might be dif-ficult to obtain experimentally. As a part’s microstructure, mechan-ical properties and residual stresses and distortions are mainly determined by its thermal history, an accurate and fast thermal simulation model for AM is of great importance[9].
Conventional thermal simulations for DED processes utilise a Lagrangian finite element model (FEM) in combination with a moving heat source formulation. This, however, limits such simu-lations to small parts due to the size dissimilarity between the heat source and part geometry, localised heating effects and highly non-linear material behaviour[6,10]. Several attempts have been made to decrease the computational time of the FE process models by reducing the size of the system of equations to be solved. Adaptive remeshing strategies utilising the character of the AM process have been proposed for this purpose. As thermal gradients decrease at increasing distance from the heat source, layer-wise coarsening strategies have been developed in which layers furthest from the heat source are being merged[11,12]. Another option to reduce
https://doi.org/10.1016/j.compstruc.2020.106463
0045-7949/Ó 2020 The Authors. Published by Elsevier Ltd.
This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).
⇑Corresponding author.
E-mail address:b.nijhuis@utwente.nl(B. Nijhuis).
Contents lists available atScienceDirect
Computers and Structures
the size of the system of equations is to coarsen the mesh of the substrate material away from the deposited material[13].
A second possibility of reducing the simulation effort is to limit the number of time steps that needs to be simulated. In[14], the WAAM deposition of a multi-layer straight wall is modelled using a steady-state approach. The material flows through a Eulerian mesh attached to the heat source, allowing to solve each layer deposition in a single step. The authors in[14]compared the per-formance of their steady-state model to a transient thermal simu-lation and reported a reduction in computational time of 99.7%. In arbitrary complex products, however, a steady-state solution may not be reached due to differences in part cross-section or boundary effects at the start or end of the layer deposition. As adapting Eule-rian models to complex geometries is difficult[15], this method is less suitable for use in a general thermal modelling framework.
Few works deal with the implementation of alternative discreti-sation techniques to predict the thermal history of an additively manufactured part. The DMD process has been modelled with the finite volume method (FVM) in[16]. A two-dimensional finite difference method (FDM) has been used to simulate the SLM pro-cess in[17]. Unlike the FEM, both the FDM and FVM can be used to obtain an explicit system of equations for the spatial discretisa-tion of the heat equadiscretisa-tion[18]. Combined with an explicit time inte-gration scheme, the thermal history can thus be updated without the necessity of solving equations. Despite this potential gain in computational efficiency, tracking the movement of the heat source with the required boundary condition updates remains an expensive task for large parts.
The amount of boundary condition updates can be alleviated by a simplified heat input modelling strategy, in which newly depos-ited elements are added at elevated temperature[19]. Even though the thermal history directly around the newly added material is incorrect, this rarely influences the subsequent thermal history after the heat source has passed[20]. It is of greater importance to input the correct amount of heat in appropriate spatial and tem-poral dimensions[20].
Prescribed-temperature heat input models in computational weld modelling (e.g. [19–22]) as well as AM process modelling (e.g. [23,24]) are typically implemented using the FEM. This method however requires continuity of the temperature field across element boundaries. Raising the temperature of all nodes of a newly activated element as in[21]will also increase the tem-perature in adjacent elements. As a result, the heat input is over-predicted. Raising the temperature of only the free nodes of the deposited element as in e.g.[24]underpredicts the heat input. This makes the FEM less suitable for a simulation which uses hot ele-ment addition as heat input.
A FDM implementation of hot element addition was proposed in[25]. The temperature of the eight nodes corresponding to the volume being deposited is raised to a predefined value. Despite sig-nificant time savings, this method will also overpredict the heat input as raising the nodal temperatures around the deposited vol-ume also affects the temperature of the already present material. Additionally, the FDM is not well suited for discretisation of com-plex geometries.
The objective of this work is to develop a simulation framework that allows efficient prediction of the thermal history of large addi-tively manufactured metal parts. A hot element addition strategy will be used as heat input model. The discontinuous Galerkin finite element method (DGFEM) will be used for the spatial discretisation of the thermal problem. By its nature, this method allows discon-tinuous temperatures across element boundaries. This removes the initial interpolation errors that are inherent to methods requir-ing a continuous temperature field across element boundaries. Unlike prescribed temperature heat input models in literature
(e.g.[19–25]), the developed method allows a predictable and con-stant heat input in simulation of multilayer deposition.
The outline of the paper is as follows. In Section2, the theoret-ical framework for the simulation model will be established by consecutively discussing the DGFEM spatial discretisation and the temporal discretisation of the transient heat equation, as well as the formulation of the heat input model. The numerical imple-mentation of the model will be described in Section3. In Section4, numerical tests will be conducted to check the validity of the dis-cretisation method and heat input model. The model will be vali-dated against experimental data from literature in Section5.
2. Theory
In this section, the spatial discretisation of the transient heat equation using the discontinuous Galerkin method is derived, fol-lowed by the temporal discretisation of the resulting system of equations and finally the discussion of a hot-element addition heat input model.
Consider a three-dimensional domain X bounded by dX and
having an initial temperature distribution T xð ; 0Þ. On dXN, the heat flux is prescribed as qN. On dXHheat is transferred to the surround-ings at T1with a heat transfer coefficient h. On dXD, the surface temperature is prescribed at a value TD. The temperature evolution inX is governed by
q
cp_T ¼$
ðj
$
TÞ8
xj x 2 X ð1Þn
j
$
T¼ qN8
xj x 2 dXN ð2Þn
j
$
T¼ h T Tð 1Þ8
xj x 2 dXH ð3ÞT¼TD
8
xj x 2 dXD ð4Þwhere
q
; cpandj
denote the material’s density, specific heat capac-ity and thermal conductivcapac-ity, respectively. These material parame-ters are generally temperature-dependent. n is the outward normal on dX. To solve Eq. (1) using the DGFEM, it is split into two first-order linear equations[26]q
cp_T ¼$
q ð5Þq¼
j
$
T ð6Þ2.1. Spatial discretisation
The domain is spatially discretised by subdivision into n non-overlapping elementsXi. Each element has a unique set of nodes. Unlike in a conventional FE-discretisation, where elements are connected through shared nodes, multiple nodes exist at the same coordinate point on shared element edges in the DGFEM. This allows a discontinuous solution field, but requires a method for inter-element communication.
To this end, numerical fluxes q and Tare introduced for the heat flux and temperature, respectively. These fluxes define a unique interface value at a coordinate point located at the connect-ing face between two adjacent elements. For a point on the face f of element i that is adjacent to element j, they are defined using the general notation[27]: T¼ ðC11þ C12ÞTifþ Cð 11 C12ÞTjf qk¼ ðC21 C22Þqki;fþ Cð 21þ C22Þqjk;f þnk;fC20 Tif T j f ð7Þ where Ti
fand qik;fare the temperature and k th
directional component of the heat flux of element i evaluated on face f, respectively. Tjf and
qjk;fare evaluated on the face of element j that is adjacent to f. nk;fis the kthcomponent of the unit outward normal on face f of element i. The flux parameters C11; C12; C20; C21 and C22, in combination with appropriate values for qj
k;f or Tjf, will be used for weak enforcement of inter-element continuity as well as to define boundary condi-tions. This will be elaborated on in Section3.1.
The temperature, heat flux components and weight functions
v
and wkare interpolated asT x; tð Þ ¼ N x½ ð ÞT tð Þ; qkðx; tÞ ¼ N x½ ð Þqkð Þ;t
v
ð Þ ¼ N xx ½ ð Þv
ð Þ; wx kð Þ ¼ N xx ½ ð Þwkð Þxð8Þ
where N½ is a row vector of nodal interpolation functions and T; qk;
v
; wkare column vectors with nodal values of T; qk;v
and wk, respectively.The Galerkin method, using the same interpolation for the pri-mary variables and weight functions, is applied to Eq.(5)and the directional components of Eq.(6), leading to a set of integral equa-tions. It is assumed that thermal parameters are constant over an element, such that these can be taken outside the integrals. Inte-grals over element boundaries dXiare written as a summation of surface integrals over the element faces dXf
i for the total number of faces nf. After substituting the definitions of Eqs. (7) and (8), the following equation for the kth directional component of the nodal heat fluxes of element i is obtained:
RRR Xi N i h iT Ni h i dV qi k¼
j
RRRX i @ N½ iT @xk N i h i dV Tij
X nf f¼1 Cf11þ C f 12 ZZ @Xf i Ni h iT Ni h i nkdS Tij
X nf f¼1 Cf11 C f 12 ZZ @Xf i Ni h iT Nj h i nkdS Tj ð9Þand the semidiscrete expression for the temperature rate of ele-ment i becomes:
q
cp RRR Xi N i h iT Ni h i dV _Ti¼ X3 k¼1 RRR Xi @ N½ iT @xk N i h i dV qi k X 3 k¼1 Xnf f¼1 Cf21 C f 22 ZZ @Xf i Ni h iT Ni h i nkdS qik X 3 k¼1 Xnf f¼1 Cf21þ C f 22 ZZ @Xf i Ni h iT Nj h i nkdS qjk X nf f¼1 ZZ @Xf i Cf 20 N i h iT Ni h i dS Ti þX nf f¼1 ZZ @Xf i Cf20 N i h iT Nj h i dS Tj ð10ÞIn Eqs.(9) and (10), the interpolation functions Nh ii belong to ele-ment i. The interpolation functions Nh ij belong to the element
adja-cent to the face f over which the surface integral is taken. The superscripted indices f for the flux parameters Cf
11; C f 12; C f 20; C f 21and Cf
22 indicate these parameters belong to face f, as each face is assigned an individual set of flux parameters depending on its func-tion (inter-element communicafunc-tion or boundary condifunc-tion enforce-ment). The terms inside the integrals are replaced by the matrices they represent to form the final system of semi-discrete equations:
Mi h i qi k¼
j
S i k h i Tij
X nf f¼1 Cf 11þ C f 12 Gii fk h i Ti þ Cf 11 C f 12 Gijfk h i Tj ð11Þ Mi h i _Ti¼ 1q
cp X3 k¼1 Sik h i qi k 1q
cp X3 k¼1 Xnf f¼1 Cf21 C f 22 Giifk h i qi k þ Cf 21þ C f 22 Gij fk h i qj k 1q
cp Xnf f¼1 Cf20 L ii fk h i Ti Lij fk h i Tj ð12ÞIt has to be noted that the number of degrees of freedom resulting from the discontinuous Galerkin spatial discretisation of the heat equation is larger compared to that of a standard finite element dis-cretisation. It will however be shown that, in combination with the temporal discretisation scheme described in Section2.2, a system of equations is obtained that can be solved very efficiently.
2.2. Temporal discretisation
The spatial discretisation of the heat equation in Section 2.1
results in a local and explicit semi-discrete system of equations. To retain these favourable properties in the temporal discretisa-tion, the explicit forward Euler time integration scheme is imple-mented. The updated nodal temperatures of element i at time step nþ 1 are calculated from the temperature rate at time step n as
Tinþ1¼ Tinþ
D
t _Tin ð13Þin whichDt is the time step size. As an explicit time integration method is chosen, the maximum time step size must satisfy the Courant-Friedrich-Lewy (CFL) condition for numerical stability. For the problem at hand, the time step limit is calculated as[18].
D
t6 C min i;k;Tq
ð ÞcT pð ÞTj
ð ÞTD
xik ð14Þ whereDxikis the dimension of an element i in the k
thdirection and the parameter C depends on the spatial discretisation settings[27].
2.3. Heat input modelling
In this work, the expensive boundary condition updates required for simulating a moving heat flux distribution, are allevi-ated by a simplified heat input model based on prescribed temper-atures. This approach has been used in computational weld modelling[28,21,22]and AM process modelling[24,25]. In these references, the temperature of nodes of newly deposited elements is set at Tdand held constant for a certain amount of time td. After this time, the material is allowed to cool according to Eq.(1). The holding time tdis typically set to the time the heat source needs to pass the deposited elements[25,21,22], but may also be fitted to match the process heat input[28,19]. Deposition temperatures are set to the liquidus temperature[28,21,24], molten zone
perature [25] or used as fitting parameter to match the process heat input[22].
Maintaining the elevated temperature for a prolonged time effectively results in two separate heating events [22]. The heat retained in the hot, deposited volume is added to the domain once the corresponding elements are being activated. Furthermore, the domain is heated by the prescribed temperature at the interface between deposited and already present material during td. The lat-ter heat input is dependent on the base malat-terial temperature and deposited bead length per activation step[22]. As both may change during an AM process simulation, replicating a constant process power is not generally possible. Therefore, the prescribed temper-ature heat input model developed in this work will be made inde-pendent of substrate temperature and activation strategy.
To eliminate prolonged surface heating at the interface between newly deposited and already present material, the deposited mate-rial is allowed to cool immediately. Newly deposited elements are assigned a constant, artificially enhanced heat capacity c
pin the temperature range between solidus (Ts) and deposition tempera-ture (Td) to be able to represent the relatively high heat input of large scale metal AM processes. After a material point has cooled down below Ts, its heat capacity is changed to the physical, possi-bly temperature-dependent, value.
The total amount of heat input per unit bead length by deposit-ing material at a deposition temperature Tdis
QA¼ A RTd T0
q
cpð ÞdTT ¼ A RTs T0q
cpð ÞdT þT RTd Tsq
c pdT ð15Þwhere T0 is a reference temperature and A is the bead’s cross-sectional area. Eq.(15)is independent of the length of the deposited bead segments and of the temperature of the material underneath. The proposed heat input model can thus be used to resemble a con-stant process power throughout the entire build. The value for cp can easily be calculated from the process heat input per unit bead length QA;pfor a given deposition temperature as
cp¼ QA;p=A RTs T0
q
cpð ÞdTTq
ðTd TsÞ ð16Þ 3. Numerical implementationThe nature of the discretised system of equations allows several simplifications in the process modelling of additive manufacturing. The definition of the numerical fluxes makes it straightforward to implement appropriate thermal boundary conditions. This will be addressed in Section 3.1. By carefully identifying solution-independent terms in Eq.(11) and (12), an efficient numerical solu-tion procedure can be devised as will be shown in Secsolu-tion3.2. It will also be discussed that the solution scheme provides additional benefits in modelling material addition. The model described in this section has been implemented in MATLAB[29].
3.1. Numerical fluxes and thermal boundary conditions
As mentioned in Section2.1, the numerical fluxes as defined in Eq.(7)are used for weak enforcement of inter-element continuity of the solution as well as for boundary condition implementation. For inter-element interpolation, a Bassi-Rebay flux with stabilisa-tion parameter
g
¼ 4 is implemented [30]. This ensures that the null-space of the semidiscrete system of equations only contains the constant solution regardless of the number of elements in the domain [30]. In this way, oscillations in the solution areavoided. For implementation of boundary conditions, qj k;f or Tjf are used to prescribe known quantities at the boundaries. The parameters for different boundary conditions are given inTable 1. Radiative heat losses to the surroundings at T1are modelled as convective boundary conditions with a temperature-dependent heat transfer coefficient. By rewriting the Stefan–Boltzmann law
qrad¼
r
B T4f T 4 1 ¼r
B T3fþ T 2 fT1þ TfT21þ T31 Tf T1 ð Þ ð17Þin which
is the material’s emissivity,r
Bis the Stefan-Boltzmann constant and Tf is the surface temperature of the material, an expression for hradð Þ is obtainedTfhradð Þ ¼Tf
r
B T3fþ T 2fT1þ TfT21þ T31
ð18Þ
Combined convective-radiative heat losses are implemented using
h Tð Þ ¼ hs radð Þ þ hTf 1 ð19Þ
where h1 is the heat transfer coefficient for convection to the surroundings.
3.2. Solution procedure
A flowchart for the developed thermal model is provided in
Fig. 1. Relevant parts of the solution algorithm will be discussed in Section3.2.1–3.2.3.
3.2.1. User input
Prior to the simulation, a mesh needs to be created for the final part. This can be done using any preferred FE-mesh generator. Each element in the DGFEM discretisation should however have a dis-tinct set of nodes (i.e. no shared nodes may exist at inter-element boundaries). Therefore, the outcome of a standard mesh generator should be postprocessed such that each node at an inter-element boundaries is duplicated. Each of the duplicated nodes is then assigned to only one element.
The deposition sequence and speed is included in the simula-tion by means of a vector with activasimula-tion times tact. The simulation starts at t¼ 0. All elements with activation times t < 0 will be regarded as substrate elements and are initialised at a predefined substrate temperature. All elements with activation times tP 0 will be deposited using the approach described in Section3.2.3.
Boundary conditions are assigned by defining planes in three-dimensional space. For each of these planes, a certain prescribed temperature, prescribed heat flux or convective boundary condi-tion can be specified as input. The solver determines whether the centroid of an element face coincides with one of these planes and assigns the respective boundary condition. Convective-radiative heat transfer tot the environment is modelled automati-cally on all remaining external faces of the part.
3.2.2. Solver initialisation
Element connectivity data is needed for inter-element commu-nication, assigning boundary conditions and calculation of the ele-ment system matrices. Connectivity data is determined on the mesh of the final, fully deposited, geometry prior to the simulation. At non-connected faces, appropriate boundary conditions need to be considered. As discussed in Section3.1, for a face f subjected to a boundary condition, known quantities at the boundary need to be stored in qj
k;f or Tjf. These are evaluated in the element adja-cent to f, which, however, does not exist for boundary faces.
There-fore, placeholder elements are introduced for all distinct boundary conditions in the model. Their nodal values qi
kand T
iare set to the respective prescribed quantities as inTable 1. All boundary faces are virtually connected to the appropriate placeholder element.
As not all elements are initially active, additional non-connected faces will appear that cannot be determined from con-nectivity data of the final geometry. These faces will be exposed to the surroundings and are assigned the appropriate convection-radiation boundary condition. As seen inTable 1, the temperature of the elements adjacent to these faces should be set to the ambi-ent temperature. By initialising all to-be-deposited elemambi-ents at T1, this requirement is automatically fulfilled. Therefore, no updates to the connectivity array or element matrices are required during the
simulation. This would have been necessary if these non-connected faces needed to be virtually non-connected to the ambient air placeholder element.
Knowing the connectivity data, the element system matrices in Eqs.(11) and (12)can be calculated. As the thermal parameters and flux constants were taken outside the integrals in Eqs.(9) and (10), these matrices are constant and only depend on the nodal coordi-nates and interpolation functions. This allows to pre-multiply the right hand sides of Eqs.(11) and (12)with the inverses of Mh ii at the element level at the start of the analysis. As a result, heat fluxes and temperature increments can be calculated by explicit matrix-vector multiplication during the solution loop.
3.2.3. Solution loop
The solution algorithm increments the time by a step sizeDt until the final time has been reached. Two functions need to be ful-filled: Simulation of material addition and thus updating the boundary conditions as well as updating the temperature distribu-tion. Where the latter action is to be performed every time-increment, simulation of material addition and boundary condition updating only has to be executed upon a change in active elements. In every time step, the updated active setAnþ1is compared to the active set in the previous step An. The difference between the current and previous active set Ad¼ Anþ1n An is the set of newly deposited elements. Nodal temperatures for all elements in Ad are set to the deposition temperature Td, directly linking the heat input modelling to the material addition procedure. Thus, no boundary condition updates other than those explained in Sec-tion3.2.3are required. This makes the simulation of the process heat input significantly less complex than when a moving heat source distribution is used.
In conventional FEM AM process models, quiet or inactive ele-ment methods are common for modelling material addition[31]. In the former, the system of equations is set up for all elements in the final geometry, but non-deposited elements are assigned material properties such that they do not influence the solution. In the latter method, elements are only introduced into the model from the moment they are deposited. As the update of the temper-ature field in this work is performed by a fully explicit, element-wise scheme on all active elements, the newly deposited elements can easily be included in the simulation. Similar to the FDM AM process model in[25], the model developed in this work combines benefits of both approaches: initialisation of the total system prior to the analysis as in the quiet element method and solving only the active part of the domain as in the inactive element method.
In the first time step or upon a change in active set, flux param-eters are assigned to the boundary faces of active elements as in
Table 1. On faces with radiative heat losses, C20 is temperature-dependent and will be calculated from Eq.(18)or(19)every time step. The surface temperatures for radiative heat transfer are approximated by element average temperatures, which are calcu-lated every time step to evaluate temperature-dependent material properties as well.
For all active elements, the nodal heat flux components are now calculated from the current temperature field using Eq.(11). With
Fig. 1. Flowchart of the DGFEM thermal solver. Table 1
Definition of numerical fluxes for various types of element boundary conditions.
Boundary condition C11 C12 C21 C22 C20 qj
k;f Tjf
None (inter-element interpolation) 0.5 0 0.5 0 g¼ 4 qj
k;f Tjf
Prescribed face temperature TD 0.5 0.5 0.5 0.5 g¼ 4 0 TD
Prescribed heat flux qN 0.5 0.5 0.5 0.5 0 qN 0
Convection to medium at T1
with heat transfer coefficient h1 0.5 0.5 0 0 h1 0 T1
this result, the nodal temperature rates are obtained using Eq.(12)
and the updated temperature field is determined using Eq.(13).
4. Numerical examples 4.1. Rod cooling
The main motivation to use the DGFEM for discretising the heat equation in this work, is its ability to resolve initially discontinuous temperature fields and the resulting heat transfer. This capability is validated in a one-dimensional analysis. A rod of 1 m length is con-sidered. All thermal parameters are set to unity. The first 70 per-cent of its length is initially at zero temperature, the remainder of the geometry is initialised at unit temperature. The fully isolated domain is discretised with 40 equal-sized, piecewise-linear ele-ments. The total simulated time is 1 s, the time step size is set to the critical time step size according to Eq. (14) with C¼ 0:08. The temperature evolution in the domain will be tracked and com-pared against the analytical solution for one-dimensional heat con-duction in an isolated rod.
The analytical solution for the initial temperature profile / xð Þ
/ð Þ ¼x 0 ; 0 6 x < bL 1 ; bL 6 x 6 L ð20Þ is given as[32] Tanðx; tÞ ¼ 1 bð Þ þX1 n¼1 2
l
nL sinl
nL sinl
nbL cosl
nxeal 2 nt ð21Þwhere b¼ 0:7; L ¼ 1, the thermal diffusivity
a
¼j
=q
cp¼ 1 andl
n¼ np
=L.Fig. 2shows the comparison between the numerical and analyt-ical solutions for the temperature distribution in the rod at t = 5 105, 0.05, 0.15 and 1 s. A perfect agreement with the ana-lytical solution is observed. This is quantified by the relative L2 errors of 10.23%, 0.11%, 0.037% and 3.96 105% at t = 5 105, 0.05, 0.15 and 1 s, respectively. These errors are calculated from the differences between the numerical (T) and analytical
tempera-ture fields, evaluated at the nodes, as
L2ð Þ ¼t ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P i T xð i; tÞ Tanðxi; tÞ 2 =PiðT xði; tÞÞ2 q
. The error at the start of the simulation is caused by the presence of initial inter-element discontinuities, which rapidly decrease as the solution proceeds.
The final temperature of the domain in the limit t! 1 equals 0.3. The results inFig. 2show that the numerical solution indeed converges to this steady-state value, indicating the conservation property of the numerical method. The good resemblance with the analytical solution at the different times shows that the correct amount of heat is transferred from the hot domain into the cold region of the model at the appropriate rate. The DGFEM thermal solver can thus be used to correctly model heat input using a hot-element addition strategy. It is emphasised that this result cannot readily be obtained using a standard finite element discreti-sation because of the initial discontinuity in temperature. 4.2. Single bead depositions
The heat input model proposed in this work (termed energy-based model in the following) was developed to add a constant heat input power to the build domain, regardless of the element activa-tion strategy and base material temperature. The invariability of the heat input under different activation strategies will be assessed by performing single bead simulations. Results will be compared against the common implementation of a prescribed-temperature heat input model that maintains the deposition temperature for the time the heat source needs to pass the deposited material as in[25,21,22](termed temperature-based model in the following).
For a single bead deposited at constant speed
v
d, multiple acti-vation strategies can be devised. The deposited bead segment length per activation step Ld and the inter-activation time Dtd should be chosen such that Ld=Dtd¼v
d. Long bead segments with longer inter-activation times can thus describe the same deposi-tion speed as shorter segments with a smaller inter-activadeposi-tion time, as shown schematically inFig. 4.Beads with a cross-sectional area of 1 1 mm2and a length of 12 mm will be deposited onto a 20 5 2 mm3
substrate. Ldis varied between 0.5 and 4 mm while settingDtdto match a deposi-tion speed of 8.33 mm/s. The domain as shown inFig. 3is meshed with piecewise-linear hexahedral elements and a uniform edge size of 0.25 mm. The domain is fully isolated. The simulated time is 5 s and a constant time step of 0.1 ms will be used in the simu-lations. The material is assigned temperature-independent thermal parameters for mild steel evaluated at room temperature, with a density
q
¼ 7860 kg/m3, thermal conductivityj
¼ 52 W/(mK) and specific heat capacity cp¼ 480 J/(kg K)[33]. The solidus tem-perature Ts¼ 1450 °C, the deposition temperature Td¼ 2000 °C to resemble the temperature in the molten zone[34]. The substrate is initially at T0¼ 0 °C.The heat input is calculated in each step n for all currently active elements as Qn¼ X i2An
q
cpVi Tin Tin1 ð22Þwhere Vi is the element volume and Ti
nand Tin1are the element average temperatures at the end of the current and previous time
Fig. 2. Numerically and analytically predicted temperature distributions at
step, respectively. Ti
0¼ 0 for all elements.
Fig. 5shows the evolution of the heat added to the domain over time for the various values of Ld for both models. The total heat input for the various simulations is plotted against Ld in Fig. 6. The energy-based heat input model was calibrated such that its heat input matches that of the temperature-based model with
Ld¼ 1 mm by calculating the artificial heat capacity cp¼ 6243 J/ (kgK) according to Eq.(16). It is observed inFig. 6that the total heat input of both models indeed agrees for this activation length. As directly visible fromFig. 6and also observed from the slopes of the curves inFig. 5, the heat input power for the temperature-based model is dependent of Ld, whereas that for the energy-based model is not. The heat retained in the bead deposited at Td, however, is constant. The size-dependent heat input of the temperature-based model is thus a direct result of holding the deposition temperature forDtd. For longer activation lengths, a lar-ger substrate surface area is heated for a lonlar-ger time compared to shorter activation lengths. This causes the heat input of the temperature-based model to increase with Ld.
Despite the large steps in heat input observed inFig. 5a for the longer activation lengths, the thermal histories at some distance from the deposited material are not significantly affected by the intermittent heat input.Fig. 7shows the temperature results of a probe point located at (6, 0, 2) for three different activation lengths. Even with the measurement location being located only two times the layer thickness from the deposited material, the thermal histories are in good agreement. This shows that the energy-based model, that cannot capture the exact distribution of the heat input, can still be used to predict the overall thermal response in a large-scale metal additive manufacturing simulation. The energy-based model is therefore well suitable for thermal modelling of additive manufacturing. Through calculation of c p using Eq.(16), the desired heat input can easily be specified. Ld
Fig. 4. Different activation strategies for a constant deposition speed. Left: short segments, short inter-activation times, right: longer segments, longer inter-activation times.
Fig. 5. Evolution of heat input versus time for both heat input models at Ld¼ 0:5, 2
and 4 mm. Fig. 6. Comparison of total heat input for both heat input models at various Ld.
andDtddo not need to be constant to obtain a constant heat input power for a given deposition speed, as is the case for the temperature-based model. For a constant ratio Ld=Dtd and fixed values of c
pand TD, constancy of the heat input is automatically fulfilled. Even if this ratio is non-constant, the desired heat input can be regained by recalculating cp. This makes that irregular meshes and varying deposition speeds can be handled in a straightforward manner.
4.3. Straight wall depositions
Straight wall deposition simulations are performed to deter-mine the effect of base material temperature on the heat input behaviour. Ten 12 mm long beads of 1 mm width and height are deposited on top of each other to create the geometry as shown inFig. 8. Deposition occurs without interpass time and in positive x-direction for all layers. Material, discretisation and heat input model parameters and boundary conditions are identical to those in Section4.2. The total simulated time is 20 s and Ldis fixed at 1 mm/ activation step.
Fig. 9shows the evolution of the heat input for both models. Opposed to the results in Section4.2, where the total heat input for Ld¼ 1 mm was equal for both models, a difference in total heat input is observed in the wall build. The energy-based model sup-plies a larger amount of total heat to the domain. The heat input power of the temperature-based model decreases with time. Only
during the first layer deposition (t 1:44 s), the behaviour of both models is similar.
The temperature evolution of the midpoint of the substrate bot-tom face is plotted inFig. 10for both heat input models. Also here, the results for both models are in good accordance during the first layer deposition, but deviate as the simulation proceeds. The final temperature for the energy-based model is significantly higher than that for the temperature-based model because of the larger heat input of the former. It has to be noted that due to the isolated boundaries, lack of interpass time and absence of latent heat, the domain is heated well above the melting point for both simula-tions. This is however of no importance for the comparison of the heat input behaviour by the two heat input models.
The results inFigs. 9 and 10show the difference between both heat input models in terms of dependency on the base material temperature. The heat input for the energy-based model is not affected by the temperature rise of the substrate. The heat input rate of the temperature-based model decreases with increasing substrate temperature. This effect was also noted in[22] and is caused by the fact that the heat input for this model depends on
the temperature difference between newly-deposited and
already-present material. As such, the temperature-based model cannot maintain a specified heat input power when heat accumu-lates in the part. The energy-based model can however represent the specified heat input even with significant heat accumulation in the already-present material.
5. Experimental validation
The ability of the DGFEM thermal solver to efficiently simulate the manufacturing of a large metal part will be evaluated by mod-elling the WAAM deposition of a 500 mm long straight wall. The proposed heat input model will be validated by comparing the simulated thermal results against experimentally obtained data from literature[14]. The model parameters are taken from[14]
and will be briefly repeated here for clarity.
Four 5 mm wide weld beads of 2 mm thickness will be sequen-tially deposited on top of each other along the centerline of a 60 mm wide substrate of 12 mm thickness. The deposition speed equals 8.33 mm/s, the interpass time between subsequent layer depositions is 400 s. Deposition occurs in positive x-direction. As the geometry is symmetric with respect to the xz-plane, only half of the domain is modelled as shown inFig. 11. The bottom face is subjected to convective heat transfer to the surroundings at 25°C with a heat transfer coefficient of h¼ 300 W/(m2K). Convective (h¼ 5:7 W/(m2K)) and radiative (
¼ 0:2) heat transfer to the
sur-Fig. 7. Thermal history probed at (6, 0,2) for the energy-based heat input model.
roundings is modelled on all remaining faces. On the symmetry plane, isolation in x-direction is included by setting h1¼ 0 and
¼ 0.The material data is taken from [33]. Fig. 12 shows the temperature-dependent thermal conductivity and heat capacity data for S355 steel. The density
q
¼ 7860 kg/m3is assumed con-stant. Convective heat transfer in the molten zone is mimickedwith an increased thermal conductivity
j
¼ 120 W/(mK) abovethe liquidus temperature Tl¼ 1500°C. The latent heat of solidifica-tion of HL¼ 270 kJ/kg is included by artificially increasing the specific heat between the solidus temperature Ts¼ 1450°C and liq-uidus temperature using the apparent heat capacity method[35]:
ceff p ¼ cpð ÞT ; T < Ts cpð Þ þ HT L= Tð l TsÞ ; Ts< T < Tl cpð ÞT ; Tl< T 8 > < > : ð23Þ
The process heat input is 269.5 J/mm. For the given thermal param-eters, the artificial specific heat c
p is calculated as 4537.9 J/(kgK) from Eq.(16)for a deposition temperature Td¼ 2000 °C.
The model is meshed with a non-uniform mesh that is coars-ened away from the wall. The mesh geometry is set to match that of[14]as closely as possible to be able to provide an estimative comparison in simulation runtime. Piecewise-linear hexahedral
elements are used. The smallest element size of
DxDyDz¼ 2:083 0:833 0:667 mm3 is used in the wall, the largest element size of DxDyDz¼ 5 5 8:33 mm3 is used furthest away from the deposition area. Every activation step, a bead segment of 2.083 mm length is deposited. In total, the model consists of 26640 elements. The time step size is set to 0.00125 s, resulting in a total number of 1.472 million time steps to simulate the 1840 s physical build time.
The process simulation is conducted using the DGFEM thermal solver on a single core of an Intel Xeon E5-2690 v3 processor with a clock speed of 2.60 GHz. Nodal temperature data is written to MATLAB native.mat-files during the solution to avoid excessive RAM usage, requiring the solution loop to only store temperature results for a single time step in the memory. Results are written every 0.05 s during a layer deposition and every 0.5 s during the interpass time.
The temperature evolution is measured at four probe points located in the substrate, P1: 250ð ; 5; 12Þ, P2: 250; 20; 12ð Þ, P3:
250; 0; 0
ð Þ and P4: 375; 0; 4ð Þ. Thermal histories of these points are compared to experimentally obtained thermocouple measure-ments at these locations[14]inFig. 13. A good agreement between numerical temperature predictions and experimentally measured values is obtained. Both peak temperatures as well as heating and cooling rates are accurately predicted by the model for all measurement locations. This indicates that the heat input resulting from the hot element addition matches well with the process heat
Fig. 10. Comparison of temperature evolution in the substrate probed at (6, 0,2) for both heat input models.
Fig. 11. Geometry used for experimental validation of the model.
Fig. 12. Temperature-dependent thermal properties of S355[33].
input and that the heat propagation through the domain is cor-rectly captured.
The total time spent on the calculations outlined in Section3.2.3
was 11 h 49 m. This is a factor 4.3 faster than the reported runtime of 51 h 24 m for an Abaqus FEM simulation of this problem using a moving heat source distribution, performed on a quad-core grid computing system in [14]. This highlights the potential of the developed model to efficiently calculate the thermal history of a large DED build. Further improvement in efficiency is expected by running the thermal solver as a compiled and optimised code, instead of the current MATLAB implementation. Reduction of the solution time is also possible by coarsening the mesh in the wall, which in this work is taken comparable to that in[14]to enable a fair runtime comparison. A coarser mesh will alleviate the time step restriction posed by Eq.(14).
6. Conclusions
In this work, a hot element addition strategy is presented that can be used for efficient thermal simulations of large scale metal additive manufacturing using directed energy deposition. The dis-continuous Galerkin finite element discretisation of the heat equa-tion results in the numerical framework required to correctly
capture the heat transfer resulting from the initially discontinuous temperature field. The resulting element-wise solution scheme enables efficient modelling of material addition. Allowing the newly activated hot elements to cool immediately after deposition, a heat input model is obtained that is independent of base material temperature and length of the deposited bead segments. This makes the model suitable for multilayer depositions of complex parts that experience significant heat accumulation. By matching the heat retained in the hot elements to the process heat input, the model accurately predicts the thermal history of additively manufactured structures.
Declaration of Competing Interest
The authors declare that they have no known competing finan-cial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgements
The authors would like to thank Dr. Jialuo Ding from Cranfield University for providing the experimental thermal data. This research was carried out under project number P16-46/S17024K,
which is part of AiM2XL program framework of the Partnership Program of the Materials innovation institute M2i (www.m2i.nl) and the Netherlands Organization for Scientific Research (www. nwo.nl). The research was conducted in collaboration with indus-trial partners and supported by the Rotterdam Fieldlab Additive Manufacturing BV (RAMLAB) (www.ramlab.com).
References
[1] Gibson I, Rosen D, Stucker B. Additive Manufacturing Technologies: 3D Printing, Rapid Prototyping and Direct Digital Manufacturing. New York: Springer; 2015.https://doi.org/10.1007/978-1-4939-2113-3.
[2] Yang L, Hsu K, Baughman B, Godfrey D, Medina F, Menon M, et al. Additive manufacturing of metals: the technology, materials, design and production, Springer Series in Advanced Manufacturing. Cham: Springer International Publishing; 2017.https://doi.org/10.1007/978-3-319-55128-9.
[3] Milewski JO. Additive manufacturing of metals: from fundamental technology to rocket nozzles, medical implants, and custom jewelry, Springer Series in Materials Science. Cham: Springer; 2017. https://doi.org/10.1007/978-3-319-58205-4.
[4] DebRoy T, Wei HL, Zuback JS, Mukherjee T, Elmer JW, Milewski JO, et al. Additive manufacturing of metallic components – Process, structure and properties. Prog Mater Sci 2018;92:112–224. https://doi.org/10.1016/j. pmatsci.2017.10.001.
[5] Gouge M, Michaleris P. Thermo-Mechanical Modeling of Additive Manufacturing. Oxford: Butterworth-Heinemann; 2017. https://doi.org/ 10.1016/C2016-0-00317-0.
[6] Francois MM, Sun A, King WE, Henson NJ, Tourret D, Bronkhorst CA, et al. Modeling of additive manufacturing processes for metals: Challenges and opportunities. Curr Opin Solid State Mater Sci 2017;21(4):198–206.https:// doi.org/10.1016/j.cossms.2016.12.001.
[7] Wu B, Ding D, Pan Z, Cuiuri D, Li H, Han J, et al. Effects of heat accumulation on the arc characteristics and metal transfer behavior in Wire Arc Additive Manufacturing of Ti6Al4V. J Mater Process Technol 2017;250:304–12.https:// doi.org/10.1016/j.jmatprotec.2017.07.037.
[8] Kok Y, Tan XP, Wang P, Nai MLS, Loh NH, Liu E, et al. Anisotropy and heterogeneity of microstructure and mechanical properties in metal additive manufacturing: A critical review. Mater Des 2018;139:565–86.https://doi.org/ 10.1016/j.matdes.2017.11.021.
[9] Stavropoulos P, Foteinopoulos P. Modelling of additive manufacturing processes: a review and classification. Manuf Rev 2018;5.https://doi.org/ 10.1051/mfreview/2017014.
[10] Ding J, Colegrove P, Williams S, Wang F, Mehnen J, Sequeira Almeida PM. A computationally efficient finite element model of wire and arc additive manufacture. Int J Adv Manuf Technol 2014;70(1–4):227–36.https://doi.org/ 10.1007/s00170-013-5261-x.
[11] Denlinger ER, Irwin J, Michaleris P. Thermomechanical modeling of additive manufacturing large parts. J Manuf Sci Eng 2014;136(6):061007.https://doi. org/10.1115/1.4028669.
[12] Martukanitz R, Michaleris P, Palmer T, DebRoy T, Liu Z-K, Otis R, et al. Toward an integrated computational system for describing the additive manufacturing process for metallic materials. Addit Manuf 2014;1–4:52–63.https://doi.org/ 10.1016/j.addma.2014.09.002.
[13] Montevecchi F, Venturini G, Grossi N, Scippa A, Campatelli G. Finite Element mesh coarsening for effective distortion prediction in Wire Arc Additive Manufacturing. Addit Manuf 2017;18:145–55. https://doi.org/10.1016/j. addma.2017.10.010.
[14] Ding J, Colegrove P, Mehnen J, Ganguly S, Sequeira Almeida PM, Wang F, et al. Thermo-mechanical analysis of Wire and Arc Additive Layer Manufacturing process on large multi-layer parts. Comput Mater Sci 2011;50(12):3315–22.
https://doi.org/10.1016/j.commatsci.2011.06.023.
[15] Megahed M, Mindt H-W, N’Dri N, Duan H, Desmaison O. Metal additive-manufacturing process and residual stress modeling. Integr Mater Manuf Innov 2016;5(1):61–93.https://doi.org/10.1186/s40192-016-0047-2.
[16] Wang J, Shi J, Wang Y, Hu Y, Dai J, Xu K. Fast computation of thermal field of direct metal deposition: a preliminary study based on quiet element method. In: ASME 2018 13th International Manufacturing Science and Engineering Conference, vol. 4.https://doi.org/10.1115/msec2018-6696.
[17] Foteinopoulos P, Papacharalampopoulos A, Stavropoulos P. On thermal modeling of Additive Manufacturing processes. CIRP J Manuf Sci Technol 2018;20:66–83.https://doi.org/10.1016/j.cirpj.2017.09.007.
[18] Hesthaven JS, Warburton T. Nodal discontinuous Galerkin methods: algorithms, analysis, and applications. In: Texts in Applied Mathematics. New York: Springer; 2007. https://doi.org/10.1007/978-0-387-72067-8.
[19] Lindgren L-E. Finite element modeling and simulation of welding part 1: increased complexity. J Therm Stresses 2001;24(2):141–92. https://doi.org/ 10.1080/01495730150500442.
[20] Lindgren L-E. Modelling options in computational welding mechanics (CWM). Cambridge: Woodhead Publishing; 2007. p. 119–63.https://doi.org/ 10.1533/9781845693558.119[chapter 9].
[21] Lindgren L-E, Runnemalm H, Näsström MO. Simulation of multipass welding of a thick plate. Int J Numer Meth Eng 1999;44(9):1301–16.https://doi.org/ 10.1002/(SICI)1097-0207(19990330)44:9<1301::AID-NME479>3.0.CO;2-K. [22] Seleš K, Peric M, Tonkovic Z. Numerical simulation of a welding process using a
prescribed temperature approach. J Constr Steel Res 2018;145:49–57.https:// doi.org/10.1016/j.jcsr.2018.02.012.
[23] Schilp J, Seidel C, Krauss H, Weirather J. Investigations on temperature fields during laser beam melting by means of process monitoring and multiscale process modelling. Adv Mech Eng 2014;6:217584.https://doi.org/10.1155/ 2014/217584.
[24] Salonitis K, D’Alvise L, Schoinochoritis B, Chantzis D. Additive manufacturing and post-processing simulation: laser cladding followed by high speed machining. Int J Adv Manuf Technol 2016;85(9):2401–11. https://doi.org/ 10.1007/s00170-015-7989-y.
[25] Stockman T, Schneider JA, Walker B, Carpenter JS. A 3D finite difference thermal model tailored for additive manufacturing. JOM 2019;71(3):1117–26.
https://doi.org/10.1007/s11837-019-03338-6.
[26] Li BQ. Discontinuous Finite Elements in Fluid Dynamics and Heat Transfer. Computational Fluid and Solid Mechanics. London: Springer; 2006.https://doi. org/10.1007/1-84628-205-5.
[27] Li P, Dong Y, Tang M, Mao J, Jiang LJ, Bagci H. Transient Thermal Analysis of 3-D Integrated Circuits Packages by the DGTD Method. IEEE Trans Compon Packag Manuf Tech 2017;7(6):862–71.https://doi.org/10.1109/TCPMT.2017.2666259. [28] Carmet A, Debiez S, Devaux J, Pont D, Leblond JB. Experimental and Numerical Study of Residual Stresses and Strains in an Electron-Beam-Welded Joint. In: Beck G, Denis S, Simon A, editors. International Conference on Residual Stresses. Dordrecht: Springer; 1989. p. 491–6. https://doi.org/10.1007/978-94-009-1143-7_82.
[29] MATLAB, version 9.5.0 (R2018b). Natick, Massachusetts: The MathWorks Inc.; 2018.
[30] Kirby RM, Karniadakis GE. Selecting the Numerical Flux in Discontinuous Galerkin Methods for Diffusion Problems. J Sci Comput 2005;22–23 (1):385–411.https://doi.org/10.1007/s10915-004-4145-5.
[31] Lundbäck A, Lindgren L-E. Modelling of metal deposition. Finite Elem Anal Des 2011;47(10):1169–77.https://doi.org/10.1016/j.finel.2011.05.005.
[32]Carslaw HS, Jaeger JC. Conduction of heat in solids. 2nd ed. London: Oxford University Press; 1959.
[33]Ding J. Thermo-mechanical Analysis of Wire and Arc Additive Manufacturing Process [Phd thesis]. Cranfield University; 2012.
[34] Kozakov R, Schöpp H, Gött G, Sperl A, Wilhelm G, Uhrlandt D. Weld pool temperatures of steel S235 while applying a controlled short-circuit gas metal arc welding process and various shielding gases. J Phys D - Appl Phys 2013;46 (47):475501.https://doi.org/10.1088/0022-3727/46/47/475501.
[35] Hu H, Argyropoulos SA. Mathematical modelling of solidification and melting: a review. Model Simul Mater Sci Eng 1996;4(4):371–96. https://doi.org/ 10.1088/0965-0393/4/4/004.
