# A Heterogeneous Parallel Algorithm for Euler-Lagrange Simulations of Liquid in Supersonic Flow
## 2 Algorithm and Code Structures
### Gas-Phase Governing Equations
\(\Delta\) is the forward difference in time, and the right-hand side residual can be written as  
\[Res_{i}^{m}=-\sum\limits_{j=1}^{N}{(\widetilde{\mathbf{F}}_{\mathrm{c},ij}^{m }-\widetilde{\mathbf{F}}_{\mathrm{v},ij}^{m})}+Q_{i}^{m}. \tag{8}\]  
Note that  
\[ \begin{array}{l}\Delta\widetilde{\mathbf{F}}_{\mathrm{c},ij}^{m}-\Delta \widetilde{\mathbf{F}}_{\mathrm{v},ij}^{m}=[\widetilde{\mathbf{F}}_{\mathrm{ c}}(\mathbf{W}_{i}^{m+1},\mathbf{W}_{j}^{m+1})-\widetilde{\mathbf{F}}_{\mathrm{c}}( \mathbf{W}_{i}^{m},\mathbf{W}_{j}^{m+1})]\\ +[\widetilde{\mathbf{F}}_{\mathrm{c}}(\mathbf{W}_{i}^{m},\mathbf{W}_{j}^{m+1})- \widetilde{\mathbf{F}}_{\mathrm{c}}(\mathbf{W}_{i}^{m},\mathbf{W}_{j}^{m})]-[ \widetilde{\mathbf{F}}_{\mathrm{v}}(\mathbf{W}_{i}^{m+1},\mathbf{W}_{j}^{m+1})- \mathbf{F}_{\mathrm{v}}(\mathbf{W}_{i}^{m},\mathbf{W}_{j}^{m+1})]\\ -[\widetilde{\mathbf{F}}_{\mathrm{v}}(\mathbf{W}_{i}^{m},\mathbf{W}_{j}^{m+1})- \widetilde{\mathbf{F}}_{\mathrm{v}}(\mathbf{W}_{i}^{m},\mathbf{W}_{j}^{m})].\end{array}  \tag{9}\]  
[PERSON]'s expansion is performed on the first and the third terms on the right-hand side of Equation (9). \[\widetilde{\mathbf{F}}_{\mathrm{c}}(\mathbf{W}_{i}^{m+1},\mathbf{W}_{j}^{m+1})- \widetilde{\mathbf{F}}_{\mathrm{c}}(\mathbf{W}_{i}^{m},\mathbf{W}_{j}^{m+1})\approx \frac{\partial\widetilde{\mathbf{F}}_{\mathrm{c},ij}}{\partial\mathbf{W}_{i}} \Delta\mathbf{W}_{i}^{m}, \tag{10}\]  
\[\widetilde{\mathbf{F}}_{\mathrm{v}}(\mathbf{W}_{i}^{m+1},\mathbf{W}_{j}^{m+1})- \widetilde{\mathbf{F}}_{\mathrm{v}}(\mathbf{W}_{i}^{m},\mathbf{W}_{j}^{m+1})\approx \frac{\partial\widetilde{\mathbf{F}}_{\mathrm{v},ij}}{\partial\mathbf{W}_{i}} \Delta\mathbf{W}_{i}^{m}. \tag{11}\]  
Substituting Equations (9)-(11) back to Equation (7), gives  
\[ \begin{array}{l}(\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}^{m}+\mathbf{W}_{ i}^{n-1}}{2\Delta t}+\sum\limits_{j=1}^{N}{(\frac{\partial\widetilde{\mathbf{F}}_{ \mathrm{c},ij}}{\partial\mathbf{W}_{i}}-\frac{\partial\widetilde{\mathbf{F}}_{ \mathrm{v},ij}}{\partial\mathbf{W}_{i}})S_{ij}}\\ +\sum\limits_{j=1}^{N}{[\widetilde{\mathbf{F}}_{\mathrm{c}}(\mathbf{W}_{i}^{m}, \mathbf{W}_{j}^{m}+\Delta\mathbf{W}_{j}^{m})-\widetilde{\mathbf{F}}_{\mathrm{c}}(\mathbf{W }_{i}^{m},\mathbf{W}_{j}^{m})]S_{ij}}\\ -\sum\limits_{j=1}^{N}{[\widetilde{\mathbf{F}}_{\mathrm{v}}(\mathbf{W}_{i}^{m}, \mathbf{W}_{j}^{m}+\Delta\mathbf{W}_{j}^{m})-\widetilde{\mathbf{F}}_{\mathrm{v}}(\bm {W}_{i}^{m},\mathbf{W}_{j}^{m})]S_{ij}}=Res_{i}^{m}.\end{array}  \tag{12}\]  
The numerical flux vectors are discretized as  
\[\widetilde{\mathbf{F}}_{\mathrm{c},ij}-\widetilde{\mathbf{F}}_{\mathrm{v},ij} \approx\frac{1}{2}[\mathbf{F}_{\mathrm{c},i}+\mathbf{F}_{\mathrm{c},j}-\lambda _{ij}(\mathbf{W}_{j}-\mathbf{W}_{i})], \tag{13}\]where \(\lambda_{ij}\) is the spectral radius of the flux Jacobian matrix, written as follows:  
\[\lambda_{ij}=\big{|}\mathbf{V}\cdot\mathbf{n}_{ij}\big{|}+a+(\frac{\mu_{ij}}{pr}+\frac{ \mu_{t,ij}}{pr_{t}})\frac{r_{ij}}{\rho_{ij}\big{\|}\mathbf{n}_{ij}\cdot(\mathbf{r}_{j}- \mathbf{r}_{i})\big{\|}}, \tag{14}\]  
\(\mathbf{n}_{ij}\) is the normal vector of the shared face, \(\mathbf{r}_{i}\), \(\mathbf{r}_{j}\) represent the position vectors of cell \(i\) and cell \(j\), \(a\) is the speed of sound, \(\mu_{ij}\), \(\mu_{t,ij}\) are the average of kinematic and turbulent viscosities of the cell \(i\) and cell \(j\). For closed-cell \(i\)  
\[\sum\limits_{j=1}^{N}\frac{\partial\mathbf{F}_{i}(n_{ij})}{\partial W_{i}}S_{ij}=0. \tag{15}\]  
Therefore,  
\[ \begin{array}{l}(\frac{3\Omega_{1}}{2\Delta t}+\frac{\Omega_{2}}{\Delta\tau }+\frac{1}{2}\sum\limits_{j=1}^{N}\lambda_{ij}S_{ij})\Delta\mathbf{W}_{i}^{m}=Re \varepsilon_{i}^{m}-\Omega_{i}\frac{3\mathbf{W}_{i}^{m}-4\mathbf{W}_{i}^{m}+\mathbf{W}_{i }^{m-1}}{2\Delta t}\\ -\frac{1}{2}\sum\limits_{j=1}^{N}[\mathbf{F}(\mathbf{W}_{j}^{m}+\Delta\mathbf{W}_{j}^{m}) -\mathbf{F}(\mathbf{W}_{j}^{m})-\lambda_{ij}\Delta\mathbf{W}_{j}^{m}]S_{ij}.\end{array}