Jul 15, 2021 A novel 3D mixed-mode multigrain model with efficient implementation of solute drag applied to austenite-ferrite phase transformations

A computational 3D model that accounts for both nucleation and interface migration is a very useful tool to monitor and grasp the complexity of microstructure formation in low-alloyed steels. In the present study we have developed a 3D mixed-mode multigrain model for the austenite-ferrite and the austenite- ferrite-austenite formation capable of following diffusional phase transformations under arbitrary thermal routes. This new model incorporates the solute drag effect of a substitutional element (in this case Mn) and ensures an automatic change in transformation direction when changing from heating to cooling and vice-versa. An analytical solution for calculating the energy dissipation of solute drag together with mul- tiple regression approximations for chemical potentials are proposed which signiﬁcantly accelerate the computation. The modelling results are ﬁrst benchmarked for an Fe-0.1C-0.5Mn (wt.%) alloy under differ- ent continuous cooling and isothermal holding conditions. The model revealed relatively large variations in transformation kinetics of individual grains as a result of interactions with neighboring grains. Then the model is applied to predict the transformation kinetics of a series of Fe-C-Mn alloys during cyclic partial phase transformations. The comparison with experimental dilatometer results nicely validates the predictions of this model regarding the change in overall transformation kinetics of the ferrite transformation as a function of the Mn content. New features of this model are its eﬃcient algorithm to compute energy dissipation by solute drag, its capabilities of predicting the microstructural state for spatially resolved grains and the minimal ﬁne tuning of modelling parameters. The code to implement this model is publicly available. Acta


Introduction
The kinetics of the austenite ( γ ) to ferrite ( α) phase transformation in steels has been studied for decades using both experimental and modelling approaches [1,2] . More recently, an increasing interest in the kinetics of the ferrite-to-austenite transformation has been observed, motivated by the design and production of advanced Mn-rich high-strength steels [3][4][5] . It is well known that common substitutional alloying elements such as Mn, Ni and Mo strongly retard the interface migration during the austenite-ing or more complex thermal routes generally lies in between the predictions for local equilibrium (LE) and para-equilibrium (PE). Understanding the effect of the slower diffusional substitutional elements on the transformation kinetics is essential to be able to tune the microstructure. So far, it has remained challenging to predict and quantify the interface migration for an arbitrary imposed thermal profile involving both austenite and ferrite.
There are two main approaches established to account for the effect of substitutional elements on the interface migration. The first approach uses the concept of an effective interface mobility that considers the effect the alloying element would have on the character of interface mobility, resulting in a decrease in the value of the interface mobility compared to the value of intrinsic γ / α interface mobility [17][18][19][20] . This approach is demonstrated to be efficient in computations, but requires a good estimation of the influence on the effective mobility beforehand. As the effective interface mobility can vary with the steel compositions, the cooling rate and the transformation direction, the value of the prefactor M 0 in the Arrhenius expression of the effective interface mobility can vary over several orders of magnitude from 10 -10 to 10 −6 molmJ −1 s −1 [21][22][23][24] . This wide range of values for the effective interface mobility generally requires repeated trials before it can be used for accurate predictions. Therefore, the general application of this approach is limited. The alternative approach is to adopt the solute drag theory to account for the effect of a substitutional element, while using the intrinsic interface mobility to account for the energy consumption by the interface friction [25][26][27][28][29][30] . In the solute drag theory, the substitutional element can segregate at the interface, causing energy dissipation by trans-diffusion inside the interface. The driving force for the phase transformation is potentially consumed by both the solute trans-diffusion inside the interface and the interface friction. As the intrinsic γ / α interface mobility can be derived experimentally from the massive phase transformation in binary Fe-X (X = Mn, Ni, Mo etc.) alloys, it could be regarded as generic for austenite-ferrite phase transformations, regardless of the changes in nominal compositions and cooling/heating rates. A recent study on a series of Fe-Ni, Fe-Mn and Fe-Co alloys provides a good reference for the parameters that determine the intrinsic interface mobility [31] . One other parameter one needs to know to calculate the energy dissipation due to solute drag is the binding energy of a specific substitutional element, which describes the affinity of the solute species to the austenite-ferrite interface. Although an accurate determination of the binding energy is not straightforward, this parameter can be estimated from detailed measurements of the composition profile across the interface by atom probe tomography [32,33] or from first-principles calculations [34] . The solute drag effect of different alloying elements can now be included and thus leads to a coupled solute drag effect, which has been demonstrated to perform well in the recent work on Fe-C-Mn-Mo alloys [13] . Therefore, the approach that uses the solute drag theory to account for the effect of substitutional elements on interface migration shows a better transferability than the one employing the effective interface mobility.
As the transformation kinetics is not only controlled by interface migration, but also by nucleation of the new phase, the evolution of the new and parent grain structure is a result of both effects. Thus, the effect of alloying elements on the interface migration cannot solely be studied by conventional continuous cooling or heating, or by isothermal holding experiments in which nucleation and growth of the ferrite phase occur simultaneously. An effective approach to circumvent this issue is to adopt the socalled cyclic partial phase transformation, where the steel is thermally cycled between two temperatures, that both lie in the twophase γα region, so that only the more energetically favored in-terface migration takes place, rather than the simultaneous occurrence of interface migration and nucleation of new grains [11] . This behavior has been observed by our recent neutron depolarization measurements, which demonstrate that new nucleation of ferrite grains is indeed absent and that the transformation kinetics is only controlled by the interface migration during cycling [35] . Although all cycling experiments of Fe-C-Mn alloys (compositions with 0.023-0.25 wt.% C and 0.17-2.1 wt.% Mn) show stagnant stages in the transformation of austenite to ferrite and vice versa, the cycling behavior is distinctly different. For example, in a lean Mnalloyed Fe-0.023C-0.17Mn (in wt.%) steel an open loop of the ferrite fraction is observed [11] , whereas a gradual net increase of the ferrite fraction is present in an Fe-0.25C-2.1Mn (in wt.%) steel during cycling [35] . Another study on austenite formation during thermal cycling in a medium-Mn steel (Fe-0.2C-4.5Mn in wt.%) showed that the fraction of the newly formed austenite decreases in each successive cycle, whilst the total amount of austenite increases during cycling [36] . These differences in behavior are due to the dynamics of the modifications of the local chemical compositions across the interface, the element chemical potentials and the relative velocities of the Mn diffusion and the moving interface. One of the best ways to clarify this complexity in dynamics is to develop models that reveal the effects of local chemistry (including C and Mn) and phase structure at the same time under a mixed-mode interface condition. To do this, in recent years extensive modelling studies have been conducted to investigate the kinetics of cyclic partial austenite-ferrite phase transformations using the solute drag theory [37][38][39][40][41][42] . However, given the heavy computational demands after incorporating solute drag, most established models are limited to 1D, with the exception of one 2D cellular automaton model [41] , and all of these models only focus on the interface migration itself. Although these studies provide a detailed insight into the transformation behavior during linear cooling, isothermal annealing and thermal cycling, two key parameters that need to be included to make the step to model transformation kinetics in real (3D) steels are still missing. Firstly, the 1D simulation provides no spatial information and gives no insight into the variation in behavior of different f errite grains. Secondly, the varying degrees of nucleation (i.e. degree of undercooling) of individual ferrite grains before entering thermal cycling has been discarded, and as a result the previous models contain no information on the grain size and the grain size distribution. However, as pointed out in our previous study [43] , the average grain size and the grain size distribution are essential information for a model to predict the microstructure and shed light on the underlying physics. Therefore, it is very desirable to develop a model that not only incorporates solute drag, but also includes the spatial information in 3D to generate a statistically relevant number of events for ferrite formation. A model of this kind will also make it possible to study how the size distribution of ferrite grains formed before cycling affects the interface migrating back and forth in the presence of solute drag.
In the present work, we have extended our previous 3D multigrain mixed-mode model by coupling it with the solute drag theory. This new model accounts for ferrite nucleation based on the classical nucleation theory and interface migration based on the balance in Gibbs free energies between the chemical driving force and the energy dissipations due to interface friction and solute drag, all calculated locally for each moving interface. The modelling results are first benchmarked for an Fe-0.1C-0.5Mn ternary alloy under continuous cooling and isothermal holding. Next, this model is applied to describe the transformation kinetics during cyclic partial phase transformations of the same steel. Finally, the model is applied to calculate the cyclic transformation behavior as a function of C and Mn concentrations in the steels. The aim is to develop a versatile and flexible modelling tool to pre-2 dict the kinetics of austenite-ferrite phase transformations with a minimal number of fitting parameters. The more efficient method to compute the local energy dissipation due to solute drag proposed in the present work is also expected to provide inspiration for further 3D simulations of multigrain solid-state phase transformations.

Model
The starting structure in the model is a fully austenitic multigrain structure produced by constructing Voronoi cells in a cubic box with a length of L b . The number density of austenite grains ( ρ γ ) is preset, resulting in an average austenite grain size d γ = ( 6 /π ρ γ ) 1 / 3 . The centers of the Voronoi cells are randomly produced with an imposed minimum distance (d min ) to control the austenite grain size distribution. The corners of the Voronoi cells are regarded as the main potential nucleation sites for ferrite grains, as ferrite nucleation is found to predominantly take place at grain corners. Besides the grain corners, also edges and faces can act as potential nucleation sites with a lower probability [44] . Once a stable ferrite nucleus forms at one of the grain corners (or alternatively at a grain edge or at a grain face), the ferrite grain is assumed to grow isotropically into the austenite as a sphere. As shown in the micro-beam X-ray diffraction during austenite-toferrite and ferrite-to-austenite transformations [45,46] , the nucleation barrier is apparently higher than the kinetic energies. Therefore, it is reasonable to use classical nucleation theory (CNT) to compute the number density of ferrite grains [47][48][49] . The migration of the austenite-ferrite interface is controlled by the diffusion of both carbon and substitutional elements and by a lattice reconstruction. This is a so-called mixed-mode approach [1,3,18] and the relevant assumptions of the model as used are described in detail later. The solute drag theory is applied to describe the energy dissipation due to the trans-diffusion of substitutional solute inside the interface. A quasi-steady state is reached when the pressure imposed on the moving interface is zero, indicating that the chemical driving force is balanced by the energy consumption by both interface friction and solute drag. The above principle has been widely adopted in various modelling approaches and currently sets the standard. However, the accuracy of the model depends on being able to calculate the concentrations of the key interstitial and substitutional alloying elements at moving interfaces as well as in the parent phases the interface is moving towards.
In a previous work we derived analytical expressions for the carbon (i.e. the interstitial element) concentration at the interface and far away from the interface for non-overlapping and overlapping of diffusion fields (soft impingement). The base model also accounted for the hard impingement of growing ferrite grains. These expressions were shown to be correct regardless of transformation directions (i.e. ferrite grains are growing or austenite grains are growing).

Ferrite nucleation
In the simulations to be presented, ferrite nucleation was assumed to take place only at the corners of austenite grains. The code has however been generalized such that it can handle different nucleation sites, such as grain corners, grain edges and grain faces or mixtures of these. Dedicated simulations showed that the choice of the dominant nucleation site had a marginal effect on the transformation and does not to affect the main conclusions of the present work (minor changes were only observed in the later stages of the transformation). The classical nucleation theory is used to describe the nucleation rate dN/dt in a selected volume. According to the CNT, the nucleation rate can be written as: where A is a pre-factor, Z is the Zeldovich factor and nearly constant (Z ≈ 0.05), N 0 represents the number of potential nucleation sites when transformation starts, f α is the volume fraction of ferrite phase, k B is the Bolzmann constant, T is the temperature in Kelvin, h is Planck's constant, Q D is the energy barrier for diffusion, is a constant that comprises all the contributions of the shape of the critical nucleus and interfacial energy between the nucleus and the surrounding parent grains, G V is the net difference in Gibbs free energy per unit volume between ferrite and austenite, G S is the misfit strain energy per unit volume. The critical nucleus size is [49] , where σ αγ is the interfacial energy per unit area for γ / α boundaries and amounts to σ αγ = 0.62 Jm −2 [50] . In Eq. 1 the constant is the most challenging parameter to be determined precisely. However, high energy X-ray diffraction studies monitoring the nucleation of ferrite or austenite suggest that in reality has a value of the order of 10 -8 J 3 m −6 [45] . Therefore, a value of ≈ 5 × 10 -8 J 3 m −6 is used in the present work. The G S is calculated by G s = 2( 1 −ν ) [51] , where ν is the Poisson's ratio, μ is the shear modulus ( ν = 0.33 and μ = 60 GPa for pure iron). The lattice parameter of ferrite a α and the lattice parameter of austenite a γ both depend on temperature and carbon concentration, and are calculated according to [52] .

γ / α interface migration
The chemical driving force per mole of atoms G chem m for the γ / α interface migration after ferrite nucleation can be expressed as: where subscript i is the element in the alloy, n is the total number of elements, x 0 i is the composition of element i transferred across the interface, the superscripts of γ α and αγ correspond to the austenite side and ferrite side on the interface, respectively, μ is the chemical potential and x is the mole fraction. For an interface moving in a quasi-steady state the chemical driving force is balanced by the interface friction G friction m and the dissipation of substitutional elements inside the interface G diff m : Such an approach to split the dissipation of the chemical driving force into two contributions was first proposed by Hillert and coworkers [25,53,54] and is now a common approach (e.g. refs. [13,30,37] ). Eq. 3 also indicates that long-range diffusion of a substitutional element (Mn in this case) in ferrite or austenite is not considered, although it may contribute to the Gibbs free energy dissipation. According to [55] , it is not evident how to separate the energy dissipation due to trans-diffusion of substitutional element inside the interface from the dissipation by the long-range diffusion. Since interstitial atoms diffuse much faster than substitutional atoms, it is reasonable to assume that there is a long-range diffusion of the interstitial element (C in this case) and a steadystate condition of the substitutional element (Mn in this case) in and across the thin interface region [56] , as is done in the current model.
The interface friction depends on the interface velocity v int and the intrinsic mobility of the γ / α interface M int : where V m is the molar volume of Fe and assumed to be the same for austenite and ferrite, M int is expressed by M int = Q int is the activation energy and R is the gas constant. The use of the intrinsic interface mobility, rather than the effective interface mobility, provides a substantial advantage in transferability of the model from one steel composition to another. The dissipation energy due to trans-diffusion of solute inside the interface G diff m is calculated using the method proposed by Purdy and Bréchet [57] . As G diff m depends on v int , which is related to the interfacial composition, the conventional way to derive the solution of v int is by calculating the values of the three energies ( G chem m , G diff m and G friction m ) as a function of v int in the expected range. The intersection point between G chem m and G diff m + G friction m is regarded as the solution of v int (when multiple intersection points are present, the minimum v int is adopted). Although this method is robust, it is computationally heavy as it needs to be calculated via numerical procedures for every moving interface while taking into account the local (thermal history dependent) chemical composition. In this work, we have derived an analytical solution for G diff m , which avoids the need to compute the Gibbs free energies as a function of v int and allows one to directly compute the solution. Thereby, the computation can be significantly accelerated enabling the model to be used in a multigrain setting. A detailed description on how to calculate G diff m is presented in Appendix A . The interfacial carbon concentrations and diffusion profiles need to be solved to derive v int . As the diffusivity of carbon in ferrite D α C is much larger than that in austenite D γ C , we only consider the bulk diffusion of carbon in austenite. This way, there is no carbon concentration gradient in ferrite and thereby the interfacial carbon concentration C α amounts to the equilibrium concentration C αγ eq . The diffusion profile of carbon in austenite is approximated by a second-order polynomial. This approximation yields mathematical simplicity and was found to give a good match with the limiting diffusional cases for which exact solutions exist. Details of the mathematical treatment of the carbon diffusion are presented in Appendix B . For non-overlapping diffusion fields the rule of mass conservation of C results in: where C 0 is the nominal carbon concentration, C αγ eq is the equilibrium concentration of carbon in ferrite, R α is the radius and V α is the volume of the ferrite grain, L is the diffusion length and r is the distance, C γ is the interfacial carbon concentration at the austenite side. As there is no carbon accumulated at the interface itself, the following equation can be derived: Now we can derive the solutions for the unknown parameters v int , C γ and the diffusion length L by solving the set of Eqs. 2 -6 . In Section 2.3 we present how we solve this problem. For the overlapping diffusion fields (soft-impingement) the carbon concentration at the soft impingement point C γ m increases. The mass conservation of carbon and the absence of accumulation of carbon at the interface still applies in this stage. Therefore, we have: where L 0 is the distance where the carbon diffusion field of a ferrite grain starts to contact with the diffusion fields of its neighbors. Now the solutions for v int , C γ , C m and L can be derived by solving the set of Eqs. 2 -4 and Eqs. 7 -9.
To distinguish whether soft impingement occurs, we first treat the diffusion of carbon without soft impingement and derive the diffusion length L. When the profile of grain j starts to intersect with that of any neighboring grain k, the following condition is fulfilled: where r jk is the distance between the centers of ferrite grains j and k. At this moment the diffusion length L is taken as the maximum distance that carbon can diffuse for that grain. With the advance of the interface migration, the ferrite grains could be in physical contact with neighboring grains, resulting in hard impingement. To account for this hard impingement effect, we use the same correction as the one applied in our previous work [43] , for which the source code now can be found at ( https://github.com/haixingfang/ 3D-mixed-mode-model ).

Simulation conditions
The model is used to compute transformation kinetics under various thermal conditions, including continuous cooling, isothermal holding and thermal cycling for a series of Fe-C-Mn alloys. This model is first benchmarked for continuous cooling at different rates and isothermal holding at a selected temperature for the Fe-0.1C-0.5Mn alloy (see temperature profile in Fig. 1 a and 1 b). Then we focus on using this model for cyclic partial phase transformations in the γ / α two-phase region between temperatures T 1 and T 2 (see Fig. 1 c). Table 1 gives these temperatures and the corresponding equilibrium ferrite fractions calculated with Thermo-Calc software using the TCFE8 database. The two cycling temperatures were chosen to be the same as in the reported experimental studies (for Fe-0.023C-0.17Mn [11] and Fe-0.1C-0.5Mn alloys [38] ) or chosen to ensure that the ferrite fraction is not too small and also not too large at the end of isothermal holding (for the other alloys examined here). For those alloys the cycling temperatures were chosen to show ortho-equilibrium ferrite fractions of f OE α (T 1 ) ≈ 0.70 and f OE α (T 2 ) ≈ 0.20. For alloys with 1.5 wt.% Mn or more the T 1 was set as f α (T 1 ) ≈ 0.50 for the end of the isothermal stage and T 2 fulfills T 2 -T 1 ≈ 50 K. The reason for setting up a different criterion for higher Mn alloys is that a substantial ferrite transformation only occurs well below the transition temperature between Fig. 1. Schematic temperature profiles for (a) continuous cooling, (b) isothermal and (c) cyclic partial phase transformations. All simulations start at A 3 (PE) and end at A 1e for continuous cooling and cyclic transformations, while the simulation does not end until there is negligible change of f α for the isothermal transformation. For simulations of the cyclic phase transformation, the primary cooling rate is set as 20 K/s, the isothermal time at T 1 is set as 30 min and the number of cycles is chosen to be 3.

Table 1
Compositions of the Fe-C-Mn alloys and the cycling temperatures T 1 and T 2 for each simulation. The ferrite fractions under ortho-equilibrium f OE α , the A 3 and A 1 temperatures under ortho-equilibrium (A 3e and A 1e , respectively), para-equilibrium temperature A 3 (PE) and PLE/NPLE temperature are also listed. partitioning local equilibrium and negligible portioning local equilibrium (PLE/NPLE) according to preliminary simulation results. The PLE/NPLE temperatures are also listed in Table 1 .
The length of the cubic sample box is set as L b = 70 μm. The average austenite grain size is chosen to be d γ = 20 μm ( ρ γ = 2.4 × 10 14 m −3 ) with d min = 12 μm. Given these dimensions and two imposed cooling rates the model describes the collective and individual behavior of about 151 parent austenite grains and about 58-204 emerging ferrite grains depending on the cooling rate. As the chemical potentials of all elements depend on both temperature and compositions, they must be calculated separately for each grain at each time step. To accelerate the computation, we use multiple nonlinear regression to approximate the dependency of the chemical potentials of Fe, C and Mn beforehand to avoid calls to the Thermo-Calc database. A detailed description is presented in Appendix B . Using this approximation, combined with the analytical expression derived for G di f f m (shown in Eq. A4 in Appendix C ), the calculation of the transformation kinetics in 3D for a multigrain setting can be significantly accelerated and be realized in realistic computing times.
The phase boundary lines of ( α+ γ )/ γ and α/( α+ γ ) are calculated with Thermo-Calc and fitted to a second-order polynomial for the temperature range of interest. The parameters for intrinsic interface mobility proposed by Zhu and coworkers [31] are used in the present work. The diffusivity of Mn in the interface is taken Except for the thermodynamic data, the values of all other parameters are kept constant independent of the alloy composition and thermal treatment. The modelling parameters are listed in Table 2 . The binding energy E 0 for Mn on the α/ γ interface is a very important parameter as it governs the Mn partitioning inside the interface, thereby strongly influencing the transformation kinetics. Slightly different values for E 0 have been derived using different approaches: scanning Auger microprobe studies on austenite grain boundaries yielded E 0 = 8 ± 3 kJmol −1 [58] ; fittings to experimental transformation kinetics yielded E 0 = 5 9 .9 kJmol −1 [37] ; calculations from the Mn profile across the in-terface measured by atom probe tomography gave E 0 = 6.0 ± 1.4 kJmol -1 [32] and first principles calculations provided a value of E 0 ≈ 12.5 kJmol −1 [34] . In the present work we have chosen an intermediate value of E 0 = 7 kJmol −1 . The choice of a slightly different value for E 0 would have affected to overall kinetics, but would not have affected the key findings in the simulations.
In this model the only adjustable parameter is the pre-factor A in Eq. 1 . This pre-factor A affects the maximum number of nuclei and the effective time (temperature) window for nucleation. The parameter depends on the rate of transformation, which depends on the imposed temperature path. It should be noted that any small uncertainty in the modelled energy barrier for nucleation (controlled by the parameter ) will directly translate into a significant change in the value of the required pre-factor A. In the current simulations for linear cooling, the value for the prefactor A is estimated from the maximum number density of ferrite grains, which was estimated from metallographic observations. For isothermal holding and thermal cycling simulations, the pre-factor A is set in the range of 0.001 ~0.0 033 (i.e. 1/10 0 0 ~1/30 0). This choice was based on reported values in synchrotron X-ray diffraction studies on isothermal ferrite-to-austenite transformations [46] .
For each time step t during the calculation of the phase transformation the number of ferrite grains is calculated according to the CNT (see equation 1 ). The interface velocity for each grain is then calculated as described in Section 2.2 . The effective radius R α, j (after corrections of hard impingement if necessary) of ferrite grain j at a time t is calculated using the velocity derived from the The overall microstructural characteristics, including the ferrite volume fraction f α , the average grain radius δ α and the standard deviation σ α for the radius distribution of ferrite grains, are computed for each time step. The programming codes for implementing the model presented here are publicly available ( https://github. com/haixingfang/3D-GEB-mixed-mode-model ).   [19] . Lines in (d) are fitting curves of lognormal distribution. Fig. 2 shows the ferrite transformation kinetics in the Fe-0.1C-0.5Mn alloy during continuous cooling at a rate of 0.4 and 10 K/s. In the simulation, the pre-factor A in Eq. 1 was set at A = 1/30 0 0 and at A = 1/35 for 0.4 and 10 K/s, respectively. The pre-factor is the only parameter required to be adjusted and it is chosen such to make sure that the maximum value for the ferrite number density was comparable to that of metallographic observations on the fully transformed sample [62] . As a result, the maximum ferrite number density for 10 K/s is about 3.5 times higher than that for 0.4 K/s, even though the nucleation starts at about the same temperature of 1095 K, as shown in Fig. 2 a. Fig. 2 b shows that the model predicted ferrite fractions (f α ) that are in good agreement with the dilatometer results for the two cooling rates. Fig. 2 c plots the evolution of the calculated average ferrite grain size ( δ α ) and standard deviation of the ferrite grain size distribution ( σ α ) as a function of temperature. The final values for δ α (10.7 and 6.9 μm for 0.4 and 10 K/s, respectively) are as intended close to the experimental data determined on metallographic images (10.7 and 6.5 μm, respectively) [61] . Fig. 2 d shows a comparison of the ferrite grain size distribution after a complete transformation between the model prediction and the experimental results for 10 K/s. The original experimental data for the ferrite grain size were presented as the equivalent circular diameter in 2D [19] . To make a comparison with our 3D data, we converted the published 2D data into 3D by R α, 3D = 4R α, 2D / π for each ferrite grain. The figure shows a reasonable agreement between the modelling results and the experimental data, both of which can be approximated by a lognormal distribution. It can also be seen in Fig. 2 d that the model predicts a narrower size distribution and a slightly larger average grain size, which could partly be attributed to an earlier transformation onset for the modelling (see Fig. 2 b). It is also worth pointing out that the actual grain size and shape distribution of the fully austenitic starting structure influences the final ferrite grain size distribution. Given a relatively uniform and isotropic starting microstructure (as presented here), the final grain structure is mainly controlled by the average grain size of the austenite. The model itself can easily be adapted to include non-equiaxed starting microstructures.

Transformations during continuous cooling
We will now try and demonstrate the capabilities of the model in predicting the transformation kinetics both at the overall sample level and at the level of individual emerging ferrite grains. Fig. 3 shows the development of two selected ferrite grains ( α A and α B ) nucleated at the same position but under different cooling rates of 0.4 and 10 K/s, respectively. Snapshots of Fig. 3 a- But once it starts at 1081 K, it has a very high initial velocity due to a large undercooling, with a larger driving force compared to α A at the starting point. Different from α A , the interface velocities of α B continuously decrease to zero with a small bump occurring at 1072 K. For cooling at a rate of 10 K/s, α A and α B show very similar behavior in terms of the starting growth temperature, the energy changes and the interfacial velocities during the whole process, which is very different from the behavior observed at a cooling rate of 0.4 K/s. The average and local concentration profiles of C and Mn also evolve distinctly differently for the two grains at different cooling rates. It can be seen from Fig. 4 b that Mn concentration at the interface builds up to a high level, especially when f α = 0.50, while the C diffusion gradient is small and even smeared out at f α = 0.50 for α A upon cooling at 0.4 K/s (see Fig. 4 a). This indicates that trans-diffusion of Mn in the interface dominates the kinetics. For the same grain nucleated at the same temperature and the same position but exposed to a higher cooling rate, the Mn concentration across the interface is much lower and the C diffusion gradient is much higher (see Fig. 4 c and 4 d), indicating that in this case (a cooling rate of 10 K/s) the transformation kinetics is mainly controlled by C diffusion.
The modelling results of the transformation kinetics during continuous cooling demonstrate that: 1) the model performs well in predicting the overall transformation kinetics during continuous cooling; 2) the development of physical and chemical characteristics can be monitored for each individual grain and 3) the contribution and dissipations of the Gibbs free energies can be calculated to understand the mechanism for the change in interfacial velocity.

Isothermal transformations
A further modelling test was made on the same alloy of Fe-0.1C-0.5Mn but imposed to isothermal phase transformations at 1058 and 1048 K to enable model validation against phase field simulations and dilatometry measurements, respectively. In both simulations the pre-factor in Eq. 1 was set at A = 1/800. Fig. 5 a shows the predicted f α as a function of isothermal time together with the result obtained from 1D phase field simulation. In both models, Mn is modelled such that trans-diffusion in the interface, i.e. solute drag, can take place. The figure shows that in both approaches f α evolves in a similar manner and ends at about the same level of 0.65, which is smaller than the para-equilibrium fraction f PE α (T = 1058 K) = 0.72. As has been reported in literature [e.g. 12,13,37,63 ], this is a case of incomplete transformation because the energy dissipation of solute drag becomes larger than the driving force at longer isothermal times. Notably, there is an obvious latent transformation in the early stage in the current modelling results (marked by the red arrows in Fig. 5 a), whilst this effect has not been observed for the 1D phase field modelling [39] . This is due to differences in nucleation modes and grain geometries: the present model implements a continuous nucleation mode Fig. 4. Carbon (a, c) and manganese (b, d) profiles across the interface for grain α A (as visualized in Fig. 3 ) at cooling rates of (a, b) 0.4 K/s and (c, d) 10 K/s. r / δ is a dimensionless quantity where r is distance and δ is the half of the interface thickness.  [39] considering the Gibbs free energy dissipation due to Mn diffusion inside the interface and (b) at 775 °C against experimentally derived f α from dilatometer measurements with a primary cooling rate of 1 K/s. The time is expressed relative to the starting moment of the isothermal holding. The experimental f α in (b) is derived from average of two separate dilatometry measurements (each starting with a fresh sample). and operates in 3D, leading to a scaling of f α with R α 3 , whereas an instantaneous nucleation and a planar geometry are assumed in the phase field and therefore, f α scales with R α . Differences in the median isothermal stage, where f α = 0.25 ~0.55, are also seen, which could be due to the fact that our model takes into account both nucleation and spatial correlations between neighboring ferrite grains.
Although the differences in the overall f α between the present 3D multigrain model and the 1D phase field model predictions presented in Fig. 5 a are not large, the 3D multigrain model shows a significant variation in transformation behavior for individual ferrite grains due to grain interactions (related to nucleation time and local environment), whereas the 1D model by default cannot give such information. Other experiments and modelling studies [45,64] also demonstrated that the transformation behavior of individual grains can vary significantly. There are a number of significant differences between 1D and 3D models. As to the 3D model, the transformation behavior for each grain is determined by the nucleation and growth behavior of the surrounding material, in addition to what has been prescribed for that particular grain. Such behavior can never be captured a priori in a 1D model. Only when one knows the full 3D behavior then the settings of the 1D simulation can be set properly. The dimensionality will also affect the way that soft and hard impingement play a role in the transfor- mation kinetics. The de-activation of nucleation sites by soft/hard impingement and the evolution of the grain size distribution can only be obtained in a 3D model. Given their geometrical simplicity, 1D models are more suitable to demonstrate the effect of various assumptions regarding solute behavior near the moving interface on the transformation kinetics. Fig. 5 b compares the predicted f α with the experimentally determined f α from dilatometry measurements by applying the lever rule [65] . By setting the average austenite grain size to < d γ > = 50 μm in accordance with experiments [66] and the sample box L b = 175 μm (using the same ratio of L b to < d γ > ), while keeping the pre-factor A the same as in the simulation shown in Fig. 5 a, we found that the predicted f α matches the experimental result closely as shown in Fig. 5 b. The model and the experiment show a latent transformation in the early stage and a similar value of f α at the final stage that is close to the para-equilibrium fraction f PE α (T = 1048 K) = 0.76). The good agreement between the simulations (using only one freely adjustable parameter and 11 independently determined parameters) and the dilatometer measurements validates the good transferability of the current model.
It is worth to note that the present model is aimed to describe bulk behavior and therefore assumes periodic boundary conditions. In the current format, the model is not appropriate to describe and predict the transformation kinetics for isothermal transformation experiments that are based on the removal of carbon via gas/sample reactions at the surface (including both mass loss and macroscopic carbon diffusion towards the outer surface of the material), e.g. the decarburization experiments as conducted in [14] .

Cyclic partial phase transformations
When analyzing a cyclic partial phase transformation the austenite-to-ferrite transformation is classified as leading to a positive interfacial migration (v int > 0) while the ferrite-to-austenite transformation leads to a negative interfacial migration (v int < 0). Similarly, the driving force for the austenite-to-ferrite transformation is connected to positive changes in the Gibbs free energy while the ferrite-to-austenite transformation is connected to negative changes in the Gibb free energy.

Comparison with dilatometer experiments
We first verify our simulation results against dilatometer experiments for two different alloys during cycling at a rate of 10 K/min. Fig. 6 shows that our model can adequately reproduce the experimental cyclic behavior of f α for both Fe-0.023C-0.17Mn and Fe-0.1C-0.5Mn alloys. The stagnant and inverse transformation stages, which are two distinct features of the cyclic behavior in low Mn steels, are captured in the simulations. Whilst it can be seen that our simulations predict a longer stagnant stage (especially during cooling) and a shorter inverse transformation stage, the simulations predict quite a good range of f α where the transformation could span during each cycle and how f α evolves in each cycle compared to the experimental data.
Interestingly, in both alloys f α starts at a higher level at the beginning of the first cycle, whilst f α never returns completely to its original level in the following cycles, which themselves are perfectly reproducible. This suggests that the transformation at the end of isothermal holding has reached an equilibrium stage, which is not reached again during cycling due to the prevailing nonequilibrium interfacial conditions caused by dynamical changes in temperature. It should be noted that these two transformation features not always occur. In our previous experimental study on a steel with a higher Mn concentration (Fe-0.25C-2.1Mn alloy [35] ), f α was found to be far from equilibrium at the end of isothermal holding. Consequently, f α progressively increased after each thermal cycle, showing an opposite behavior from what is seen in Fig. 6 . Fig. 7 shows the behavior of the first nucleated grain during cycling from the simulations of Fe-0.1C-0.5Mn alloy. Similar to f α , R α also experiences a loop consisting of two stagnant, two normal and two inverse stages during each cycle shown in Fig. 7 a. Each stage is characterized by different f eatures of the interfacial velocities as plotted in Fig. 7 b and the features can be understood from the difference between the chemical driving force G chem m and the energy dissipation of solute drag G di f f m (see Fig. 7 c-e). To assist the analysis, we marked points of interest by A-G in each graph. The character of the transformation stage marked by these points is summarized as follows: 2) Normal stage where α → γ proceeds during heating (from point C to D). An increase in temperature above point C makes the driving force sufficiently larger than the solute drag force, To investigate the evolution of the cyclic behavior over successive cycles, we plot the radius and the interfacial velocities as a function of temperature in Fig. 8 for six ferrite grains (including the one shown in Fig. 7 ). Five of the grains behave similarly in terms of both R α (forming a loop) and v int , whereas one grain that nucleated the latest shrinks during heating after a stagnant period and finally disappears, showing an abnormal behavior. This abnormal behavior is also found for some other grains that nucleate late and are relatively small before cycling begins. The reason is that the overlapping of carbon concentration fields from their neighboring larger grains, together with the temperature increase during heating, decreases the chemical driving force for ferrite formation to stop the continuous shrinkage. In principle, the curvature is also expected to play a role here, but this effect is not considered in the current modelling.

Cyclic behavior of individual grains
For the five grains showing a normal transformation behavior, their radii span different ranges and their velocity peaks at different temperatures, showing a significant variation in behavior, as has also been observed in in-situ EBSD and TEM experiments [67,68] . Two of the five grains ( t nuc = 21.3 and 106.1 s in Fig. 8 a) show a considerably longer inverse transformation stage than the rest. For one of them ( t nuc = 106.1 s) the interface does not go back to the original position, which coincides with the observation that f α does not recover to the starting value upon the completion of the first cycle (see Fig. 6 b). The reason is that the ferrite grains which nucleated late in the cycle have higher remote carbon concentration in their parent austenite grains because of impingement of the carbon diffusion fields from earlier nucleated ferrite grains.
As a result, α → γ is more favored than γ → α for these late nucleated grains. It can be seen in Fig. 8 b that the v int for the grain with t nuc = 106.1 s becomes negative earlier (i.e. α → γ starts earlier) and positive later (i.e. γ → α starts later) than other earlier nucleated ferrite grains. For the grain with t nuc = 113.4 s it even disappears during the first cycle.
Since these 6 grains selected include the earliest nucleated one, the latest nucleated one and the some nucleated in between, together they are expected to represent the spectrum of cyclic behavior for all grains. It can be seen from Fig. 8 b that the interface velocity always shows a peak and its behavior agrees well with the recent in-situ TEM observation on two single γ / α interfaces [68] .
The peaks of interface velocities for the normal behaved grains are located in the range 803 -808 °C for α → γ and ~840 °C for γ → α, which is in close agreement with the TEM results.
It should be pointed out that a full grain shrinkage during cycling cannot be predicted in any previous models for cyclic transformations because they do not consider ferrite nucleation before the cycling starts. However, in-situ neutron depolarization experiments showed that the actual number density of ferrite grains during cycling over certain temperature windows can decrease slightly, due to the occurrence of coarsening (resulting in the disappearance of some small grains) [35] . To unambiguously track the behavior of individual grains during cyclic transformations, grain resolved 3D/4D techniques based on synchrotron high-energy Xray diffraction using micro beam [69,70] seems to be the most promising approach.

Cyclic behavior as a function of C and Mn concentrations
The model can also be extended to predict the cyclic transformation kinetics as a function of wide range concentrations for C and Mn, which will allow us to elucidate their interplay. Fig. 9 shows the evolution of f α for a series of Fe-C-Mn alloys with temperature in comparison with the equilibrium fraction. Different behaviors for f α are clearly observed. As the Mn concentration increases from 0.17 to 1.0 wt.% for Fe-0.1C-Mn alloys ( Fig. 9 a-f), an open loop for f α can clearly be observed, while the stagnant stage (measured by the temperature range where interface is immobile) becomes longer. It is interesting to note that when the Mn concentration increases to 1.0 wt.%, f α spans a much smaller range compared to the lower-Mn alloys. With further increasing Mn concentrations, the amplitude of the changes in f α in each cycle is so small that the open loop is hardly seen. This means the interfacial Mn partitioning largely controls the cyclic kinetics and limits any substantial interfacial movements in these alloys. However, carbon partitioning also plays a role here. Even for a relatively high Mn alloy, when C concentration decreases, the transformation kinetics of both γ → α and α → γ could be accelerated. This acceleration is observable when comparing Fig. 9 e and 9 g. Conversely, when the C concentration is increased in alloys with a fixed Mn concentration, then the kinetics is gradually suppressed and f α cycles over a smaller range. This is well illustrated when Fig. 9 a is compared with Fig. 9 h. Similar behavior is seen when Fig. 9 e is compared with Fig. 9 i.
Previous studies have focused on the effect of a single alloying element (either C or Mn) on the cyclic transformation behavior [66,71,72] . To characterize the stagnant behavior during cycling, the length of the stagnant stage (the temperature range where f α does not change) during cycling or the length of the post cyclic stagnant stage (the temperature range where f α does not change after the thermal cycling is finished, which is due to residual Mn spikes formed during cycling in front of the interfaces) are used. It is shown that increasing the C or Mn concentration leads to more stagnant movement of interfaces during cycling [71,72] . The post cyclic stagnant stage, present right after the cycling finishes, became evident ( > 2.5 K) when the Mn concentration increased above 1.0 wt.% Mn [66] . It can be seen from Fig. 9 that C and Mn have a coupled effect on the cyclic transformation behavior. To quantify this effect, it is not very suitable to use the length of stagnant stage as such, as has been done in previous studies, because the equilibrium states at the beginning of the cycling and the interval of equilibrium fractions between T 1 and T 2 are different. Here, we integrate f α over the temperature normalized by integration of the net para-equilibrium faction f α PE surplus from that of T 2 over the same temperature range to describe the stagnant transformation behavior (denoted as η stag ): A high value of η stag corresponds to a more stagnant behavior during the cyclic phase transformation. Intuitively, the numerator in Eq. 12 represents the area of the cyclic loop in the f α -T plot and the denominator represents the area of f α over T as if it was formed by quenching an alloy from T 2 to T 1 , then holding at T 1 for a long time and finally heating back to T 2 under para-equilibrium conditions. Since the transformation behavior becomes repeatable from the second cycle on, we used the data for the second cycle to calculate η stag for all alloys. The metric η stag can be linked to the evolution of the C and Mn concentration profiles during cyclic phase transformation. In Fig. 11 we plot the radius of the first nucleated ferrite grains and their C and Mn concentration profiles across the interface at five selected positions in the first cycle for the Fe-0.1C-0.1Mn alloy with a relatively low η stag and for the Fe-0.1C-1Mn alloy with a relatively high η stag (as marked in Fig. 10 ), respectively. Although the two grains both show characteristic cyclic behavior (stagnant, inverse and normal transformations) as described in Section 3.3.2 , the amplitudes of the change in R α and the resulting loop outlined by R α are significantly different (see Fig. 11 a and 11 d). Plots of the C profiles at the five selected positions indicate that the evolution of the C profile is much more pronounced for the Fe-0.1C-0.17Mn alloy than for the Fe-0.1C-1Mn alloy (see Fig. 11 b and 11 e). However, in both alloys the C profile in austenite has a very small concentration gradient, suggesting that C diffusion is not the rate limiting step for both alloys.   consistent with the observation that η stag for the Fe-0.1C-1Mn alloy is larger than that for the Fe-0.1C-0.17Mn alloy.

Conclusions
A novel 3D mixed-model has been developed by considering ferrite nucleation and growth governed by a balance in the Gibbs free energy contributions between chemical driving force and dissipation by interface friction and by solute drag for austeniteferrite transformations in Fe-C-Mn alloys. A new implementation of the solute drag dissipation calculations is proposed and shown to be computationally very efficient. The capabilities of this model have been illustrated. Based on the model application to austeniteferrite transformations in Fe-C-Mn alloys, the following conclusions can be drawn: 1) The derived analytical solution for calculating the energy dissipation due to solute drag is computationally much more efficient than the conventional approach that has to compute the Gibbs free energies as a function of interface velocity in a predefined range. This new method proposed here is general and can be used in other 3D simulations. 2) The model performs well in predicting the overall transformation kinetics during continuous cooling, isothermal holding and inter-critical thermal cycling. The grain size distribution can be predicted reasonably well for continuous cooling.
3) The transformation behavior of each individual grain can be monitored. A wide variation in behavior was found. The local Gibbs free energy of the chemical driving force and the dissipations can be calculated to understand the mechanism for the change in interfacial velocity. The carbon partitioning and the solute drag effect of Mn complement each other depending on the specific conditions. 4) Small ferrite grains that nucleated relatively late may fully shrink upon heating during cyclic partial phase transformations. This phenomenon has not been predicted by models that do not consider ferrite nucleation before the beginning of thermal cycling.

5)
A new metric is proposed to describe the stagnant interfacial behavior during cyclic phase transformations. It is shown that the coupled effect of C and Mn on the cyclic transformation behavior can adequately be described. A transition region as a function of C and Mn concentrations is identified to distinguish different types of cyclic behavior.
Although we only demonstrate the application of the model for Fe-C-Mn alloys, the model can also be used for calculating the transformation kinetics for other ternary alloys once the corresponding thermodynamic information is available. Since a more efficient implementation of computing solute drag dissipation is proposed, the present work is expected to encourage more advanced 3D simulations using the principle of the Gibbs free energy balance. To this aim the code is made publicly available at https://github.com/haixingfang/3D-GEB-mixed-mode-model .

Declaration of Competing Interest
We declare no conflict of interest. considered as a potential well. The solute segregation thereby imposes a drag effect on the interface migration. Assuming a triangular potential well within the interface, Purdy and Brechét [57] proposed that the equilibrium potential of solute μ(r) over a distance r depends on the temperature T, solute concentration profile C(r) and the interaction energy E(r): where R is the gas constant. The interaction energy E(r) is given by The substitutional concentration profile for an interface moving with a quasi-steady velocity v int must fulfil: where D int is the interfacial diffusivity of the substitutional solute. The concentration profile can be derived by solving Eq. A2. The solution at four different regimes can be expressed by the following equations containing four dimensionless quantities a, b, V and S: where C 0 is the nominal concentration of the substitutional solute, The concentration profile is thus asymmetric, as shown in Fig. A1b. Now the energy dissipation G di f f m due to trans-diffusion of substitutional alloying element inside the interface can be calculated using Cahn's equation [73] : where β 1 and β 2 are two parameters: β 1 = a 2 RT v int C 0 δ DV ( 1+2 a + a 2 ) and β 2 = bRT v int C 0 δ DV ( 1+ a +2 b+2 ab + b 2 + b 2 a ) . The solution expressed in Eq. A4 can considerably accelerate the computation by avoiding the integration in Eq. A3. Like G f riction m , G di f f m also depends on v int , which makes the solutions for the two energy dissipation terms interdependent and thereby requires the solution method described in Section 2.2 .
Eq. A4 also applies for v int < 0 , i.e. the ferrite to austenite phase transformation. For v int < 0 the dissipation energy due to transdiffusion of solute inside the interface G di f f m has a negative value. The chemical driving force G chem m and the dissipation energy due to interface friction G f riction m when v int < 0 are obtained in the same way as for v int > 0 . Therefore, a change in the transformation direction is automatically accounted for when Eq. 3 is used.

Appendix B: Carbon diffusion profiles
For an early stage of ferrite growth, where the carbon diffusion fields surrounding the ferrite grains do not overlap (see Fig. B1 a), the carbon concentration far away from the γ / α interface, C ∞ , equals the nominal concentration: C ∞ = C 0 . The carbon concentration profile surrounding the ferrite grain, C(r), can be approximated by a second-order polynomial: where r is the distance from the interface (r = 0 at the γ / α interface), C 0 , C γ and L have the same meaning as described in Section 2.2 . This concentration profile fulfills the following boundary conditions: Combining Eqs. B1-5 yields the expression of Eq. 5 as given in Section 2.2 .
As the austenite-to-ferrite transformation proceeds, the carbon diffusion field surrounding a growing ferrite grain may start to overlap with that of a neighboring one (see Fig. B1 b). When this occurs, the carbon concentration profile can be written as: where C m is the carbon concentration at the soft impingement point m. Eq. B6 fulfills the following conditions: Similar to the non-overlapping stage, the mass of carbon must be conserved. Replacing C(r) in Eq. B5 by the expression of Eq. B6 results in Eq. 7 as presented in Section 2.2 .
Assuming that the carbon can diffuse fast enough outside of the soft impingement region, the carbon concentration in the matrix of the austenite grain becomes homogenous: C ∞ = C m . This means that the soft impingement causes an increase in C ∞ , which slows down the growth of any other ferrite grains nucleated at other sites of the same austenite grain even if they have not shown soft impingement with any other grains (e.g. grain α 3 in Fig. B1 b).

Appendix C: Approximations for the chemical potentials
As the transformation proceeds, the interfacial concentrations of the ferrite grains are changing continuously. Thus, the chemical potentials are changing with the transformation time for all of the grains. To calculate the driving force for individual grains, it is necessary to calculate the chemical potentials as a function of temperature and composition. For the Fe-C-Mn system, this dependence for austenite or ferrite can be described as follows: μ C = μ 0 , C + RT ln ( x C ) + RT ( e CC x C + e CMn x Mn where R is the gas constant, T is temperature (in K), e CC , e CMn , e MnC and e M nM n are interaction coefficients. As μ F e mainly depends on temperature, with a weak dependence on the composition, we rewrite Eq. C3 to μ F e = μ 0 , F e + RT ln ( where ω is a constant. We employed Thermo-Calc with the TCFE8 database and calculated the chemical potentials, which were fitted to Eqs. C1-2 and Eq. C4. The obtained fitting parameters for the chemical potentials in γ are given by: The relative errors are about 11.2%, 3.3% and 0.05% for Eqs. C8-10, respectively.