Dislocation density-based modeling of the yield drop phenomenon in nickel-based single crystal superalloy

The phenomenon of yield drop, characterized by a decrease in flow stress after initial yield, has been observed in various nickel-based superalloys. Despite numerous proposed physical mechanisms, there is still a lack of a meso-mechanism-based constitutive model to explain this phenomenon. In this study, the tensile behavior of a nickel-based single crystal superalloy (DDX), was investigated at different strain rates and a temperature of 900 °C. It was observed that the yield drop phenomenon in DDX became more pronounced with increasing strain rate. To predict the yield drop phenomenon during tensile processing, an improved strength law based on continuum dislocation density theory was considered in the crystal plasticity framework. The proposed constitutive model was implemented using nonlinear iteration and incorporated into a finite element analysis software. The simulation results exhibited a good agreement between the experimental data and the stress–strain curve in the vicinity of the yield drop region, affirming the predictive aptitude of the proposed model in elucidating the yield drop phenomenon at various strain rates.


Introduction
Superalloys, particularly nickel-based single crystal superalloys (NSCS), play a vital role in aeroengine and gas turbine applications [1][2][3].These alloys are subjected to demanding service conditions that necessitate high temperature strength, corrosion resistance, and oxidation resistance, ultimately ensuring prolonged service life.These desirable properties are attributed to the presence of a cubic strengthening phase, γ′ (Ni3(Al, Ti)), which constitutes a significant volume fraction (over 60%) within the matrix.The high temperature tensile behavior of various NSCS is governed by the volume fraction of the strengthening phase and the concentration of rare earth elements [4][5][6].At elevated temperatures, the motion of dislocations and the energy of phase interfaces are significantly influenced by the morphology and distribution of the γ′ strengthening phase [7].
Despite extensive research on the fundamental high temperature tensile properties of nickel-based single crystal superalloys (NSCS) over the past few decades [8][9][10][11][12][13][14][15][16], there is limited literature addressing the yield drop phenomenon specifically in the tensile process.He et al [9] investigated the tensile strength of different NSCS and derived a unified strength prediction law using the energy equivalence method.Kim et al [16] examined the tensile properties of an NSCS employed in turbine disks, but their proposed model is purely numerical and lacks practical physical significance.While the proposed model accurately predicts the tensile strength of NSCS, it fails to account for the yield drop phenomenon.Additional studies have explored the influence of crystal orientation [17][18][19], alloy composition [20][21][22], and heat treatment techniques [23][24][25] on NSCS behavior.
Understanding the underlying physical mechanisms of deformation during the tension process is of paramount importance for optimizing the performance of nickel-based single crystal superalloys (NSCS).The yield drop phenomenon, characterized by a decrease in flow stress after initial yield, has been observed in various NSCS and other metallic materials [26][27][28].Two theoretical frameworks, namely static theory [29,30] and dynamic theory [31,32], have been proposed to explain this phenomenon.According to the static theory, the presence of immobilized solute atoms creates locking effects that require high stress to overcome.On the other hand, the dynamic theory suggests that the yield drop phenomenon is a result of an increasing density of mobile dislocations following the macroscopic yield.Kim et al [26] argue that the dynamic evolution of dislocation density is a more appropriate explanation for the physical mechanism of the yield drop phenomenon.
While a formula provided in [26] can estimate the extent of the yield drop phenomenon, it does not offer guidance on how to incorporate this phenomenon into a constitutive model.However, it does confirm the correlation between the yield drop phenomenon and the dynamic evolution of dislocation density.This perspective encourages the adoption of a crystal plastic constitutive model coupled with the dynamic evolution of dislocation density to elucidate the yield drop phenomenon.The development of hardening models considering the thermal activation evolution of dislocations and dislocation substructure formation was initially proposed by Beyerlein [33].This constitutive model, which incorporates multiple slip and deformation twinning, has been employed to predict the anisotropic behavior of zirconium alloys under various loading paths and ambient temperatures.Over the past decade, the dislocation density-based constitutive model has demonstrated applicability to a range of metallic materials, including HCP materials such as zirconium [33], magnesium [34], and beryllium [35]; nickel-based single crystal superalloy [36][37][38][39]; FCC materials like aluminum alloys [40]; BCC materials such as tantalum-tungsten alloys [41]; and composite metal materials like Mg-Li-(Al) alloys [42] and Sn-Ag-Cu alloys [43][44][45].
Typically, the hardening law or slip resistance model incorporates various mechanisms, including initial slip resistance, forest dislocation term, dislocation substructure (debris), Hall-Petch effect term, and twin boundaries.These mechanisms cover the common hardening stages from stage I to stage IV [46].While these proposed constitutive models effectively account for different slip modes, non-Schmidt effect, anisotropy, and other factors, they do not consider the yield drop phenomenon.In the case of metallic glasses, the yield drop can be successfully demonstrated by introducing the concept of free volume concentration [47][48][49].Defects resulting from strain and thermal annihilation are interpreted as changes in the free volume concentration, and the equilibrium concentration term controls the occurrence of the yield drop [47].However, the visco-plasticity mechanism of metallic glasses differs from that of metallic materials, as the stress-strain response in metallic glasses exhibits strong thermally activated influences, and viscosity plays a critical role in shaping stress-strain curves.Metallic glasses are essentially amorphous metal alloys that exhibit mechanical behavior more akin to traditional glasses.Therefore, the free volume concentration model cannot be directly applied to nickel-based superalloys.
It is widely accepted that the yield drop phenomenon (YDP) is closely related to the evolution of dislocation density [50].Consequently, the dislocation density is influenced by factors such as forest dislocations and ambient temperature.In this study, we specifically focus on the dynamic evolution of forest dislocations as the primary deformation mechanism during the yield drop stage.To account for the yield drop, we incorporated an additional term representing forest dislocation into the slip resistance model.
The structure of this paper is outlined as follows: section 2 presents the results of tensile tests conducted on a nickel-based single crystal superalloy (NSCS) DDX at 900 °C under various strain rates.The fracture characteristics and microstructure features are analyzed using a scanning electron microscope.Section 3 provides a comprehensive review of crystal plasticity theory with consideration for finite deformation, elucidating the explicit forms of shear stress and shear strain increment on the slip system.Furthermore, a slip resistance model that incorporates the yield drop phenomenon is proposed.In section 4, the constitutive equation is derived and implemented, taking into account nonlinear iteration.Section 5 validates the proposed slip resistance model and constitutive model using the experimental data presented in section 2, also considering the key parameters within the slip resistance model.

Tensile test specimen of NSCS
In this study, a second-generation nickel-based single crystal superalloy (NSCS) DDX was utilized to investigate the YDP.Table 1 presents the composition of major chemical elements in DDX.The original DDX bar was produced using the vacuum directional solidification method with the selector technique.The intended solidification orientation was [001], although the actual cast bars may exhibit orientation deviations ranging from 0.5°to 13°.In this study, the orientation deviation of the bars was controlled within 5°.The standard heat treatment procedure applied to all original materials consisted of the following steps: 1290 [1].Dog boneshaped specimens were obtained by cutting DDX bars with an orientation along the [001] direction.The tensile test specimens had a length of 56mm, a gauge length of 12.5mm, and a gauge diameter of 4mm.The surface roughness of the specimens after turning and precision polishing was measured to be 0.4 μm. Figure 1 illustrates the scanning electron microscope (SEM) image of DDX after electrochemical corrosion, where the cubic strengthening phase (γ' phase) is surrounded by the matrix phase (γ phase).

Experiment and results
The tensile tests under constant strain rates were performed using a tensile testing machine, which can automatically record real-time stress and strain data through the controller.To maintain a consistent high temperature environment, a high temperature furnace was employed, ensuring an ambient temperature of 900 °C during the tensile tests.The temperature was monitored by three thermocouples in contact with the specimen's surface, with an error range of ±3 °C.Tensile tests were conducted at different strain rates: 10 −3 s −1 , 5 × 10 −4 s −1 , 10 −4 s −1 , and 10 −5 s −1 .
The results obtained from the tensile tests conducted at different strain rates are depicted in figure 2. It can be observed that the yield drop phenomenon (YDP) occurs at all strain rates when tested at 900 °C, and the extent of the YDP becomes more pronounced with increasing strain rate.The difference between the upper yield point (s y up ) and lower yield point (s y low )) is 195MPa for a strain rate of 10 −3 s −1 , while it decreases to 33.1MPa for a strain rate of 10 −5 s −1 .The fracture surfaces of the DDX specimens after tensile testing, as shown in figure 2, were examined using a metallurgical microscope.Due to the high temperature environment, severe oxidation is observed on the fracture surfaces.From the micrographs, it is evident that the fracture surfaces of the DDX specimens exhibit crystallographic characteristics, with the surfaces being parallel to at least one of the {111} planes.
The purpose of this section is to present the incorporation of the yield drop phenomenon into a crystallographic slip resistance model and its subsequent integration into the crystal plasticity finite element model (CPFEM).The forthcoming section outlines the framework of CPFEM, which is coupled with the crystallographic slip resistance model.

Kinematic and constitutive model
In this section, we provide a concise derivation of the constitutive equations that describe the elastic-viscoplastic deformation of single crystals under large strains.While these equations are widely available in the existing  literature, it is crucial to include them in order to incorporate the constitutive model involving yield drop phenomenon based on dislocation density, which will be discussed in the following section.
The complete deformation of a face-centered cubic (FCC) single crystal can be represented using the deformation gradient tensor F and the velocity gradient tensor L. These tensors provide a comprehensive description of the overall deformation of the crystal.
where the plastic deformation gradient, denoted as F , p represents the contribution of crystallographic slip to the overall deformation.On the other hand, the elastic deformation gradient, F , e characterizes the crystal rotation and lattice distortion.The elastic part of the velocity gradient tensor, denoted as L , e is calculated as the product of the time derivative of F e and the inverse of F .
e Similarly, the plastic part of the velocity gradient tensor, L , p is obtained by multiplying F , e the time derivative of F , p the inverse of F , p and the inverse of F .
e Following the deformation process, the sliding system experiences a change in the normal direction, represented as a n , and the sliding direction, denoted as a m .
The velocity gradient can be decomposed into two components: a symmetric part and an asymmetric part.
where the deformation rate tensor = + D D D , e p which is the sum of the elastic and plastic deformation rate tensors, can be expressed as follows: where the shear strain rate in slip α, denoted as g a ,  contributes to the formation of the spin tensor W = W + W . e p The spin tensor W = W + W , e p which is the sum of the elastic and plastic spin tensors, can be expressed as follows: The Kirchhoff stress t K is determined based on the following definition [51][52][53]: where the Jacobian matrix J, which quantifies the change in volume, is used in defining the Kirchhoff stress t .

K
Additionally, the Jaumann rate of Kirchhoff stress rate, t  , e K is associated with the symmetric rate of lattice stretching D e [51,53]: where the elastic stiffness tensor C, a fourth-order tensor, exhibits symmetric properties: K stress rate on axes that rotate with the material.
where the material derivative of Kirchhoff stress, denoted as t , K  can be obtained from equations (10) and (12).As a result, the Jaumann rate of t , K which undergoes rotation with the crystal lattice, can be expressed as follows: where the corotational Cauchy stress rate, denoted as s  , e represents the stress rate on axes that rotate with the crystal lattice.
Neglecting the volume change, with the assumption that only slip occurs from the initial configuration to the intermediate configuration, we have an approximation J≈1, leading to the following expression: : By combining equations (9), (15), and (16), we can express the difference between the two stress rates as follows: By combining equations (7), (17), and (18), we can express the Jaumann rate of Cauchy stress as follows: The resolved shear stress (RSS) on slip system a, taking into account finite elastic distortions, is defined as follows, as described in previous studies [51,54,55]: The inclusion of the additional terms a n* and a m* in equation (20) accounts for the influence of finite elastic distortions.To obtain the rate form of τ α , it is necessary to perform partial differentiation of equation (20): Through a series of mathematical deductions, we can express the rate form of t a as follows: :

 * *
To capture the aforementioned relationship, constitutive functions are required to define the shear rate g a  on slip system a.In the present study, A widely used power-law was used to describe the relationship between slip rate and the RSS on slip system a [ [43][44][45][56][57][58].
where g a 0  represents the reference slip rate, while a R c signifies the slip resistance, which will be elaborated on in the subsequent section.The rate sensitivity index, denoted as n, describes the rate sensitivity of materials.As n increases, the material exhibits a greater degree of rate independence.Furthermore, the function sign() is a symbolic function that maintains the same sign as the value enclosed in parentheses.Although it has been demonstrated that equation (23) can describe the deformation curves of many materials, Wu and Zaiser [59] have shown that the adoption of equation ( 23) can lead to a violation of the first and second laws of thermodynamics.
The resistance originating from the initial friction stress, denoted as a R , for is related to the initial critical resolved shear stress (CSS) of the slip system and the initial dislocation density [33].Additionally, a R HP 0, represents the initial slip resistance arising from the Hall-Petch effect [60].
The symbol μ represents the average shear modulus of the material, and it is commonly acknowledged that the effect of temperature on μ is more significant than its dependence on shear direction [33].The term f HP corresponds to the Hall-Petch coefficient, a b denotes the Burgers vector, and D g,0 represents the initial grain size.
The term a R for characterizes the resistance arising from the interaction between forest dislocations and mobile dislocations.Experimental evidence has demonstrated that a R for is proportional to the square root of the total forest dislocation density r a for [33,60].

R b 27
for for The dislocation interaction coefficient χ, which varies between 0.1 and 1.0, is set to a value of 0.9 in this study [33,60,64].The evolution of r a for can be understood as a competition between the increase in dislocation density resulting from dynamic loading and the decrease in dislocation density due to dynamic recovery processes [33].
for gen for rem for for for The parameter a k 1 represents the forest dislocation density, which is independent of the deformation rate and represents the accumulation of dislocations due to the statistical interaction between mobile dislocations and forest obstacles [61].On the other hand, the coefficient ( ) is rate-sensitive and dependent on the deformation rate and temperature.This coefficient is associated with dislocation dynamic recovery and accounts for processes such as dislocation cross-slip and climb through thermal activation [60].The following form is commonly recommended for expressing this relationship: The activation energy a Q nor is independent of stress, and it is multiplied by the Boltzmann constant k. a D represents the drag stress, and e = s 10 0 7 1


is the reference strain rate.Equation (29) is valid for strain rates ranging from - 10 5 to /s 10 . 3 The term a R deb in equation (24) represents the contribution of dislocation substructures to the slip resistance.This term is commonly associated with stage IV hardening.Madec et al [65] introduced the following formula to determine a R : The expression for a R deb in equation (30) involves a material constant k , deb which is commonly assigned a value of 0.086.r deb represents the dislocation density of the substructure and can be determined using functions that account for the dynamic recovery of all active dislocations: The rate coefficient a q is associated with the process of debris formation through local thermal activation.The final term in equation (24) accounts for the reduction in slip resistance that gives rise to the phenomenon of yield drop.The origin of yield drop can be attributed to the Johnston-Gilman mechanism [28], which assumes that the yield drop phenomenon occurs when the primary density of mobile dislocations is too low [27].The primary dislocations are immobilized by other dislocations and obstacles, resulting in a rapid increase of forest dislocation density before the yield drop occurs.When applied stress reaches a relatively higher level, the dislocations can break away from their pinning point and generate large amounts of mobile dislocations, causing the stress to drop suddenly [27].In addition, upon analyzing the experimental results, it can be found that the yield drop occurs when the strength of the specimens reaches certain values.This suggests that the yield drop phenomenon could potentially be a consequence of the evolution in dislocation density.Therefore, it is assumed that a R dro in equation ( 24) evolves with the forest dislocation density.
where the payield rameter Θ quantifies the extent of the yield drop phenomenon, while r a for denotes a reference dislocation density that is influenced by the strain rate.The value of r a for governs the rate at which the yield drop occurs.Additionally, κ is an index that accounts for the impact of forest dislocation density evolution on the yield drop phenomenon.
By combining equations ( 24) and (32), the comprehensive formulation for the hardening law can be expressed as follows:

Numerical implementation
In this section, an algorithm is developed to implement the proposed dislocation density-based slip resistance model coupled with the crystal plasticity model.To calculate the nonlinear shear strain increment on the slip system for a given strain increment, the Newton-Raphson (N-R) iterative process is employed.However, the initial value of the iteration is determined using the exhibit method.The stresses and state variables are then evaluated at the end of the time increment using an implicit time integration scheme.

Determination of initial value of iteration
At a given time t and time t+Δt, the strain increment ∆g a of the slip system can be determined through linear interpolation.
The parameter q varies between 0 and 1.By performing a first-order Taylor expansion of  we obtain: From equations (34) and (35) gives: In equation (36), the partial derivative of g a t  is easy to determine: From equation (22), the incremental form of t a  is given as: To incorporate the proposed hardening model, it is necessary to provide an explicit expression for the incremental form of a R .c Since a R c 0, remains invariant according to equation (33 Differentiating a R c with respect to the corresponding variables yields: where the ab H , for ab H deb and ab H dro are the hardening modulus due to the forest dislocation evolution, substructure evolution and yield drop, i.e.
When a b = : When a b ¹ : From equations (34) to (44), the iterative initial value of shear strain increment on slip system is determined by the following formula:

Implicit solution of N-R iteration
As the initial value of ∆g a is obtained through the first-order Taylor expansion, it provides an approximate solution with some deviation.To achieve a more accurate solution, the Newton-Raphson iterative method is employed to solve the nonlinear incremental equation (10).


The equation mentioned above is derived from the equation specific to the α slip system, but it is dependent on the strain increment ∆g b for all slip systems (where β ranges from 1 to n).Hence, equation (15) forms a system of n nonlinear equations with n unknowns.To solve this system, the Newton-Raphson iterative method is employed in the following form: In the equation provided, ∆g a + t k , 1 represents the result of the kth calculation, ∆(∆ ) g a + t k , 1 denotes the increment computed through the Newton-Raphson iterative process, and (∆ ) , with respect to the strain increment.Referring to equation (46), we have: The partial derivative of ∆ a R c to ∆g b is given by equation (40):  From equation ( 49) to (54), the differential operator of (∆ ) g a F t k , with respect to strain increment is given by: The solution for ∆g b is obtained using the Newton-Raphson iterative method, as described in equation (47), with the initial value obtained from the exhibit solution given by equation (45).The remaining increments of the state variables are determined using the same iterative procedure.

Simulation and discussion
The enhanced crystal plasticity constitutive model was implemented as a user-defined material subroutine.A cubic finite element model (FEM) was employed for the simulations, with the C3D8R element chosen as the material element.To prevent hourglass modes, the enhanced hourglass stiffness control method was applied.The boundary conditions for the FEM are depicted in figure 3. The fixed point (O) was selected as the origin of the FEM, with no displacement in any direction (dx = dy = dz = 0).The following boundary conditions are imposed to prevent the surface ADBO from rotating about point O, but not to limit the stretching of the individual edges.Constraints were applied to three edges connected to point O, with specific displacement constraints set for each edge (OA: dy = dz = 0; OB: dx = dy = 0; OC: dx = dz =0).Similarly, constraints were applied to the edges of CF (dy = dz = 0) and CE (dx = dy = 0) only to prevent rotation of the surface CEFG.At point D in figure 3, the displacement in the y direction was fixed.A displacement of 7% was applied to the surface CEFG.The duration of the analysis step was determined based on different strain rates (70s, 350s, 700s, 7000s).
The integration parameter q and the reference slip rate g a 0  are set to constant values based on the literature sources [43-45, 51, 52].The mechanical properties of the material, including E, v and μ, are adopted from [66].The reference strain rate e 0  is set to 10 7 s −1 following the work of Beyerlein et al [33].The dislocation interaction coefficient c is chosen as the recommended value of 0.9, as suggested by previous studies [33,60,64].It is known that at the same temperature, the extent of yield drop is closely associated with the strain rate.Hence, the dislocation interaction index k remains constant for different strain rates, as it plays a role in controlling the rate of yield drop. Increasing k results in a narrower range of plastic strain over which yield drop occurs.The stressindependent activation energy a Q nor can be determined by fitting experimental data obtained at various strain rates [44,45].The drag stress a D is a proportionality factor that shares the same units as stress [33].Its value is determined through fitting procedures based on experimental data obtained during the hardening stage.The Burgers vector a b and the initial dislocation density r a for,0 for DDX with [001] orientation can be found in reference [66].The trapping rate parameter a k 1 is typically obtained from the slope of the initial stress-strain curve following the yield point [60].However, for the materials investigated in this research, determining a k 1 using the initial slope method is not feasible due to the presence of the yield drop phenomenon.Therefore, in this study, the value of a k 1 is determined by fitting the initial slope of the stress-strain curve after the lower yield point.
The relevant parameters utilized in this crystal plasticity finite element model (CPFEM) are provided in table 2.
The comparison between the simulation results and experimental data is presented in figure 4, demonstrating the capability of the proposed model to accurately predict the stress-strain curves and yield drop phenomenon of DDX at different strain rates.The accurate prediction of the yield drop is primarily attributed to the formulation in equation (32), where the strain rate-dependent parameter ( ) r e a ref  governs the extent of the yield drop and its determination requires fitting the yield drop range at various strain rates.The relationship between ( ) r e a ref  and the applied strain rate e  is depicted in figure 5.It is evident from figure 5 that a strong logarithmic linear relationship exists between the reference dislocation density and the applied strain rate:   The net hardening rate proposed by Mecking et al [46] can be used to examine the change of hardening modulus in the process of yield drop stage, it can be defined as the gradient of macro stress-strain curve as a generation of initial plastic: When the parameter q H approaches zero, a significant range of plastic flow occurs under nearly constant stress, indicating a stable state of plastic flow.The influence of the reference dislocation density r a ref and yield drop strength Q on the net hardening rate is presented in figures 6(a) and (b), respectively.As mentioned earlier, the value of r a ref determines the occurrence range of the yield drop phenomenon.With an increase in r a , ref the evolution of dislocation density is accelerated, leading to a narrower range of strain where the yield drop occurs (figure 6(a)).The yield drop strength Θ affects the magnitude of the yield drop, but it also impacts the range of the yield drop phenomenon (figure 6(b)).Therefore, the determination of Θ needs to be considered in conjunction with the reference dislocation density to ensure accurate modeling.

Conclusion
This study investigates the rate-dependent tensile behavior of nickel-based superalloy DDX at 900 °C, with a particular focus on the observed yield drop phenomenon across different strain rates ranging from 10 −5 s −1 to 10 −3 s −1 .To accurately simulate this phenomenon, an enhanced slip resistance model incorporating dislocation density is proposed and integrated into a finite deformation crystal plasticity framework.The performance of the model is verified and analyzed by comparing the simulation results with experimental data.The main conclusions drawn from this research are as follows: 1. Experimental observations reveal that the yield drop phenomenon in nickel-based superalloy DDX becomes more pronounced as the strain rate increases.Higher strain rates lead to wider yield drop ranges following the initial yield point.
2. The developed slip resistance model, which considers dislocation density, is successfully incorporated into the finite element program, enabling nonlinear iterations within the framework.
3. The simulation results of the tensile stress-strain curve exhibit excellent agreement with experimental data, indicating the model's strong predictive capability for different strain rates of nickel-based superalloy DDX.
4. The proposed reference dislocation density governs the extent of the yield drop phenomenon, while the yield drop strength determines its magnitude.Both parameters collectively describe the entire process of the yield drop phenomenon.

Figure 1 .
Figure 1.The scanning electron microscope of DDX after electrochemical corrosion.

Figure 2 .
Figure 2. True stress-strain curve of tensile test with different strain rate and the corresponding fracture surface of DDX at 900 °C.

Figure 3 .
Figure 3. Schematic of the boundary conditions used in the finite element model.

Figure 4 .
Figure 4.The uniaxial tensile true stress-strain curve of DDX is depicted in both experimental data and simulation results.The test data are represented by points with different shapes, while the solid lines correspond to the simulation results.

Figure 5 .
Figure 5.The relation between reference dislocation density and strain rate.

Table 2 .
The relevant parameters used in this CPFEM.