Effect of Vacancies on Dynamic Response and Spallation in Single-Crystal Magnesium by Molecular Dynamic Simulation

: The effect of vacancies on dynamic response and spallation in single-crystal magnesium (Mg) is investigated by nonequilibrium molecular dynamics simulations. The initial vacancy concentration ( C v ) ranges from 0% to 2.0%, and the shock loading is applied along [0001] and [10–10] directions. The simulation results show that the effects of vacancy defects are strongly dependent on the shock directions. For shock along the [0001] direction, vacancy defects have a negligible effect on compression-induced plasticity, but play a role in increasing spall damage. In contrast, for shock along the [10–10] orientation, vacancy defects not only provide the nucleation sites for compression-induced plasticity, which mainly involves crystallographic reorientation, phase transition, and stacking faults, but also signiﬁcantly reduce spall damage. The degree of spall damage is probably determined by a competitive mechanism between energy absorption and stress attenuation induced by plastic deformation. Void evolution during spallation is mainly based on the emission mechanism of dislocations. The {11–22} <11–23> pyramidal dislocation facilitates the nucleation of void in the [0001] shock, as well as the {1–100} <11–20> prismatic dislocation in the [10–10] shock. We also investigated the variation of spall strength between perfect and defective Mg at different shock velocities. The relevant results can provide a reference for future investigations on spall damage.


Introduction
Spallation is an important subject in dynamic damage research, which can be observed in the local tensile stress region caused by the interaction of the unloading compression wave and rarefaction wave formed from the free surface of the target. Experimentally, the most widely used methods to investigate spall failure are plate impact experiments [1][2][3] or laser pulse [4][5][6]. For example, Wang et al. [3] compared the spall strength of four aluminum materials and discussed the damage mechanism by plate impact experiments. By using a laser shock, the spall strength and fragmentation of single-crystal magnesium at very high strain rates of about 10 7 s −1 were provided [6]. In experiments, the study of spall failure was to obtain effective information about spall strength (σ sp ) and dynamics by measuring and analyzing the velocity history of free surface of the samples [7]. However, it was difficult to describe and identify the individual contributions of various microscopic processes to the overall material reaction because of the very short time and length scale of the processes involved. Molecular dynamics (MD) simulations make it possible to study atomic level evolution of microstructure, and have been widely used to investigate the dynamic failure behavior in different materials [8][9][10][11][12][13][14]. For example, Luo et al. [10] studied spall damage of single-crystal Cu under square and Taylor shock wave loading using MD simulations. The results showed that the spall strength of Taylor wave loading was higher than that of square shock loading at the same shock velocities. Agarwal and Dongare [14] investigated the effects of loading direction, impact pressure, and initial temperature on spallation of single-crystal magnesium with MD simulations, and found that the spall strength decreased with the increase in impact strength and initial temperature. Gunkelmann et al. [11] performed MD simulations of shock-induced spallation in polycrystalline iron and found that spall damage was strongly affected by the phase transition, leading to the change of fracture surface morphology. Chen et al. [13] studied the effect of grain boundaries (GBs) on spall behavior in the metal Ta by MD simulations. They demonstrated the correlation between GB structure, GB plasticity mechanism, and the resulting spall strength. Wang et al. [8] found that release melting and shock melting occurred when the shock velocities increased in the metal tin via MD. The spall strength decreased rapidly when release melting occurred, while the downtrend of the spall strength became gentle during shock melting. In addition to the effects of pressure, temperature, and shock velocity, the influence of defect structures on the spallation of materials has also been given attention in the abovementioned studies. However, considering the wide variety of defects in actual metals and alloys, individual contributions are difficult to separate. Vacancy defects are not only ubiquitous in materials, but also one of the simplest forms of defects. Under environmental conditions, the vacancy concentration of the material is very low, such that its effects are considered to be negligible [15]. However, materials with high vacancy concentration can be obtained by ion radiation [16], rapid quenching [17], and thin film growth [18]. Luo et al. [19] studied spallation of Cu with C v up to 2.0%, and the results showed that the shear flow strength and spall strength decreased with C v increasing. Lin et al. [20] further studied and found that the effect of vacancy defects on the spallation damage of Cu was anisotropic. Spallation in Ni with C v ranging from 0% to 2.0% was explored by Qiu et al. [21]. They showed that the vacancy concentration had two competing effects on the spallation, one was to dissipate some of the applied stress, and the other was to provide sites for void nucleation and accelerate dislocation motion. However, most MD simulations of the effect of high vacancy concentration on spallation are carried out on FCC metals, there is a clear lack of reports on HCP metals. As typical HCP materials, magnesium and magnesium alloys have been widely used in the aerospace, automobile, defense, and medical fields due to their excellent properties of low density and high strength [22][23][24]. In practical engineering materials, the effects of unavoidable or artificial defects, such as voids [25][26][27][28], cracks [29], grain boundaries [30,31], and dislocations [32], on the plasticity of Mg materials, have received extensive attention. For example, Li et al. [28] revealed two typical void collapse mechanisms of nanoporous magnesium under shock loading by MD simulations: the plasticity mechanism and the internal jetting mechanism. Tang et al. [29] investigated the fatigue crack propagation behavior of single-crystal Mg by using MD simulations and found that the fatigue crack growth rate and the crack path vary significantly with the crystallographic orientations of the initial crack. However, there has been little attention paid to the influence of the defect structure in Mg single crystals on spallation, or the microscopic mechanism of spallation evolution. Therefore, it is of great significance to study the spallation mechanism of magnesium with vacancy defects under dynamic loading.
In this work, the main purpose is to analyze the effect of vacancy concentration on spallation of single-crystal magnesium under different loading directions, with an initial concentration ranging from 0% to 2.0%. Large scale MD simulations of shock loading of single-crystal Mg were first performed to investigate the anisotropic deformation and spallation response at a piston velocity of 0.5 km/s along the [0001] and  loading directions. The void evolution of initial spallation was also analyzed and discussed in detail. In addition, to investigate the effect of shock strength, MD simulations were carried out at piston velocities ranging from 0.5 km/s to 2.0 km/s. Meanwhile, the coupling between anisotropic deformation and initial spallation was also analyzed and discussed in detail under weak shock loading. The outline of the paper is organized as follows. The

Models and Simulation Methods
Large-scale atomic/molecular massively parallel simulator (LAMMPS, version 2013, copyright Sandia Corporation) codes [33] were used for the molecular dynamics simulations with a Mg embedded atom method (EAM) potential in the form of a Finnis-Sinclair (FS) function [34]. Compared with the existing empirical or semi-empirical interatomic potentials, this potential effectively describes the plasticity and phase transition properties of Mg at high pressure. Vacancy-free single-crystal Mg samples were first prepared with the Cartesian (x, y, z) axes corresponding to the (   ) crystal orientations, and they had sizes of 20.0 nm × 19.9 nm × 200.1 nm and 19.9 nm × 19.8 nm × 200.1 nm, respectively. The vacancies were achieved by randomly removing atoms from the perfect configuration, and the vacancy concentration in our simulations ranged from 0% to 2.0%. The samples were first equilibrated at 300 K and zero pressure with the constant pressure-temperature (NPT) ensemble in three-dimensional (3D) periodic boundary conditions for 100 ps. After being quenched, vacancies existed predominantly as monovacancies, followed by divacancies, and there were only a few vacancy clusters in all samples. During dynamic loading, the shock wave propagated along the z direction with a microcanonical ensemble. Periodical boundary conditions were applied only in the x and y directions, while the surface along the z direction was regarded as the free surface. The shock simulation was performed by using a rigid piston composed of several layers of atoms on one end of samples along z direction [35]. The moving infinite massive piston impacted the rest of the atoms at a constant loading velocity, V p , to generate the shock pulse. The time step for integrating the equation of motion was set to 1 fs and the shock loading period was up to 25 ps. Then, the atomic piston was removed, generating a release wave from the side of piston. This work mainly discusses the dynamic response and spallation under the piston velocity of 0.5 km/s.
The microstructures (such as crystal structure types and local crystal defects) and voids evolution during the simulations were analyzed by adaptive common neighbor analysis (ACNA) [36] and "construct surface mesh" within OVITO [37]. Dislocation extraction algorithm (DXA) [38] was used to identify various dislocation types and the directions of dislocation slip. A binning analysis [39] was used to obtain the wave profiles by evenly dividing the simulated system into many slices along the z direction, and averaging the thermodynamics quantities, such as local stress, σ ij (i, j = x, y, z), temperature, T, and particle velocity in each slice. Local stress is an important reference quantity, and the expression is defined as follows [35]: where Λ α → r is 1 if α lies within the average volume Ω centered on; otherwise, it is 0. Function is the fraction length of bond that lies within the average volume and satisfies The free surface velocity V fs is measured by averaging the total atomic velocities of the two atomic layers at the end of the sample.

Shock Spallation
For a more intuitive understanding of the shock compression, release, tension, spall, and recompression stages, a position-time (z-t) diagram in terms of local stress, σ zz , as an example, is given in Figure 1a. Under shock loading, a compressive shock wave was produced in the sample that travels towards the rear surface. At 25 ps, the unloading process began, generating a release wave from the side of piston. When the shock wave reached the rear surface (~30 ps), the rear surface expanded freely, yielding a reflected rarefaction wave. This rarefaction wave continued to propagate to the left and then interact with the tail of the release wave at about 40 ps to induce an evolving tensile region. When the tensile stress was strong enough, this interaction generated a high tri-axial tensile wave, as shown by the intersecting dashed-black arrows, which further lead to nucleation, growth, and collapse of voids, eventually forming a spall region. The peak value of the tensile stress at the spall plane is defined as the spall strength (σ sp ) of the material. We calculated the spall strength of each vacancy concentration system according to the maximal tensile stress during tensile. The method to determine the accurate spall strength is displayed in Figure 1b. We plotted the spatial distribution diagram of stress before and after spall damage nucleation, finding that the tensile stress reached its peak at 45 ps and the corresponding tensile stress was 3.75 GPa. This maximal tensile stress is called spall strength. The formation of spallation would induce recompression in the sample, resulting in a pullback hump in the free surface velocity history ( Figure 2). wave. This rarefaction wave continued to propagate to the left and then interact with the tail of the release wave at about 40 ps to induce an evolving tensile region. When the tensile stress was strong enough, this interaction generated a high tri-axial tensile wave, as shown by the intersecting dashed-black arrows, which further lead to nucleation, growth, and collapse of voids, eventually forming a spall region. The peak value of the tensile stress at the spall plane is defined as the spall strength (σsp) of the material. We calculated the spall strength of each vacancy concentration system according to the maximal tensile stress during tensile. The method to determine the accurate spall strength is displayed in Figure 1b. We plotted the spatial distribution diagram of stress before and after spall damage nucleation, finding that the tensile stress reached its peak at 45 ps and the corresponding tensile stress was 3.75 GPa. This maximal tensile stress is called spall strength. The formation of spallation would induce recompression in the sample, resulting in a pullback hump in the free surface velocity history ( Figure 2).  The free surface velocity history can be used to describe the various stages of spallation and measure the spall strength. Figure 2 shows the free surface velocity Vfs (t) evolution results along the two orientations at different Cv. When the shock front reaches the wave. This rarefaction wave continued to propagate to the left and then interact with the tail of the release wave at about 40 ps to induce an evolving tensile region. When the tensile stress was strong enough, this interaction generated a high tri-axial tensile wave, as shown by the intersecting dashed-black arrows, which further lead to nucleation, growth, and collapse of voids, eventually forming a spall region. The peak value of the tensile stress at the spall plane is defined as the spall strength (σsp) of the material. We calculated the spall strength of each vacancy concentration system according to the maximal tensile stress during tensile. The method to determine the accurate spall strength is displayed in Figure 1b. We plotted the spatial distribution diagram of stress before and after spall damage nucleation, finding that the tensile stress reached its peak at 45 ps and the corresponding tensile stress was 3.75 GPa. This maximal tensile stress is called spall strength. The formation of spallation would induce recompression in the sample, resulting in a pullback hump in the free surface velocity history ( Figure 2).  The free surface velocity history can be used to describe the various stages of spallation and measure the spall strength. Figure 2 shows the free surface velocity Vfs (t) evolution results along the two orientations at different Cv. When the shock front reaches the The free surface velocity history can be used to describe the various stages of spallation and measure the spall strength. Figure 2 shows the free surface velocity V fs (t) evolution results along the two orientations at different C v . When the shock front reaches the free surface, the free surface undergoes an acceleration from zero to about twice the velocity of the piston (V max ≈ 2 V p ). As the tensile state begins, the free surface velocity gradually decelerates to the first minimum value of V fs (t). The pullback velocity ∆V fs is the difference between the maximum and the first minimum value. When the spallation occurs, the tensile stress in the damaged area suddenly drops to zero, and then a compression wave forms in the spall plane. As a result, the reverberation of stress waves causes the formation of pullback hump (spall signal).
We can find that the effect of vacancy defects on V fs (t) is highly dependent on the crystal orientations. For the [0001] direction (Figure 2a), there was no difference in dynamic response when C v ≤ 0.5%, and no pullback hump appeared. Spallation occurred for C v > 0.5%, and the sharper pullback hump with increasing C v in V fs (t), indicating that the degree of spall damage increased gradually. However, for the [10-10] direction (Figure 2b), the unstable shock plateaus were formed during the release stage in the case of C v > 0.5%, which means that the vacancies collapsed, and plasticity may have absorbed part of the shock energy during the previous compression stage. No spallation occurred at C v = 0% because of its low tensile stress level. Interestingly, the pullback hump started to decline when C v increased to 0.5%, implying a gradual decrease in the degree of spall damage.
To further explain the anisotropic effect of C v on spall damage, some parameters of simulation results are listed in Table 1. The temperature of the spall plane, at the moment right before the incipient spall, is defined as spall temperature, T sp . The shock temperature, T H , is gradually released to the spall temperature, T sp , during adiabatic expansions. Vacancies collapse and plasticity can absorb part of shock energy, resulting in shock heating state. For [0001] direction, the shock stress, σ zz,H , and shock temperature, T H , remain nearly constant even at C v = 2.0%, while the spall temperature, T sp , increases by 7.5% for C v = 2.0% relative to C v ≤ 1.0% due to the heating effect of tensile plasticity. It is noted that the spall strength, σ sp , decreases with the increase in T sp when C v ≥ 1.0%, so that the spall damage is more likely to occur. However, for  direction, the σ zz,H reduces by about 16% for  direction when C v is increased to 2.0%. The stress attenuation will weaken the stress field around the spallation region, making it difficult to achieve the tensile stress required for void nucleation, so that the spall damage no longer easily occurs. Obviously, it is reasonable that the spall damage would reduce at higher C v for the  direction. Therefore, spallation damage may be determined by the competition mechanism between energy absorption and stress attenuation generated by plasticity. Table 1. Shock and spall parameters of Mg with different initial vacancy concentrations C v (%). The steady shock or Hugoniot state is denoted by the subscript H. Values listed are the shock stress, σ zz,H (GPa), the shock temperature, T H (K), the spall temperature, T sp (K) and the spall strength, σ sp (GPa).

Directions
C

Microstructures under Compression and Tension
Vacancies have an appreciable effect on the plasticity and spall damage of singlecrystal Mg if C v is above a critical value, and the effects highly depend on the shock loading direction, as shown by Figures 3 and 4. For shock along the [0001] direction, only elastic deformation occurred during compression for C v ≤ 1.0% (Figure 3a-c). When C v = 2.0% (Figure 3d), a small number of partial dislocations nucleated and propagated in the compression region of the sample. During tension, no spallation occurred at C v ≤ 0.5%. However, it is interesting that spallation still occurred at C v = 1.0% despite the absence of compression plasticity, which also occurred in the direction of  with C v = 0.5% (Figure 4b,f). This phenomenon is same as that in Cu [19,20] and Ni [21], implying that tension plasticity plays a key role in spall damage. For shock along the [10-10] direction, severe plastic deformation was found in the sample of C v ≥ 1.0% (Figure 4c,d) under compression, which mainly manifested as crystallographic reorientation, phase transition, stacking faults, and slips. The compression-induced plasticity will be partially or completely reversed during the release process, and plays an important role in the subsequent damage development. During tension, spallation occurred at C v ≥ 0.5% (Figure 4f,h), and the degree of spall damage decreased with the increase in C v , which is consistent with the degree of spall damage obtained from the evolution history of free surface velocity.
Metals 2022, 12, x FOR PEER REVIEW 6 of 13 deformation occurred during compression for Cv ≤ 1.0% (Figure 3a-c). When Cv = 2.0% (Figure 3d), a small number of partial dislocations nucleated and propagated in the compression region of the sample. During tension, no spallation occurred at Cv ≤ 0.5%. However, it is interesting that spallation still occurred at Cv = 1.0% despite the absence of compression plasticity, which also occurred in the direction of  with Cv = 0.5% ( Figure  4b,f). This phenomenon is same as that in Cu [19,20] and Ni [21], implying that tension plasticity plays a key role in spall damage. For shock along the [10-10] direction, severe plastic deformation was found in the sample of Cv ≥ 1.0% (Figure 4c,d) under compression, which mainly manifested as crystallographic reorientation, phase transition, stacking faults, and slips. The compression-induced plasticity will be partially or completely reversed during the release process, and plays an important role in the subsequent damage development. During tension, spallation occurred at Cv ≥ 0.5% (Figure 4f,h), and the degree of spall damage decreased with the increase in Cv, which is consistent with the degree of spall damage obtained from the evolution history of free surface velocity.  As we all know, the fracture mechanism of spall fracture in ductile metals is dominated by nucleation, growth, and coalescence of voids [40,41]. Figures 5 and 6 display the void evolution of magnesium with different Cv shocked along the [0001] direction. In order to identify the dislocations, the DXA was applied to track the trajectories of dislocation slip. For sample with Cv = 1.0%, the plastic deformation process at the initial spallation direction under V p = 0.5 km/s. Atoms are colored by CNA values, red atoms correspond to vacancies as well as disordered atoms, blue and yellow atoms correspond to the HCP structure and BCC structure, respectively, while the green atoms correspond to stacking faults (FCC structure), and the pink atoms represent the surface of void (the following atomic colors represent the same structure).
Metals 2022, 12, x FOR PEER REVIEW 6 of 13 deformation occurred during compression for Cv ≤ 1.0% (Figure 3a-c). When Cv = 2.0% (Figure 3d), a small number of partial dislocations nucleated and propagated in the compression region of the sample. During tension, no spallation occurred at Cv ≤ 0.5%. However, it is interesting that spallation still occurred at Cv = 1.0% despite the absence of compression plasticity, which also occurred in the direction of  with Cv = 0.5% ( Figure  4b,f). This phenomenon is same as that in Cu [19,20] and Ni [21], implying that tension plasticity plays a key role in spall damage. For shock along the [10-10] direction, severe plastic deformation was found in the sample of Cv ≥ 1.0% (Figure 4c,d) under compression, which mainly manifested as crystallographic reorientation, phase transition, stacking faults, and slips. The compression-induced plasticity will be partially or completely reversed during the release process, and plays an important role in the subsequent damage development. During tension, spallation occurred at Cv ≥ 0.5% (Figure 4f,h), and the degree of spall damage decreased with the increase in Cv, which is consistent with the degree of spall damage obtained from the evolution history of free surface velocity.  As we all know, the fracture mechanism of spall fracture in ductile metals is dominated by nucleation, growth, and coalescence of voids [40,41]. Figures 5 and 6 display the void evolution of magnesium with different Cv shocked along the [0001] direction. In order to identify the dislocations, the DXA was applied to track the trajectories of dislocation slip. For sample with Cv = 1.0%, the plastic deformation process at the initial spallation As we all know, the fracture mechanism of spall fracture in ductile metals is dominated by nucleation, growth, and coalescence of voids [40,41]. Figures 5 and 6 display the void evolution of magnesium with different C v shocked along the [0001] direction. In order to identify the dislocations, the DXA was applied to track the trajectories of dislocation slip. For sample with C v = 1.0%, the plastic deformation process at the initial spallation stage manifested in the pattern of the nucleation, growth, and multiplication of dislocation emitted out from vacancy region, as shown in Figure 5. The local tensile stress reached the maximum value at 45 ps. However, it is not until 53 ps that we clearly observed the formation of the vacancy cluster defects. Meanwhile, 2 leading partial dislocations (1/6 [0-22-3] and 1/3 [-1100]) were emitted from the vacancy region, which can coordinate the deformation in the C-axis direction. The slip of dislocations merged the surrounding vacancies to form a strip-like vacancy cluster defect microstructure. Subsequently, the generated vacancy cluster defect microstructure collapsed and eventually formed a void at about 56 ps. As the dislocation slipped, {11-22} <11-23> shear bands formed at the edges of void, and the void propagated along the direction of the shear bands. Basal stacking faults are observed near the shear bands with the slip of 1/3 <1-100> partial dislocations to accommodate the plastic strain. Due to the limited number of slip systems in the HCP structure, twinning deformation plays a crucial role in the plastic deformation of Mg [42]. At 60 ps, {11-21} <-1-126> twin appeared at the bottom left of the void. It may be that the prismatic and pyramidal slips are suppressed by periodic boundary conditions, and only twinning deformation is available near the voids [29]. The voids grow slowly due to the formation of the twinning which moves with the void tip. Therefore, it is difficult for the macroscopic fracture to occur. stage manifested in the pattern of the nucleation, growth, and multiplication of dislocation emitted out from vacancy region, as shown in Figure 5. The local tensile stress reached the maximum value at 45 ps. However, it is not until 53 ps that we clearly observed the formation of the vacancy cluster defects. Meanwhile, 2 leading partial dislocations (1/6 [0-22-3] and 1/3 [-1100]) were emitted from the vacancy region, which can coordinate the deformation in the C-axis direction. The slip of dislocations merged the surrounding vacancies to form a strip-like vacancy cluster defect microstructure. Subsequently, the generated vacancy cluster defect microstructure collapsed and eventually formed a void at about 56 ps. As the dislocation slipped, {11-22} <11-23> shear bands formed at the edges of void, and the void propagated along the direction of the shear bands. Basal stacking faults are observed near the shear bands with the slip of 1/3 <1-100> partial dislocations to accommodate the plastic strain. Due to the limited number of slip systems in the HCP structure, twinning deformation plays a crucial role in the plastic deformation of Mg [42]. At 60 ps, {11-21} <-1-126> twin appeared at the bottom left of the void. It may be that the prismatic and pyramidal slips are suppressed by periodic boundary conditions, and only twinning deformation is available near the voids [29]. The voids grow slowly due to the formation of the twinning which moves with the void tip. Therefore, it is difficult for the macroscopic fracture to occur.  For the sample with Cv = 2.0% (Figure 6), the plastic deformation pattern at the initial spallation stage was similar to the sample with Cv = 1.0%, with some key differences. At stage manifested in the pattern of the nucleation, growth, and multiplication of dislocation emitted out from vacancy region, as shown in Figure 5. The local tensile stress reached the maximum value at 45 ps. However, it is not until 53 ps that we clearly observed the formation of the vacancy cluster defects. Meanwhile, 2 leading partial dislocations (1/6 [0-22-3] and 1/3 [-1100]) were emitted from the vacancy region, which can coordinate the deformation in the C-axis direction. The slip of dislocations merged the surrounding vacancies to form a strip-like vacancy cluster defect microstructure. Subsequently, the generated vacancy cluster defect microstructure collapsed and eventually formed a void at about 56 ps. As the dislocation slipped, {11-22} <11-23> shear bands formed at the edges of void, and the void propagated along the direction of the shear bands. Basal stacking faults are observed near the shear bands with the slip of 1/3 <1-100> partial dislocations to accommodate the plastic strain. Due to the limited number of slip systems in the HCP structure, twinning deformation plays a crucial role in the plastic deformation of Mg [42]. At 60 ps, {11-21} <-1-126> twin appeared at the bottom left of the void. It may be that the prismatic and pyramidal slips are suppressed by periodic boundary conditions, and only twinning deformation is available near the voids [29]. The voids grow slowly due to the formation of the twinning which moves with the void tip. Therefore, it is difficult for the macroscopic fracture to occur.  For the sample with Cv = 2.0% (Figure 6), the plastic deformation pattern at the initial spallation stage was similar to the sample with Cv = 1.0%, with some key differences. At For the sample with C v = 2.0% (Figure 6), the plastic deformation pattern at the initial spallation stage was similar to the sample with C v = 1.0%, with some key differences. At about 41 ps, the reflected wave and the compression wave met, and then the compressioninduced plasticity gradually disappeared during tension. When the local tensile stress reached the maximum value at 45 ps, we observed that the tensile dislocations were emitted from the residual compression plasticity region. It means that the compression-induced plasticity has a significant effect on the nucleation sites of voids. Because of the weakening effects of local plastic deformation, the void formed early and its distribution was relatively concentrated for the sample with C v = 2.0%, so the macroscopic fracture occurred easily. Figure 7 displays the void evolution of magnesium with C v = 0.5% shocked along the  direction. Similar to the [0001] direction with C v = 1.0%, the dislocation nucleated preferentially at the vacancy cluster defects, and then the vacancy cluster defect microstructure collapsed to form a void. Due to the dislocation emission caused by the activation of two prismatic slip systems {1-100} <11-20>, the void was significantly blunted and almost stopped growing after a small amount of expansion. This phenomenon is similar to the simulation results of Tang et al. [29]. Subsequently, new voids were generated in the vacancies with the prismatic dislocation emissions at about 61 ps. Since the activation of prismatic slip systems were insufficient to release the concentrated stress, the dislocations in the basal and pyramidal planes near the voids were activated, as shown at 64 ps, followed by stacking faults. The nucleation and emission of dislocations eventually achieved the growth and coalescence of the voids. about 41 ps, the reflected wave and the compression wave met, and then the compressioninduced plasticity gradually disappeared during tension. When the local tensile stress reached the maximum value at 45 ps, we observed that the tensile dislocations were emitted from the residual compression plasticity region. It means that the compression-induced plasticity has a significant effect on the nucleation sites of voids. Because of the weakening effects of local plastic deformation, the void formed early and its distribution was relatively concentrated for the sample with Cv = 2.0%, so the macroscopic fracture occurred easily. Figure 7 displays the void evolution of magnesium with Cv = 0.5% shocked along the  direction. Similar to the [0001] direction with Cv = 1.0%, the dislocation nucleated preferentially at the vacancy cluster defects, and then the vacancy cluster defect microstructure collapsed to form a void. Due to the dislocation emission caused by the activation of two prismatic slip systems {1-100} <11-20>, the void was significantly blunted and almost stopped growing after a small amount of expansion. This phenomenon is similar to the simulation results of Tang et al. [29]. Subsequently, new voids were generated in the vacancies with the prismatic dislocation emissions at about 61 ps. Since the activation of prismatic slip systems were insufficient to release the concentrated stress, the dislocations in the basal and pyramidal planes near the voids were activated, as shown at 64 ps, followed by stacking faults. The nucleation and emission of dislocations eventually achieved the growth and coalescence of the voids. (b) Mechanism of void growth based on dislocation emission is analyzed using DXA: 1/3 <1-100> partial basal dislocations, 1/3 <1-210> partial prismatic dislocations, other dislocations (mainly 1/6 <20-23> and 1/9 <01-13> partial dislocations), etc., are presented.
For the samples with Cv = 1.0% loading along the [10-10] direction, crystallographic reorientation was the main plastic deformation mode during shock compression. We described the reorientation in more detail by selecting a representative region of the sample with Cv = 1.0%, as shown in Figure 8a. Crystallographic reorientation occurred around the [-12-10] orientation with a rotation angle of 90°, and the basal stacking fault appeared in the rotated region. The mechanism is discussed relatively clearly in a qualitative manner, displayed in Figure 8b. Similar to the atomic rearrangement of HCP single crystals, such as Mg [43][44][45], under the compression experiment and simulation, and Ti [46] under impact loading, the rotation was the result of a shear and shuffle mechanism. However, the lattice reorientation is reversible, and returned to the parent lattice gradually after unloading. In Figure 9, at 50 ps, the rotated lattice in the spall region gradually recovered to the parent lattice, and stacking faults can be observed on both sides of the grain boundary. When the tensile stress reached the maximum value at 53 ps, only the parent lattice exists in the spall region. There were intensive defects in the region of basal stacking faults interaction, which made the voids nucleate rapidly along the faults at a lower tension stress. As in the sample with Cv = 0.5%, voids nucleation and growth were accompanied by the (b) Mechanism of void growth based on dislocation emission is analyzed using DXA: 1/3 <1-100> partial basal dislocations, 1/3 <1-210> partial prismatic dislocations, other dislocations (mainly 1/6 <20-23> and 1/9 <01-13> partial dislocations), etc., are presented.
For the samples with C v = 1.0% loading along the [10-10] direction, crystallographic reorientation was the main plastic deformation mode during shock compression. We described the reorientation in more detail by selecting a representative region of the sample with C v = 1.0%, as shown in Figure 8a. Crystallographic reorientation occurred around the [-12-10] orientation with a rotation angle of 90 • , and the basal stacking fault appeared in the rotated region. The mechanism is discussed relatively clearly in a qualitative manner, displayed in Figure 8b. Similar to the atomic rearrangement of HCP single crystals, such as Mg [43][44][45], under the compression experiment and simulation, and Ti [46] under impact loading, the rotation was the result of a shear and shuffle mechanism. However, the lattice reorientation is reversible, and returned to the parent lattice gradually after unloading. In Figure 9, at 50 ps, the rotated lattice in the spall region gradually recovered to the parent lattice, and stacking faults can be observed on both sides of the grain boundary. When the tensile stress reached the maximum value at 53 ps, only the parent lattice exists in the spall region. There were intensive defects in the region of basal stacking faults interaction, which made the voids nucleate rapidly along the faults at a lower tension stress. As in the sample with C v = 0.5%, voids nucleation and growth were accompanied by the emission of {1-100} <11-20> prismatic dislocations, so the void shapes were non-spherical for shock along the [10-10] direction. Basal and pyramidal dislocations around the void were also activated, Metals 2022, 12, 215 9 of 13 but the density was much lower than that in the sample with C v = 0.5%. The reason is that the rapid growth of the voids greatly released the concentrated stress in the spall region. In addition, no twinning was observed in any simulations along the  direction. emission of {1-100} <11-20> prismatic dislocations, so the void shapes were non-spherical for shock along the  direction. Basal and pyramidal dislocations around the void were also activated, but the density was much lower than that in the sample with Cv = 0.5%. The reason is that the rapid growth of the voids greatly released the concentrated stress in the spall region. In addition, no twinning was observed in any simulations along the  direction.  By visualization, dislocation emission, void nucleation, and coalescence of voids and micro-cavities lead to the formation of spall plane. Void nucleates at vacancies, and compressive plasticity has an important effect on void nucleation sites.

Effect of Shock Velocity on Spall Strength
The spall strength is an important parameter for estimating the resistance to spall damage in a metal. In this work, the effect of vacancies on spall strength of single-crystal magnesium at different shock velocity Vp was also investigated, and the results are shown in Figure 10a,b. For both the [0001] and  directions, σsp generally decreased with the increase in Vp, which was consistent with Agarwal's MD results [14]. However, the pre-existence of vacancies further reduced the σsp of single-crystal Mg, making the spall occur at a lower level of tension stress. It is well known that the wave structure is a single elastic wave when a given loading strength is lower than the Hugoniot elastic limit (HEL) of a perfect single crystal [14]. In this work, the Hugoniot elastic limit velocities (VpHEL) of emission of {1-100} <11-20> prismatic dislocations, so the void shapes were non-spherical for shock along the  direction. Basal and pyramidal dislocations around the void were also activated, but the density was much lower than that in the sample with Cv = 0.5%. The reason is that the rapid growth of the voids greatly released the concentrated stress in the spall region. In addition, no twinning was observed in any simulations along the  direction.  By visualization, dislocation emission, void nucleation, and coalescence of voids and micro-cavities lead to the formation of spall plane. Void nucleates at vacancies, and compressive plasticity has an important effect on void nucleation sites.

Effect of Shock Velocity on Spall Strength
The spall strength is an important parameter for estimating the resistance to spall damage in a metal. In this work, the effect of vacancies on spall strength of single-crystal magnesium at different shock velocity Vp was also investigated, and the results are shown in Figure 10a,b. For both the [0001] and  directions, σsp generally decreased with the increase in Vp, which was consistent with Agarwal's MD results [14]. However, the pre-existence of vacancies further reduced the σsp of single-crystal Mg, making the spall occur at a lower level of tension stress. It is well known that the wave structure is a single elastic wave when a given loading strength is lower than the Hugoniot elastic limit (HEL) of a perfect single crystal [14]. In this work, the Hugoniot elastic limit velocities (VpHEL) of (b) Mechanism of void growth based on dislocation emission is analyzed using DXA: 1/3 <1-100> partial basal dislocations, 1/3 <1-210> partial prismatic dislocations, other dislocations (mainly 1/6 <20-23> and 1/9 <01-13> partial dislocations), etc., are presented.
By visualization, dislocation emission, void nucleation, and coalescence of voids and micro-cavities lead to the formation of spall plane. Void nucleates at vacancies, and compressive plasticity has an important effect on void nucleation sites.

Effect of Shock Velocity on Spall Strength
The spall strength is an important parameter for estimating the resistance to spall damage in a metal. In this work, the effect of vacancies on spall strength of single-crystal magnesium at different shock velocity V p was also investigated, and the results are shown in Figure 10a,b. For both the [0001] and  directions, σ sp generally decreased with the increase in V p , which was consistent with Agarwal's MD results [14]. However, the pre-existence of vacancies further reduced the σ sp of single-crystal Mg, making the spall occur at a lower level of tension stress. It is well known that the wave structure is a single elastic wave when a given loading strength is lower than the Hugoniot elastic limit (HEL) of a perfect single crystal [14]. In this work, the Hugoniot elastic limit velocities (V pHEL ) of perfect single-crystalline Mg in the [0001] direction and the  direction were about 1.1 km/s and 0.6 km/s, respectively. Vacancy defects promoted the nucleation of compression plasticity, which made the V pHEL in the defective sample lower than that in perfect single crystal. It can be found that when the impact velocity is slightly greater than V pHEL , the spall strength decreases significantly. This means that the plastic deformation in the compression stage has a significant effect on σ sp . For the [0001] direction, vacancy defects were found to have negligible effect on the spall strength of single-crystal Mg at higher shock velocity (V p ≥ 1.0 km/s). A similar trend can also be observed in Li's simulations for small size void magnesium [28]. Experimentally, the values of σ sp for magnesium and magnesium alloys reported were generally in the range of 1-2 GPa [2,[47][48][49][50], and the values were higher in the direction of  than that in the direction of [0001]. Our simulation results are close to the experimental results at higher shock velocities. Moreover, it can be observed that the fracture toughness of magnesium metal is much lower than that of BCC or FCC metals, owing to the very limited slip systems available for dislocation glide.
perfect single-crystalline Mg in the [0001] direction and the  direction were about 1.1 km/s and 0.6 km/s, respectively. Vacancy defects promoted the nucleation of compression plasticity, which made the VpHEL in the defective sample lower than that in perfect single crystal. It can be found that when the impact velocity is slightly greater than VpHEL, the spall strength decreases significantly. This means that the plastic deformation in the compression stage has a significant effect on σsp. For the [0001] direction, vacancy defects were found to have negligible effect on the spall strength of single-crystal Mg at higher shock velocity (Vp ≥ 1.0 km/s). A similar trend can also be observed in Li's simulations for small size void magnesium [28]. Experimentally, the values of σsp for magnesium and magnesium alloys reported were generally in the range of 1-2 GPa [2,[47][48][49][50], and the values were higher in the direction of  than that in the direction of [0001]. Our simulation results are close to the experimental results at higher shock velocities. Moreover, it can be observed that the fracture toughness of magnesium metal is much lower than that of BCC or FCC metals, owing to the very limited slip systems available for dislocation glide. Consequently, the effect of vacancy defects on plastic deformation and spallation response was found to be anisotropic. Since the  direction had the lowest activated slip system, the effect of vacancy defects can be expected to be most significant in this direction.

Conclusions
Molecular dynamics simulations are used to study the effect of vacancies on dynamic response and spallation in single-crystal magnesium in the shock directions [0001] and . The mechanisms of plastic deformation, spall damage, and void evolution of single-crystal magnesium with various vacancy concentrations are investigated. The main findings are summarized as follows.
Vacancies provide nucleation sites for plasticity, and the increasing vacancy concentration greatly promotes the heterogeneous nucleation of plasticity. Our research shows that the plastic deformation is various for different loading directions. The shock compression stage exhibits the deformation modes of crystallographic reorientation, phase transition, stacking faults, and slips when Cv is above a critical value. The crystallographic reorientation can only take place in samples loaded along the  direction. Under tension, the basal, prismatic, and pyramidal slips occur at the spallation region and are responsible for the plasticity deformation. The compression plasticity has a significant effect on the nucleation sites, while the tensile plasticity plays a key role in the evolution of void shape. Consequently, the effect of vacancy defects on plastic deformation and spallation response was found to be anisotropic. Since the  direction had the lowest activated slip system, the effect of vacancy defects can be expected to be most significant in this direction.

Conclusions
Molecular dynamics simulations are used to study the effect of vacancies on dynamic response and spallation in single-crystal magnesium in the shock directions [0001] and . The mechanisms of plastic deformation, spall damage, and void evolution of single-crystal magnesium with various vacancy concentrations are investigated. The main findings are summarized as follows.
Vacancies provide nucleation sites for plasticity, and the increasing vacancy concentration greatly promotes the heterogeneous nucleation of plasticity. Our research shows that the plastic deformation is various for different loading directions. The shock compression stage exhibits the deformation modes of crystallographic reorientation, phase transition, stacking faults, and slips when C v is above a critical value. The crystallographic reorientation can only take place in samples loaded along the  direction. Under tension, the basal, prismatic, and pyramidal slips occur at the spallation region and are responsible for the plasticity deformation. The compression plasticity has a significant effect on the nucleation sites, while the tensile plasticity plays a key role in the evolution of void shape.
Vacancy collapse and subsequent plastic deformation mainly result in two competing mechanisms-energy absorption and stress attenuation. For the [0001] direction, energy absorption dominates, so the degree of spall damage increases with the increase in C v . For the  direction, when C v ≥ 0.5%, stress attenuation dominates, so the degree of spall damage decreases with the increase in C v . Whether spallation damage occurs or not, eventually depends on the competition between the two mechanisms.
By visualizing the evolution of spallation, we found that the voids usually nucleate at defects, and the void evolution is mainly based on the emission mechanism of dislocations. The {11-22} <11-23> pyramidal dislocation facilitates the nucleation of void in the [0001] shock direction, as well as the {1-100} <11-20> prismatic dislocation in the  shock direction. In addition, the {11-21} <-1-126> tensile twin is observed near the voids in the orientation of [0001], and no tensile twin participates in the deformation of the  direction.
Author Contributions: C.J. and Z.J. concept this work; C.J. performs simulations, collects data, analyzes data and writes the original draft; all authors analyze results; X.L. and S.X. contribute to writing-review and editing manuscript. All authors have read and agreed to the published version of the manuscript.