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 (nx,ny,nz) is a unit normal vector across the boundary Oil , q = unx + vny + wnz and Ng = [0,0,0,пж,пу,п2,(/]т.
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 7if. 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, p2,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
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).
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