**4. Conclusions**

The literature review of multilayer simulations for the EBM process and more in general for AM processes is poor and did not allow the prediction of the local quality. In this work, an accurate but simple 3D multilayer model was developed for the optimization of process parameters. The aim of the model was to rapidly identify a set of potentially faulty process conditions. To speed up the runtime, we developed a new tailor-made procedure. This procedure was based on the progressive element activation under a specific temperature, which was used not only to simulate the layer by layer production but also to limit the number of elements that participated in the calculation during the single-layer melting. The activation temperature was determined by accurately analyzing the HAZ simulation using a microscale model that predicted the temperature distribution and the actual melt pool. The activation temperature was calibrated, and a small di fference in heat transfer was observed between the two models.

To identify the quality inherent in the lack of fusion, our methodology included the definition of a material state variable (ID), which allowed us to distinguish between melted and unmelted material. The graphical visualization of ID evaluated the roughness profile of a sample surface along a build direction, as well as the lack of fusion appearance. For validation purposes, the comparison between the microscale model used for the activation temperature calibration and the proposed one was shown. Successively, the model's ability to predict surface texture was verified against experimental results. The model response was tested under a critical condition that, in real production, generated a process failure. The right activation temperature led to a small deviation between the microscale and the proposed models in terms of heat transfer with a savings time of about 40%. The model showed good agreemen<sup>t</sup> with the experimental and numerical measurements in forecasting the 2D surface profile and the roughness descriptors (Ra and Rq). The little deviation between the experimental and the numerical demonstrates that the proposed approach will be an e fficient model for simulating the roughness profile. However, a more detailed numerical analysis will be carried out for the verification of the model capability in roughness prediction at di fferent process conditions. Additionally, further work will consider implementing the e ffect on the powder layer due to account both for thermal expansion of powder particles and for porosity reduction of the powder bed during the heating [23].

The model also demonstrated the ability to predict a lack of fusion associated with thermal effects due to, for example, a large line o ffset. In such a case, with respect to the experimental trial, the simulation allowed for hugely significant time-saving process optimization. From this point of view, the presented results demonstrated that the proposed approach could be easily applied to the process development stage in order to reject faulty process conditions. From an industrial point of view, this approach can save time and resources with respect to the experimental routine. Further model validations and developments with the aim of simulating complex parts may help in the production of defect-free parts. Additionally, an accurate prediction of surface profile roughness may allow for the identification of minimum allowances, which are needed to add for the subsequent finishing operations. Moreover, there is a need for microcrack identification, which can jeopardize a product's fatigue limits.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4352/10/6/532/s1, **S1**: Video: Multilayer process simulation of the EBM process for Ti-6Al-4V: the size of the simulated domain of material addition is 2.2 mm × 2 mm × 1.5 mm. The domain contains 10 layers with a layer thickness equal to 0.050 mm. The process parameters are: beam diameter = 0.780 mm; Line o ffset = 0.200 mm; scan speed = 471.69; Heat flux = 347.920 W/mmˆ2; Preheat temperature = 923 K. The movie shows the temperature distribution evolution over the time (on the right) and the corresponding melted material (on the left). **S2**: Thermal material properties for Ti-48Al-2Cr-2Nb implemented in the simulation.

**Author Contributions:** Conceptualization, M.G.; data curation, M.G. and O.D.M.; investigation, M.G. and O.D.M.; methodology, M.G.; resources, L.I.; software, M.G. and O.D.M.; supervision, L.I.; validation, M.G. and O.D.M.; visualization, M.G. and O.D.M.; writing—original draft, O.D.M.; writing—review and editing, M.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Experimental and Numerical Roughness Profile and Ra and Rq Calculation**

Figure A1 depicts all surface roughness collected during the experiments. For each surface of the sample, the roughness profiles have been acquired along the build direction and replicated three times. The first profile was acquired along the surface center, while the other two profiles were equally spaced from the first one and the edges of the sample. These measurements were collected and indexed as follows: the first number (from 1 to 4) indicates the index of the surface; the second number (from 1 to 3) specifies the measure replica. The surface roughness profiles of the sample were acquired using an RTP-80 profilometer, equipped with a TL90 drive unit. The cut-o ff length and sampling length were chosen according to ISO 4288:1997.

**Figure A1.** Surface roughness profiles from each surface. The measurements have been replicated three times. The first number (from 1 to 4) indicates the index of the surface; the second number (from 1 to 3) specifies the measure replica.

For each measurement, the parameters Ra and Rq and the surface profile were also collected. The average value of the Ra and Rq was equal to 28 μm and 33 μm, respectively. All collected measurements and roughness descriptors can be found in the database [41].

A portion of the samples has been simulated numerically according to the proposed methodology. In order to have a significant and reliable number of roughness values to calculate the predicted Ra and Rq values, the addition of 17 layers over a solid substrate has been simulated. According to the methodology presented in the paper, the boundary between the unmelted and the melted materials has been extracted from the simulation results by using the ID variable (Figure A2a). According to the standard [45], Ra is calculated, over the evaluation length, as the arithmetical mean of the absolute values of the profile heights (yi) from the mean line of the e ffective profile. The mean line is the line having the shape of the geometrical profile and dividing the e ffective profile (Figure A2) so that within the sampling length, the sum of the squares of the distances between the e ffective profile points and the mean line is minimum. Analogously, Rq is the root mean square average of the profile heights over the evaluation length. From the numerical model, Ra and Rq resulted in being equal to 35 μm and 42 μm, respectively.

**Figure A2.** Numerical measurements of the surface roughness: (**a**) cross-section of the model that shows the ID state variable map (**b**) Numerical profile and its mean line for Ra calculation.

#### **Appendix B. Source Code Listing of DFLUX and UEPACTIVATIONVOL Subroutines as Implemented in Abaqus**

Appendix B lists the user code as written in Abaqus to implement the modeling as discussed in the main paper. The subroutines are in FORTRAN language.

SUBROUTINE DFLUX(FLUX,SOL,KSTEP,KINC,TIME,NOEL,NPT,COORDS,

```
JLTYP,TEMP,PRESS,SNAME)
```

```
INCLUDE 'ABA_PARAM.INC'
DIMENSION FLUX(2), TIME(2), COORDS(3)
CHARACTER*80 SNAME
c
ltot = 10 !scan track lenght
D = 0.78 !beam diameter
V = 471.69 !beam speed
Q = 347920 !heat flux magnitude
LO = 0.2 !line o
               ffset
LT = 0.05 !layer thickness
n_layer = 10 !number of layers
c
```
x1 = COORDS(3) !spatial coordinate x1

```
x2 = COORDS(1) !spatial coordinate x2
       x3 = COORDS(2) !spatial coordinate x3
       t = TIME(1) !total time
     c
     x1c0 = 0 !coordinates of the first position of the beam
       x2c0 = −D/2 !coordinates of the first position of the beam
     c
     t_track = (ltot-D)/v !time to perform a single line
     t_layer = ltot/LO*t_track !time to scan a single layer, tlayer
     c
     n = int(t/t_layer) !layer index
     m = int((t-n*t_layer)/t_track) !line index
     c
     a11 = −sin(n*(pi/2)) !matrix coefficients
       a22 = −sin(n*(pi/2))
       a12 = cos(n*(pi/2))
       a21 = −cos(n*(pi/2))
     c
     x1_inc = x1c0 + a12*(v*(t-n*t_layer)-m*t_track) + a11*m*LO !incrementalandthelastpointofthescanningpath
```

```
x2_inc = x2c0 + a22*(v*(t-n*t-layer)-m*t_track) + a21*m*LO !incremental coordinates from x2c0
and the last point of the scanning path
```
 "quiet"

coordinates

 from x1c0

c

c

if(x3.EQ.-(n\_layer + n − 1)\*LT)then !check the x3

if ((x1-x1\_inc)\*\*2 + (x2−x2\_inc)\*\*2.LE.(D/2)\*\*2)then

FLUX(1) = q !Heat source application. Since only a portion of the whole cube was simulated, end if !was considered the idle times due to the non-simulated portion of the layer end if

```
RETURN
END
c
```
subroutine UEPACTIVATIONVOL(lFlags, epaName, noel, nElemNodes, iElemNodes, mcrd, coordNodes, uNodes, kstep, kinc, time, dtime, temp, npredef, predef, nsvars, svars, sol, solinc, volFract, nVolumeAddEvents,volFractAdded,csiAdded)

 include 'aba\_param.inc' character\*80 epaName c ltot = 10 !scan track lenght D = 0.78 !beam diameter V = 471.69 !beam speed LO = 0.2 !line offset LT = 0.05 !layer thickness n\_layer = 10 !number of layers Tact = 1000 !activation temperature c At the beginning of the simulation the status of all elements must be set asvolFractAdded(1) = 0.d0 t = TIME(1) !total time x1 = coordNodes(3,1) !spatial coordinate x1 x2=coordNodes(1,1)!spatialcoordinatex2

x3 = coordNodes(2,1) !spatial coordinate x3 c x1c0 = 0 !coordinates of the first position of the beam x2c0 = −R !coordinates of the first position of the beam n = int(t/t\_layer) !layer index c Activation of the initial set of elements if(x3.EQ.-(n\_layer+n-1)\*LT)then if ((x1-x1c0)\*\*2+(x2-x2c0)\*\*2.LE.(D/2)\*\*2)then volFractAdded(1) = 1.d0 end if end if c Thermal activation if (x3.LE.-(n\_layer+n-1)\*LT)then if(sol(1).GE. Tact) then end if end if c return end
