Multiple elastic shock waves in cubic single crystals

Abstract

Multiple elastic shock waves carry the information on elastic properties under dynamic extreme conditions, but may complicate the interpretation of wave structure including the elastic–plastic transition. On the basis of the acoustic wave-equation analysis, we predict the absence or presence of multiple elastic shock waves in a single crystal subjected to shock loading along a specific crystallographic orientation. Typical FCC and BCC single crystals are taken as validation and application cases. Large-scale molecular dynamics simulations are performed for Cu and Ta; double-wave or triple-wave structures of elastic shock waves (quasilongitudinal and quasitransverse) are observed in the simulations, and the multi-wave structures are in excellent agreement with the wave-equation analysis. Also, the acoustic wave-equation analysis is used to analyze MD calculations, as well as the complex structure of the shock wave during plastic deformation. Free-surface velocity history, transverse velocity history of free surface, and ultrafast X-ray diffraction are explored as experimental means to resolve multiple elastic shock waves.

Appendices

Appendix 1: More details of wave-equation analysis

Substituting (2) into (1) yields

$$\begin{aligned} \begin{aligned} (c_{ijkl}n_jn_l-\delta _{ik}\rho v_{\textrm{c}}^2)S_k&= 0, \end{aligned} \end{aligned}$$

(7)

where $\delta _{ik}$ is the Kronecker delta function. For a plane wave propagating along the x-axis, ${\textbf {n}} = (1,0,0)$; $c_{ijkl}$ at a new coordinate system can be obtained as

$$\begin{aligned} \begin{aligned} {c_{ijkl} = Q_{im}Q_{jn}Q_{kp}Q_{lq}c_{mnpq},} \end{aligned} \end{aligned}$$

(8)

where Q is the rotation matrix from new coordinate system to the original coordinate system, $c_{mnpq}$ is the elastic stiffness tensor at the original coordinate system.

Equation (7) can be simplified to [43]

$$\begin{aligned} \begin{aligned} \begin{bmatrix}C_{11}\!-\!\rho v_{\textrm{c}}^2&{}\quad C_{16}&{}\quad C_{15}\\ C_{16}&{}\quad C_{66}-\rho v_{\textrm{c}}^2&{}\quad C_{56}\\ C_{15}&{}\quad C_{56}&{}\quad C_{55}-\rho v_{\textrm{c}}^2\end{bmatrix}\begin{bmatrix}a_x\\ a_y\\ a_z \end{bmatrix}&\!=\! \begin{bmatrix}0\\ 0\\ 0\end{bmatrix}, \end{aligned}\nonumber \\ \end{aligned}$$

(9)

where ${\textbf {a}}=(a_x,\,a_y,\,a_z)$ is the unit direction vector for ${\textbf {s}}$.

To obtain the nonzero solutions for (9), the determinant of the matrix should equal zero, i.e.,

$$\begin{aligned} \left| \begin{array}{ccc} C_{11}-\rho v_{\textrm{c}}^2 &{}\quad C_{16} &{}\quad C_{15} \\ C_{16} &{}\quad C_{66}-\rho v_{\textrm{c}}^2 &{}\quad C_{56} \\ C_{15} &{}\quad C_{56} &{}\quad C_{55}-\rho v_{\textrm{c}}^2 \end{array}\right| = 0. \end{aligned}$$

(10)

According to the Cardano’s formula, three solutions for $v_{\textrm{c}}$, $v_{\mathrm{c_1}}$, $v_{\mathrm{c_2}}$, and $v_{\mathrm{c_3}}$ can be obtained from (10), and here $v_{\mathrm{c_1}}>v_{\mathrm{c_2}}\ge v_{\mathrm{c_3}}$. According to (9), the unit direction vector, ${\textbf {a}}_1$, ${\textbf {a}}_2$, and ${\textbf {a}}_3$ can be obtained as (3).

For $C_{15}=0$ and $C_{56}=0$, the Cardano’s formula is not applicable, and (10) is written as

$$\begin{aligned}{} & {} (C_{55}-\rho v_{\textrm{c}}^2)[(\rho v_{\textrm{c}}^2)^2-(C_{11}+C_{66})\rho v_{\textrm{c}}^2\nonumber \\{} & {} \quad +C_{11}C_{66}-(C_{16})^2] = 0. \end{aligned}$$

(11)

Three solutions for $v_{\textrm{c}}$ can be obtained as well, and $v_{\mathrm{c_1}}>v_{\mathrm{c_2}}\ge v_{\mathrm{c_3}}$.

Appendix 2: Rotation matrix

Given Euler angles $(\varphi _1,\Phi ,\varphi _2)$ [67, 68], the rotation matrix ${\textbf {R}}$ can be written as

$$\begin{aligned} \begin{aligned}&{\textbf {R}} =\begin{bmatrix}\cos {\varphi _1}\cos {\varphi _2}-\cos {\Phi }\sin {\varphi _1}\sin {\varphi _2} &{}\quad -\sin {\varphi _1}\cos {\Phi }\cos {\varphi _2}-\cos {\varphi _1}\sin {\varphi _2} &{}\quad \sin {\varphi _1}\sin {\Phi }\\ \sin {\varphi _1}\cos {\varphi _2}+\cos {\varphi _1} \cos {\Phi }s_3&{}\quad \cos {\varphi _1}\cos {\Phi }\cos {\varphi _2}-\sin {\varphi _1}\sin {\varphi _2} &{}\quad -\cos {\varphi _1}\sin {\Phi }\\ \sin {\Phi }s_3&{}\sin {\Phi }\cos {\varphi _2}&{}\cos {\Phi }\end{bmatrix}. \end{aligned} \end{aligned}$$

(12)

The unit vector ${\textbf {p}}$ $(p_x,p_y,p_z)$ (corresponds to the [100] direction closest relative to the x-axis) can be calculated as

$$\begin{aligned} \begin{aligned} {\textbf {R}}\cdot {\textbf {n}}={\textbf {p}}. \end{aligned} \end{aligned}$$

(13)

Here, ${\textbf {n}}=(1,0,0)$ is the unit vector along the x-axis. As shown in Fig. 1, the angle $\beta$ (between ${\textbf {n}}$ and ${\textbf {p}}$) and $\phi$ (between the projection vector of ${\textbf {p}}$ onto the yz-plane and the y-axis) are

$$\begin{aligned} \begin{aligned} \beta&= \arccos {({\textbf {p}}\cdot {\textbf {n}})}, \end{aligned} \end{aligned}$$

(14)

and

$$\begin{aligned} \begin{aligned} \phi&= \arctan {\frac{p_z}{p_y}}. \end{aligned} \end{aligned}$$

(15)

Appendix 3: Supplementary results of MD simulation for single-crystal Cu at $u_{\textrm{p}}=0.1$ km s$^{-1}$

Upon the arrival of a QL wave, shear stress $\sigma _{xy}$ is nonzero (Fig. 7e) and drives the transverse movement of atoms and shear strain $\gamma _{xy}$ in single-crystal Cu. The profiles of particle velocity component $u_y$ and transverse displacement $s_y$ for [311] Cu at 45 ps are shown in Fig. 18a. In the $\mathrm F_2F_1$ segment, $s_y$ decreases with increasing x. Particles move transversely at $u_y=16.88$ $\mathrm m\,s^{-1}$. In the $\mathrm OF_2$ segment, the transverse motion stops; $s_y$ increases from 0 to peak P. Shear strain is calculated as $\gamma _{xy} = \frac{\partial s_y}{\partial x}+\frac{\partial s_x}{\partial y}$. In shocked Cu, $\frac{\partial s_x}{\partial y}=0$, so $\gamma _{xy}=\frac{\partial s_y}{\partial x}$; $\gamma _{xy}^{(1)}=-3.62\times 10^{-3}$ and $\gamma _{xy}^{(2)}=4.56\times 10^{-3}$ for S1 and S2, respectively. Different from the [311] case, $u_z$ (Fig. 7c) and shear stress $\sigma _{xz}$ (Fig. 7f) for [321] Cu are not zero and drive transverse movement of atoms in the xz-plane and giving rise to shear strain $\gamma _{xz}$. Figure 18b shows the profiles of transverse displacement $s_y$, $s_z$, and particle velocity component $u_z$. Shear strain $\gamma _{xy} = \frac{\partial s_y}{\partial x}+\frac{\partial s_x}{\partial y}$ and $\gamma _{xz} = \frac{\partial s_z}{\partial x}+\frac{\partial s_x}{\partial z}$. However, partial derivative $\frac{\partial s_x}{\partial z}=0$ and $\frac{\partial s_x}{\partial y}=0$, so $\gamma _{xy}=\frac{\partial s_y}{\partial x}$ and $\gamma _{xz}=\frac{\partial s_z}{\partial x}$. In Fig. 18b, the profiles of $s_y(x)$ and $s_z(x)$ indicate that $\gamma _{xy}$ and $\gamma _{xz}$ of S1, S2, and S3 achieve a constant value. For S1, $\gamma _{xy}^{(1)}=-2.06\times 10^{-3}$ and $\gamma _{xz}^{(1)}=0$. For S2, $\gamma _{xy}^{(2)}=-1.58\times 10^{-3}$ and $\gamma _{xz}^{(2)}=-1.37\times 10^{-3}$. For S3 ($\mathrm OF_3$), $\gamma _{xy}^{(3)}=3.32\times 10^{-3}$ and $\gamma _{xz}^{(3)}=0.73\times 10^{-3}$.

Fig. 18

a Profiles of transverse particle displacement $s_y$ and particle velocity component $u_y$ for [311] Cu at 45 ps. b Profiles of transverse particle displacements $s_y$ and $s_z$, and particle velocity component $u_z$ for [321] Cu at 45 ps; O: impact location; $\textrm{F}_1$: wavefront of S1; $\textrm{F}_2$: wavefront of S2; $\textrm{F}_3$: wavefront of S3

Fig. 19

Position–time or $x{-}t$ diagram for flyer plate impact loading at $u_{\textrm{p}}=0.1$ km s$^{-1}$. Color coding is based on particle velocity component $u_x$; U: unshocked region; S1: elastic shock state 1; S2: elastic shock state 2; R1: release wave 1; N: a state induced by the interaction of S2 and R1

Appendix 4: Possible experimental measurements

For [311] single-crystal Cu, the wave propagation and interaction are detailed in the position–time ($x{-}t$) diagram in terms of particle velocity component $u_x$, as shown in Fig. 19. The impact yields double elastic shock waves (correspond to S1 and S2) propagating into the target. The wavefront of S1 first reaches the free surface (D in Fig. 19) and then reflects at the free surface as a centered simple release fan (R1) traveling toward the opposite direction. When R1 meets S2, the interaction of R1 and S2 gives rise to a new state N. The wavefront of N travels faster than that of S2 and reaches the free surface first (E); later, the wavefront of S2 reaches the free surface (F). Thus, a triple-platform structure is observed on the $u_{\textrm{fs}}(t)$ curve (Fig. 14a).

In XRD simulations with GAPD [72, 73], the diffraction intensity $I({\textbf {q}})$ is [76]

$$\begin{aligned} \begin{aligned} I({\textbf {q}}) = \frac{F({\textbf {q}})F^{*}({\textbf {q}})}{G}. \end{aligned} \end{aligned}$$

(16)

Here, ${\textbf {q}}$ is the scattering vector. $F({\textbf {q}})$ is the structure vector, and its complex conjugate $F^{*}({\textbf {q}})$. G is the number of atoms in the selected area. The structure vector $F({\textbf {q}})$ is

$$\begin{aligned} \begin{aligned} F({\textbf {q}}) = \sum _{j=1}^G f_j \textrm{exp}(2\pi i{\textbf {q}}\cdot {\textbf {x}}_j), \end{aligned} \end{aligned}$$

(17)

where f is the atomic scattering factor for X-rays [77] and ${\textbf {x}}_j$ is the position of the jth atom in the real space. According to the Bragg’s law [76], ${\textbf {q}}$ can be obtained with

$$\begin{aligned} \begin{aligned} |{\textbf {q}}|= \frac{2\sin {\theta }}{\lambda } = \frac{1}{d_{hkl}}, \end{aligned} \end{aligned}$$

(18)

where 2$\theta$ is the diffraction angle, $\lambda$ is the wavelength, and $d_{hkl}$ is the interplanar spacing. For a polychromatic beam, the intensity at a specific position on a 2D detector, $I(2\theta ,\gamma )$, is the weighted integration over the incident beam wavelength range, [$\lambda _1, \lambda _2$], i.e.,

$$\begin{aligned} \begin{aligned} I(2\theta ,\gamma )=\frac{\int _{\lambda _1}^{\lambda _2} \ I(2\theta ,\gamma ,\lambda )\,\textrm{d}\lambda }{\int _{\lambda _1}^{\lambda _2} \ w(\lambda )\, \textrm{d}\lambda }. \end{aligned} \end{aligned}$$

(19)

Here, $\gamma$ is the azimuthal angle (Fig. 16) and $w(\lambda )$ is the flux fraction of the incident beam. Each set of $(2\theta ,\gamma ,\lambda )$ corresponds to a scattering vector ${\textbf {q}}$. More details including intensity projection onto a detector were presented elsewhere [72, 73].