# An All-Mach Number HLLC-Based Scheme for Multi-Phase Flow with Surface Tension
## 3 Numerical Methodology
### Inviscid HLLC-CICSAM Edge Flux
#### 3.5.4 Surface Tension Term  
HLLC solvers rarely take into account surface tension effects for multi-phase compressible flow. The first mention of an HLLC-based scheme with surface tension was that of [PERSON] et al. [1], where a compression technique was used to remove any numerical diffusion present in the VoF field. In contrast, in this work, the conservative VoF-CICSAM method is combined with the HLLC solver to allow for curvature to be computed accurately. Similar to [1], this requires a well-balanced implementation with the surface tension source term. A well-balanced implementation implies that the surface tension source term is balanced by the computed pressure gradient in a steady flow field with zero velocity and zero gravity. This may be mathematically expressed by considering a static bubble configuration [4; 12; 50; 51] where the Euler equation in volumetric form with a CSF surface tension discretisation reduces to  
\[-\
abla p+\sigma\kappa\
abla\alpha=0.\]  
The above is discretised for a node \(l\) as  
\[\frac{1}{\mathcal{V}_{l}}\sum_{f\in\partial\Omega_{l}}(p\mathbf{n}\mathcal{A})|_{ f}+\varepsilon_{p}=\frac{\sigma\kappa_{l}}{\mathcal{V}_{l}}\sum_{f\in\partial \Omega_{l}}(\alpha\mathbf{n}\mathcal{A})|_{f}+\varepsilon_{\sigma}. \tag{35}\]  
Here, \(\varepsilon_{p}\) and \(\varepsilon_{\sigma}\) are the discretisation errors associated with pressure and alpha, respectively.  
For the case where \(s^{+}\) and \(s^{-}\) are non-zero, the pressure face as per Equation (31) is given as  
\[p_{f}=\frac{1+\text{sign}(s_{*})}{2}[p_{L,M}+\chi_{L}s_{L}\rho_{L,M}s_{*}]+ \frac{1-\text{sign}(s_{*})}{2}[p_{R,M}+\chi_{R}s_{R}\rho_{R,M}s_{*}]. \tag{36}\]  
with \(s_{L}=\min(-\tilde{c},-c_{L})\) and \(s_{R}=\max(\tilde{c},c_{R})\). The contact wave-speed simplifies to  
\[s_{*}=\frac{p_{L,M}-p_{R,M}-\sigma\kappa_{f}(\alpha_{L,\sigma}-\alpha_{R, \sigma})}{\rho_{R,M}s_{R}-\rho_{L,M}s_{L}}. \tag{37}\]  
In the case of a static bubble, as the pressure field converges to the analytical solution, the consistency pressure condition (Equation (24)) is satisfied and hence the contact wave-speed tends to zero (or machine precision, \(\epsilon\)). Equation (36) may therefore be expressed as  
\[p_{f}=\frac{1+\text{sign}(s_{*})}{2}p_{L,M}+\frac{1-\text{sign}(s_{*})}{2}p_{ R,M}, \tag{38}\]  
where \(s_{*}\to\epsilon\).  
As explained by [PERSON] et al. [51], to recover the discrete equilibrium of Equation (35), it is required for \(\varepsilon_{p}\) to cancel \(\varepsilon_{\sigma}\) to zero. To achieve this, the same operator used for the pressure face must also be employed for the volume fraction field in the surface tension term. Hence, \(\alpha_{f,\sigma}\) is interpolated for the left or right state using MUSCL and up-winded using HLLC as follows:  
\[\alpha_{f,\sigma}=\frac{1+\text{sign}(s_{*})}{2}\alpha_{L,M}+\frac{1-\text{ sign}(s_{*})}{2}\alpha_{R,M}, \tag{39}\]  
Similarly, the contact wave-speed now reads for the general case,  
\[s_{*}=\frac{p_{L,M}-p_{R,M}-\sigma\kappa_{f}(\alpha_{L,M}-\alpha_{R,M})}{\rho_ {R,M}(s_{R}-\mathbf{u}_{R,M}\cdot\mathbf{n})-\rho_{L,M}(s_{L}-\mathbf{u}_{L,M}\cdot\mathbf{n})}, \tag{40}\]  
where \(\kappa_{f}\) is averaged as per [1], i.e., \(\kappa_{f}=\frac{\kappa_{L}+\kappa_{m}}{2}\). This approach for the computation of \(\kappa_{f}\) is also widely used for solving the pressure Poiseuille equation in semi-implicit projection solvers [50; 51]. The approximation is valid for a sufficiently fine mesh and for cases where the curvature is uniform and constant.  
Finally, the energy surface tension term is re-expressed in its conservative form and discretised as follows:  
\[\frac{\sigma\kappa}{\mathcal{V}_{l}}\int_{\Omega}\
abla\alpha\cdot ud\mathcal{ V}\approx\frac{\sigma\kappa_{l}}{\mathcal{V}_{l}}\Bigg{(}\sum_{f\in\partial \Omega_{l}}\alpha_{\sigma,f}u_{ff}-\alpha_{l}\sum_{f\in\partial\Omega_{l}}u_{ ff}\Bigg{)}. \tag{41}\]  
This guarantees a well-balanced method as demonstrated numerically via the static bubble test-case in Section 4.2.4.