**Numerical Simulation of Three-Dimensional Mesoscopic Grain Evolution: Model Development, Validation, and Application to Nickel-Based Superalloys**

### **Zhao Guo, Jianxin Zhou \*, Yajun Yin, Xu Shen and Xiaoyuan Ji**

State Key Laboratory of Materials Processing and Die & Mould Technology, Huazhong University of Science and Technology, Wuhan 430074, China; zg2014@hust.edu.cn (Z.G.); Yinyajun436@hust.edu.cn (Y.Y.); shenxu@hust.edu.cn (X.S.); jixiaoyuan@hust.edu.cn (X.J.)

**\*** Correspondence: zhoujianxin@hust.edu.cn; Tel./Fax: +86-027-87541922

Received: 27 November 2018; Accepted: 5 January 2019; Published: 9 January 2019

**Abstract:** The mesoscopic grain model is a multiscale model which takes into account both the dendrite growth mechanism and the vast numerical computation of the actual castings. Due to the pursuit of efficient computation, the mesoscopic grain calculation accuracy is lower than that of dendrite growth model. Improving the accuracy of mesoscopic grain model is a problem to be solved urgently. In this study, referring to the calculation method of solid fraction in microscopic dendrite growth model, a cellular automata model of 3D mesoscopic grain evolution for solid fraction calculated quantitatively at the scale of cell is developed. The developed model and algorithm validation for grain growth simulation is made by comparing the numerical results with the benchmark experimental data and the analytical predictions. The results show that the 3D grain envelopes simulated by the developed model and algorithm are coincident with the shape predicted by the analytical model to a certain extent. Then, the developed model is applied to the numerical simulation of solidification process of nickel-based superalloys, including equiaxed and columnar dendritic grain growth. Our results show good agreement with the related literature.

**Keywords:** numerical simulation; cellular automaton; dendritic grain growth; quantitative prediction

### **1. Introduction**

Grain evolutions including nucleation and growth during solidification process have a profound effect on the material mechanical properties of alloys. In particular, the competitive growth of equiaxed and columnar dendritic grains determine the microstructure of alloys. Over the past few decades, many computer models and methods have been developed to study theoretically and experimentally the grain growth in alloy solidification, such as the volume-averaged, front tracking, cellular automata (CA), and phase field models [1–5]. Among them, due to the advantages of clearly describing the physical phenomena and high computational efficiency of solidification transition, the CA model has been widely used and developed for prediction of grain growth in additive manufacturing [6,7], welding [8,9], and casting [10,11].

In 1993, the CA model was applied for the first time in alloy solidification by Rappaz and Gandin [12]. It mainly dealt with nucleation and growth of grains in a uniform temperature field. At this point, the basic framework of the CA model for grain growth simulation was constructed, although it still had some drawbacks. Subsequently, in order to satisfy the adaptability of non-uniform temperature field at the same time, the "2D rectangular algorithm" [13] and "decentered square algorithm" [14,15] coupled with finite element (FE) calculations were proposed. Similarly, the CA model could also be coupled with the finite difference (FD) method to form the CA-FD model [16],

or coupled with the finite volume (FV) method to form the CA-FV model [17]. As one of the effective methods to simulate the growth of grain, the CA-FE model integrated into the casting commercial software of ProCAST have been widely used [18,19].

However, the growing octahedron in these models is no longer considered if the corresponding cell is fully surrounded by the mushy cells [15]. In this way, this CA-FE model is still a weak coupling method and overestimates the heat release of cell transformation in the mushy zone [20]. In the CA model, the solid fraction of growing cell or interface cell in the mushy zone cannot be calculated quantitatively at the scale of CA grid. And it is calculated approximately at the scale of FE grid by the Scheil's equation and the Scheil microsegregation model. Recently, Guillemot et al. [20] have employed a front tracking method proposed by Gandin [21] to deal with the coupling problem. Similar works have also been conducted by Carozzani et al. [22] and Liu et al. [23,24].

Additionally, the CA model of grain growth, ignoring the dendritic morphology of the grain itself, is generally referred to as mesoscopic model. Corresponding CA models of dendritic growth called microscopic model have been proposed and developed, typically works like Nastac's [25], Zhu's [26] and Pan's [27]. In this regard, it is a full coupling model that accounts for the time evolution of the fraction of solid predicted at the scale of CA grid in the CA model of dendritic growth. It is further observed in the CA models of two scales that both are faced with the same problems of grid anisotropy [13,28], transformation rules [14,29], and coupling methods [30]. From the point of view of the simulation method, the method for dealing with the same problem in the CA models of two scales are interrelated [16,31]. Nevertheless, there is a problem that the CA model of dendritic growth requires sufficiently higher mesh resolution than the CA model of grain growth [32].

The purpose of this paper is to achieve the full coupling simulation of grain growth so as to improve the simulation accuracy. A 3D modified cellular automaton coupling model of grain growth for solid fraction calculated quantitatively at the scale of CA gird is proposed. The relevant models and algorithms are presented and validated. Then, the developed model is applied to simulate the formation of equiaxed and columnar dendritic grains during solidification process of nickel-based superalloys.

### **2. Model Description and Numerical Algorithm**

A 3D computation domain is firstly divided into a series of regular hexahedral cells. The cell spacing is equal to the corresponding mesh size. Each cell is characterized with such variables as solid fraction, grain growth rate, grain orientation, etc. All the cell states are defined as three types: liquid (*fs* = 0), solid (*fs* = 1), or interface (0 < *fs* <1). And the cell neighbors are characterized by the first-order neighbors (6 neighborhoods), second-order neighbors (12 neighborhoods), and third-order neighbors (8 neighborhoods). At the beginning of the numerical simulation, the cells in computation domain are initialized with given initial and boundary conditions. When the undercooling reaches or exceeds the critical nucleation undercooling, the grain begins to nucleate and grow. As the dendritic grain grows, the solid fraction of growing cell or interface cell is calculated quantitatively at the scale of CA grid. The governing equations for calculating the nucleation distribution and growth rate of dendritic grains are presented. The numerical algorithm about capture rules for new interface cells and the calculation of solid fraction at the scale of CA grid are described in detail below.

### *2.1. Grain Nucleation*

Only heterogeneous nucleation is considered in this study, since high undercooling is not achieved due to the presence of foreign phases. For the sake of practicability, the continuous nucleation model proposed by Rappaz and Gandin et al. [12] is used to describe the grain nucleation. The total grain density *n*(Δ*T*) is calculated by:

$$m(\Delta T) = \int\_0^{\Delta T} \frac{dn}{d\Delta T'} d\left(\Delta T'\right) \tag{1}$$

where ΔT is the undercooling.

The continuous nucleation distribution *dn*/*d*Δ*T* is given by:

$$\frac{dn}{d\Delta T'} = \frac{n\_{\text{max}}}{\sqrt{2\pi}\Delta T\_{\sigma}} \exp\left[-\frac{1}{2}\left(\frac{\Delta T - \Delta T\_N}{\Delta T\_{\sigma}}\right)\right] \tag{2}$$

where *n*max is the total density of grains, Δ*T<sup>σ</sup>* is the standard deviation of undercooling, and Δ*TN* is the mean undercooling in nucleation, respectively. Additionally, the total density of grains on the two-dimensional section (*nmax*)2D can be obtained by the stereological relationships [12], (*nmax*)2D = (6·*nmax*2/π) 1/3.

Unlike the numerical realization of nucleation process determined by random number method, in this present study, the nucleation position and critical undercooling are randomly assigned to the nucleation region in advance, and then the actual undercooling and the critical undercooling of each cell are compared in each time step to determine the nucleation.

### *2.2. Grain Growth*

Grain growth rate is mainly determined by the growth behavior of dendrite tips. According to the Lipton-Glicksman-Kurz (LGK) growth model [33,34], the dendrite tips growth of alloys are affected by thermal undercooling (Δ*Tt*), constitutional undercooling (Δ*Tc*) and curvature undercooling (Δ*TR*). With the growth of dendrite tips, latent heat and solute are released from tips. At the interface of dendrite tips, heat and solute transfer lead to changes in thermal and constitutional undercooling. In addition, the growth of dendrite tips also lead to the change of curvature undercooling. At a given total undercooling, the grains will grow under the interaction of various factors until they grow up completely or encounter obstacles to stop growth.

The thermal undercooling, constitutional undercooling and curvature undercooling are defined as:

$$
\Delta T\_t = I \upsilon(P\_t) \frac{\Delta H}{c\_p} \tag{3}
$$

$$
\Delta T\_c = m\_L \mathbb{C}\_0 \left( 1 - \frac{1}{1 - (1 - k\_0)Iv(P\_\natural)} \right) \tag{4}
$$

$$
\Delta T\_R = \frac{2\Gamma}{R} \tag{5}
$$

where *Iv*(*P*) is the Ivantsov function of the Peclet number; *Pt* = *RV*/(2*α*) and *Pc* = *RV*/(2*D*) are the heat Peclet number and solutal Peclet number, respectively; Δ*H* is the latent heat; *cp* is the specific heat; *mL* is the liquidus slope; *C*<sup>0</sup> is the initial components concentration of alloys; *k*<sup>0</sup> is the solute partition coefficient; Γ is Gibbs-Thompson coefficient; *R* is the tip radius; *V* is the interface growth rate; *α* is thermal diffusion coefficient; and *D* is solute diffusion coefficient.

The 3D Ivantsov function is expressed as:

$$I\upsilon(P) = P\exp(P)E\_1(P) \tag{6}$$

where *E*1(*P*) is the exponential integral function.

The relationship between the tip radius and the interface growth rate can be established by combining the above equations, but it cannot be solved, and other more definite conditions are needed.

Based on the stability criterion of a dendrite tip, the tip radius is equal to the shortest wavelength *λi*, which can be evaluated by:

$$
\lambda\_i = \left(\frac{\Gamma}{\sigma^\*(m\_L\mathbf{G}\_c - \mathbf{G})}\right)^{0.5} \tag{7}
$$

where *Gc* is the concentration gradient; *G* is the temperature gradient; *σ\** is a stability parameter, it should be noted here that the stability parameter according to the microsolvability theory varies

with the degree of anisotropy of the interface energy, and its value can be obtained by the linearized solvability theory [35].

For ease of numerical calculations, the above equations are realized and solved using Matlab. The growth rate is simplified as a polynomial function of local undercooling.

### *2.3. Capture Rules for New Interface Cells and Calculation of Solid Fraction at the Scale of the CA Grid*

A single regular octahedron is defined as a constituent element of a grain. Each cell corresponds to a regular octahedron. The sizes of regular octahedrons are characterized by the distance (*SL*) between any two adjacent vertices. Its three orthogonal axes ([100]/[010]/[001]) represent the directions of dendrite tips. Meanwhile, in order to avoid the infinite growth problem of the single regular octahedron, the distance (*SL*) is limited to 2·Δ*x* (Δ*x* is the cell size). The grain envelope, which consists of a number of regular octahedrons, is also shown as a larger regular octahedron in the uniform temperature field. A single grain whose crystallographic orientation is assumed to be the Euler angles (*ϕ*, *θ*, *ψ*) with respect to the axis of CA, as shown in Figure 1.

**Figure 1.** Schematic diagram of a single regular octahedron with arbitrary orientation (*ϕ*, θ, *ψ*).

In this present work, the nucleated cell is defined as a solid cell. And all its neighbor cells are defined as interface cells with a solid fraction of *ε* (*ε*→0). When the solid fraction of the interface cells reaches 1 or greater than 1, all the first-order, second-order, and third-order neighbors' liquid cells of the former interface cell transform to interface cells. At this point, it is guaranteed that at least one of the 26 neighbor cells of the interface cells evolves to a solid cell. The Schematic diagram of cell state is shown in Figure 2. For the octahedron corresponding to the nucleated cell, the center of octahedron coincides with the center of nucleated cell. While the center of octahedron corresponding to the interface cell also coincides with the center of nucleated cell. Similar to 3D decentered square algorithm, the tips length of octahedron (*L*) from its center to any of its six vertices are given by:

$$L(t\_2) = L(t\_1) + \int\_{\tilde{t}\_1}^{t\_2} V(\Delta T(\tau))d\tau \tag{8}$$

where *t*<sup>2</sup> represents the next time, *t*<sup>1</sup> represents the current time. The initial tip length of the grain corresponding to the nucleated cell is taken as the maximum value (√2·Δ*x*). The initial tip length of the grain corresponding to the interface cell is taken as ε (ε→0).

**Figure 2.** Schematic diagram of the cell state.

During the evolution of grains, due to the displacement between the center of the interface cell and the center of the growing octahedron, the growing octahedron corresponding to the interface cell will capture the cell center of the neighbors of solid cell. When the interface cell is captured, it is then switched to the solid cell. At this moment, all the remaining liquid cells being neighbors of this new solid cell will be automatically transferred from liquid to interface cells. The new octahedron corresponding to this new interface cell is initialized. Once the tip length of the growing octahedron is equal to or greater than the maximum value (√2·Δ*x*), the growing octahedron will cease to grow. The corresponding interface cell will be switched to the solid cell regardless of the solid fraction of the interface cell.

In order to realize the self-consistency between the tip length of the growing octahedron and the solid fraction of interface cell, a quantitative method for calculating the two is described below.

Figures 3 and 4 sketch the capture rule for solid cell in a time step calculation. When the center of growing octahedron corresponding to the interface cell (for example, IJ) coincides with the center of nucleated cell, the fraction increment of the interface cell IJ is calculated by:

$$
\Delta f\_s = \frac{\Delta L}{L\_1} with \qquad L\_1 \le \sqrt{2} \Delta x \tag{9}
$$

where Δ*L* represent the increment of the tip length of the growing octahedron, *L*<sup>1</sup> represent the final tip length (OA1) of the growing octahedron.

**Figure 3.** Schematic diagram of the capture rule for solid cell in a time step calculation, the center of IJ is in the final plane A1D1E1 of the growing octahedron with *<sup>L</sup>*<sup>1</sup> <sup>≤</sup> <sup>√</sup>2·Δ*x*.

Subsequently, as soon as the center of the interface cell IJ is captured by the growing octahedron, this interface cell is immediately switched to solid cell and the remaining liquid cells adjacent to this new solid cell will be transformed into new interface cells (for example, IJK), as shown in Figure 4. As mentioned above, the new growing octahedron corresponding to the new interface cell IJK is initialized. Especially, the center (O0) of the new growing octahedron is also located at the center of the octahedron corresponding to the original interface cell IJ. And the tip length of the new growing octahedron is initialized with the maximum value *<sup>L</sup>*<sup>m</sup> (*L*<sup>m</sup> <sup>=</sup> <sup>√</sup>2·Δ*<sup>x</sup>* <sup>−</sup> *<sup>ε</sup>*). As a result, before the solid fraction of the new interface cell IJK reaches 1, the center (O2) of the new growing octahedron has moved gradually along the direction of growth so as to maintain its continuous growth. In this case, the fraction increment of the new interface cell IJK is calculated by:

$$
\Delta f\_{\circ} = \frac{\Delta L}{L\_2 - L\_m} \tag{10}
$$

where *L*<sup>2</sup> represent the final tip length (O0A3) of the growing octahedron.

**Figure 4.** Schematic diagram of the capture rule for solid cells in a time step calculation, the center of IJK is in the final plane A3D3E3 of the growing octahedron with *Lm* <sup>=</sup> <sup>√</sup>2·Δ*x*−*ε*.

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

### *3.1. Comparison with the Benchmark Experimental Data*

Firstly, to validate the numerical calculation program for mesoscopic grain growth in this study, the dendritic tip growth rate of a succinonitrile-acetone mixtures are calculated and compared with benchmark experimental data [34]. According to the microsolvability theory, the stability parameter remains unchanged within a certain range of Peclet number. Here, two stability parameters are considered. Figure 5 presents the dendritic tip growth rate as a function of initial concentration in succinonitrile-acetone mixtures.

**Figure 5.** Dendritic tip growth rate as a function of initial concentration in succinonitrile-acetone mixtures.

As shown in Figure 5, under both undercooling conditions of 0.5 K and 0.9 K, the dendrite tip growth rate calculated by numerical methods increases first and then decreases with the initial concentration, which is consistent with the experimental results. For the initial concentration range of 0–0.25 (mol. %), the curve of the dendritic tip growth rate with a stability parameter of 0.025 is approximately matched with the experimental results. For the initial concentration range of 0.25–0.45 (mol%), the curve of the dendritic tip growth rate with a stability parameter of 0.0195 is approximately matched with the experimental results. It is confirmed that the stability parameters match the experimental results only in a certain range of Peclet numbers.

### *3.2. Comparison with the Analytical Model of Grain Growth*

The model and algorithm validation for grain growth simulation is made by comparing the numerical results with the analytical model predictions [14,36]. In the analytical growth model, it was assumed that the grain growth was controlled by the solute diffusion and only applicable to the low Peclet number regime. This is basically similar to the assumption of the LGK model. There are two kinds of grain branches predicted by the analytical model, here the external envelope is assumed to be the grain envelope. Figure 6 presents the schematic diagram of two-dimensional dendrite growth envelope trajectory under directional temperature. In the process of two-dimensional dendritic growth under directional solidification, the coordinates of the external envelope trajectory can be calculated by Equation (11) [36]:

$$\begin{cases} \begin{aligned} \mathbf{x} &= \mathbf{x}\_{0} + (y - y\_{0}) \cdot f(\theta) \\\\ \Delta T &= \frac{\Delta T\_{L}}{\sqrt{\mathbf{g}(\theta)}} \cdot \frac{h\left(\mathbf{G} \cdot \mathbf{A} \cdot \Delta T\_{L} \cdot \sqrt{\mathbf{g}(\theta)} \cdot (t - t\_{0})\right) + \frac{\Delta T\_{\overline{0}} \cdot \sqrt{\mathbf{g}(\theta)}}{\Delta T\_{L}}}{\Delta T\_{L}} \\\\ \mathbf{y} &= \mathbf{v}\_{L} \cdot t - \frac{\Delta T}{\mathbf{G}} \end{aligned} \end{cases} \tag{11}$$

where *VL* is the moving rate of liquid isotherm, Δ*TL* is the dendrite tip undercooling, Δ*T*<sup>0</sup> is the undercooling at initial time *t*0, *θ* is the dendrite orientation angle, *A* is a constant coefficient in LGK model, *f*(θ), *g*(θ) and *S*<10> are directional functions related to the orientation angle respectively, which are given by Equations (12)–(14).

**Figure 6.** Schematic diagram of two-dimensional dendrite growth envelope trajectory under directional temperature.


$$\begin{cases} \ S\_{\langle 10 \rangle} = 1 & \text{for the } [10] \text{ and } [01] \text{ branches} \\ \ S\_{\langle 10 \rangle} = -1 & \text{for the } [\overline{1}0] \text{ and } [0\overline{1}] \text{ branches} \end{cases} \tag{14}$$

Since the above analytical model is based on two-dimensional dendrite growth, it is necessary to extend this method to three-dimensional grain growth. In this present study, the projection method and the three cosine theorem are used to calculate the outermost envelope trajectory of a three-dimensional grain. Besides, since the analytical model was restricted to a quadratic approximation of the velocity-undercooling relationship, the growth rate was temporarily simplified as the same quadratic approximation relationship of local undercooling in the numerical simulation. Figure 7 presents the 3D envelopes and tips length of a single grain growing in a uniform temperature field and in a Bridgman-like situation. In Figure 7a, the cooling rate is set to 0 K/s and the thermal gradient is set to 0 K/m, the tips length (*L*1, *L*2, *L*3, *L*4, *L*5, *L*6) of the grain corresponding to Figure 7a are shown in Figure 7b. In Figure 7c, the cooling rate is set to −0.1 K/s and the thermal gradient along the *Z*-axis is set to 250 K/m, the tips length (*L*1, *L*2, *L*3, *L*4, *L*5, *L*6) of the grain corresponding to Figure 7c are shown in Figure 7d. The cell size Δ*x* is set to 0.0001 m.

*Metals* **2019**, *9*, 57

**Figure 7.** 3D envelopes and tips length of a single growing grain at timet=7s after nucleation, where the grain is assumed to nucleate and grow with an initial undercooling ΔT0 = 2 K, the crystallographic orientation (20◦, 40◦, 60◦) and the growth constant A = 0.0001 m/(s·K2). (**a**) 3D envelopes in a uniform temperature field; (**b**) tip length in a uniform temperature field; (**c**) 3D envelopes in a Bridgman-like situation; and (**d**) tip length in a Bridgman-like situation.

As shown in Figure 7a,c the 3D views of analytical and numerical predictions, for a given crystallographic orientation, are approximately regular octahedrons in a uniform temperature field. In Bridgman-like situations, due to the effects of deep undercooling, the tips length (*L*4, *L*6) of the grain closest to the negative thermal gradient direction are the longer distances. Altogether, the 3D grain envelopes (a2, c2) simulated by the developed model and algorithm in both cases are visually reproduced and coincident with the shape (a1, c1) predicted by the analytical model. Furthermore, from the Figure 7b,d, it can be noted that there are only a few mesh size gaps among the tips length of the grain in both predictions, which also approximately present the same distribution.

### *3.3. Model Application to Nickel-Based Superalloys*

The present model is used to simulate the formation of equiaxed and columnar dendritic grains during nickel-based superalloys solidification process. The elemental composition (mass fraction, wt.%) of the material is Ni–7.82Cr–5.34Co–2.25Mo–4.88W–6.02Al–1.94Ti–3.49Ta. For multicomponent alloys, a large number of simplifications and calculations are required due to the complex multicomponent

coupling solidification process. Here, a multicomponent pseudo-binary alloy method [37] is used. The material properties and model parameters used in these simulations are listed in Table 1.


**Table 1.** Material properties and model parameters.

### 3.3.1. Equiaxed Dendritic Grains

The calculation domain is uniformly divided into 200 × 200 × 200 meshes with mesh size of 0.0001 m. It is assumed that the temperature field is homogeneous and cooled down from liquidus temperature at different cooling rates leading to nucleation and growth process. The cooling rates are set to 1 K·s−1, 10 K·s−1, and 100 K·s−1, respectively. Here, the nucleation process only takes into account the nucleation inside the melt. It must be stressed that the nucleation parameters are affected by many factors, such as random nucleation, solidification conditions, etc. It cannot be uniquely determined or calculated in advance. For the ease of numerical calculations, the nucleation parameter is assumed to be known in advance. In the bulk of the melt, the total density of grains is 1 × <sup>10</sup><sup>10</sup> <sup>m</sup><sup>−</sup>3, the median undercooling is 3 K, and the standard deviation of undercooling is 0.5 K. Figure 8 presents the size distribution of equiaxed dendritic grains at different cooling rates.

**Figure 8.** *Cont*.

**Figure 8.** Size distribution of equiaxed dendritic grains at different cooling rates. (**a**) Simulation results of three-dimensional equiaxed dendritic grains; and (**b**) size distribution of equiaxed dendritic grains in a two-dimensional section.

(**b**)

Figure 8 shows that, with the increase of cooling rate, the grain size of equiaxed grains decreases obviously, and the grain density increases from 0.13 mm−<sup>2</sup> to 4.08 mm−2. From the data analysis of the two-dimensional section, it can be seen that the number of grains approximately presents a normal distribution curve at d/<d> from 0 to 2.0. With higher cooling rates the maximum of the grain number as function of d/<d> grows, while its dispersion in grain-size lowers. The highest grain number-density is always given close to the median grain-size. It is indicated that the increase of cooling rate can significantly promote the grain refinement. These results are consistent with the experimental results of solidification microstructures [38], which preliminarily validate the effectiveness of the proposed model and algorithm.

### 3.3.2. Columnar Dendritic Grains

Cooling rate and temperature gradient are two important parameters for directional solidification microstructure. Since the thermal diffusion rate is more than several orders of magnitude larger than the solute diffusion rate, it can be assumed that the temperature is uniform at the scale of a grain. At a fixed temperature gradient (*G*) and cooling rate (*C*) the temperature isotherms of the whole calculation domain will change at the same drawing speed (*R* = *C***/***G*). The calculation domain is uniformly divided into 60 × 180 × 600 cubic cells with size of 0.0001 m. Two types of nucleation parameters are considered. At the bottom surface, the total density of grains is 5.0 × <sup>10</sup><sup>5</sup> <sup>m</sup><sup>−</sup>3, the mean undercooling is 0 K, and the standard deviation is 0 K. In the bulk of the melt, the total density of grains is 1 × <sup>10</sup><sup>10</sup> <sup>m</sup><sup>−</sup>3, the mean undercooling is 3 K, and the standard deviation is 0.5 K.

Equiaxed-to-Columnar Transformation (ECT) for Different Cooling Rates

The temperature gradient is set to 1000 K·m<sup>−</sup>1. And the cooling rates are set to 1 K·s<sup>−</sup>1, 10 K·s<sup>−</sup>1, and 100 K·s<sup>−</sup>1, respectively. Figure <sup>9</sup> presents the grain size distributions for different cooling rates.

**Figure 9.** *Cont*.

**Figure 9.** Grain size distributions for different cooling rates. (**a**) *<sup>C</sup>* =1K·s−1; (**b**) *<sup>C</sup>* = 10 K·s−1; (**c**) *<sup>C</sup>* = 100 K·s<sup>−</sup>1; and (**d**) grain number along the *<sup>Z</sup>* axis.

Figure 9 shows, that for different cooling rates, the developed solidification microstructures of the calculated domain gradually changes (with respect to the size distribution of grains). At the cooling rate of 1 K·s−1, some thick columnar dendritic grains are shown in the calculation domain, and the number of grains decreases gradually along the *<sup>Z</sup>* axis. At the cooling rate of 10 K·s−1, the grain structure in the calculation domain is composed of both columnar and fine equiaxed dendritic grains, and the number of grains increases first and then decreases. At the cooling rate of 100 K·s−1, almost all grains in the calculation domain are small equiaxed dendritic grains, and the number of grains increases suddenly and decreases slowly with increasing Z. From the above results, it can be inferred that at low cooling rate the grain growth rate of directional solidification is faster than the grain nucleation rate in the supercooled region of melt. In addition, due to the boundary limitation of the calculation domain, some grains deviating from the direction of temperature gradient will be eliminated at the boundary, resulting in the gradual decrease of the number of grains along the Z axis.

Equiaxed-to-Columnar Transformation (ECT) for Different Temperature Gradients

The temperature gradient is set to 500 K·m−1, 1000 K·m−1, and 1500 K·m−1, respectively. The cooling rates are set to 10 K·s−1. Figure <sup>10</sup> presents the grain size distributions for different temperature gradients.

**Figure 10.** Grain size distributions for different temperature gradients. (**a**) *<sup>G</sup>* = 500 K·m−1; (**b**) *<sup>G</sup>* = 1000 K·m<sup>−</sup>1; (**c**) *<sup>G</sup>* = 1500 K·m<sup>−</sup>1; and (**d**) grain number along the *<sup>Z</sup>* axis.

Figure 10 shows for the calculated domain, that with increasing temperature gradient, the solidification microstructure develops to the opposite of the results presented in Figure 8. In other words, at the temperature gradient of 500 K·m<sup>−</sup>1, almost all grains in the calculation domain are small equiaxed dendritic grains, and the number of grains increases suddenly and decreases slowly. At the temperature gradient of 1500 K·m<sup>−</sup>1, some thick columnar dendritic grains are shown in the calculation domain, and the number of grains decreases gradually along the Z axis. It can also be deduced that a rapid growth of columnar dendritic grains is promoted more by an increasing temperature gradient than by an increasing nucleation rate.

### 3.3.3. Columnar-to-Equiaxed Transformation (CET)

The ECT and CET are common phenomena of the grain evolution in solidification process. Under different solidification conditions, the transformations rule of the two influences the grain size distribution in the casting. During directional solidification, the ECT transformation is beneficial to the

growth of columnar dendritic grains, and it is also affected by the nucleation parameters and cooling rate. Since the cooling rate and the temperature gradient are difficult to control, the casting of the metallic alloy may exhibit a non-uniform distribution of the temperature field. Especially, when the temperature gradient gradually decreases to a negative temperature gradient, this will lead to the CET transformation. To characterize the positive-to-negative transition in the temperature gradient, the temperature gradient is set to 1000 K·m−<sup>1</sup> from mesh point 0 to 300 in the direction of the *<sup>Z</sup>* axis, and the temperature gradient is set to −300 K·m−<sup>1</sup> from mesh point 300 to 600 in the direction of the *<sup>Z</sup>* axis. The cooling rate is set to 10 K·s−1. The other grain nucleation and growth parameters set here are the same as those above Section 3.3.2. Figure 11 presents the simulation results of the CET transformation.

**Figure 11.** CET transformation. (**a**) t = 0.5 s; (**b**) t = 2.0 s; (**c**) t = 2.7 s; (**d**) t = 3.3 s; (**e**) grain number (0 ≤ Z ≤ 300); and (**f**) grain number (300 ≤ Z ≤ 600).

As shown in Figure 11, in the calculation domain, the grains first nucleate and grow at the bottom surface. From mesh point 0 to 300 in the direction of the *Z* axis, under the positive temperature gradient, grains begin to compete with each other along the *Z* axis, which is similar to the ECT transformation process mentioned above. From mesh point 300 to 600 in the direction of the *Z* axis, the grains nucleate and grow at the top surface under the negative temperature gradient. Compared with the ECT transformation process, the grains grow in the form of a large number of equiaxed

crystals. Finally, the two parts of the grain stop growing at intermediate position, thus forming the CET transformation process. The above results show a good agreement with the results of related literature [23,39,40], which further validates the applicability of the proposed model and algorithm.

### **4. Conclusions**

A cellular automata model of 3D mesoscopic grain evolution is developed. In this model, the 3D grain envelope trajectory is coupled with the solid fraction in a cell. The solid fraction in grain growth is calculated quantitatively at the scale of a CA grid. The detailed model and algorithm validations are performed by comparing the numerical results with the analytical model predictions. The simulated 3D grain envelopes are found to be in agreement with the analytical predictions. The developed model was then applied to simulate the formation of equiaxed and columnar dendritic grains during the solidification process of nickel-based superalloys. The results of this simulation could well explain the evolution law of grain growth in solidification microstructures, including equiaxed dendritic grain growth, columnar dendritic grain growth, and microstructure transformation.

Although the mesoscopic grain growth in a simple temperature field is studied in this present work, the proposed grain envelope trajectory coupled with the cell solid fraction method can be easily extended to complex multiscale and multiphysical field calculations. Such a 3D mesoscopic grain model coupled to a macroscopic temperature field and solute field is in progress and can be used to simulate the solidification process of actual castings.

**Author Contributions:** Z.G. and J.Z. wrote the paper; Z.G. and Y.Y. did the simulations; Z.G. and X.S. analyzed the data; Z.G. and X.J. revised the paper.

**Funding:** This work was supported by the National Natural Science Foundation of China (No. 51605174, No. 51775205) and AECC Beijing Institute of Aeronautical Materials.

**Conflicts of Interest:** The authors declare no conflicts 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* **Analysis of Forming Parameters Involved in Plastic Deformation of 441 Ferritic Stainless Steel Tubes**

**Orlando Di Pietro 1, Giuseppe Napoli 1, Matteo Gaggiotti 1, Roberto Marini <sup>2</sup> and Andrea Di Schino 1,\***


Received: 19 June 2020; Accepted: 23 July 2020; Published: 28 July 2020

**Abstract:** A welded stainless steel tube is a component used in several industrial applications. Its manufacturing process needs to follow specific requirements based on reference standards. This calls for a predictive analysis able to face some critical issues affecting the forming process. In this paper, a model was adopted taking into account the tube geometrical parameters that was able to describe the deformation process and define the best industrial practices. In this paper, the effect of different process parameters and geometric constraints on ferritic stainless steel pipe deformation is studied by finite element method (FEM) simulations. The model sensitivity to the input parameters is reported in terms of stress and tube thinning. The feasibility of the simulated process is assessed through the comparison of Forming Limit Diagrams. The comparison between the calculated and experimental results proved this approach to be a useful tool in order to predict and properly design industrial deformation processes.

**Keywords:** stainless steels; plastic deformation; mechanical properties

### **1. Introduction**

Stainless steels are nowadays used in many applications targeting corrosion resistance coupled to mechanical resistance [1,2]. In particular, they are employed in automotive [3], construction and building [4–6], energy [7–11], aeronautics [12], food [13–15] and three-dimensional (3D) printing [16] applications. Stainless steels find application in automotive field, thanks to their ability to be worked into complex geometries [17–19]. Compared to the most common techniques, the metal forming is used to manufacture complex geometries in efficient way. In order to demonstrate compliance with the most relevant standards, many quality tests are commonly carried out in the industrial sector: this implies a strong increase in terms of cost and time-consumption.

The modelling and simulation of manufacturing processes is a well-established prerequisite for cost-effective manufacturing and quality-controlled manufacturing plant operation. Therefore, the management and control of manufacturing, which have the potential to be game changers, are a key issue. A game-changing technology is one which either directly provides a new and advantageous method of manufacture or else provides an important missing link, thereby enabling a method of manufacturing which would not otherwise be viable or so significantly advantageous. This should be based on the development and application of models able to support manufacturing process stability and product quality.

In the case of ferritic steels, for example, the plastic processing of welded pipes is characterized by a poor behavior homogeneity [20]; because of the nature of this material, characterized by an intrinsic inhomogeneity, a certain percentage of tests will not be reliable. The quality control, performed after the manufacturing process, is generally carried out by means of tensile tests according to specifications. In many cases, the measured tensile properties are found to be not sufficient to guarantee compliance with the required standards. The numerical simulation of sheet metal forming processes has become an indispensable tool for the design of components and their forming processes [21,22]. This role was attained due to the strong impact in reducing the time to market and the cost of developing new components in industries ranging from automotive to packing, as well as enabling an improved understanding of the deformation mechanisms and their interaction with process parameters. Despite being a consolidated tool, its potential for application continues to be discovered with the continuous need to simulate more complex processes, including the integration of the various processes involved in the production of a sheet metal component and the analysis of the in-service behavior. The request for more robust and sustainable processes has also changed its deterministic character into stochastic to be able to consider the scatter in mechanical properties induced by previous manufacturing processes.

In this framework, many tube bending approaches have been developed in response to the different demands of tube specifications, shapes, materials and forming tolerances [23–25]. Many mathematical models able to describe the steel plastic deformation behavior have been developed, considering the type of material and application field [26]. All of these approaches require the analysis of the mechanical properties, considering the steel in its macroscopic and microscopic structure. Another relevant characteristic concerns the anisotropic behavior caused by the plastic deformation due to the pipe manufacturing. He et al. [27] reported about the occurring mechanisms, accurate prediction, and efficient controlling of multiple defects/instabilities during bending forming. Wu et al. [28] investigated the effects of the temperature, bending velocity, and grain size on the wall thickness change, cross-section distortion, and spring-back of a tube in rotary draw bending. Liu et al. [29,30] analyzed the effect of dies and process parameters on the wall thickness distribution and cross-section distortion of a thin-walled alloy tube. Tang [31] deduced several bending-related equations to predict the stress distributions of bent tubes, wall thickness variation, bending moment, and flattening based on the plastic-deformation theory. Al-Qureshi et al. [32] derived the approximate equations of tube bending to predict spring-back and residual stress quantitatively with assumptions of the ideal elastic-plastic material, plane strain condition, absence of defects, and "Bauschinger effect". E et al. [33] derived an analytical relation to predict stress distributions, neutral layer shift, wall thickness change, and cross-section distortion based on in-plane strain assumption in tube bending with the exponential hardening law. Li et al. [34], based on the plastic-deformation theory, established the analytical prediction model for spring-back, taking into account the tube specifications and material properties. Jeong et al. [35] proposed equations able to calculate the bending moment and spring-back of bent tubes, considering the work hardening in the plastic-deformation region. The cross-section distortion formula was deduced based on the virtual force principle by Liu et al. [36]. It is noted that, though many factors cannot be considered, such as contact conditions and unequal stress/strain distributions, the analytical models are still very useful to estimate and predict the forming quality for tube bending.

Considering the complexity of the tube bending forming process, the finite element numerical simulation method has been widely used to explore the bending deformation under various forming conditions [37–39]. The main goal of simulations is to predict the behavior of different tubes' geometries during the bending process or the cold metal forming of steel sheets. In fact, a proper procedure of tube bending and a correct simulation of pipe yielding, after bending, turns out to be critical. Several approaches devoted to ferritic and austenitic stainless steel grades are based on the Von Mises and Johnson–Cook criteria [40,41]. Such criteria describe the elastic-plastic behavior of isotropic materials, defining the stress induced as a function of the deformation, strain rate, and temperature. A step forward is based on Hill's criterion, which introduces different equations for orthotropic and anisotropic materials [42,43].

In this paper, the deformation process of a ferritic stainless steel tube is simulated by means of a commercial software package adopting Hill's criterion. The results of the simulations are compared to those coming from experimental tests, with the aim to develop a tool able to predict the bending behavior of stainless steel tubes easily applicable in industrial applications. As a matter of fact, even if many researchers are active in this field and fundamental papers have been published, stainless steel tube bending is commonly performed in industries just based on an empirical approach (trial and error). The main novelty of the paper is to present an approach able to cover the gap between fundamental theoretical research and industrial application, putting in evidence the effect of bending process parameters. A proper forming process can therefore be designed by the comprehension of the effect of different pipe geometric parameters on the plastic deformation operation, thus allowing an optimum process fitting to the component under examination.

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

### *2.1. Materials*

The steel chemical composition of the selected material (EN 1.4509 441) and the geometrical parameters of the tubes considered in this paper are reported in Tables 1 and 2, respectively.


**Table 1.** Chemical analysis of 441 (main elements, mass. %).

**Table 2.** Materials selected for the simulations with their geometric characteristics.


### *2.2. Methods*

The steel, for each combination of diameter and thickness, has been tested according to the UNI EN 69892−1 and UNI EN 6892−2 standard. A Formability Limit Diagram (FLD) determination has been carried out to describe the sample deformation paths. This type of diagram contains the Formability Limit Curve (FLC), which shows the maximum capacity of a material to be deformed, calculated by carrying out repeated Nakazima tests and measuring the deformation along the two perpendicular directions [44]. Based on this test, the strains are measured with the conventional grid method, plotting a pattern of circles on the specimen. Such circles are deformed into ellipses during the deformation process. On these ellipses, strains are measured on minor and major sizes, thus identifying on the FLC diagram the deformation state points of the material. The materials properties needed as input parameters for the subsequent modelling are reported in Table 3. They include the steel density, Young modulus, Poisson ratio, Lankford value, and strain hardening coefficient.

**Table 3.** Steel properties [2].


### *2.3. Model*

The FEM model used in such an investigation is based on the rigid-plastic variational principle following:

$$
\pi = \int\_{V} \overline{\sigma^{\varepsilon}} dV - \int\_{S\_F} F\_i u\_i dS\_\prime \tag{1}
$$

which defines the allowed deformation rates. In Equation (1), the incompressibility condition can be driven by adding a penalty constant α, as reported in Equation (2):

$$
\pi = \int\_{V} \overline{\sigma \dot{\varepsilon}} dV + \int\_{V} \frac{\alpha}{2} \dot{\varepsilon}\_{V} \overline{\ }^{2} dV - \int\_{S\_{\overline{\mathbb{F}}}} F\_{i} u\_{i} dS\_{\ast} \tag{2}
$$

thus:

$$
\delta\pi = \int\_{V} \overline{\sigma} \delta\overline{\overline{\varepsilon}} dV + \alpha \int\_{V} \dot{\varepsilon}\_{V} \delta\dot{\varepsilon}\_{V} dV - \int\_{S\_{\overline{\varepsilon}}} F\_{l} \delta u\_{l} dS\_{\star} \tag{3}
$$

where σ, . ε, ε*V*, *Fi*, *ui* are the average stress, deformation coefficient, volume deformation coefficient, surface stress, and speed of the pipe, respectively. The 8-nodes element is adopted to simulate the pipe bending process. The position (*x*,*y*,*z*) and speed vectors (*ux*, *uy*, *uz*) of node-i are therefore defined as follows:

$$\begin{aligned} \mathbf{x} &= \sum\_{i=1}^{8} N\_i(\xi, \eta, \zeta) \mathbf{x}\_i \\ y &= \sum\_{i=1}^{8} N\_i(\xi, \eta, \zeta) y\_{i\prime} \\ z &= \sum\_{i=1}^{8} N\_i(\xi, \eta, \zeta) z\_{i\prime} \\ u\_x &= \sum\_{i=1}^{8} N\_i(\xi, \eta, \zeta) u\_{x\prime} \\ u\_y &= \sum\_{i=1}^{8} N\_i(\xi, \eta, \zeta) u\_{y\prime} \\ u\_z &= \sum\_{i=1}^{8} N\_i(\xi, \eta, \zeta) u\_{z\prime} \end{aligned} \tag{4}$$

where ξ, η, ζ are the co-ordinates in the cell space.

The model takes into account of the friction subjected to the pipe during bending according to the Tresca model [20]:

$$f = -\frac{2}{\pi}mk \tan^{-1} \left(\frac{u\_s}{u\_0}\right)t,\tag{5}$$

where *m* is the friction factor, *k* is the effectve shear stress, *us* is the speed of the pipe, and *t* is the unit vector in the *us* direction. Therefore, a FEM model has been implemented, considering the tube as discretized with a 4 × 4 mm × mm mesh. The pressure die, wiper die, and bending die tools are discretized, according to the theory of plastic deformation, such as rigid body.

The input related to the shape of tools are described in Figure 1.

**Figure 1.** Tools for the bending simulation.

The input tools sizes are reported in Table 4.

**Table 4.** Adopted bending tool sizes.


The other related tools' geometric parameters are calculated by simulation. Internal spherical elements supporting the bending process (so called balls) are present even if they are not reported in Figure 1. In this paper, four balls with 3.3 mm decreasing diameter sizes are considered. The starting diameter of the first ball is calculated by decreasing the tube diameter by the same measure. It is worth mentioning that an additional support machinery element (called a booster) is not taken into account in the calculations.

The numerical calculation was performed using the Altair HyperWorks™ (2017 version, Troy, MI, USA, commercial software). Such software is able to take into account the material properties reported in Table 3 together with the true stress–true strain steel curves. As far as concerns the tensile material properties, five true stress–true strain curves have been considered, aiming to obtain an average curve. The above materials properties, in particular, are collected in a file readable by the commercial package solver (namely, RADIOSS™). The solver is able to take into account the strain rate by reading the tensile curves performed at different deformation rate conditions. This software allows the adoption of the Hill 48 yield function. Such a function is known to be ideal for small-sized tubular geometries [45] as a constitutive equation for stainless steel behavior, taking into account the following parameters in order to simulate the bending process:


The analysis of the simulation outputs is carried out through mappings of the values calculated by the solver (e.g., internal stress, thinning, and deformation). An example of the obtained maps is reported in Figure 2.

**Figure 2.** Stress mapping on bent tube.

In order to analyze the parameters' effects on the final process and its feasibility, the maximum stress values obtained on the mappings will be considered in order to consider the critical points of the geometry.

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

The effects of the geometric and operational inputs parameters described for the tubes' plastic deformation are reported below.

### *3.1. E*ff*ect of Tube Diameter*

In this section, the effect of the tube's diameter on the bending process is analyzed, keeping the *R*/*D* ratio as a fixed value. *R*/*D* = 1 is considered being such a value significative in order to simulate what is commonly performed in the industrial bending process. The tube's stress behavior as a function of the diameters is reported in Figure 3.

**Figure 3.** Mean maximum stress behavior as a function of diameter size for 441 steel with a 1.5 mm thickness.

The results allow us to evaluate an effect of about 5% of the tube's diameter on the internal stress. The same effect is also found for the tube's thinning (Figure 4). Such results well fit with what is expected from common industrial experience. In addition, they allow us to estimate the stress and thinning loss as a function of the tube's diameter, thus allowing a process optimization for different tube geometries.

**Figure 4.** Maximum thinning as a function of the diameter size for 441 steel with a 1.5 mm thickness.

In order to evaluate the deformation capacity and support the above calculation outcomes, the results from FLD diagrams were considered. Figure 5 shows the formability limit (%) as reached during deformation for 441 stainless steel (diameter sizes ranging from 40 to 60 mm). The dashed red line represents the sample breakage. The FLD diagrams confirms that the deformation of the various geometric elements is affected by the diameter size. The results show that an increase in the diameter size will allow a decrease in the breakage risk. Such a result allows us to design the forming process for different pipes sizes, avoiding pipe breakage.

**Figure 5.** Formability limit (%) as a function of the tube's diameter for 441 with a 1.5 mm thickness.

### *3.2. Tubes Thickness E*ff*ect*

The stress and the tube's thinning dependence on the thickness are reported in Figures 6 and 7, respectively. The results show that a not-significative stress distribution was found. In particular, the total variation is lower than 2%. On the other hand, the thinning (%) increases as the tube's thickness increases (Figure 7). Although there have been not relevant variations in terms of stress and thinning, the trend reported in Figure 7 is significant. As a matter of fact, while in the case reported in Section 3.1 the thickness reduction (which is kept constant) has the possibility of being distributed over ever larger diameters, in this case the diameter was fixed as a constant parameter; following that, an increase in thickness led to an increase in the tube's thinning.

**Figure 6.** Mean maximum stress behavior as a function of thickness for 441 steel with a 50 mm diameter.

**Figure 7.** Maximum thinning as a function of thickness size for 441 steel with a 50 mm diameter.

An FLD plot as a function of the tube's thickness is reported in Figure 8. The results show that the success of the manufacturing process depends on the tube's thickness. Figure 8 clearly shows the strong effect exercised by the initial thickness on the bending process success.

**Figure 8.** Formability limit (%) as a function of the thickness for 441 with a 50 mm diameter.

### *3.3. Speed and Bend Angle E*ff*ects*

The effect of the speed and bending angle are reported in the following Figure 9.

As far as concerns the speed, the influence of its variation for every bending angle (in the range of 30◦ and 90◦ for a common manufacturing process) has been considered. The formability limit (%) for the combination of angle and thickness is shown in Figure 9. In particular, the geometric parameters were fixed, together with the relationship between the diameter and the thickness.

**Figure 9.** A 2D plot of the formability limit (%) for different speed and angle combinations.

In order to analyze the data, the lines shown in Figure 9 have been interpolated and reported in Figure 10 to better evaluate the influence of the speed. For this reason, the percentage variations between the percentages of formability limits obtained at minimum and maximum feed speeds for each angle, calculated according to Equation (6), are reported in Figure 11.

$$
\Delta FLD = FLD\_{\upsilon\_{\text{max}}} - FLD\_{\upsilon\_{\text{min}}} \,. \tag{6}
$$

**Figure 10.** Formability Limit Diagram (FLD) dependence on speed for different bending angles (linear interpolation).

**Figure 11.** Goodness of the simulation output beyond the breaking of the worked piece (red dotted line).

Δ*FLD* is closely related to the slope of the interpolated line. Figure 11 shows that in the angle region ranging between 30◦ and 90◦ degrees (the most interesting for common industrial processes), the Δ*FLD* varies almost linearly with respect to the bending angle. If higher angle values are considered the results show that such a parameter reaches a maximum at about 120 degrees of bending, then decreases. The motivations leading to this particular behavior need to be deeply analyzed, but currently we can hypothesize that this is due to the stress concentration being more localized in the first 90◦ of bending.

In order to obtain more consistent results, the analysis has been repeated using conditions better reproducing the industrial process. Higher curvature radius values (and the *R*/*D* ratio, consequently) were considered. After that, the data were newly interpolated (Figure 12) and the FLD deltas were calculated for the new data set (Figure 13).

**Figure 12.** Linear interpolation of the formability limit (%) for each combination of speed and angle.

**Figure 13.** FLD dependence on the bending angle.

Figure 13 shows that the better choice of the *R*/*D* ratio leads to an improvement in the sample formability. ΔFLD increases in the bending angle range 30◦–90◦; it then tends to be stable, keeping away from breaking conditions. Such results suggests that a focus should be put on the *R*/*D* ratio on the process. As a matter of fact, this is one of the key issues in the setup of an industrial bending process.

### *3.4. R*/*D E*ff*ect*

The *R*/*D* ratio is commonly adopted as a feasibility index for the bending industrial process. In common industrial practice, this value ranges between 1.0 and 1.5. In fact, *R*/*D* < 1 values increase the breakage risk; on the other hand, an *R*/*D* > 1.5 is not commonly adopted in the automotive field. The results reported in Figure 14 show a negligible effect of *R*/*D* on stresses. On the other hand, a marked *R*/*D* ratio effect is found on the tube's thinning (Figure 15). The results show that the tube's thinning decreases as the *R*/*D* increases. Additionally, in this case the results well fit with what is expected from common industrial experience. In addition, they allow us to estimate the stress and thinning loss as a function of the tubes' *R*/*D* parameter.

**Figure 14.** Maximum equivalent stress dependence on the *R*/*D* ratio.

**Figure 15.** Maximum thinning dependence on the *R*/*D* ratio.

### *3.5. Experimental Validation*

Six samples of 441 steel tubes (50 mm diameter and 1.2 mm thickness size) were considered for validation. All of the considered samples have been chosen with sizes corresponding to the above reported simulations. One of the six samples (for each group) was subjected to tensile tests in order to measure the related stress–strain curve and reduce the uncertainty due to the use of a mean curve. Experimental tests have been carried out according to operational parameters (e.g., rotational velocity and bending radius), as reported in Table 5. The chosen values are those commonly adopted by industries operating in pipe forming for the typical finishing process. Both the thicknesses reached during the bend along the backbone at specific angles, as shown in Figure 16, were measured.

**Table 5.** Testing conditions.


**Figure 16.** Thickness measuring grid on the backbone.

The thickness values were measured for each marked angle, and the average value was considered. The obtained values are shown in Table 6.


**Table 6.** Thickness for the 441 steel samples as measured at different considered angles (mm).

The measured thickness variation with the angle and the percentage discrepancy between the measured and calculated thickness as a function of the angle are reported in Table 7. In Table 7, ΔThickness is defined as the difference between the simulation thickness and the sample mean thickness, as conventionally adopted by companies operating in the tube forming sector.

**Measurement Angle Simulation Thickness [mm] Sample Mean Thickness [mm] <sup>Δ</sup>Thickness [mm] Percentage Variation [%]** 0◦ 1.168 1.193 −0.025 −2.10 22.5◦ 1.014 1.049 −0.035 −3.35 45◦ 0.876 1.035 −0.159 −15.34 67.5◦ 0.892 1.067 −0.175 −16.36 90◦ 1.171 1.170 0.001 0.06

**Table 7.** Comparison between the calculated and measured thickness at different angles.

The deviation of the model from the experimental results (maximum underestimation of 16.36%) corresponds to 67.5 degrees. Such a deviation is probably related to the presence in the experimental tests of an additional support machinery element (the so-called booster), which was not considered in the simulation model. The booster effect in the industrial deformation process is that of pushing the tube during bending in order to avoid the deformations or failure caused by the friction between the element and the machine or by stress concentration. Its action also affects the distribution of the thinning caused by the deformation. As a matter of fact, the tube being pushed by the booster will have more evenly distributed stresses. As a consequence, the deformations and the thinning will take place on a wider area and will not lead to the failure of the piece. Anyway, the accuracy between the modelling and experimental tests is considered good and quite promising, since the discrepancy between the calculated and experimental values is lower than 15% for an angle position lower than 45◦. Moreover, the fact that the model appears to underestimate the thickness values allows us to consider its results conservative with respect to the real behavior.

### **4. Conclusions**

This paper analyzes the influence of geometric and operational parameters on the bending process of 441 ferritic stainless steel pipes. Experimental investigations combined with simulations highlighted the effect of each parameter, both operational and geometric, on the final results. The main novelty of the paper is to present an approach able to bridge the gap between fundamental theoretical research and industrial application, putting in evidence the effect of bending process parameters.

Based on the reported results, the following conclusions can be drawn:

• Pipe diameter effect on forming: If the same parameters are considered and the R/D ratio between the radius and diameter is kept constant, an increase in the diameter size (in the typical automotive range) will result in a −9% variation in terms of internal stresses, evaluated according to the von Mises criterion. The thinning of the tube will decrease by −4%, and the feasibility characteristics of the process will improve; in fact, the percentage reaching the formability limit will drop by −20%. The diameter size choice appears therefore to be a key issue in terms of process feasibility. In this sense, the reported approach is a useful tool aimed to properly design, also in quantitative terms, the forming process for different pipe sizes.


The comparison between the calculation and experimental results proved the reported approach to be a useful tool in order to predict and properly design industrial deformation processes. The experimental validation showed a deviation of the model from the experimental case, with a maximum underestimation of 16.36%. Such discrepancy is probably related to the presence in the experimental tests of a booster, which was not included in the simulation model.

**Author Contributions:** Conceptualization, O.D.P. and R.M.; methodology, O.D.P. and R.M.; formal analysis, O.D.P., G.N., and M.G.; paper preparation, O.D.P. and A.D.S.; supervision, R.M. and A.D.S. All authors have read and agreed to the published version of the manuscript.

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

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

### **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/).
