# A Heterogeneous Parallel Algorithm for Euler-Lagrange Simulations of Liquid in Supersonic Flow
## 2 Algorithm and Code Structures
### Gas-Phase Governing Equations  
The governing equations based on the integral form for the gas-phase can be written as follows [31; 37; 38]:  
\[\frac{\partial}{\partial t}\int\limits_{\Omega}\textbf{W}d\Omega+\oint \limits_{\partial\Omega}\left(\textbf{F}_{\text{c}}-\textbf{F}_{\text{v}} \right)dS=\int\limits_{\Omega}\textbf{Q}d\Omega, \tag{1}\]  
\[\textbf{W}= \begin{pmatrix}\rho\\ \rho u\\ \rho v\\ \rho v\\ \rho E\end{pmatrix} ,\textbf{F}_{\text{c}}= \begin{pmatrix}\rho V\\ \rho v+n_{x}p\\ \rho v+n_{y}p\\ \rho wV+n_{z}p\\ V(\rho E+p)\end{pmatrix} ,\textbf{F}_{\text{v}}= \begin{pmatrix}0\\ n_{x}\tau_{xx}+n_{y}\tau_{xy}+n_{z}\tau_{xz}\\ n_{x}\tau_{yx}+n_{y}\tau_{yy}+n_{z}\tau_{yz}\\ n_{x}\tau_{zx}+n_{y}\tau_{zy}+n_{z}\tau_{zz}\\ n_{x}\Theta_{x}+n_{y}\Theta_{y}+n_{z}\Theta_{z}\end{pmatrix} ,\textbf{Q}= \begin{pmatrix}\dot{Q}_{\text{m}}\\ \dot{Q}_{\text{p},x}\\ \dot{Q}_{\text{p},y}\\ \dot{Q}_{\text{p},z}\\ \dot{Q}_{\text{E}}\end{pmatrix} , \tag{2}\]  
where  
\[V=n_{x}u+n_{y}v+n_{z}w, \tag{3}\]  
\[\Theta_{x}=u\tau_{xx}+v\tau_{xy}+w\tau_{xz}+k\frac{\partial T}{ \partial x},\] \[\Theta_{y}=u\tau_{yx}+v\tau_{yy}+w\tau_{yz}+k\frac{\partial T}{ \partial y}, \tag{4}\] \[\Theta_{z}=u\tau_{zx}+v\tau_{zy}+w\tau_{zz}+k\frac{\partial T}{ \partial x}.\]  
Here \(\textbf{W},\textbf{F}_{\text{c}},\textbf{F}_{\text{v}}\) are the vector of conservative variables, convective fluxes, and viscous fluxes, respectively, \(\Omega\) is a control volume, and \(V\) is defined as the scalar product of the velocity vector and the unit normal vector. \(\dot{Q}_{\text{m}}\) is the mass flow rate, \(\dot{Q}_{\text{p},x^{\prime}}\)\(\dot{Q}_{\text{p},y^{\prime}}\)\(\dot{Q}_{\text{p},z}\) are the momentum in the \(x\), \(y\), and \(z\) directions, respectively, and \(\dot{Q}_{\text{E}}\) is the energy source terms of the droplet donation, and the functions will be introduced in Section 2.2. For the vector of viscous fluxes, the components of viscous stress \(\tau_{ij}\) are as follows [37]:  
\[\tau_{ij}=\mu(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{ \partial x_{i}}-\frac{2}{3}\frac{\partial u_{k}}{\partial x_{k}}\delta_{ij}). \tag{5}\]  
The spatial discretization of Equation (1) is based on the cell-centered finite volume method. The upwind Advection upstream splitting method + UP (AUSM + UP) and central scheme are used to solve the convective fluxes and the viscous fluxes, respectively. A Monotone upstream-centered scheme for conservation laws (MUSCL) reconstruction scheme is used, and the van Leer limiter is employed to make the scheme monotone. The \(k-\omega\) shear stress transport (SST) turbulence model is adopted for the turbulent model. A liquid in supersonic flow is a typical unsteady compressible flow problem. A dual time-scheme based on Data Parallel Lower-Upper Relaxation Method (DP-LUR) is used as follows. Equation (1) is discretized spatially in the \(i\)th cell to yield. \[\Omega_{i}\frac{\partial\mathbf{W}_{i}}{\partial t}+\sum\limits_{j=1}^{N}{( \widetilde{\mathbf{F}}_{\mathrm{c},ij}-\widetilde{\mathbf{F}}_{\mathrm{v},ij} )S_{ij}}=Q_{i}, \tag{6}\]  
where \(N\) is the total number of neighboring cells that share a face with cell \(i\), \(S_{ij}\) is the area of the face shared by cell \(i\) and cell \(j\). \(\widetilde{\mathbf{F}}_{\mathrm{c}}\) and \(\widetilde{\mathbf{F}}_{\mathrm{v}}\) are the second-order-accurate numerical inviscid and viscous flux vectors at the face, respectively. The time derivative is discretized with a 3-point backward-difference dispersion in Equation (6), and a pseudo time derivative term is added to the left, resulting in  
\[(\frac{3\Omega_{i}}{2\Delta t}+\frac{\Omega_{i}}{\Delta\tau})\Delta\mathbf{W}_{i}^ {m}+\Omega_{i}\frac{3\mathbf{W}_{i}^{m}-4\mathbf{W}_{i}^{n}+\mathbf{W}_{i}^{n-1}}{2\Delta t }+\sum\limits_{j=1}^{N}{(\Delta\widetilde{\mathbf{F}}_{\mathrm{c},ij}}^{m}- \Delta\widetilde{\mathbf{F}}_{\mathrm{v},ij}^{m})S_{ij}}=Res_{i}^{m}, \tag{7}\]  
where the superscript \(m\) represents the \(m\)th pseudo time.