# A Heterogeneous Parallel Algorithm for Euler-Lagrange Simulations of Liquid in Supersonic Flow
## 2 Algorithm and Code Structures
### Liquid-Phase Equations  
The Lagrange approach, which is an efficient method that supports parallel computing for full-field simulation is used to model the spray dynamics. In addition, the computational droplet is used in this paper. One computational droplet represents a collection of real droplets with the same properties (diameter, velocity, temperature, etc.) and in the same location [40]. The number of real droplets carried by each computational droplet is represented by a parameter \(\omega_{k}\). Since the density of the liquid phase is much larger than that of the gas phase, it is possible to ignore the unsteady drag forces in the Basset-Boussinesq-Oseen equation (BBO equation), and only the Stokes drag force is included in the liquid momentum equation. Because of the small liquid volume fraction, the gas-liquid flow can be treated as dilute two-phase flow. As a result, the collision and coalescence among droplets are also ignored [13].  
[PERSON]'s second law determines the motion of each computational droplet as follows:  
\[\frac{dx_{k,i}}{dt}=u_{k,i}, \tag{20}\]  
\[\frac{du_{k,i}}{dt}=\frac{F_{k,i}}{m_{k}}, \tag{21}\]  
where \(x_{k,i}\), \(m_{k}\), \(u_{k,i}\) are the position vector, mass, and velocity of the droplet, respectively, \(i=x\), \(y\), \(z\). The surface drag \(F_{k,i}\), which is a force acting on the droplet body, can be modeled by:  
\[\frac{F_{k,i}}{m_{k}}=\frac{3}{4}\frac{\rho_{\rm g}C_{k}}{\rho_{\rm d}d_{k}}|u_ {\rm g,i}-u_{k,i}|\big{(}u_{\rm g,i}-u_{k,i}\big{)}, \tag{22}\]  
where \(d_{k}\), \(\rho_{\rm d}\), \(u_{k}\) represent the droplet diameter, density, and velocity, the \(\rho_{\rm g}\), \(u_{\rm g}\) are the gas-phase density and velocity, which are calculated by weighting the inverse distance between the control volume containing the droplet and its adjacent control volumes. \(C_{k}\) is the drag coefficient, which can be calculated by \(C_{k}=f_{kf}C_{k}(M_{k},{\rm Re}_{k})\). \(M_{k}\) and \(f_{kf}\) are the relative Mach number and the droplet deformation correction factor, respectively. The calculation of the modified drag coefficient \(C_{k\rm s}\) is detailed in reference [22].  
The effects of the droplet on the gas phase can be calculated as follows:  
\[\dot{Q}_{\rm m}=-\frac{\omega_{k}}{\Omega}\sum_{k}\frac{dm_{k}}{dt}, \tag{23}\]  
\[\dot{Q}_{\rm p,x}=-\frac{\omega_{k}}{\Omega}\sum_{k}\frac{dm_{k}u_{k}}{dt}, \tag{24}\]  
\[\dot{Q}_{\rm p,y}=-\frac{\omega_{k}}{\Omega}\sum_{k}\frac{dm_{k}v_{k}}{dt}, \tag{25}\]  
\[\dot{Q}_{\rm p,z}=-\frac{\omega_{k}}{\Omega}\sum_{k}\frac{dm_{k}w_{k}}{dt}, \tag{26}\]  
\[\dot{Q}_{\rm E}=-\frac{\omega_{k}}{\Omega}\sum_{k}\frac{d(m_{k}C_{pl}T_{k}+m_{ k}\mathbf{U}_{k,i}^{2})}{dt}, \tag{27}\]  
where, \(m_{k}\), \(T_{k}\), are the mass and temperature of the \(k\)th droplet respectively, \(u_{k}\), \(v_{k}\), \(w_{k}\) are the velocities in \(x\), \(y\), and \(z\) directions of the \(k\)th droplet, \(C_{pl}\) is the specific heat of the liquid phase. \(\Omega\) is the volume of the cell containing the droplet. In this paper, evaporation is not considered, therefore,  
\[\dot{Q}_{\rm m}=0, \tag{28}\]  
\[\dot{Q}_{E}=-\frac{\omega_{k}}{\Omega}\sum_{k}\frac{d(m_{k}\mathbf{U}_{k,i}^{2})} {dt}. \tag{29}\]  
In the present study, the droplet stripping process is modeled by Kelvin-Helmholtz model (KH model). Then, the droplet deformation and oscillation are modeled by the Taylor Analogy Breakup model (TAB model), and the TAB model and Rayleigh-Taylor model (RT model) compete to simulate the droplet secondary breakup process. It is worth noting that both the TAB and RT models are explosion-type breakup modes, where the original droplets are treated as if they disappeared while smaller droplets are newly produced. More information can be found in reference [22]. The three-step Runge-Kutta method with third-order accuracy is used for the temporal discretization of droplet iteration.