A computationally efficient thermo-mechanical model with temporal acceleration for prediction of residual stresses and deformations in WAAM

ABSTRACT

The present paper introduces a novel temporal acceleration strategy for computationally efficient prediction of residual stresses and deformations in wire arc additive manufacturing (WAAM) components. It employs a semi-analytical approach, dividing temperature into the analytical and the complementary fields. The analytical field is obtained by a closed-form solution, and the complementary field is employed to solve the boundary conditions. A temporal acceleration factor for the heating period is applied. Meanwhile, the diffusion time in the analytical field is manipulated to guarantee accurate temperature prediction. Validation via WAAM experiments, including both a thin-wall structure and a practical engineering component with characteristic dimensions on the order of metres, indicates that the predicted stresses is 50 MPa lower than the experimental values. Moreover, the discrepancy of between the predicted deformations and the experimental measurements is less than 10%, demonstrating reasonable accuracy can be achieved with attractive computational efficiency.

KEYWORDS:

• Wire arc additive manufacturing
• process modelling
• thermo-mechanical model
• residual stresses and deformations

1. Introduction

Wire Arc Additive Manufacturing (WAAM) is a process that operates by utilising an electric arc as the heat source to melt and deposit the filler wire, layer-by-layer, to build intricate 3D structures [1]. WAAM has gained considerable attention for its ability to fabricate large-scale metallic components with efficiency and precision. One of the significant advantages of WAAM is its remarkable deposition rate, which can range from 2 to 10 kg per hour [2]. This rate far exceeds some other additive manufacturing techniques, such as Laser Powder Bed Fusion (LPBF). Due to its potential to create sizeable objects while maintaining reasonable production rates, WAAM has found applications in various industries, including aerospace, automotive, and marine [3,4].

However, residual stresses and deformations are inherent challenges in WAAM, stemming from the complex thermal and mechanical interactions that occur during the fabrication of metallic components [5,6]. As molten material solidifies and cools, internal stresses develop due to uneven temperature gradients. Additionally, the layer-by-layer deposition characteristic of WAAM introduces unique mechanical interactions that influence the final shape and structural integrity of the manufactured object. Understanding and managing residual stresses and deformations is crucial, as they can adversely affect the dimensional accuracy, mechanical properties, and overall performance of WAAM produced parts [7]. Uncontrolled stresses can lead to warping, distortion, and even cracking of the fabricated components [8].

Simulation plays an important role to address these issues, because it is impractical to optimise the residual stresses and deformations merely by experimental trial and error [9–11]. The part-scale residual stresses and deformations are typically obtained through coupled thermo-mechanical finite element (FE) analysis [12–14]. The transient thermal analysis is performed first to acquire the temperature history, and then the thermal strain that is related to the thermal history is calculated and employed as the thermal load for the subsequent mechanical analysis. The double-ellipsoidal Goldak heat source, which was first developed for welding simulation, is the most commonly used heat source model in the thermal calculation of WAAM [15,16]. However, other types of heat sources, including point [17,18], line [19], Gaussian distributed surface [20] and volumetric heat sources [21–23] were also employed in some numerical models, and the results suggested it may achieve similar accuracy by choosing appropriate heat source parameters.

Lu et al. [24,25] conducted a thermo-mechanical analysis for the metal AM process, encompassing nonlinear effects such as the solid phase transformation and creep behaviour at high temperatures. The results suggested the cooling period contributes more to the residual stresses than the heating period. Weisz et al. [26] simulated the building process for a thin-wall structure with a height of 100 mm. This study suggested solid phase transformation and creep behaviour may be neglected without significantly affecting predicted stresses and deformations for 316L stainless steel. Sahoo et al [27,28] studied influence of hatch spacing through simulation and indicated that increasing hatch spacing may decrease temperature and the residual stresses.

However, a primary challenge in simulation lies in obtaining the thermo-mechanical response for the entire fabricated part within a reasonable computational cost. The main issue involves resolving the temporal and spatial mismatch between the local melting/cooling behaviour and the overall fabricated structure [9]. Closed-form solutions for a moving heat source on a semi-infinite space are commonly employed in various studies [29]. To address boundary conditions (BCs), Li et al. [30] utilised mirror heat sources, which can also be expressed analytically, while Ning et al. [31,32] characterised the BCs using heat sinks, for which the analytical thermal solution was developed. However, it appears these analytical solutions are applicable only to simple geometries.

Based on the analytical thermal solution, some semi-analytical approaches are developed. Daniel et al. [33] proposed a semi-analytical solution of the transient heat conduction problem for components with rotational features. A more general semi-analytical approach is proposed in Refs. [33–36], in which the real temperature is divided into two distinct fields, i.e. the analytical field, which provides a closed-form solution to describe the temperature evolution without considering BCs, and the complimentary field, which accounts for the BCs and is solved numerically. For solving the complimentary field, there is no need to explicitly considering the heat source term because of the assistance of the analytical field. Consequently, the mismatch in spatial domain between the heat source and the overall part is alleviated in solving the complementary field, and thus it allows for the employment of a coarser mesh that scales with the dimension of the fabricated part.

While the semi-analytical approach can mitigate the spatial mismatch problem, the issue of temporal mismatch persists. Therefore, in the present study, based on the semi-analytical model we developed previously [34,35], it is attempted to solve the temporal mismatch issue by proposing a temporal acceleration strategy. The outline of the paper is as follows. Section 2 introduces the semi-analytical thermo-mechanical model and the temporal acceleration strategy. In Section 3, to assess the impact of the proposed acceleration strategy on temperature, stresses, and deformations, a numerical example regarding a simple thin-wall structure is conducted. The obtained results are compared with corresponding experimental data in the literature in Section 4. Following validation of the accuracy of the temporal acceleration strategy, a simulation and corresponding experiment conducted in-house for fabricating a practical engineering WAAM component are reported in Section 5. The article concludes with a reiteration of the most salient points of the study in Section 6.

2. Model description

The semi-analytical thermal model and subsequent mechanical model are first explained in this section. Then the proposed temporal acceleration strategy is illustrated. The numerical implementations are also briefly introduced at the end of this section.

2.1. Thermal model

In the thermal model, based on the Fourier’s law, the governing equation is expressed as: (1) $\rho c_p{\displaystyle \frac{\mathrm{\partial }T}{\mathrm{\partial }t}=\mathrm{\nabla }\cdot \left(k\mathrm{\nabla T}\right)+{\dot{Q}}_v,}$(1) where the temperature is denoted as T, and ${\dot{Q}}_v$ is the rate of volumetric heat generation, i.e. the source term. The thermal parameters ρ, c_p and k are the density, the constant-pressure specific heat and the thermal conductivity, respectively.

The semi-analytical thermal approach is based on the assumption that considerable accuracy can be ensured by choosing appropriate temperature independent thermal properties [34,37,38]. Additionally, it assumes the solid phase transformations and phase transitions between the solid, liquid and vapour can be neglected, as they are considered to cause only a second order effect on stresses and deformations [34,35]. With these assumptions, Equation (1)(1) $\rho c_p{\displaystyle \frac{\mathrm{\partial }T}{\mathrm{\partial }t}=\mathrm{\nabla }\cdot \left(k\mathrm{\nabla T}\right)+{\dot{Q}}_v,}$(1) can be written as a linearised form: (2) ${\displaystyle \frac{\mathrm{\partial }T}{\mathrm{\partial }t}=\alpha {\mathrm{\nabla }}^{2}T+{\dot{Q}}_v,}$(2) where $\alpha =k/\rho c_p$ is the thermal diffusivity.

Then considering the temperature T to be divided into two parts, i.e. (3) $T=\tilde{T}+\hat{T}$(3) where $\tilde{T}$ is the analytical temperature field, and $\hat{T}$ is the complimentary field. Substituting Equation (3)(3) $T=\tilde{T}+\hat{T}$(3) into Equation (2)(2) ${\displaystyle \frac{\mathrm{\partial }T}{\mathrm{\partial }t}=\alpha {\mathrm{\nabla }}^{2}T+{\dot{Q}}_v,}$(2) , we obtain: (4) ${\displaystyle \frac{\mathrm{\partial }\tilde{T}}{\mathrm{\partial }t}-\alpha {\mathrm{\nabla }}^2\tilde{T}-{\displaystyle \frac{{\dot{Q}}_v}{\rho c_p}={\displaystyle \frac{\mathrm{\partial }\hat{T}}{\mathrm{\partial }t}-\alpha {\mathrm{\nabla }}^2\hat{T}}}}$(4) If an analytical solution exists to let the left part of Equation (4)(4) ${\displaystyle \frac{\mathrm{\partial }\tilde{T}}{\mathrm{\partial }t}-\alpha {\mathrm{\nabla }}^2\tilde{T}-{\displaystyle \frac{{\dot{Q}}_v}{\rho c_p}={\displaystyle \frac{\mathrm{\partial }\hat{T}}{\mathrm{\partial }t}-\alpha {\mathrm{\nabla }}^2\hat{T}}}}$(4) be zero, then only the following equation needs to be solved: (5) ${\displaystyle \frac{\mathrm{\partial }\hat{T}}{\mathrm{\partial }t}=\alpha {\mathrm{\nabla }}^2\hat{T}}$(5) Upon comparing Equations (5)(5) ${\displaystyle \frac{\mathrm{\partial }\hat{T}}{\mathrm{\partial }t}=\alpha {\mathrm{\nabla }}^2\hat{T}}$(5) to (2)(2) ${\displaystyle \frac{\mathrm{\partial }T}{\mathrm{\partial }t}=\alpha {\mathrm{\nabla }}^{2}T+{\dot{Q}}_v,}$(2) , the advantage is the heat source term ${\dot{Q}}_v$ disappears in Equation (5)(5) ${\displaystyle \frac{\mathrm{\partial }\hat{T}}{\mathrm{\partial }t}=\alpha {\mathrm{\nabla }}^2\hat{T}}$(5) . Consequently, there is no need to employ a considerably fine mesh to describe the heat source, thus allowing for the use of a coarser mesh that scales with the size of the component.

Now the problem is how to find a closed-form solution to let the left part of Equation (4)(4) ${\displaystyle \frac{\mathrm{\partial }\tilde{T}}{\mathrm{\partial }t}-\alpha {\mathrm{\nabla }}^2\tilde{T}-{\displaystyle \frac{{\dot{Q}}_v}{\rho c_p}={\displaystyle \frac{\mathrm{\partial }\hat{T}}{\mathrm{\partial }t}-\alpha {\mathrm{\nabla }}^2\hat{T}}}}$(4) be zero. For a moving heat source, it can be considered as the combination of a large number of instantaneous point heat sources. For any given point heat source I, if it is imposed on the boundary surface of a semi-infinite space, a closed-form solution of the temperature caused by the point heat source is written as [34]: (6a) ${\tilde{T}}^{\left(I\right)}\left(x_i,t\right)={\displaystyle \frac{P\mathrm{\Delta tA}}{4\rho c_p{\left(\mathrm{\pi \alpha }\left(t-t_0^{\left(I\right)}+t_r\right)\right)}^{3/2}}\exp \left({\displaystyle \frac{-R^2}{4\alpha \left(t-t_0^{\left(I\right)}+t_r\right)}}\right)}$(6a) where the superscript I means the information for the I^th point heat source. The $t_0^{\left(I\right)}$ is the activated moment of the I^th point heat source and P is the corresponding power. The Δt is the duration of each point heat source. The absorptivity is denoted as A, and R is the distance between the location of the heat source $x_i^{\left(I\right)}$ and the point of interest x_i, which is expressed as: (6b) $R^2=(x_1-x_1^{\left(I\right)}{)}^2+(x_2-x_2^{\left(I\right)}{)}^2+(x_3-x_3^{\left(I\right)}{)}^2$(6b) The parameter t_r is called the diffusion time and proposed to be equal to $r_l^2/8\alpha$, where $r_l$ is the characteristic size of the heat source in the real situation. The diffusion time is on one hand to avoid singularity when t = t₀, and on the other hand to characterise the size effect of the heat source in the real situation. It is considered that the material is heated after the heat energy has diffused a certain distance from the centre of the heat source, and $t_r=r_l^2/8\alpha$ corresponds to the time that the heat energy has diffused a distance of $r_l$. This has been proved to be an accurate characterisation of the thermal response in metal additive manufacturing [34].

Remind Equation (6a)(6a) ${\tilde{T}}^{\left(I\right)}\left(x_i,t\right)={\displaystyle \frac{P\mathrm{\Delta tA}}{4\rho c_p{\left(\mathrm{\pi \alpha }\left(t-t_0^{\left(I\right)}+t_r\right)\right)}^{3/2}}\exp \left({\displaystyle \frac{-R^2}{4\alpha \left(t-t_0^{\left(I\right)}+t_r\right)}}\right)}$(6a) only holds true when Δt is small enough, which means each point heat source is instantaneous. When Δt approaches 0, the $\tilde{T}$ of the moving heat source is expressed as: (7) $\tilde{T}\left(x_i,t\right)=\int {\tilde{T}}^{\left(I\right)}\left(x_i,t\right)\mathrm{d}t,$(7) If the heat source moves straightly along a line, then the closed-form expression of Equation (7)(6b) $R^2=(x_1-x_1^{\left(I\right)}{)}^2+(x_2-x_2^{\left(I\right)}{)}^2+(x_3-x_3^{\left(I\right)}{)}^2$(6b) for this moving line heat source is derived in Ref. [36], and the specific expression is given in Appendix A.

Finally, after obtaining $\tilde{T}$, the $\hat{T}$ can be calculated with corresponding BCs, which are expressed as a function of $\tilde{T}$. According to the real conditions, both Dirichlet and Neumann conditions can be considered in solving $\hat{T}$, and the details are given in Appendix B. The relationship between the real temperature T, the analytical temperature $\tilde{T}$ and the complementary temperature $\hat{T}$ can be illustrated in Figure 1.

Figure 1. Illustration of the semi-analytical thermal approach.

2.2. Mechanical model

Once the temperature T is obtained, the thermal and mechanical models are connected through the thermal strain ${ϵ}_{\mathrm{ij}}^{\mathrm{th}}$, which is calculated as: (8) ${ϵ}_{\mathrm{ij}}^{\mathrm{th}}={\delta }_{\mathrm{ij}}{\int }_{T_{\mathrm{ref}}}^T{\text{β}}_{\mathrm{th}}\left(T\right)\mathrm{d}T$(8) where δ_ij is the Kronecker delta and β_th is the coefficient of thermal expansion. The reference temperature T_ref is usually set as the room temperature. The total strain ${ϵ}_{\mathrm{ij}}$ is expressed as: (9) ${ϵ}_{\mathrm{ij}}={ϵ}_{\mathrm{ij}}^e+{ϵ}_{\mathrm{ij}}^p+{ϵ}_{\mathrm{ij}}^{\mathrm{th}}$(9) where ${ϵ}_{\mathrm{ij}}^e$ and ${ϵ}_{\mathrm{ij}}^p$ are the elastic and the plastic strains, respectively.

The constitutive law in an elasto-plastic analysis is written as: (10) $\mathrm{d}{\sigma }_{\mathrm{ij}}=D_{\mathrm{ijkl}}^{\mathrm{ep}}(\mathrm{d}{ϵ}_{\mathrm{kl}}^e+{ϵ}_{\mathrm{kl}}^p)=D_{\mathrm{ijkl}}^{\mathrm{ep}}(\mathrm{d}{ϵ}_{\mathrm{kl}}-\mathrm{d}{ϵ}_{\mathrm{kl}}^{\mathrm{th}})$(10) where ${\sigma }_{\mathrm{ij}}$ is the stress, and $D_{\mathrm{ijkl}}^{\mathrm{ep}}$ is the elasto-plastic constitutive matrix which is expressed as: (11) $D_{\mathrm{ijkl}}^{\mathrm{ep}}=D_{\mathrm{ijkl}}^e\left(1-{\displaystyle \frac{S_{\mathrm{ij}}{S}_{\mathrm{ij}}{D}_{\mathrm{ijkl}}^e}{{\displaystyle \frac{4H_i{\sigma }_Y^2}{9}+S_{\mathrm{ij}}{D}_{\mathrm{ijkl}}^e{S}_{\mathrm{kl}}}}}\right)$(11) The $D_{\mathrm{ijkl}}^e$ is the hook elastic tensor, ${\sigma }_Y$ is the yield stress, S_ij is the deviatoric stress, and H_i is the plastic tangent modulus. The von Mises yield criterion with isotropic hardening is employed.

Small deformation is assumed, and thus the relationship between strain and displacement is written as: (12) ${ϵ}_{\mathrm{ij}}={\displaystyle \frac{1}{2}\left({\displaystyle \frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }{x}_j}+{\displaystyle \frac{\mathrm{\partial }{u}_j}{\mathrm{\partial }{x}_i}}}\right)}$(12) where x_i (i = 1, 2, 3) represents the coordinates, and u_i is the corresponding displacement component.

For every material point, the force equilibrium equation is expressed as: (13) ${\displaystyle \frac{\mathrm{\partial }{\sigma }_{\mathrm{ij}}}{\mathrm{\partial }{x}_j}+b_i=0}$(13) where b_i is the body force components along the x_i direction. The body force, such as gravity, is commonly neglected [39–41], and hence the influence of body force is not considered in the present study. The stress/strain and displacement for every material point can be finally obtained by solving Equations (8)(7) $\tilde{T}\left(x_i,t\right)=\int {\tilde{T}}^{\left(I\right)}\left(x_i,t\right)\mathrm{d}t,$(7) –(13)(12) ${ϵ}_{\mathrm{ij}}={\displaystyle \frac{1}{2}\left({\displaystyle \frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }{x}_j}+{\displaystyle \frac{\mathrm{\partial }{u}_j}{\mathrm{\partial }{x}_i}}}\right)}$(12) .

2.3. Acceleration strategy in temporal domain

The fabrication process in WAAM can be divided into the heating period, for which the heat source is imposed on the part, and the cooling period, for which no heat source is imposed and the temperature of the entire component decreases. During the heating period in the simulation, the velocity of the moving heat source is accelerated by a factor of κ. This means the moving speed of the heat source in the simulation v_sim equals to κv_real, where v_real is the real speed of the heat source in the experiment. To maintain the same power density P/v, the power employed in the simulation P_sim is also increased to be P_sim = κP_real, where P_real is the power employed in the experiment.

By the above treatment, the time for the heating period will be decreased with a factor of κ. For the cooling period, the real cooling time is adopted without any acceleration. The main reason for this is to ensure the simulated part can still be cooled down to a reasonable temperature. Besides, the temperature gradient for the cooling period is normally much smoother than that of the heating period, thus larger time increment can be applied during the cooling period to guarantee computational efficiency and accuracy.

It also can be found later in Section 4 that the predicted temperature will be higher with more focused heat energy distribution. To compensate for this effect, it is proposed to manipulate the diffusion time t_r shown in Equation (6a)(6a) ${\tilde{T}}^{\left(I\right)}\left(x_i,t\right)={\displaystyle \frac{P\mathrm{\Delta tA}}{4\rho c_p{\left(\mathrm{\pi \alpha }\left(t-t_0^{\left(I\right)}+t_r\right)\right)}^{3/2}}\exp \left({\displaystyle \frac{-R^2}{4\alpha \left(t-t_0^{\left(I\right)}+t_r\right)}}\right)}$(6a) . Remind t_r is the diffusion time corresponding to the time that the heat energy has diffused from the centre of the heat source to the point of interest. Hence, if larger t_r is employed, the temperature distribution will be more evenly distributed with reduced peak value, and it is possible to obtain similar temperature profiles to those without temporal acceleration. More discussions regarding this aspect can be found in Section 4.

2.4. Numerical implementations

The proposed model is realised using the commercial FE code Abaqus with several user subroutines. The subroutine umatht is utilised to calculate the $\tilde{T}$ and $\hat{T}$, and the subroutines dflux and disp are employed to impose the thermal BCs. The subroutine uexpan is used to calculate the thermal strain from the temperature field, and the subroutine uexternaldb is employed to read the moving path of the heat source.

3. Problem description

A WAAM experiment conducted by Ding et al. [42] is employed for the purpose of validation in this section. A four-layer structure as shown in Figure 2(a) was fabricated in Ref. [42]. Each layer was built with only a single track. The power was 2245 W, and the velocity of the heat source was 8.33 mm/s. The cooling time after building each layer was 400 s. In Figure 2(a), the temperature history of point G1 and G2, and the temperature distribution along path 1 are presented later in Section 4 to investigate the proposed temporal acceleration strategy.

Figure 2. (a) The geometry of the simulated structure [42]. (b) Half of the geometry modelled by taking advantage of the symmetry of the problem.

Previous studies [36,43] have shown that a thermo-mechanical model based on the semi-analytical approach without adopting any acceleration in the temporal domain can achieve considerable accuracy in prediction of the thermal and mechanical response for this experiment in Ref. [42]. Hence, the same thermal and mechanical properties, mesh sizes and BCs are employed in the present paper. The material is mild steel S355, and the thermal and mechanical properties are listed in Appendix C. The mesh size for the built part is 5 × 0.3125 × 2 mm, and 5 × 5 × 2 mm for the baseplate.

The mechanical BCs are shown in Figure 2(a), i.e. the baseplate is fastened to a large backing plate by clamping it at six specific locations. By taking advantage of the symmetry of the problem, only half of the geometry is modelled, as shown in Figure 2(b). The backing plate is modelled as a rigid body and is in contact with the baseplate. In the thermal model, Neumann BCs are applied for all the surfaces. A large convection coefficient of 300 W/m²K is applied on the bottom surface of the baseplate to represent the effect of the cooling system placed below that surface, and for the other boundary surfaces the convection coefficient is 5.7 W/m²K.

As suggested in Ref. [42], a cut-off temperature of 1000°C can be applied in the simulation to avoid difficulties in numerical convergence. The cut-off temperature means any variations of mechanical properties that are above this cut-off temperature are disregarded, and thus the cut-off temperature serves as the upper limit for temperature considerations in the mechanical analysis.

4. Results and discussions

The predicted temperature profile is first investigated in Figure 3. In Figure 3, the reference means the temperature profile that is calculated without acceleration on the temporal domain, and the original diffusion time $t_r=r_l^2/8\alpha =0.01\mathrm{s}$ is applied. Figure 3(a) plots the analytical temperature $\tilde{T}$ along the x₂ axis when the moving heat source arrives at point P. At this stage, no BCs is involved, i.e. it is the condition that the moving heat source is imposed on the boundary surface of a semi-analytical space. This describes the situation where the heat source is located far from the boundary. The red solid line is the reference, and the other solid lines are the values calculated with various acceleration factor κ. The original diffusion time $t_r=r_l^2/8\alpha =0.01\mathrm{s}$ is applied for all the calculations in Figure 3(a).

Figure 3. (a, b) The analytical temperature $\tilde{T}$ distribution along x₂ axis. In (a), various acceleration factor κ with the original diffusion time t_r = 0.01 s are employed, and in (b), the acceleration factor κ is set as 200 with various diffusion time t_r. (c) The real temperature T along path 1 (see Figure 2(a)).

It can be found that with increasing acceleration factor κ, the peak temperature also increases with narrower profile, which means a more concentrated energy distribution around the heat source. However, the temperature profile tends to converge when the κ is greater than 50. This is because when the acceleration factor κ is greater than a certain value, it means the heat source moves almost instantaneously from one location to another. It is also of interest to note that the peak temperature is unrealistically high which is due to the neglect of the phase transitions. It should be noted that because all the mechanical properties are set as the same when it is above the cut-off temperature of 1000°C, the temperature profiles that are below the cut-off temperature mainly affect the subsequent mechanical analysis.

In Figure 3(b), for the acceleration factor of κ  = 200, various values of t_r are tested to achieve a similar temperature profile to the reference. It can be observed that with increasing t_r, the temperature is more evenly distributed with reduced peak value. As indicated before, only the temperature below the cut-off temperature (indicated as the blue dashed line in Figure 3(b)) is vital for the subsequent mechanical calculation, it is found when t_r = 0.22 s, a similar temperature profile to the reference for the values that are below the cut-off temperature can be obtained.

When the BCs are considered, the real temperature profiles are plotted in Figure 3(c). The temperatures are obtained along path 1 (see Figure 2(a)) at the moment when the temperature at the centre of the path (x₂ = 0) reaches the maximum. Because a large convection coefficient of 300 W/m²K is applied on the bottom surface of the baseplate to represent the effect of the cooling system placed beneath that surface, it can be observed that the peak temperature in Figure 3(c) is lower than that in Figure 3(b). In Figure 3(c), it can be seen that when the acceleration factor κ is 200 without any manipulation of t_r, more concentrated temperature profile is obtained. In comparison, when t_r is increased to be 0.22 s, the corresponding temperature profile agrees well with the reference, especially for the values that are below 1500°C.

Figure 4 further shows investigation results regarding temperature history with the acceleration factor κ = 200 and t_r = 0.22 s. The temperature evolutions of points G1 and G2 (see Figure 2(a)) are plotted. The blue dashed lines are obtained from the experimental measurements [42], and the black dashed lines are from the simulation. Since in the simulation the heating period is accelerated, the corresponding time range of the heating period is different from that in the experiment. Hence, only the temperature evolutions during the cooling period are compared, and the correspondingly starting moment for the simulated cooling period is also shifted to match that of the experiment. It also should be noted that the stresses and deformations are more closely related to the cooling rates compared to the heating rates [24,44]. It can be observed that for point G1 (see Figure 4(a)) the predicted peak temperatures are slightly higher than that of the experiments (∼50°C), which is due to the acceleration strategy applied. This overestimation diminishes rapidly and the predicted peak temperatures closely match the experimental data when it is for point G2, which is only 20 mm away from the moving heat source line. With the accurately predicted peak temperature, the simulated cooling curve agrees well with the experimental measurements as the cooling stage is not affected by the acceleration factor κ in the simulation.

Figure 4. The temperature as a function of time for point (a) G1 and (b) G2 (see Figure 2(a)). For the simulated results, only the cooling period is displayed, and the correspondingly starting moment of the simulated cooling period is shifted to match that of the experiment.

It should be noted that the acceleration factor κ cannot be increased without limitation. This is because as the acceleration factor κ increases, the heating period decreases, resulting in a corresponding decrease in the time increment. A Convergence problem arises when the time increment becomes too small, and it is found the simulation becomes difficult to converge when the acceleration factor κ is greater than 200. Therefore, the acceleration factor κ is set as 200 for the subsequent mechanical calculation.

With accurate estimation of the temperature history, the determined acceleration parameters (κ = 200, t_r = 0.22 s) are employed to further investigate the prediction on stresses and deformations. Figure 5 displays the von Mises stresses ${\sigma }_{\mathrm{VM}}$ calculated along path 2 which is at the interface between the built part and baseplate. Figure 5(a) corresponds to the moment when the fourth layer is built and the baseplate is still clamped on the backing plate, while Figure 5(b) corresponds to the moment that the clamps on the baseplate are removed. The red solid lines are calculated from the model without temporal acceleration and manipulation of t_r, and hence taken as the reference. The black lines are calculated with the acceleration parameters (κ = 200, t_r = 0.22 s).

Figure 5. The von Mises stress along path 2 (a) before and (b) after removal from the backing plate.

For the results of the reference, the maximum ${\sigma }_{\mathrm{VM}}$ is observed at both ends, while at other locations, the predicted ${\sigma }_{\mathrm{VM}}$ is ∼450 MPa before the baseplate is removed from the backing plate, and this value decreases to ∼350 MPa after removal from the backing plate. The entire distribution of ${\sigma }_{\mathrm{VM}}$ obtained based on the acceleration strategy is in line with the reference results, although the predicted values are ∼50 MPa lower than those in the reference findings. This small discrepancy is considered as the consequence of the slight mismatch of the predicted temperature with the reference as shown in Figures 3(c) and 4(a).

Figure 6 shows the deformation along the x₃ direction of the bottom surface of the baseplate. The blue squares are obtained from the experiment [42], while the black solid and red dashed lines are the simulated results with and without temporal acceleration, respectively. The results in Figure 6 correspond to the moment that the clamps on the baseplate are removed after the fourth layer is built. It can be seen that the predicted deformations with and without temporal acceleration both closely match with the experimental measurements. This implies that temporal acceleration has little influence on the accuracy of the simulated deformations, and based on the predicted temperature, stresses and deformations from Figures 4 to 6, it provides confidence in the proposed temporal acceleration strategy.

Figure 6. The displacement component U₃ after building the fourth layer and removing the clamps.

Regarding the computational efficiency, all the calculations are performed on Abaqus using four CPUs (13th Gen Intel Core i7-1365U) with a processor speed of 1.80 GHz and a memory of 32 GB RAM. The total CPU time for the reference model without temporal acceleration is 10,551 s, and the total CPU time for the model with the acceleration factor κ = 200 is 231 s. It can be seen that by implementing the acceleration strategy, it exhibits attractive computational efficiency without sacrificing much accuracy.

5. Simulation for a practical engineering component

In this section, to further validate the effectiveness of the proposed acceleration strategy, the WAAM process for a more complicated part is simulated. The dimension of the part is illustrated in Figure 7(a), while Figure 7(b) shows the corresponding experimental setup conducted in-house. The simulated deformation of the fabricated part is compared with the experimental measurements to assess the accuracy of the numerical model.

Figure 7. (a) Geometry of the WAAM component. (b) Experiment setup of the WAAM process.

There are 11 layers for the part shown in Figure 7(a), and the layer thickness is 2.5 mm. These 11 layers correspond to a full plate of metal wire, for which the weight is ∼50 kg. The metal wire was melted using the plasma arc. The material employed was TC11, and the thermal and mechanical properties are listed in Appendix C. This fabricated part belongs to a practical engineering component used in rocket bundling.

The power of the heat source was 4200 W and the moving speed was 3 mm/s. The moving path of the heat source for a certain layer is shown in Figure 8. It can be seen that multiple moving tracks were applied for each layer. Upon completion of each track, the heat source moved to the starting location of the subsequent track, with the moving speed of 20 mm/s. Consequently, the cooling time between adjacent tracks varied depending on the location of each track. Files containing this information about the moving heat source can be generated from the WAAM machine. In the simulation, these files are read by the numerical model to determine the location of the heat source at any given time.

Figure 8. Moving path of the heat source for a certain layer.

The total real fabrication time for the 11 layers were ∼63.6 h, which is difficult to model without using a temporal acceleration. Hence, the proposed temporal acceleration strategy is employed in the simulation, and the acceleration factor κ for the heating period is set as 200. The same approach as indicated in Figure 3(b) is employed to determine the appropriate diffusion time t_r, and the t_r is finally determined as 0.5 s. The initial temperature is set as 25°C. The characteristic mesh size is ∼30–50 mm, and the results show little difference when further decreasing the mesh size. For the thermal BCs, the Dirichlet BC is applied for bottom surface of the baseplate with a constant temperature of 25°C during the whole process, and Neumann BCs are applied for the other boundary surfaces with the convection coefficient of 5.7 W/m²K. The mechanical BCs are shown in Figure 7(a). The bottom surface of the baseplate is in contact with the backing plate, and the displacements of the four clamps are fixed. The backing plate is again modelled as a rigid body. The simulation is performed on the same conditions as in Section 3, and the total CPU time for fabricating the 11 layers is ∼8 h.

In the experiment, once the 11 layers were fabricated, the clamps were removed and the part along with the baseplate were subsequently placed in an oven to conduct an annealing process aimed at stress relief. The oven was initially heated to 850°C at a rate of 200°C/h, and then maintained at 850°C for 5 h. Finally, the oven was allowed to cool down to the room temperature over a period of 11 h. This annealing process is also simulated in this section. As the material during the annealing process is mainly governed by the creep behaviour, the creep behaviour of the material is described as: (14) ${\dot{ϵ}}_{\mathrm{eq}}^{\mathrm{creep}}=\left\{B{\sigma }_{\mathrm{VM}}^{\lambda }\right[(m+1){ϵ}_{\mathrm{eq}}^{\mathrm{creep}}{]}^m{\}}^{1/(m+1)}$(14) where ${ϵ}_{\mathrm{eq}}^{\mathrm{creep}}$ and ${\dot{ϵ}}_{\mathrm{eq}}^{\mathrm{creep}}$ are the equivalent creep strain and equivalent creep strain rate. The B, λ and m are the material constants. The creep data of TC11 in Ref. [45] is used and how the material constants are determined is explained in Appendix C.

The simulated deformation along the x₃ direction, i.e. the building direction, of the fabricated component before removal from the backing plate is depicted in Figure 9. The distortion at the left end is 8.8 mm, while it is 14.4 mm at the right end. This is close to the experimental observations of 8.0 mm at the left end and 13.2 mm at the right end. This good agreement again verifies the proposed temporal acceleration strategy is capable of achieving satisfied accuracy.

Figure 9. The deformation after the part is fabricated. A scale factor of 5 is applied.

Figure 10 compares von Mises stress of the component before (Figure 10(a)) and after (Figure 10(b)) the annealing process. The component is removed from the backing plate for both the situations in Figure 10(a) and (b). Before annealing, the maximum stresses can be observed at the interface between the part and the baseplate, while after the annealing process, the stresses of the component considerably decrease to a level that is less than 50 MPa (see Figure 10(b)). This reduction is attributed to the relaxation of the elastic strain components induced by the creep strain.

Figure 10. The von Mises stresses of the component (a) before and after (b) the annealing process. The clamps are removed for both situations.

Figure 11 shows the simulated deformation along the building direction before and after the annealing process along path 3 (inset of Figure 11). The displacements along the building direction at right and left ends are set as 0. Before the annealing process, when the component is removed from the backing plate, the warpage at the right and left ends comparing to the middle of the part is of ∼10 mm, while after the annealing process, the warpage decreases to ∼4 mm. For the experiment, a warpage of 4.2 mm at the right and left ends after annealing was observed. Hence, the numerical predictions agree well with the experimental measurement, which results from the accurate calculations of the stress/strain states of the WAAM process.

Figure 11. The displacement component U₃ along path 3 before and after the annealing. The clamps are removed for both situations.

6. Conclusions

The present study proposed an acceleration strategy on the temporal domain based on the semi-analytical approach for computationally efficient prediction of stresses and deformations of the WAAM component. With the acceleration of the moving speed for the heat source, the applied power is also increased to maintain the same power density. This treatment results in more concentrated heat energy distribution around the heat source. To compensate for this effect, the diffusion time in the semi-analytical approach is manipulated to achieve similar temperature profiles to those without temporal acceleration.

By comparing with the experiments for a simple thin-wall structure, it shows the predicted stresses are only 50 MPa lower than the experimental values, while the predicted deformations agree considerably well with those of the experiment. For a more complicated large WAAM component, the discrepancy between the predicted deformations and the experimental measurements is less than 10%. Overall, the proposed temporal acceleration strategy offers attractive computational efficiency without sacrificing much accuracy. This enables the modelling of the WAAM process for a practical engineering component. The moving path of the heat source can explicitly be considered in the calculation, and the subsequent post-heat treatment can also be integrated in the model.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability

The raw/processed data required to reproduce these findings is available from the corresponding authors on reasonable request.

Additional information

Funding

Y. Yang is sponsored by National Natural Science Foundation of China (Grant Number 52105421) and Guangdong Provincial Natural Science Foundation (Grant Number 2022A1515011621). Q. Li is sponsored by Beijing Nova Program (Grant Number 20230484294) and Young Elite Scientist Sponsorship Program by CAST (Grant Number 2019QNRC001).

Appendix A.

Analytical solution of the moving line heat source

The closed-form solution of the analytical field $\tilde{T}$ for a moving line heat source is expressed as: (A1) $\begin{array}{rl}& \tilde{T}\left(x_i,t\right)={\displaystyle \frac{\mathrm{PA}}{4\mathrm{k\pi }\sqrt{Y}}\exp \left(-{\displaystyle \frac{\left(\sqrt{Y}+C\right)\nu }{2\alpha }}\right)}\\ & \left[\mathrm{erfc}\left(z_1\right)-z_5\mathrm{erfc}\left(z_2\right)-\mathrm{erfc}\left(z_3\right)+z_5\mathrm{erfc}\left(z_4\right)\right]\end{array}$(A1) where (A2) $\begin{array}{rl}& Y=R_s^2-2\left(x_1-x_1^s\right)\nu t^{\ast }\cos \theta -2\left(x_2-x_2^s\right)\nu t^{\ast }\sin \theta +\nu t^{\ast }{}^2,\\ & C=\left(x_1-x_1^s\right)\cos \theta +\left(x_2-x_2^s\right)\sin \theta -\nu t^{\ast },\\ & z_1={\displaystyle \frac{\nu }{2\sqrt{\alpha }U}-{\displaystyle \frac{\sqrt{Y}U}{2\sqrt{\alpha }},}}\\ & z_2={\displaystyle \frac{\nu }{2\sqrt{\alpha }U}+{\displaystyle \frac{\sqrt{Y}U}{2\sqrt{\alpha }},}}\\ & z_3={\displaystyle \frac{\nu }{2\sqrt{\alpha }L}-{\displaystyle \frac{\sqrt{Y}L}{2\sqrt{\alpha }},}}\\ & z_4={\displaystyle \frac{\nu }{2\sqrt{\alpha }L}+{\displaystyle \frac{\sqrt{Y}L}{2\sqrt{\alpha }},}}\\ & z_5=\exp \left({\displaystyle \frac{\sqrt{Y}\nu }{\alpha }}\right),\\ & L=(t^{\ast }+t_r{)}^{-1/2},\\ & U=\{\begin{array}{c}{t}_r^{-1/2},\\ {\left(t-t_f+t_r\right)}^{-1/2},\end{array}\begin{array}{c}t\le t_f.\\ t>t_f.\end{array}\end{array}$(A2) The moving speed of the heat source is denoted as v. The termination instance of the moving line heat source is denoted by t_f. The t* = t − t₀, where t₀ is the starting moment of the moving line heat source. The θ represents the angle of the moving line heat source with respect to the x₁ axis. The R_s is the distance between the point of interest x_i to the origin of the moving line heat source x_i^s, i.e. (A3) $R_s^2=(x_1-x_1^s{)}^2+(x_2-x_2^s{)}^2+(x_3-x_3^s{)}^2$(A3) The error function erfc(φ) is defined as: (A4) $\mathrm{erfc}\left(\phi \right)={\displaystyle \frac{2}{\sqrt{\pi }}{\int }_{\varphi }^{\mathrm{\infty }}\exp \left(-{\xi }^2\right)\mathrm{d}\xi }$(A4)

Appendix B.

The BCs for solving the complimentary field

The BCs in thermal problem can normally expressed as: (B1) $T=\overline{T},\mathrm{Dirichlet}\mathrm{BC}$(B1) (B2) $\mathrm{\nabla T}\cdot \mathit{n}={\overline{q}}_f,\mathrm{Neumann}\mathrm{BC}$(B2) where $\overline{T}$ and ${\overline{q}}_f$ are the prescribed temperature and heat flux, respectively, and n is the outward facing normal of the boundary surface. According the superposition of the temperature shown in Equation (3)(3) $T=\tilde{T}+\hat{T}$(3) , Equations (B1) and (B2) yield: (B3) $\hat{T}=\overline{T}-\tilde{T},\mathrm{Dirichlet}\mathrm{BC}$(B3) (B4) $\mathrm{\nabla }\hat{T}\cdot \mathit{n}={\overline{q}}_f-\mathrm{\nabla }\tilde{T}\cdot \mathit{n},\mathrm{Neumann}\mathrm{BC}$(B4) The expression of $\mathrm{\nabla }\tilde{T}$ is given as: (B5) $\begin{array}{rl}& {\displaystyle \frac{\mathrm{\partial }\tilde{T}}{\mathrm{\partial }{x}_1}=W_x\left(-\tilde{T}{W}_1\left(W_1+{\displaystyle \frac{v}{2\alpha }}\right)+W_0{W}_2{W}_1^2{W}_T\right)+\tilde{T}\left(-{\displaystyle \frac{\mathrm{vcos}\theta }{2\alpha }}\right),}\\ & {\displaystyle \frac{\mathrm{\partial }\tilde{T}}{\mathrm{\partial }{x}_2}=W_y\left(-\tilde{T}{W}_1\left(W_1+{\displaystyle \frac{v}{2\alpha }}\right)+W_0{W}_2{W}_1^2{W}_T\right)+\tilde{T}\left(-{\displaystyle \frac{\mathrm{vcos}\theta }{2\alpha }}\right),}\\ & {\displaystyle \frac{\mathrm{\partial }\tilde{T}}{\mathrm{\partial }{x}_3}=W_z\left(-\tilde{T}{W}_1\left(W_1+{\displaystyle \frac{v}{2\alpha }}\right)+W_0{W}_2{W}_1^2{W}_T\right)}\end{array}$(B5) where (B6) $\begin{array}{rl}& W_x=x_1-x_1^s-\nu t^{\ast }\cos \left(\theta \right),\\ & W_y=x_2-x_2^s-\nu t^{\ast }\sin \left(\theta \right),\\ & W_z=x_1-x_3^s,\\ & W_0={\displaystyle \frac{\mathrm{PA}}{4\mathrm{k\pi }},}\\ & W_1={\displaystyle \frac{1}{\sqrt{Y}},}\\ & W_2=\exp \left(-{\displaystyle \frac{\left(\sqrt{Y}+C\right)\nu }{2\alpha }}\right),\\ & W_T={\displaystyle \frac{U}{\sqrt{\mathrm{\pi \alpha }}}\left[\exp \left(-z_1^2\right)+\exp \left({\displaystyle \frac{\sqrt{Y}\nu }{\alpha }-z_2^2}\right)\right]}\\ & \begin{array}{cc}& \end{array}-{\displaystyle \frac{L}{\sqrt{\mathrm{\pi \alpha }}}\left[\exp \left(-z_3^2\right)+\exp \left({\displaystyle \frac{\sqrt{Y}\nu }{\alpha }-z_4^2}\right)\right]}\\ & \begin{array}{cc}& \end{array}+z_5\left[\mathrm{erfc}\left(z_4\right)-\mathrm{erfc}\left(z_2\right)\right]{\displaystyle \frac{\nu }{\alpha }.}\end{array}$(B6)

Appendix C.

Material properties

The thermal and mechanical properties of mild steel S355 are listed in Tables C1 and C2.

Table C1. Thermal properties of mild steel S355 employed in the simulation [27].

A (−)	ρcp	k
0.78	3.9 MJ/m^3K	48.6 W/mK

Table C2. Mechanical properties of mild steel S355 employed in the simulation.

T (°C)	E (GPa)	βth (°C^−1)	ν (−)	${\sigma }_Y$ $\left({ϵ}_{\mathrm{eq}}^P=0\right)$ Part (MPa)	${\sigma }_Y$ $\left({ϵ}_{\mathrm{eq}}^P=0.01\right)$ Part (MPa)	${\sigma }_Y$ $\left({ϵ}_{\mathrm{eq}}^P=0\right)$ Baseplate (MPa)	${\sigma }_Y$ $\left({ϵ}_{\mathrm{eq}}^P=0.01\right)$ Baseplate (MPa)
20	206	1.2 × 10^−5	0.29	450	520	350	420
100	203	1.2 × 10^−5	0.39	450	520	330	400
200	201	1.2 × 10^−5	0.29	420	500	305	380
300	200	1.2 × 10^−5	0.29	390	450	270	350
400	165	1.2 × 10^−5	0.30	320	370	230	290
500	100	1.2 × 10^−5	0.30	260	300	180	230
600	60	1.2 × 10^−5	0.32	170	215	125	160
700	40	1.2 × 10^−5	0.32	60	100	60	100
800	30	1.2 × 10^−5	0.35	50	80	60	100
900	20	1.2 × 10^−5	0.35	50	50	60	60
1000	10	1.5 × 10^−5	0.39	50	50	60	60

The E is the elastic modulus, υ is the Poisson’s ratio, and ${ϵ}_{\mathrm{eq}}^P$ is the equivalent plastic strain [32].

Table C3. Thermal properties of TC11 employed in the simulation.

A (−)	ρcp	k
0.8	3.8 MJ/m^3K	17.2 W/mK

Table C4. Mechanical properties of TC11 employed in the simulation.

T (°C)	E (GPa)	ν (−)	βth (°C^−1)	${\sigma }_Y$ (MPa) $\left({ϵ}_{\mathrm{eq}}^P=0\right)$	${\sigma }_Y$ (MPa) $\left({ϵ}_{\mathrm{eq}}^P=0.004\right)$
20	123.6	0.3	9.3 × 10^−6	1004	1004
200	105.6	0.3	9.3 × 10^−6	736.2	764.6
300	99	0.3	9.5 × 10^−6	676.7	719
400	94.3	0.3	9.7 × 10^−6	656.6	693.6
500	86.1	0.3	10 × 10^−6	628.4	664.4
600	73	0.3	10.2 × 10^−6	552.2	567.6

The thermal and mechanical properties of TC11 are listed in Tables C3 and C4. The creep data of TC11 is listed in Table C5, and the material constants for creep behaviour are listed in Table C6.

Table C5. Creep behaviour of TC11 [36].

T (°C)	${ϵ}_{\mathrm{eq}}^{\mathrm{creep}}$ (%)	${\dot{ϵ}}_{\mathrm{eq}}^{\mathrm{creep}}$ (×10^−7s^−1)	${\sigma }_{\mathrm{VM}}$ (MPa)
450	0.335	1.31	350 MPa
450	0.445	1.40	400 MPa
450	0.538	1.67	450 MPa
500	0.410	2.70	350
500	0.461	5.06	400
500	0.526	6.94	450
550	0.428	7.81	300
550	0.449	17.3	350
550	0.704	29.8	400

Table C6. Creep parameter employed in the simulation (in MPa-mm unit).

T (°C)	B	λ	m
450	5.20e−9	1.57	−0.66
500	1.06e−7	1.29	−0.89
550	1.61e−15	3.89	−0.26

According to Equation (14)(13) ${\displaystyle \frac{\mathrm{\partial }{\sigma }_{\mathrm{ij}}}{\mathrm{\partial }{x}_j}+b_i=0}$(13) , a logarithmic form of the equation can be obtained as: (C1) $\log \left(\dot{ϵ}\right)=f_1+f_2\log \left({\sigma }_{\mathrm{VM}}\right)+f_3\log \left(ϵ\right)$(C1) where (C2) $\begin{array}{rl}& f_1={\displaystyle \frac{\log \left(B\right)}{m+1}+{\displaystyle \frac{m}{m+1}\log (m+1)}}\\ & f_2={\displaystyle \frac{\lambda }{m+1}}\\ & f_3={\displaystyle \frac{m}{m+1}}\end{array}$(C2) It can be seen that the parameters of f₁, f₂ and f₃ are a function of B, λ, m. By substituting the data of Table C5 in Equation (C1)(C1) $\log \left(\dot{ϵ}\right)=f_1+f_2\log \left({\sigma }_{\mathrm{VM}}\right)+f_3\log \left(ϵ\right)$(C1) , the parameters of f₁, f₂ and f₃ can be determined, thereby enabling the determination of the creep constants B, λ, and m.