## 3 Numerical Simulation of Shockwave Soil Loosening
### Governing Equations and the Preconditioning System  
This research employs a multi-block grid technique, utilizing a CFD algorithm that integrates the control volume method and the preconditioning approach to address Navier-Stokes equations for modeling interactions between compressible and incompressible flows. The inviscid flux vector \(\overset{\rightarrow}{F}\) is expressed in the standard conservation form as follows:  
\[\overset{\rightarrow}{F}= \begin{bmatrix}\rho\overset{\rightarrow}{V}\\ \rho\overset{\rightarrow}{V}v_{x}+p\overset{\rightarrow}{i}\\ \rho\overset{\rightarrow}{V}v_{y}+p\overset{\rightarrow}{i}\\ \rho\overset{\rightarrow}{V}v_{z}+p\overset{\rightarrow}{k}\\ \rho\overset{\rightarrow}{V}E+p\overset{\rightarrow}{V}\end{bmatrix}  \tag{1}\]  
Figure 4: Design and production of a shock output tube: (**a**) design of a shock output tube; (**b**) physical entity of a shock output tube. In this context, \(\rho\) represents the density, while \(\overset{\rightarrow}{V}\) denotes the velocity vector, with \(v_{x}\), \(v_{y}\), and \(v_{z}\) corresponding to the velocity components in the \(x\), \(y\), and \(z\) directions, respectively. The parameter \(p\) signifies pressure, and \(E\) denotes the total energy per unit mass. The relationship between \(\mathbf{E}\) and the total enthalpy \(H\) is given by \(H\) is \(H=E-p/\rho\), where \(H=C_{p}T+V^{2}/2\), with \(C_{p}\) being the specific heat at a constant pressure and \(T\) representing the temperature. The system is governed by the equation of state, typically expressed as \(\rho=\rho(p,T)\). In this study, the AUSM+ method, which employs a low-diffusion flux splitting technique, is applied [27]. The inviscid interface flux \(F_{i+1/2}\) in the \(x\) direction can be separated into two distinct components: the convective contribution \(F_{1/2}^{c}\) and the pressure contribution \(F_{1/2}^{p}\). Here, the subscript \(1/2\) denotes the midpoint of the computational grid, and both the convection term \(F_{1/2}^{c}\) and the pressure term \(F_{1/2}^{p}\) are formulated as follows:  
\[F_{1/2}^{c}=\left(\rho v_{x}\right)_{1/2}\left[ \begin{array}{c}1\\ v_{x}\\ v_{y}\\ v_{z}\\ H\end{array} \right]_{i/i+1},\;\;F_{1/2}^{p}=\left(\rho v_{x}\right)_{1/2} \left[ \begin{array}{c}0\\ p_{1/2}\\ 0\\ 0\end{array} \right] \tag{2}\]  
In the convection term \(F_{1/2}^{c}\), if interface mass flux \((\rho v_{x})_{1/2}\) is non-negative, then state \(i\) is selected for column vector (\(1\), \(v_{x}\), \(v_{y}\), \(v_{z}\), \(H\))\({}^{T}\); if \((\rho v_{x})_{1/2}\) is negative, then state \(i+1\) is selected. The interface quantities \(\left(\rho v_{x}\right)_{1/2}\) and \(p_{1/2}\) are defined as follows:  
\[\left(\rho v_{x}\right)_{1/2}=a_{1/2}\left(\rho_{i}m_{1/2}^{+}+\rho_{i+1}m_{1/ 2}^{-}\right) \tag{3}\]  
\[p_{1/2}=P_{(5)}^{+}(M_{i})p_{i}+P_{(5)}^{-}(M_{i+1})p_{i+1} \tag{4}\]  
Among them, \(a_{1/2}\) is the interface sound speed, \(M_{i}\) is the MACH number, \(M_{i}=v_{xi}/a_{1/2}\), and \(m_{1/2}^{+}\) and \(m_{1/2}^{-}\) are defined as follows:  
\[m_{1/2}^{+}=M_{(4)}^{+}(M_{i})+M_{(4)}^{-}(M_{i+1}) \tag{5}\]  
\[m_{1/2}^{-}=\frac{1}{2}(M_{1/2}+|M_{1/2}|) \tag{6}\]  
The terms \(P_{(5)}^{+}\) and \(P_{(5)}^{-}\) in Equation (4) are both fifth-degree polynomials of \(M\), and they are defined as follows:  
\[P_{(5)}^{\pm}\left\{ \begin{array}{c}\frac{(M\pm 1)^{2}(2\mp M)}{4}\pm \frac{3M\left(M^{2}-1\right)^{2}}{16},\;if\;|M|<1\\ \frac{M+|M|}{2M},\;otherwise\end{array} \right. \tag{7}\]  
Similarly, \(M_{(4)}^{+}\) and \(M_{(4)}^{-}\) in Equation (5) represent fourth-degree polynomials of \(M\), with their definitions given below:  
\[M_{(4)}^{\pm}\left\{ \begin{array}{c}\frac{\pm(M\pm 1)^{2}}{4},\;if\;|M|<1\\ \frac{M\pm|M|}{2},\;otherwise\end{array} \right.