A three-dimensional (3D) cellular automaton model describing the austenite grain growth process of Fe-1C-1.5Cr alloy steel is established. This model proposes a discrete representation of the three-dimensional grain boundary curvature, and calculates the 26 neighbor lattice points to the center lattice point in three-dimensional space. The strength of neighbor lattice point’s influence is determined by the distance from the center lattice point. The paper presents the calculation steps of the three-dimensional cellular automaton model and the calculation methods of key parameters such as grain boundary curvature and grain growth speed, and validates the model calculation results. The results show that the model established in this paper can well describe the austenite grain growth process of Fe-1C-1.5Cr alloy steel at different temperatures. The model can be used to study the grain growth process of Fe-1C-1.5Cr alloy steel under three-dimensional conditions.

Fe-1C-1.5Cr alloy steel is a typical bearing steel, it has high fatigue resistance, high ductility and good wear resistance, and it has wide application in industry. The austenite growth process of Fe-1C-1.5Cr alloy steel has an important impact on the quality of products.

The austenite matrix is composed of crystal grains as a basic unit. Austenite grains tend to grow under the driving of grain boundary energy. During grain growth process, the grain boundary area decreases and the energy of the entire system decreases. There are two trends in the understanding of the mechanism of grain growth. One believes that grain growth is carried out by a combination mechanism, that is, the grain boundaries of adjacent grains disappear due to decomposition of large-angle grain boundaries. Another reason is that the grain growth process is achieved by grain boundary migration. For the grain boundary migration, scholars believe that the grain boundary has two structures of dislocation and step. The movement of the grain boundaries is accomplished by dislocation motion or dislocation step motion.

Many scholars have conducted a large number of experimental studies on the process of grain growth. Kurtz et al. [1] and Antonione et al. [2] studied the grain growth phenomenon during isothermal annealing, and obtained the grain size distribution in accordance with the log-positive distribution. Hu et al. [3] systematically studied the grain growth kinetics of pure iron, and obtained the grain growth index with the increase of temperature, and in general, lower than the theoretical value of 0.5. Gil et al. [4] studied the grain growth process of pure titanium and found that the grain growth index is close to the theoretical value of 0.5. Hu1l [5] studied the grain size and morphology of β-brass and found that the grain size is linear with the number of grains. In recent years, the automatic serial sectioning for three-dimensional analysis of microstructures introduced by Spowart et al. [6] can quickly obtain the three-dimensional structure of the grain. Ji et al. [7] indicate that the austenite grains of both steels grow up gradually with increasing the heating temperature and holding time; the size and growth rate of austenite grain of Nb–Ti-bearing high carbon steel are much lower than those of Nb–Ti-free high carbon steel. Xu et al. [8] proposes a coupled multiphysical reaction-transport (c-MPRT) model that is based on the implicit finite difference method for the promising bisolid porous pellet to investigate its dynamic characteristics concerning heat transfer, chemical reaction, and phase change. This model successfully reveals a negative influence on the heat conduction from the outer layers to the inner layers of the pellet, due to the enormous heat absorption of the chemical reacting and physical melting process.

In addition to experimental research, many scholars have proposed theoretical models of grain growth. The models in this area mainly include two types, one is the classical dynamic model, including the Beck model [9] and the Sellars model [10]. The other type is a model that simulates the microstructure of a material, including the Vertex or Network model [11], the front tracking model [12,13], and the phase field model [14,15], Monte-Carlo Potts model [16] and cellular automata model [17]. The classical kinetic model can describe the average size of the grains, but for the morphology of the grains, the distribution of grain sizes, and the process of abnormal grain growth, there is no power.

The cellular automaton model is a research method that discretizes both space and time. By superimposing the local evolution process, the change law of the whole system is obtained. There have been many studies on grain growth based on two-dimensional cellular automata [18–26]. However, the study of grain growth based on three-dimensional cellular automata is still rare [27–30]. The three-dimensional model is closer to the actual grain growth process and the austenite grain growth process can be well described by developing a precise 3D cellular automaton model.

2Three-dimensional simulation of austenite grain growth processBased on the Kevin Kremeyer’s two-dimensional model, this paper proposes a three-dimensional method of discretization of grain boundary curvature. Based on this method, a three-dimensional cellular automaton model of austenite grain growth is established.

2.1Mathematical model of austenite grain growth processIn the cellular automaton model, the three-dimensional computational region is discretized into square lattice points, each lattice point representing a small region in the actual material. Each lattice point is given a crystallographic orientation, and a grain is composed by adjacent lattice points that have the same metallurgy. The grain boundaries are exist between adjacent grains, and the growth process of the grains is determined by the migration of grain boundaries. The cellular automaton model is asynchronous updating in time and local evolution in space. In the three-dimensional cellular automaton model, the state of the specific lattice point in the next time is determined by the state of the adjacent lattice point at this moment. Fig. 1 is a schematic diagram of the specific lattice point and its three-dimensional neighbors. It can be seen from the figure that there are a total of 26 neighbor lattice points for a certain lattice point. There are 3 types of neighbors. The neighbor I is the lattice point that is in face contact with the center lattice point O, and we call those neighbors the Neumann type neighbors that are shown in yellow in the figure, there are 6 neighbor I lattice points in total. Neighbor II is in line contact with the center lattice point O. The lattice points are shown in black in the figure, and there are 12 neighbor II lattice points in total. The neighbors III are neighbors that are in lattice point contact with the center lattice point O. The lattice points is shown in gray in the figure, and there are 8 neighbors III lattice points in total.

This paper considers the influence of all 26 neighbor lattice points on the center lattice point, and the strength of its influence is determined by the distance from the center lattice point. The grain growth process is controlled by the grain boundary moving speed, the frequency of the crystallographic orientation transform into a neighboring lattice point can be descript as follow:

Where i, j represents the coordinates of the center lattice point; m, n is the relative coordinate of the neighbor lattice point and the center lattice point; Vi, j is the moving speed of the grain boundary adjacent to the lattice point of the coordinate (i, j) ; Lm, n is the distance between the neighbor lattice point and the center lattice point. If the size of the lattice point is A, the distances of the neighbors I, II, and III from the center lattice point are A,2A and 3A, respectively.In the absence of strain, the grain growth process is driven by the grain boundary curvature. According to Burke and Turnbull’s theory [31], the grain growth rate can be expressed by the following equation:

Where M is the mobility of the grain boundary; P is the driving force of the grain boundary migration. M is a parameter related to the crystallographic angle and can be expressed by the following equation [32,33]:Where M0 is the pre-factor. Qb is the grain boundary moving activation energy. θ is the angle between adjacent crystallographic orientations. θm is the crystallographic orientation difference of the grain boundary of the large crystallographic angle. T is the absolute temperature. R is a general gas constant.The driving force P of the grain boundary migration can be expressed by the following equation:

Where γ is the grain boundary energy. κ is the grain boundary curvature. γ is also related to the angle of the crystallographic orientation θ and can be calculated using the Read–Shockley equation [34].Where γm is grain boundary energy of large crystal grain boundary.2.2Calculation of grain boundary curvatureIn the calculation of two-dimensional grain boundary curvature, the equation proposed by Kevin Kremeyer [35] in calculating the solidification process is often used:

Where a is the side length of the square lattice point, A = 1.28 is the shape factor, N = 24 is the total number of the first adjacent lattice point (neighbor I) and the second adjacent lattice point (neighbor II), and Ni is the number of lattice points that are the adjacent lattice points of the grain i. Kink = 15 is the number of lattice points belonging to gain that includes the center lattice point when the grain boundary is a flat grain boundary (κ = 0).In the above two-dimensional local curvature calculation model, it is considered that there is a linear relationship between the dimensionless curvature κa and the outer lattice point fraction Kink−NiN+1, and the linear parameter is the shape factor A. The reciprocal of the dimensionless curvature 1/κa is the equivalent radius of the circle in which the surface is located. A suitable shape factor A can be obtained to obtain a equation for calculating the two-dimensional curvature. This method has achieved good results in the two-dimensional CA model [23]. This paper assumes that κa can still be expressed as the product of Kink−NiN+1 and shape factor A in the three-dimensional cellular automaton model.

The first adjacent lattice point and the second adjacent lattice point of the center lattice point are also taken in the three-dimensional cellular automaton model curvature calculation, as shown in Fig. 2. The object in the figure labeled F is the first adjacent lattice point, and the object labeled S is the second adjacent lattice point. N=124 is the total number of three-dimensional neighboring lattice points, Ni is the number of lattice points of the three-dimensional neighboring lattice points belonging to the center lattice point i, and Kink = 75 is the number of lattice points of the adjacent lattice points belonging to gain that includes the center lattice point when the grain boundary is a flat grain boundary (κ = 0). When the side length of the lattice point is determined, the only unknown parameter is the shape factor A.

In this paper, the shape factor A at different equivalent radii is obtained by the following simulation calculation steps:

- 1)
Generating a sphere with equivalent radius R in a discrete three-dimensional lattice point region;

- 2)
Set a shape factor A, and use the Eq. (7) to calculate the dimensionless curvature κa1;

- 3)
Compare κa1 with the actual dimensionless curvature 1/R. If κa1 = 1/R, the shape factor set in step (2) is an appropriate value, otherwise, the shape factor A is adjusted, and the calculation process of step (2) is continued.

As shown in Fig. 3, the values of the shape factor A under different equivalent radii (represented by the dimensionless curvature in the figure) are calculated separately. It can be seen from the figure that in the three-dimensional case, there is a quadratic function relationship between the dimensionless curvature and the shape factor:

Bringing the Eq. (8) into the Eq. (7), the three-dimensional dimensionless curvature satisfies the following quadratic equation:

Where M=Kink−NiN+1.By solving Eq. (9), we can get the three-dimensional curvature to satisfy the following equation:

In Eq. (10), it is guaranteed that the solution of the quadratic equation needs to satisfy M > 0.0524. Also ensure κ > 0 and M > 0.0643. Taking the intersection of two conditions gives a range of M > 0.0643. When M < 0.0643, κ = 0 is considered.

2.3Calculation step of Cellular automaton modelIn this paper, the C#.net is used to program the three-dimensional cellular automaton model of the austenite grain growth of Fe-1C-1.5Cr Alloy Steel, and the simulation is carried out using the self-compiled program. The initial grain size of the austenite structure is generated by the Monte Carlo method. The calculation steps of the cellular automaton model established in this paper are as follows:

- (1)
Calculate the curvature of the grain boundary and the grain boundary moving speed based on the curvature. The frequency vi,j,m,n that the specific lattice point is changed to a neighbor lattice point that’s crystallographic orientation is different are respectively calculated. The total frequency at which the specific lattice point transitions can be calculated by following equation:

- (2)
Take the maximum value vmax among the transformation frequencies of all the grain boundary points, and set the calculated time step Δt = 1/vmax. This time step setting ensures both the accuracy of the calculation and the calculation speed.

- (3)
Calculate the probability that a specific lattice point crystal orientation will change to a neighboring lattice point with a different crystallographic orientation fi,j,m,n=ΔtVi,j/Lm,n. All grain boundary points are updated according to this probability.

- (4)
Update the calculation time and return to step (1).

The austenite grain growth process of Fe-1C-1.5Cr alloy steel (The composition of Fe-1C-1.5Cr alloy steel is shown in Table 1) was calculated using the above model in this paper.

The grain boundary energy γm of the large-angle grain boundary is 0.56 J/m2 mentioned in the references [36]. The mobility parameters were obtained by summarizing the experimental data, taking M0 = 2.82 × 10−8 m4J−1s−1, and the growth activation energy was 142.58 kJ mol−1. The comparison between the calculated results and the experimental results [37] is shown in Fig. 4. It can be seen that the calculation results are in good agreement with the experimental results, the maximum absolute error of the simulation and experimental results is about 7.2 μm, the average absolute error is about 2.7 μm, the maximum absolute error is about 9%, and the average relative error is about 4%. The CA model of this paper can better describe the growth process of the grains.

4Comparison with other modelsIn general, the Burke dynamics relationship is satisfied between the grain radius and the growth time:

Where r and r0 are the grain radius and the initial grain radius, respectively. t is the grain growth time, k is a constant, and m is the grain growth index.The grain growth index m is different in different experiments and simulations [14]. In general, when the influence of the second phase particles on grain growth is neglected, the grain growth index is approximately 2. As shown in Fig. 5, the r2 value obtained at different temperatures and the grain growth time were linearly fitted. The fitting accuracy was found to be above 0.995, indicating that the grain growth index of the simulation results is approximately 2 in this paper. The results are consistent with the simulation results of most two-dimensional and three-dimensional models, indicating that the proposed model can better describe the grain growth process.

In two-dimensional and three-dimensional simulations, the three most commonly used functions for describing grain size distribution are the log-normal function [38], the Louat function [39], and the Weibull function. Fig. 6a shows the simulation results of grain size distribution at 1200 °C after holding different time (1 h, 2 h and 3 h) and Fig. 6b shows the simulation results of grain size distribution at different temperatures (1000 °C, 1100 °C and 1200 °C) after holding 1 h. It can be seen from the figure that the Weibull function can better fit the grain size distribution with respect to the other two functions in the case where the tissue is uniform and not affected by the second phase particles. The holding time and the temperature have little effect on the grain size distribution. The average α = 1.15, β = 2.79 obtained by fitting the Weibull function. This result is closer to the results of the study by He et al. In their study, β = 2.74 [20] in the two-dimensional case and α = 1.08, β = 2.82 [29] in the three-dimensional case. It can also be seen from the figure that the peak value of the crystal frequency appears at a relative size r/ra ≈ 1.0, indicating that most of the grain size is concentrated near the average grain size. The maximum r/ra obtained in this paper’s model is about 2.0, which is close to the results of the two-dimensional model 2.2 [40] and 2.14 [41]. And the simulation results are quite different with the traditional MC method 2.7 [42] and the simulation results of the three-dimensional cellular automaton model which is based minimum energy principle 2.6 [28].

The kinetic equation of grain growth studies the growth process of all grains. However, the growth process of individual grains is related to grain size and grain shape. The study of average grain growth dynamics cannot be a substitute for the dynamics of a single grain growth. Hillert [43] proposed a relationship between the velocity of a single grain growth and the grain radius:

Where c is a constant, r is the grain radius, and rcr is the critical grain radius.In general, in the case of two-dimensional, the critical grain radius is equal to the average grain radius ra, and in the case of three-dimensional, the critical grain radius is equal to 9ra/8. As shown in Fig. 7, the growth process of three grains is calculated in this paper, and it is found that the grain growth relationship proposed by Hillert is generally satisfied, especially when the grain radius is close to the critical grain radius. When the grain radius is larger than the critical grain radius, the growth velocity is substantially greater than 0, and when the grain radius is less than the critical radius, the growth velocity is substantially less than zero. This phenomenon indicates that the grain growth process is a process in which large grains engulf small grains.

Ding et al. [28] summarized the relationship between grain volume and the number of grain faces in the calculation results of a three-dimensional cellular automata model based on the principle of minimum energy.

Where Nf is the number of faces of a certain grain, Vn is the volume of the grain, Va is the average volume of all grains, and α is a constant.In this paper, the simulation results of 1 h, 2 h and 3 h respectively at 1100 °C were counted, and the results are shown in Fig. 8. It can be seen from the figure that the number of gain faces Nf and (Vn/Va)2/3 are the basic satisfying linear relationship and are the same as the results of Ding. The intercept of the linear equation is 6.84, which is basically equal to the calculation result of Ding (6.81), indicating that the average number of grain faces is between 6 and 7 when grains is disappearing. In addition, when the grain volume is equal to the average grain volume, the corresponding number of grain faces is 16 in this paper, which is slightly different from the calculation result 14 of Ding.

5ConclusionIn this paper, a new discretization expression of three-dimensional curvature is established, and a three-dimensional cellular automaton model of austenite grain growth is established. In order to verify the reliability of the model, the experimental results was used to fit the calculation results. It was found that the grain growth index calculated by the model was 2 and the precision was above 0.995, which was consistent with the previous research results: the Weibull grain size distribution function, the single crystal length dynamic equation proposed by Hillert and the three-dimensional grain topological relationship proposed by Ding et al. The results show that the three-dimensional cellular automaton model proposed in this paper has a high consistency with some previous theoretical research results and experimental results on austenite grain growth, and can successfully simulate the three-dimensional austenite grain growth process.

Conflicts of interestThe authors declare no conflicts of interest.

This work is supported by National Natural Science Foundation of China (51706017).

*et al*.

*et al*.

*et al*.

*et al*.

*et al*.