*2.1. Modeling*

Numerical modeling was based on the fundamentals presented in a previous study [15] in which a 3D microscale thermal model for the EBM process was developed. That work presented a new heat source model tailored according to the actual beam impact on the powder bed and new modeling for the powder thermal proprieties. The beam impact was simulated by Monte Carlo simulations. These simulations determined the effective beam diameter DE that mimics the electron distribution on the powder surface. The energy source model was, therefore, defined as a uniform heat distribution on a circular area in which a coefficient η depends on DE. This was used to couple the electron beam size with the uniform heat source. The material was modeled as a continuum, and the calculation of actual thermal conductivity of the powder considered both thermal radiation in the pores and heat transfer through the powder necks formed during preheating. The powder bed's specific heat and latent heat of fusion were assumed to be equal to those of the solid bulk material [15,16]. The porosity of the powder material, ϕ, was considered equal to 0.32, which corresponds with a particle arrangemen<sup>t</sup> in body cubic centered (BCC) structure. Because of the thermal properties of the powder differ from the bulk (material after the melting) one, the powder was handled as different material, and during the melting simulation, a material change procedure updated the material properties from powder to bulk according to the achieved temperature. In that and subsequent works [15,16], the model was used to predict and validate the temperature distribution and melt pool size during the melting of a single track. Although the accuracy of this kind of modeling, the modeling of the actual EBM process should consider the melting of the entire layer and the addition of the material layer by layer. With this aim, the heat source developed in Reference [15] was adapted to consider the electron beam path across the powder bed. The heat source was, therefore, calculated according to the following equations:

$$\mathbf{q} = \mathfrak{n} \mathbf{U} \mathfrak{l} \mathbf{S} \left( \mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3 \right) \in \Gamma\_{\mathfrak{i}} \tag{1}$$

$$\mathbf{q} = \mathbf{0} \ (\mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3) \in (\mathbf{D} - \Gamma\_i) \tag{2}$$

where η is calculated according to the procedure [15], D is the whole simulation domain, and Γi is the domain of the flux application. The shape of Γi is circular and corresponds to the cross-section of the beam diameter. x1 and x2 are the coordinates of the nodes in the plane of the layer, while x3 is

the coordinate along the build direction, which is the axis normal to the building platform (Figure 2). If nlayer is the total number of layers to be simulated, x3 varies layer after layer according to Equation (3).

$$
\Delta \chi\_3 = - (\mathbf{n}\_{\text{layer}} + \mathbf{n} - 1) \Delta s \tag{3}
$$

The index n is the index of actual layer and Δs is the layer thickness and thus an integer variable, which is calculated (Equation (4)) as the ratio between the total time increment, t, and time to scan a single layer, tlayer.

$$\mathbf{n} = \text{int}(\mathbf{t}/\mathbf{t}\_{\text{layer}}) \tag{4}$$

The beam movement over the x1x2 plane is governed by Equation (5) that describes the selection of the nodes that belong to the domain Γi in the x1x2 plane.

$$(\mathbf{x}\_1 - \mathbf{x}\_{1, \text{inc}})^2 + (\mathbf{x}\_2 - \mathbf{x}\_{2, \text{inc}})^2 \le (\mathcal{O}/2)^2\tag{5}$$

where Ø is the beam diameter and x1,inc and x2,inc are the incremental coordinates from the first position of the beam on the plane (x1c0, x2c0) and the last point of the scanning path. The beam is, therefore, moved point by point along a predefined scanning path. In this work, a standard single direction hatching strategy (Figure 2) was implemented. According to this strategy, the beam is moved with a constant speed, v, across the powder bed along parallel lines. The line offset parameter (LO) establishes the distance between the centerline of two adjacent lines, and all hatch melt lines are in the same direction. At each layer, the scanning path is rotated by 90 degrees clockwise. Therefore, x1,inc and x2,inc are calculated as follows:

$$\mathbf{x}\_{1,\text{inc}} = \mathbf{x}\_{1\text{c}0} + \mathbf{a}\_{12}[\mathbf{v}(\mathbf{t} - \mathbf{n}t\_{\text{layer}}) - \mathbf{m}t\_{\text{track}}] + \mathbf{a}\_{11}\mathbf{m}\mathbf{L}\mathbf{O} \tag{6}$$

$$\mathbf{x}\_{2,\text{inc}} = \mathbf{x}\_{2\text{c0}} + \mathbf{a}\_{22} \left[ \mathbf{v}(\mathbf{t} - \mathbf{n}\mathbf{t}\_{\text{layer}}) - \mathbf{m}\mathbf{t}\_{\text{track}} \right] + \mathbf{a}\_{21}\mathbf{m}\mathbf{L}\mathbf{O} \tag{7}$$

where aii are the coefficients of an antisymmetric matrix that considered the single direction hatching and the scanning path rotation (Figure 2). The aii elements were calculated as follows:

$$\mathbf{a}\_{11} = \mathbf{a}\_{22} = -\sin(\mathbf{n}\pi\mathbf{/2}),\tag{8}$$

$$\mathbf{a}\mathbf{a}\_{12} = -\mathbf{a}\_{21} = \cos(\mathbf{n}\pi\mathbf{/2}).\tag{9}$$

In Equations (6) and (7), m is the line index and represents the actual line to be melted. The line index, m, was calculated from the time for the single line, ttracks, as follows:

$$\text{Im} = \text{int}[(\text{t} - \text{rt}\_{\text{layer}})/\text{t}\_{\text{track}}] \tag{10}$$

Since the simulation here considered the melting of a square of side ltot, ttrack was calculated as follows:

$$\mathbf{t}\_{\text{track}} = (\mathbf{l}\_{\text{tot}} - \mathcal{Q})\mathbf{v} \tag{11}$$

The addition of the material layer by layer was emulated using a new element activation procedure based on the quiet element approach [33]. This method is a numerical approach used in FE simulation for processes that need to account for addition of material during the process. All elements of the model must be modeled from the beginning. Therefore, all elements within the workspace are assembled into the global thermal conductivity matrix, but the elements ye<sup>t</sup> to be deposited are assigned 'void' material properties. Practically, the quiet elements are elements in which the thermal properties are adequately scaled and do not affect the calculation. The domain of analysis is progressively increased according to specific conditions. In this work, the quiet element approach was used not only to simulate the layer addition but also to keep the number of elements involved in the calculation into the current melted layer as low as possible. In this sense, the quiet elements confine the heat exchange only in the HAZ. This strategy is expected to increase the computational speed with respect to the literature models that consider the entire layer at the beginning of the simulation. This method uses a constant mesh and is simple. However, incorrect selection of elements to be eventually activated can lead to ill heat transfer. For this reason, a new criterion was implemented that aims to determine the lowest number of elements belonging to the HAZ without decreasing the accuracy of the heat transfer. Practically, this corresponds to determining the minimum temperature that refers to the location of the HAZ boundary at which the element must be "activated". This temperature is called "activation temperature". Element activation means changing the element material properties to the real ones. This methodology is called thermal activation because it depends exclusively on heat transfer. Since the nodal temperatures of a quiet element are di fferent from those of an active element, body heat flux is numerically generated in some integration points to compensate for the temperature di fference [39]. In this way, when the heat source moves, heat transfer can occur. This propagation generates a temperature field in the nodes of the quiet elements, which are constantly compared with the predefined activation temperature.

The determination of HAZ during the melting involves temperature monitoring and, therefore, experimental trials. To overcome these issues, the microscale model [15,16] was used to determine the HAZ, and thus the activation temperature. Firstly, few melting tracks are simulated using the microscale model. During the simulation, the minimum temperature that corresponds to the boundary of HAZ is identified. This temperature is set as activation temperature in the macroscale model, and the same simulation is performed. The results from the two models are compared in terms of temperature distribution and melt pool size. The macroscale model is then calibrated by changing the activation temperature iteratively until a good match between the two models is observed. With this approach, a trade-o ff between runtime, numerical convergence, and simulation accuracy is guaranteed. Several simulations at di fferent process conditions showed that the activation temperature is slightly higher than the preheating temperature used in the model.

To further reduce the running time, material change was not considered. Therefore, the material properties weremodeled directly as bulkmaterial within the solid part are simulated. However, thematerial surrounding the external surfaces was modeled as powder material in order to account for the heat transfer between the melted and surrounding materials.

For each node of the model, a material state variable ID was defined to indicate the material's temperature history. The material state was defined as "unmelted" (ID = 0) before the melting point was reached and as "melted" (ID = 1) after the melting point was reached. At the beginning of the process, all elements had an ID equal to zero. This variable, therefore, identified the material melted during the EBM process simulation from the one not melted. The melted material state (ID = 1) was used for both the liquid state and for all solid material in the model. The elements for which the ID was still equal to zero at the end of the process were identified as having a lack of fusion areas.

#### *2.2. Model Implementation*

The thermal model was implemented in Abaqus/Standard using user code subroutines for the heat source motion, element activation, and molten material representation. The whole domain is rather simple and consists of a solid substrate on which a set of layers are laid. To reduce the calculation time, the analysis domain used was a small portion of a specimen produced by EBM process with an external lateral surface in contact with a certain volume of loose powder. This powder volume was included to consider the heat exchange between the melted part and the surrounding area. The whole domain was divided into four subdomains (Figure 2): the bulk substrate, the powder substrate, the bulk quiet domain, and the quiet powder domain. The bulk substrate represented material already processed and not simulated in the analysis. The powder substrate was adjacent to the solid substrate. Therefore, a vertical surface of the bulk domain was numerically coupled with one surface of the powder domain. The other surfaces were considered adiabatic. The quiet bulk domain represented the domain scanned by the electron beam. It contained quiet elements that were activated once they achieved the activation temperature. The quiet powder domain represented a portion of material adjacent to the quiet bulk domain that was not scanned. This domain contained quiet elements that were automatically activated at the beginning of the layer and simulated heat transfer from the electron beam track and the surrounding powder. The quiet powder and bulk domains share a lateral surface. The other surfaces were considered adiabatic.

The mesh included DC3D8 (heat transfer brick) elements with a size equal to 0.05 mm in the quiet bulk and powder domains. The rest of the domain was discretized with a progressive coarsening mesh of up to 0.5 mm with DC3D4 elements (heat transfer tetrahedron).

**Figure 2.** Finite element (FE) model configuration.

The overall domain was subjected to a predefined field corresponding to the preheating temperature. To simulate the preheating step, the elements activated during the simulation had a preheating temperature.

According to the EBM melting strategy, the heat source swiped across the layer using a hatching mode. A specific user code, DFLUX, controlled the electron beam movement according to a user code that implemented the Equations from (1) to (11). Initially, the heat source moved at a constant speed along a certain axis. At each time increment, the code identified the nodes with the domain Γi and applied the heat flux according to Equations (1) and (2). After the scanning of a track, the heat source was repositioned at the origin of the previously scanned line. Then, it was translated perpendicularly to the previously scanned line and was of a quantity equal to LO. Hence, a new track was melted. The process was repeated until the layer was completed. Since only a portion of the whole cube was simulated, the code also considered the idle times due to the non-simulated portion of the layer.

To consider the melting strategy rotation, the scan lines were tilted on every layer by 90 degrees (Figure 2). Therefore, the code considered the heat source path rotation and its translation along the build direction. Both the path rotation and the translation were controlled by the time and the scanning speed.

The addition of the material during the simulation was implemented using the UEPACTIVATION VOL subroutine. This subroutine allows for the element activation control using rules predefined by the user. According to the thermal activation described in the previous paragraph, the elements are activated if the nodal temperature (T) is higher than the activation temperature (Tact). Owing to the Abaqus limitation, at the beginning of each layer, a few elements must be activated to trigger the subsequent thermal activation. Therefore, a predefined set of elements corresponding to the beam cross-section was activated by default at the beginning of each layer. This initial set was concentrically positioned at the center of the first beam spot. During subsequent increments, the elements were activated via the achieved temperature. Because of the heat propagation, it was possible that subsequent melting would determine the activation of elements belonging to the underlying layer if the temperature was higher than the activation temperature.

At the beginning of the simulation, the ID value at each node was equal to zero. At each increment, the ID was updated using the USDFLD subroutine. This user code identified and collected the nodes that had achieved the solidus temperature (TSolidus), while the temperature at the material point for each increment (T) was read using the GETVRM subroutine. Practically, all ID values were read at each node at the end of each increment, and the ID values equal to zero were updated (from 0 to 1) if the node achieved the melting point.

Figure 3 summarizes the FE whole structure implementation. Appendix B reports the source code of DFLUX and UEPACTIVATIONVOL subroutines as implemented in Abaqus.

**Figure 3.** FE model structure and the subroutine for the material addition (UEPACATIVATION VOL), the material state variable for the melted material (USDFLD), and heat flux application (DFLUX) implementations.

#### *2.3. Model Validation*

As mentioned above, the proposed model needs proper calibration before it can be run. The calibration is performed using the microscale model [15], named the reference model. Practically, the activation temperature in the macroscale model is iteratively changed until similar heat transfer conditions between the models are obtained.

First, to demonstrate the reliability of the heat transfer obtained via the proposed model after the calibration, a preliminary validation shows the comparison between the proposed and the reference models in terms of temperature distribution and melt pool size. This analysis considered the melting of a portion of a single layer. Then, further validation evaluates the ability of the model to predict the lateral surface roughness. In this case, a multilayer simulation was performed and compared with the experimental data.

In order to give a measure of the computational time, all simulations were performed on a personal computer with the following specifications: Intel® Core ™ i7-8700K 3.70 GHz processor, 32 GB RAM, and 1 TB SSD. A six-core parallelization was used according to the computational resources used to perform the analysis.

#### 2.3.1. Comparison with the Reference Model: Single Layer Validation

This validation was carried out by simulating the melting of a single layer, 30 mm × 30 mm, of Ti-48Al-2Cr-2Nb alloy. The layer was 0.090 mm thick. The material properties for the bulk material were extracted from Reference [40], while the powder proprieties were calculated from the bulk one according to Reference [15]. The material properties implemented in the model were reported as Supplementary Materials. The simulation by the proposed and reference models consisted of the

melting of six unidirectional parallel scan tracks. Table 1 lists the process parameters used for both simulations. The comparison considered both the melt pool size and temperature distribution.


**Table 1.** Simulation parameters for both analyses.

#### 2.3.2. Multilayer Experimental Validation

The experiment consisted of producing a 10 mm × 10 mm × 20 mm-sized parallelepiped sample made on Ti6Al4V. The sample was produced using an Arcam A2X, standard Arcam powder (average particle size equal to 75 μm) with a layer thickness equal to 0.050 mm. The process parameters are listed in Table 2. According to the standard melting strategy, the melting strategy consisted of unidirectional parallel lines with a clockwise rotation of 90 degrees every layer (Figure 2). After the production, the sample was cleaned by loose powder. The 2D surface roughness profiles of the lateral surface samples were acquired using an RTP-80 profilometer, equipped with a TL90 drive unit. The cut-off length and sampling length were chosen according to ISO 4288:1997. For each lateral surface of the sample, the roughness profile was measured along the build direction. Additionally, for each surface, three profiles were collected: one on the center of the specimen and the other two equally spaced from the first one and from the sample's edges. The collected profiles for each surface are reported in Appendix A, while all collected data and roughness descriptors can be found in the database [41]. The average roughness (Ra) and root mean square (Rq) values were equal to 28 μm and 33 μm, respectively.

**Table 2.** Process parameters for the experiment and simulation.


The numerical simulation by the proposed model emulated the experimental setup. The material properties for Ti6Al4V were extracted from a prior study [15]. The layer was fixed to be equal to 0.050. Figure 2 shows the portion of the simulated specimen. The simulated portion of the cube was 2.2 mm × 2 mm × 1.5 mm. As described above, the analysis domain consisted of four subdomains. The quiet bulk domain contained 10 layers (soft green with black dotted edges in Figure 2) that were added over a solid substrate (soft blue with black edges in Figure 2). The powder substrate (soft blue with orange edges) shared a lateral surface with the solid substrate, while the quiet powder domain (soft green with orange dotted edges in Figure 2) was located over the powder substrate and adjacent to the quiet bulk domain. According to the proposed procedure, at the beginning of each layer, the initial activated area was cylindrical and concentric to the initial position of the beam. The height of the

cylinder was set to be equal to the layer thickness, and its diameter was equal to the beam diameter. The video of the simulation after the model calibration is provided as Supplementary Materials.

#### **3. Results and Discussion**

#### *3.1. Model Validation*

#### 3.1.1. Single-Layer Simulation Validation

As mentioned above, the reliability of the proposed model response was compared with the reference one in terms of melt pool dimensions and maximum temperature (Tmax). The length (l), the width (w), and the depth (d) of the melt pool were measured. The activation temperature was fine-tuned to a convergence in the heat transfer between the two models. The model was considered calibrated when small deviations were observed between the two models. For an activation temperature equal to 1700 K, the deviations from the reference model are reported in Table 3. This deviation can be attributed to the di fferent material properties attributed to the layer. The reference model, in fact, considered the simulation of the entire layer with the powder thermal material properties that were then updated from powder to bulk according to the achieved temperature during the process. In the proposed model, the portion of the layer to be melted was modeled directly as bulk material. Because of that, deviations were more significant in the first track. In this case, the beam should melt a track surrounded completely by material with low conductivity (powder), while in the proposed model, the powder was modeled only on one side of the beam track. In subsequent passages, the beam overlapped with a portion of the already solidified bulk material, and the deviation between the two models was, therefore, lower, as shown in the third scanning line (Table 3). The di fference in the melt pool depth remained, implying that there was an overestimation of the heat penetration in the proposed model. However, the penetration depth was also a ffected by uncertainty due to the substrate temperature after melting, which could not have been forecasted with a single-layer simulation.


**Table 3.** Single-layer validation: the deviations between the reference model against the proposed model in terms of the melt pool dimension and maximum temperature.

Overall, the numerical results obtained adopting the new modeling approach can be considered to be consistent with the reference model.

Additionally, using the same computational resources, the total simulation with the reference model took around 8 CPU hours, while the same simulation with the proposed model after calibration took around 5 h. Therefore, the proposed approach to the EBM modeling with proper calibration of the activation temperature allows to obtain an accurate enough heat transfer simulation and a significant runtime saving of about 40% with respect to the reference model.

#### 3.1.2. Lateral Surface Roughness Prediction

According to the proposed approach, firstly, few melting lines were simulated using both the reference and the proposed models. Then, the proposed model was calibrated by changing the activation temperature iteratively until a heat transfer convergence between the two models. For an activation temperature equal to 1000 K, the obtained numerical results from the calibrated model showed a small deviation in comparison with the corresponding reference model measurements. For this reason, the activation temperature was set equal to 1000 K.

To be numerically consistent with the experiment, the surface roughness was described by the surface profile between the quiet bulk and powder domains, considering the melted material. This information was extracted by the material state variable ID. As mentioned above, the ID was equal to 1 if the node achieved the solidus temperature, otherwise, it remained equal to zero. Graphically, this result can be represented by a colorful map in which the elements in red were an ID equal to 1 and those in blue were an ID equal to zero. Figure 4 shows the experimental 2D roughness profile (blue line) overlapped with the numerical one. The numerical profile was obtained using a cross-section of the model. The surface roughness profile of the sample (blue line) was characterized by peaks and valleys, which were repeated with a certain frequency. This frequency seems to be linked to the rotation layer and, therefore, to a distance equal to four layers. This frequency was also observed in the numerical results. The peaks of the numerical simulation appeared to be representative of the experimental ones, while the numerical valleys appeared to be more profound. This could be explained by the presence of powder particles attached to the experimental surface because of the high sintering level [42,43]. These particles could not be removed during cleaning because they were partially melted into the surface. Since the simulation approach used the FE method and a continuum domain, the presence of these particles was not detected.

**Figure 4.** Comparison between the numerical and experimental roughness profiles.

The di fferences of the wave width between the experimental and the numerical profiles can be explained by the e ffective layer thickness [12], which was around 2 or 3 times larger than the lowering of the build platform [23] and was not considered in this simulation.

A certain e ffect may also be caused by the material change, which was not considered in the simulation and may increase the heat transfer toward the powder domain because of the local thermal conductivity change.

Overall, the numerical simulation appeared to represent the actual surface profile well and established the capability of the model to give information about the surface texture resulting from the EBM process.

Figure 5 shows a portion of the 3D profile of the lateral surface obtained using numerical simulation. This profile was obtained by manipulating the numerical output visualization to show the boundary surface between the melted and unmelted material. Figure 5 detects the peaks and valleys along the surfaces previously observed in both 2D experimental and numerical profiles.

To estimate the average roughness (Ra) and the root mean square (Rq) of the surface and to increase the sampling length in the model, a longer simulation was performed up to 17 layers, which corresponds to the height of the simulated sample equal to 0.85 mm. From this simulation, Ra and Rq resulted in being equal to 35 μm and 42 μm, respectively (see Appendix A for details). The obtained values are comparable with the experimental one (Ra = 28 μm and Rq = 33 μm). Therefore, the results obtained from the model are believed to serve as good indicators of the lateral surface profiles. However, as will be discussed in the Conclusion, a more detailed investigation of both experimental and simulated melting is planned.

#### *3.2. Lack of Fusion in EBM Process*

Given the results presented in Section 3.1.2, further analysis was carried out to evaluate the model response in predicting the lack of fusion under a critical process condition. Lack of fusion occurs easily if the EBM conditions are not appropriate. For instance, an increased line o ffset toward beam diameter results in lower energy densities, which causes a lack of fusion [44]. With this scope, the line offset of the previous multilayer simulation was increased up to 0.800 mm, which was approximately equal to the beam diameter (Table 2). The other parameters remained constant, including the activation temperature.

**Figure 5.** 3D numerical surface profile: peaks and valleys can be observed.

Numerically, the lack of fusion was represented by those elements named "unmelted material" and, therefore, did not achieve the melting point during the simulation. The ID variable gathered the state of the element. As before, within the quiet bulk domain, the unmelted elements (ID = 0) were plotted in blue (Figure 6), while the melted elements (ID = 1) were plotted in translucent. As can been seen, Figure 6 shows the presence of numerous lack-of-fusion areas entrapped in the solid material. The locations of these unmelted areas look to be dependent on the layer rotation strategy. The simulation, therefore, identified this process condition as faulty. In fact, larger line offset values that did not allow for the full melting of the material, causing the lack of fusion appearance, voids, and porosities, thus inhibiting the continuation process [18].

**Figure 6.** Simulation of a process condition that causes the appearance of a lack of fusion. The translucent elements indicate the material that achieved the melting point during the process (ID = 1), while blue ones indicate the material that has never achieved the melting point (ID = 0). After the cooling phase, the unmelted material during the process represents a lack-of-fusion defect in the final part.

#### *3.3. Runtime Considerations*

The total duration of the 10-layer simulation was around 20 h. The process simulation with a larger line offset took around 2 h.

A convergence analysis was carried out on mesh size during the activation temperature calibration. The mesh size was finer than 0.050 mm and significantly increased the runtime without increasing the

output accuracy. Lower activation temperatures lead to the activation of a higher number of elements, therefore, worsening the runtime. Higher temperatures caused the non-convergence of the model because of the incorrect heat transfer.
