Mesoscale modeling and simulation of microstructure evolution during dynamic recrystallization of a Ni-based superalloy

Microstructural evolution and plastic flow characteristics of a Ni-based superalloy were investigated using a simulative model that couples the basic metallurgical principle of dynamic recrystallization (DRX) with the two-dimensional (2D) cellular automaton (CA). Variation of dislocation density with local strain of deformation is considered for accurate determination of the microstructural evolution during DRX. The grain topography, the grain size and the recrystallized fraction can be well predicted by using the developed CA model, which enables to the establishment of the relationship between the flow stress, dislocation density, recrystallized fraction volume, recrystallized grain size and the thermomechanical parameters.


Introduction
The mechanical properties of a metallic component such as high-temperature alloy are largely depended on their chemical composition and microstructures, resulting from processing and heat treatment history in industrial production. Accurate prediction and precise control of the microstructural evolution are of great importance for designers to achieve excellent mechanical properties for metals and alloys through hot working. Up to date, DRX has been recognized as a powerful mechanism for grain refinement under hot plastic deformation. In this context, in-depth understanding of microstructural evolution during DRX can help develop more accurate prediction methods to simulate microstructural characterization and mechanical properties in the hot deformation of metals and alloys [1][2][3].
With the rapid development of computer technology and a thorough understanding of microstructural evolution thermodynamics and kinetics, mesoscale methods, i.e., the CA method, Monte Carlo (MC) method, vertex models, phase-field approach and level-set modeling, have been successfully used to simulate microstructural evolution in hot working. Mesoscale models can not only predict the general average microstructure properties (the kinetics of recrystallization and grain size), but also able to visually provide the virtual microstructural evolution during hot deformation. Compared with the other methods, CA method can be used as an alternative approach to describe the spatial and temporal microstructure evolution on a mesoscale. There are some shortcomings of CA method. For example, CA has poor representation of grain boundary curvature as compared to the phase-field model and levelset model. The first attempt to model of microstructural evolution during DRX by CA method was made by Goetz and Seetharaman [4]. Up to now, there are many cases of successful application of CA in modeling microstructural evolution. In these studies, Ding et al. [5], Raabe [6], Kugler and Turk [7], Xiao and Zheng [8] and Jin and Cui [9] contributed greatly to the development of the CA method for the simulation of DRX from the aspects of optimization algorithm, grain morphology change and parameter identification. Recently, based on the work of Ding et al. [5] and Kugler and Turk [7], Yazdipour et al. [10] simulated DRX during thermomechanical deformation of 304 austenitic stainless steel by the CA method. In this case, initial grain size, initial grain orientation and dislocation density are used as input data to the CA model; the flow curve, dislocation density, finial grain size and orientation and recrystallized fraction are the output data. Hallberg et al. [11] employed the CA model with probabilistic cell switches in simulating DRX of pure Cu. Ding et al. [12] investigated the hot compression behavior of the Mg-9Al-1Zn alloy using EBSP analysis and the CA method. Chen et al. [13,14] predicted microstructural evolution during DRX of an ultra-super-critical rotor steel. In addition, the DRX behavior of TRIP steel [15], C-Mn micro-alloy steel [16], Ti-alloy [17] and Ni-based alloy [18] was studied by using improved CA models. However, in order to perfect the classical CA model for DRX, the following two aspects still need further investigation. Firstly, it is highly necessary to develop reasonable and practicable method to calculate CA simulation constants such as nucleation parameter. Secondly, it is desirable to introduce the mechanisms of grain morphology change on the nucleation and grain growth during DRX in either twodimensional or three-dimensional cases.
Due to its excellent mechanical characteristics at elevated temperatures up to 650°C, good fatigue and corrosion resistance, Ni-based superalloy named IN 718 has been widely used in aerospace industry [19]. Usually, IN718 superalloy components are processed by hot forging. The mechanical properties of forgings are extremely sensitive to the dynamic and interconnected variables, such as strain, strain rate and temperature [20]. Therefore, in recent years, the flow behaviors and microstructural evolution of IN718 have been studied experimentally by many researchers [21][22][23][24][25][26]. However, there have been few investigations on modeling and simulation of microstructure evolution for IN718 superalloy during hot working. In order to deeply understand DRX mechanism and to predict and control the microstructure evolution during DRX for IN718, there is a need to develop new model to simulate the grain refinement during hot working from mesoscale aspect.
The main purpose of this article is to develop a 2D CA model to predict the microstructure evolution for IN718 superalloy during hot working. The CA model has the following advantages. Firstly, reverse algorithm was adopted for parameter identification. Secondly, grain morphology deformation technique was used in the CA model, in which a cellular coordinate system (CCS) and a material coordinate system (MCS) were established separately. The CCS remains unchangeable in the whole simulation. The MCS and the corresponding grain boundary shape change with deformation in the developed DRX module. By using the developed model, a wide range of deformation temperature and strain rate over DRX were examined to investigate the quantitative and topographic prediction of the microstructural evolution during hot working for IN718 superalloy. As a result, the appropriate hot working parameters for IN718 were recommended based on the CA simulations.

Experiment materials and procedures
The chemical composition (wt%) of IN718 superalloy used in this study is as follows: 18.07Cr-17.74Fe-5.41Nb-2.86Mo-1.02Ti-0.9Si-0.5Al-0.025C-(bal.)Ni. The diameter and length of the specimens are 8 mm and 12 mm, respectively. According to ISO 6892-2, suitable test procedures were designed. All the specimens were solution treated at 1100°C for 0.75 h and then cooled to room temperature in the air. The isothermal compression tests were carried out at the temperatures of 940, 980, 1020 and 1060°C and at the strain rates of 0.001, 0.01, 0.1 and 1 s -1 . The tests were performed on a Gleeble-1500 thermomechanical simulator. Before loading, the specimens were heated to the test temperature with the heating rate of 10°C/min. The maximum true strain is 1.0. In order to study the effects of forming process on the microstructure, the specimens were immediately quenched by water after hot compression. The deformed specimens were sliced along the compression axis section for microstructural analysis. The tested specimens were prepared for optical microscopy (OM) using standard metallographic techniques of mounting, coarse and fine grinding, final polishing and etching. The specimens were immersion etched with HCL (20 ml) ? 40 % HF (50 ml) ? H 2 O 2 (2 ml) to reveal the microstructures. Etching times varied with test conditions from 2 to 3 min at room temperature. As shown in Fig. 1, the optical microstructures in the section plane are examined. The average grain sizes for the nine positions at axial section are measured using the method described in ASTM standards. Figure 2 shows the true stress-true strain curves of IN718 superalloy obtained at the strain rates of 0.001, 0.01, 0.1 and 1 s -1 and temperatures from 940 to 1060°C. As shown in Fig. 2, it is easily found that at the beginning of deformation, the flow stress increases rapidly. It is attributed to the high work hardening rate caused by the rapid generation and multiplication of dislocation. Meanwhile, the dynamic recovery resulted from the dislocation climbing, and sliding is too weak to overcome the work hardening effect. With the increase in strain, the flow stress curves exhibit similar trends, i.e., a single peak at a critical strain followed by a strain softening stage and then sometimes a steady stage at high strain zone. The characteristics of the flow stress curves are typical that can be observed in many metals and alloys. This implies the occurrence of DRX phenomenon during hot deformation.

Flow characteristics
However, at temperatures from 940 to 1060°C and strain rate of 1 s -1 , there is no obvious steady stage in the test strain range.
3 Modeling and simulation on DRX by CA method 3.1 Mechanism of CA simulation of DRX When CA method is adopted, the evolution of the microstructure can be represented by some basic microstructural constituents such as the grain orientation. This means that the actual microstructures can be predicted at mesoscale and the statistical information of internal variables used in the macroscale models are just the inevitable results. The interaction mechanism between mesoscopic and macroscopic evolution for CA simulation is illustrated in Fig. 3. Actually, the dislocation density provides a link of multi-scale modeling to bridge the macroscopic thermomechanical parameters and the mesoscopic microstructures. Figure 3 also schematically shows the variation of flow stress as well as the mean dislocation density of the material and instantaneous dislocation density at a single cell of CA. In the case of none recrystallization, the flow stress depends on the combined effects of work hardening and dynamic recovery, showing a monotone increase toward a saturated stress. But in fact, DRX takes place once the dislocation density reaches a critical value, resulting in a significant decrease in the dislocation density in the recrystallized zone. The difference between the dislocation density of recrystallized grain and the matrix supplies the thermodynamic driving force for the growth of recrystallized grain. In continuous deformation, the dislocation density of a particle changes as the strain increases and the recrystallization happens, which shows a decreased mean dislocation density and a decreased flow stress. Therefore, when the CA method is adopted, one of the most important issues is to track the changes in dislocation density.

Modeling of dislocation evolution
During plastic deformation, dislocations are created by work hardening (WH) and annihilated by dynamic recovery (DRV). To describe the deformation behavior of austenite, the Kocks-Mecking model [27] is employed, in which the flow stress is proportional to the square root of dislocation density as given in: where a is a dislocation interaction term, which is around 0.5 for most metals, b is Burger's vector, l is the shear modulus, q mean is the mean dislocation density. In CA simulation, q mean can be expressed as: where N total is the total number of cells in a CA model and q i,j is the dislocation density of site (i, j). The dislocation hardening and recovery are modeled by the change in the dislocation density exactly as proposed by Estrin [28]. The variation of the dislocation density of site (i, j) with respect to strain can be expressed as: where the two terms represent a competition between storage of mobile dislocations and dynamic annihilation of dislocations, k 1 represents work hardening, k 2 is the softening parameter that represents the recovery of dislocation. k 1 and k 2 in Eq. (3) can be evaluated from the measured flow stress and can be expressed as functions of temperature and strain rate.

Modeling of nucleation and grain growth
The onset of nucleation of DRX is usually associated with the accumulation of dislocations. When the dislocation density of a deforming matrix reaches a critical value Fig. 3 Interaction mechanism between mesoscopic and macroscopic evolution for CA simulation (q critical ) during thermomechanical processing, nucleation of DRX is initiated. q critical is calculated by using Eq. (1). Nucleation of DRX is handled by randomly choosing a certain fraction of the total cells located on grain boundaries (including primary grain boundaries and recrystallized grain boundaries) as the new nucleus at each CA simulation step. Each nucleus consists of one cell. The orientations of the nucleus are randomly set in the program and are represented by random integer in the range of 1-180. When such a nucleus is placed in the system, it is initially assigned a stored energy of zero and increases with increasing strain based on Eq. (3). The nucleation rate per unit grain boundary area for DRX can be expressed as [5]: where C dynamic is a constant, m is set to 1.0 in the present simulation, Q act is the activation energy, _ e is the strain rate, R is the universal gas constant (8.31 J mol -1 K -1 ) and T is the absolute temperature. There are two basic methods to calculate the constant (C dynamic ). First, the percentage of DRX (g) can be experimentally measured for a specific deformation condition [8]; then, the constant C dynamic in Eq. (4) can be determined for a deforming material. Second, the constant C dynamic can be determined using an inverse analysis method [9,14]. The dislocation density of recrystallized grain evolves from zero to a saturated value. The difference between the dislocation densities of recrystallized grain and the matrix supplies the thermodynamic driving force for the growth of recrystallized grain. It is generally assumed that the velocity of grain boundary movement (v) fulfills the following equation [29]: where Df is the driving force and M is the grain boundary mobility, which can be calculated by the following equation: where d is the characteristic grain boundary thickness, D ob is the boundary self-diffusion coefficient, Q b is the grain boundary diffusion activation energy, and K is the Boltzmann's constant. It is assumed that the recrystallized grain is spherical. The driving force Df can be expressed as: where q matrix is the dislocation density of the matrix, q R-grain is the dislocation density of the ith recrystallized grain, s is the dislocation line energy (s = clb 2 ) (c is a constant equal to 0.5). The 2D simulation metallography is viewed as a cross section of the corresponding 3D microstructure. An equivalent average radius of a grain (r i_recrystallized ) is used to define the recrystallized grain size and is calculated by using the following equation: where N i is the total number of cells belonging to the recrystallized grain and where A cell is the area of one cell. Besides, the grain boundary energy (c i ) can be calculated by the Read-Shockley equation. The time step (Dt) is defined as the ratio of the cell diameter (d 0 ) to maximum grain boundary velocity (v max ).

Consideration of grain geometry of real deformation
The deformation of grain geometry makes change of grain boundary and consequently affects nucleation and growth of nuclei during DRX. To describe the effect of compression on the grain shape more accurately, a 2D CA model with the incorporation of a grain morphology deformation technique was adopted in this study. In the modified CA model, a cellular coordinate system and a material coordinate system were established separately. The cellular coordinate system was used to describe the nucleation of DRX and the growth of new recrystallized grains; the size of the cell remains unchangeable. The material coordinate system and the corresponding grain boundary shape change with deformation during a CA simulation. A detail description of grain morphology mechanism in CA simulation is shown in Fig. 4. As shown in this figure (step 1), the original matrix was compressed in the compression direction and elongated in the perpendicular direction during DRX. Then, the grain boundary will be mapped from the material coordinate system to the cellular coordinate system before the appearance of DRX (step 2). The dislocation density increases gradually, and when the value reaches a critical density (q critical ), new nucleus appears on the grain boundaries and continues to grow with the model of equiaxed growth in the cellular coordinate system (step 3). Finally, the grain boundary will be mapped back to the material coordinate system for the next round of the simulation.

Implementation of recrystallization in CA model
In this study, the CA framework for recrystallization is composed of the following parts: (1) Generation of initial microstructure; (2) Definition of the cell neighborhood; (3) Assignment of state variables to each cell; (4) Definition of transition rules for CA models.
Every cell has four state variables: (1) The dislocation density variable represents the stored energy. The initial dislocation density of all primary grains is uniform and identical. During the progress of DRX, the dislocation density increment (Dq) is calculated based on Eq. (3). When a cell recrystallized, the dislocation density of the cell is equal to q initial and increases with increasing strain. indicates the increment of the distance at t À Dt When a cell recrystallized, the recrystallized fraction of the cell equals to 1.
In the CA model for recrystallization, the transformation of a cell from the unrecrystallized state to the recrystallized state can be initiated when the following conditions are fulfilled simultaneously: (1) The driving force is positive; (2) The cell is at the grain boundary; meanwhile, at least one cell in the von Neumann neighborhood is already in the recrystallized state; (3) The recrystallized fraction variable is equal to 1; (4) The computer-generated random number is less than the transformation probability P, which is defined as P = i/4, where i is the number of neighboring cells, which have the same orientation in the von Neumann neighborhood. Figure 5 gives the operation of the CA model. Generally, it includes the following steps: (1) The initial grain is generated by running the initial microstructure generation CA module based on the state transition rules. In this module, two of four state variables are used: the orientation variable and grain boundary variable. As shown in Fig. 5, the initial orientation of the matrix is randomly set as Q (1^Q^180). In each cellular automaton steps (CASs), the entire cells are investigated in accordance with the transition rules. More details on the transition rules can be found in previous work [30]. (2) Calculation of the time step by the following Eq. (9).
The time step is defined as the ratio of the cell diameter (d 0 ) to maximum grain boundary velocity (v max ).
The total number of cycles can be determined by e= _ eDt in the whole deformation.
(3) Calculation of dislocation density increment with the deformation by Eq. (3). When the dislocation density of the matrix reaches a critical value (q critical ) for DRX, the recrystallized nuclei form on the grain boundaries as shown in Fig. 5. (4) Once recrystallization occurs in a cell, the density of the cell comes back to the relative low density, which is equal to q initial in this CA simulation.

CA simulation parameters 4.1 Parameter categories
The parameters used in the modified CA simulation are divided into two categories: (1) The material parameters calculated from the flow stress obtained by compression tests, such as k 1 and  Table 1.

Determination of k 1 and k 2
Equation (3) describes the variation of dislocation density without the occurrence of DRX. Therefore, k 1 and k 2 can be determined from the measured stress-strain data before the onset of DRX. By integrating Eq. (3), the following expression is obtained, where q 0 is calculated using Eq. (1) and r 0 is obtained  from the flow stress. Once dq de ¼ 0, it can be obtained, where q saturation is the saturation dislocation density. Substituting Eq. (11) into Eq. (10), it gives where r WH is the flow stress in the working hardening stage. For a certain deformation condition, substituting r WH , r 0 and r saturation into Eq. (12), it gives the value of k 2 . From Eq. (11), the value of k 1 can be obtained.
To calculate r saturation and q critical on the flow curve, the approach based on the analysis by Kocks and Mecking was employed [27]. Using this method, the critical and saturation stress for dynamic recrystallization can be indirectly determined from the curve of working hardening rate (h ¼ dr=de). Figure 6 gives the stress dependence of the experimental work hardening rate at different deformation conditions. As an example (940°C and 1 s -1 ), the key points have been marked by red circle in Fig. 6.

Calculation of work hardening rate and the activation energy
The work hardening rate can be calculated from the slope of the experimental true stress-true strain curves in the work hardening stage. Recrystallization activation energy can be expressed as [27], where R is the gas constant, a is constant, r p is the peak stress, T is the temperature. The value of n can be derived from the slope in a plot of the value of ln[sinh(ar p )] as a function of ln _ e. Q act can be derived from the slope in a plot of ln[sinh(ar p )] as a function of 1/T. The value of Q act is 430 kJ/mol.

Reverse method of determination of nucleation parameter
The intrinsic information of DRX such as nucleation and growth of the new recrystallized grain is embedded in the material flow stress. The optimal values of nucleation parameter can be identified when the difference between the measured and calculated stress-strain curves is minimized. Essentially, it is a minimization problem and the minimization objective function can be defined as: where r exp k and r cal k x ð Þ are experimental and calculated flow stresses at the strains e k (k = 1, …, N), respectively. N is the number of sampling points on the stress-strain curve and x is the nucleation parameters (C dynamic ) to be identified. Table 2 shows the values obtained from flow stress curves.

Simulation results and discussion
In the developed CA model, the simulation mesh is 512 9 512 square lattice, and periodic boundary conditions are used. The size of every cell represents 2 lm of the real dimension of the material. In order to balance the accuracy and computational efficiency, the increment of strain is 0.2 during coupling with the effect of grain morphology change on DRX. Figure 7 shows a comparison of predicted true stress-strain curve and experimental one. It also gives the microstructural evolution of DRX at a strain rate of 1 s -1 and temperature of 1060°C with consideration of grain morphology deformation. The different color regions represent the recrystallized grains, and the orange regions denote the deformed matrix, and the black lines denote the grain boundaries (GBs). The initial microstructure (D 0 = 210 lm) is obtained by running the initial microstructural generation module. As shown in the figure, the grains that have not completed the DRX in each round are compressed in the vertical direction with the continuous deformation, and the newly formed recrystallized grains grow with the type of equiaxed growth in the cellular coordinate system. The recrystallized fraction increases with the increase in the strain from 0.62 to 0.93.  Figure 8 gives the simulated microstructures without considering the grain morphology change of the initial matrix at the temperature of 1060°C and the strain rates of 1 and 0.1 s -1 , respectively. Figure 9 shows the simulated results by using the CA model with support of the grain morphology technology under the same deformation conditions. Additionally, Fig. 10 illustrates the simulated average grain size and recrystallized fraction using the topological CA model at a strain rate of 1 and 0.1 s -1 and temperature of 1060°C in comparison with the conventional CA simulation. Overall, there is a reasonable agreement between the final average grain size in CA simulations and experiments from Figs. 8 and 11. Meanwhile, the results also show that the grain morphology deformation has an accelerated effect on the recrystallization process. The recrystallized fractions are 89 and 75 % for the strain rates of 0.1 and 1 s -1 , respectively, when the grain morphology deformation technique was adopted as compared to only 58 and 46 % of recrystallized fractions without the consideration of grain morphology deformation. The simulation results also indicate that the simulated average grain size decreases slightly due to the effect of matrix deformation. The reason for such a difference may be because that the elongated grain boundaries provide more potential nucleation sites for recrystallization due to a higher ratio of the grain boundary volume to grain volume.  At the same time, the newly formed recrystallized grains impinge earlier with support of the morphology change of the matrix, resulting in a smaller average grain size as compared to the conventional CA simulation [5][6][7][8][9][10][11][12][13][14][15][16][17][18]. In addition, it can be also found that the grain morphology deformation has a decreased effect on the critical strain for DRX. The reason may be due to that the reduction in the initial grain size/matrix accelerates the progress of DRX [5]. In conclusion, the CA model with the incorporation of the morphology deformation technique can better reflect the topological characteristics of the deformed grains in the simulation of DRX. Figure 11 shows the simulated microstructures of IN718 superalloy under different thermomechanical conditions. The average size of the DRX grains increases from 2 to 12 lm with the temperature increases from 940 to 1020°C. The average grain size increases from 5 to 8 lm with strain rate decreasing from 0.01 to 0.001 s -1 . It is obvious that decrease in temperature and/or increase in strain rate can refine the DRX grain size. The influences of the strain rate and the temperature on the DRX kinetics are presented in Fig. 12. Figure 12a shows the DRX kinetics at different temperatures under the strain rate of 0.01 s -1 , while Fig. 12b shows the DRX kinetics at 980°C under   Fig. 12 a Simulated kinetic curves of DRX at the strain of 1 and strain rate of 0.01 s -1 , and b simulated kinetic curves of DRX at the strain of 1 and temperature of 980°C different strain rates. It is clear that only when the strain reaches a critical value, DRX can occur at a certain temperature and strain rate under hot deformation. In other words, there is a need for time-dependent deformation to initiate DRX. Also, this time-dependent deformation would account for the incubation period that was observed by Fisher et al. [33]. From these figures, we can also observe that increasing the temperature and/or decreasing the strain rate can reduce the incubation period. This phenomenon is mainly attributable to the rate of dislocation accumulation and the critical strain. The faster the rate of dislocation accumulation and the smaller the critical strain, hence the shorter the incubation time.  Fig. 13a, f). Meanwhile, it can be found that number of DRX grains obviously decreases under the strain rate of 1 s -1 , and the elongated original grains is shown easily in Fig. 13b. Also, it reveals that some fine DRX grains preferentially nucleate at the vicinities of grain boundaries. From Fig. 13a, b, d and f, it can be found that the DRX  Figure 13c, d and e shows the grain microstructures in the deformed specimens under 940, 980 and 1020°C, respectively. Here, the true strain and strain rate are 1 and 0.01 s -1 , respectively. It is clearly seen that the grain microstructures are greatly dependent on the deformation temperature. Under the relatively low temperature (940°C), there are many elongated original grains, but only a few DRX grains, as shown in Fig. 13c, and the recrystallized grain size is less than 2 lm. When the deformation temperature increases to 980°C, some equiaxed DRX grains appear and the average recrystallized grain size can be estimated as 6-8 lm, as shown in Fig. 13d. Under the deformation temperature of 1020°C, the fraction volume of recrystallized DRX grains increases greatly and the grain microstructure is notably reined (shown in Fig. 13e). As compared with the CA simulated microstructures (shown in Figs. 9 and 11), the simulated recrystallized DRX grain size values are slightly higher than the experimental values. The main reason may conclude that the modified CA simulation did not consider the effect of precipitation phases such as c 0 , c 00 and d on the DRX behavior of Ni-based superalloy [34]. Figure 14 shows the microstructural evolution during large strain deformation. As shown in the figure, when the strain reaches 1.0, the recrystallized fraction is less than 50 %, when the strain reaches 1.5, complete recrystallization takes place, and the initial grains are refined greatly. It is well known that one of the typical characteristic of severe plastic deformation (SPD) is refining the grain size through the large deformation [35]. Therefore, CA method is capable of predicting microstructural evolution for DRX of IN718 superalloy in SPD process such as multi-direction forging.

Conclusions
In this paper, the 2D CA model was developed to study on the DRX behaviors of IN718 superalloy. The main conclusions can be summarized as follows: (1) By introducing the mechanism of morphology change of grain, a CA model was established to describe the DRX for IN718 superalloy. At the same time, the reverse method was adopted to identify the nucleation parameter. The CA model is shown to be a new and powerful tool for quantitative and explicit simulation of microstructural evolution during DRX. The established CA model can well describe the DRX behaviors of the studied superalloy. Furthermore, the developed CA model can also be used as an efficient simulation method to investigate the microstructural evolution in large material deformation on a mesoscale under the conditions such as multi-direction forging for IN718 superalloy. (2) The effects of deformation temperature and strain rate on the recrystallized grain size are significant. Increasing the deformation temperature or decreasing the strain rate coarsens the recrystallized grain size. DRX is difficult to take place under the relative low temperature (900°C). At a relative low strain rate of 0.001 s -1 and temperature higher than 980°C, it is easier for full DRX to take place.