**Performance of Optimization Algorithms in the Model Fitting of the Multi-Scale Numerical Simulation of Ductile Iron Solidification**

### **Eva Anglada 1, Antton Meléndez 2, Alejandro Obregón 2, Ester Villanueva <sup>2</sup> and Iñaki Garmendia 3,\***


Received: 1 July 2020; Accepted: 5 August 2020; Published: 8 August 2020

**Abstract:** The use of optimization algorithms to adjust the numerical models with experimental values has been applied in other fields, but the efforts done in metal casting sector are much more limited. The advances in this area may contribute to get metal casting adjusted models in less time improving the confidence in their predictions and contributing to reduce tests at laboratory scale. This work compares the performance of four algorithms (compass search, NEWUOA, genetic algorithm (GA) and particle swarm optimization (PSO)) in the adjustment of the metal casting simulation models. The case study used in the comparison is the multiscale simulation of the hypereutectic ductile iron (SGI) casting solidification. The model fitting criteria is the value of the tensile strength. Four different situations have been studied: model fitting based in 2, 3, 6 and 10 variables. Compass search and PSO have succeeded in reaching the error target in the four cases studied, while NEWUOA and GA have failed in some cases. In the case of the deterministic algorithms, compass search and NEWUOA, the use of a multiple random initial guess has been clearly beneficious.

**Keywords:** model fitting; optimization; FEM; metal casting; SGI; numerical simulation; compass search; NEWUOA; genetic algorithm; particle swarm optimization

### **1. Introduction**

Today, numerical simulation of metal casting processes at macro scale is a well-known technology. Although new approaches based for example in dynamic mesh techniques are being currently investigated [1], the use of commercial programs based on the finite element method (FEM) or on the finite volume method (FVM) is quite widespread in the industry. However, the agreement between the obtained predictions and the actual results obtained at the foundry plant, is not always the desired one. This agreement is more difficult even in the case of multiscale simulation. For this reason, the model fitting is of high interest to improve the confidence in the predictions. Even more in the present context where, in one side the integrated computational materials engineering and in other side the hybrid prognostic models combining numerical simulation and data-driven models applied to the Industry 4.0 movement of the manufacturing industry, appear as two of the most promising and challenging research trends [2].

The obtaining of reliable predictions in metal casting processes is not a trivial task. The results do not only depend on numeric codes in which they are based. The use of the appropriate values of the thermophysical parameters along with the proper establishment of boundary conditions, is essential to avoid poor results or even erroneous ones. In fact, the ASM Committee [3] remarks that each manufacturing process has unique boundary conditions that can be even equipment specific, which must be identified and characterized for the specific application being simulated.

The characterization of the material properties for the whole range of temperatures involved in the metal casting processes is very difficult. Additionally, the determination of the values of some boundary conditions, as for example the heat transfer coefficients (HTC) between the alloy and the mold, is even more complex. The limitations of bibliographic data and the extreme difficulty of the direct determination of these values have inspired the application of inverse methods to avoid these difficulties.

The inverse methods can be considered a subset of the field of numerical optimization where the objective is to minimize one objective function, which represents the error between the simulation results and the experimental measurements. The challenge is that the corresponding objective function can be highly nonlinear or non-monotonic, may have a very complex form or its analytical expression may be unknown. Moreover, in the case of metal casting inverse problems, the effect of changes in boundary conditions are normally damped or lagged, i.e., the varying magnitude of the interior temperature profile lags behind the changes in boundary conditions and is generally of lesser magnitude. Therefore, such a problem would be a typically ill-posed and would normally be sensitive to measurement errors. Additionally, the uniqueness and stability of the solution are not generally guaranteed [4–6].

This type of optimization problems is usually solved using iterative strategies. The methodology consists basically in an iterative process where several parameters of the model are modified until the simulation predictions agree with the experimental measurements. Different strategies can be applied to advance in the iterative process depending on the selected optimization method.

In general, it is not possible to establish which optimization methods are the best between the high number of alternatives that exist. The performance of the optimization methods is completely related with the problem type to be tackled. Therefore, for each problem type, the most appropriate method should be searched.

The use of optimization algorithms to correlate the numerical models with experimental values has been applied in other fields different from the metal casting. The sector where possibly more efforts have been done in this sense, is the space industry. The correlation of the spacecraft thermal mathematical models with the results of the thermal tests is mandatory, which has encouraged the evaluation of different solutions. Some authors have explored the use of deterministic methods, as Klement, who used several quasi-newton type algorithms based on the Broyden methods [7,8] or Torralbo et al. who used a generalized iterative pseudo-inverse algorithm [9,10]. But most of them have studied and compared different stochastics methods. For example Dudon [11] use the Latin hypercube sampling (LHS), the self-adaptive evolution (SAE) and the branch and bound methods; Beck [12] use the adaptive particle swarm optimization (APSO); Van Zijl [13] use the Monte Carlo, the genetic algorithm (GA) and the APSO methods; Trinoga and Frey et al. [14,15] use the simulated annealing (SA), the threshold accepting (TA) and the GA; Anglada et al. work mainly with GA [16,17] but also made some comparisons with Klement algorithms [18] and with the deterministic algorithms TOLMIN, NEWUOA, BOBYQA and LINCOA [19]. In addition, works about the use of hybrid algorithms combining deterministic and stochastic methods can also be found, as the works of De Palo [20] where the downhill simplex is combined with the LHS or the work of Cheng et al. [21] where the Broyden–Fletcher–Goldfarb–Shanno (BFGS) is combined with the Monte Carlo method.

In contrast, the efforts done in the field of metal casting are much more limited. The number of bibliographic references devoted to the application of optimization algorithms to the adjustment of metal casting simulation models is quite reduced. Most of them [22–25] use the maximum A posteriori (MAP) algorithm, which is a variation of the least-square technique explained in reference [26]. Others perform the adjustment manually [27] or combining the manual adjustment with the MAP

algorithm [28–32]. Moreover, some authors have tested different methods as Ilkhchy et al. [33] which use a nonlinear estimation technique similar to that stated by Beck in [34]; Dong et al. [35] that use a combined function specification and regularization method; Zhang et al. [36] with a neural network or Zhang et al. [37] that use a globally convergent method (GCM).

The authors of the present study believe that a deeper study in the performance of different optimization methods may contribute to give a step forward in the state of art of the model fitting applied to the metal casting. It would make possible to obtain adjusted models in less time and moreover will open the door to the future integration of automatic model fitting algorithms in prognostic models to be used in the framework of the industry 4.0.

The work presented henceforward compares the performance of several algorithms, deterministic and stochastic, in the adjustment of the multiscale simulation model of the solidification of a ductile iron casting, with the objective of contributing to the development of this research line in the search of a methodology solid enough to be used in the adjustment of metal casting simulation models.

### **2. Methodology**

The methodology used for the adjustment of the simulation model is summarized in Figure 1.

**Figure 1.** Flowchart of the model adjustment procedure.

The multiscale simulation model of the metal casting process was done with the commercial finite element software ProCAST 2018.0 (ESI Group, Paris, FRANCE), focused on metal casting simulation [38]. The fitness function represents the error between the simulation model predictions and the experimental results, that in this case are the alloy tensile strength. The parameters modification is guided by the optimization method used in each case, which was implemented by means of the PyGMO (Parallel Global Multiobjective Optimizer) library. PyGMO is a scientific library for massively parallel optimization developed in the context of the European Space Agency to evolve interplanetary spacecraft trajectories as well as designs for spacecraft parts and more [39].

### *2.1. Simulation Model*

### 2.1.1. The Physics

Different approximations have been studied for the prediction of the tensile strength of cast alloys, as the use of machine learning algorithms [40] or analytical approximations as the Equation (1) stated by [41]. In this case, this last approximation, Equation (1), has been used by the simulation model to predict the tensile strength at room temperature (σ*ULT*), where *fg* is the graphite fraction, *n* is the shape factor for graphite nodules which is related with their sphericity (maximum value, that is a perfect sphere, is unity), *f*<sup>α</sup> is the ferrite fraction, *fp* the perlite fraction, and (*dT*/*dt*) <sup>850</sup> is the cooling rate at 850 ◦C. Therefore, the microstructure and the cooling rate must be previously calculated.

$$
\sigma\_{ULT} = \left(1 - f\_{\mathcal{g}}^n\right) \left(482.2f\_a + 991.5f\_p\right) + 50\left[(dT/dt)\_{850} - 0.5\right],\tag{1}
$$

In cast irons the carbon and silicon contents are high, compared with steels, which origins a rich carbon content phase in their structure. When the alloy graphitization potential is high, which is determined by its composition and melt treatment, the iron solidifies according to the stable Fe-graphite system and the rich carbon phase is graphite, a hexagonal-close-pack form of carbon. In addition, the presence of magnesium makes the graphite structure be spheroidal instead of lamellar. One of the main melt treatments is the inoculation, which consists in the deliberate addition of elements that promote more active graphite nucleation sites (usually ferrosilicon enriched with Sr, Ba, Al, Zr, Ca and with very reduced and well calculated quantities of rare earths). The matrix structure is essentially determined by the composition and the cooling rate through the eutectoid temperature range. The eutectoid reaction leads to the decomposition of austenite into ferrite and graphite for the case of the stable eutectoid and to pearlite for the metastable eutectoid transformation. Slower cooling rates result in more stable eutectoid structure promoting a more ferritic matrix. If the complete transformation of austenite is not achieved when the metastable temperature is reached, pearlite forms and grows in competition with ferrite.

In the metal casting model used in this case, the influence of the chemical composition is managed by means of the thermodynamic database for multicomponent iron-based alloys developed by CompuTherm, PanFe 2017 (CompuTherm LLC, Middleton, WI, USA) [42]. This thermodynamic database is based on the calculation phase diagrams (CALPHAD) methodology, also called computer coupling of phase diagrams. CALPHAD is a materials genome method which plays a central and fundamental role in materials design, see reference [43] for more information. This method develops models to represent thermodynamic properties for various phases. The thermodynamic properties of each phase are described through the Gibbs free energy, based on the experimental and theoretical information available on phase equilibria and thermochemical properties in a system. Following this, it is possible to recalculate the phase diagram, as well as the thermodynamic properties, of all the phases and the system as a whole. It allows to reliably predict the set of stable phases and their thermodynamic properties in regions without experimental information and for metastable states during simulations of phase transformations together with the properties of multicomponent systems from those of binary and ternary subsystems. This is accomplished by considering the physical and chemical properties of the system in the thermodynamic model. In this case, the "lever" micro-segregation model, which assumes a complete mixing of the solute in the solid, was applied.

The influence of the melt treatment is managed by means of the nucleation parameters which must be specified/adjusted for each case, as they are not an intrinsic property of the material.

• The nucleation of the primary dendritic phase, austenite, is modeled following the gaussian distribution model proposed by [44]. This model defines the relationship between the number of nuclei and the undercooling (Δ*T*) following Equation (2), where *nmax* is the maximum grain nuclei density, Δ*T*<sup>δ</sup> the undercooling standard deviation and Δ*Tn* the average undercooling (see Figure 2):

$$m(\Delta T(t)) = \frac{n\_{\text{max}}}{\sqrt{2\pi} \cdot \Delta T\_{\delta}} \int\_{0}^{\Delta T(t)} \exp\left(-\frac{\left(\Delta T(t) - \Delta T\_{n}\right)^{2}}{2\Delta T\_{\delta}^{2}}\right) d(\Delta T(t)),\tag{2}$$

• The nucleation of the graphite nodules is calculated as a power law of the undercooling following Equations (3) and (4) of the model proposed by [45], where *Ae* and *n* are the nucleation constants. This model assumes bulk heterogeneous nucleation at foreign sites which are already present within melt or intentionally added to the melt by inoculation:

$$N\_{\rm cut} = \mathcal{A}\_{\rm c} (\Delta T)^{\rm u},\tag{3}$$

$$\frac{dN\_{\rm cut}}{dt} = -nA\_{\rm ef} (\Delta T)^{n-1} \frac{dT}{dt} \,\,\,\,\,\tag{4}$$

• The graphite nodules growth is calculated with a quadratic power of the undercooling following Equation (5), where μ*<sup>e</sup>* is the eutectic growth coefficient:

$$\frac{d\mathcal{R}\_{\mathcal{S}}}{dt} = \mu\_{\mathcal{E}} (\Delta T)^2,\tag{5}$$

**Figure 2.** Undercooling definition.

Relating to the cooling rates, the governing equations of the physics involved in the mold filling, and in the alloy cooling and solidification, are the energy conservation, Equation (6), the mass conservation or continuity, Equation (7), and the Boussinesq form of the Navier–Stokes equation for incompressible Newtonian fluids, Equation (8), where ρ is the density; *h* is the specific enthalpy; *t* is the time; <sup>υ</sup> is the velocity; *<sup>k</sup>* is the thermal conductivity; . *Rq* is the heat generation per unit mass; ρ<sup>0</sup> is density at reference temperature and pressure; *p*ˆ is the modified pressure (*p*ˆ = *p* + ρ<sup>0</sup> *gh*); μ*<sup>l</sup>* is the shear viscosity; *g* is the gravity; β*<sup>T</sup>* is the volumetric thermal expansion coefficient; *T* is the temperature and *T*<sup>0</sup> is the boundary temperature. More detailed information about them can be found at references [46,47].

$$
\rho \frac{\partial h}{\partial t} + \rho \nu \cdot \nabla h = \nabla \cdot (k \nabla T) + \rho \dot{R}\_{q\prime} \tag{6}
$$

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \nu) = 0 \tag{7}$$

$$
\rho\_0 \left( \frac{\partial \upsilon}{\partial t} + \upsilon \cdot \nabla \upsilon \right) = -\nabla \mathfrak{P} + \mu\_l \nabla^2 \upsilon - \rho\_0 \mathfrak{g} \mathcal{J}\_T (T - T\_0) \tag{8}
$$

### 2.1.2. Model Reduction

A preliminary detailed 3D model has been developed. As can be observed in Figure 3a, it represents the complete geometry of the cast parts. The two similar, but different cast parts, the feeding systems formed by two risers, one for each cast part, the filling system and the mold. This 3D model is formed by 1,148,327 tetrahedral elements and 199,031 nodes. In this first attempt the complete simulation has been performed, that is, the part filling and the cooling and solidification process together with the microstructural calculations. The time step has taken values between 1 <sup>×</sup> <sup>10</sup>−<sup>3</sup> (during filling) and 50 s (during the last cooling) depending on the calculation evolution. The resolution of all the equations, shown in last section, for a transient analysis from pouring temperature (1410 ◦C) down to values below the eutectoid transformation temperature (738 ◦C), with the mentioned mesh and time steps, requires a high CPU time. In this case, the simulation has required about 5 h using a computer with 64 Gb of RAM using 8 cores Intel-Xeon CPU E5-2630 v3 @ 2.4 GHz (Dell Inc., Round Rock, TX, USA).

**Figure 3.** (**a**) Mesh used in the 3D simulation model; (**b**) mesh used in the pseudo-2D model.

As it has been previously introduced, the model adjustment is an iterative procedure where the simulation model must be solved each time to get the predictions needed to calculate the fitness function. Therefore, it is highly advised to employ simulation models whose calculation times are as reduced as possible.

For this reason, due to the long calculation times of this model, a model reduction strategy was followed to reduce the high CPU time for the fitness function calculation. In first place the mold filling was suppressed from the calculation. As in this case the mold used was manufactured in resin bonded sand, the temperature drop during filling is not very pronounced. Therefore, the assumption of a mold completely filled of alloy at uniform temperature at the initial step of the simulation, is admissible. In second place the 3D geometry was replaced by a slice of the central area where the tensile test bar was extracted (see Figure 3b and Figure 5). This pseudo-2D geometry assumes that the heat flows mainly across the mold walls.

This model reduction presents some differences in the cooling rate (see Figure 4) which implies an error of 3.6% in the tensile strength prediction (433.64 MPa in the full 3D model and 418.14 MPa in the 2D model). However, the CPU time required to solve this 2D model is of only 29 s (620 times less). Therefore, regarding the CPU time saving, the error level is considered acceptable.

### *2.2. Fitness Function*

The fitness function—also called cost function or objective function—is the mathematical representation of the aspect under evaluation which will be minimized. Therefore, it is the function that represents the difference between the hypothesis (model predictions) and the real results (experimental measurements), that is, the error. The expression used in this case, is the relative error between the predicted tensile strength (σ*predicted*) and the measured tensile strength (σ*re f erence*), following Equation (9). Note that absolute value of the difference is used.

$$fitness = \frac{\left| \sigma\_{predicted} - \sigma\_{reference} \right|}{\sigma\_{reference}} \,\tag{9}$$

**Figure 4.** Cooling curve comparison between the pseudo-2D model and the 3D model.

### *2.3. Variables Selection*

The variables selection that is the parameters of the model that will be modified by the algorithm to adjust the model, was selected in a supervised mode based on the expert's knowledge about the physics involved.

It is advised to select those variables which have more influence in the fitness function. Moreover, among them, those whose values are more uncertain. As it has been previously introduced, the tensile strength prediction is mainly related with the chemical composition of the alloy, the melt treatment and the cooling rate.

In this case the influence of the chemical composition is managed by the thermodynamic database, which calculates the material properties and the different phases that appears during the solidification. CALPHAD is a well-known methodology widely used in computational materials engineering, therefore, although new updated databases are continuously developed, we can assume its predictions as reliable.

The melt treatment can be modeled by means of the parameters that control the nucleation of the austenite: the maximum grain nuclei density *nmax*, the undercooling standard deviation Δ*T*<sup>δ</sup> and the average undercooling Δ*Tn*; those that control the nucleation of the graphite nodules: *Ae* and *n*; and, finally, the parameter that controls the graphite nodules growth, the eutectic growth coefficient μ*e*. All these parameters are very difficult to estimate, so they are good candidates to be included in the parameters group to be modified by the optimization method.

The cooling rate is mainly determined by the material properties of the alloy and the mold and for the heat transfer coefficient (HTC) between them. As said, the material properties of the alloy are calculated by the thermodynamic database, so they can be considered reliable. Resin bonded sand molds are manufactured mixing together clean and dry sand, with a binder and a catalyst in a continuous mixer. Different granulometries and types of sand can be used together with different binders and catalysts from several manufacturers, which can be also mixed in different proportions, so the resin bonded sand molds may present differences in their properties. Therefore, they can be also included in the parameters group to be modified by the optimization method. The heat transfer coefficient (HTC) between the alloy and the mold is difficult to estimate. When the metal enters in the mold, the macroscopic contact between the alloy and the mold is good because of the conformance of the molten metal. In this period the HTC presents high values, typically between 100 and 3000 W/m2 K depending on different factors (surface roughness, pressure, wettability, etc.). Once a solidified layer with enough strength has grown, the contraction suffered by the alloy, together with the mold distortion caused by the mold heating, are likely to start a gap between them. In this phase, the HTC will mainly depend on the gas conductivity and the gap size, and its values may drop to values 10 times

lower (10–100 W/m2 K). In the ductile iron case, the expansion that takes place in the alloy during the graphite precipitation tends to reduce that gap but makes even more difficult the estimation of the HTC value. Therefore, the HTC is also a good candidate to be modified during the adjustment. However, in order to facilitate the implementation of the algorithm designed to perform the model fitting, a constant value was used for the HTC.

### *2.4. Optimization Methods*

In the model fitting procedure followed, the parameters modification is guided by the optimization method. The performance of four different methods was compared. Two of them are deterministic methods, the compass search and the NEWUOA. The other two are heuristic methods, the particle swarm optimization (PSO) and the genetic algorithm (GA).

Compass search is a deterministic local optimization algorithm. It uses a simple and slow, but very sure procedure. The algorithm advances varying one of the parameters that form the solution vector at a time by steps of the same magnitude and when no such increase or decrease in any one parameter further improves the fit to the experimental data, they halve the step size and the process is repeated until the steps are deemed sufficiently small. For example, in a minimization problem with only two variables, the algorithm can be summarized as follows: Try steps to the East, West, North and South. If one of these steps yields a reduction in the function, the improved point becomes the new, iterate. If none of these steps yields improvement, try again with steps half as long. More details about this algorithm can be found in [48].

NEWUOA is also a deterministic local optimization algorithm. It seeks the least value of a function *F* without constraints and without derivatives. A quadratic approximation *Q* of the *F* function is used to obtain a fast convergence ratio, using a trust region strategy (see Figure 5). In this type of strategy, a model which represents a behavior similar to the objective function in a point *xk* is generated. The search of the optimum is restricted to the surroundings to that point, which is called trust region. As the optimization progress, successive approximations of *Q* are used. More details about this algorithm can be found in [49].

**Figure 5.** Quadratic approximation of a function

Particle swarm optimization (PSO) is a heuristic bioinspired algorithm for global optimization. A number of particles are placed in the search space of some problem or function and each evaluates the objective function at its current location. Each particle then determines its movement through the search space by combining some aspects of the history of its own current and best (best-fitness) locations with those of one or more members of the swarm, with some random perturbations. The next iteration takes place after all particles have been moved. Each individual in the particle swarm is composed of three D-dimensional vectors, where D is the dimensionality of the search space. The first vector is the current position *xi*, which can be considered as a set of coordinates describing a point in space that is evaluated as a problem solution. Second vector is formed by the coordinates of previous best position *pi*. The value of the best global position so far is stored in a variable for comparison on

later iterations. The third vector is the velocity *vi*, which can effectively be seen as a step size. In this way, new points are chosen by adding *vi* coordinates to *xi* and the algorithm operates by adjusting *vi*. Each particle communicates with some other particles and is affected by the best point found by any member of its neighborhood. More details about this type of algorithm can be found in reference [50].

Genetic algorithms (GA) are heuristic search algorithms of general purpose inspired in the genetic processes that take place in biologic organisms and in the principles of the natural evolution of the populations. The basic idea consists in the generation of a random population of individuals or chromosomes and the advance to better individuals by means of the application of genetic operators, which are based in the genetic processes that takes place in the nature. Each of the individuals represents a possible solution to the problem, so the population represents the solutions search space. During the successive iterations, called generations, the individuals of the population are evaluated depending on their adaptation as solution. To do that, each one has associated a fitness value which measures how of good or bad is the possible solution that represents. The new population of individuals is created by selection mechanisms and by crossover and mutation operators. In this way, with the successive iterations or generations the algorithm converges to better solutions until reaching the optimum or the maximum number of iterations. More information about this type of algorithms can be found in reference [51].

### **3. Case Study**

The case study is a part with a stepped geometry, cast in a sand mold bonded with resin where two similar parts are cast simultaneously. The reference value used for the model fitting is the tensile strength measured experimentally. Three tensile test bars were extracted from the central area of the 35-mm high step of the longer part, shown in Figure 6. The mean tensile strength obtained is equal to 552 MPa and the standard deviation equal to 4 MPa. The material is ductile iron, also called spheroidal gray iron (SGI), whose slightly hypereutectic composition is shown in Table 1, together with the tensile strength value measured experimentally. The mold was manufactured with silica sand mixed with a phenolic urethane resin (Pep-Set) and the corresponding catalyst. Figure 7 shows the manufactured cast part.

**Figure 6.** Top and bottom view of the part geometry.


**Table 1.** Ductile iron chemical composition.

**Figure 7.** Manufactured cast part.

One important aspect in the performance of optimization algorithms is the dimensionality problem—that is, the number of variables that can be modified by the algorithm. Typically, classical optimization algorithms, as compass search or NEWUOA in this case, performs better in low dimensional problems, but when the number of variables increases it may require too long calculation times. Instead, the advantage of the heuristic algorithms, as PSO and GA in this case, takes place when the problem dimensionality is high. For this reason, the performance of the four optimization algorithms was evaluated for different numbers of variables included in the model fitting. The cases of 2, 3, 6 and 10 variables were studied. The corresponding variables considered in each case are detailed in Table 2. Table 3 summarize the default values used for these variables together with their ranges.


**Table 2.** Variables included in the model fitting cases.


**Table 3.** Default values and ranges of the variables included in model fitting.

The optimization algorithms were used as they are implemented in the PyGMO library previously mentioned, using quite standard values for their internal parameters. Compass search was defined with a start range equal to 0.1, a stop range equal to 0.001 and a reduction coefficient equal to 0.5. In the case of NEWUOA, no additional parameters were setup. In order to compare both algorithms in similar conditions, the same seed were used in both cases.

Compass search and NEWUOA are local optimization algorithms, for this reason the selection of the initial guess may have influence in the final solution of the algorithm. For this reason, two different configurations were used for both. In the first configuration the algorithms are initialized with only one random guess. In the second one, they are initialized with 20 random guesses. In this last case, after the first iteration with those 20 individuals, the best of them is selected and the algorithms continue iterating from that point.

PSO and GA are stochastic population-based algorithms. In both cases a population size of 20 individuals were used. The canonical version of PSO based on the constriction coefficient velocity update formula was used. In addition, in this case quite standard values were used: constriction factor equal to 0.7298, social component equal to 2.05, a cognitive component equal to 2.05, maximum allowed particle velocities normalized respect to the bounds width equal to 0.5 and 4 individuals considered as neighbors. In the GA case, an exponential crossover operator with a probability equal to 0.90, a polynomial mutation operator with a probability equal to 0.02 and a tournament selection operator were used.

More details about the algorithm parameters can be found in [39].

### **4. Results and Discussion**

The simulation model solved for the default parameters shown in Table 3, predicts a tensile strength value equal to 418.14 MPa, which corresponds to a relative error equal to 24% (see Equation (10)). The objective is to evaluate the ability of the algorithms to reduce that error to values below 0.1%, modifying the variables shown in Table 2.

$$\text{Error} = \frac{|418.14 - 552|}{552} = 0.2425,\tag{10}$$

Figure 8 collects the convergence evolution of each algorithm for the different cases studied. The behavior of the deterministic algorithms was uneven. The performance of the compass search was very good, reducing the error to values below the 0.1% in all cases with a reasonable number of calculations. The use of 20 initial random points has improved the algorithm performance reducing the number of calculations, especially in the 10 variables case where the difference is clear.

**Figure 8.** Convergence evolution.

NEWUOA algorithm instead, failed in reaching the error target in several cases. Using an initialization based on only one point, the algorithm has failed in the 6 and 10 variables cases with error values equal to 22.62% and 5.77%, respectively. The use of 20 random initial points clearly has improved the behavior of the algorithm making possible to reach the target in the 6 variables case and to reduce the error value to 2.07% in the 10 variables case.

The heuristic algorithms, GA and PSO, have succeeded in reaching the error target in all the cases, with the exception of the GA in the 3 variables case, where the minimum error value was equal to 0.124%. The number of calculations needed to reach the target error was lower for PSO than for GA, except for the 2 variables case.

Therefore, the best options were the compass search algorithm with 20 random initial points, followed by the PSO algorithm.

A priori, it can be thought that as higher the number of variables, larger the search domain and higher the difficulty of the algorithm to reach the target. However, a higher number of variables may also imply a bigger number of solutions that fulfil the target value, making easy to find one of them. On the other side, a reduced number of variables may be too restrictive limiting the number of possible solutions and making difficult to reach the target. The results obtained are not clear in this sense. There is a tendency in the deterministic algorithms to find more difficulties in cases with higher number of variables in contrast with the heuristics, which have found more difficulties in cases with a lower number of variables. However, the obtained results are not conclusive in this sense.

Relating to the values assigned to the variables, as it was explained in the introduction, it does not exist a unique solution of the problem. Therefore, the model may be fitted assigning different combinations of the values assigned to the variables. The risk of this fact is that sometimes, the values assigned to the variables by the algorithm may cause a certain loss in the physical sense of the model although the results are accurate from the mathematical point of view. However, this is a common difficulty of this type of approximation that has not been solved by now.

Figure 9 shows the distribution of the values assigned to the different variables in boxplot form. Black points represent the values assigned for each variable jittered, the blue line corresponds to the value assigned for each variable in the initial model and red lines to the bounds fixed for each of them (In the case of NEWUOA no bounds are used). As can be observed the distribution is quite different depending on the adjusted variable.

**Figure 9.** Distribution of the values assigned to the different variables in boxplot form.

Figure 10 collects the values assigned grouped by the number of variables included in the fitting and colored by the algorithm used in each case. Colored bars represent the values assigned to the variables by the different optimization algorithms, the dotted red lines represent the bounds fixed for the corresponding variable and the dotted blue line the value initially assigned.

As expected, the higher the number of variables included in the adjustment, the higher the variability in the values assigned. This can be observed in the cases of the maximum grain nuclei density, the undercooling standard deviation and the heat transfer coefficient. Otherwise, the variability of the values assigned to the variable seems also related with the algorithm used. For example, compass search presents a bigger tendency to adjust the variable to higher values. However, as the "true" value of each variable is actually unknown, it is difficult to evaluate if any of the algorithms has a better performance in assigning more reliable values to the variables.

**Figure 10.** Values assigned to variables grouped by the number of variables included in the fitting and colored by the algorithm.

### **5. Conclusions**

Only compass search and PSO have succeeded in reaching the error target in the four cases studied (2, 3, 6 and 10 variables). The NEWUOA and GA algorithms have failed in some case.

In the case of the deterministic algorithms, compass search and NEWUOA, the use of a multiple random initial guess was clearly beneficious.

Compass search with an initial 20 random points has presented the best performance, reaching the error target in all cases with the lower number of calculations.

The known problem related with the absence of a unique solution was found, obtaining different combinations of the values assigned to the variables. However, as the 'true' value of each variable is actually unknown, it is difficult to evaluate if any of the algorithms has a better performance in assigning more reliable values to the variables.

Future works. The authors believe that the research line focused in the fitting of metal casting models by means of optimization algorithms is very interesting, as the advances in this area may contribute to reduce the time needed to adjust the metal casting models improving the confidence in their predictions and contributing to reduce tests at laboratory scale. In fact, some of the future works planned are:



**Author Contributions:** Conceptualization, E.A.; data curation, A.O. and E.V.; formal analysis, E.A.; investigation, E.A., A.M., A.O. and E.V.; methodology, E.A.; resources, A.M.; software, E.A.; validation, A.M.; writing—original draft, E.A.; writing—review & editing, E.A., A.M. and I.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Basque Government under the ELKARTEK Program (ARGIA Project, ELKARTEK KK-2019/00068) and by the HAZITEK Program (CASTMART Project, HAZITEK ZL-2019/00562).

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

### **Nomenclature**


### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Analysis of Melt-Pool Behaviors during Selective Laser Melting of AISI 304 Stainless-Steel Composites**

**Daniyal Abolhasani 1,2, S. M. Hossein Seyedkashi 2, Namhyun Kang 3, Yang Jin Kim 1, Young Yun Woo <sup>1</sup> and Young Hoon Moon 1,\***


Received: 29 June 2019; Accepted: 6 August 2019; Published: 8 August 2019

**Abstract:** The melt-pool behaviors during selective laser melting (SLM) of Al2O3-reinforced and a eutectic mixture of Al2O3-ZrO2-reinforced AISI 304 stainless-steel composites were numerically analyzed and experimentally validated. The thermal analysis results show that the geometry of the melt pool is significantly dependent on reinforcing particles, owing to the variations in the melting point and the thermal conductivity of the powder mixture. With the use of a eutectic mixture of Al2O3-ZrO2 instead of an Al2O3 reinforcing particle, the maximum temperature of the melt pool was increased. Meanwhile, a negligible corresponding relationship was observed between the cooling rate of both reinforcements. Therefore, it was identified that the liquid lifetime of the melt pool has the effect on the melting behavior, rather than the cooling rate, and the liquid lifetime increases with the eutectic ratio of Al2O3-ZrO2 reinforcement. The temperature gradient at the top surface reduces with the use of an Al2O3-ZrO2 reinforcement particle due to the wider melt pool. Inversely, the temperature gradient in the thickness direction increases with the use of an Al2O3-ZrO2 reinforcement particle. The results of melt-pool behaviors will provide a deep understanding of the effect of reinforcing particles on the dimensional accuracies and properties of fabricated AISI 304 stainless-steel composites.

**Keywords:** selective laser melting; additive manufacturing; SLM; FEM; Al2O3; reinforced; Al2O3-ZrO2; 304; stainless; composite

### **1. Introduction**

Austenitic stainless steels can find more applications if their strength is improved. Hence, some studies have investigated the reinforcement of stainless steel matrices via the selective laser melting (SLM) process [1,2]. Reinforcing can be achieved through the use of Al2O3 to produce a metal matrix composite that exhibits high mechanical properties, in addition to high corrosion and wear resistance. However, Al2O3 particle has very high melting temperatures and very low thermal conductivities, which result in a very high thermal gradient during the process, causing high local stresses and crack formation [3]. The use of a eutectic mixture of Al2O3 and ZrO2 powders decreases the melting temperature of Al2O3 from 2313 K to about 2133 K [3,4]. If an optimum volume can be identified for the samples reinforced with a eutectic mixture of Al2O3-ZrO2, it can be a good alternative to strengthen the AISI 304 stainless steel.

To provide a deep understanding in relation to composites reinforced with alumina particles, this paper presents an innovative numerical study dealing with melt-pool behaviors during selective laser melting (SLM) of Al2O3-reinforced and a eutectic mixture of Al2O3-ZrO2-reinforced AISI 304 stainless-steel composites. As is well known, the dimensional accuracies and properties of fabricated composite parts are significantly dependent on the melt-pool behaviors [2–4].

SLM is a manufacturing technique to construct three-dimensional parts in which a high-power density laser is used to melt and fuse metallic powders [4–7]. Compared with other traditional techniques, laser processing typically does not require mechanical tooling and therefore exhibits high flexibility [8–12]. In SLM, the energy density should be sufficient to ensure both the melting of the powder and the bonding of the underlying substrate. If the bonds between the scan tracks and layers are weak, defects such as cracks may be generated [13–16]. The trapped gas in the melt pool can increase the number of pores in SLM-fabricated parts. Additionally, the reduction in solubility of the element during the rapid melting and cooling process can cause detrimental defects [17]. Schleifenbaum et al. [18] found that with the increase in the laser power, the rate of material evaporation and the incidence of spattering increased. These defects should be prevented by identifying appropriate process parameters of the SLM process.

The prediction of the temperature evolution in SLM is commonly performed using the finite-element (FE) method or statistical approaches. The temperature and stress fields in samples of 90W-7Ni-3F and 316L stainless steels were predicted by Zhang et al. [19] and Hussein et al. [20]. Dai and Shaw [21] simulated the transient temperature, as well as the thermal and residual stress fields. Kundakcioglu et al. [22] performed transient thermal modeling of laser-based Additive Manufacturing (AM) for 3D freeform structures.

In the case of composites, AlMangour et al. [23] presented a simulation model for predicting the temperature evolution of the melt pool of TiC/SUS316L samples. Shi et al. [24] investigated the effects of the laser power and scan speed on the melting and solidification mechanisms during SLM of TiC/Inconel718 via a simulation approach. Li et al. [1] investigated the microstructural and mechanical properties of Al2O3/316L stainless-steel metal-matrix composites (MMCs) fabricated via SLM through experimental and numerical methods.

The aforementioned studies focused on a single reinforcement particle in the case of SLM of a metal matrix composite. They did not consider an SLM composite with binary-phase reinforcement, such as a eutectic mixture. The thermal cycle in the SLM process increases the complexity of the material melting and thermal behavior, if new combinations of the reinforcements are applied [25]. The melting temperature of the powder mixture is an important microstructural feature in the thermal process, as it affects the final temperature of the surface during the heating or cooling of the melt pool. No detailed analysis involving the 3D FE modeling of a selective-laser-melted part, including the metal matrix and binary-phase reinforcement and its relationship with the geometric features of the melt pool during the process, has been performed.

In this study, a eutectic mixture of Al2O3 and ZrO2 powder particles was added to AISI 304 stainless steel as the metal-matrix media in a set of numerical simulation runs, incorporating experimental runs as validation tests. The use of a eutectic mixture of Al2O3 and ZrO2 powders reduces the melting temperature of Al2O3 particles. Hence, for a better understanding of the results, both sets of parts modeled with various weight percentages of Al2O3 and Al2O3-ZrO2 particles were considered.

For this investigation, the melt-pool dimensions and the thermal evolution of AISI 304-Al2O3 and AISI 304-Al2O3-ZrO2 SLM composites were simulated, and the predicted thermal results were described.

### **2. Materials and SLM System**

In this study, the eutectic mixture was prepared using powders containing 58.5 wt% Al2O3 and 41.5 wt% ZrO2, corresponding to the eutectic point of the Al2O3-ZrO2 binary phase diagram [4], with the reinforcement content within the AISI 304 stainless-steel powder as the metal matrix. An experiment was performed using an SLM system with a continuous-wave, ytterbium fiber laser (IPG YLR-200, IPG Photonics, Burbach, Germany), as shown in Figure 1. The maximum available power of 200 W at 6 A was used, and the laser scanning was controlled using a scanner (hurrySCAN®20, SCANLAB, Puchheim, Germany). Argon gas was used as a shielding gas to prevent oxidation in the melt pool. The layering bar spread the powder, and the build cylinder moved vertically to control the powder-bed height. During the numerical and experimental runs, a 30-μm-thick layer of powder

mixture was used. The laser spot diameter was fixed at 80 μm. The laser power and scanning velocity were selected after conducting a large number of preliminary experiments associated with the single-line formation tests.

**Figure 1.** Selective laser melting (SLM) system: (**a**) Experimental setup; (**b**) Schematic illustration.

The laser power and scan speed were 200 W and 732 mm/s, respectively. For comprehensive comparison, various weight percentages of Al2O3 and the eutectic mixture of Al2O3-ZrO2 were employed in the reinforcing experiments, as shown in Table 1.


**Table 1.** Numerical and experimental runs.

### **3. FE Modeling**

### *3.1. Physical Description of Model*

A 3D model was developed, and the ABAQUSTM commercial software (version 6-14, Dassault Systems, Providence, RI, USA) was utilized to simulate volumetric laser energy deposition with a Gaussian distribution. A schematic of the SLM process, which shows the interactions between the laser beam and powder bed, as well as the melting-pool dimensions and solidification, is presented in Figure 2a. A square composite powder layer with dimensions of 7 mm × 7 mm × 0.030 mm was placed on a stainless-steel substrate. The powder bed was meshed with eight-node linear hexahedral elements with a size of 70 <sup>×</sup> <sup>10</sup>−<sup>6</sup> m that were distributed uniformly throughout the powder-bed model, as shown in Figure 2b. The sweep method was used to mesh the substrate and accurately capture the flux distribution of the moving laser beam.

**Figure 2.** (**a**) Schematic showing the thermal behavior of the powder bed under laser irradiation. (**b**) 3D finite-element (FE) model. (**c**) Schematic showing the transfer of laser radiation into the powder bed.

### *3.2. Initial Governing Equation*

In the powder-bed additive layer manufacturing system with a moving laser source, the heat-transport equation for 3D Cartesian coordinate systems is expressed as follows [22]:

$$\rho \cdot \mathbf{C}\_{\mathcal{V}} \frac{\partial T}{\partial t} = \frac{\partial}{\partial \mathbf{x}} (\mathbf{k} \cdot \frac{\partial T}{\partial \mathbf{x}}) + \frac{\partial}{\partial y} (\mathbf{k} \cdot \frac{\partial T}{\partial y}) + \frac{\partial}{\partial \mathbf{z}} (\mathbf{k} \cdot \frac{\partial T}{\partial \mathbf{z}}) + \rho \cdot \mathbf{C}\_{\mathcal{V}} \mathbf{v} \cdot \frac{\partial T}{\partial \mathbf{x}} + Q\_{\prime} \tag{1}$$

where *<sup>T</sup>* is the temperature (K), <sup>ρ</sup> is the density (kg/m3), *Cp* is the specific heat (J/(kg·K)), *<sup>k</sup>* is the thermal conductivity (W/m K), and ∇ is the differential or gradient operator. *Q* is the rate of internal energy conversion per unit volume (referred to as the source term, W/m3).

The corresponding boundary condition was [20]:

$$-\frac{\partial T}{\partial z} = \frac{h}{k}(T - T\_0) \ + \frac{\varepsilon}{k}(T^4 - T\_0^4),\tag{2}$$

where *T* is the surface temperature of the powder bed, *z* is the axis perpendicular to the powder surface, *h* = 200 [26] is the forced-convection heat-transfer coefficient (W/m2 K) due to argon gas, ε = 0.8 is the emissivity of the heated surface, and <sup>б</sup> <sup>=</sup> 5.6703 <sup>×</sup> <sup>10</sup>−<sup>8</sup> <sup>W</sup>/m2K4 is the Stefan–Boltzmann constant.

### *3.3. Heat-Source Model*

The travel distance of the laser beam into the powder media (*z*), which is schematically presented in Figure 2c, can be modeled using a volumetric heat source (*Q*(*x*, *y*, *z*)) [27]:

$$Q(\mathbf{x}, \mathbf{y}, z) = (1 - R) \, n \, A \, \frac{P}{\pi r\_0^2} \exp(\frac{-n(\mathbf{x}^2 + \mathbf{y}^2)}{r\_0^2}) \exp(-\frac{1}{d} z),\tag{3}$$

where *R* is the surface reflectivity of the powder mixture, and ݅ is a shape parameter of the heat-flux distribution, ݅ = 2. *A* is the laser-absorption coefficient of the irradiated surface and it was modified during numerical validation. *P* is the laser power (W), *r*<sup>0</sup> is the laser-beam radius (m), *x* and *y* are the Cartesian coordinates at the surface, and *z* is the Cartesian coordinate perpendicular to the powder bed. *d* is the penetration depth equal to 40 μm for the SLM process [27].

The surface reflectivity of the powder mixture (*R*) is expressed as follows [25,28]:

$$\mathcal{R} = \frac{\sum\_{i=1}^{n} \frac{\mathbb{R}\_{(i)} \: \mathbb{R}\_{(i)} \le \mathbf{w}\_{(i)}}{\rho\_{(i)}}}{\sum\_{i=1}^{n} \frac{\mathbb{R}\_{(i)} \: \mathbf{w}\_{(i)}}{\rho\_{(i)}}},\tag{4}$$

Here, in accordance with the three components utilized in this study (AISI 304, Al2O3, and ZrO2), *n* is 3; *i* is a subscript referring to a specific mixture component; *R*(*i*) is the surface reflectance of component *i*; β(*i*) is the extinction coefficient of component *i*, ρ(*i*) is the density of component *i*, and w(*i*) is the weight percentage of component *i* in the mixture.

The extinction coefficient β is defined as [29]:

$$\beta \, = \, \frac{3 \, (1 - \phi)}{2 \, \phi \, D},\tag{5}$$

where φ is the porosity of the bed, and *D* is the average particle diameter of the powder. The powder density was calculated using the following equation [26]:

ρpowder = (1 - φ) ρbulk, (6)

The laser-absorption coefficient of a powder bed with a Gaussian distribution (*A*powder) was estimated using the absorption coefficient of the bulk material (*A*bulk) [30]:

$$A\_{pvmler} = 0.0413 + 2.89 A\_{bulk} - 5.36 A^2 \, \_{bulk} + 4.50 A^3 \, \_{bulk} \tag{7}$$

The values of *A* in Equation (3) were calculated according to the linear rule of mixtures [31]:

$$A = w\_1 \cdot A\_1 + w\_2 \cdot A\_2 + w\_3 \cdot A\_3 \tag{8}$$

where *w*1, *w*2, and *w*<sup>3</sup> represent the volumetric contents of the three components in the mixture, and *A*1, *<sup>A</sup>*2, and *<sup>A</sup>*<sup>3</sup> represent the absorption coefficients of the components ( <sup>3</sup> *i*=1 *wi* = 1, 0 < *A* < 1). The final values of *A* applied in Equation (3) were obtained through numerical validations, as discussed in Section 4.1.

The parameters presented in Table 2 were obtained from the material supplier and previous works [26,32], as well as an online source. Then, the corresponding values were calculated using Equations (4)–(8).

**Table 2.** Parameters used for obtaining the values of the surface reflectivity (*R*) and absorption coefficient (*A*) of the powder mixture.


### *3.4. Physical Properties of Materials*

The values of *k* are expressed by the following equations [29,31,33], where *T* is the temperature and the superscripts *s* and *m* correspond to "solidus" and "melting", respectively.

$$k = \begin{cases} \text{ for AISI 304, } T < T\_s: k\_{\text{pound}} = k\_{\text{bulk}} (1 - \phi)^n, n = 4\\ \text{ for AI}\_2 \text{O}\_3 \text{ and } \text{ZrO}\_2, \text{ } T < T\_s: \\ \frac{k\_{\text{pound}}}{k\_{\text{dyn}}} = 1 - \sqrt{1 - \Phi} \left( 1 + \frac{\Phi\_k \frac{k\_{\text{pd}}}{k\_{\text{dd}}}}{k\_{\text{dm}}} \right) + \sqrt{1 - \Phi} \left( \frac{2}{1 - \frac{k\_{\text{dm}}}{k\_{\text{dd}}}} \ln \left( \frac{k\_{\text{DL}}}{k\_{\text{dd}}} \right) - 1 \right) + \frac{k\_{\text{md}}}{k\_{\text{dd}}} \text{y}, \\ \text{for all components, } T\_s < T < T\_m: \\ k = \begin{pmatrix} k\_0 \ln \left( \frac{(T\_m) - k\_0 \text{ quadr}}{T\_m - T\_s} \right) (T - T\_s) + k\_0 \text{ pound} & (T\_s) \right) \times 10^{-3} \\ \text{for mixture, } T < T\_s, T\_s < T < T\_m: k\_{\text{mickure}} = w\_1 k\_1 + w\_2 k\_2 + w\_3 k\_3 \end{cases} (9)$$
 
$$\text{for mixture, } T \ge T\_m:$$
 
$$k\_{\text{mickure}} = k\_1 \left( \frac{k\_2 \left( 1 + 2w\_2 \right) - k\_1 (2w\_2 - 1)}{k\_2 (1 + w\_2) - k\_1 (2w\_2)} \right), \text{ } w\_2 < w\_1, k\_1 = k\_1 \text{ bulk}, k\_2 = k\_2 \text{ hulk}$$

Here, *Katm* is the thermal conductivity of the ambient atmosphere and is 0.018 for argon gas; φ is the porosity of the powder bed; *k*<sup>0</sup> is the thermal conductivity at the ambient temperature; *Krad* is the thermal conductivity due to the radiation among the particles in the powder bed and is evaluated as in [29]:

$$K\_{rad} = 4F \cdot 6T^3 \cdot D\_\prime \tag{10}$$

where *F* = 1/3 is the view factor; б is the Stefan–Boltzmann constant; *T* is the temperature; and *D* is the average particle diameter of the powder.

$$\mathbf{C}\_{p} = \begin{cases} \text{for all components,} & T\_s < T < T\_m \ : \mathbf{C}\_p = \mathbf{C}\_{p0} + \left( \frac{1}{(T\_m - T\_s)\sqrt{\pi}} e^{-\left(\frac{T}{T\_m} - \frac{T\_m}{T\_s}\right)} \right) \times \mathbf{L} \\ & \text{for mixture,} \ \mathbf{T} < T\_s, \ T\_s < T < T\_m \ : \mathbf{C}\_{p\text{mixture}} = w\_1 \mathbf{C}\_{p1} + w\_2 \mathbf{C}\_{p2} + w\_3 \mathbf{C}\_{p3} \\ & \text{for mixture,} \ \mathbf{T} \ge T\_m \ \mathbf{C}\_p \text{ of mixture was not available} \end{cases} \tag{11}$$

Here, *Cp*<sup>0</sup> is the specific heat capacity at the ambient temperature, *L* is the latent heat of melting, and *Cp*1, *Cp*2, and *Cp*<sup>3</sup> are the specific heat capacities of the three components [33].

$$\rho = \begin{cases} \text{ for AISI 304, } T < T\_s: \rho = \rho\_1 \\ \text{ for AI}\_2\text{O}\_3 \text{ and } \text{ZrO}\_2, \; T < T\_s: \; \rho = \rho\_0 \\ \text{for all components, } T\_s < T < T\_m: \; \rho = (1 - \frac{T - T\_s}{T\_m - T\_s})\rho\_0 + (\frac{T - T\_s}{T\_m - T\_s})\rho\_0 \text{ bulk} \\ \text{ for mixture, } \; T < T\_s, \; T\_s < T < T\_m: \; \rho\_{\text{mixture}} = w\_1\rho\_1 + w\_2\rho\_2 + w\_3\rho\_3 \\ \text{ for mixture, } \; T \ge T\_m: \; \rho\_{\text{mixture}} = w\_1\rho\_{\text{bulk}} \; 1 + w\_2\rho\_{\text{bulk}} + w\_3\rho\_{\text{bulk}} \end{cases} \tag{12}$$

Here, ρ1, ρ<sup>2</sup> and ρ<sup>3</sup> are the densities of the 304 stainless-steel, Al2O3 and ZrO2 powder components, respectively. ρ<sup>0</sup> is the density of the powder component at the ambient temperature [31,33,34]. Note that ρ*powder* = ρ*bulk* when φ = 0 at the melting point.

The phase change was included in the model by using the definition of the latent heat of melting. The *Ts* of the component with the lowest value (1615 K for AISI 304 stainless steel) and the *Tm* of the component with the highest value were used for the mixture when simulating composite parts. The solution-dependent state variables of the three fields (powder, solid, and liquid) were based on the density of the mixture (ρ) using a User Defined Field (USDFLD) subroutine.

### **4. Results and Discussion**

To analyze melt-pool behaviors, cross-sectional views of the melt pools obtained from experimental and numerical models are compared first. Additional simulation results are then presented to exhibit the temperature effects on the liquid lifetime of the melt pools. The effect of temperature gradient on microstructural features is also discussed.

### *4.1. Numerical Validation and Melt-Pool Characterization*

To reduce the simulation time, the single track of conditions presented in Table 1 was modeled. The numerical model was validated using the calculated dimensions of the melt pool and compared with the experimental results of the single-line formation test. Image processing was performed on the cross-sectional optical micrograph to evaluate the melt-pool morphology, as shown in Figure 3.

**Figure 3.** Schematic showing the dimensions of the melt pool. W and D represent the width (μm) and depth (μm), respectively.

The melt-pool boundary was distinguished by the melting-temperature line of the numerical-simulation thermal field. Thus, the width and depth of the numerical melt pool were compared with the average width and depth of the experimental melt pool.

Through numerical validation, the best results for the melt-pool dimensions were obtained by changing the absorption coefficient (*A*) of the mixture. The experimental and numerical cross-sections of three samples—pure AISI 304, 3 wt% Al2O3, and 3 wt% Al2O3-ZrO2 (eutectic mixture)—are shown in Figure 4a–f. The black lines in the cross-section of the simulated conditions represent the melt-pool boundaries, which are the melting temperatures, i.e., 1670 K for pure AISI 304, 2313 K for Al2O3 reinforced, and 2133 K for the eutectic mixture reinforced samples [4].

**Figure 4.** *Cont.*

**Figure 4.** Cross-sectional views of the experimental and numerical models: (**a**,**b**) pure AISI 304; (**c**,**d**) 3 wt% Al2O3; (**e**,**f**) 3 wt% Al2O3-ZrO2. (NT11: Nodal Temperature).

The experimental and numerical melt pools exhibited similar shapes. The width and depth of the melt pools under the conditions listed in Table 1 are plotted in Figure 5a,b, respectively. The two reinforced composites exhibited slight differences in width, and the width of the melt pool was smaller for the pure AISI 304 sample. As thermal conduction is the most influential factor for the melt pool [29], the results are attributed to the lower thermal conductivity in the reinforced samples compared with the pure AISI 304, as shown in Figure 5c at 800 K. The Al2O3-reinforced composite exhibited the smallest melt-pool depth among the samples, as it had the highest melting point, making it difficult for the molten powders to penetrate deep inside the melting pool. The larger melt-pool width and depth for the eutectic-reinforced sample compared with the Al2O3-reinforced sample at each weight percentage are attributed to the reduction in the melting point, along with the reduction in the thermal conductivity.

**Figure 5.** (**a**) Cross-section width and (**b**) depth of the experimental and numerical models. (**c**) Thermal conductivity (*k*) of the mixtures at 800 K.

### *4.2. Temperature Distribution and Liquid Lifetime*

The variation of the temperature with the time needed for scanning a single path during the process, with different weight percentages of reinforcement contents, is explored. The center point of the top surface shown in Figure 6a is considered for plotting the corresponding temperature-time profiles shown in Figure 6b. Figure 6c,d also exhibit the maximum temperature of the top surface and melt-pool lifetime reaching to 300 K for each condition, respectively. As the eutectic mixture of Al2O3-ZrO2 is replaced with the Al2O3 reinforcing particle, the maximum temperature is increased significantly due to the reduction in the thermal conductivity and latent heat of the mixture. Also according to Figure 6b, concerning the slope of the cooling curve as the cooling rate, this rate decreased especially in the cases with 7 wt% of reinforcement particles. The simulation results show that as the reinforcement content increases, the cooling rate decreases. The main factor affecting the maximum temperature in various weight percentages of a specific reinforcement in the simulations is identified to be the absorption coefficient. The absorption coefficient of the mixtures was calculated using Equation (8) and exhibited a decreasing trend when increasing a reinforcement content, as is shown in Figure 6e. It is important to note that the cooling rate, i.e., the slope of the cooling curve of Al2O3 and Al2O3-ZrO2 systems upon each weight percentage of the reinforcing particle, showed a negligible discrepancy. From the cooling step in simulation runs, the liquid lifetime of the melt pool is distinguished. From Figure 6d, it can be seen that the liquid lifetime increases with the use of a eutectic ratio of Al2O3-ZrO2 reinforcement. Meanwhile, it also rises as the reinforcement content increases. Because the gas atoms are released from the lattice in the heat-affected zone, yielding a longer liquid lifetime in samples with a high weight percentage of reinforcement, e.g., 7 wt%, may result in the formation of a higher amount of gas. It is reported that this phenomenon may cause the formation of porosity in SLM-produced parts [35]. This was explored by multi-line formation tests shown in Figure 7, where the presence of cracks is evident in the reinforced sample with 7 wt% Al2O3-ZrO2, seen in Figure 7f. On the other hand, a short liquid lifetime of the melt pool combined with a lower temperature, observed in the sample with 3 wt% Al2O3, is detrimental for wetting the farther gaps of powders, which generated some inter-gaps, as observed in Figure 7a. This can also be observed in Figure 5b, in which in the sample with 3 wt% Al2O3 is shown, and the depth of the melt pool cannot reach to the powder bed thickness indicated by the dashed line, i.e., 30 μm. This can be expanded to other reinforced samples with Al2O3 particles, where lower temperatures and shorter liquid lifetimes are seen compared with reinforced samples with eutectic mixtures. An appropriate weight percentage of the reinforcement particle of Al2O3 or the eutectic mixture of Al2O3-ZrO2 plays an important role in the SLM process of AISI 304 stainless-steel composites. In the case of 3 wt% Al2O3-ZrO2, the formation of an uneven surface, as observed in Figure 7d, is the reason for which it is not regarded as an optimum condition. A proper temperature to melt the particles with appropriate liquid lifetime and melt-pool depth was achieved using a 3 wt% eutectic mixture. Considering a more moderate condition for liquid lifetime, however, the sample of 5 wt% eutectic mixture shows a reasonable trend, as indicated by Figures 6b–d and 7e. From the above observation and discussion, it is evident that the melt-pool lifetime has an effect on the melting behavior, rather than the cooling rate.

### *4.3. Temperature Gradient*

As SLM is an unsteady-state process, temperature gradients along surface direction and thickness direction are unavoidable. The heating and solidification during the process will alter the temperature gradient. This affects the microstructural features, such as grain morphology, grain size and its size distribution, as well as residual stress accumulation due to a large temperature gradient. In general, the unstable material flow, warpage and delamination between fabricated layers are detrimental results of process-induced stresses [36–38].

**Figure 6.** (**a**) Center point on the top surface chosen for plotting the (**b**) Temperature-time curves. (**c**) Results for the maximum temperatures of the reinforced samples. (**d**) Liquid lifetimes. (**e**) The final estimations of the absorption coefficients.

The slope of the curve of temperature distribution presents the temperature gradient of the material [20]. Figure 8a,b show the effects of reinforcement contents on the temperature distribution along the Y-direction at the top surface and Z-direction, or thickness direction, during SLM of Al2O3/Al2O3-ZrO2, AISI 304 systems. With the use of an Al2O3-ZrO2 reinforcement system, the temperature gradient at the top surface reduces. This can be recognized in Figure 8a by reducing the slope of the curve of temperature distribution compared with those of Al2O3 systems with similar weight percentages of reinforcement content. The generation of a wider melt pool, seen in Figure 5a, or the observation of a larger temperature distribution, seen in Figure 6b, are responsible for this observation. In Figure 8b, the temperature gradient in the thickness direction (*Z*-direction) of the melt pool is larger in Al2O3-ZrO2-reinforced composites than those of Al2O3-reinforced composites. The temperature at the bottom area of the powder bed, shown on a vertical dashed line in Figure 8b, yields higher magnitudes in Al2O3-ZrO2 reinforced samples, but the slope of the plot or temperature gradient in the powder bed is higher compared with Al2O3 samples. When a large temperature gradient in the thickness direction is accompanied by a larger depth of the melt pool, as seen in

Figure 5b, this may imply a larger dissipation of laser energy through the pre-fabricated layers or metal substrate. So, it can exert a larger liquid flow in the melt pool in Al2O3-ZrO2 samples.

**Figure 7.** Optical micrographs showing the multi-line formation test results. Al2O3-reinforced samples of: (**a**) 3 wt%., (**b**) 5 wt%., and (**c**) 7 wt%. Al2O3-ZrO2 reinforced samples of: (**d**) 3 wt%., (**e**) 5 wt%., and (**f**) 7 wt%.

Liquid flow in the melt pool is an important issue in processes related to molten metals. The free surface energy, which is changed by local heating or cooling, is used to drive liquid metal. However, the term "free surface energy" is not commonly used in regard to liquids and refers to the "gradient in the surface tension" [39]. Thus, the thermal creep in different directions can be obtained by applying a temperature gradient to the surfaces, because heating a surface causes the surface tension to decrease or increase. Inhomogeneities in the gradient of the temperature of a liquid surface generate forces related to the Marangoni effect. This effect can typically be defined as a dimensionless number (*M*) in the characterization of flows, as indicated by Equation (13) [39].

$$M = \left| \frac{d\gamma}{dT} \right| \times \left| \frac{dT}{di} \right| \times \frac{L^2}{\mu \,\alpha'} \tag{13}$$

Here, *d*γ/*dT* is the surface-tension gradient or surface-tension temperature coefficient (N m−<sup>1</sup> K<sup>−</sup>1), *dT*/*di* is the temperature gradient (K m<sup>−</sup>1) along a specific direction of *i*, *L* is the characteristic length (m), μ is the dynamic viscosity (N m−<sup>2</sup> s), and α is the thermal diffusivity (m<sup>2</sup> s<sup>−</sup>1). The characteristic length (*L*) can be considered in processes involving melting, such as welding and AM, for generating deep penetration of surfactants, joints, or additives. For evaluating the convection of the minimum length between the highest temperature (*Tmax*) and the lowest temperature (*Tmin* = 300) along the melt-pool depth, as schematically shown in Figure 8c, the horizontal dashed line drawn at 300 K in Figure 8b can illustrate this length. Following the points of minimum temperatures in the *Z*-direction on the horizontal dashed line in this figure, for all the models the distance between the points of the maximum and minimum temperatures (*L*) decreased with the increase in the reinforcement content within the metal matrix. The decrease in the maximum temperature at the top surface observed in Figure 6c may be related to the decrease in the distance between the points of the maximum and minimum temperatures, i.e., *L*. However, a significant discrepancy is observed between, e.g., 3 wt% Al2O3-reinforced and 3 wt% eutectic-reinforced samples and so on. From a numerical viewpoint, this may be due to the higher melting temperature and hence the smaller melt-pool depth in the Al2O3-reinforced sample, as previously discussed. Because the thermocapillary (Marangoni) convection always flows from a lower-surface tension region to a higher-surface tension region [39], the larger

temperature gradient in each direction causes a larger Marangoni convection in that direction. As a result, when applying the Al2O3 reinforcement particle, the Marangoni flow at the top surface increases, which may benefit processes such as surface alloying/hardening. While using the Al2O3-ZrO2 reinforcement particle, the Marangoni flow at the thickness direction increases, which can benefit processes such as 3D-printing. The micro-hardness measurements of the multi-layered samples were investigated at different locations, as shown in Figure 8d, based on the American Society for Testing and Materials (ASTM) E384-16 standard. Figure 8e illustrates that with the use of reinforcement particles, an improvement in the micro-hardness of the fabricated parts is achieved. However, it can be seen that the discrepancy between the micro-hardness of Al2O3-reinforced samples and the eutectic ratio of Al2O3-ZrO2 reinforced samples becomes narrower at the top layer, indicating an improvement in reinforcement particles' distribution towards the surface in Al2O3-reinforced samples.

**Figure 8.** Temperature distribution. (**a**) Along the *Y*- direction at the top surface. (**b**) Along the *Z*-direction, i.e., thickness direction. (**c**) Definition of the distance between the points of the maximum and minimum temperatures (*L*). (**d**) The locations for measuring the micro-hardness. (**e**) Measurements of micro-hardness at different conditions and locations.

### **5. Conclusions**

Melt-pool behaviors during selective laser melting (SLM) of Al2O3-reinfored and a eutectic mixture of Al2O3-ZrO2-reinforced AISI 304 stainless-steel composites were analyzed both numerically and experimentally, and the following conclusions are drawn.


**Author Contributions:** Conceptualization, Y.J.K.; Investigation, D.A. and Y.Y.W.; Methodology, S.M.H.S.; Supervision, Y.H.M.; Validation, N.K.

**Funding:** This work was supported by a National Research Foundation of Korea (NRF) grant (No. 2012R1A5A1048294) funded by the Korean Government (Ministry of Science, ICT and future Planning).

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

### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Atomistic Simulations of Pure Tin Based on a New Modified Embedded-Atom Method Interatomic Potential**

### **Won-Seok Ko 1,\*, Dong-Hyun Kim 2,\*, Yong-Jai Kwon <sup>1</sup> and Min Hyung Lee <sup>3</sup>**


Received: 16 October 2018; Accepted: 31 October 2018; Published: 3 November 2018

**Abstract:** A new interatomic potential for the pure tin (Sn) system is developed on the basis of the second-nearest-neighbor modified embedded-atom-method formalism. The potential parameters were optimized based on the force-matching method utilizing the density functional theory (DFT) database of energies and forces of atomic configurations under various conditions. The developed potential significantly improves the reproducibility of many fundamental physical properties compared to previously reported modified embedded-atom method (MEAM) potentials, especially properties of the *β* phase that is stable at the ambient condition. Subsequent free energy calculations based on the quasiharmonic approximation and molecular-dynamics simulations verify that the developed potential can be successfully applied to study the allotropic phase transformation between *α* and *β* phases and diffusion phenomena of pure tin.

**Keywords:** tin alloy; modified embedded-atom method; molecular dynamics simulation; phase transformation; diffusion

### **1. Introduction**

The materials systems with tin (Sn) have received great attention owing to their importance in many modern technologies. In the electronics industry, soldering is the most important technique of connecting the substrate and electronic devices. Among various kinds of solders, Sn-Pb alloys are the most popular and traditional solders. Recently, Sn-based lead-free solders such as Sn-Ag-Cu alloys have replaced the traditional solders [1]. In the electrochemical industry, Sn-Li alloys are notable materials due to their applications to electrodes [2]. Moreover, liquid Sn-Li alloys are currently being considered as candidates for plasmafacing materials in tokamak fusion reactors [3,4].

For the understanding of detailed phenomena in Sn and its alloys, atomic scale simulations are effective to complement experiments providing detailed insights into the atomic scale processes. Among those simulations, density functional theory (DFT) calculation is a widely used method with many advantages such as a straightforward selectivity of alloy components and the accuracy. However, the DFT calculation is computationally demanding and cannot be applied in all applications. If the material phenomena are related to a comparably large system which cannot be covered by the DFT calculation, large-scale atomistic simulations such as molecular dynamics (MD) and Monte Carlo (MC) are highly desirable. Because the predictive accuracy of such atomistic simulations is critically dependent on the accuracy of interatomic potentials, the availability of reliable interatomic potentials

of target systems is of crucial importance. In particular, a relevant description of pure Sn should take precedence because it can be a basis for the development of interatomic potentials for concerned Sn-based multi-component alloy systems.

Until now, the development of an accurate interatomic potential for pure Sn has been a difficult task because of the complexity of the system due to the presence of allotropic phase transformation. In the periodic table, Sn is located at the borderline between covalent and metallic bonding elements. Above 286 K, Sn crystallizes into a body-centered tetragonal crystal structure (*β*-Sn) having a characteristic of metallic bonding [5]. Below 286 K, Sn crystallizes into a diamond cubic structure (*α*-Sn) having a characteristic of covalent bonding [5]. So far, several interatomic potentials for pure Sn have been developed in previous works. There are available potentials based on pair potential models [6–8], but the practical utilization of these potentials is greatly limited because these models ignore many-body bonding characteristics of pure Sn. There is an available many-body interatomic potential [9] based on embedded-atom method (EAM) model [10], but this potential cannot predict the *β*-Sn phase as a stable phase at the ambient condition because target properties of this potential are confined to phases at a high-pressure. Alternatively, there are available many-body interatomic potentials based on the modified embedded-atom method (MEAM) model [11] by Ravelo and Baskes [5], by Vella et al. [12], and by Etesami et al. [13]. It seems reasonable to apply the MEAM model to the Sn system since this model was originally devised to well describe various types of atomic bonds, including the metallic and covalent bonds, in a single formalism. However, these MEAM potentials were developed mostly focusing on properties of liquid phase and showed deficiencies in reproducing physical properties of solid phases, especially for properties of the *β*-Sn phase.

In the present study, we provide a new interatomic potential based on the second nearest-neighbor modified embedded-atom method (2NN MEAM) model [14–16]. The 2NN MEAM is a model that overcomes several shortcomings of the original MEAM model by partially considering 2NN interactions of atoms. In addition, the present study considered a relevant method of the parameter optimization to improve the reliability of the developed potential. In previous works for MEAM potentials [5,12,13], potential parameters were optimized focusing on specific physical properties obtained by experimental and density-functional theory information. In contrast, the present optimization was performed based on the force-matching method proposed by Ercolessi and Adams [17]. This method considers not a physical property itself, but forces and energies related to various atomic configurations including configurations at finite temperatures expected from the DFT calculation. It has been confirmed that this method can greatly improve the accuracy and transferability of developed interatomic potentials [18–20].

The remainder of the present article is organized as follows. Section 2 describes detailed processes of the DFT calculations and optimization of potential parameters. In Section 3, the accuracy and transferability of the developed potential is presented with suitable examples. The conclusion is finally drawn in Section 4.

### **2. Methods**

The optimization of the present interatomic potential was performed based on the force-matching method which considers both the structural energy and the forces acting on each atom as a fitting target. The overall procedure follows previous works [20,21] in which 2NN MEAM potentials for pure Ni, Ti and Li were similarly determined based on the same method in a systematic manner. First, a database of atomic energies and forces of various atomic configurations was established based on the DFT calculation. Then, an optimum set of potential parameters was determined by minimizing errors between the expectation by potential parameters and the DFT database.

### *2.1. Construction of a DFT Database*

A series of DFT calculations were performed using Vienna Ab initio Simulation Package (VASP) code [22–24] based on the projector augmented wave (PAW) method [25]. For the exchange-correlation functional, the Perdew-Burke-Ernzerhof generalized-gradient approximation (GGA) [26] was used. A plane-wave kinetic-energy cutoff of 400 eV and the Methfessel-Paxton smearing method with a width of 0.1 eV were used. All calculations were performed with a density of the *k*-point mesh equivalent to a 21 × 21 × 21 mesh for the face-centered cubic (fcc) primitive unit cell, and the corresponding similar density of the *k*-point mesh were employed for other unit cells and supercells. In the present DFT calculation, atomic configurations resulting from various conditions were considered for the construction of the DFT database to ensure the sufficient transferability of the developed interatomic potential to various possible applications. Detailed information on cell structures and corresponding *k*-point meshes is summarized in Table 1.

**Table 1.** Atomic configurations considered in the present density functional theory (DFT) calculations. In the Stability column, "stable" indicates that the structure is reported in an equilibrium phase diagram. Other phases are labeled as "hypothetical". The Strain column indicates the strain applied to the supercells, where *H*, *O*, and *M* denote hydrostatic, orthorhombic, and monoclinic strains, respectively.


The equilibrium lattice constants and bulk modulus at 0 K were calculated by employing the Birch-Murnaghan equation of state [27,28]. For the calculation of properties involved with defects, positions of each atom were relaxed at a constant volume and cell shape. The vacancy migration energy was calculated with a suitable saddle-point configuration utilizing the nudged elastic band (NEB) method [29,30]. For the calculation of the surface energy, rectangular cells with a stacking of 12 Å to 13 Å thick slab and a vacuum region of 10 Å were employed without the relaxation of cell dimensions parallel to free surfaces. The convergence criteria for energies and forces of all defect calculations were 10−<sup>6</sup> eV/atom and 10−<sup>2</sup> eV/Å, respectively. Phonon calculations were performed using the "Phonopy" code [31,32] based on the direct force constant approach [33]. Supercells of 64 atoms were used for the phonon calculation for *α* and *β* phases with the convergence criteria for energy and forces of 10−<sup>8</sup> eV and 10−<sup>4</sup> eV/Å, respectively. To calculate the cohesive energy, a reference energy with a single atom was obtained by considering the spin-polarized calculation.

To obtain a sufficient number of effective force data at finite temperatures, two-step DFT calculations were conducted. First, ab initio MD simulations [22] were performed for a total of 2000 steps with a timestep of 1.5 fs using relatively low convergence parameters such as a single *k*-point and a default value of the cutoff energy. To determine accurate energies and forces of each configuration, well-converged calculations with a higher cutoff energy (400 eV) and denser *k*-point mesh were then followed using randomly extracted configurations from ab initio MD simulations.

### *2.2. Optimization of Potential Parameters*

For the development of a unary potential based on the 2NN MEAM formalism, an optimization of 14 independent potential parameters is required. Four parameters [the cohesive energy (*Ec*), the equilibrium nearest-neighbor distance (*re*), the bulk modulus (*B*) of the reference structure and the adjustable parameter *d*] are involved with to the universal equation of state. Seven parameters [the decay lengths (*β*(0), *β*(1), *β*(2), and *β*(3)) and the weighting factors (*t* (1), *t* (2), and *t* (3))] are involved with the electron density. The parameter *A* is required for the embedding function and two parameters [*C*min and *C*max] are required for the many-body screening. Detailed explanations on these potential parameters can be found in literature [14–16].

The target configurations in the DFT database consist of configurations of stable phases (*α*, *β*, and liquid) and hypothetical phases (fcc, bcc, and hcp) at various temperatures. In particular, the inclusion of configurations at finite temperatures is indispensable for the force-matching because these configurations can provide effective atomic force data for the optimization. A total of 63 configurations resulting from various temperatures and strain conditions as well as various defect configurations were considered for the fitting. Detailed information on atomic configurations used for the fitting process is listed in Table 1.

The optimization was performed by comparing energies of target configurations and forces on each atom expected by a candidate set of potential parameters and those expected by the DFT calculation. The optimization started with specifying a reference structure of the potential parameters, a radial cutoff distance, and fitting weights of each configuration. An optimizer based on the genetic algorithm then adjusted the candidate set of parameters to reduce the sum of energy errors of each configuration and force errors of each atom. If the derived set of potential parameters does not satisfactorily reproduce overall physical properties, another optimization was performed with a different set of the reference structure, radial cutoff value, and fitting weights of configurations. The trial reference structures were diamond cubic, fcc and bcc structures, and the trial radial cutoff values were 4.5, 5.0, 5.5 and 6.0 Å. During the optimization process, we realized that the consideration of similar weighting between configurations of *α* and *β* phases results in a general worsening of the reproducibility of various physical properties. Considering the importance of the *β* phase stable at ambient condition, the final optimization was performed using increased fitting weights for atomic configurations of the *β* phase.

Figure 1 shows results of the final optimization represented by statistical correlations between the calculated energies and forces from the developed potential and the DFT calculation. Each correlation was obtained using atomic configurations expected by the DFT calculation without further atomic relaxations. The resultant root mean square (RMS) errors of energies and forces are 0.044 eV/atom and 0.097 eV/Å, respectively. These values are significantly higher than the results of previous optimization for pure Li [21]. It seems that this is because the handling of phases with different types of atomic bonding (covalent and metallic) simultaneously is significantly more difficult than the handling of phases with single type of atomic bonding (metallic) based on a single potential formalism.

**Figure 1.** Scatter plots for (**a**) energies and (**b**) forces of pure Sn with respect to the DFT database. The dashed lines indicate a perfect correlation between the present 2NN modified embedded-atom method (MEAM) and the DFT values.

We finally confirmed that a reference structure of fcc provides the optimum result in reproducing various physical properties of pure Sn. A cutoff value of 5.0 Å was confirmed to be large enough to reproduce various physical properties with an acceptable computational efficiency. In the subsequent section, all calculations based on the new MEAM potential were performed using this radial cutoff distance. Table 2 lists the final set of potential parameters for pure Sn.

**Table 2.** Optimized 2NN MEAM potential parameter set for the pure Sn system. The following properties are dimensionful: the cohesive energy *Ec* (eV/atom), the equilibrium nearest-neighbor distance *re* (Å), and the bulk modulus *B* (1012 dyne/cm2). The reference structure is fcc Sn.


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

In this section, results of atomistic simulations performed by the developed potential are presented. The results are compared to corresponding experiments and DFT calculations to evaluate the accuracy and transferability of the developed potential. In addition, the performance of the developed potential was compared with previous MEAM potentials by Ravelo and Baskes [5], by Vella et al. [12], and by Etesami et al. [13]. We used the potential of Ravelo and Baskes [5] in a form adopted by Vella et al. [12]. For the clear comparisons between the present potential and the previous potentials, all physical properties were recalculated in the same manner, whether or not some of properties were already reported in previous studies [5,12,13]. All atomistic simulations were performed based on the LAMMPS code [34]. If not specially designated as MD simulations, obtained properties represent results of molecular statics simulations at 0 K using atomic configurations of at least 4000 atoms. All MD runs were performed starting from initial configurations optimized by corresponding molecular statics simulations at 0 K using a timestep of 0.002 ps. Nosé-Hoover thermostat and barostat [35,36] were used for controlling temperature and pressure, respectively.

In Table 3, calculated bulk, elastic and defect properties of pure Sn at 0 K using the present potential are compared with experimental data, DFT calculations and results using previous MEAM potentials. The experimental cohesive energy, which is also in good agreement with the DFT expectation, is well reproduced by the present potential. The ground state structure of pure Sn expected by the DFT calculation is the *α* structure. The present potential well reproduces such tendency and the structural energy differences between various stable and hypothetical solid phases as well. Among structural energy differences, the difference between the low temperature stable *α* phase and the high temperature stable *β* phase is important for the expectation of the allotropic transformation. The present potential accurately reproduces this trend while previous potentials by Vella et al. [12] and by Etesami et al. [13] show significantly higher energy differences. The structural parameters such as lattice constants of the *α* phase and lattice constants and *c*/*a* ratio of the *β* phase are also well reproduced by the present potential in closer agreement with experimental data.

In the elastic properties, all MEAM potentials generally present a difficulty in reproducing the *C*<sup>44</sup> of the *β* phase. Especially, previous potentials by Ravelo and Baskes [5] and by Etesami et al. [13] indicate the deviation more than one order of magnitude. The present potential significantly improves the reproducibility of the *C*<sup>44</sup> and other elastic details as well. For example, it was reported that the *β* phase lattice is stiffer in the *c*-direction compared to *a*- or *b*-direction (*C*<sup>33</sup> > *C*11). The present potential correctly reproduces this trend as well as the absolute values of each elastic constant.

The defect properties of the *β* phase such as the properties related to the vacancy and the free surface are also examined and listed in Table 3. In the present study, the activation energy of vacancy diffusion (*Q*vac) is defined as a sum of the vacancy formation energy (*E*vac *<sup>f</sup>* ) and the vacancy migration energy (*E*vac *<sup>m</sup>* ) considering the reported substitutional diffusion mechanism of the *β* phase [37]. Based on previous experimental study [38] and the present DFT calculation, the self-diffusivity of the *β*

phase is anisotropic along crystal directions, and the activation energy of vacancy diffusion along the *c*-direction ( *c*) is higher than that along the *a*- or *b*-direction (⊥*c*). The present potential successfully reproduces this trend. In the calculation of the vacancy migration energy and resultant activation energy of vacancy diffusion, we realized that the calculation along the *a*- or *b*-direction by previous potentials indicates unphysical results such as the negative migration energy and the instability of the saddle point position. These values are thus not presented in Table 3. The surface energies of the *β* phase are also correctly reproduced only by the present potential while previous potentials show a significant overestimation.

**Table 3.** Calculated bulk, elastic and defect properties of pure Sn using the present 2NN MEAM potential, in comparison with experimental data, DFT data, and the calculation results using previous MEAM potentials by Ravelo and Baskes [5], Vella et al. [12], and Etesami et al. [13]. The following quantities are listed: the cohesive energy *Ec* (eV/atom), the lattice constant *a* (Å), the bulk modulus *B* and the elastic constants *C*11, *C*12, *C*13, *C*33, *C*<sup>44</sup> and *C*<sup>66</sup> (1012 dyne/cm2), the structural energy differences Δ*E* (eV/atom), the vacancy formation energy *E*vac *<sup>f</sup>* (eV), the vacancy migration energy *<sup>E</sup>*vac *m* (eV), the activation energy of vacancy diffusion *Q*vac (eV), and the surface energies *E*surf (erg/cm2) for the orientations indicated by the superscript. For the *E*vac *<sup>m</sup>* and *Q*vac, values considering diffusional paths along the *a*- or *b*-direction (⊥*c*) and along the *c*-direction ( *c*) of the *β* phase are presented. The DFT and MEAM calculations were performed at 0 K while the experimental data were obtained at finite temperatures.


<sup>a</sup> Ref. [39]. <sup>b</sup> Ref. [40]. <sup>c</sup> Ref. [41]. <sup>d</sup> Ref. [42]. <sup>e</sup> Ref. [38]. <sup>f</sup> Present DFT calculation. <sup>g</sup> Present calculation using the reported potential parameters. <sup>h</sup> The potential parameters adopted by Vella et al. [12] are used.

We further examined the transferability of the developed potential by comparing a group of properties at finite temperatures. Table 4 lists the thermal properties of the *β* phase (the thermal expansion coefficient and the specific heat), in comparison with experimental data. These properties were obtained based on the MD simulations using an isobaric-isothermal (*NPT*) ensemble at the target temperature and zero pressure. For the thermal expansion coefficient of the *β* phase, the potential by Ravelo and Baskes [5] shows a better agreement with the experiment than other potentials. Instead, other potentials show a better reproducibility than the potential by Ravelo and Baskes [5] in the case of the specific heat of the *β* phase.

**Table 4.** Calculated thermal properties of pure Sn using the present 2NN MEAM potential, in comparison with experimental data and the calculation results using previous MEAM potentials by Ravelo and Baskes [5], Vella et al. [12], and Etesami et al. [13]. The listed quantities correspond to the thermal expansion coefficient *ε* (10−6/K), the heat capacity at constant pressure *CP* (J/mol K), the melting temperature *Tm* (K), the enthalpy of melting Δ*Hm* (kJ/mol), and the volume change upon melting Δ*Vm*/*V*solid (%).


<sup>a</sup> Ref. [41]. <sup>b</sup> Ref. [43]. <sup>c</sup> Ref. [44]. <sup>d</sup> Present calculation using the reported potential parameters. <sup>e</sup> The potential parameters adopted by Vella et al. [12] are used.

Figure 2 shows the temperature dependence of the atomic volume of pure Sn. Initially, the *β* phase with 16,224 atoms was equilibrated at 5 K, and the temperature was then gradually increased 1000 K with a heating rate of 1.0 K/ps using the *NPT* ensemble at zero pressure. Each potential indicates a discontinuous jump in the volume at a certain temperature. This jump indicates the occurrence of the melting (*β* → liquid). Because the heating simulation was performed without any heterogeneous nucleation sites, the melting occurs at a temperature much higher than the equilibrium melting temperature. Therefore, this temperature can be regarded as an overheated melting point. The equilibrium melting temperature at zero pressure was further calculated using the interface method [45,46] employing a simulation cell consisting of liquid and solid phases and listed in Table 4. The table also lists the enthalpy change and the volume change upon melting which were obtained at the calculated equilibrium melting temperature. For the previous potential by Ravelo and Baskes [5], it was reported that the liquid part of the interface simulation is crystallized into a structure different from the *β* phase [12], and thus the equilibrium melting point and resultant enthalpy and volume changes are not presented. The developed potential shows discrepancies in the properties associated with the melting while the potential by Etesami et al. [13] indicates the best agreement in the melting point with the experiment. This can be interpreted by the fact that the potential [13] was developed focusing mostly on the melting phenomenon without consideration to the reproducibility of many other important properties.

As listed in Table 4, the present potential is somewhat deficient in describing properties related to the melting compared to previous potentials. The present potential underestimates the enthalpy of melting and overestimates the volume change upon melting compared to experimental values. Considering that the present potential satisfactory describes physical properties of *β* phase, it is necessary to further investigate the properties associated with the structure of the liquid phase. Figure 3 shows the calculated radial distribution function of liquid structures at different temperatures. Although the present potential shows a generally acceptable quality to reproduce the height and position of the first peak, it shows a deficiency in reproducing characteristics of other peaks. We attribute these deficiencies to the use of increased fitting weights of the *β* phase compared to other phases during the fitting process.

**Figure 2.** Atomic volume dependence of pure Sn, calculated using the present 2NN MEAM potential, in comparison with the calculation results using previous MEAM potentials by Ravelo and Baskes [5], Vella et al. [12], and Etesami et al. [13]. The heating simulations were started with the *β* phase at 5 K. The discrete jumps represent the occurrence of the melting.

**Figure 3.** Calculated radial distribution function of liquid Sn at (**a**) 573 K and (**b**) 1373 K using the present 2NN MEAM potential, in comparison with experimental data (Ref. [7]) and the calculation results using previous MEAM potentials by Ravelo and Baskes [5], Vella et al. [12], and Etesami et al. [13].

We further examined the reproducibility of phonon properties and resultant properties at finite temperature based on the harmonic approximation (HA) and quasiharmonic (QHA) approximation. Figure 4 shows the phonon density of states (DOS) of *α* and *β* phases expected by the present and previous potentials compared to the DFT expectation. As shown in Figure 4a, the phonon DOS of the *α* phase is closely reproduced by the present potential. In contrast, previous potentials indicate a significant overestimation of the phonon DOS at high frequencies compared to the DFT expectation. A similar trend is also presented for the phonon DOS of the *β* phase as shown in Figure 4b. Moreover, previous potentials cannot reproduce the stability of the *β* phase under perturbative forces. As shown in the magnified figure of the phonon DOS of the *β* phase (Figure 4c), previous potentials indicate dynamical instability (imaginary phonon frequencies) in contrast to expectations by the present potential and DFT calculations. Even though this problem of previous potentials is expected to interfere with various possible applications of previous potentials, the present potential is free from this problem.

**Figure 4.** Phonon density of states (DOS) of (**a**) *α* phase and (**b**) *β* phase, calculated based on the harmonic approximation (HA). The results using the present 2NN MEAM potential are compared with the present DFT results and the results using previous MEAM potentials by Ravelo and Baskes [5], Vella et al. [12], and Etesami et al. [13]. (**c**) A magnified figure of Figure (**b**) is presented to clarify the existence of imaginary phonon modes.

Therefore, only the present potential can be further utilized to expect the physical properties at finite temperatures based on the QHA. One simple example is calculations of the thermal expansion coefficient and the heat capacity as presented in Figure 5. The present potential exhibits an acceptable quality comparable to the DFT calculation when reproducing such properties. Another important example is the expectation of the allotropic phase transformation between *α* and *β* phases. Because the present potential can accurately reproduce the phonon properties without suffering from the imaginary phonon modes, a vibrational contribution on the free energy can be calculated based on the QHA. The Gibbs free energy at the given temperature *T* and pressure *P* can be obtained from the Helmholtz free energy *F*(*T*, *V*) at the given temperature and volume *V* through the transformation,

$$G(T, P) = \min[F(T, V) + PV],\tag{1}$$

where the right-hand side of this equation represents a value when the minimum value is found in square brackets by varying the volume. In the present study, the *F*(*T*, *V*) was computed by a sum of the DFT total energy *EDFT*(*V*) at the given volume, and the vibrational contribution to the free energy *Fvib*(*T*, *V*) at the given temperature and volume.

$$F(T, V) = E\_{DFT}(V) + F\_{vib}(T, V) \tag{2}$$

A detailed explanation on the determination of the free energy based on the phonon calculation is given in Ref. [47].

**Figure 5.** (**a**) Thermal expansion coefficient (*ε*) and (**b**) heat capacity (*CP* ) of *β*-Sn, calculated based on the quasiharmonic approximation (QHA). The results using the present 2NN MEAM potential are compared with those using the present DFT calculation and experimental data (Refs. [43,48]).

Figure 6 shows the calculated temperature dependence of Gibbs free energies of *α* and *β* phases when the reference state is set to the *β* phase at 0 K. As stated, both DFT and the present MEAM potential expect that the *α* phase is a ground state of pure Sn. The Gibbs free energy of each phase decreases with increasing temperature due to the increasing contribution of the vibrational entropy. Because the vibrational contribution is more significant in the *β* phase than in the *α* phase, the *β* phase starts to become more stable than the *α* phase at a certain temperature. We define this temperature as the phase transformation temperature between *α* and *β* phases. The calculated transformation temperature by the present MEAM potential (460 ± 5 K) agrees well with that by the present DFT calculation (470 ± 5 K) while these values are higher than experimental data (286 K) [5]. Although this discrepancy seems to be caused by the limitation of the QHA, the present potential at least reproduces the DFT expectation.

**Figure 6.** Relative Gibbs free energies of *α* and *β* phases at various temperatures, calculated based on the quasiharmonic approximation (QHA). The results using (**a**) the present 2NN MEAM and (**b**) the present DFT calculation are illustrated. Each energy value was calculated with respect to the reference state of the *β* phase at 0 K.

As a final confirmation of the transferability of the developed potential, the self-diffusion of liquid Sn was examined using the MD simulation. In the atomic scale, the MD simulation provides a time dependency of the mean square displacement (MSD), and this is related to the bulk diffusivity based on the following equation (Einstein relation),

$$D = \frac{\left< R^2(t) \right>}{6t},\tag{3}$$

where *D* is the diffusivity, - *R*2(*t*) is the MSD of atoms, and *t* is the time.

The self-diffusivity of liquid Sn was calculated using *NPT* ensemble MD runs starting with the *β* phase of 16,224 atoms. First, the initial structure was melted by maintaining the cell at a high temperature (2000 K) for a sufficient time (40 ps). The temperature was then changed to each target temperature within a range of 600 and 2000 K, and the relaxation was performed for 40 ps. The simulation cell was further maintained at the target temperature for a total simulation time of 2000 ps, and the time evolution of the MSD was counted at every 4 ps. We confirmed that the MSD of atoms accumulated through these conditions shows a clear linear relationship with the time. Figure 7 shows resultant self-diffusivities of the liquid phase at various temperatures (600–2000 K) compared with experimental data. Despite a deviation between experimental results, the present potential successfully reproduces the general trend of an experimental result by Bruson et al. [49]. The calculated self-diffusivities using the present potential are similar to those using the previous potential by Vella et al. [12] even though the values are slightly underestimated at low temperatures. Compared to the previous potentials by Ravelo and Baskes [5] and by Etesami et al. [13], the present potential indicates significantly improved results, especially at high temperatures. It is interesting to note that the self-diffusivities of the liquid phase are sufficiently reproduced despite the discrepancy of properties related to the liquid structure. The result seems to emphasize the effectiveness of the force-matching method to obtain a robust potential that demonstrates a good transferability to kinetic properties such as the migration energy and the attempt frequency of diffusing atoms. However, it should be noted that the transferability of the present potential to kinetic properties of low-coordinated atoms on surfaces (surface diffusion) is not guaranteed because the present fitting was mostly performed by using atomic configurations in the bulk state.

**Figure 7.** Calculated self-diffusivity of liquid Sn using the present 2NN MEAM potential, in comparison with experimental data (Refs. [49–52]) and the calculation results using previous MEAM potentials by Ravelo and Baskes [5], Vella et al. [12], and Etesami et al. [13].

To summarize, we have shown that the present potential generally reproduces various physical properties of pure Sn, even though there is a difference in the performance for each property. The developed potential performs very well in describing structural, elastic and thermal properties of the *β* phase, phonon properties of solid phases, and diffusion properties of solid and liquid phases. It is

less suited to describing the melting behavior and the liquid structure compared to previous MEAM potentials [5,12,13], which were developed to target specific physical properties of the liquid phase. This fact should be kept in mind in future applications of the present potential. The LAMMPS files implementing the developed potential can be obtained from an online repository [53].

### **4. Conclusions**

We provide a new robust interatomic potential for pure Sn on the basis of the 2NN MEAM model. The potential is developed based on the force-matching method utilizing the DFT database of energies and forces of atomic configurations under various conditions. The developed potential significantly improves the reproducibility of various physical properties compared to previously reported MEAM potentials, especially for properties of the *β* phase. The present potential exhibits superior transferability to the phonon properties and the properties related to diffusion phenomena. As a possible example, the allotropic phase transformation between α and *β* phases is investigated based on the QHA. The results indicate that the present potential can be successfully used for the expectation of the transformation temperature at least at the level of the DFT expectation. The present potential also successfully reproduces the self-diffusivity of liquid Sn. The developed potential can be a suitable basis for implementing atomistic simulations of technically important Sn-based alloy systems.

**Author Contributions:** Conceptualization, W.-S.K. and D.-H.K.; Methodology, W.-S.K. and D.-H.K.; Validation, W.-S.K., D.-H.K, Y.-J. K, and M.H.L.; Formal analysis, W.-S.K. and D.-H.K.; Investigation, W.-S.K. and D.-H.K.; Resources, Y.-J.K.; Data curation, W.-S.K.; Writing—original draft preparation, W.-S.K.; Writing—review and editing, D.-H.K, Y.-J.K, and M.H.L; Visualization, W.-S.K.; Project administration, D.-H.K.; Funding acquisition, M.H.L.

**Funding:** This work is supported by KITECH research fund (Grant No. UR180012), and National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP) (Grant No. NRF-2017R1C1B5015038).

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

### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Identification of the Flow Properties of a 0.54% Carbon Steel during Continuous Cooling**

**Christoph Rößler 1,2, David Schmicker 1,2, Oleksii Sherepenko 3, Thorsten Halle 3, Markus Körner 2,3, Sven Jüttner 3,\* and Elmar Woschke <sup>1</sup>**


Received: 16 December 2019; Accepted: 4 January 2020; Published: 9 January 2020

**Abstract:** The determinination of material properties is an essential step in the simulation of manufacturing processes. For hot deformation processes, consistently assessed Carreau fluid constitutive model derived in prior works by Schmicker et al. might be used, in which the flow stress is described as a function of the current temperature and the current strain rate. The following paper aims to extend the prior mentioned model by making a distinction, whether the material is being heated or cooled, enhancing the model capabilities to predict deformations within the cooling process. The experimental identifaction of the material parameters is demonstrated for a structural carbon steel with 0.54% carbon content. An approach to derive the flow properties during cooling from the same samples used at heating is presented, which massively reduces the experimental effort in future applications.

**Keywords:** flow stress, hot deformation, carbon steel, continuous cooling, phase transformations

### **1. Introduction**

The input of manufacturing process simulations has to include information about the actual process as well as the geometries and the materials being processed. Geometry and process information tend to be easier to obtain, because those are constantly monitored for quality assurance, for instance. The assessment of the material information on the other hand is more challenging, because several factors have to be considered.

A manufacturing process involving a broad range of temperatures and strain rates is rotary friction welding (RFW). In this process, the energy to form a permanent bond is directly introduced in the joining zone in form of frictional heat, which is generated by pressing the parts together while a relative motion is performed. Therefore, through this process temperatures close to the liquidus temperature of the softer material are achieved [1,2] and the strain rates reach one-digit values [3,4]. To simulate RFW processes, Schmicker et al. [4–7] elaborated a non-linear fluid model to describe the material flow. Although there are material models more capable in representing the dependencies of the yield stress on the deformation temperature, the strain rate, the degree of deformation, the phase composition and more, this model is suitable for RFW and for other hot deformation processes, too. A major advantage are the fairly expedient material parameters that can be gathered with only a few samples in a tensile testing routine, in which the test specimen is continuously heated. Besides conducting own experiments, the model can also be parametrized using data from encyclopedias such as [8,9], which significantly lowers the necessary efforts for the application of the process simulation.

As Schmicker et al. define the flow stress solely by steady-state values, the material point history is not taken into account. Similarly, aforementioned encyclopedias may also contain only the first three factors listed previously.

During heating, the microstructure of carbon steels starts to change to austenite by diffusion processes, if the temperature is greater than *Ac*<sup>1</sup> and which finish at *Ac*3. Both, *Ac*<sup>1</sup> and *Ac*3, depend on the rate of heating, which is documented in time-temperature-austenitization (TTA) diagrams as found in [10]. During cooling, the reverse transformation deviates from this and primarly depends on the cooling rate, which is documented in continuous cooling transformation (CCT) diagrams as depicted in Figure 1. Further information in the CCT include the microstructural composition as well as the hardness *H*, typically in Vickers hardness (HV),which occur at different cooling paths.

The microstructural changes have an impact on the flow properties of the material and the properties during cooling therefore differ from the properties during heating [11]. This should be taken into account in the process simulation if the material is still being deformed during cooling and independently of this in the residual stresses analysis as the yield strength limits the stress formation.

**Figure 1.** CCT diagram of a 0.55% carbon steel austenitized at 840 ◦C [12].

The current paper aims to quantify the flow properties of a 0.54% carbon steel during continuous cooling and to link the gathered results to the properties at heating, which are assumed to match the steady-state reference values, because the investigated temperature range is limited around and below *Ac*<sup>1</sup> temperature, at which the diffusional transformations just start.

To experimentally determine the material properties, a prior developed testing routine by Schmicker is extended by tests with continuous cooling. Since there exist an infinite number of cooling rates, a practical approach is proposed, which involves the hardnesses and transformation temperatures found in CCT diagrams as shown in Figure 1. This approach based on reference flow properties and hardness values aims to ensure that all model parameters can also be determined using the data available in literature collections. Concerning hot forming processes, data found in deformation CCT diagrams [13] indicate that compared to the cooling rate, the effect of deformations on the hardness is significantly smaller.

### **2. Materials and Methods**

### *2.1. Mathematical Approach*

The description and evaluation of the plastic material properties presented below is based on the following assumptions:

• The material behaves isotropic regarding all mechanical properties.

*Metals* **2020**, *10*, 104


In the Norton–Bailey constitutive model [14], in the double logarithmic depiction the stress–strain rate relation is linear, which can be formulated as

$$
\dot{\varepsilon} = A \sigma^n \tag{1}
$$

or as *<sup>ε</sup>*˙

$$\frac{\dot{\varepsilon}}{\dot{\varepsilon}\_0} = \left(\frac{\sigma}{\sigma\_0}\right)^n \tag{2}$$

in which *σ*0(*θ*) and *n*(*θ*) are material parameters specified for a reference strain rate *ε*˙0 and as a function of temperature *θ*. Equation (2) is valid for a wide range of strain rates and materials and documented in several data collections [8,9]. Under the premise that all isothermal flow curves *σiso*(*θiso*,*ε*˙) intersect in one characteristic point *C*(*σC*,*εC*), Schmicker et al. [7] bypass the determination of *n*(*θ*) for every temperature. Therefore, knowing the point *C*, the Norton–Bailey exponent can be expressed as

$$m = \frac{\log\left(\frac{\ell\_C}{\ell\_0}\right)}{\log\left(\frac{\sigma\_C}{\sigma\_0}\right)}\tag{3}$$

as a feature in the so-called consistently assessed Carreau fluid model. The strain rate sensitivity typically increases with increasing temperatures and to produce higher strain rates, higher stresses are required, which is ensured by

$$
\sigma\_{\mathbb{C}} > \max(\sigma\_0(\theta, \varepsilon\_0))
$$

$$
\varepsilon\_{\mathbb{C}} > \varepsilon\_0
$$

To distinguish heating and cooling, the flow properties during continuous cooling are subsequently denoted by an apostrophe. Two effects are taken into account for.

Firstly, a heat treatment effect, in which the material either hardens or softens due to microstructural changes. Rapid quenching causes the formation of martensite, which achieves more than twice the hardnesses than ferrite and pearlite. In annealing processes on the other hand, the cooling is typically slow to avoid this transformation.

Secondly, a transformation inertness that causes the transformation to shift to other temperature ranges depending on the cooling rate. The quicker the cooling process, the less is the time for the carbon diffusion processes. If the diffusion can not take place at all, the carbon become trapped in a body-centered tetragonal lattice configuration below martensite start temperature. In Figure 1 it is also seen that even for very slow cooling, the transformation starts well below *Ac*3.

Assuming that *σ* <sup>0</sup> during cooling is not necessarily identical to *σ*<sup>0</sup> at heating, but similarly shaped, the two curves are linked by adding offset parameters *κ* and *θκ*

$$
\sigma\_0'(\theta, \dot{\theta}) = \kappa \sigma\_0(\theta + \theta\_\mathbb{K}) \tag{5}
$$

for prior discussed transformation effects, depending on *θ*, ˙ *θ*.

The hardening factor *κ* can be interpreted as a vertical scaling of *σ*0(*θ*) to account for the heat treating effect as prior applied by Rößler et al. [15]. For its evaluation the linear relation in-between the yield strength *σ<sup>y</sup>* and the Vickers hardness *H* [16,17]

$$
\sigma\_y = aH + b \tag{6}
$$

is utilized, in which it is physically reasonable that *b* is zero. It should be mentioned that for other hardness scales the correlation can be non-linear. At room temperature, *κ* can be expressed as the proportion

$$\kappa(\theta\_0, \dot{\theta}) = \frac{\sigma\_0'(\theta\_0, \dot{\theta})}{\sigma\_0(\theta\_0 + \theta\_\kappa)} = \frac{H\_0'(\dot{\theta})}{H\_0} \tag{7}$$

in which *H*<sup>0</sup> is the initial hardness corresponding to *σ*0(*θ*0) and *H* <sup>0</sup> the hardness after cooling with a specific rate. To couple the hardening to the actual transformation, sigmoid function

$$\kappa(\theta,\theta) = \kappa(\theta\_0,\theta) - \frac{\kappa(\theta\_0,\dot{\theta}) - 1}{2} \tanh\left(\frac{2}{\theta\_r(\dot{\theta}) - Ac\_1}(\theta - Ac\_1)\right) \tag{8}$$

limits hardening to the lower temperature range. The parameters *H* <sup>0</sup> and *θ<sup>r</sup>* are either found in the CCT diagram (Figure 2) or can be experimentally determined by indentation and dilatometric testing. The use of the *Ac*<sup>1</sup> temperature is a recommendation for a free value of this equation, because it guarantees (*θ*( ˙ *θ*) − *Ac*1) < 0, which must always be fulfilled for mathematical reasons.

**Figure 2.** Extraction of *κ*, *θκ* for fast and slow cooling rates.

To shift *σ*0(*θ*) horizontally due to a transformation inertness, the delay temperature

$$\theta\_{\mathbf{k}}(\theta,\dot{\theta}) = \frac{Ac\_1 - \theta\_r(\dot{\theta})}{2} \tanh\left(\frac{2}{\theta\_r(\dot{\theta}) - Ac\_1}(\theta - Ac\_1)\right) \tag{9}$$

is similar defined as *κ*. It is worth mentioning that *κ* and *θκ* actually start to raise above *θ<sup>r</sup>* to compensate discontinuities of *σ* <sup>0</sup> introduced by large *θκ* for martensitic transformations, for instance.

Concerning Equation (3), the characteristic intersection point *C* has to be re-evaluated, too, to satisfy (4) for *κ* 1. To maintain the same high strain rate senstivity at high temperatures around the melting point *θ<sup>M</sup>* and assuming that the low sensitivity for low temperatures will not change, point *C* is identified using

$$\begin{aligned} n'(\theta\_M) &= n(\theta\_M) \\ n'(\theta\_0) &= n(\theta\_0) \end{aligned} \tag{10}$$

To estimate hardnesses not documented in the CCT diagram, law of mixture

$$H = \sum\_{i} \xi\_{i} H\_{i} \tag{11}$$

as presented by Ion et al. [18] can be applied, in which *ξ* are the phase fractions of ferrite, pearlite, bainite, martensite, and austenite. The calculation of the phase fraction might be done numerically using evolution equation

$$\dot{\xi}(t) = \frac{\ddot{\xi}\_{\infty}(\theta) - \ddot{\xi}(t)}{\tau(\theta)} \tag{12}$$

by Leblond and Devaux [19].

### *2.2. Materials*

A detailed summary of the chemical composition of the investigated steel is provided in Table 1. The round tensile specimens used in tensile testing are machined from rods measuring 10 mm in diameter.



### *2.3. Experimental Setup*

The uniaxial tensile tests were carried out using Gleeble 3500 testing machine at Otto von Guericke University Magdeburg. Part of its capabilities was the execution of displacement controlled, uniaxial tensile tests, while the temperature of the test specimen could be varied using conductive heating at the same time. Figure 3 provides an overview of the setup. Within the machine, the force *F*(*t*) required to produce the stroke *s*(*t*) was captured and an additionally installed extensometer monitored the diameter *d*(*t*) in the axial center of the test specimen. At the same location, the temperature *θ*(*t*) was measured by welded-on K-type thermo couples, which were used to control the amperage for the conductive heating. The initial diameter in the monitored section was 6 mm based on the permissible force of the testing machine. For cooling upto 30 K/s, two air cooling jets were installed, which aimed at the center of the test specimen. The rate of 30 K/s could be maintained down to 250 ◦C, which was well below the martensite start temperature.

**Figure 3.** Gleeble 3500 machine setup (**a**) installed test specimen, (**b**) positions of the air cooling jets.

In both routines, the specimens were exposed to linear temperature and displacement profiles in the evaluation range (Figure 4). This allowed us to gather information about the flow properties for a number of temperatures in a single test. In order to examine different ranges of *θ* and *ε*˙, the start temperature *θin* was varied as well as the stroke rate *s*˙. A drawback of altering these parameters simultaneously was that factors such as crystal recovery, recrystallization, and grain growth remained inseparable from the thermal dependence and were therefore averaged in the evaluation of the tensile tests. To overcome elastic deformations, pre-stroke *sin* was set to 0.3 mm, before the actual evaluation range begins. A comparitive summary of the test conditions during heating and cooling is provided in Table 2.

**Figure 4.** Schematic thermal and axial stroke profiles in the tests during (**a**) heating and (**b**) cooling.



In the cooling tests, all test specimens were heated up to 1200 ◦C to trigger the phase transformation to austenite. The controlled cooling process was then started by turning on the air jets, where the central thermocouple was used to ensure constant cooling rates. The cooling rates ˙ *θ* were exemplarly set to 5 K/s, at which it was assumed that the microstructure was similar to the initial one, and to 30 K/s, at which martensite is formed. For heating, test data for about 40 specimens were available by Jüttner and Körner [20] and additional 20 tests were performed for both examined ˙ *θ*, because of the narrowed temperature evaluation range.

To evaluate Equation (7), Vickers hardness measurements were carried out across the center of the lateral cross-sections of test specimens, tested with the same cooling routine without applying any axial forces. The indentation locations were arranged in three lines, each of which got 10 points spaced over the diameter, with 2 mm axial spacing.

### **3. Results**

Based on the measured data for *F*(*t*) and *d*(*t*), the longitudinal true stress

$$
\sigma = \frac{4F}{\pi d^2} \tag{13}
$$

and the longitudinal true strain rate

$$
\dot{\varepsilon} = -2\frac{d}{dt}\tag{14}
$$

were derived.

In order to fit the model according to Equation (3), the generated *σ*-*ε*˙-*θ* triples were classified by *θ* and least squares method was applied to identify *C*(*σC*,*εC*) and a monotonous *σ*0(*θ*,*ε*0) in a one-step optimization. All measured data and the fitted models are depiced in Figure 5a for the heating experiments respectively Figure 5b for cooling at 5 K/s and Figure 5c at 30 K/s. Tabular summaries of the fitted models are given in Tables 3 and 4. The coefficients of determination

$$R^2 = \frac{\sum\_{i} (\log \sigma(\theta\_{i\prime} \dot{\epsilon}\_i) - \log \vartheta)^2}{\sum\_{i} (\log \sigma\_i - \log \vartheta)^2} \tag{15}$$

are listed in Table 3, confirming the suitablity of the consistently assessed Carreau fluid model for the investigated steel.

Table 5 lists the averaged hardness measurements, which show that at 30 K/s more than twice as high hardnesses arises, which is linked to the formation of martensite. Cooling by 5 K/s produces hardnesses about 10 HV below the initial hardness.


**Table 3.** Correlation of the experimental data and fitted model.


**Table 5.** Vickers hardness measurements (averaged) and transformation temperatures taken from Figure 1.


**Figure 5.** Flow properties examined (**a**) during heating, (**b**) at 5 K/s cooling, (**c**) at 30 K/s cooling; The surface plots represent the fitted consistently assessed Carreau fluid models for each routine.

### **4. Discussion**

Using the presented testing routine during cooling, it becomes difficult to evaluate data for *θ* < 200 ◦C, because the samples start to fracture, due to the increasing brittleness, especially when martensite is formed. Another impediment is the strain hardening at such temperatures. Narrowing down the temperature evaluation ranges in the tests fixes both problems, but involves an increased number of tests, which contradicts the original idea of this routine.

A visualization of *σ*0(*θ*) of Table 4 is presented in Figure 6. The results at heating show that the slope of *σ*0(*θ*) and *n*(*θ*) is small at around 300 ◦C, justifying to extrapolate both values to *θ*0.

Comparing *σ*<sup>0</sup> in the range of 200 ◦C < *θ* < 800 ◦C, at the same temperature, differences up to 500 MPa occur, which can be explained by the amount of austenite that is still present at 300 ◦C at sufficient quick cooling. Below 300 ◦C, the slope starts to increase as well, though. The values during cooling at 5 K/s on the other hand almost correspond to the heating values shifted by about −50 K.

Besides the optimzation results from the measurements, Figure 6 also contains the estimated *σ* 0 using Equation (5) for both cooling rates, which reproduces the experimentally determined curves very well. The parameters in use are listed in Table 5.

The hardness measurements greater than 730 HV confirm the martensite formation in the experiments. Therefore, *θ<sup>r</sup>* for cooling at 30 K/s is assumed to match the martensite start temperature. In this context it is necessary to mention, that CCT diagrams are sensitivite to fluctuations in the chemical composition as well as the specific test conditions for the identification of the CCT diagram. Compared to the tensile tests, the used CCT diagram in Figure 1 has been determined for a lower austenitizing temperature of 840 ◦C. Analogous to deformation CCT diagrams, special weld CCT diagrams exist, which are characterized by high austenitisation temperatures and short holding phases. The described conditions increase the transformation inertness during cooling, which leads to the formation of microstructures of higher hardness at lower cooling rates. Relative to conventional CCT diagrams for heat treatment, the transformation points in the weld CCT diagrams are shifted to the bottom right [21].

**Figure 6.** Estimated *σ* <sup>0</sup> based on *σ*<sup>0</sup> and parameters from Table 5.

### **5. Conclusions**

The flow properties of a 0.54% carbon steel are identified during continuous cooling and compared to results at heating. The findings are summarized as follows.


**Author Contributions:** Conceptualization, C.R., D.S.; methodology, C.R., D.S., T.H. and E.W.; software, C.R., D.S.; validation, O.S. and M.K.; formal analysis, C.R. and M.K.; investigation, C.R. and O.S.; resources, E.W. and S.J.; data curation, O.S.; writing—original draft preparation, C.R.; writing—review and editing, D.S., O.S. and M.K.; visualization, C.R.; supervision, T.H., S.J. and E.W.; project administration, C.R. and S.J.; funding acquisition, D.S. and S.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by AIF as part of IGF 18.966 B.

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

### **Abbreviations**

The following abbreviations are used in this manuscript:


### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **A Phenomenological Mechanical Material Model for Precipitation Hardening Aluminium Alloys**

### **Hannes Fröck 1, Lukas Vincent Kappis 1, Michael Reich 1,\* and Olaf Kessler 1,2**


Received: 25 September 2019; Accepted: 21 October 2019; Published: 29 October 2019

**Abstract:** Age hardening aluminium alloys obtain their strength by forming precipitates. This precipitation-hardened state is often the initial condition for short-term heat treatments, like welding processes or local laser heat treatment to produce tailored heat-treated profiles (THTP). During these heat treatments, the strength-increasing precipitates are dissolved depending on the maximum temperature and the material is softened in these areas. Depending on the temperature path, the mechanical properties differ between heating and cooling at the same temperature. To model this behavior, a phenomenological material model was developed based on the dissolution characteristics and experimental flow curves were developed depending on the current temperature and the maximum temperature. The dissolution characteristics were analyzed by calorimetry. The mechanical properties at different temperatures and peak temperatures were recorded by thermomechanical analysis. The usual phase transformation equations in the Finite Element Method (FEM) code, which were developed for phase transformation in steels, were used to develop a phenomenological model for the mechanical properties as a function of the relevant heat treatment parameters. This material model was implemented for aluminium alloy 6060 T4 in the finite element software LS-DYNA (Livermore Software Technology Corporation).

**Keywords:** aluminium alloy; EN AW-6060; precipitation hardening aluminium alloys; material model; heating; cooling; flow cures; LS-DYNA

### **1. Introduction**

Heat treatments are used to influence the mechanical properties of metallic materials. The material properties during and after heat treatment are dependent on different heat treatment parameters. Numerical heat treatment simulation offers the possibility of predicting process results, such as temperature gradients, mechanical properties, residual stress or distortion, in the heat-treated component. By using numerical simulation, process understanding can be further improved by simulating intermediate states, which are very difficult or impossible to measure experimentally. At the same time, the experimental effort for optimizing process parameters can be reduced by numerical simulation.

Precipitation-hardening aluminium alloys achieve their strengths through fine particles in the aluminium matrix. Different stable and metastable secondary phases can be formed or dissolved as a function of the temperature during the short-term heat treatment of precipitation-hardened states. Thus, the precipitation state, and therefore, the mechanical properties, change during short-term heat treatment. In heat treatment simulation, modelling the quenching of steels [1–4] and aluminium alloys [5–8] has been the subject of intensive research and is now used in many industrial applications [9]. In addition, numerical descriptions of the heat treatment processes are also used to simulate various manufacturing processes such as welding steels [10], aluminium alloys [11] or induction hardening [12]. Heat treatment simulation of the short-time heat treatments used for tailored heat-treating blanks (THTB) [13,14] and profiles (THTP) for aluminium has not yet been performed. These short-term heat treatments, which locally heat material areas using a laser, are characterized by very high heating rates of several 100 Ks<sup>−</sup>1, almost no soaking at peak temperature and subsequent cooling at a few 10 Ks−1. These short-term heat treatments can be used to intentionally soften local areas of semi-finished components and, thus, create a tailored strength layout. By tailoring the strength layout, the forming limits can be extended for subsequent forming. To further advance the development of THTP, a linked heat treatment and bending simulation can be used.

The mechanical properties of an aluminium alloy during a short-time heat treatment depend on the current temperature and the maximum heat treatment temperature. As a result, the mechanical properties during heating differ from those during cooling [15]. Currently, there is no material model for the mechanical properties of aluminium alloys depending on temperature and temperature path. For steels, this problem is solved using the mixture rules of different coarse phases like ferrite, pearlite, austenite, martensite, bainite, etc. [16]. For precipitation-hardened aluminium alloys with fine precipitates, this approach does not work at first glance. For aluminium alloys, some material models have also been developed. However, these are designed for special areas such as quench simulation [5–7] or ageing simulation [17,18]. However, these quench simulation material models can provide only temperature-dependent flow curves, while the ageing models only describe the change in mechanical properties during ageing. The mechanical properties as a function of the relevant parameters during short-time heat treatment cannot provide both types of models. Material models were also developed for producing tailored heat-treated semi-finished products [13,14]. However, these material models were developed for forming simulations and provide the mechanical properties at room temperature depending on the maximum temperature during heat treatment. The mechanical properties at elevated temperatures, as required in a heat treatment simulation, cannot be provided.

In this work, an efficient phenomenological material model is developed, which is based on the mixture rules of imaginary phases. Empirical-phenomenological model approaches provide simple descriptions of the material behavior as a function of the process parameters.

### **2. Materials and Methods**

### *2.1. Examined Aluminium Alloy*

The material model was developed using the precipitation hardening aluminium alloy EN AW-6060 in the natural aged state (T4). As a base material, a 20 mm × 20 mm × 2 mm hollow quadratic extrusion profile was used. The chemical composition of the investigated alloy is given in Table 1.


**Table 1.** Mass fraction of the alloying elements in the investigated alloy.

### *2.2. Database Used*

The precipitation and dissolution behavior of the alloy was recorded in previous work by direct [19] and indirect [20] differential scanning calorimetry (DSC) over a wide range of heating rates. The continuous time-temperature dissolution diagram of the alloy EN AW-6060 T4 can be derived from the DSC results. This diagram shows the temperature ranges at which a certain precipitation or dissolution reaction dominates depending on the heating rate. Figure 1 shows the continuous time-temperature dissolution diagram of the investigated alloy, EN AW-6060 T4.

**Figure 1.** Continuous heating dissolution diagram of EN AW-6060 T4 [20].

Figure 2 shows a DSC heating curve at 1 Ks−1, together with 0.2% yield strength and tensile strength after heating tensile specimens to different maximum temperatures followed by overcritical quenching [21]. The results indicate that the endothermic peak B, which is interpreted as the dissolution of clusters and GP-Zones [22], is accompanied by a significant decrease in strength. The exothermic peak c + d, which is considered to be the precipitation of the β" and β' phases [23], leads to a renewed increase in strength. The subsequent endothermic peak E, is considered to be the dissolution of the β" and β' precipitates [24], with further material softening.

**Figure 2.** Continuous heating DSC curve of EN AW-6060 T4 at 1 Ks−<sup>1</sup> correlated with the 0.2% yield and tensile strength after heating tensile specimens to different maximum temperatures followed by overcritical quenching. Experiments described in [21].

Short-term laser heat treatment, as well as recording time-temperature profiles, was carried out at the Institute of Manufacturing Technology of the University of Erlangen-Nürnberg [25]. Laser heat treatment is characterized by a high heating rate of up to several 100 Ks−1, with no soaking at the maximum temperature and a relatively slow, non-linear, cooling with a few 10 Ks<sup>−</sup>1. The investigated time-temperature courses are shown in Figure 3.

**Figure 3.** Recorded temperature courses of different laser heat treatments, described in [25].

These time-temperature profiles were imitated in a quenching and deformation dilatometer, interrupted at defined temperatures and tensile tests were carried out immediately at these temperatures in the same device [15]. The schematic measurement plan of the thermo-mechanical analysis, with the major parameters of the tensile test, is shown in Figure 4.

**Figure 4.** Schematic measurement plan of the thermo-mechanical analysis, with the major parameters of the tensile test. Described in more detail in [15].

It can be seen in Figure 5 that the strength of the alloy decreases during heating with increasing temperature. During subsequent cooling, the strength increases again with decreasing temperature. The mechanical properties during heat treatment with a 180 ◦C maximum temperature are identical during heating and cooling and only dependent on temperature, as shown in Figure 5A. Up to this temperature, no permanent softening of the material has taken place. In short-term heat treatment with a 250 ◦C maximum temperature, the strength during cooling does not increase to the same extent as it decreased during heating. The mechanical properties at a certain temperature are, therefore, not identical during heating and cooling, as shown in Figure 5B. The same behavior is evident in the short-term heat treatments at 300 ◦C and 400 ◦C maximum temperatures, as seen in Figure 5C,D. For a purposeful simulation of short-term heat treatment, these mechanical properties dependencies must be implemented in a material model.

**Figure 5.** Mechanical properties (0.2% yield and tensile strength) during various short-term heat treatments. Experiments described in [15].

### **3. Material Model Development**

### *3.1. The Basic Idea*

We developed a phenomenological mechanical material model for precipitation hardening aluminium alloys, which can describe flow curves as a function of the actual temperature and the peak temperature of the heat treatment. The basic idea was to define the precipitation-hardened initial state as an imaginary hardened phase (A) and the state in which the strength-increasing precipitates were dissolved as an imaginary softened phase (B). These individual phases can be assigned temperature-dependent flow curves.

From a defined starting temperature (Tstart), the imaginary hardened phase (A) is transformed into the imaginary softened phase (B) as a function of the temperature during heating. This transformation is completed at a defined finish temperature (Tfinish). At temperatures between the start and finish temperature, a phase mixture of an imaginary hardened and an imaginary softened phase exists, as shown in Figure 6. At temperatures above the final temperature, the material consists of the softened phase (B). Through developing the material model, it has been found that the accuracy of the model can be increased by introducing an intermediate phase (Z) between the hardened phase (A) and the softened phase (B), see Figure 7. During cooling, no phase transformations are permitted. The resulting mechanical properties of the material state are defined by linearly mixing the mechanical properties of the occurring phases. The schematic of the implemented phase transformation during heat treatment is shown in Figure 7.

**Figure 6.** Schematic approach for an empirical-phenomenological material model of precipitation hardening aluminium alloys.

**Figure 7.** Schematic of the implemented phase transformation during heat treatment.

The following simplifications are assumed in the model:


### *3.2. Determination of the Phase Transformation Temperatures*

To adapt the material model to the investigated alloy, the phase transformation temperature ranges must be determined. This was performed based on the continuous heating dissolution diagram (Figure 1) and mechanical results (Figure 5). Figure 5A shows that the mechanical properties of the heating and cooling step did not differ for short-term heat treatment up to a peak temperature of 180 ◦C. Figure 1 shows that the clusters and GP-zones of the initial state start to dissolve between 200 ◦C and 225 ◦C at a high heating rate of 100 Ks<sup>−</sup>1. Thus, it was determined that the imaginary hardened phase (A) was completely present up to a maximum temperature of 212.5 ◦C. At higher temperatures, phase (A) continuously transformed into the intermediate phase (Z). At a maximum temperature of 250 ◦C,

the material consisted only of the imaginary intermediate phase (Z). It can be seen in Figure 1 that at a heating rate of 100 Ks<sup>−</sup>1, the dissolution of the cluster and GP-zones was completed at 300 ◦C. This is also reflected in the mechanical properties, as shown in Figure 5C,D. For this reason, it is assumed that the simulated phase transformation Z to B is completed at 300 ◦C. The properties above this peak temperature are described by phase (B). Figure 8 shows the implemented temperature transformation ranges of the imaginary phases.

**Figure 8.** Schematic temperature transformation ranges of the imaginary phases for alloy 6060 T4.

### *3.3. Calculating the Resulting Flow Curves*

The basis for the mechanical properties of individual phases is experimental stress–strain curves after different short-term heat treatments (Figure 5). For the thermo-mechanical simulation, flow curves are required as input data (true stress; *kf* vs. logarithmic plastic strain; ϕ). Experimentally obtained flow curves were converted into a temperature-dependent mathematical description. The flow curves of face-centered cubic metals, like aluminium alloys, can be described well by the Hockett–Sherby [26] relationship, see Equation (1).

$$k\_f(qr) = k\_s - \varepsilon^{-mq^{\nu}}(k\_s - k\_0) \tag{1}$$

This equation contains four parameters, the initial flow stress (*k*0), the saturation stress (*ks*) and the hardening exponents, *m* and *P*. These parameters were adapted to the experimentally determined flow curves, using OriginPro 2018. By adapting the parameters to the flow curves at different temperatures, the Hockett–Sherby parameters of a phase can be derived over the existence range of this phase.

The Hockett–Sherby parameters of the imaginary phase (B) can be obtained over the entire temperature range directly from the experimental flow curves. The experimental basis for the temperature-dependent flow curves of the imaginary softened phase (B) is the experimentally obtained flow curves during cooling after the 400 ◦C maximum temperature, see Figure 9, and an additional tensile test at 550 ◦C. The flow curves are shown in the following figures as true stress *kf* vs. logarithmic plastic strain ϕ. The points in Figure 10 show the Hockett–Sherby parameters for phase (B), which were derived from the experimental data. These parameters from the experimental data can be fit via linear or nonlinear mathematical functions.

0 5 10 0 50 100 150 200 kf in Nmm-2 M in % 25 °C 175 °C 225 °C 275 °C 325 °C 375 °C 400 °C Experimentally obtained flow curves of EN AW-6060 T4 at different temperatures during cooling after TMax: 400 °C

**Figure 9.** Experimentally obtained flow curves of EN AW-6060 T4 alloy at different temperatures during cooling after the 400 ◦C maximum temperature, phase (B).

**Figure 10.** Temperature-dependent Hockett–Sherby parameters of phase (B).

In the temperature range between 212.5 ◦C and 300 ◦C, the mechanical properties were represented by an imaginary phase mixture, A + Z respectively Z + B, as shown in Figure 8. The resulting flow curve (*kf*) can be described using Equation (2), taking into account the phase proportions of the existing phases (*xi*) and the flow curves of the individual phases (*kfi*). The phase proportions were calculated via a phase transformation model as a function of temperature (Section 3.3), while the individual phases were assigned to temperature-dependent flow curves.

$$k\_f = \sum\_{i=1}^{n} \mathbf{x}k\_{fi} \tag{2}$$

The resulting flow curves during the heating of the imaginary phase mixture A + Z respectively Z + B can be determined experimentally by the tensile tests performed during heating. Figure 11 shows the experimentally obtained flow curves of the imaginary phase mixture A + Z respectively Z + B at various temperatures during heating. The points in Figure 12 show the Hockett–Sherby parameters while heating alloy EN AW-6060, which were derived from the experimental data. The points were fit by mathematical formulas, which allows the Hockett–Sherby parameters to be determined for all temperatures during heating.

Experimentally obtained flow curves of EN AW 6060-T4 at different temperatures

**Figure 11.** Experimentally obtained flow curves of imaginary phase mixture A + Z respectively Z + B for alloy EN AW-6060 T4 at different temperatures during heating at 100 Ks<sup>−</sup>1.

**Figure 12.** Temperature-dependent Hockett–Sherby parameters of imaginary phase mixture A + Z respectively Z + B during heating.

The simulated microstructure after a 250 ◦C maximum temperature consists entirely of phase (Z). The temperature-dependent flow curves and the Hockett–Sherby parameter of phase (Z) can thus be determined from tensile tests during cooling after a short-time heat treatment with a 250 ◦C maximum temperature, as seen in Figures 13 and 14.

As can be seen in Figure 8, phase (Z) is transformed into phase (B) up to a 300 ◦C maximum temperature. Thus, temperature-dependent flow curves of phase (Z) up to 300 ◦C are necessary. By definition, the experimentally determined flow curves during heating between 250 ◦C and 300 ◦C represent the resulting mechanical properties of phase mixture (Z) and (B). As already described, the temperature-dependent flow curves of phase (B) can be calculated over the entire temperature range. The flow curves of phase (Z) in the interval between 250 ◦C and 300 ◦C can thus be determined by the experimentally determined resulting flow curves (*kfres*) from the phase mixture of (Z) and (B) during heating, as well as the experimentally determined flow curves of phase (B) (*kfB*) using Equation (3).

$$k\_{fZ}(T) = \frac{k\_{f\_{ra}}(T) - \mathbf{x}\_{\mathcal{B}}(T) \cdot k\_{fB}(T)}{\mathbf{x}\_{\mathcal{Z}}(T)}\tag{3}$$

The phase proportions of the phases are given by the temperatures and can be calculated from the phase transformation (Section 3.3). Figure 14 shows the Hockett–Sherby parameters determined for phase (Z) in its entire existence temperature range from 300 ◦C to room temperature.

Experimentally obtained flow curves of EN AW-6060 T4 at different temperatures

**Figure 13.** Experimentally obtained flow curves of alloy EN AW-6060 T4 at different temperatures during cooling after a 250 ◦C maximum temperature (phase Z).

**Figure 14.** Temperature-dependent Hockett–Sherby parameters of phase (Z).

The experimental database for phase (A) provided the tensile tests performed during the short-time heat treatment with a 180 ◦C maximum temperature, see Figure 15. The same procedure for determining temperature-dependent Hockett–Sherby parameters (Figure 16) was performed for phase (A) from the imaginary phase mixture A + Z.

**Figure 15.** Experimentally obtained flow curves of alloy EN AW-6060 T4 at different temperatures during cooling after a 180 ◦C maximum temperature (phase A).

**Figure 16.** Temperature-dependent Hockett–Sherby parameters of phase (A).

The presented procedure makes it possible to adjust the Hockett–Sherby parameters for each phase. As a result, flow curves can be calculated for each phase at any temperature in the temperature range of phase existence.

### *3.4. Phase Transformation*

According to the model in Figure 6, we have chosen simplified linear transformations of imaginary phases A → Z and Z → B in the specified temperature ranges. Corresponding to the continuous heating dissolution diagram in Figure 1, we have further assumed no dependence of phase transformations on heating rate in the relevant range of approximately 10 Ks−<sup>1</sup> to 100 Ks−1. Figure 17 shows the phase fractions during a simulated heating up to 400 ◦C. It becomes clear that an entire linear transformation of the individual phases was achieved as a function of temperature, considering the defined start and end temperatures of the transformations.

**Figure 17.** Simulated phase transformations at different heating rates up to 400 ◦C.

### *3.5. Simulation Model*

The simulation model was implemented using the finite element software LS DYNA (Livermore Software Technology Corporation). The material model MAT\_GENERALIZED\_PHASECHANGE was released in software version R9.0.1 of LS DYNA (August 2016). It is used to model phase transformations in metallic materials and the associated changes in material properties. Up to 24 individual phases that transform into each other can be defined. The transformation can be defined for heating or cooling [27]. The phase transformation laws according to Koistinen–Marburger, Johnson–Mehl–Avrami–Kolmogorov (JMAK), and Kirkaldy and Oddy are predefined in the keywords.

A simple geometric model was used to verify the phenomenological mechanical material model for precipitation hardening aluminium alloys. The simulation model consisted of a single cube-shaped solid element with a 1 mm edge length. This element was subjected to a defined temperature profile to simulate heat treatment. The element was firmly clamped on one side by boundary conditions. The nodes on the other side were loaded with a predetermined displacement. Thermal expansion was not considered during this simulation to prevent thermally induced deformation.

Figure 18 illustrates an example time-temperature profile and the time-displacement course used. The time-temperature curve was varied in both the maximum temperature and in the temperature during deformation, in different individual simulations. The time-displacement curve remained constant during all individual simulations. The simulated flow curves were then compared with experimental data.

**Figure 18.** Exemplary time-temperature and time-displacement graph for checking the simulated transformation behavior and the resulting properties.

In the keywords in LS-DYNA, the Hockett–Sherby parameters were not entered directly. The representation of the flow curves via the Hockett–Sherby parameters serves only to calculate an associated flow curve for each phase at any temperature. In LS-DYNA, the flow curves were assigned to a temperature in tabular form for each phase. In the presented material model, a flow curve was tabulated every 25 K for each phase.

### **4. Results**

The experimental and the simulated flow curves during different short-time heat treatments of alloy EN AW-6060 are shown in Figure 19. The experimentally determined flow curves shown represent the average of at least three individual measurements. The true stress over the plastic strain was plotted. The flow curves clearly show that our phenomenological mechanical material model for precipitation hardening aluminium alloys agrees very well with the experiments.

Figure 20 highlights some applications of the model. Figure 20A shows that the material model can provide temperature-dependent flow curves during heating. Strength decreases with increasing temperature, comparable to Figure 5. Figure 20B shows that the simulated flow curves during heating and cooling run differently. By heating to 275 ◦C, the material is softened. The flow curves after cooling to 200 ◦C and 25 ◦C are therefore significantly below the flow curves from heating to the same temperatures. Figure 20C shows that the flow curves at the same current temperature (here 25 ◦C) depend on the maximum temperature of the previous short-time heat treatment. A short heating to 200 ◦C does not cause softening and conforms to the flow curve of the initial state. Up to 300 ◦C, the material is softened, and the strength decreases with the increasing maximum temperature. Above the 300 ◦C maximum temperature, the material is maximally softened, and the flow curves do not change any further.

Measured and simulated flow curves according to different short-time heat treatments of alloy EN AW-6060 T4

**Figure 19.** Measured and simulated flow curves of alloy EN AW-6060 T4 during different heat treatments.

**Figure 20.** Simulated flow curves of alloy EN AW-6060 T4 during various short-time heat treatments. (**A**) During heating to different temperatures; (**B**) Distinction heating/cooling to 25 ◦C and 200 ◦C after *TMax*: 275 ◦C; (**C**) At 25 ◦C according to different maximum temperatures.

The developed material model can provide flow curves that depend on both the current temperature and the previous maximum temperature, and distinguish between heating and cooling. The developed material model thus takes into account the essential influencing factors of a short-term heat treatment on the mechanical properties and is therefore, suitable for the coupled thermal, metallurgical, and mechanical simulation of such a heat treatment. In this case, metallurgical does not mean the real nm-sized precipitation processes in Al-alloys, but a feasible macroscopic approach based on imaginary phases.

### **5. Conclusions**

The two main influencing factors for softening aluminium alloys during short-term heat treatment are the current temperature and the previous maximum temperature [21]. Until now, no material model can describe the mechanical properties as a function of the current temperature and the previous maximum temperature during the short-term heat treatment and can also differentiate between heating and cooling. In this work, a phenomenological mechanical material model was developed for precipitation hardening aluminium alloys. In this model, the real hardening behavior of aluminium alloys, which depends on nm-sized precipitates, was not considered. Instead, the mechanical properties were defined by a mixture of different imaginary phases comparable to coarse phase mixtures in steels. The initial precipitation-hardened state was defined as an imaginary hardened phase (A). During short-term heat treatment, precipitates were dissolved, and the material was softened. This state was defined as an imaginary softened phase (B). Between the hardened phase (A) and the softened phase (B), an imaginary intermediate phase (Z) was defined. These phases transformed into each other during heating depending on the temperature. During cooling, no further phase transformation was assumed. The individual phases were assigned to temperature-dependent flow curves. Due to the phase composition as a function of the maximum temperature and the temperature-dependent flow curves of the individual phases, the developed material model can calculate a flow curve for each current temperature and each previous maximum temperature.

The material model was developed for calculating the mechanical properties during typical short-term heat treatments. For this reason, the following factors deviating from short-term heat treatment conditions were not considered in the material model—heating and cooling rate, strain rate and subsequent natural ageing. In principle, these factors can be considered in the model, but need an even broader experimental database.

The presented material model with the imaginary phases can be adapted to different age hardening aluminum alloys. For this purpose, the flow curves of the corresponding phases, as well as the transformation temperatures of the individual phases must be adapted to the alloy.

The model was implemented in the finite element software LS-DYNA using the keyword MAT\_GENERALIZED\_PHASECHANGE. The laser short-time heat treatment of Al-alloys can be simulated with this material model. Residual stresses and distortions can be calculated and handed over to a subsequent forming simulation, to realize a through process simulation.

**Author Contributions:** H.F., L.V.K., M.R. and O.K. conceived and designed the experiments; H.F. performed the experiments and analysed the data, L.V.K. implemented the model; H.F., L.V.K., M.R., and O.K. discussed and interpreted the results together; H.F. wrote the paper.

**Funding:** This research was funded by the German Research Foundation (DFG), within the scope of the research project Improvement of formability of extruded aluminium profiles by a local short-term heat treatment (DFG KE616/22-2).

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

### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
