# Numerical Method

## Treatment of the advection equation

System (14) has a volume fraction transport equation which is not conservative. As indicated by Johnsen and Colonius (2006), directly using the non-conservative formula will present difficulties when dealing with the material interface due to an inconsistency between the wave speeds (shock wave, rarefaction and contact discontinuity) and the velocity vector *V. *They advise transforming the advection equation into a conservative formulation to overcome this obstacle for one-dimensional multi-fluid problems. Based on their work, here, we move forward to re-construct the volume fraction equation appropriately for three-dimensional multiphase flow problems.

Introducing a vector / defined as

then the divergence of / is given by Consequently, we can obtain

Substituting Equation (17) into Equation (12) to obtain the quasi-conservative form

then replacing the non-conservative equation in system (14) by Equation (18). the integral formulation of the overall system can now be expressed as

in which О represents the flow domain, *()Q* is its boundary; the vectors U. F and G are given by

where *(n _{x},n_{y},n_{z})* is a unit normal vector across the boundary

*Oil*,

*q*=

*un*+

_{x}+ vn_{y}*wn*and N

_{z }_{g}= [0,0,0,п

_{ж},п

_{у},п

_{2},(/]

^{т}.

## Spatial discretisation

For three-dimensional problems, hexahedrons are used to fill the domain Q. Without loss of generality, if 17 is a cuboid, we may simply partition it using a Cartesian mesh with / x *J* x *К *cells. Within the mesh, any cell can be indicated by the subscripts *i, j* and *k.* As shown in Figure 2, a hexahedral volume cell is closed by six quadrilateral faces. For this cell, six neighbouring cells are selected to form a computational stencil.

Supposing the mesh does not vary with time, the discrete form of Equation (19) for a mesh cell mcan be written as

where the subscript *m* stands for the mesh cell itself, the subscript / represents a mesh cell face and *A* is the area of the face. Introducing the residual defined as

Equation (21) can be written as

As flow discontinuities (shock waves and material interfaces, etc.) may be expected to occur in multi-phase multi-component compressible flows, attention should be paid to the computation of surface flux terms. In the present work, the surface flux term F across any mesh cell face / is evaluated by an approximate Riemann solution based on the numerical flux function

where the symbols ^{+} and ~ indicate the left and right sides of the face considering a normal vector 7*if*. The conservative variables at neighbouring cell centres may be assumed piecewise

Figure 2: A mesh cell and its computational stencil.

constant and assigned directly to UJ and Uj respectively, and this very simple but diffusive treatment is first order accurate in space. To improve the accuracy, the solution data is formally reconstructed using a third order MUSCL scheme (monotone upstream-centred schemes for conservation laws Brain van Leer, 1979) based on the primitive variables W = («1,*p, p**2**,u,v,w,*p)^{T}, written as

where the subscripts L and R represent the left and right neighbouring mesh cells respectively,

in which the subscript LL indicates the left neighbour of the left cell, and RR represents the right neighbour of the right cell. The parameter *к* provides different options of upwind or centred schemes

In the present work, *к =* 1/3 is adopted. To prohibit spurious oscillations introduced by the high order interpolation near discontinuities, a slope limiter function *ф* is utilised to modify the reconstruction procedure as follows

In this study, we apply van Albada’s limiter (van Albada, 1997) defined as

where

In smooth regions, the slope limiter doesn’t operate on spatial derivatives and the MUSCL scheme provides high resolution. The slope limiter operates in critical regions of flow discontinuities to prevent extrema and ensure stability and monotonicity of the solution (van Albada, 1997). The reconstructed conservative variables at the left and right sides of a mesh cell face can be easily obtained as

The numerical flux term represented by Equation (24) is calculated by employing the HLLC approximate Riemann solver (ARS) for the homogeneous mixture.

## The HLLC Riemann solver

The HLLC ARS originates from the work of Harten et al. (1983) and Toro et al. (1992). Its applications for separated multiphase multicomponent compressible flows can be found in the work of Hu et al. (2009) and Johnsen and Colonius (2006), etc. However, to the authors’ knowledge, no work has been reported thus far that solves dispersed multi-phase multi-component compressible flows governed by the one-dimensional Kapila model or three- dimensional system (14) with the HLLC ARS. The numerical flux term represented by Equation (24) is calculated by the HLLC ARS defined as

in which the middle left and right states are evaluated by

the intermediate wave speed can be calculated by

and the intermediate pressure may be estimated as The left and right state wave speeds are computed by

where *q* is an averaged velocity component evaluated by the component Roe-averages (Batten et ah. 1997). The averaged speed of sound c is calculated by the formulation proposed by Hu et al. (2009)

More detail of Equation (37) can be found in reference (Hu et ah, 2009).

## Temporal discretisation

In the present work, a third order total variation diminishing Runge-Kutta scheme (Gottlieb and Shu, 1998) is applied to update the numerical solution from time level *n* to *n* + 1

For a mesh cell *m,* the partial differential operator L is defined as