Review of γ ’ Rafting Behavior in Nickel-Based Superalloys: Crystal Review of ’ Rafting Behavior in Nickel-Based Superalloys: Crystal Plasticity and Phase-Field Simulation Plasticity and Phase-Field Simulation

: Rafting is an important phenomenon of the microstructure evolution in nickel-based single crystal superalloys at elevated temperature. Understanding the rafting mechanism and its effect on the microstructure evolution is of great importance in determining the structural stability and applications of the single crystal superalloys. Phase-field method, which is an excellent tool to analyze the microstructure evolution at mesoscale, has been gradually used to investigate the rafting behavior. In this work, we review the crystal plasticity theory and phase-field method and discuss the application of the crystal plasticity theory and phase-field method in the analysis of the creep deformation and microstructure evolution of the single crystal superalloys.


Introduction
Turbine blades are critical structural components in modern aircraft engines and function under extreme service conditions, such as high temperature, high pressure, and high stress level.The performance of turbine blades plays an important role in determining the stability and service life of aircrafts.There is a great need to develop high-performing materials for turbine blades.
Nickel-based single crystal superalloys are currently used in turbine blades due to their excellent mechanical properties at elevated temperature, including high mechanical strength and creep resistance.The excellent mechanical properties of the superalloys are attributed to the two-phase microstructure, i.e., γ' precipitate phase and γ matrix phase, with L12 ordered γ' phase being coherently embedded in γ phase of face-centered cubic structure (Figure 1a).The ordered γ′ phase with high volume fraction (about 60-70% in most cases) serves as the strengthening phase in the superalloys and hinders the motion of dislocations during creep deformation [1][2][3].Under mechanical loading at high temperature, the morphology of the γ' precipitate evolves from cubic shape to plate-like shape, as shown in Figure 1, which is referred to as directional coarsening or rafting.The rafting is a complex process, which can contribute to the hardening of the materials, since the raft structure inhibits dislocation climb around the γ' phase and restricts the motion of dislocations in γ channels [2][3][4].However, there are reports that the rafting behavior can also lead to the softening of materials [5,6].The onset of rafting widens the width of some γ channels which reduces the blocking effect of Orowan stress and causes the bowing of dislocations in the γ channels (Figure 2).This behavior results in the increase of the plastic strain in the corresponding γ channels and the nearby γ/γ' interfaces.Note that the plastic deformation around the γ/γ' interfaces plays a key role in the rafting process.Therefore, the rafting can lead to the strengthening and softening of nickel-based superalloys during creep, which necessitates the understanding of the evolution of the raft structure under mechanical loading at elevated temperature.In the field of view, some dislocations are present in the γ' phase, which is attributed to the dislocation climb [7].(Reproduced with permission from [8] © 2020 Elsevier.) There are a variety of methods available to investigate the rafting behavior, which can be divided into two categories: experimental methods [1,9] and numerical simulation [10][11][12][13].In general, it is very costly and impractical to use conventional experimental methods to study the rafting behavior due to the need to "record" the changes of the microstructures during creep deformation.Numerical methods, such as finite element method [10], Monte Carlo method [11], and phase-field method [12], have become a major approach to understand the microstructure evolution associated with the rafting process.Among the numerical methods, the phase-field method has become an important technique to analyze the microstructure evolution, which involves phase transition at mesoscale for a long-time period without explicit interface tracking [14][15][16][17][18][19][20], and has been used efficiently to simulate the rafting process.Realizing the presence of plastic deformation during the rafting, the theory of crystal plasticity is usually incorporated in the phase-field method to investigate the plastic deformation concurrently presented during the rafting.
In this review, we briefly introduce the rafting behavior in nickel-based single crystal superalloys, together with the crystal plasticity theory and the phase-field method.The models, which are based on the crystal plasticity theory and the phase-field method, respectively, for the analyses of plastic deformation and microstructure evolution are summarized.The applications of the crystal plasticity theory and the phase-field method in studying the rafting behavior are then discussed.Finally, a brief summary is given.

Types of Rafting
Rafting is an important phenomenon of the microstructure evolution in L12 hardened superalloys, such as nickel-and cobalt-based superalloys [21][22][23].There are extensive studies on the formation mechanism in single crystal superalloys, which are also applicable to individual grains in polycrystalline superalloys [24].Here, we are mainly focused on the rafting behavior in single crystals.
Under mechanical loading in <001> direction, raft structures can form in single crystals in the direction either parallel (P-type) or normal (N-type) to the direction of external loading.The factors controlling the rafting orientation mainly include the sign of applied stress, σ A (positive for tensile loading and negative for compressive loading), the γ/γ' lattice misfit, δ, and the difference in the mechanical behavior between γ and γ' phase, ∆G [25,26].Here,  = 2  −   +  , with  and  being the lattice constants of γ' phase and γ phase, respectively.∆ is usually used to represent the hardness ratio of the γ' precipitate to the γ matrix with ∆ > 1 for hard precipitates and ∆ < 1 for soft precipitates.The P-type raft structures form under  (1 − ∆) < 0, and the Ntype raft structures form under  (1 − ∆) > 0, as shown in Figure 3. Schmidt and Gross [25] summarized several cases for the raft structures presented in regions (1-7) in Figure 3, in which regions (5-7) represent extreme cases.Note that the raft structures presented in region (5) are only observed in very soft precipitates under compressive loading, and the raft structures presented in regions ( 6) and ( 7) are observed in two-phase materials with different Poisson's ratio and/or Zener anisotropy parameter [26].The rafting direction plays an important role in determining the mechanical behavior of superalloys.In most cases, the N-type raft structures can reduce the low-cycle fatigue strength of materials and the P-type raft structures generally possess high creep resistance and low-cycle fatigue strength [27][28][29].Note that the N-type raft structures can also increase the creep resistance by effectively hindering the motions of dislocations at high temperature under the action of small stress [28] in comparison to the corresponding alloys with cubic γ' precipitates.Mughrabi and Tetzlaff [28] reported that the prerafting to form the P-type raft structures in superalloys with negative γ/γ' lattice misfit also resulted in the improvement in fatigue strength at elevated temperature.Compromise needs to be made in controlling the rafting behavior of superalloys with negative γ/γ' lattice misfit.
Commercial nickel-based superalloys possess negative lattice misfit at room temperature, whose magnitude decreases with the increase of environmental temperature [30,31].To better describe the effect of lattice misfit on the rafting behavior, cobalt-based superalloys are taken as an example for comparison.Cobalt-based superalloys possess positive lattice misfit at room temperature, whose magnitude also decreases with the increase of environmental temperature, but the lattice misfit remains positive even at temperature of 1000 °C [22,31].Negative lattice misfit in nickel-based superalloys and positive lattice misfit in cobalt-based superalloys at such a high temperature result in the formation of N-type and P-type raft structures under tension, respectively.The opposite behavior is observed under compression.
It is known that a superalloy of single crystal exhibits different mechanical behavior from the corresponding superalloy of polycrystal, and its creep and fatigue properties are extremely sensitive to crystal orientation [32,33].The complex geometry of single-crystal turbine blades and the multiaxial stress state under service necessitate careful characterization and analyses of the mechanical behavior and γ/γ' microstructure evolution of single crystal superalloys in typical crystal orientations, including <011> and <111>.Under a tensile loading in <011> direction, directional coarsening of γ' phases was found parallel to (010) plane, while equiaxed coarsening behavior was observed in (100) plane, which exhibits a 45-degree rafting [34,35].Under a tensile loading in <011> direction, no rafting behavior was present [34,36].However, there always exists a slight misorientation between actual loading direction and <011> or <111> crystal orientation, which can lead to the formation of raft structures different from theoretical results.If the loading direction slightly deviates from <011> direction, see Figure 4 as an example, directional coarsening of γ' phases is still found parallel to (010) plane, but the coarsening behavior in (100) plane is not equiaxed, tending to align along [010] or [001] direction.In the case of loading direction slightly deviating from <111> direction, rafting behavior along one main direction is also observed (Figure 5).These unexpected phenomena were also reported in some works [37][38][39], which is believed to be related to asymmetric stress distribution in γ/γ' microstructures, caused by a slight misorientation from the loading direction and microdefects randomly distributed in the material.Therefore, the rafting in a superalloy of single crystal is sensitive to the effective loading direction.

Kinetics of Rafting
In the heart of rafting is a stress-limited diffusion process, depending on applied stress, σ A , lattice misfit, δ, and the difference between the mechanical property of γ and γ' phases, ∆G [3,40].Both applied stress and inherent misfit stress play important roles in the distribution of internal stress in the grains/crystals consisting of γ and γ' phases due to the dependence of chemical potential on local stress state [40].Under tensile loading, the difference of the chemical potentials between the γ' phase and the γ phase (γ channels) drives the forming atoms of the γ' phase (Al, Ti, Ta, etc.) in horizontal γ channels into vertical γ channels and drives nearly insoluble atoms (Co, Cr, Mo, Re, W, etc.) away from vertical γ channels.This results in the formation of raft structures perpendicular to the tensile loading, as shown in Figure 6 [41,42].Assuming that the stress-limited diffusion process is the dominant mechanism controlling the rafting, Fan et al. [34] proposed a von-Mises stress-based criterion in determining the rafting direction.Under concurrent action of external stress and misfit stress, the von-Mises stress in different γ channels can be calculated according to the formulations in Table 1.Dislocations can easily accumulate in γ channels under a large stress, leading to the relief of the misfit stress through the migration of the γ/γ' interface and the directional coarsening in the γ channels with maximum von-Mises stress [34,43].
Table 1.Calculation of the von-Mises stress in γ channels under uniaxial loading.Here, σ0 is applied stress in the global coordinate system, σi is the magnitude of the misfit stress, and α is a reduction factor.Adapted from [34].Geometrical characteristics of the microstructures have been used to determine the degree of rafting, such as the width of γ channels [34,44], the dimension of γ' phase [45], and the shape of γ' phase [46].The images of the microstructures can be captured through scanning electron microscopy (SEM) and transmission electron microscopy (TEM), and the geometrical characteristics of the microstructures can be then analyzed via image-processing algorithms.There are a few algorithms commonly used to analyze the SEM and TEM images, including Connectivity number method [46], AutoCorrelation Method [47,48], Rotational Intercept Method [48], Fourier analysis [9], and Moment invariants method [49].
Fedelich et al. [44] introduced a dimensionless variable ξ for the analysis of rafting as where (t) is the width of the γ channel at time t,  and  are the channel widths of the cubic structure and the raft structure, respectively.Their numerical values are correlated to the volume fraction of γ' phase,  , and the microstructure periodicity, λ.The dimensionless variable, ξ, varies in a range of 0 (initial cubic morphology) to 1 (complete raft morphology).
Tinga et al. [45] proposed an evolution law for the dimension of γ′ phase under the action of a multiaxial stress as where  are the dimensions of γ′ phase along three orthogonal directions,  are the diagonal components of deviatoric stress tensor,  is von-Mises stress, and  is thermal shear-activation volume, Q, R, and T are activation energy, gas constant, and absolute temperature, respectively,  * and  are constants, and  is the dimension of γ′ precipitate in cubic shape.Desmorat et al. [50] used the width of γ channel and the dimensions of γ' phase as internal variables in the framework of thermodynamics with elasto-visco-plasticity and Orowan stress (dislocation mechanics).They calculated the Orowan stress  , which acts as a resistance to the dislocation motion in γ matrices, during the rafting in the following formulation, where  or is a material constant ranging from 0.238 to 2.15 [51], and μ and b are the shear modulus and the magnitude of Burgers vector, respectively.Their simulation results provided quantitative description of the rafting behavior in nickel-based single crystal superalloys.There exists the interaction between rafting and dislocation motion.Rafting leads to the accumulation of dislocations in the γ channels, resulting in the relief of the misfit stress through the migration of the γ/γ' interface and the directional coarsening of γ' precipitates.Without mechanical loading, rafting can also prevail at high temperature if the plastic strain in superalloys reaches a threshold value [52,53].The deformation field (plastic strain) is associated with the presence of dislocations in crystal, even though there is no external loading.The dislocation motion in superalloys at high temperature relieves the misfit stress and promotes the directional coarsening of γ' precipitates through rafting.It is the plastic deformation in the γ matrices (channels) that determines the rafting process and the microstructure evolution.

Crystal Plasticity Theory
Crystal plasticity theory is based on the work from Taylor and his coworkers [54].Their work suggested that plastic deformation of crystalline metals is correlated to the change of crystallographic structure.There are two major aspects contributing to the change of crystallographic structure: one is the lattice distortion, and the other is the plastic deformation from the gliding of dislocations in slip planes.There is a significant amount of dislocations in crystalline metals, which cannot be easily described by traditional continuum mechanics [54].The large amount of dislocations in crystalline metals (≈10 7 per cm 2 at annealed state) makes it reasonable to use the concept of continuum mechanics in the analysis of plastic deformation controlled by dislocation motion.
Assuming that dislocation gliding occurs uniformly in a grain/single crystal, Rice and Hill [55,56] proposed a kinematic formalism for the plastic deformation of crystals.They divided the deformation gradient tensor F for a crystal into F e , representing the lattice distortion and rigid rotation, and F p , representing dislocation gliding (Figure 7) as The plastic deformation rate tensor due to the dislocation gliding is calculated as Here, L p is the velocity gradient for plastic deformation, consisting of the contribution of the shear strain rate  ˙( ) on all active slip systems as where  ( ) and  ( ) are the unit slip direction and normal vectors of the slip plane for the -th slip system, respectively, and N is the number of slip systems.L p can be further divided into a symmetric part, i.e., plastic deformation rate tensor, D p , and an antisymmetric part, i.e., spin tensor, W p , as The plastic deformation rate tensor, D p , represents the incremental change in the deformation behavior of the material, and the spin tensor, W p , represents the gliding-induced change in the crystal orientation.Equations ( 5)-( 8) lay the foundation to establish the relationship between the deformation at macro-scale and the shear strains of individual slip systems in a grain/crystal during plastic deformation.
In addition to the kinematic relations, constitutive equations, which capture the microstructure evolution, such as the rafting and dislocation activities, and correlate stresses to strains, need to be developed in order to completely describe the deformation behavior of a grain/crystal.In the following, we present two classes of constitutive models: one is referred to as phenomenological constitutive models, and the other is referred to as physics-based constitutive models.

Phenomenological Constitutive Models
Several internal state variables are used in the development of phenomenological constitutive models.The deformation state of a material is determined by the variations of the internal state variables with thermal and mechanical loading.The phenomenological constitutive models can be categorized into power-law type (medium stress condition), hyperbolic-sine type (medium to high stress condition) and linear type or combination of linear and power-law types (low to medium stress condition).Chowdhury et al. [57] compared the applications of these constitutive models in the study of crystal plasticity.For engineering practices, the analysis of the plastic deformation of turbine blades is usually based on the power-law constitutive models.The power-law constitutive models have been widely used in the stress analysis due to that it is cost-effective in determining material parameters in the constitutive models and they are applicable in the analysis of the steady-state creep deformation of turbine blades at elevated temperature over a wide range of stresses.However, the pre-factor of the power-law constitutive models is a function of temperature due to different rate mechanisms.
For nickel-based superalloys, dislocation motion in both γ and γ' phases are the dominant mechanism for plastic deformation at elevated temperature [58,59].In the power-law constitutive models, the ratio or difference between the resolved shear stress,  ( ) , and the slip resistance,  ( ) , determines the slipping activity of the slip systems.A general mathematic relationship between the shear strain rate, the resolved shear stress, and the slip resistance can be expressed as [60]  ˙ =   ( ) ,  ( ) There are three different types of power-law constitutive models for creep deformation as [57]  ˙( ) =  ˙  ( )  ( )   ( ) (10) where  ˙ and n are the shear rate at the reference state and the rate sensitivity parameter, respectively,  ( ) and  ( ) are two threshold stresses (both evolve with the activity of crystallographic slips [57]),  ( ) is internal stress or the back stress, and K is temperature-dependent material parameter.The angle bracket "<y>" in Equation ( 12) represents the positive part of y, and the "sgn(•)" is the sign function.
For the power-law constitutive model of Equaton (10), the activity of crystallographic slips is determined by the resolved shear stress,  ( ) .The slip resistance,  ( ) , is the only internal state variable as with ℎ as the hardening matrix, which determines the hardening effect of the α-th slip system on the β-th slip system [61].For α = β, the hardening parameter, ℎ , represents the self-hardening modulus as [61] ℎ = ℎ() = ℎ ℎ ℎ   −  (14) for  ≠ , the hardening parameter, ℎ , represents the latent hardening modulus as Here, ℎ is the initial hardening modulus,  and  are the saturated shear stress and initial yield stress, respectively, and q is a constant.The cumulative shear strain, , is calculated as Note that the hardening matrix, ℎ , is also reported in a power-law form as [60,62] with n0 as a material constant.For the self-hardening modulus,  = 1, and for the latent hardening modulus,  = 1.4.
For the constitutive models of Equations ( 11) and ( 12), the activity of crystallographic slips is determined by the effective stress  ( ) −  ( ) and  ( ) −  ( ) −  ( ) , respectively.The internal stress,  ( ) , can be calculated by the following equation [63] In Equation (18),  ( ) is the determinate internal stress at a temperature, and  ( ) is a control variable [57].The constant d is a recovery parameter.The flow accumulation function of   ( ) is calculated from the accumulated shear strain of the -th slip system of  ( ) and the constants of  and  in Equation (20).It needs to be pointed out that some works [12,19] replaced the term of   ( )  ˙( ) with  ˙( ) in Equation (19) and did not consider the flow accumulation function of   ( ) .
In general, the use of the effective stress instead of the resolved shear stress for the constitutive models of Equations ( 11) and ( 12) makes it possible to analyze the plastic deformation in materials consisting of complex microstructures involving multiphases [64,65].This is because the internal stresses from the complex microstructures do not directly contribute to the activities of crystallographic slips and needs to be deducted.Additionally, the increments of the internal state variables in Equations ( 11) and ( 12) have provided the basis to accurately describe crystallographic slips and determine the threshold stress.All of these suggest that the power-law constitutive models of Equations ( 11) and ( 12) can likely provide better correlations between stresses and strains for the analysis of the plastic deformation of the nickel-based superalloys.

Physics-Based Constitutive Models
In the phenomenological constitutive models, the threshold stresses,  ( ) and/or  ( ) , whose evolution follows a hardening rule, are used to represent the contribution of dislocation motion.However, it is very difficult, if not impossible, to experimentally determine the hardening rule and to validate the hardening rule under service-like conditions.Additionally, these constitutive models fail to capture the orientation dependence of the mechanical behavior of single crystals [66] and the hardening rule of materials at the micron scale [67].There is a great need to develop physics-based constitutive models, which use dislocation density as an important internal state variable.
There are various physics-based constitutive models, which incorporate dislocation density in the theory of plasticity for the analysis of the plastic deformation of crystalline materials [68][69][70][71].In the heart of the physics-based constitutive models is the evolution of statistically stored dislocation (SSD) and geometrically necessary dislocation (GND) during plastic deformation.The SSD is associated with the "homogeneous" deformation, while the GND is associated the "inhomogeneous" deformation in single crystals only at small length scale [68].
SSDs are quantified by dislocation density, ρ, (line length per unit volume), Burgers vector, b, and unit vector of dislocation-line segment, t.The magnitude of the Burgers vector is discrete and related to the lattice constant of crystal.For simplification, the unit vector of dislocation-line segment is sometimes limited to a finite set, which makes it easy to include dislocations in numerical calculation.
The Orowan relationship instead of the phenomenological constitutive relationship is used to correlate the plastic shear rate with dislocation density in a physics-based constitutive model as [60] where  is the density of mobile dislocations and is calculated from the SSD density,  , and b and  are the magnitude of Burgers vector and average velocity of the mobile dislocations, respectively.One specific form of Equation ( 22) is [68,69]: where the subscripts, e and s, represent edge and screw dislocations, respectively, and the symbols, "+" and "−", represent the polarity of the SSD density.The density of mobile dislocation,  , in Equation ( 22) can also be divided into two portions:  for the dislocations parallel to slip planes, and  for the dislocations perpendicular to slip planes as [70,71] where kB is Boltzmann constant, c1 is a constant for the evolution of dislocation density, and  ( )( ) represents the interaction strength between different slip systems.The dislocation interaction plays an important role in determining the plastic deformation in single crystal superalloys.There are nucleation and annihilation of dislocations during plastic deformation, which determine the evolution of dislocation densities [68,70].The SSD alone is not enough to describe the plastic deformation of crystalline materials.If heterogeneous two-phase microstructure exists in a material, a local strain gradient related to the activity of GNDs is always generated between the strengthening precipitates, such as γ' phases in nickel-based superalloys [72].In crystal plasticity models, GNDs can be obtained from slip gradients, and they are further divided into the edge dislocation,  , , along the slip direction and the screw dislocation,  , , perpendicular to the slip direction as [67]  , ( )  = − ( ) •  ( ) (28) where  ( ) =  ( ) ×  ( ) ,  is a unit vector parallel to the Burgers vector, and  is the unit normal to the slip plane.
With the SSD density and GND density, the internal stress or back stress,  ( ) , and the slip resistance,  ( ) , in phenomenological constitutive models can be calculated as [67]  ( ) =   ( )( )  ( ) (30) where  is Poisson's ratio, R is a length scale in the calculation model [73],  ( )( ) is the interaction matrix between two slip systems of α and β with six types of interaction in the matrix [74]. ( )( ) is the interaction coefficient as [73]  ( )( ) = 1 (, ) = (4,13), (6,18), (8,17), (9,15), (10,16), (11,14), (1,16), (2,17), (3,18), (5,14), (7,13), (12,15) 0 ℎ Here, α (= 1 to 12) represents the edge dislocation, and α (= 13 to 18) represents the screw dislocation. ( ) in Equation ( 30) is the total dislocation density on the β-th slip system.Note that Tinga et al. [51] considered the contributions of the SSD density and GND density to the total dislocation density, respectively.Some works incorporated dislocation dynamics, such as discrete dislocation dynamics (DDD) [16,17] and continuum dislocation dynamics (CDD) [18], in the analysis of plastic deformation and microstructural evolution.The DDD models are based on discrete description of dislocation motion and require sufficiently fine grid spacing and great computational cost in the simulation of dislocation activities.The CDD models are based on a continuum quantity (dislocation density) instead of individual dislocations and need much less computational cost [18].However, the simulation with either type of dislocation dynamics models still costs more computational effort than that with the crystal plasticity models, which incorporate the dislocation activity in a phenomenological or empirical form [19,20]. Thus, crystal plasticity models are widely used to account for the microstructure evolution and plastic deformation and to provide valuable information for engineering applications [57,75].Table 2 presents the comparisons of the plasticity models used in the rafting analysis.

Phenomenological constitutive models
Pros: Be cost-effective in determining material parameters and applicable in engineering calculations.Cons: Fail to capture the orientation dependence of the mechanical behavior of single crystals; difficult to experimentally verify the hardening rule used in the constitutive models.
Rafting with creep damage [8,19] Physics-based constitutive models Pros: Be able to model the microstructure evolution and include the contribution of dislocations.Cons: Fail to explicitly capture the motion of dislocations.
Coupling between rafting and crystal plasticity with dislocation densities [67] Discrete dislocation dynamics models Pros: Explicitly describe the dislocation distribution during microstructural evolution.Cons: Require sufficiently fine grid spacing and great computational cost in simulation.
Distribution of plastic strain in γ-channels and its effect on rafting [16] Continuum dislocation dynamics models Pros: Consider average distribution of dislocations and need less computational cost.Cons: Difficult to be compared with phenomenological constitutive models in engineering calculations.
Effect of initial dislocation density on rafting [18] 4. Phase-Field Models

Phase-Field Method
Phase-field method is based on the description of diffuse-interface model, and the concept of diffuse interface is derived from van der Waals [76] and Cahn and Hilliard [77].In contrast to conventional sharp-interface models (field variables are discontinuous at interface), the interface in the diffuse-interface model is represented by a series of consecutive values [18].Assigning different values to different phases (for example, γ' strengthening phase and γ matrix phase in nickel-based superalloys), the interface region is implicitly given [78].Such a method makes it possible to analyze the evolution of complex microstructures, including precipitation, dissolution, coarsening, connection, and topological inversion of γ' strengthening phases [79,80].Additionally, the physical characteristics, i.e., the γ/γ' lattice misfit, elastic constants and dislocation densities, and the mechanical behavior, i.e., heterogeneous elasticity and plastic deformation, of superalloys can be well incorporated in the mathematical framework of the phase-field method.The phase-field method has become an excellent technique to study internal stresses and plastic deformation and illustrate the mechanisms for the microstructure evolution in nickel-based superalloys.
The application of the phase-field method in the simulation of microstructural evolution was initially based on an elastic framework, where the mechanical stress, σ A , the γ/γ' lattice misfit, δ, and the difference of mechanical properties between γ and γ' phases, ∆, were incorporated in the model to reveal the rafting behavior [14,15,26].The elasto-plastic frameworks, which take into account the plastic deformation in γ matrices during creep, were later developed [8,12,19].
Superalloys are generally multiphase and/or multicomponent materials, which require large number of field variables or physical quantities in the phase-field method in addition to mechanical deformation (elastic and elasto-plastic deformation).They provide a practical application of the phase-field theory in the understanding of the microstructure evolution of multiphase, polycrystal, and multicomponent materials [81].The combination of the phase-field model and the crystal plasticity likely opens a new avenue for the study of the mechanical deformation and microstructure evolution of superalloys under concurrent action of thermal and mechanical loading.

Ni-Al Binary System
Nickel-based superalloys consist of multielements, including Ni, Al, Re, Co, etc., while Ni and Al are two major elements determining the microstructure evolution [12].A Ni-Al binary system has been used in different phase-field models to simulate the γ/γ' microstructure evolution in the nickelbased superalloys.
(, ) is defined as the concentration of Al, which is used to distinguish γ' phase (mainly Ni3Al) from γ phase (mainly Ni).For a given concentration, C, it is equal to  in the γ phase and  in the γ' phase, and the region with the Al concentration changing from  to  corresponds to the γ/γ' interface.Note that the Ni-Al phase diagram is used to determine the numerical values of  and  at a given temperature [82].A dimensionless parameter,  , is sometimes used instead of the concentration, C, to distinguish the γ' phase from the γ phase as [26]  =  −   −  (33) with C' = 1 representing the γ' phase and C' = 0 representing the γ phase.
The principle of maximum entropy has been mainly used in the solidification analysis of alloys.It is very difficult to calculate the change of entropy during plastic deformation.The principle of minimum free energy is more common in the study of the rafting process.With the principle of minimum free energy, the concentration field, (, ), and the field parameters,  (i = 1, 2, 3), satisfy the Cahn-Hilliard and Allen-Cahn equations as , (i=1, 2, 3) The local volume fraction of the γ' phase, f(, ), is used sometime instead of (, ) in Equation ( 34) [20,84].Note that the local γ'-volume fraction is equivalent to the dimensionless parameter, C'.The use of f(, ) makes it easy to extend the phase-field models for the Ni-Al binary system to multicomponent systems.Here, m0 and l0 are the mobility coefficient and the kinetic coefficient, respectively, and F is the total free energy consisting of chemical free energy,  , and strain energy,  .The chemical free energy can be expressed by Ginzburg-Landau functional approximation as [12] The gradients of the concentration field and field parameters in Equation ( 36) define the numerical resolution and interface thickness in the simulation [18].The gradient energy coefficients of  and  are related to interfacial energy, and their values are adjusted to ensure that the interface thickness of two-phase microstructure can represent real situation [12,85].The use of  is to distinguish the γ phase from the γ' phase and to assure the stability of four different γ' variants at  {1, 1, 1},  {1 , 1 , 1},  {1 ´, 1, 1 }, and  {1, 1 , 1 }.The strain energy,  , is calculated as where  is local elasticity modulus tensor depending on the concentration field, (, ). is a homogeneous term, depending on loading condition.For strain-control condition,  = 0; for stresscontrol condition,  = − ̅ .Here  is the components of applied stress, and ̅ is the average strain components [13].
The contribution of the hardening free energy, F vp , or plastic energy, F pl , to the total free energy during plastic deformation was also considered in some works [12,19].However, the partial derivative of F vp or F pl to (, ) in Equation ( 34) is negligible since most studies have been using the same viscoplastic parameters for both the γ phase and the γ' phase, i.e., the hardening or plastic energy function is homogeneous in two phases and independent of the concentration field, (, ).
It is worth mentioning that calculations are performed for small deformation in the phase-field simulation.In this situation, the change of the crystal dimensions is negligible and the total strain rate tensor,  , can be divided into three parts as [26,86] Here,  ,  , and  are elastic strain, eigenstrain from the γ/γ' lattice misfit , and plastic strain rate tensors, respectively.The plastic strain rate tensor,  , is obtained through the theory of crystal plasticity as [86] Then the shear strain rate  ˙( ) is calculated through phenomenological or physics-based constitutive models.However, the use of the theory of crystal plasticity for small deformation indicates that the calculated strains are much less than the strains measured in experiments.
The equilibrium equations must be satisfied all the time during the microstructure evolution.The local stress equilibrium is expressed as [12]   ~= 0  =  ~ ∈ Minimizing the strain energy function,  , with respect to the displacement or strain components, ui or̅ , under given boundary conditions yields [13] The Cahn-Hilliard and Allen-Cahn equations and the equilibrium equations constitute the main framework of the phase-field models, which involve the two-way interaction between the concentration field and the energy function (stresses).

Multiphase-Field Model
Phase-field models have evolved from simple binary models with only one order parameter to multiorder or multiphase models recently.A typical example of the multiorder phase-field model is to have four different γ' variants with three ordered parameter fields,  (i = 1, 2, 3), as discussed in the above section.Multiphase-field models are mainly used for polycrystalline materials [87] and single crystal with multiple components [41].The following discussion is focused on the multicomponent models for the analysis of the rafting behavior in single crystals.
The total free energy, F, in multicomponent models can be calculated from the integration of the strain energy density, f el , the interface energy, f it , and the chemical free energy, f ch , as [87]  = ( +  +  ) where  denotes γ matrix,  (=1, 2, 3, 4) denote different γ' variants,  is the interface energy,  is the width of interface,  =   is the repulsive potential function, and  is the concentration vector.In Equation (44),  ( ) is the bulk free energy of each phase, and  =  is chemical potential vector.
The strain energy density,  , in multiphase-field models has a similar expression to that in binary-field models as where  * and  * are the components of effective stiffness tensor and effective eigenstrain tensor, respectively [41,87].The kinetic equations for the microstructure evolutions are [41,87] Here,  is the components of the mobility-coefficient tensor, and  is the chemical mobility tensor.Comparing Equations ( 47)-( 49) of the multiphase-field model with those in the Ni-Al binaryfield model, we note similarities between both models.The driving force for the microstructure evolution in both models is the variation of the energy functions of individual phases, and the equilibrium equation needs to be satisfied all the time.
Figure 8 summarizes the basic process to numerically solve the Ni-Al binary-field models and the multiphase-field models.Step 1: Obtain initial microstructures.This can be achieved by the heat treatment of superalloys (precipitation and coarsening) at high temperature and imaging (SEM and/or TEM) [79] and/or by setting field variables from external files [8,42].
Step 2: Calculate plastic strain.The plastic strain is calculated by solving the related equations, which are based on phenomenological models or physics-based models.
Step 3: Calculate strain energy.The total strain field,  ~, and the stress filed,  , are firstly obtained by solving equilibrium equations [12,26].The elastic strains are obtained by subtracting the eigenstrain and plastic strains from the total strains.Finally, the strain energy is calculated from the elastic strains and elastic constants.
Step 4: Update the field variables by solving the Allen-Cahn and Cahn-Hilliard equations.

Uniaxial Tension
Uniaxial tension with constant strain rates has been used to experimentally study the rafting behavior in superalloys.Using strain control in the phase-field simulation, we can achieve the tensile mode in solving equilibrium equations.Figure 9 presents the microstructure evolution and plastic strain during the tension of a nickel-based single crystal superalloy.It is evident that there is no significant change in the microstructure, and there is plastic deformation near the γ/γ' interface and the center of the γ' phases.Comparing the numerical results with experimental results, we can determine the mechanical properties of the γ/γ' microstructure and the bulk phase (if two phase structures are not distinguishable) from the input parameters for the numerical calculation.Additionally, we can evaluate the strengthening effect of the γ' phase and determine the effect of the γ' phase on the mechanical properties of nickel-based superalloys on the macroscale.Cottura et al. [67,88] introduced a characteristic plastic length, , in 1D configuration to illustrate the size effect on the plastic deformation of the γ phase and the mechanical properties of the material.Their simulation results suggest that plastic strain is homogeneous in the γ phase if the characteristic plastic length is negligible.Increasing the characteristic plastic length led to the inhomogeneous distribution of the plastic strain and the decrease of the peak value of .For the value of  being larger than the width of the γ channels, no plastic deformation would occur in the γ channels.Increasing the characteristic plastic length also increased the flow stress during the tensile deformation.
Wang et al. [8] included the Orowan stress in a visco-plasticity phase-field model to analyze the size effect in a two-dimensional system.Their results show that the flow stress during the tensile deformation decreases with the increase of the γ channel width and the flow stress remains unchanged after the γ channel width reaches a critical value.Zhang et al. [7] also reported similar results for the tension of nickel-based superalloys at high temperature.They pointed out the dependence of the yielding stage during tensile deformation on the penetration of dislocations into γ' precipitates.

N-Type/P-Type Rafting
There are reports on the use of the phase-field models in the analysis of the rafting during creep deformation in the frameworks of elastic deformation [14,26] and elasto-plastic deformation [8,12,19,67].In the framework of elastic deformation, the rafting is determined by the sign of applied stress, σ A , the γ/γ' lattice misfit, δ, and the difference of mechanical properties between γ and γ' phases, ∆.The results from the phase-field simulation reveal the formation of N-type and P-type raft structures in consistence with the predictions shown in Figure 3.However, the simulation cannot capture the rafting features associated with plastic deformation [12,26] despite the good replication in the types of the raft structure.For example, Gaubert et al. [12] compared the simulated microstructures in elastic framework with the corresponding ones in elasto-plastic framework.They found raft structures with straight edges in elastic framework and raft structures with wavy-type edges in elasto-plastic framework.In addition, the plastic deformation increased the rafting rate and promoted the formation of raft structures.Wang et al. [8] presented raft structures and the distribution of plastic strain during creep deformation in elasto-plastic framework.As shown in Figure 10, the raft direction was not strictly perpendicular to the loading direction, and the plastic strain almost appeared in horizontal γ channels and around the γ/γ' interface.They pointed out that further plastic strain that occurred around the γ/γ' interface might result in the instability of raft structures.Zhou et al. [89] pointed out that the kinetics for rafting are not solely determined by the difference of mechanical properties between γ and γ' phases, ∆.They revealed that N-type raft structures can form under tensile loading even for homogeneous γ/γ' microstructures (∆ = 1) if there are dislocations present in horizontal γ channels (normal to the loading direction) at initial state.The time needed to form stable P-type raft structures is less than that to form N-type raft structures, which can be attributed to the difference in the diffusion paths of solute atoms.
According to Equations ( 34) and ( 36), the atomic migration in γ channels is controlled by the difference of chemical potentials, where forming elements of the γ' phase (mainly Al in Ni-Al binary system) diffuse from the γ channels with large chemical potential to the matrix with small chemical potential.Under compression in <001> direction, the forming atoms of the γ' phase migrate from two types of vertical γ channels (channel 1 and channel 2, see Table 1) to the horizontal γ channel (channel 3) leading to the formation of P-type raft structure.Under tension in <001> direction, atoms migrate from channel 3 to channels 1 and 2, leading to the formation N-type raft structure.
The phase-field models for Ni-Al binary system have been focused on the diffusion of Al, which is likely not enough for multielement superalloys.Multiple-component phase-field models, which contain the main chemical components of nickel-based superalloys, have been developed in order to systematically investigate the effects of alloy elements on the rafting kinetics [42].As shown in Figure 6, the results reveal that, under tensile loading, the coalescence of γ' precipitates to form raft structures causes the diffusion of Al, Ti, and Ta from the horizontal γ channels to the vertical γ channels and removes atoms of Co, Cr, and Mo from the vertical γ channels.However, some refractory elements, such as Re and W, likely accumulate near the γ/γ' interface.

Complex Types of Rafting
The types of raft structures formed are largely dependent on the loading condition (direction, tension/compression, etc.), such as the 45-degree rafting in <011> direction [37] (Figure 4).Kamaraj [40] reported the formation of raft structures with 45° angle to the shear stress in a double shear creep test.Such a type of raft structure was also observed experimentally during the creep of MC2 nickelbased superalloy at high temperature and attributed to the highly localized creep strain [90].However, the angle between the direction of the raft structure and the mechanical loading varied from the region near the fracture with highly localized creep strain to the region far away from the fracture.
Gaubert et al. [91] analyzed the microstructure evolution of <011>-oriented Ni-based superalloys using a 3D mean-field visco-plasticity model.Their results demonstrated the directional coarsening of raft structures along <100> direction.A slight deviation between the loading direction and <011> direction drove initial cubic γ' precipitates to coarsen firstly along <100> direction and extended then along <001> or <010> direction.Yang et al. [92] performed a similar study on the microstructure evolution of nickel-based superalloys with the loading direction deviated from <001> orientation.They referred to this loading mode as monoclinic loading, which can be equivalent to a shear-loading mode and a tension-loading mode.They analyzed the microstructure evolution under monoclinic, shear, and tension loadings, respectively, and found nearly no rafting behavior under the shear loading and synchronous N-type rafting for the monoclinic and tension loading under small stress.Increasing applied stress led to the formation of 45-degree rafting in the overall region under the shear loading and in the partial region under the monoclinic loading.A N-type rafting was still found under the tensile loading, and its rafting rate was smaller than that under the monoclinic loading.Increasing applied stress increased the role of the shear stress component in the rafting process under the monoclinic loading.
Ali et al. [93] combined the phase-field method with physics-based crystal plasticity model to explain the formation of a 45-degree rafting in the region with local creep strain larger than 10% under a tensile loading in <001> direction.Their results showed that creep strain in some regions was significantly higher than average creep strain.A large amount of geometrically necessary dislocations was found in these regions, which caused the change of the original direction of the raft structures and led to the formation of the 45-degree rafting.

Collapse and Topological Inversion
It is known that rafting usually occurs under mechanical loading at high temperature.Continuous creep deformation can cause the change of initial γ' precipitates of cubic shape to stable raft structures.However, further creep deformation may distort the stable raft structure and lead to topological inversion of γ/γ' phases, i.e., the γ' strengthening phase gradually surrounds the γ phase [94].Incorporating damage parameters [8,19,80] and the change of misfit strain [84,95] in the phasefield models allows for the determination of the major factors contributing to the collapse and topological inversion of the γ/γ' microstructures.
The damage parameters can represent the change of the elastic constants of the γ' phase and/or the increase of the plastic shear rate during creep deformation.Some slip systems can cut into the γ' phase and distort stable raft structures [8,80] when plastic strain accumulated near the γ/γ′ interface reaches a critical value.The results from the phase-field models with the variation of misfit strain reveal the decrease of misfit strain during the deformation, which makes the γ/γ' interface become semi-incoherent or incoherent.This trend destroys the stabilization mechanism in the γ channels and reduces the disjoining potential between adjacent γ' precipitates [95].This is an essential step to achieve topological inversion, as observed experimentally.
Currently, the phase-field models with damage parameters are mainly under the framework of elasto-plastic deformation.It is a challenge to fully describe the three creep stages for the tensile creep at elevated temperature in addition to the analysis of the collapse and topological inversion of the γ/γ' microstructures and the microstructure evolution.The incorporation of the change of misfit strain in the phase-field models also enables the observation of the topological inversion of the γ/γ' microstructures, and the simulation results are in good agreement with those observed in long-term aging experiments [95].

Summary
Phase-field method is an excellent technique to study the rafting process in nickel-based superalloys.Incorporating the theory of crystal plasticity in the phase-field method makes it possible to analyze the rafting kinetics with internal physical characteristics (the γ/γ' lattice misfit, elastic constants, and dislocation densities) as well as the mechanical behavior (heterogeneous elasticity and plastic deformation) of superalloys.Raft structures with N-type/P-type or other complex types can be formed during creep deformation under different loading modes.The creep deformation may lead to the collapse of raft structures and the topological inversion of the γ/γ' microstructures.
Introducing damage parameters or variation of lattice misfit in the phase-field models under elasto-plastic framework makes it possible to illustrate the individual processes for the microstructure evolution.The development of multiphase-field models provides the foundation to study the contribution of solute atoms to the rafting behavior and to better design and optimize superalloys for industrial applications.Future efforts are expected to focus on complex situations for industrial applications of single crystal superalloys, which likely include orientation-dependent, temperature-dependent, and composition-dependent mechanical properties.The evolution of the mechanical properties needs to be correlated to the microstructure evolution, such as the rafting.Additionally, the interactions between raft structures and defects, such as dislocations, micro-voids, micro-cracks, etc., need to be incorporated in the models.Further study is needed to expand the capability of the phase-field simulation in the analysis of the effects of voids and cracks of millimeter size.

Figure 1 .
Figure 1.SEM images of the typical microstructure of nickel-based single crystal superalloy (DD5): (a) initial microstructure without deformation, and (b) raft structure after the creep deformation at 980 °C.

Figure 2 .
Figure 2. Dislocations bowing through γ channels after long-term aging at 1050 °C for 100 h: (a) TEM image, and (b) schematic of dislocation bowing.In the field of view, some dislocations are present in the γ' phase, which is attributed to the dislocation climb [7].(Reproduced with permission from [8] © 2020 Elsevier.)

Figure 3 .
Figure 3. Types of raft structures under mechanical loading in <001> orientation of single crystals and positive lattice misfit.∆ represents the hardness ratio of the γ' precipitate to the γ matrix.Seven different cases are summarized and marked with corresponding numbers in the figure.1-N-type structure under compression for ∆ > 1; 2-P-type structure under tension for ∆ > 1; 3-P-type structure under compression for ∆ < 1; 4-N-type structure under tension for ∆ < 1; 5-N-type structure under compression for ∆ < 1; 6-N-type structure under tension for ∆ > 1; 7-P-type structure under tension for ∆ < 1. Adapted from [25].

Figure 6 .
Figure 6.Schematic for the migration of solute atoms under tension in <001> direction.

Figure 7 .
Figure 7. Deformation of a single crystal: (a) example of pure dislocation deformation F = F p , and (b) example of pure lattice deformation F = F e .

Figure 8 .
Figure 8. Flowchart showing the basic process to numerically solve the Ni-Al binary-field models and the multiphase-field models.

Figure 9 .
Figure 9. Numerical results of the microstructure and plastic strain at different instants during the tension of a nickel-based single crystal superalloy at a strain rate of 10 −3 per second and 1253 K: (a1-d1) microstructure, and (a2-d2) plastic strain.(Reproduced with permission from [8] © 2020 Elsevier.)

Figure 10 .
Figure 10.Microstructure evolution (top row) and the corresponding plastic strain field (bottom row) at different instants at 1253 K under 330 MPa from the phase-field simulation: (a) t = 0, (b) t = 5 h, (c) t = 16.83h, (d) t = 32 h, and (e) t = 51.39 h.Dark blue in the top row denotes γ matrix, and other four colors represent four different γ' variants.(Reproduced with permission from [8] © 2020 Elsevier.)

Table 2 .
Comparisons of different plasticity models used in the rafting analysis.