## 3 Numerical Model
### Governing Equations  
Three-dimensional compressible Navier-Stokes equations were established based on conservation of mass, momentum, and energy [24; 25; 26]. The CFD software SC/Tetra was used to simulate the gas-solid two-phase flow of solid quartz sand particles driven by the shock wave. The air was described by the Euler method; the solid particles were described by the Lagrange method. The governing equations are written as follows.  
Mass conservation equation:  
\[\frac{\partial\rho}{\partial t}+\frac{\partial\rho u_{i}}{\partial x_{i}}=0 \tag{1}\]  
Momentum conservation equation:  
\[\frac{\partial\rho u_{i}}{\partial t}+\frac{\partial u_{j}\rho u_{i}}{\partial x _{j}}=-\frac{\partial\alpha_{ij}}{\partial x_{j}}+\rho g_{i} \tag{2}\]  
Energy conservation equation:  
\[\frac{\partial\rho H}{\partial t}+\frac{\partial u_{j}\rho H}{\partial x_{j}}= \frac{\partial P}{\partial t}+\frac{\partial u_{j}P}{\partial x_{j}}+\sigma_ {ij}\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial}{\partial x_{j}}K \frac{\partial T}{\partial x_{j}}+Q \tag{3}\]  
Gas state equation:  
\[P=\rho RT \tag{4}\]  
where \(\rho\) is the density of air; \(t\) is the time; \(u_{i}\) is the velocity of the \(x_{i}\) direction; \(\sigma_{ij}\) is the stress tensor; \(\mu\) is viscosity; \(H\) is the specific enthalpy; \(g_{x}\) is the gravity; \(T\) is the temperature; \(K\) is the thermal conductivity; \(Q\) is the heat source; \(P\) is the pressure; and \(R\) is the gas constant.  
In addition, during the momentum exchange between the shock wave and the solid particles, a drag force is generated, which manifests that the particles get an initial velocity and disperse outwards. To accurately calculate the drag force of shockwave on particles, it was defined as a function of the Reynolds number (Re), and the coupling effect between the gas and solid phases considered simultaneously. The equilibrium equations of particle force during the interaction process are as follows:  
\[\frac{du_{p}}{d_{t}}=\frac{u-u_{p}}{t_{\rm s}}+\frac{gx\Big{(}\rho_{p}-\rho \Big{)}}{\rho_{p}}+F_{x} \tag{5}\]  
\[t_{\rm s}=\frac{\rho_{p}d_{p}^{2}}{18\mu}\frac{24}{C_{\rm D}\rm R_{e}} \tag{6}\]  
Figure 3: Jacobian matrix of mesh quality.  
\[\mathrm{R_{e}}=\frac{\rho d_{\rho}\big{|}u_{p}-u\big{|}}{\mu} \tag{7}\] \[C_{D}=a_{1+}\frac{a_{2}}{\mathrm{R_{e}}}+\frac{a_{3}}{\mathrm{R_{e }^{2}}} \tag{8}\]  
where \(u\) is the velocity of air; \(u_{p}\) is the velocity of particle; \(\rho_{p}\) is particle density; \(d_{p}\) is the particle diameter; \(C_{D}\) is the drag coefficient; \(t_{s}\) is the relaxation time; \(F_{x}\) is the additional acceleration (force/unit particle mass) term; and \(a_{1}\), \(a_{2}\), and \(a_{3}\) are constants.  
The \(k\)-\(\omega\) turbulent model was adopted to describe the turbulence phenomena in the numerical simulation. It is an empirical model based on model transport equations for the turbulence kinetic energy \(k\) and the specific dissipation rate \(\omega\)[27].  
The turbulence kinetic energy \(k\) and specific dissipation rate \(\omega\) can be obtained from the following transport equations:  
\[\frac{\partial}{\partial t}\rho k+\frac{\partial}{\partial x_{i}}(\rho ku_{i} )=\frac{\partial}{\partial x_{j}}\bigg{(}\Gamma_{k}\frac{\partial k}{ \partial x_{j}}\bigg{)}+G_{k}-Y_{k}+S_{k} \tag{9}\]  
\[\frac{\partial}{\partial t}\rho\omega+\frac{\partial}{\partial x_{i}}(\rho \omega u_{i})=\frac{\partial}{\partial x_{j}}\bigg{(}\Gamma_{\omega}\frac{ \partial\omega}{\partial x_{j}}\bigg{)}+G_{\omega}-Y_{\omega}+S_{\omega} \tag{10}\]  
where \(G_{k}\) is the generation of turbulence kinetic energy due to mean velocity gradients; \(G_{\omega}\) is the generation of \(\omega\); \(\Gamma_{k}\) and \(\Gamma_{\omega}\) represent the effective diffusivity of \(k\) and \(\omega\); \(Y_{k}\) and \(Y_{\omega}\) are the dissipation of \(k\) and \(\omega\) due to turbulence; and \(S_{k}\) and \(S_{\omega}\) are user-defined source terms.