Efficiency and Accuracy in Thermal Simulation of Powder Bed Fusion of Bulk Metallic Glass

Additive manufacturing by powder bed fusion processes can be utilized to create bulk metallic glass as the process yields considerably high cooling rates. However, there is a risk that reheated material set in layers may become devitrified, i.e., crystallize. Therefore, it is advantageous to simulate the process to fully comprehend it and design it to avoid the aforementioned risk. However, a detailed simulation is computationally demanding. It is necessary to increase the computational speed while maintaining accuracy of the computed temperature field in critical regions. The current study evaluates a few approaches based on temporal reduction to achieve this. It is found that the evaluated approaches save a lot of time and accurately predict the temperature history.


INTRODUCTION
This article summarizes the applicability of approaches in computational welding mechanics (CWM) to simulate additive manufacturing (AM) and pinpoints the renewed necessity for performing efficient simulations to manage extremely large amounts of 'welds' in the case of additive manufacturing.
The study has a special emphasis on efficiency with retained accuracy for the powder bed fusion (PBF) of amorphous powder in creating bulk metallic glass (BMG). The rapid cooling of an appropriate multicomponent system provides it with an amorphous structure. This metallic glass can have special electric, magnetic, corrosion, and mechanical properties. 1 Powder bed fusion processes can be used to produce BMGs, as they give substantially high cooling rates of the solidifying material. 2,3 However, there is a risk that this material, when reheated at a temperature not sufficiently high to cause melting, may crystallize, i.e., become devitrified. Therefore, it is advantageous to simulate the process to fully comprehend it and design it to avoid the aforementioned risk. However, a detailed simulation is computationally demanding. Thus, simplifications enabling faster simulations are advantageous. The accuracy of these models in critical regions is important to be able to predict crystallization and devitrification. The studied amorphous alloy is AMZ4, Zr 59.3 Cu 28.8 Al 10.4 Nb 1.5 , and preliminary data are used to evaluate the modelling approaches.

COMPUTATIONAL WELDING MECHANICS
The overall aim of CWM 4 is to establish methods and models that are usable for the control and design of welding processes so as to obtain the appropriate mechanical performance of the welded component. It is focused on subjects ranging from the modelling of heat generation, weld pool phenomena, and heat flow to thermal stresses and deformations. Models in CWM can be combined with other models to obtain microstructure evolution and other features that enable the prediction of microstructure cracking and other phenomena. The centrepiece in modelling is the heat generation process, which is very complex for fusion welding processes. Some models may include fluid flow but exclude the physics of heat generation. However, the classical, i.e., common, approach is to also exclude this flow. The introduced heat input model must be calibrated with respect to experiments or obtained from more detailed weld process models. Initially, researchers applied various techniques for JOM https://doi.org/10.1007/s11837-018-2919-8 Ó 2018 The Author(s) reducing computing time, 4 whereas it is now possible to solve multipass welds in three dimensions (3D). However, the necessity for efficiency remains important when applying CWM to AM processes. 5,6 APPROACHES FOR EFFICIENT MACROSCALE SIMULATIONS OF PBF Approaches for reducing the wall-clock time when performing CWM simulations can be divided into three groups: spatial and temporal reductions and numerical techniques. The first two groups always bring some degree of approximation into the problem, whereas it is easier to control the degree of accuracy in the last group.
Spatial reduction explores the possibility of reducing a 3D problem to one that is two-dimensional. 7 This also implies a type of temporal reduction similar to the axisymmetric model of a circumferential pipe, and the use of a plane deformation model for a weld in a plate corresponds to infinite welding speed.
Temporal reduction ignores some of the history of the process. The use of analytic solutions, as well as solutions ignoring phase changes and similar simplifications, are examples of temporal reduction. One temporal reduction used in this study is the lumping of several welds into one, which can be combined with infinite welding speed. The latter means that, at least, parts of the weld are laid instantaneously. Ding et al. 8 based the stress calculation on peak stress during a cycle. This appears to be similar to envelope methods used early in CWM. Afazov et al. 6 used a combination of remeshing to group layers and an analytic temperature field for increased speed. The inherent strain approach replaces the weld process by an estimated final residual strain field that is supposed to yield the same residual stresses and deformations as the weld. It was previously published in 1970 9 and developed and promoted further by Ueda and coworkers. 10,11 This approach has recently been applied to AM and is available in commercial software.
The current study focuses on temporal reduction techniques for accelerating simulations of the PBF process (see Fig. 1). For each case shown, heating advances according to the arrows. The case in Fig. 1a has no temporal reduction, whereas the whole width of a string is heated instantaneously for the case in Fig. 1b, and the entire layer is heated for the case in Fig. 1c.

COMPUTATIONAL SETUP
Thermal analyses were performed using MSC.Marc with full Newton-Raphson iterations solving the resulting non-linear system of coupled equations. Adaptive time stepping was used to adjust the length of the time step to the nonlinearity of the solution during cooling. A cylinder with a diameter of 15 mm on a 10-mm-thick base plate was modelled with eight-node solid brick elements. The mesh size was such that one element represents the thickness of a layer, i.e., 40 lm. Local mesh refinement was used to resolve the heat source, and coarsening was used to reduce the computational domain where the temperature gradient is small. A 5-s pause was included between each added layer. The laser power was set to 340 W at an efficiency of 30%. The spot diameter was 100 lm and scanning speed was 1000 mm/s. More details pertaining to heating in different models are provided in the following section.
The powder was not modelled explicitly. Instead, elements were activated with bulk material properties at the position of the heat source. A face film boundary condition with a film coefficient of 100 W/ m 2 was used on all faces interacting with the powder bed.

Heat Input Variants
A cylinder with a diameter of 15 mm is built on a 10-mm-thick base plate. Each added layer has a thickness of 40 lm. The laser power, spot diameter, and scanning speed are 340 W, 100 lm, and 1000 mm/s, respectively. The most detailed model, denoted M1, simulates the moving heat source along the actual scanning path, as shown in Fig. 1a. The classical double ellipsoid heat input model with Gaussian distribution proposed by Goldak et al. 12 is used in this case. where Q and g are the laser power and its efficiency, respectively. The variables a, b, and c are geometrical dimensions limiting the region of heat input in a local moving coordinate system. The ellipsoid is a circle in the xz plane with the parameters a and c replaced by the spot diameter d, and b is the penetration depth of 0.16 mm. Goldak et al. 12 suggested that at least four elements should be covering each semi-axis of the ellipsoid. However, we were limited to four elements over the heat source width after mesh refinement because of limited computer memory. Thus, the model does not experience a Gaussian distribution, but the total heat acquired by the model is controlled as described in Ref. 13. This method makes the total heat input independent of discretisation. The first model with temporal reduction is denoted M2, as shown in Fig. 1b. The entire length of a hatch is briefly heated with a Gaussian distribution only in the depth direction. The heat was applied as a box with no overlap or spacing between hatches. The power density Eq. 1 was then reduced to one dimension and rewritten using the above assumptions.
Consequently, the power density becomes independent of the scanning area.
The total supplied energy is the same for all cases. Thus, the integral of the heat source along one scan in model M1 is the energy supplied for each heating line in model M2. The sum of all heating lines in model M2 yields the amount of heating in model M3.
As described above, the net time spent on heating is different for each model, but the total time spent on one layer was set to 5 s for all models.
The total scanning time for each model is a result of chosen process parameters and geometry as illustrated in Fig. 2. The characteristic heating time of any point is the same for all models and is given by The time to sweep the entire area in model M1 is t 1 = A/dv, as shown in Fig. 2a. The time to cover one layer in this model then becomes 1.767 s. However, the scanning path was limited to exclude the edge at the beginning and end of each hatch. This approach enabled numerical stability and avoided extremely high temperatures at the contour. This reduced the total scan time to 1.736 s. Model M2 has an infinite speed for each hatch, and the heating line moves across the area with the scanning speed to maintain the characteristic heating time. Thereafter, the time to cover the whole cross section is t 2 = D/dv. The simplest model, M3, briefly heats the whole area with the characteristic heating, t 3 = t*.

Thermal Properties
The glass-forming alloy AMZ4 is the focus of an on-going research project. This alloy has the capability of becoming amorphous, provided the cooling rate from the liquid state is greater than approximately 20°C/s. The material model is based on thermal properties from the literature. The required high cooling rate limits the possibility to measure material properties at higher temperatures for the amorphous state. However, a complete transformation diagram for crystallization is currently not available. The thermal conductivity of a similar glass-forming alloy, Zr 55 Al 10 Ni 5 Cu 30 , 14 is used because of the lack of published data for AMZ4. Data are available at temperatures of up to the crystallization temperature and are extrapolated to higher temperatures using the slope of available data. Above the liquidus temperature, the conductivity is assumed to increase linearly to imitate convective heat transport. It is considered constant in the superheated liquid state (see Fig. 3).
The heat capacity of AMZ4 has been measured by Heinrich 15 and fitted to Kubaschewski equations 16 for the supercooled liquid region and solid glass state, where the fitting parameters are given as a = 5.224 9 10 À4 J 1 mole À1 K À2 , b = 1.031 9 10 7 J 1 K 1 mole À1 , c = 0.00623 J 1 mole À1 K À2 , and d = À 6.047 9 10 À7 J 1 mole À1 K À3 in Ref. 15. The transition region between these two equations is based on measurements at a gradual heating rate, as indicated in Fig. 3.

Finite Element Models
Simulations with M1-3 start with the same initial state. The initial structure was built on the base plate and consists of 16 layers, with a total height of 0.64 mm. The simulation of the initial state was performed on M3. The simulation time between each layer is 5 s, including heating and subsequent cooling. The mesh was coarsened at the end of the analysis, leaving two layers of elements corresponding to the previous 16 layers (see Fig. 4). The model is cut in a plane perpendicular to the building direction showing the temperature profile.
The computational efficiency of each model was evaluated according to the computer wall-clock time. All models were solved on 16 multiple threads for assembly, recovery, and matrix solving. Model M1 required 381,000 s (approximately 4.5 days) to build one layer. Models M2 and M3 required 1660 s (30 min) and 770 s (13 min), respectively. Thus, only 0.2% of the computational time of M1 is used by M3.  Furthermore, the time increments to heat the entire surface were reduced from 34,700 for M1 to 150 for M2 and only a few for M3.

RESULTS
The cooling rate of the molten material is in the order of magnitudes greater than the critical cooling rate of 20°C/s. Therefore, the material will remain in an amorphous state with certainty. The critical locations are at the lower layers that are merely reheated but not molten. However, they run the risk of devitrification, i.e., crystallization. The choice of the critical location, four layers down, is based on preliminary evaluations. The peak temperature computed by the most detailed model is just above 900°C at this point. Figure 5 shows the temperature history at the selected point versus time for the three models. The cooling rate for the models is evaluated and shown in Fig. 6.

DISCUSSIONS
Model M1 has no temporal reduction and is considered as a reference with respect to simplified cases because measured temperatures are lacking. Its numerical accuracy has been evaluated by changing the temporal and spatial discretisation. Doubling the time step and/or mesh size affect the temperature field near the heat source. However, the effect is negligible at the evaluated critical location, which is four layers down. The heat conductivity used in the model had to be employed for another BMG. Consequently, a numerical experiment was conducted to evaluate whether changing this property will affect the relative behaviour of the models, thereby the conclusion below. Reducing the value of the conductivity by ten times did not change the relative performance of the models.
As further shown in Fig. 5, it is evident that the number of cyclic heating sequences differs among the models. The high temperature region also differs but solutions approximate each other at lower temperatures. However, this is unimportant as they all exhibit extremely high cooling rates. Therefore, this accuracy is not a problem. More relevant for the current application are the predicted cooling rates (Fig. 6). It is possible to define a quantitative error measurement on this rate, but what certainly is important is whether crystallization or devitrification occurs. However, this requires a model that is currently not available. Therefore, a sufficiently accurate model is employed as a conservative model, i.e., one that has a lower cooling rate than the real case. The simplified approaches are conservative in this respect compared with M1.
The critical cooling rate for AMZ4 can be estimated as approximately 20°C/s, with the nose at an undercooling of approximately 80°C. 15 The cooling rates at higher temperatures (see Fig. 6) are in the order of magnitudes greater than this critical value. The risk for devitrification is for reheated layers located further down that experience lower temperatures during the reheating cycle. Any thermal  iii ii i Fig. 6. Cooling rates computed at a location four layers down. Multiple curves for M1 result because of subsequent heating cycles that introduce negative cooling rates. These heating rates are excluded in the graph, as the logarithm of a negative number is not defined. The first three peaks are indicated, and corresponding times can be seen in Fig. 5. history of remelted material can be ignored. According to Fig. 7, all models predict temperatures at these locations that seem to be away from the region for the nucleation of a crystal structure. It should be noted that the figure has no complete data set. The solid line in the upper right corner is based on differential scanning calorimetry measurements of materials cooled at low cooling rates. The symbols at the lower right are from heating experiments. Temperature history lines for models M1 and M3 at the selected point are well to the left of these symbols, indicating a negligible risk of devitrification.

CONCLUSION
The presented work describes an approach to significantly reduce the computational effort to merely 0.2% of a fully detailed model, but one that still ensures sufficient accuracy in the thermal simulation of PBF. The approach includes the temporal reduction of the heating sequence and obtains a conservative estimate of the risk of devitrification of bulk metallic glass.