Three-dIMenSIonal phaSe fIeld SIMulaTIon for rafTIng of MulTIparTIcle precIpITaTe In elaSTIc InhoMogeneouS alloy under exTernal STreSS

Nickel-based super alloys are the main candidate materials for aero-engines, gas turbine blades etc. This paper focuses on the simulation of nucleation and growth kinetics of γ' phase, and stress response mechanism of γ’ phase particles during their preferential coarsening (rafting) in elastic inhomogeneous system. A phase-field model is employed in the present study, which incorporates chemical, interfacial, and elastic energies, and it couples essentially to externally imposed mechanical field. Due to the limitations of the 2D model on analyzing the shape and size of the precipitate particles, the process of γ' phase particles growing and coarsening is further modeled by performing 3D simulation. The results show that the average particle size is linearly related to the evolution time and satisfies the Lifshitz-Slyozov-Wagner (LSW) classical coarsening theory when the external stress is not applied. Particles exhibit a strong special orientation under tensile stress, and the orientation is in excellent agreement with previous studies. In the nucleation stage, the collision and coalescence between particles promote rafting significantly, and the number of soft particles is obviously larger than that of hard particles. In the coarsening stage, the growth rate of soft particles is higher than that of hard particles. Three-dimensional simulation results show that the effect of final characteristic size of precipitated particles is not significant by external loads. The morphology evolution and coarsening mechanism of the precipitated particles are of great significance for studying the strengthening mechanism of super-alloy.


Introduction
Nickel based super-alloys, as one of the most important high-temperature structural materials, exhibit unique performance at elevated temperatures. The nickel based super alloys are typically a mixture of disordered FCC phase (γ) and the ordered L1 2 phase (γ') [1]. Precipitation hardening by the ordered phase is the major strengthening mechanism for Nibase super alloys [2]. Therefore, the microstructure control of γ' and γ phase is the key technique in the development of advanced Ni-base super alloys.
Rafting (which is a preferential coarsening) in nickel alloy is essentially a phenomenon about morphological instabilities, which leads to the destruction of an initially periodical arrangement of precipitates during service and eventually causes the formation of regular coarsening structure. Rafting can bring about either hardening or softening of the microstructure [3][4][5][6][7]. The kinetics of directional coarsening is affected by the type of the applied stress (tensile or compression condition), the sign and magnitude of the γ / γ' lattice misfit, and the difference in the elastic constants of the two phases [7,9]. More generally, the application of a macroscopic external load triggers an anisotropic misfit relaxation, which is ascribed to the dislocations generated during creep test [10].
The mechanical properties (e.g., strength, fatigue, and creep) are strongly dependent on the volume fraction, size, microstructure and morphological change of the γ' phase. Therefore, extensive efforts have been devoted to understanding the γ' precipitate evolution experimentally. Tien and Gamble [11] found that γ' phase particles distributed uniformly in the matrix phase along the elastic soft direction after the standard solution and aging treatment. Due to its significant advantages in the characterization of elemental composition, the three-dimensional atom probe (3D-AP) is employed to investigate the microstructure and the evolution of the composition distribution of a series of Ni-based super alloys [12,13] to obtain a better understanding about the coarsening behavior of the precipitated phase. Murakumo et al [14] found that the creep fracture life is strongly dependent on the volume fraction of γ' phase. A close examination of the experimental results has shown that strain accumulation in different regimes depends on the value of the applied stress [15]. 3D-AP provides an effective way for the fine regulation of the type, quantity, morphology, size, and orientation of precipitated strengthening phases. However, the factors that influence the microstructure formation are extensive, which increases the difficulty and complexity of controlling the microstructure. Many trial-and-error experiments are performed for searching the best combination of desired microstructures and material properties. One way of overcoming these shortcomings is to model the microstructural evolution by the simulation method, which captures the image of the microstructure evolution of the phase transition process.
Many attempts have been made to simulate the microstructure evolution and coarsening kinetics of precipitated phases, especially by means of the phase field method (PFM). PFM is a well-developed and powerful tool, in which a set of conserved or nonconserved physical and artificial fields are employed to describe the microstructure of crystalline. And the evolution of microstructure is predicted by solving a set of mathematical equations which governs the field's evolution. It has been demonstrated in quite a lot of areas, such as grain growth, dendrite growth, crack propagation, spinodal decomposition, and so on [16][17][18]. For example, Wang and Khachaturyan investigated the morphology of a single coherent precipitate during the growth within two dimensions (2D) [19]. Based on the inhomogeneous elastic theory, Li and Chen [20] studied the rafting behavior of γ' phase under external loading. The results showed that the effect of the direction of the applied force field on the direction of precipitation rafting is very important, which is consistent with the results of precipitation rafting under tensile and compression conditions. In addition, the direction of precipitation rafting under external loading was studied by Wang et al [21] and Yang et al [22], and the multi-particle effects, the elastic and diffusional interactions among precipitates, were studied [23]. The experimentally observed morphologic patterns, such as the change of precipitate shape and its alignment along the elastic soft direction of the matrix, have been successfully predicted by 2D phase field simulation. The effects of precipitate volume fraction on the coarsening kinetics are also modeled by Vaithyanathan and Chen [24].
Despite the fact that the phase field simulation has successfully predicted the evolution of various forms of precipitated particles and provided some basic thermodynamic and kinetic mechanisms, most of the present phase field simulations are qualitative. In order to describe the quantitative information of the characteristic size of the precipitated phase in elastic inhomogeneous alloy systems under stress, the elastic difference between the precipitated particles and the matrix phase was defined by model parameters. An elastic contribution is added to the existing phase field model to describe the kinetic process of precipitated particles. Moreover, the precipitation phase in elastic inhomogeneous alloy is mostly carried out in 2D phase field simulation, and the three-dimensional (3D) simulation of precipitated phase is seldom involved. The 2D information to estimate particle shape and particle volume fraction is not rigorous enough, resulting in errors in the research results. Therefore, the quantitative information of the volume fraction of the precipitation phase under stress is obtained by three-dimensional simulation, and the internal mechanism of the precipitation phase rafting is explained by comparing the results of twodimensional simulation and three-dimensional simulation. A comparison with experimentally observed coarsening rate would require 3D simulations due to solute precipitates in all the directions.
In this paper, this PFM approach is applied to the specific case of γ / γ' system. First, 2D simulations are performed to investigate the differences in microstructure evolution with homogeneous and inhomogeneous elasticity, and to analyze directional coarsening on a sufficiently large scale. The classical theory is verified for predicting rafting from the view of two-dimensional simulation, and the variation of nucleation particle number with respect to evolution time is analyzed. Secondly, we present 3D simulations of γ / γ' microstructural evolution under external load in order to investigate the influence of the elastic driving force on the directional coarsening. The encounter and collision process of precipitated particles in bulk materials are described in detail, and the alloy strength with the evolution time is preliminary revealed by the change of the average size of precipitated particles. We will analyze the obtained microstructure under external load in terms of a meanfield approach of inhomogeneous elasticity in order to better understand the underlying physical processes.
We briefly describe, in Section 2, the phase field model and the numerical algorithm for the study of external stress and elastic inhomogeneous systems. Then the Schmitt and Gross's theory [25] (hereafter referend to as SG) is confirmed by two-dimensional phase field simulation in Section 3. In addition, we provide the three-dimensional simulation results of multi-particle system rafting in Section 3. We end the paper with a discussion and a brief summary in Section 4.

The phase-field model with elastic energy 2.1 Formulation
In this study, we have adopted the phase-field  [26]. A binary phase separating Ni-Al alloy system exhibiting a miscibility gap [27,28] is considered in this model. In equilibrium state, the alloy consists of two phases (matrix and precipitate) in the miscibility gap, and the compositions of these phases are expressed by and , respectively. We rescale the composition variable C' to yield the scaled composition [26]: (1) By defining the reduced concentration c, a single phase field variable can be used to describe the microstructure of the whole simulation region. In particular, the value of concentration field variable in matrix phase is zero, and that for the precipitate phase equals to 1. Given an initial microstructure , its future evolution can be studied by solving the following version of the Cahn-Hilliard equation [26]: Where M is the mobility, t is evolution time,  ch and  el are the chemical and elastic contributions to the chemical driving force, respectively, and are defined as the variational derivatives of chemical (F ch ) and elastic (F el ) free energies, respectively.
The chemical free energy F ch is given by the following expression: Where  is the gradient energy coefficient related to surface energy and f o (c) is the bulk free energy density.
We assume the mobility M in Eq. (2) and the gradient energy coefficient  in Eq. (3) to be constants; thus, the interfacial energies and the diffusivities are isotropic.
The elastic contribution to the total free energy is given by (4) Where is the elastic stress and is the elastic strain, which is given by: The eigenstrain is assumed to be dilatational and an explicit function of composition [26]: Where  T is a constant that determines the strength of the eigenstrain, is the Kronecker delta and (c) is a scalar function of the composition.
The elastic modulus tensor is also an explicit function of the composition. It is convenient to describe C ijkl using the following expression: Where (c) is a scalar function of the composition. In addition, the entire macroscopic system could be subjected to a homogeneous stress state  A . To obtain the elastic energy, and hence the elastic chemical potential, we have to solve the equation of mechanical equilibrium, in (8) The fast Fourier transform (FFT)-based iterative method is applied to solve the above set of partial differential equations and the detailed solution can be found in Ref. [26].

System parameters
We use 256△x×256△y computational grids for 2D simulations. The simulated system size is 256×256 nm 2 and the periodic boundary condition is assumed in all dimensions. A thermal fluctuation [-0.001, 0.001] is introduced into the initial composition in the simulation. The system is strained in the axial direction with an applied strain of  11 =0.01. The misfit strains (eigenstrains) are assumed to be o = ± 0.01 , in which is the Kronecker delta function. All the parameters used in our computation are listed in Table 1. The elastic constants are presented in terms of the average shear modulus G, the Poisson's ratio and the anisotropy parameter , which are related to the circular averages of the Voigt constants, , and :

H. Mao et al. / JMM 55 (1) B (2019) 101 -110
We define an inhomogeneity ratio as the ratio of the shear modulus of the two-phase p and m: = / . If is equal to unity, the system is elastically homogeneous. If it is greater (less) than unity, the precipitate phase is harder (softer) than the matrix. In all our calculations, we use: v p = v m and A z = A z ;

Two-dimension (2D) simulation of rafting
The 2D rafting behavior in SG theory is summarized schematically in Fig. 1. The applied stress is along the x-axis, and the eigenstrain is assumed to be positive. Regions 1 and 3 correspond to a hard particle under tensile stress and a soft particle under compressive stress respectively, which leads to P-type rafting; in a similar manner, a hard particle under compressive stress (region 2) and a soft particle under tensile stress (region 4), lead to N-type rafting. The rafting process of elastic inhomogeneous system subjected to external uniaxial tensile stress is simulated by the phase field method. The results are in agreement with the theoretical predictions, as shown in Fig. 2.
The results of Fig. 2b, c show that the coarsening process of precipitated phase has a very obvious elongation orientation under the external stress, the particles are arranged regularly in the matrix phase, and the rafting orientation is basically along <10> direction. The tensile stress direction perpendicular to the particle elongation direction in Fig. 2b is N-Type rafting and the tensile stress parallel to the particle elongation direction in Fig. 2c is P-Type rafting. The results show that the elongation of soft precipitates in [0 1] direction is larger than that of hard particles in [1 0] direction, which means it is easier to coalescence the soft precipitated phase particles into large particles and form a long and narrow raft structure in the process of raft formation.
As shown in Fig. 3, it is seen that the number of precipitated particles has passed through two stages: the number of particles firstly increases steeply from 0 to the maximum in the first stage, and then the number of particles decreases in the second stage but with a slower rate than the first stage. These two stages of rafting phenomenon have been experimentally observed in Al-Cu alloys by Eto et al [29][30]. The authors found that the effects of stress on precipitation direction originated from nucleation  stage by two-stage aging process. The first stage is nucleation process. The precipitated particles appear due to the composition fluctuation, and the number of particle increases to the maximum after a short time.
It can be observed from the diagram that the peak value of the particle number in EH system is greater than that of soft and hard particles. The main reason for this phenomenon is that the precipitated particles rarely collide and fuse in the elastic homogeneous system. On the other hand, the external stress and the inhomogeneous elastic modulus of the system have important influences on the evolution morphology of precipitated particles. There are many soft elastic field channels in the elastic inhomogeneous system, and it is easier for the particles to choose the soft channel to merge and grow up. The peak value, particularly of the number of soft particles is higher than that of hard particles, demonstrating that the critical nucleation radius of soft particles is smaller than that of hard particles. As a result, the number of soft particles increased more rapidly than that of hard particles. This conclusion can also be verified from the histogram of the distribution of soft and hard particles. As shown in Fig. 4a (t=500 time steps), the number of small soft particles is more concentrated than that of hard particles. The second stage is coarsening. As particle surface area increases, the competition between bulk free energy and surface energy causes the impingement and coalescence between adjacent particles, which leads to the decrease of the particle number as shown in Fig. 4a (After t=500 time steps). The histogram of particle distribution of two specified evolution time (t =500 time steps and t =5000 time steps) also shows that the number of particles decreases with the coarsening process (from the maximum value of 170 reduced to 40), and the particle size increases significantly with evolution (from the maximum value of 400 increased to the maximum of 1000). During the process of coarsening, the number of hard particles gradually begin to exceed that of soft particles. It shows that the coalescence rate of hard precipitated particles is smaller than that of soft particles. In another word, the performance of hard particles to resist raft is better than that of soft particles. In general, in two-dimensional simulation, the average particle size changes during the coarsening. It can be seen from Fig. 5a that the average particle size increases with the evolution time, implying that there exists Ostwald ripening that the large particles engulf  small particles. In addition, the results show that the characteristic size of precipitated phase particles under tensile stress is larger than that of coarsening under non external stress state, which indicates that external load significantly promotes the formation of rafting. In the early stage of nucleation, as shown in Fig. 5b, due to component fluctuations, the growth rate of the nucleation particles is increased and the particle feature size grow up with parabolic relation. The rate of coarsening is stable during late stage, and the cubic power of the mean particle radius is proportional to the evolution time . The classical phase coarsening theory was proposed by Lifshitz, Slyozov, and Wagner [31] who investigated the precipitative phase coarsening behavior of diffusion control. This theory was referred to as LSW theory. According to the LSW theory, the evolution of the characteristic scales of the precipitates satisfies:

H. Mao et al. / JMM 55 (1) B (2019) 101 -110
and t are the characteristic length and precipitation time of the precipitated phase at the beginning of the coarsening; is the roughening rate constant, and its value depends on the diffusion coefficient, interface energy and molar volume.
The result of 2D simulation has a similar linear relationship with Lifshitz-Slyozov-Wagner (LSW) classical coarsening theory , the value of at this point is 0.004. Nevertheless, a well quantitative match cannot be expected in 2D simulation since the particle characteristic size depends on the area parameter of a certain section. In fact, the two-dimensional simulation only shows that rafting has a strong orientation characteristic. This study can't obtain the morphology of precipitated particles accurately from the 2D simulation results. This result is only a display of two-dimensional plane, but can't determine the shape of the space is needleshaped, Strip, flat, or columnar structure? Quantitative information about the volume fraction of concern can't be obtained exactly. This method lacks the comprehensive statistics for complete information of the overall particle size, especially the calculation of the characteristic size of strongly anisotropic growth particles, such as rafting.

Three-dimension (3D) simulation of rafting
In order to overcome the limitation of twodimensional phase field simulation of alloy during the aging process, and compare the predictions of our model with the experimental results, the 3D phase field method is used to model the precipitation process of multi-particles without external stress. All simulations were conducted in the domain of 128 △x × 128 △y × 128 △z, which corresponds approximately to 128×128×128 nm 3 . The classical theoretical model for growth describes the general law of the second phase particle coarsening through diffusion with uniform size. According to Gibbs Thomson equation, the solubility of small particle is always larger than that of large particle. Figure 6 shows the diffusion growth of the multiparticle during the aging precipitation without external stress. Fig. 6a, c correspond to aging conditions at t=1000, 3000, 5000 time steps, respectively, and all the precipitated particles are in the irregular arrangement state. It can be seen from Theoretically, only considering the interfacial energy affect, the interfacial energy and the surface area of particles increase in proportion with the particle growth, and particles with spherical or small shape are more favored. As the coarsening goes on, the growth rate of particles, which are larger than critical size is positive, and that of particles smaller than the critical size is negative. As shown in Fig. 6b, some adjacent particles begin to collide and fuse during the growth process. Particles impinge and coalesce to form bulky particles, and a few long striped particles appear, while some smaller particles are gradually reduced in volume to dissolve and disappear because of their high solubility. It is seen from Fig. 6c that several large particles have emerged. At the same time, the small particles continue to dissolve and disappear. It is seen from the evolution morphology that the coarsening process of the alloy is motivated by the conventional Ostwald aging process and the merging of adjacent particles.
Since the precipitated particle is not completely spherical, the average particle size is used as the characteristic scale for measurement. It can be seen from Fig. 7 that R 3 of the particle increases with evolution, whether in the stress state or non-stress state. Under the action of stress, the increasing rate of the characteristic size is greater than that in the nonstress state, and there is a linear relationship between the cubic power of characteristic size and the evolution time in no stress system which is consistent with the two-dimensional simulation. But, the value of slope =0.0025 in 3D simulation is slightly less than the result of the 2D simulation =0.004. The growth rate of the particles is lower than the rate of 2D simulation when the volume effect is taken into account.
In particular, as shown in Fig. 7, the particle growth rate is lower than that in free growth after the time step t=4×10 3 , and the characteristic size of the particles is eventually the same, which is different from the 2D simulation. The composition content of the final precipitated phase is determined by the thermodynamic equilibrium property. With the same thermodynamic equilibrium condition, the final component content of the precipitated phase is rarely affected by the external stress in elastic homogeneous system.
Essentially, the shape of the precipitation phase in the aging process depends on the competition and coordination between the interface free energy and the elastic distortion energy. In the elastic constrained system, the finite lattice mismatch can cause significant internal stress in the system, which has a strong impact on the coarsening kinetics. In addition, the external stress will change the stress distribution in the system and thus affect the morphology of the precipitated phase particle size, precipitates volume fraction, etc. This study attempts to explore the rafting mechanism in the coarsening process by studying the interactions among the precipitated phase, the matrix phase and the external stress in elastic inhomogeneous system.
It can be observed that the precipitate particles are suspended irregularly in the matrix phase from Fig.  8a, c, and the shape of the precipitated phase particles is mostly round, oval, and lath. From Fig. 8a, c, the number of soft particles is lower than that of the hard particles, but the volume of soft particles is larger than that of the hard particles. Accurate precipitate particle information cannot be obtained from plane (2D) data. In order to observe the evolution morphology of multi-particle precipitates more directly, the shape patterns of precipitate particles after removing the matrix phase are shown in Fig. 8b, d. It can be seen from the two images that the precipitation behavior is different from the coarsening process in the non-stress state studied above. Under the action of external stress, the aging process of both soft and hard particles shows a strong orientation along <100>. It suggests that the coarsening of soft particles shows a large elongation ratio along <100> under the interaction from both external stress and elastic distortion stress, and most of the soft particles are flat in shape, the coarsening of the hard particles has a greater elongation in the <001> direction, and most of which exhibit a cylindrical shape. Moreover, some particles form large flat plate-shaped particles along <100>, which is intuitively similar to many "rafts" made of precipitated phases in the matrix. Compared with the growth process of precipitation phase in the non-stress state, the effect of external stress on the elastic distortion energy is considered. With the growth of particles, the elastic energy and interfacial energy increase proportionally with the volume of particles. However, the elastic distortion energy gradually dominates in the process of phase separation, and with the release of elastic distortion energy the precipitation phase is distributed regularly along the elastic "soft" direction. It is found from Fig. 8b that the elongation direction of soft particles is perpendicular to the stress direction, which is N-Type rafting, and the coarsening orientation of hard particles is parallel to the applied stress called P-Type rafting. The numerical simulation results are in good agreement with SG theory. ) of the particles increases, and it reveals that the average volume of the soft particles in each time step is larger than that of the hard particles. It leads to the same conclusion that hard particles are stronger in resisting rafting and retaining their compact shapes than the soft ones. In the multi-particle system, the deformation resistance of the interface is very significant, the grain boundary hinders the dislocation movement during the plastic deformation, and the deformation of each grain is restrained by the surrounding grain. Therefore, the strength of polycrystals always increases with the grain refinement (that is, the increase of grain boundary area).
Compared with the Ostwald ripening process in the non-stress state, it is found that the rafting process of precipitation phase no longer follows the classical LSW roughening theory due to the external stress.
When the anisotropic parameter Az is equal to 1, the elastic constant of precipitated phase is isotropic. If the elastic constant is greater or less than 1, the elastic constant has cubic anisotropy. For the system in which interface energy and distortion energy interact, Gururajan and Abinandanan established the theory of shape transition of breaking symmetry by minimizing the two theoretical formulas [26], which involves the loss of normal mirror symmetry with x axis and y axis. For two-dimensional precipitated particles, when Az > 1, the < 10 > direction is the softest, and the < 11 > direction for Az < 1 system is the softest. It can be seen from Fig. 10 that when Az = 0.3, whether hard or soft precipitated particles, the rafting orientation is no longer along the < 100 > direction, but it is along the < 111> direction under the action of external stress. It can be concluded that this elongation direction is the softest rafting direction according to Gururajan and Abinandanan theory [26].

conclusions
In this paper, 2D and 3D phase field models are applied to simulate the coarse kinetic process of precipitation phase particles in elastic inhomogeneous alloy. In the elastic homogeneous system, only the interfacial diffusion relationship between precipitated phase particles and matrix phase is considered, and the results show that the Ostwald ripening process accords with the classical LSW coarsening theory. Compared to the 2D simulation, 3D simulation results have more physical significance when the volume effect is taken into consideration. Under the action of external stress of precipitated phase particles, if the rafting and shape stability of precipitated phase are predicted based on simple elastic analysis, the predicted results are not always consistent with LSW coarsening theory. The applied stress changes the distribution of elastic strain energy in the parent phase and precipitated phase. It leads to a rafting structure with a particular extension direction and a regular arrangement by the collisions and mergers between precipitated particles perform in the softest channel of the elastic stress field. Most of the soft particles undergo the transformation from spherical to flat shape under stress, while the hard particles undergo the transformation from spherical to long cylindrical shape. The simulation results are in good agreement with the experimental results, and the particle distribution, morphology and spatial correlation are also in good agreement. Quantitative simulation results are expected to play a significant role in the field of material performance prediction.