Spallation is a typical failure mode associated with the tensile stress generated by the interaction of a reflected rarefaction wave (RRW) and an unloading rarefaction wave (URW) when a material is exposed to intense impact.1 For ductile metals, this dynamic performance is characterized by void nucleation, growth, and coalescence,2 accompanied by elastic–plastic deformation, phase transition, melting, and other nonlinear behaviors.3,4 Spallation is usually categorized into classical spallation, multiple spallation, and micro-spallation.5 Although some investigations of multiple spallation with several spall planes have been reported,6–8 complex microstructural evolution occurs during multiple spallation, which renders it challenging to explore the formation micromechanism of this phenomenon.

Although significant progress has been made in the study of the spall mechanism in the case of flat-plate impact, existing techniques are not well suited to in situ experimental exploration of material shock performance during nonplanar loading (such as cylindrical and spherical impacts). Implosion is one of the typical impact types, in which a cylindrical converging detonation wave is induced to propagate in a cylindrical shell. In this case, materials always experience complex stress states, exhibiting a unique manner of fracture.9,10 Xue et al.11 conducted an implosion impact experiment on ductile metal. They reported that multiple adiabatic shear bands occurred in a thick-walled cylinder, which were initiated by momentum diffusion from stress unloading, stress/strain/temperature perturbations, and microstructural inhomogeneities. Shirinkina et al.12 also experimentally investigated multi-spall behavior under implosion loading and found that it involved intensive deformation, fragmentation, and melting. On the basis of these observed phenomena, Zhang et al.13 proposed a unified phase field theory to describe the complex failure coupled with multiple spallation and adiabatic shear banding in ductile metals under implosion, revealing the dynamic fracture mechanism of metal shells. Vishnu et al.14 performed finite element simulations for electromagnetically collapsing thick-walled cylinders, demonstrating that voids contributed to nucleation and the development of dynamic shear localization. Owing to the intricacy of the damage patterns and their transient nature, the formation mechanism of multiple spallation under implosion impact has rarely been investigated and requires simultaneous consideration of the interaction of multiple damage events. Moreover, there is a lack of mechanistic explanations of the penetration of multiple spall planes. Owing to their remarkable temporal and spatial resolution, molecular dynamics (MD) simulations have proved to be an effective way to obtain a deep understanding of the physical mechanism of dynamic damage combined with transient microstructural evolution at the atomic level.15–17 For instance, Chen et al.18 proposed a MD-based approach to explore the anisotropy of melting and spallation behavior of copper under implosion loading, highlighting that spallation occurred in premelting regions. Tan et al.19 observed twinning of hexagonal close-packed variants under cylindrical impact, revealing coupling between plastic slip and phase transitions at the atomic scale.

It has been shown that there is a strong dependence of mechanical properties on grain size in polycrystalline metals, involving a competitive mechanism of plastic flow initiation.20,21 Cordero et al.22 argued that a decrease in grain size facilitated dislocation blockage for grain-size strengthening; however, for the smallest sizes, a softening with decreasing grain size was observed, which was triggered by activation of grain boundary deformation. Clearly, spall performance possesses a strong sensitivity to grain size. Flanagan et al.23 pointed out that there was an inverse Hall–Petch relationship for spall strength at the micrometer scale, while the spall strength exhibited a Hall–Petch relationship with grain size as the size decreased to the nanoscale level, owing to excessive grain boundary content (up to 50%). Zhu et al.24 reported that there was a minimum spall strength when the grain size was in the ranged 10–14 nm, indicative of the existence of a critical size. In fact, as one of primary defects occurring in materials, grain boundaries have been found to dominate plastic deformation (including grain boundary sliding, migration, grain rotation, and dislocation motion).23,25 Bronkhorst et al.26 proposed that local stress profiles were affected by grain boundaries and triple junctions in tantalum polycrystals during shock loading. Fensin et al.27 found that voids were more likely to nucleate at grain boundaries, especially those perpendicular to the direction of impact. Furthermore, Minich et al.28 proposed that as the grain size decreased, an increase in grain boundary surface area provided a greater number of potential void nucleation sites, thereby decreasing the spall strength. As yet, however, the relationship between grain size and this dynamic failure during implosion loading remains unclear, and there is therefore an urgent need for a deeper understanding of the spallation of nanocrystalline metals from the perspective of the initial microstructures.

Among binary shape memory alloys, nickel–titanium (NiTi) alloys are promising candidates for structural components that are subjected to extremely high strain rates in engineering applications such as shielding against hypervelocity impacts in outer space, owing to their excellent properties such as shape memory effect, superelasticity, and good damping characteristics.29,30 There has been considerable interest in the dynamic response of NiTi alloys to intense impact. According to Nemat-Nasser et al.,31 the NiTi deformation mechanism is stress-induced martensitic transition at low strain rates, whereas dislocation-induced plastic slip of the austenite phase became the dominant factor determining deformation at high strain rates. Zhang et al.32 found that the dynamic behavior of Ni₅₂Ti₄₈ subjected to high strain rate could be attributed mainly to elastic–plastic transition, because there was almost no martensitic phase observed in recovered samples. Li et al.33 proposed that such a difference in NiTi dynamic response mechanisms was strongly dependent upon alloy components. Xi and Su34 even discovered that a grain size at the nanoscale led to disappearance of the martensitic transition due to grain boundary inhibition. As reported in previous work,35–37 there are significant differences in dynamic damage behavior between NiTi alloys and classical metals, because these alloys exhibit more complex microstructural evolution (involving, for example, dislocation motion and plastic response) under strong shocks. To date, very little work has been done on the dynamic behavior of cylindrical nanocrystalline NiTi alloys subjected to high strain rate loading, and so little information is available regarding complex dynamic failure modes coupled with multiple spallation and shear localization. However, such phenomena have been experimentally observed in the expanding fracture behavior of an explosively driven steel cylinder.10 Therefore, the physical mechanism underlying the dynamic behavior of NiTi alloys under extreme loading conditions deserves to be explored more comprehensively.

Hence, in the present study, atomistic simulations are conducted to investigate the dynamic behavior of a cylindrical shell composed of nanocrystalline NiTi alloy subjected to implosion loading. In contrast to other studies, the primary objective of this work is to explore the micromechanisms of multi-spall formation and spall plane penetration, paying particular attention to the interaction between spallation and shear localization. In conjunction with dislocation motion, void statistics are employed to reveal the effects of grain size on spall strength and void nucleation mode. These results are expected to provide novel insights into the nature of complicated dynamic failure of NiTi-based cylindrical shells, facilitating the further development of shape memory alloys for practical applications.

With the help of the Large-Scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) open source software,38 nonequilibrium molecular dynamics (NEMD) simulations were conducted to investigate the dynamic performance of nanocrystalline NiTi alloy subjected to implosion loading. The atomic interactions occurring in the material were described by the embedded atom method (EAM) potential proposed by Lai and Liu39 and improved by Zhong et al.,40 which has been confirmed as being suitable for the description of Ni₅₂Ti₄₈ dynamic performance.5 Voronoi tessellation41 was employed to generate several cylindrical models with average grain sizes in the range of 2–50 nm, in which the grains, with random orientations, were distributed in the form of nanopillars. The systems centered on the z axis contained ∼8 × 10⁶ atoms, with simulated dimensions of 300 nm external diameter × 200 nm internal diameter × 3 nm height. All the initial Ni₅₀Ti₅₀ models were set as the B2 austenite phase with lattice constant 3.008 Å. Subsequently, some of the Ti atoms were randomly replaced by Ni atoms to produce Ni₅₂Ti₄₈. Since atomic overlaps were present at grain boundaries in the systems, achieved through Voronoi tessellation,42 the conjugate gradient method was used for energy minimization to sufficiently stabilize atoms. The simulated system with periodic boundaries in all three directions (x, y, and z) was first relaxed for 30 ps under an isobaric–isothermal ensemble and eventually reached an equilibrium state with 300 K and zero pressure to eliminate the effect of residual stress. Any point in the polar coordinate system could be equivalently represented by a corresponding point in the rectangular coordinate system.

As presented in Fig. 1, implosion loading was accomplished by rapid shrinkage of a cylindrical potential wall of infinite height at a velocity of 1.0 km/s, in which atoms interacted with the wall, resulting in the formation of a cylindrical convergent shock wave. This method of loading has been proved to be effective for the introduction of a cylindrical wave.18,19 It is worth noting that the cylindrical potential wall had the same axis as the cylinder sample. Since there was a slight change in the external radius after relaxation, the radius of the potential wall was set to 151 nm so that all the atoms were kept inside the wall. When the atoms were forced along the radial direction by the nearest potential wall, the interaction energy E was given by the Lennard-Jones potential38$$E=4\epsilon \left[{\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right],r<r_{\text{c}},$$(1)where ɛ and σ determine the interaction strength and distance, respectively, r is the distance from an atom to the nearest potential wall, and r_c is the cutoff distance, beyond which there was no interaction between atoms and potential wall. After 10-ps loading, the system was unloaded by removing the potential wall. A microcanonical ensemble was applied to the whole simulated system so that energy transfer under high strain rate was ignored. To reduce boundary effects, the x and y axes were maintained with shrink-wrapped boundary conditions, while a periodic boundary condition was applied to the z axis. All the simulations were conducted using a time step of 1 fs.

The 2D Lagrangian binning method was employed to trace local stress, density, and temperature profiles in samples during progressive damage. The area of each bin was set as 1 × 1 nm² in the x–y plane. Thus, instantaneous spatial distributions of those physical quantities were obtained by averaging over all atoms in each bin. The stress tensor for each atom was given according to the virial definition43$${\mathbf{\sigma }}_{\text{ab}}=-m{\mathbf{v}}_{\mathrm{a}}{\mathbf{v}}_{\mathrm{b}}-{\mathbf{W}}_{\text{ab}},$$(2)where a and b take values x, y, z to generate the six components of the symmetric tensor, and m and v are the atomic mass and thermal vibration velocity of atoms, respectively. The first term on the right-hand side of Eq. (2) is the kinetic energy contribution for each atom. The second term is the virial contribution due to intra- and intermolecular interactions.

According to classical mechanics, the radial stress component σ_ρ and azimuthal stress component σ_φ of the entire bin are given by$${\sigma }_{\mathrm{\rho }}={\sigma }_{\text{xx}}{\cos }^2\theta +{\sigma }_{\text{yy}}{\sin }^2\theta +2{\sigma }_{\text{xy}}\sin \theta \cos \theta ,$$(3a)$${\sigma }_{\mathrm{\phi }}={\sigma }_{\text{xx}}{\sin }^2\theta +{\sigma }_{\text{yy}}{\cos }^2\theta -2{\sigma }_{\text{xy}}\sin \theta \cos \theta ,$$(3b)where σ_xx, σ_yy, and σ_zzare the mean stresses in the x, y, and z directions, respectively, σ_xy is the shear stress, and θ is the angle between the radial direction and the positive x axis. Positive and negative values of σ_ρ correspond to compressive and tensile states, respectively. Here, the spall strength has been taken to be given approximately by the maximum mean tensile stress, which is numerically equivalent to ${\left({\sigma }_{\mathrm{\rho }}\right)}_{\max }$. It must be emphasized that removal of the local center-of-mass velocity along the radial direction provides an accurate determination of stress and temperature during simulations.44 In addition, the atomic-level deformation gradient and the strain tensor at each particle from the relative motion of its neighbors were determined using a modifier from the OVITO visualization software.45

To clearly observe microstructural evolution, voids were identified by the construct surface mesh (CSM) method based on the alpha-shape algorithm,46 where the probing sphere radius was chosen as 3.6 Å in the coordinated analysis. The dislocation statistic was achieved through the dislocation extraction algorithm (DXA).47 In addition, the polyhedral template matching (PTM) method was used to identify atomic arrangements at high strain rate, because it possesses an identification capability with higher accuracy48 than common neighbor analysis. The cutoff value of the root-mean-square deviation used for PTM was selected as 0.18.48 On the basis of the PTM results, the grain size was determined through the grain segmentation method, which has an advantage over the grain tracking algorithm.49 All the visualizations under various configurations were analyzed using OVITO.

First, implosion loading was applied to the 50-nm sample to explore unique dynamic response of NiTi alloy. It is well known that spallation is characterized by a tensile failure caused by the interaction of a URW and RRW reflected from the material surface. Figure 2 plots the radial stress profiles at various time points to display stress wave propagation, in which compressive and tensile states are represented in red and blue, respectively. After a continuous loading of 8.5 ps, the shock wave reaches the inner surface and is reflected. Subsequently, it propagates toward the outer surface. In contrast to the state at 13 ps, we can see that a tensile stress is present at 14 ps, which can be attributed to the collision of the RRW with the URW. The tensile stress increases to the peak at 14.5 ps, indicative of void nucleation. As the damage process continues, the stress wave is decomposed into two parts, as shown at 16.5 ps, with one wave propagating toward the inner surface and the other outward. At 21.5 ps, cavitation regions (represented in white) can clearly be observed between the inner and outer surfaces, suggesting the formation of spallation in the material.

Figure 3 presents void evolutions during implosion loading. As observed in Fig. 3(a), the variation in void number exhibits a distinct bimodal trend. Considering Figs. 3 and 4 together, it can be seen that when the number of voids reaches its first peak at 16.5 ps, large voids are concentrated in those zones where the RRW and URW have interacted, indicative of the formation of the first spall plane. When the number of voids reaches its second peak at 21.5 ps, many individual voids are found to nucleate in a strip-like manner within the inner region [between the first spall plane and the inner surface, as marked by the red dashed line in Fig. 3(b)]. These voids are distributed in some lath-shaped zones and oriented at an angle to the radial direction. It can also be seen from Fig. 3(a) that new void nucleations occur in the inner region only after the total number of voids has reached its first peak. The void nucleation in the inner region promotes the second peak of the total number of voids, indicating that multistage void nucleation processes occur during spallation.

In Fig. 5(a), a specific position is chosen to provide a better observation of the spalling process. Figure 5(b) reveals several regions with lower density (represented in blue), suggesting the presence of multiple spallation. Figure 5(c) presents a time–position–radial stress plot to facilitate exploration of the mechanism of multiple spallation. It is clear that the tensile stress is generated by interaction of the RRW reflected from the inner surface with the URW. Chen et al.18 also conducted MD simulations to investigate the spall failure of single-crystal copper under cylindrical converging impact, and they found that the wave reflected from the outer surface interacted with the unloading wave, resulting in a tensile stress state. In the present study, some stress-release regions can be identified in the material, corresponding to low-density zones. One is characterized by the first spall plane at the maximum tensile stress, with ∼10.6 GPa being generated in the sample at 14.5 ps. As mentioned earlier, the spall strength is taken to be approximately given by the maximum tensile stress. Other regions are located near the inner surface along the propagation of the tensile stress wave and corresponding to the secondary spall plane. As shown in Fig. 5(d), there is another tensile stress peak corresponding to about 6.5 GPa at 17 ps, at which time new voids have nucleated in the inside region, indicating that the secondary spallation can still be regarded as a dynamic tensile failure. The significant difference between the two stress peaks indicates that a lower maximum tensile stress is responsible for the generation of the second spall plane, but some influential factors contributing to the reduced spall strength need to be further studied.

In fact, a cylindrical sample shrinks inward during implosion loading, and so the whole system is exposed to a more complex stress state in comparison with the case of plate shock. Figure 6 shows shear strain profiles at various time points. Initially (5 ps), before the shock wave has reached the internal free surface, the shear strain is primarily concentrated at grain boundaries. In polycrystalline materials, when adjacent grains are subjected to shear stress, plastic deformation such as grain boundary sliding occurs, leading to shear strain localization.50 In the present work, we find that the inner surface begins to shrink inward when the shock wave reaches it (8.5 ps), and this is followed by shear-deformation enrichment in the inner part. Subsequently, some stripe-like zones covered by high shear strain form at the inner surface and grow toward the outside, which is consistent with experimental observations of metal-based cylindrical shells under explosive loading.51

Figure 7(a) plots the evolution of the azimuthal stress component, corresponding to Fig. 5(d). From the profiles of azimuthal stress and shear strain in Figs. 7(b) and 7(c), it can be seen that after the formation of the first spallation (15 ps), the high-shear-strain zones are occupied by atoms with high azimuthal stress (represented in red). Subsequently, the azimuthal stress propagates toward the inner and outer surfaces, centered on the spall plane. At 17 ps, two high-azimuthal-stress rings (represented in red) are observed in the cylindrical sample. Hence, it is believed that azimuthal stress contributes to the presence of high shear strain.

In the present case, the deformation is at a high rate, and so, following shell convergence, these localized deformation zones become pronounced. Here, such events are taken to be shear deformation bands (SDBs). From the shear strain and dislocation profiles in Fig. 8, it can be seen that over time, dislocations are activated and continuously assemble within the deformed region, prompting a shearing effect.52 This suggests that the survival of shear localization is closely associated with plastic deformation, and dislocation activity is conducive to the origin and propagation of SDBs. As can also be seen from Fig. 8, there is a significant increase in temperature during shear localization. This is consistent with experimental findings by Zhu et al.,53 who showed that temperature increase was not primarily responsible for the initiation of shear banding in metal-based samples, and, conversely, shear instability induced by microstructural effects may be the trigger of shear band formation. This viewpoint is consistent with crystal plasticity modeling carried out by Mayer and co-workers,54,55 who reported that perturbations of dislocation density were more effective than temperature perturbations in causing plastic localization.

More SDBs are distributed in the inner region at 16.5 ps compared with 14.5 ps, as presented in Fig. 6. Voids belonging to the second spall plane then nucleate at 17 ps. Figures 4 and 6 reveal that voids in the inner region are concentrated along the bands, and thus the SDBs can be expected to play a role in subsequent secondary spallation.

Figure 9 presents distribution patterns of voids and SDBs at 21.5 ps, with initial grain boundaries and voids colored in black and orange, respectively. It should be mentioned that the atoms with shear strain exceeding 0.5 (colored in yellow) in the inner region are retained while atoms that do not satisfy this condition are deleted. As can be seen, all the voids in the inner region have nucleated at positions of high shear strain. Careful examination reveals that some of the void nucleations have occurred on initial grain boundaries, while the other voids have nucleated on the shear-strain-accumulation regions inside the grains. In contrast to the first spall plane, the formation of the secondary spall plane is most likely to be closely associated with SDBs. Similarly, previous work56 has revealed that multiple spallation in typical metals is strongly governed by shear failure.

In ductile metals, the formation of shear bands under high strain rates is usually accompanied by a localized adiabatic temperature rise, which is responsible for a decrease in spall strength.5,13 Evidently, temperature plays a vital role in affecting multi-spall performance. Figure 10 presents temperature profiles at various time points. Overall, the higher temperatures are distributed on grain boundaries with higher shear strain, as compared with other regions. At the first spallation, the temperature does not significantly increase, and in fact drops before void nucleation, owing to the presence of tensile stress.57 Subsequently, the temperature increases as voids nucleate and become agglomerated. At the secondary spallation, there is a localized temperature rise in the strip zones where the SDBs are born. Such high-temperature regions were in fact present before the formation of new voids (16.5 ps), as can be seen in Fig. 10.

Figure 11 compares shear strain and temperature profiles at 16.5 ps for a circular region with a radius of 95 nm. According to Fig. 5(d), this radius is determined as the median between the inner surface (80 nm) and the position of the tensile stress peak (110 nm). Here, the x axis is at 0° and the counterclockwise direction is taken to be positive.

To further explore the dependence of spall strength on temperature, Table I presents the spall strength of a NiTi sample at various initial temperatures. It can be seen that the spall strength for the first spallation is close to the value predicted for planar impact.5 These results demonstrate that an increase in initial temperature leads to a degradation of spall strength, causing a softening effect. It is believed that along with dislocation activities, SDBs are produced in the secondary spall region, accompanied by a localized temperature rise. Thus, high temperatures induced by shear localization promote material softening, contributing to a decrease in spall strength. This is the reason why the maximum tensile stress for the secondary spallation is lower than that for the first one. Figure 12 presents an enlarged view of the temperature profile of the simulated system at 16.5 ps. Clearly, there are some high-temperature zones distributed on the grain boundaries and SDBs, in which the maximum temperature is mostly beyond 1000 K. In fact, a higher temperature is most likely to bring about material melting, leading to a more significant decrease in spall strength. Accordingly, the spall strength can be considered to be sensitive to temperature.

Figures 13 and 14 display the shock response of the NiTi-based shell during the formation of spallation under different loading velocities. Multiple dynamic events, such as multiple spallation and shear localization, are still pronounced for a loading velocity of 0.75 km/s (Fig. 13). Before the formation of the first spallation (15 ps), significant shear strain concentrations are present in the material. At 18 ps, the first spallation has been generated, accompanied by more severe shear localization induced by plastic deformation. Along with the formation of secondary spallation (corresponding to 21 and 24 ps, respectively), a large number of SDBs have nucleated and grown within the sample, leading to a wider distribution zone. Importantly, voids are seen to nucleate on SDBs, ultimately facilitating the occurrence of multiple spallation. When the loading velocity is decreased to 0.5 km/s (Fig. 14), there is a weakening of the severe localized shearing effect in the sample. In this case, fewer SDBs are observed than for a loading velocity of 0.75 km/s. Thus, fewer void nucleation sites are provided in the material, and so it is difficult to produce multiple spallation. These findings reveal that the dynamic damage performance is sensitive to loading velocity, conforming to experimental observations for other alloys.58,59 It is most likely that there is a critical loading velocity responsible for activating complex dynamic failure modes in NiTi alloy subjected to implosion loading; however, to simplify the analysis, the influence of loading velocity on NiTi dynamic performance has not been systematically investigated in this work,

Figure 15(a) displays the void evolution at several average grain sizes to reveal the formation of multiple spallation. For these cases, two peaks in the number of voids are observed, indicating multiple void-nucleation. Figures 15(b)–15(d) are time–position–density plots at the position shown in Fig. 5(a). The atomic configurations colored according to shear strain in the range from 0 to 1 at 30 ps are also shown on the right side of each plot. For grain sizes of 5 and 10 nm [Figs. 15(b) and 15(c)], it can be seen that the material is divided by multiple blue regions with low density, corresponding to the cavitation zone in atomic configurations. The blue region closest to the outer surface represents the first spall plane, in which the first void-nucleation sites are located. Similarly, the secondary spallation corresponds to the remaining blue zones. At a grain size of 30 nm [Fig. 15(d)], there are no distinct multiple blue regions, but plenty of voids are distributed between the first spall plane and the inner surface, promoting the occurrence of the secondary spallation.

Figure 16 shows the dependence of spall strength on grain size for multiple spallation. For the first spallation, as the average grain size decreases, the spall strength undergoes a transition from an inverse Hall–Petch to a Hall–Petch relationship at a critical size of 10 nm. At larger grain sizes (>10 nm), the content of grain boundaries decreases with increasing grain size, providing fewer void-nucleated sites. As a consequence, a higher spall strength is produced. Regarding the Hall–Petch relationship for grain sizes in the range from 10 to 2 nm, a reverse increasing trend for the spall strength is observed in the NiTi alloy. This suggests that an excessive number of grain boundaries enhance superelasticity, leading to superior spall strength.24 For the second spallation, a similar transition is observed at 8 nm. By contrast, as the average grain size increases from 8 to 50 nm, there is another transition in spall strength at a grain size of 30 nm. Compared with Fig. 15(a), the corresponding void number also shows a transformation from an increasing to a decreasing trend. This is likely to be related to shear localization, which will be discussed in the subsequent analysis.

Figure 17 shows the void and SDB distributions in the inner region for various grain sizes. In the case of a grain size of 5 nm, for the first spallation, all the voids have nucleated along the initial grain boundaries. For nanocrystals, it has been shown that as grain boundaries subjected to shear strain become weak regions,23 voids preferentially nucleate along them rather than in grains. The shear strain distributions before the first void nucleation (14 ps) for various grain sizes are shown in Fig. 18. Clearly, for a grain size of 5 nm, a large number of grain boundaries with high shear strain are densely distributed, providing enough potential sites for the first void nucleation. In this case, the dominant mechanism of plastic flow in the system is grain boundary deformation. As the grain size is increased from 10 to 50 nm, there are insufficient nucleation sites on the grain boundaries, and consequently the voids nucleate in a coexisting intergranular/transgranular manner. As a consequence, there is an increase in void number with increasing grain size, which is similar to what happens in the case of one-dimensional planar impact.35 Furthermore, an improved spall strength is expected, because fewer voids nucleate and become aggregated at grain boundaries.

For the secondary spallation, voids individually nucleate on the initial grain boundaries when the grain size is small (5 nm). Additionally, the high-shear-strain atoms cluster near the grain boundaries, and no SDBs are formed in the grains. This behavior is completely different to that found for a grain size of 50 nm (Fig. 9), which is to be expected, since the grain boundary sliding becomes the dominant mechanism of plastic deformation at smaller grain sizes.60 At a grain size of 10 nm, there are almost no voids nucleating within grains, in spite of the presence of some generated SDBs. As shown in Figs. 18 and 19, the shear strain in these bands is not as high as that at grain boundaries, and intergranular nucleation is maintained. In addition, the void distribution for secondary nucleation exhibits a strip pattern along the grain boundaries, in contrast to the scattered distribution at 5 nm. This is the reason for the greater number of voids for a grain size of 10 nm, compared with 5 nm. As the grain size is increased to 30 nm and then to 50 nm, a greater number of voids nucleating along SDBs are distributed with a strip-like manner in the grains. In this case, shear strain concentrations at fewer grain boundaries cannot be completely accommodated, and so grain boundary sliding loses its dominance of plastic deformation, and in turn a large quantity of SDBs are formed within grains. In other words, shear deformations became strongly localized in the larger grains, demonstrating the sensitivity of plastic deformation to grain size. These results deduced from the MD modeling reveal the correlation of the development of shear localization with initial grain size. However, owing to rapid aggregation induced by high void concentrations at grain boundaries, the peak number of voids for a grain size of 50 nm is lower than that for 30 nm.

Figure 20 displays atomic configurations colored by PTM at the initial time and after secondary void nucleation. It can be seen that for a grain size of 5 nm, the grain boundaries became thicker after secondary nucleation, because those atoms located along the boundaries are converted to disordered structures by sliding. At larger grain sizes (especially for 30 and 50 nm), the initiation and propagation of SDBs leads to a localized disorder of atoms within grains, and so new grain boundaries are produced. As consequence, a single grain is divided into multiple parts.

On the basis of the PTM results, grain sizes before and after void nucleation are counted by the grain segmentation method, as displayed in Fig. 21. The ordinate represents the number of atoms in a grain, which is equivalent to grain size. The grains are arranged and labeled according to size from large to small. It is evident that as the grain size increases, the number of grains in the material grows significantly after nucleation, promoting grain refinement. A similar behavior was also captured by experiments on other ductile metals,6 demonstrating that smaller grain size in multiple spallation is likely to be caused by high shear localization under shock loading followed by heating. Thus, grain refinement is expected to be closely associated with the presence of SDBs.

Figure 22 shows the evolution of dislocation length for various grain sizes, demonstrating that an increase in grain size leads to a decrease in initial total dislocation length. This is because for larger grain size, there are fewer initial grain boundaries to provide primary places for accumulation of dislocations. During the loading stage, the 1/2⟨111⟩ dislocation length remains almost unchanged for a grain size of 5 nm, whereas a significant increase is found at other grain sizes, demonstrating that there is a unique plastic deformation mechanism for 5 nm. In fact, after void nucleation, nearly all the dislocations are annihilated (but do not disappear). Similar behavior was observed in the case of piston impact,35 in which a large number of ⟨100⟩ dislocations were found to be generated and accumulated after the peak in void number, leading to local strengthening for larger grain sizes. By contrast, for implosion loading, only a few ⟨100⟩ dislocations are generated, even at larger grain sizes. Importantly, the distribution of such dislocations is concentrated almost solely in the region between the secondary spall plane and the outer surface, and they are almost absent from the inner region, as can be seen in Fig. 23. This can be attributed to the presence of SDBs producing many refined grains with smaller sizes in the inner region and thus preventing the proliferation of ⟨100⟩ dislocations. For the secondary spallation, grain refinement in the sample with a grain size of 50 nm, which is a consequence of the formation of new grain boundaries, results in a lower spall strength than in the case of a 30-nm grain size, indicative of a reduced dislocation-induced strengthening effect.

Figures 4 and 9(b) reveal a large number of voids located on a specific grain boundary, causing a through penetration between the first and secondary spall planes. As already mentioned, there is a higher shear strain as well as a higher temperature on the linking grain boundary, compared with other positions. To explore the penetration mechanism, the Simon function determined by the two-phase method is employed as follows to evaluate the melting state of the material:5,61$$T_{\text{m}}\left(P\right)=1226{\left(1+\frac{P}{7.226}\right)}^{0.368},$$(4)where T_m(P) is the melting point at pressure P.

The evolution of the temperature state at the position corresponding to Fig. 9(b) is displayed in Fig. 24, where the regions with temperatures above T_m and below T_m are colored in red and blue, respectively, and voids are represented by white areas between the inner and outer surfaces. Clearly, there is no melting when voids began to nucleate (14.5 ps). When the number of voids reaches its first peak (16.5 ps), some red areas appear in the system, indicative of the formation of high-temperature zones (outlined by the yellow dashed box). This can be attributed to the intense increase in temperature accompanying the initiation and propagation of shear localization, as displayed in Figs. 8 and 12. Indeed, this phenomenon was experimentally observed by Shirinkina et al.,12 who confirmed the existence of localized deformation zones with elevated temperature exceeding the melting point of V95 alloy during the collapse of hollow cylindrical shells. Subsequently, these specific zones extend forward to the secondary spall plane. High temperatures promote material softening, facilitating void nucleation. Ultimately, penetration between the two spall planes is produced in the material.

NEMD simulations have been performed to explore the dynamic response of NiTi alloy under implosion loading and its grain-size dependence. For larger grain size, we have found that voids nucleate in a coexisting intergranular/transgranular manner. Moreover, the void evolution exhibits bimodal behavior, indicative of sequential void nucleation. Interestingly, secondary nucleated voids are aggregated in lath-shaped zones concentrated by high shear strain at a certain angle to the radial direction. The void distribution and density profiles reveal that multiple spallation occurs in the sample. According to the radial stress profiles, some spall planes are found to sequentially appear as the tensile stress wave propagates. The maximum tensile stress for promoting the secondary spallation is less than that for the first spallation, which can be attributed to the presence of high shear localization in the second spall region. The resulting temperature rise is primarily responsible for a decrease in spall strength. Moreover, the results show that higher shear strains are concentrated on specific grain boundaries, and in turn the temperature exceeds the melting point, causing penetration between the two spall planes.

For the first spallation, a decrease in grain size leads to a lower spall strength, which can be attributed to the greater number of void-nucleation sites provided by increasing grain boundary content. Surprisingly, at a grain size of 10 nm, there is a transition of the spall strength behavior from an inverse Hall–Petch to a Hall–Petch relationship. For smaller grain size, voids nucleate in an intergranular manner. As a result, the spall strength is increased. In contrast to the first spallation, for the secondary spallation, the spall strength at a grain size of 50 nm is lower than that at 30 nm, which can be attributed to the fact that the large number of SDBs that are formed lead to localized disorder of atoms within grains, promoting grain refinement. The refined grains inhibit the movement of ⟨100⟩ dislocations, leading to a reduction in the dislocation-induced strengthening effect. By contrast, no SDBs are observed in samples with a smaller grain size, and new voids nucleate only at grain boundaries. These results reveal a correlation of the development of shear localization with initial grain size.

Thus, this study has demonstrated the existence of an interaction between multiple spallation and shear localization in NiTi alloys under implosion loading. In addition, the results of the MD modeling may provide useful guidance for experimental work on such systems. As the dynamic response of a material is significantly affected by its microstructure (e.g., grain boundary patterns and initial defects) and by the loading conditions applied, a more precise calibration of the model to enable it to explain the influence of these factors on dynamic failure will be conducted in our future work.

Acknowledgment. The authors gratefully acknowledge the financial support of the National Natural Science Foundation of China under Grant Nos. 12372367 and 12202081, and the Special Foundation from the Institute of Fluid Physics of CAEP under Grant No. 2022-YCHT-0641.