Crystal Plasticity Model Analysis of the Effect of Short-Range Order on Strength-Plasticity of Medium Entropy Alloys

: Numerous studies have demonstrated the widespread presence of chemical short-range order (SRO) in medium and high entropy alloys (M/HEAs). However, the mechanism of their inﬂuence on macroscopic mechanical behavior remains to be understood. In this paper, we propose a novel dislocation-based model of crystal plasticity, by considering both the dislocation blocking and coplanar slip induced by SRO. The effect of SRO on the plastic deformation of CoCrNi MEAs was investigated. We found that the yield strength increases monotonically with increasing SRO-induced slip resistance, but the elongation ﬁrst appeared to increase and then decreased. Further analysis suggested that the plastic elongation is a result of the competition between grain rotation-induced deformation coordination and stress concentration, which depends on the slip resistance of the SRO.


Introduction
Due to their remarkable mechanical properties, medium and high entropy alloys (M/HEAs), which include multiple primary elements in almost equal amounts, have recently received a lot of attention [1][2][3][4][5][6][7][8][9][10][11][12][13]. M/HEAs were initially conceived as random solid solutions, with a perfect chemical disorder organization. The mismatched atomic radii and intricate interactions between the constituent elements, however, provided evidence that the arrangement of atoms in M/HEAs is not entirely random [14][15][16]. Chemical elements can be rearranged to lower their free energy during solidification or/and heat treatment operations. This results in the formation of specific preferential inter-atomic neighborhoods and the production of chemical short-range orders (SRO) [17][18][19][20]. Accordingly, the presence of SROs may be a general, yet distinctive, structural characteristic of M/HEAs, and which can be used to adjust properties even more effectively.
From many atomic-level simulations and experimental studies, it is known that the SRO interacts with dislocations, the SRO impedes dislocation motion, and when sufficient dislocations pass through an SRO region (severe plastic deformation), the order degree of the SRO region decreases (towards a random solid solution) [18,[21][22][23]. In addition, the increasing degree of SRO in an alloy material leads to a transformation of the dislocation slip, from a wavy to coplanar configuration [19,[24][25][26]. Although in different ways, all these results highlight the strong influence of SRO on mechanical properties. However, due to the chemical and structural complexity of M/HEAs, understanding their full impact on mechanical behavior is challenging. Further efforts are needed, to characterize the extent and macroscopic processes that influence the mechanical behavior of M/HEAs through the SRO mechanism.
In this paper, a modified crystal plasticity finite element model is developed to analyze the influence of SRO on deformation behaviors, in which the effects of both dislocation impedance and coplanar slip by SRO are considered. The typical CoCrNi MEAs, in which SRO is well observed in both experiments and simulations, was taken as the model alloy [4,19,22]. The modified model was used to simulate a tensile experiment of the polycrystalline CoCrNi MEAs and confirm its corresponding parameters. On this basis, the effect of deformation is discussed, by adjusting the degree of SRO.

Kinematics of Crystal Plasticity
The total deformation gradient tensor F (bold font indicates second-order tensor) can be decomposed as: where the superscript * denotes lattice deformation and superscript P denotes plasticity. The total velocity gradient is defined as [27]: The total velocity gradient can be decomposed into the plasticity part L P and the lattice deformation part L * . The strain rate D * is the symmetric part of L * . Using the Jaumann objective stress rate for Kirchhoff stress, the instanton equation is expressed as: where the stiffness matrix is considered to be a linear combination of the stiffness matrix of the matrix and the twin, Kirchhoff stress τ = Jσ, J = det(F), σ is the Cauchy stress, f β is the total volume of the twin, f β (regular italic font indicates scalar) and C β tw (Euclid math font indicates fourth-order tensor) are the volume fraction and the stiffness of twins for the twin system β, respectively. C m is the stiffness of the matrix material. N tw is the total number of twin systems. The velocity gradient of plasticity consists of dislocation slip and twinning: where N s is the total number of slip systems.
. γ α and . γ β are the individual shear strain rates of the slip system α and twin system β. m α and n α (bold italic font indicates vector (first-order tensor)) are the unit vectors along the shear and normal directions of the slip plane, while m β tw and n β tw are those of the twin plane.

Dislocation Slipping
In a dislocation-based model, the Orowan equation [28] gives the shear rate on the slip system α as: where ρ α and τ α e f f are the dislocation density and the effective resolved shear stress on the slip system α, τ sol is the solid solution strength, b is the length of the Burgers vector for slip, v 0 is the dislocation glide velocity, Q s is the activation energy for dislocation slip, k B is the Boltzmann constant, T is the temperature, and p and q are fitting parameters.
The effective shear stress τ α e f f is modified on the basis of Xiaochong Lu [29]: where τ α is the corresponding resolved shear stress of the slip system α, τ α b is the back stress, τ α p is caused by the obstacles of forest dislocation and Hall-Petch stress, and τ α SRO is dislocation resistance of SRO. The regions of SRO can be thought of as a smaller-sized and more densely distributed eutectic lattice precipitation [16,19], which has a certain pegging effect on the dislocation motion and can improve the yield strength at the macroscopic level [22]. Therefore, adding a negative term to the effective stress for the resistance of SRO increases the yield strength and decreases the slip velocity in the slip system.
The first term τ α in the right of Equation (6) is the resolved shear stress is [27]: The second term τ α b in the right of Equation (6) is the back stress, which is caused by the dislocation pile-ups at the grain and twin boundaries [28]: .
where G is the shear modulus, Λ α s is the geometrical length scale of the microstructure, l is the average distance between slip bands, N α is the number of piled-up dislocations at boundaries, and N * is the saturated number of N α .
The microstructural geometrical length Λ α s in Equation (8) is defined by the combined effect of the grain and twin boundaries [29]: 1 λ α where d grain is the grain size, ξ αβ characterizes the interaction between the slip system α and twin system β, and t is the average thickness of twin lamellas. The passing stress τ α p in Equation (6) is caused by the obstacles of forest dislocation and Hall-Petch stress [29]: where σ 0 , k HP are coefficients of the Hall-Petch law, and ξ αα is the interaction coefficient between the slip system α and α , including the self-hardening, coplanar, collinear, orthogonal, glissile, and sessile interactions [30]. According to the modified Kocks-Mecking rule, the evolution rate of dislocation density can be expressed by [31]: where k is the forest dislocation hardening constant, and d anni is the annihilation distance for dislocations. In order to describe the dislocation blocking and coplanar slip induced by SRO, τ α SRO can be decomposed into two terms: one for the initial resistance τ 0 SRO , and one for the evolution coefficient D α SRO . The evolution coefficient D α SRO decreases rapidly with cumulative slip, and since the coefficient D α SRO should tend to zero, as the cumulative slip tends to infinity, an exponential form is used. The specific form is as follows: where τ 0 SRO is the initial value of τ α SRO . γ 0 is the reference strain. ξ αη is the dislocation resistance coefficient of SRO between the slip system α and η, including coplanar and noncoplanar interactions.

Deformation Twinning
Twining causes plastic deformation and serves as a grain-refining process. The twinning evolution adopts a similar form to that of Xiaochong Lu [29]. The equation of twinning evolution is briefly explained in this section. The twin system's contribution to the shear rate can be represented as follows: where γ tw is the characteristic twinning shear strain as √ 2/2 for FCC (Face Center Cubic) lattice. Λ β t is the mean spacing between two obstacles caused by deformation twins. w is the available width of twin lamellas, when w > Λ β t , w = Λ β t . The addition of w results from the fact that twins usually do not penetrate the entire grain when the grain size is large.
The equation of Λ β t is similar to Equations (10) and (11): where ξ ββ is the interaction coefficient between the twin systems β and β . The twin nucleation rate . N β t in Equation (15) is influenced by the dislocation cross-slip and the resolved shear stresses on twin systems: N 0 p ncs p tn (18) where . N 0 is the reference twin nucleation rate. p tn is the probability of forming twins. The probability that no cross-slip occurs p ncs is calculated by: where V cs is the activation volume for dislocation cross-slip, τ r is the stress needed to bring two partial dislocations from the equilibrium distance x 0 to the critical distance x c in which the twin nucleation can be facilitated. τ r is expressed as: The equilibrium distance x 0 between two partial dislocations is determined by: where Γ s f is the value of stacking fault energy and ν is the Poisson ratio. b p is the length of the Burgers vector of Shockley partial dislocation, and τ β is the resolved shear stress on the twin system β: The resolved shear stress on a twin plane affects the twin nucleation rate. The probability of forming a twin is calculated by: where A is a fitting parameter. τ t is the critical stress for twinning, equal to the activation stress for a Frank-Read dislocation source plus the critical stress for generating a wide stacking fault: In summary, using a dislocation-based crystal plasticity model that considers back stress, forest dislocation hardening, and twinning, the model introduces the impacts of SRO. The effect of the SRO impeding dislocation motion is reflected by adding τ α SRO to the Orowan equation (Equations (5) and (6)), τ α SRO decreases τ α e f f , causing the dislocation slip velocity to decrease. Since the impeding effect of SRO varies with the slip of the material [22], the evolutionary coefficient D α SRO of τ α SRO is given. Meanwhile, the interaction coefficient of the coplanar slip system is positive and the interaction coefficient of the non-coplanar slip system is negative, making the coplanar slip occurs easily, while the non-coplanar slip is suppressed.

Polycrystalline Finite Element Model
Simulation of polycrystalline CoCrNi MEAs used ABAQUS software (2020, Dassault Systemes Simulia Crop., Johnston, RI, USA). The material model of Section 2 was implemented using the Umat subroutine. A bone-shaped flat tensile sample is shown in Figure 1a. To reduce the calculation, an equivalent polycrystalline model was used for the clamping end (white part), and a crystal plasticity model was used for the samples in the parallel part, with different colors representing different grains (Figure 1a). The C3D8 (8-node linear brick) element was used in the simulation, with an element size of~1 µm. The parallel segment size was 16 × 40 × 100 µm 3 and contained a total of 55 grains (55 random points in parallel segment space generate 55 Voronoi polyhedra), so the grain size was about 13 µm.
Three groups of terms used in the paper are clarified prior to the analysis. The first group consists of macroscopic strain and local strain; the macro-strain is defined as the difference in displacement below and above the parallel segment divided by the length of the parallel segment, denoted by E; and the local strain is defined as integral point strain, denoted by ε. The second group consists of macroscopic stress and local stress; the macroscopic stress is equal to the sum of the restraint reaction forces divided by the cross-sectional area of the parallel section; and local stress is the stress of the integral point. The third group is the elongation, which is the macro-strain that corresponds to the highest point on the macro-stress. vs. macro-strain curve. difference in displacement below and above the parallel segment divided by the length of the parallel segment, denoted by E ; and the local strain is defined as integral point strain, denoted by  . The second group consists of macroscopic stress and local stress; the macroscopic stress is equal to the sum of the restraint reaction forces divided by the crosssectional area of the parallel section; and local stress is the stress of the integral point. The third group is the elongation, which is the macro-strain that corresponds to the highest point on the macro-stress. vs. macro-strain curve.

Parameter Validation
To verify the validation of the developed constitutive model, the tensile responses of CoCrNi MEAs [32] were simulated in this section. The material parameters can be divided into three categories: the first category is the parameters that can be clarified, such as the lattice parameter, stacking fault energy [33], elastic properties [34], and dislocation activation energy [35]. The second category is parameters with ambiguous reference values, such as the average distance between slip bands, interaction coefficient between the twin systems, solid solution stress, the thickness of twin lamellas, etc. The values of these parameters were adjusted based on the parameters of TWIP (twinning induced plasticity) steel [29], which is similar to CoCrNi MEAs. Hall-Petch coefficients obtained in tensile experiments needed to be divided by the Taylor factor and adjusted appropriately according to the simulation results. The third category is the fitting parameters, such as, p and q Equation (5), A in Equation (23). In addition, the parameters related to SRO were also found by fitting with the experiment data. The specific parameters are shown in Table 1,

Parameter Validation
To verify the validation of the developed constitutive model, the tensile responses of CoCrNi MEAs [32] were simulated in this section. The material parameters can be divided into three categories: the first category is the parameters that can be clarified, such as the lattice parameter, stacking fault energy [33], elastic properties [34], and dislocation activation energy [35]. The second category is parameters with ambiguous reference values, such as the average distance between slip bands, interaction coefficient between the twin systems, solid solution stress, the thickness of twin lamellas, etc. The values of these parameters were adjusted based on the parameters of TWIP (twinning induced plasticity) steel [29], which is similar to CoCrNi MEAs. Hall-Petch coefficients obtained in tensile experiments needed to be divided by the Taylor factor and adjusted appropriately according to the simulation results. The third category is the fitting parameters, such as, p and q Equation (5), A in Equation (23). In addition, the parameters related to SRO were also found by fitting with the experiment data. The specific parameters are shown in Table 1, and the fitting results is the red curve in Figure 1b, which fits well with the experimental results of Miao et al. [32].
To further confirm the feasibility of these parameters, the feature length of the model was adjusted from 1 µm to 69 µm and 77 µm (the average grain size from 13 µm to 900 µm and 1000 µm) to compare with the results in Zhang R and Zhao S [19] (the true stress-strain curve was transformed into engineering stress-strain curves, while taking the strain rate as 4 × 10 −3 s −1 ). The blue curve corresponds to the curve of low SRO (τ 0 SRO = 10 MPa), and the purple curve is the curve of high SRO (τ 0 SRO = 50 MPa) (Figure 1b). These calculated results and the experiments match well, which shows the applicability of the parameters to CoCrNi MEAs.

Influence of SRO on Deformation Behavior
Based on the model developed in Section 2 and the parameters in Section 3, the effect of the resistance of SRO on the deformation behavior was evaluated by adjusting the value τ 0 SRO . Figure 2a,b show the macro-stress vs. macro-strain curves at different τ 0 SRO . By observing the curves, it can be seen that the yield strength increases monotonically when τ 0 SRO increases (Figure 2(a 1 ,b 1 )). Furthermore, when observing the end of the curves, the elongation of the material rises and then falls as τ 0 SRO increases (Figure 2(a 2 ,b 2 )). The observed phenomenon is plotted in Figure 2c. From Figure 2c, when τ 0 SRO ≈ 29 MPa, the elongation takes the maximum value.  To analyze why the elongation increases and then decreases with increasing τ 0 SRO , we plotted the pole figures at τ 0 SRO = 0, 30, 60 MPa and macro-strains E 22 = 0.1, 0.3 ( Figure 3). Figure 3a shows the distribution of 55 grains before deformation. As the micro-strain increases, the crystallographic orientation on the integration points unfold on the pole figures (Figure 3b-d). When E 22 = 0.1, the orientations of the integration points from grain #2 are dispersed on the pole figures (the red points in (Figure 3(b 1 ,c 1 ,d 1 )). When E 22 = 0.1 and τ 0 SRO = 0 MPa, the integration points from grain #2 are scattered in the region ∆l 0 = 0.077 and ∆θ 0 = 11.17 • (Figure 3(b 1 )). Similarly, when E 22 = 0.1 and τ 0 SRO = 30 MPa, the integration points are scattered in ∆l 30 = 0.083 and ∆θ 30 = 12.78 • (Figure 3(c 1 )). Moreover, when E 22 = 0.1 and τ 0 SRO = 60 MPa, ∆l 60 = 0.099, ∆θ 60 = 13.12 • (Figure 3(d 1 )). The crystallographic orientation of the material changes drastically as τ 0 SRO increases. It can be further inferred that the increase in elongation with increasing τ 0 SRO , when τ 0 SRO < 29 MPa, is due to the intense local rotation, which can provide additional macro-strain.  To illustrate the effect of local crystal rotation on macroscopic deformation, we chose a single crystal with ϕ 1 ≈ 9.74 • , Φ = 45 • , and ϕ 2 = 0 • , for the simulation (Figure 4a). The grain was in the following orientation: the slip plane {111} was 45 • to the tensile direction, and the slip plane {111} was easily activated, with a clear coplanar slip effect under the effect of SRO. As shown in (Figure 4b), the blue curve (τ 0 SRO = 0 MPa) and the red curve (τ 0 SRO = 30 MPa), both below the diagonal line (the black), indicate that the macro-strain was not only contributed to by the local (element) deformation, the rotation of the elements also provided a certain amount of macro-strain. When the macro-strains    As seen from the results in Figure 4, the local rotation of the material becomes more intense when τ 0 SRO is larger. In polycrystalline materials, due to the different rotational directions generated by different grains, a more intense rotation causes a stronger stress concentration ( Figure 5(a 1 ,b 1 ,c 1 )). From Equation (12), it is known that high local stress implies a high dislocation density; thus, quickly attaining a saturated dislocation density locally and causing the material to lose its capacity to harden ( Figure 5(a 2 ,b 2 ,c 2 )).
Metals 2022, 12, x FOR PEER REVIEW 11 of 13 As seen from the results in Figure 4, the local rotation of the material becomes more intense when 0 SRO  is larger. In polycrystalline materials, due to the different rotational directions generated by different grains, a more intense rotation causes a stronger stress concentration ( Figure 5(a1,b1,c1)). From Equation (12), it is known that high local stress implies a high dislocation density; thus, quickly attaining a saturated dislocation density locally and causing the material to lose its capacity to harden ( Figure 5(a2,b2,c2)). The impact of SRO on the material's macroscopic deformation was explored using the developed model. The yield strength of the material rises monotonically as the amount increases, yet the elongation exhibits a trend of rising and then falling. The rise in dislocation resistance is the cause of the rise in yield strength. The intense grain rotation brought The impact of SRO on the material's macroscopic deformation was explored using the developed model. The yield strength of the material rises monotonically as the amount increases, yet the elongation exhibits a trend of rising and then falling. The rise in dislocation resistance is the cause of the rise in yield strength. The intense grain rotation brought on by coplanar slip is the source of the initial rise and subsequent fall in elongation. When the τ 0 SRO < 29 MPa, the local grain rotation causes additional macro strain; while when the τ 0 SRO > 29 MPa, the intense rotation causes severe stress concentration, resulting in a decrease in elongation. There is a competitive relationship between the contribution of grain rotation to elongation and the decrease in elongation due to the stress concentration. When the resistance of SRO is small, this contribution dominates, and when the resistance is large, the degradation of stress concentration dominates. This model is a good explanation for the fact that the SRO causes non-uniform dislocation motion, which leads to plasticity improvement. On the other hand, with the increase of the order degree or density of SRO, the dislocation plug is aggravated during the deformation, which reduces the deformation ability of M/HEAs [16,19]. The monotonic increase in strength predicted by the model is consistent with the experimental results of Zhang R et al. [19], Ding et al. [16], and D. Han [24]. However, the results of elongation require more precise control of the SRO in experiments. The results of the simulation of Zhang R et al. [19], to a certain extent, show that the elongation decreases when the SRO is outside of the appropriate range.

Conclusions
Through the incorporation of the resistance of SRO into the dislocations and the coplanar slip effect into the dislocation-based model, this study offers a crystal plasticity model that is better suited to the M/HEAs. The main conclusions are as follows: (1) A set of parameters consistent with CoCrNi MEAs was determined and can be used to discuss the influence of various factors on a material's deformation behavior. (2) Adjusting the resistance of SRO at a certain range increases both the yield strength and elongation simultaneously, but beyond this range, the yield strength increases but the elongation decreases. (3) As the resistance of SRO increases, the elongation increases and then decreases, which is attributed to the more intense local rotation with coplanar slip. Local rotation can increase the additional macro strain, while also causing a more intense stress concentration; when the resistance of SRO is low, the additional macro strain dominates the elongation increase; when the resistance is high, the stress concentration dominates the elongation decrease.
The results predicted by the model indicate that the effect of the SRO on the uniform deformation of the M/HEAs is not monotonic; therefore, in the subsequent material design and treatment, whether the SRO is regulated by annealing treatment, or by changing or adding elements, the SRO should be controlled within a suitable range, to achieve a simultaneous improvement of strength and plastic deformation.  Data Availability Statement: Available from the corresponding author upon request.

Conflicts of Interest:
The authors declare no conflict of interest.