Atomistic Investigation on the Strengthening Mechanism of Single Crystal Ni-Based Superalloy under Complex Stress States

: Single crystal Ni-based superalloy, with excellent mechanical properties in high temperature, always works under complex stress states, including multiaxial tension and compression, which results in various strengthening mechanisms. In this paper, the atomistic simulation is applied to investigate the microstructure evolution under complex mechanical loading conditions, including uniaxial, equibiaxial, and non-equibiaxial tensile–compressive loadings. By comparison of the strain– stress curves and analysis of dislocation motion, it is believed that the tension promotes the bowing out of dislocations into the channel at loading direction, while compression limits it. Moreover, the dislocation analysis shows that the initial dislocation network, comprised of Lomer dislocations, will dissociate to form Lomer–Cottrell lock upon loading, which acts as a barrier to the further glide of dislocations. The mechanism of dislocation evolution is analyzed in detail by combining Schmid factor analysis and the comparison of energy density difference between γ and γ (cid:48) phases.


Introduction
Single crystal (SC) Ni-based superalloy, known for its high strength, oxidation resistance, corrosion resistance, and creep resistance under high temperature, is a vital aerospace material, widely used in crucial components, such as aero-engine turbine blades, rocket engine nozzles, and aerospace aircraft thermal coatings [1,2]. With the development of material manufacture technology, the current SC Ni-based superalloy is composed of the face-centered cubic (FCC) Ni (γ) matrix phase and cubic morphology ordered L1 2 Ni 3 Al precipitates embedded in γ phase, namely γ' phase, at as-received state [3]. The special microstructure of Ni-based superalloy requires a unique strengthening mechanism, which has a great influence on the mechanical property of the alloy during service [4][5][6].
In the past few decades, numerous experimental tests have been carried out to reveal the relation between the microstructure and mechanical property of the Ni-based superalloy [7][8][9]. Luo et al. [10] studied the variation of the lattice parameter with the temperature of SC Ni-based superalloy and found that the lattice misfit between γ and γ phases increases with temperatures below 1100 • C. Zhang et al. [11] and Chang et al. [12] reported that the interfacial dislocation networks can prevent the dislocation motion effectively, which leads to the prolonging of the creep rupture time. Pei et al. [13] investigated the evolution and influence mechanism on creep failure of the oxide layer, heat-affected layer, and substrate microstructure. They showed that creep properties are seriously affected by pre-oxidation. Ding et al. [14] displayed the mechanical behavior of SC Ni-based superalloy by tensile testing at various temperatures and revealed the deformation substructures and mechanisms from micrometer to atomic scales. In recent decades, plentiful worthwhile discoveries have been published by experiment researchers. However, due to the limitation of the experimental method and complexity of loading conditions, the in situ observations are difficult to conduct.
On the other hand, molecular dynamic (MD) simulation can capture dynamic deformation in real time and provide atomistic details, such as the interaction between dislocation and γ /γ structure, atomistic stress, and energy distribution [15,16]. Zhu et al. [16] and Xiong et al. [17] performed the matrix dislocation motion in SC Ni-based superalloy and provided the interplay mechanism between the matrix dislocation and interfacial dislocation network by MD simulation. Wu et al. [18][19][20] established 2D and 3D MD models to observe the interfacial dislocation network evolution under uniaxial loading. Most of these MD investigations regarding the Ni-based superalloy adopt a uniaxial loading and overlook the fact that the multiaxial loading will result in a considerable difference in the microstructure evolution of SC Ni-based superalloy during the MD simulation.
In general, the SC Ni-based superalloy always works in a complex mechanical loading environment, resulting in numerous degeneration behaviors, such as rafting and precipitation of TCP phase [7]. Gebura and Lapin [21] studied the effect of multiaxial stress state, which is caused by a notch effect, on the microstructure degradation of CMSX-4 and showed that notches affect the degradation process of alloy significantly. Caccuri et al. [22] provided the orientation and type of rafting that are governed by the stress type and triaxiality in the initial and middle states during service of SC Ni 3 Al superalloys with a small strain. Cao et al. [23] investigated how the evolving multiaxial stress states affect the kinetics of rafting during creep of SC Ni-based superalloy. Their results supported the fact that the maximum principal stress is the key component of the overall stress state, which governs rafting.
Meanwhile, due to the L1 2 atomistic structure of γ phase, the SC Ni-based superalloy clearly illustrates the tension-compression asymmetry, which leads the compression loading to play an important role in the study of complex stress conditions [7,19,[24][25][26][27]. The tension-compression asymmetry of yield and creep strengths was studied by Tsuno et al. [28], and they believed that the twin boundary dominated the tension-compression asymmetry mechanism. Wang et al. [29] examined the tension-compression asymmetry of SC Ni-based superalloy DD6 during low cycle fatigue. They found that the tensile-compressive stresses are related to the dislocation microstructures in the materials. Yin et al. [30] investigated the effect of orientation on tensile-compressive properties through 3D MD simulation and found that the [001] orientation performed more evidently. In these researches, MD simulation fully demonstrates its advantages in capturing the microstructure evolution in real time. However, there is still a lack of systematic study of Ni-based superalloy deformation under complex stress states.
Based on the above discussion, the simulation results of tension and compression under uniaxial, equibiaxial, and non-equibiaxial loadings are compared and analyzed systematically, which can provide guidance for the design of working loading states of SC Ni-based superalloy. In addition, the strengthening mechanism of SC Ni-based superalloy under different loading conditions has been concluded in this paper. The MD model and loading conditions are constructed in detail in Section 2. Then, Section 3 presents the strainstress relation and the dislocation evolution during the loading process for all loading modes. In Section 4, the dissociation of the initial dislocation network is explained, and the dislocation bowing mechanism is illustrated by analyzing the energy density distribution and calculation of the Schmid factor. Finally, this paper is summarized in Section 5.

Model and Method
To investigate the microstructure evolution of SC Ni-Al superalloy under multiaxial loading, large-scale MD simulations are performed using the classical MD code LAMMPS [31]. For the MD method, the input parameters are the atomic potential, atomic structure, ensemble, boundary conditions, loading states, temperature, etc. The atomic interactions are governed by the embedded atom method (EAM) potential developed by Mishin [32], which describes lattice properties of γ and γ fields on the Ni-Al phase diagram and compares fairly well with first-principles calculations. From the EAM potential, the lattice parameters of γ and γ phases can be obtained as α γ = 0.352 nm and α γ = 0.357 nm, respectively. The lattice mismatch δ is then calculated as: As the bulk lattice parameters of γ and γ phases nearly commensurate with the ratio n = α γ : α γ ≈ 70 : 69, the lattice number ratio at γ /γ interface is set as 69:70 to reduce the misfit stress and form the misfit dislocations network. The 3D MD model of SC Ni-Al superalloy is built by embedding a 69 × 69 × 69 cubical γ precipitate phase in a 100 × 100 × 100 γ matrix phase [19], as shown in Figure 1. The peripheral red atoms represent the Ni atoms in γ matrix phase, whereas the blue and white atoms in the center represent the Ni atoms and Al atoms in γ precipitate, respectively. A relatively wide matrix phase width is defined in this model to reduce the influence of periodic boundary conditions [33]. Accordingly, the initial box sizes of our typical MD models are l x = l y = l z = 31.7 nm, where the x, y, z axis directions represent [100], [010], and [001] orientations, respectively. Moreover, throughout this paper, the γ channels normal to the x, y, z axes denote X channel, Y channel, and Z channel, respectively. The length of the internal cubic γ precipitate phase is 24.46 nm, which indicates that the volume fraction of γ phase is 46% less than the fifth-generation SC Ni-based superalloy [7]. As this work focuses on the microstructure evolution during the loading process, this difference will not have a considerable impact on our study. The entire model is comprised of 2,858,036 atoms, of which 2,529,527 are Ni atoms and 328,509 are Al atoms.
To investigate the microstructure evolution of SC Ni-Al superalloy under multiaxial loading, large-scale MD simulations are performed using the classical MD code LAMMPS [31]. For the MD method, the input parameters are the atomic potential, atomic structure, ensemble, boundary conditions, loading states, temperature, etc. The atomic interactions are governed by the embedded atom method (EAM) potential developed by Mishin [32], which describes lattice properties of γ and γ′ fields on the Ni-Al phase diagram and compares fairly well with first-principles calculations. From the EAM potential, the lattice parameters of γ and γ′ phases can be obtained as = 0.352 nm and ′ = 0.357 nm, respectively. The lattice mismatch is then calculated as: As the bulk lattice parameters of γ and γ′ phases nearly commensurate with the ratio = ′ : ≈ 70: 69, the lattice number ratio at γ′/γ interface is set as 69:70 to reduce the misfit stress and form the misfit dislocations network. The 3D MD model of SC Ni-Al superalloy is built by embedding a 69 × 69 × 69 cubical γ′ precipitate phase in a 100 × 100 × 100 γ matrix phase [19], as shown in Figure 1. The peripheral red atoms represent the Ni atoms in γ matrix phase, whereas the blue and white atoms in the center represent the Ni atoms and Al atoms in γ′ precipitate, respectively. A relatively wide matrix phase width is defined in this model to reduce the influence of periodic boundary conditions [33]. Accordingly, the initial box sizes of our typical MD models are = = = 31.7 nm, where the X, Y, Z axis directions represent [100], [010], and [001] orientations, respectively. Moreover, throughout this paper, the γ channels normal to the X, Y, Z axes denote X channel, Y channel, and Z channel, respectively. The length of the internal cubic γ′ precipitate phase is 24.46 nm, which indicates that the volume fraction of γ′ phase is 46% less than the fifth-generation SC Ni-based superalloy [7]. As this work focuses on the microstructure evolution during the loading process, this difference will not have a considerable impact on our study. The entire model is comprised of 2,858,036 atoms, of which 2,529,527 are Ni atoms and 328,509 are Al atoms. The MD simulations are performed with an NPT ensemble. Periodic boundary condition is applied on the three axes and the simulations are conducted with an integration time step of 1.0 × 10 -15 s. First, the model is heated up to a temperature of 300 K by velocity rescaling and a Nose-Hoover thermostat is applied to maintain the temperature for 10,000 timesteps (10 picoseconds, ps) as a relaxing process. Meanwhile, the MD The MD simulations are performed with an NPT ensemble. Periodic boundary condition is applied on the three axes and the simulations are conducted with an integration time step of 1.0 × 10 −15 s. First, the model is heated up to a temperature of 300 K by velocity rescaling and a Nose-Hoover thermostat is applied to maintain the temperature for 10,000 timesteps (10 picoseconds, ps) as a relaxing process. Meanwhile, the MD box dimensions are adjusted by a Berendsen barostat to maintain traction-free periodic boundary conditions in all three directions. The loadings are applied by deforming the simulation box at a fixed strain rate of ε 3) for tension and compression, respectively. It is worth mentioning that the loading is negative for compression, while maintaining zero  Figure 2. The loading lasts 500 ps until ε zz = 5% for each simulation. The ε and ε • represent the engineering strain and engineering strain rate in this paper, which are calculated by the ratio of length at the moment to the original length and its derivative with time, respectively. The loading states in this paper are considerably similar to the authors' previous work [34]. The dislocation evolutions are analyzed by the dislocation analysis (DXA) proposed by Stukowski et al. [35], while the blue lines are perfect dislocations, green lines are Shockley partials, red lines are stair-rod partials, and yellow lines are stacking faults (SFs) atoms. The atomic configurations are visualized using visualization software OVITO in the subsequent analysis [36]. The output strain and stress for strain-stress curves are calculated by the box size evolution and global pressure from the built-in algorithm in LAMMPS. box dimensions are adjusted by a Berendsen barostat to maintain traction-free periodic boundary conditions in all three directions. The loadings are applied by deforming the simulation box at a fixed strain rate of • = 10 8 s -1 , • = • = 10 8 s -1 ( • : • = 1: 1), and • = 1 3 ×10 8 s -1 , • = 10 8 s -1 ( • : • = 1:3) for tension and compression, respectively. It is worth mentioning that the loading is negative for compression, while maintaining zero box pressures in the other directions, as shown in Figure 2. The loading lasts 500 ps until = 5% for each simulation. The and • represent the engineering strain and engineering strain rate in this paper, which are calculated by the ratio of length at the moment to the original length and its derivative with time, respectively. The loading states in this paper are considerably similar to the authors' previous work [34]. The dislocation evolutions are analyzed by the dislocation analysis (DXA) proposed by Stukowski et al. [35], while the blue lines are perfect dislocations, green lines are Shockley partials, red lines are stair-rod partials, and yellow lines are stacking faults (SFs) atoms. The atomic configurations are visualized using visualization software OVITO in the subsequent analysis [36]. The output strain and stress for strain-stress curves are calculated by the box size evolution and global pressure from the built-in algorithm in LAMMPS.  Figure 3a displays the strain-stress diagrams for uniaxial loading. Here, the absolute value of strain-stress is taken for compressive loading to show the difference between tension and compression. The marked dot in Figure 3a represents the moments of the initial stage of dislocation emission, the stage at which dislocations are about to cut through γ channel, and the stage when γ channel is already penetrated by dislocations, respectively. Figures 3b1-b3 and 3c1-c3 show the snapshots of these marked moments. These snapshots are analyzed by DXA and the transparent green cube at the center represents γ′ phase, while the stacking faults (SFs) atoms are coloured in yellow and the FCC atoms are hidden.  Figure 3a displays the strain-stress diagrams for uniaxial loading. Here, the absolute value of strain-stress is taken for compressive loading to show the difference between tension and compression. The marked dot in Figure 3a represents the moments of the initial stage of dislocation emission, the stage at which dislocations are about to cut through γ channel, and the stage when γ channel is already penetrated by dislocations, respectively. Figure 3(b1-b3,c1-c3) show the snapshots of these marked moments. These snapshots are analyzed by DXA and the transparent green cube at the center represents γ phase, while the stacking faults (SFs) atoms are coloured in yellow and the FCC atoms are hidden.

Uniaxial Loading
In Figure 3a, the uniaxial strain-stress curve for compression has a higher gradient at the early stage. By linear fitting with the elastic segment of the strain-stress diagram, the slopes represent the elastic modulus which are 96.5 and 149.3 GPa for tension and compression, respectively. This indicates that the SC Ni-based superalloy presents a larger elastic modulus for [100] orientation compression. Meanwhile, the maximum stress for compression is greater than the tensile stress and the tensile curve turns at a later time. In Figure 3a, the uniaxial strain-stress curve for compression has a higher gradien at the early stage. By linear fitting with the elastic segment of the strain-stress diagram the slopes represent the elastic modulus which are 96.5 and 149.3 GPa for tension and compression, respectively. This indicates that the SC Ni-based superalloy presents larger elastic modulus for [100] orientation compression. Meanwhile, the maximum stress for compression is greater than the tensile stress and the tensile curve turns at later time.
In Figure 3b1-b3, it can be seen that the decomposed Shockley dislocations slip int [001] γ matrix channel, which follows the same direction of uniaxial loading at = 2.2%. The primary dislocation reaction has already been discussed in the authors' previ ous work [34]: The original 1 2 〈110〉 pure edge misfit dislocation dissociates into one stair-rod dislo cation locked in γ′/γ interface &#x¯; and two different Shockley partials with stackin faults bowing out into [010] γ matrix channel, which is perpendicular to the externa loading direction [37]. Then, the motion of decomposed dislocations is resisted by the γ precipitate phase in X and Y directions. As a result, the bowing-out dislocations ar mainly located in the channel of Z direction at = 3.5%. When reaches 3.6%, the dis locations emitted from upper and lower interfaces start to intersect in γ phase, while th original pure edge dislocations at the interface turn into mixed dislocations. At this mo In Figure 3(b1-b3), it can be seen that the decomposed Shockley dislocations slip into [001] γ matrix channel, which follows the same direction of uniaxial loading at ε Z = 2.2%. The primary dislocation reaction has already been discussed in the authors' previous work [34]: The original 1 2 110+ pure edge misfit dislocation dissociates into one stair-rod dislocation locked in γ /γ interface and two different Shockley partials with stacking faults bowing out into [010] γ matrix channel, which is perpendicular to the external loading direction [37]. Then, the motion of decomposed dislocations is resisted by the γ precipitate phase in x and y directions. As a result, the bowing-out dislocations are mainly located in the channel of z direction at ε Z = 3.5%. When ε Z reaches 3.6%, the dislocations emitted from upper and lower interfaces start to intersect in γ phase, while the original pure edge dislocations at the interface turn into mixed dislocations. At this moment, the dislocations are still largely in the channel of Z direction, but a few dislocations start to bow out at x and y directions since the dislocation networks begin to dissociate.
In Figure 3(c1-c3), the snapshots of uniaxial compression are displayed and the motion of dislocation is very similar to the equibiaxial tension, as shown in Figure 4(b1-b3). Evidently, the main dislocation reaction during compression in Equation (2) and the dislocations cut through the γ phase near the stress peak of strain-stress curves. However, the dislocation evolution is considerably different from the tensile loading. At ε Z = −2.0%, the dislocations primarily dissociate at Z channel and then bow out toward the x and y directions. Therefore, when ε Z reaches −2.1%, the dislocations are mainly distributed in X and Y channels, similar to the equibiaxial loading on x and y directions. In Figure 3(b1-b3,c1-c3), the dislocation network has been severely damaged, due to the fact that the strength of material rapidly decreased, as shown in Figure 3a.
tions start to bow out at X and Y directions since the dislocation networks begin to dissociate.
In Figure 3c1-c3, the snapshots of uniaxial compression are displayed and the motion of dislocation is very similar to the equibiaxial tension, as shown in Figure 4b1-b3. Evidently, the main dislocation reaction during compression in Equation (2) and the dislocations cut through the γ phase near the stress peak of strain-stress curves. However, the dislocation evolution is considerably different from the tensile loading. At = −2.0%, the dislocations primarily dissociate at Z channel and then bow out toward the X and Y directions. Therefore, when reaches −2.1%, the dislocations are mainly distributed in X and Y channels, similar to the equibiaxial loading on X and Y directions. In Figure 3b,c, the dislocation network has been severely damaged, due to the fact that the strength of material rapidly decreased, as shown in Figure 3a.

Equibiaxial Loading
In Figure 4, the diagrams show the strain-stress lines and microstructure evolution under equibiaxial loading, while the snapshots are conducted in Figure 3. Figure 4a illustrates the strain-stress relation under equibiaxial loading. Here, the tensilecompressive curves have the same gradient at the early stage. As the strain reaches = 0.5%, the compressive graphs show a higher elastic modulus, which is 202.4 GPa compared with 157.4 GPa for tensile stress and eventually turns at almost the same strain. The X and Z direction curves are in practical unanimity until the turning point for both loading states.
In Figure 4b1-b3, the interfacial dislocations in the Y direction, which bear no loading, are first decomposed by the reaction above at = 1.0%. As increases to 1.7%, the

Equibiaxial Loading
In Figure 4, the diagrams show the strain-stress lines and microstructure evolution under equibiaxial loading, while the snapshots are conducted in Figure 3. Figure 4a illustrates the strain-stress relation under equibiaxial loading. Here, the tensile-compressive curves have the same gradient at the early stage. As the strain reaches ε zz = 0.5%, the compressive graphs show a higher elastic modulus, which is 202.4 GPa compared with 157.4 GPa for tensile stress and eventually turns at almost the same strain. The x and z direction curves are in practical unanimity until the turning point for both loading states.
In Figure 4(b1-b3), the interfacial dislocations in the y direction, which bear no loading, are first decomposed by the reaction above at ε Z = 1.0%. As ε Z increases to 1.7%, the decomposed dislocations cut into the channel at x and z directions, which have the loading, while the dislocations that cut toward the y direction are obstructed by the γ phase at the same time. When the dislocations cut through the X and Z channels at ε Z = 1.8%, there are still no dislocations in the Y channel. As shown in Figure 4(c1-c3), the dislocation evolution under equibiaxial compressive loading of x and z directions is the same as the uniaxial tension at y direction. The difference is that the dislocation networks decompose rapidly when the dislocations cut through the Y channel. Similar to uniaxial loading, the moment when the dislocations cut through the channel is also around the turning point in Figure 4a, which confirms that the completeness of the dislocation networks significantly contributes to the material strength. Figure 5 shows the strain-stress relation and microstructure change under nonequibiaxial loading, and the pictures are modified as previously mentioned. In particular, for comparison, the strain in the x direction for non-equibiaxial loading is magnified by three times. Figure 5a shows that the non-equibiaxial strain-stress curves, such as equibiaxial loading, also bend at a similar strain. In contrast to equibiaxial loading, the strain-stress curves of z direction are almost the same for tension and compression, while the elastic modulus for tension and compression in z direction are 144.9 and 149.2 GPa, respectively. Moreover, the curves of x direction show a higher compressive stress compared with tensile stress and the elastic modulus are 231.6 and 335.4 GPa for tension and compression in x direction. Compared with Figure 5(b2,c2), the dislocations in X channel are more stable than in Z channel. This is the reason the elastic modulus in x direction is larger than in z direction. From these strain-stress diagrams, it can be seen that the maximum stress for compression is slightly larger than tension at the same strain for all of the loading conditions. phase at the same time. When the dislocations cut through the X and Z channels at = 1.8%, there are still no dislocations in the Y channel. As shown in Figure 4c1-c3, the dislocation evolution under equibiaxial compressive loading of X and Z directions is the same as the uniaxial tension at Y direction. The difference is that the dislocation networks decompose rapidly when the dislocations cut through the Y channel. Similar to uniaxial loading, the moment when the dislocations cut through the channel is also around the turning point in Figure 4a, which confirms that the completeness of the dislocation networks significantly contributes to the material strength. Figure 5 shows the strain-stress relation and microstructure change under nonequibiaxial loading, and the pictures are modified as previously mentioned. In particular, for comparison, the strain in the X direction for non-equibiaxial loading is magnified by three times. Figure 5a shows that the non-equibiaxial strain-stress curves, such as equibiaxial loading, also bend at a similar strain. In contrast to equibiaxial loading, the strain-stress curves of Z direction are almost the same for tension and compression, while the elastic modulus for tension and compression in Z direction are 144.9 and 149.2 GPa, respectively. Moreover, the curves of X direction show a higher compressive stress compared with tensile stress and the elastic modulus are 231.6 and 335.4 GPa for tension and compression in X direction. Compared with Figure 5b2,c2, the dislocations in X channel are more stable than in Z channel. This is the reason the elastic modulus in X direction is larger than in Z direction. From these strain-stress diagrams, it can be seen that the maximum stress for compression is slightly larger than tension at the same strain for all of the loading conditions. For non-equibiaxial loading in Figure 5(b1-b3), the dislocation evolution at the initial stage is the same as equibiaxial loading. As the strain increases, the dislocations also appear at the X and Z channels, but the dislocations move more slowly in the x direction with lower loading evidently. In Figure 5(c1-c3), the dislocations bow out into the X and Y channels, which sustain a smaller loading compared with Z channel and no loading, respectively. The velocity of dislocation bowing is higher at y direction in comparison with x direction. In contrast to the uniaxial and equibiaxial loadings, the models still maintain the elastic performance after the first dislocation cuts through γ channel, as shown in Figure 5a. In Figure 5(b3,c3), the dislocation networks still have a relative complete morphology compared with Figure 4(b3,c3). As the dislocations in the x direction move more slowly than in other channels, when dislocations in other channels have penetrated the channel, they still move in the γ channel, which leads the dislocation network to remain more intact.

Non-Equibiaxial Loading
In summary, the dissociated dislocations always move into the tensile loading direction channels, and the bowing velocity shows a positive correlation with the magnitude of loading, while the compressive loading shows the opposite results. From the comparison of dislocation motion under tensile-compressive loadings, it can be concluded that tensile loading promotes the bowing out of dislocation into the loading channel, while compressive loading suppresses the bowing out of dislocation into the loading direction. In addition, the correlation of the turning point on the strain-stress curves to the moment in which the decomposed dislocations cut through the γ phase indicates that the accumulation of dislocation in γ channel is the starting point of the plastic deformation of SC Ni-based superalloy [38].
Moreover, compared with these strain-stress lines, an evident tension-compression asymmetry for SC Ni-based superalloy is observed during the loading process. This is due to the fact that the alloy crystals are subjected to different stresses as a result of different intergranular slip directions during tension and compression [39]. From our simulation, it can be concluded that the compression requires greater stress in the 100 direction. From the experimental results by Wen et al. [25] and Leidermark et al. [40], the tension-compression asymmetry for [001] and [111] orientations is small, while the strength of [011] orientation is stronger in compression than in tension. Leidermark et al. found that the tensioncompression asymmetry significantly relates to the critical resolved shear stress (CRSS) and the CRSS for compression in [001] direction is only slightly larger than tension. However, it has to be pointed out that the atomistic simulation is distinguished from the macroscopic experiment, as the real materials have many defects and the MD models are generally perfect in crystal structure. As a result, the value of CRSS can be magnified many times in atomistic simulation leading to an exaggerate tension-compression asymmetry in [001] direction in this paper. On the other hand, by comparing these graphs, it can be observed that the maximum stress of tension only changes from 10~20% under different loading states, while the maximum stress of compression at z direction transforms around 45% between uniaxial and equibiaxial loadings. This result demonstrates that the maximum stress for compression is more sensitive with biaxial loading.

Discussion
To explain the microstructure evolution, we first investigate the initial dislocation structure. Figure 6 illustrates a detailed view of the initial dislocations, where the dislocation line has been emphasized by an arrow, red atoms represent Ni atoms in γ phase, blue atoms represent Ni atoms in γ phase, and white atoms represent Al atoms in γ phase. The slip directions of [100] and [111] orientations have been indicated in dashed green line. From the DXA algorithm, the dislocation is considered as a perfect pure edge dislocation, pointed out by a black arrow in Figure 6. Nevertheless, the perfect dislocation is surrounded by both the γ and γ phases and the interface at (010) plane, which indicates that the movements of the dislocation into different phases belong to different slip systems at the general (111) slipping plane. In fact, the initial dislocation is known as Lomer dislocation, which is formed by dislocation reaction given below: Moreover, it acts as a barrier to the glide of further dislocations on the two {111} planes, which is known as Lomer-Cottrell lock, as shown in Figure 7 [37]. When the Shockley partials completely cut through the γ channel, the stress of alloy rapidly decreases. This is because the interaction between the dislocations breaks the Lomer-Cottrell lock and leads to a degeneration of material strength. The strengthening mechanism is considerably different from the normal metal materials, in which the yield stress is related to the dislocation nucleation [38,42].  To investigate the mechanism of the Lomer dislocation dissociation in the interface, the energy density distribution of initial dislocation networks along all three directions Under all of the loading states, the Lomer dislocation scarcely glides on any of the four {111} planes and is always sessile as a result of its location in the (100) plane [41]. Therefore, the reaction mentioned above represents the classical Lomer dislocation dissociation with a reduction in b 2 , according to: where b is the module of Burgers vector and a is the lattice parameter. Therefore, the stair-rod partial exerts a repulsive force on the two Shockley partials and these three partial dislocations form a stable, sessile arrangement under incipient loading conditions. Moreover, it acts as a barrier to the glide of further dislocations on the two {111} planes, which is known as Lomer-Cottrell lock, as shown in Figure 7 [37]. When the Shockley partials completely cut through the γ channel, the stress of alloy rapidly decreases. This is because the interaction between the dislocations breaks the Lomer-Cottrell lock and leads to a degeneration of material strength. The strengthening mechanism is considerably different from the normal metal materials, in which the yield stress is related to the dislocation nucleation [38,42].
as Lomer dislocation, which is formed by dislocation reaction given below: Under all of the loading states, the Lomer dislocation scarcely glides on any of the four {111} planes and is always sessile as a result of its location in the (100) plane [41]. Therefore, the reaction mentioned above represents the classical Lomer dislocation dissociation with a reduction in 2 , according to: where is the module of Burgers vector and is the lattice parameter. Therefore, the stair-rod partial exerts a repulsive force on the two Shockley partials and these three partial dislocations form a stable, sessile arrangement under incipient loading conditions. Moreover, it acts as a barrier to the glide of further dislocations on the two {111} planes, which is known as Lomer-Cottrell lock, as shown in Figure 7 [37]. When the Shockley partials completely cut through the γ channel, the stress of alloy rapidly decreases. This is because the interaction between the dislocations breaks the Lomer-Cottrell lock and leads to a degeneration of material strength. The strengthening mechanism is considerably different from the normal metal materials, in which the yield stress is related to the dislocation nucleation [38,42].  To investigate the mechanism of the Lomer dislocation dissociation in the interface, the energy density distribution of initial dislocation networks along all three directions To investigate the mechanism of the Lomer dislocation dissociation in the interface, the energy density distribution of initial dislocation networks along all three directions is shown in Figure 8. The energy density is calculated by summing up the potential energy of each atom and then dividing by the volume. In addition, only the central part of the atom is selected in the calculation to avoid the influence of the interface at the other sides. In Figure 8, two dashed lines mark the location of the interface. From the curves, it is evident that the energy density of γ phase is considerably larger than γ phase and the atoms in γ phase have a higher energy density as they get closer to the interface. For all of the three directions, the energy density distribution curves are almost identical. As the dislocation tends to slip toward the lower energy states, the difference in energy density between γ and γ phases acts as an energy barrier, which allows the Lomer dislocations to dissociate toward γ in most of our simulation cases. the atom is selected in the calculation to avoid the influence of the interface at the other sides. In Figure 8, two dashed lines mark the location of the interface. From the curves, it is evident that the energy density of γ phase is considerably larger than γ′ phase and the atoms in γ phase have a higher energy density as they get closer to the interface. For all of the three directions, the energy density distribution curves are almost identical. As the dislocation tends to slip toward the lower energy states, the difference in energy density between γ and γ′ phases acts as an energy barrier, which allows the Lomer dislocations to dissociate toward γ′ in most of our simulation cases. Traditionally, the Schmid factor of decomposed Shockley partials, which indicates the magnitude of the resolved shear stress, is calculated as Equation (5): SF = cos cos (5) where is the angle between applied force and the norm of the slip plane, is the angle between applied force and the slip direction. The calculation details in each channel under different tension loading states have been shown in Table 1. For the Lomer dislocations in each channel, they all dissociated into two types of different Shockley partials, as shown in Table 1. The Schmid factors have been listed corresponding to the Burgers vector of slipping Shockley partials, which is calculated by combining the Burgers vectors and resultant force direction vector with Equation (5). The resultant force direction vector is calculated by summation of forces along X, Y, and Z directions. It is worth mentioning that the resultant force direction vector for non-equibiaxial loading, which is not on account of the loadings, is conducted by the strain rate. The actual stress ratio between the X and Z directions is 1:2 rather than 1:3. For Z direction uniaxial tension, the Schmid factor of dislocations in Z direction is 0.47, which is greater than 0.24 in X and Y directions. As a result, the dissociation of Lomer dislocations first bows out in Z direction and then the dislocations in X and Y channels start to decompose. Nevertheless, for the dislocations in non-loading channels, they prefer to bow toward γ′ phase due to the stress caused by necking. Under biaxial loading, the mechanism is more complex. Under Traditionally, the Schmid factor of decomposed Shockley partials, which indicates the magnitude of the resolved shear stress, is calculated as Equation (5): where φ is the angle between applied force and the norm of the slip plane, λ is the angle between applied force and the slip direction. The calculation details in each channel under different tension loading states have been shown in Table 1. For the Lomer dislocations in each channel, they all dissociated into two types of different Shockley partials, as shown in Table 1. The Schmid factors have been listed corresponding to the Burgers vector of slipping Shockley partials, which is calculated by combining the Burgers vectors and resultant force direction vector with Equation (5). The resultant force direction vector is calculated by summation of forces along x, y, and z directions. It is worth mentioning that the resultant force direction vector for non-equibiaxial loading, which is not on account of the loadings, is conducted by the strain rate. The actual stress ratio between the x and z directions is 1:2 rather than 1:3. For z direction uniaxial tension, the Schmid factor of dislocations in z direction is 0.47, which is greater than 0.24 in x and y directions. As a result, the dissociation of Lomer dislocations first bows out in z direction and then the dislocations in X and Y channels start to decompose. Nevertheless, for the dislocations in non-loading channels, they prefer to bow toward γ phase due to the stress caused by necking. Under biaxial loading, the mechanism is more complex. Under equibiaxial loading, the Schmid factor is 0.24 and 0 for the two decomposed Shockley partials in X and Z channels, while it is 0.47 at Y channel. In x and z directions, the Lomer dislocations cannot dissociate by virtue of only one Shockley partial, which has resolved shear stress, as shown in Figure 4(b1-b3). On the other hand, although the dislocations in Y channel are blocked by γ phase, the Shockley partials which are close to γ phase can still cut into γ phase in the X and Z channels under the continuous loading, resulting in numerous dislocations in X and Z channels, as shown in Figure 9. For non-equibiaxial loading, the Schmid factors of dislocation in X channel are 0.189 and 0, thus the interfacial dislocations in this direction do not decompose. In addition, since the Schmid factors of the dislocations in z direction are 0.42 and 0.24, while the ones in the y direction are both 0.42, the dislocations in the y direction interface are decomposed and cut through the γ channel first due to the fact that the dislocation decomposition is limited by the minimum Schmid factor. ited by the minimum Schmid factor.  For compressive loading, the energy barrier at the interface plays an important role in the bowing direction of dislocations. Under uniaxial compressive loading, the Lomer dislocation in z direction interface decomposes directly and is caused by the resolved shear stress. Meanwhile, the dislocations in x and y direction interface are affected by the tension stress as a result of the Poisson effect, and the dislocations tend to decompose directly into the γ channel phase. However, due to the action of the energy barrier between the two phases, the dislocations cannot be decomposed toward the γ phase, resulting in a situation similar to equibiaxial tension. The dislocations under equibiaxial and nonequibiaxial loading work in the same way. While under equibiaxial compressive stress states, the Lomer dislocation dissociates and bows out toward the γ channels different from compressive direction, while the dislocations maintain the original state in the loading direction interface. On the whole, the mechanism of dislocation bowing is controlled by the resolved shear stress together with the energy barrier at the interface.

Conclusions
In this paper, MD simulation is applied to investigate the microstructure evolution of SC Ni-based superalloy under multiaxial tension and compression. On the basis of our simulation results and dislocation analyses, several significative conclusions can be drawn: (1) The strain-stress response demonstrates a significant tension-compression asymmetry, and the strength of compression is larger than that under tension, which is consistent with the experimental observation. In addition, contrary to the tensile loading states, in which the maximum stress does not significantly change between uniaxial and multiaxial loadings, the maximum stress under compressive loading varies evidently. (2) From the atomistic simulation results, we conclude that tensile loading can promote the bowing of dislocations into the channel, while compressive loading restrains it. This phenomenon is believed to be caused by a coupling effect of both the difference of resolved shear stress and energy barrier at the interface. (3) The initial dislocation networks are comprised of Lomer dislocations, which dissociate and form the Lomer-Cottrell lock, acting as a barrier to the further glide of dislocations. Meanwhile, as the Lomer-Cottrell lock disappears, the tensile stress of the alloy declines rapidly, which indicates that the main strengthening mechanism of the SC Ni-based superalloy is the formation of the Lomer-Cottrell lock. To initiate the decomposition of one Lomer dislocation, resolved shear stress applied in both dissociated Shockley partials are required to exceed the critical value, that is why a large amount of Lomer dislocations remain undissociated under loading. (4) During the loading process, the dissociated direction of Lomer dislocation is affected by the energy barrier between γ and γ phase. Nevertheless, in the subsequent dislocation movement, the resolved shear stress plays a dominant role and the slip of Shockley partials is consistent with the calculation results of Schmid factors.