An Eulerian interface sharpening algorithm for compressible two-phase flow: The algebraic THINC
approach
Keh-Ming Shyuea, Feng Xiaob
aInstitute of Applied Mathematical Sciences, National Taiwan University, Taipei 10617, Taiwan
bDepartment of Energy Sciences, Tokyo Institute of Technology, 4259 Nagatsuta Midori-ku, Yokohama 226-8502, Japan
Abstract
We describe a novel interface-sharpening approach for efficient numerical resolution of a compressible homogeneous two-phase flow governed by a quasi-conservative five-equation model of Allaire et al. (J. Comput. Phys.
181 (2002) 577-616). The algorithm uses a semi-discrete wave propagation method to find approximate solution of this model numerically. In the al- gorithm, in regions near the interfaces where two different fluid components are present within a cell, the THINC (Tangent of Hyperbola for INterface Capturing) scheme is used as a basis for the reconstruction of a sub-grid dis- continuity of volume fractions at each cell edge, and it is complemented by a homogeneous-equilibrium-consistent technique that is derived to ensure a consistent modeling of the other interpolated physical variables in the model.
In regions away from the interfaces where the flow is single phase, standard reconstruction scheme such as MUSCL or WENO can be used for obtaining high-order interpolated states. These reconstructions are then used as the initial data for Riemann problems, and the resulting fluctuations form the basis for the spatial discretization. Time integration of the algorithm is done by employing a strong stability-preserving Runge-Kutta method. Numeri- cal results are shown for sample problems with the Mie-Gr¨uneisen equation of state for characterizing the materials of interests in both one and two space dimensions that demonstrate the feasibility of the proposed method
Email addresses: shyue@ntu.edu.tw(Keh-Ming Shyue), xiao@es.titech.ac.jp (Feng Xiao)
for interface-sharpening of compressible two-phase flow. To demonstrate the competitiveness of our approach, we have also included results obtained using the anti-diffusion interface sharpening method.
Key words: Compressible two-phase flow, Five-equation model, Interface sharpening, THINC reconstruction, Mie-Gr¨uneisen equation of state, Semi-discrete wave propagation method
2000 MSC: 65M06, 65M50, 65Y15, 35L60, 58J45
1. Introduction
Our goal is to describe a novel Eulerian interface-sharpening approach for the efficient numerical resolution of problems with material interfaces arising from inviscid compressible two-phase flow. We consider an unsteady, inviscid, homogeneous two-phase flow that is governed by a five-equation model system of the form
∂
∂t (α1ρ1) + ∇ · (α1ρ1~u) = 0,
∂
∂t (α2ρ2) + ∇ · (α2ρ2~u) = 0,
∂
∂t(ρ~u) + ∇ · (ρ~u ⊗ ~u) + ∇p = 0,
∂E
∂t + ∇ · (E~u + p~u) = 0,
∂α1
∂t + ~u · ∇α1 = 0,
(1)
as an example, for the principal motion of the state variables such as the partial densities, momentum, total energy, and volume fraction, respectively (cf. [1]). Here ρk and αk ∈ [0, 1] denote in turn the kth phasic density and volume fraction for k = 1, 2; α1 + α2 = 1. We have ρ = α1ρ1 + α2ρ2
representing the total density, ~u the vector of particle velocity, and p the mixture pressure.
To close the system, the constitutive law for each of the fluid phases of interest is assumed to satisfy a Mie-Gr¨uneisen equation of state,
pk(ρk, ek) = p∞,k(ρk) + ρkΓk(ρk) (ek− e∞,k(ρk)) , (2) where ek is the phasic specific internal energy, Γk = (1/ρk)(∂pk/∂ek)|ρk the Gr¨uneisen coefficient, and p∞,k, e∞,k are the properly chosen states of the
pressure and internal energy along some reference curve in order to match the experimental data of the material being examined. For simplicity, each of the expressions Γk, p∞,k, and e∞,k is taken as a function of the density only, see Section 5 for an example. As usual, E = ρe + ρ~u · ~u/2 is the total energy with the total internal energy defined as a volume-fraction average of the form ρe = P2
k=1αkρkek.
If the isobaric closure is assumed also in this model, where we have p1 = p2 = p in a region that contains more than one fluid component, from the total internal energy with (2), it is easy to derive the expression for mixture pressure as
p = ρe −
2
X
k=1
αkρke∞,k(ρk) +
2
X
k=1
αk
p∞,k(ρk) Γk(ρk)
! 2 X
k=1
αk
Γk(ρk). (3) With that, it can be shown that this five-equation model is hyperbolic when each physically relevant value of the state variables of the flow are defined in the region of thermodynamic stability, see [1] for the detail.
As reported in the literature, conventional numerical approach originally developed for single phase compressible flows can be used to solve this sys- tem. However, particular attention must be paid when computing the volume fraction. It is well known that even high order Eulerian transport scheme on fixed grids cannot completely remove numerical dissipation which then tends to continuously smear out the initial jump in the volume fraction func- tion. It is obviously problematic for immiscible interfaces. Even worse, it will eventually lead to the failure in identifying the moving interface which separate different materials for long term computations. So, we have to keep the interface sharp or at least recognizable throughout the simulation. To this end, numerical techniques, such as front tracking [5, 23, 45], anti diffu- sion [51, 47, 48], interface compression [42, 54], ALE (Arbitrary Lagrangian- Eulerian) and its variant [3, 4, 18, 25, 27, 46, 49], have been proposed.
In this work, our approach to interface-sharpening is a variant of the THINC (Tangent of Hyperbola for INterface Capturing) scheme that was proposed previously as an advection solver for the sharp resolution of a sub- grid discontinuity of volume fraction in the simulation of incompressible two- phase flow (cf. [2, 12, 60, 61, 62]). Shown in these works, the improved THINC method can get numerical accuracy comparable to the VOF schemes that use explicit geometrical reconstruction, such as PLIC(piecewise linear interface calculation) reconstruction, to retrieve the moving interfaces in in-
compressible multi-phase flows. Without the geometrical reconstruction, a THINC scheme is just a pure advection scheme, and is thus perhaps the simplest interface-capturing scheme in multi-phase flow simulations.
In the present case with the five-equation model for compressible two- phase flow, in regions near the interfaces, the original THINC scheme is used without any change for the reconstruction of volume fractions at each cell edge. Using the particular reconstruction function, THINC scheme effectively recovers the non-oscillatory sharp resolution of the volume fractions at every time step, and thus completely avoids the appearance of negative volume fractions that always accompanies the high-order Eulerian type schemes over fixed computational grids. With that, assuming the retain of constant cell averages of ρ1, ρ2, ~u, and p within a cell, a new reconstruction procedure is derived so as to ensure a consistent modeling of the remaining state variables such as α1ρ1, α2ρ2, ρ~u, and E at the cell edges, see [51] for a similar procedure that was employed in an anti-diffusion type interface-sharpening algorithm.
In regions away from the interfaces, where the flow is single phase, standard MUSCL (Monotone Upstream-centered Schemes for Conservation Laws) or WENO (Weighted Essentially Non-Oscillatory) reconstruction may be ap- plied for the interpolation of state variables to high order (cf. [20, 41, 57]);
this will be described further in Section 4.
To find approximate solutions of (1) with (2) and (3) numerically, we use a semi-discrete wave propagation method developed by Ketcheson et al. [15, 17] for general hyperbolic systems. In this method, the spatial dis- cretization is constructed by computing fluctuations from the solutions of Riemann problems with the proposed interface-sharpened initial data at cell edges. We employ the strong stability-preserving Runge-Kutta scheme [6, 7]
for the integration in time, yielding a simple and efficient implementation of the algorithm in the framework of the Clawpack software [16, 22].
It should be mentioned that our THINC-based interface-sharpening ap- proach is based on fixed underlying grids which is different from time-varying grid approaches such as the Lagrangian or ALE [3, 4, 18, 25, 27, 49], interface tracking [32, 38, 45], and adaptive moving grid method [53] for sharpening interfaces. In addition, our method is simpler than the ones proposed re- cently by Shukla et al. [42] and So et al. [51] (see also Appendix B) in the sense that there is no need to introduce additional interface-compression and anti-diffusion source terms to the compressible two-phase flow model, respectively, for the purpose of sharpening interfaces numerically (this may be difficult to do for more complex multiphase flow system and with phase
changes, see [40, 63], for example). Numerical results presented in Section 5 and Appendix A show the feasibility of the method for practical compressible two-phase flow and moving interface problems.
The format of this paper is outlined as follows. In Section 2, we review the basic idea of a semi-discrete finite volume method in wave-propagation form for hyperbolic conservation laws. In Section 3, we describe the THINC scheme for the sharp reconstruction of interfaces based on a given set of dis- crete volume fractions. The interface-sharpening algorithm for compressible two-phase flow is discussed in Section 4. Results of some sample validation tests as well as application of the method to practical problems are shown in Section 5. In Appendix A, additional advection benchmark tests are consid- ered, and in Appendix B, a brief review of an anti-diffusion based interface sharpening scheme is described.
2. Semi-discrete wave propagation method
We begin our discussion by reviewing a semi-discrete finite volume method in wave-propagation form [15, 17] that is essential in our interface-sharpening algorithm for solving the present quasi-conservative five-equation model nu- merically for compressible two-phase flow. For the ease of the latter descrip- tion, let us restate (1) in the form
∂q
∂t +
N
X
j=1
∂fj(q)
∂xj +
N
X
j=1
Bj(q) ∂q
∂xj = 0, (4a)
where with ~u = (u1, u2, . . . , uN) the terms q, fj, and Bj are defined by q = (α1ρ1, α2ρ2, ρu1, . . . , ρuN, E, α1)T , (4b) fj = (α1ρ1uj, α2ρ2uj, ρu1uj+ pδ1j, . . . , ρuNuj + pδN j, Euj+ puj, 0)T , (4c) Bj = diag (0, 0, . . . , 0, u1δ1j, . . . , uNδN j) , (4d) respectively, j = 1, 2, . . . , N. Here δij is the Kronecker delta, and N denotes the number of spatial dimension.
2.1. One-dimensional case
Consider the one-dimensional case N = 1 with a uniform grid with fixed mesh spacing ∆x1 and the number of cell in total M1 that discretize a spatial domain in the x1 coordinate. We use a standard finite-volume formulation
in which the value Qi(t) approximates the cell average of the solution over the grid cell Ci : x1 ∈ [x1i−1/2, x1i+1/2] at time t,
Qi(t) ≈ 1
∆x1
Z x1i+1/2 x1i−1/2
q(x, t) dx.
Let Q(t) = (Q1, Q2, . . . , QM1)T(t) be the vector of the approximate solu- tion of (4) at time t, and L1(Q(t)) = (L1(Q1), L1(Q2), . . . , L1(QM1))T (t) be the associated spatial-discretization vector. Then the semi-discrete version of the wave-propagation method is a method-of-lines discretization of (4) that can be written as a system of ODEs (Ordinary Differential Equations)
∂Q(t)
∂t = L1(Q(t)) , (5a)
where each component of L1(Q(t)) is defined by L1(Qi(t)) = − 1
∆x1 A+1∆Qi−1/2+ A−1∆Qi+1/2+ A1∆Qi
(5b) for i = 1, 2, . . . , M1. Here A+∆Qi−1/2 and A−∆Qi+1/2, are the right- and left-moving fluctuations, respectively, that are entering into the grid cell, and A1∆Qi is the total fluctuation within Ci. To determine these fluctuations, we need to solve Riemann problems.
Considering the fluctuations A±1∆Qi−1/2 arising from the edge (i − 1/2) between cells i − 1 and i, for example. This amounts to solving a Cauchy problem that consists of
∂q
∂t + ∂f1(q)
∂x1 + B1(q) ∂q
∂x1 = 0 (6a)
as for the equations and
q x1, t0 =
qi−1/2L if x1 < x1i−1/2
qi−1/2R if x1 > x1i−1/2, (6b) as for the initial condition at a time t0. Here qi−1/2L = limx→x1
(i−1/2)−q˜i−1(x) and qi−1/2R = limx→x1
(i−1/2)+ q˜i(x) are the interpolated states obtained by tak- ing limits of the reconstructed piecewise-continuous function ˜qi−1(x) or ˜qi(x)
(each of them are determined based on the set of discrete data {Qi(t0)}, see Sections 3 and 4) to the left and right of the cell edge at x1i−1/2, respectively.
If an approximate solver is used for the numerical resolution of the above Riemann problem (cf. [1, 55]), this would result in three propagating discon- tinuities that are moving with speeds s1,k and the jumps W1,k across each of them for k = 1, 2, 3, yielding the expression for the fluctuations as
A±1∆Qi−1/2 =
3
X
k=1
s1,k qi−1/2L , qRi−1/2±
W1,k qLi−1/2, qi−1/2R . (7a) In a similar manner, we may define the fluctuation A1∆Qi based on the Riemann problem with the initial data qi−1/2R and qi+1/2L at the cell center which gives
A1∆Qi =
3
X
k=1
s1,k qi−1/2R , qi+1/2L ±
W1,k qi−1/2R , qLi+1/2 ; (7b) this completes the definition of the fluctuations and also the spatial discretiza- tion operator L1(Qi(t)) in (5b). As usual, the notations for the quantities s± are set by s+ = max (s, 0) and s− = min (s, 0).
To integrate the system of ODEs (5a) in time, we employ the SSP (Strong Stability-Preserving) multistage Runge-Kutta scheme [6, 7]. That is, in the first-order case we use the Euler forward time discretization as
Qn+1 = Qn+ ∆tL1(Qn), (8a) where we start with the cell average Qn ≈ Q(tn) at time tn, yielding the solution at the next time step Qn+1 over ∆t = tn+1− tn. In the second-order case, however, we use the classical two-stage Heun method (or called the modified Euler method) as
Q∗ = Qn+ ∆tL1(Qn), Qn+1 = 1
2Qn+ 1
2Q∗+ 1
2∆tL1(Q∗). (8b)
It is common that the three-stage third-order scheme of the form Q∗ = Qn+ ∆tL1(Qn) ,
Q∗∗ = 3
4Qn+1
4Q∗+1
4∆tL1(Q∗) , Qn+1 = 1
3Qn+2
3Q∗∗+ 2
3∆tL1(Q∗∗) .
(8c)
is a preferred one to be used in conjunction with the fifth-order WENO scheme that is employed for the reconstruction of ˜qi(x1) during the spatial discretization (cf. [41]).
2.2. Multidimensional case
To extend the one-dimensional method described above to more space di- mensions, here we take a simple dimension-by-dimension approach. Namely, assuming a uniform Cartesian grid with fixed mesh spacing ∆x1 and ∆x2, in the x1-, and x2-direction, respectively, for example, the semi-discrete wave propagation method in two dimensions can be written of the form
∂Q(t)
∂t = L2(Q(t)) . (9a)
Here both Q(t) = {Qij(t)} and L2(Q(t)) = {L2(Qij(t))} are two-dimensional arrays for i = 1, 2, . . . , M1, j = 1, 2, . . . , M2, with Qij(t) denoting the ap- proximate cell average of the solution for the (i, j)th grid cell over the region Cij : (x1, x2) ∈ [x1i−1/2, x1i+1/2] × [x2j−1/2, x2j+1/2] at time t as
Qij(t) ≈ 1
∆x1∆x2
Z x1i+1/2 x1i−1/2
Z x2j+1/2 x2j−1/2
q(x, y, t) dy dx,
and L2(Qij(t)) representing the approximate spatial derivatives of the equa- tions (4) for Cij as
L2(Qij(t)) = − 1
∆x1 A+1∆Qi−1/2,j + A−1∆Qi+1/2,j+ A1∆Qij − 1
∆x2 A+2∆Qi,j−1/2+ A−2∆Qi,j+1/2+ A2∆Qij .
(9b)
Note that the fluctuations A+1∆Qi−1/2,j, A−1∆Qi+1/2,j and A1∆Qij are obtained by solving the one-dimensional Riemann problems in the direc- tion normal to the x1-axis (see (6)), while the fluctuations A+2∆Qi,j−1/2, A−2∆Qi,j+1/2 and A2∆Qij can be computed in a similar manner by solving the one-dimensional Riemann problems of in the direction normal to the x2-axis with
∂q
∂t + ∂f2(q)
∂x2 + B2(q) ∂q
∂x2 = 0
as for the equations and with chosen reconstructed data as for the initial condition. The SSP Runge-Kutta scheme can be employed also for the inte- gration of the semi-discrete scheme (9) in time.
It should be mentioned that the class of dimension-by-dimension semi- discrete method considered here would only give second order results when it is applied to solve general nonlinear hyperbolic systems, even if the fifth-order WENO reconstruction scheme is used in the spatial discretization, and the third order Runge-Kutta method is employed in time integration, see [64] for the details. We will not discuss the fully multidimensional spatial discretiza- tion in wave-propagation form here which is beyond the scope of this paper (see [20, 41] and references cited therein for more discussions on how this can be done), but will devote our attention to the devise of a new reconstruc- tion procedure for the purpose of sharpening interfaces from compressible two-phase flows.
2.3. Include source terms
To end this section, we comment that if x1 is the axisymmetric direction and u1 is the radial velocity, an axisymmetric version of the five-equation model (4) in one space dimension can be written as
∂q
∂t + ∂f1(q)
∂x1 + B1(q) ∂q
∂x1 = ψ(q), (10)
where ψ is the source term derived directly from the geometric simplification, yielding
ψ = −κ
x1 α1ρ1u1, α2ρ2u1, ρu21, Eu1+ pu1, 0T
.
Here in the case of a 2D radially or 3D spherically symmetric flow, the quantity κ takes 1 or 2, respectively.
Now, to find approximate solutions of (10) numerically, in this work, we take a simple unsplit approach (cf. [20]) in that the semi-discrete scheme (5) is employed again, but the term L1(Q(t)) on the right-hand side of the equation is modified to include the effect of source term ψ as
L1(Qi(t)) := L1(Qi(t)) + ψ(Qi(t)). (11) We then continue by employing an SSP Runge-Kutta method as usual to integrate the resulting ODEs in time for the updated solution at the next time step, see Section 5.1 for a sample result obtained using the method.
3. THINC reconstruction scheme
The reconstruction for any physical field to identify a sub-grid discon- tinuity can be derived using the THINC formulation for volume fraction function as a building block. To describe the basic idea, we start with the one-dimensional case, and for simplicity rather than using x1 for the spatial variable and α1(x1, t) for the volume fraction as in Section 2, we take the often-employed symbols x and φ(x, t) for that matter instead. Then as be- fore we define the approximate value of the cell-average of φ(x, t) over the grid cell Ci at time t by
φ¯i(t) ≈ 1
∆xi
Z xi+1/2
xi−1/2
φ(x, t) dx,
where ∆xi = xi+1/2 − xi−1/2 is the mesh size of the cell; assuming a non- uniform discretization of the spatial domain. Recall that being a volume fraction of a specified fluid, say fluid F , we have the values of ¯φi(t) as follows,
φ¯i(t) =
1 if Ci is filled with fluid F , α ∈ (0, 1) if an interface lies in Ci, 0 if there is no fluid F in Ci.
(12)
In practice, we identify an “interface cell” where the volume fraction satisfies, (i) ǫ < ¯φi(t) < 1 − ǫ,
(ii) ( ¯φi+1(t) − ¯φi(t))( ¯φi(t) − ¯φi−1(t)) > 0, (13) where ǫ is a small positive parameter (e.g., 10−4). Note that the latter condition (ii) is a monotonicity constraint on the data near the interfaces;
this is introduced as a measure to prevent from incurring spurious oscillations in the method.
Now given the volume fraction for the interface cell Ci, we seek a recon- struction that mimics the sub-grid structure of the jump between 0 and 1 in the volume fraction function. For a fixed time t and x ∈ [xi−1/2, xi+1/2], the following hyperbolic tangent function is well suited for this purpose,
Φi(x) = 1 2
1 + σitanh
β
x − xi−1/2
∆xi
− ˜xi
, (14a)
where σi = sgn(∆ ¯φi(t)) is the sign function of the variable ∆ ¯φi(t) = ¯φi+1(t)−
φ¯i−1(t), and β is a prescribed parameter to control the slope and thickness of
the jump which can be determined flexibly in practice (here we take β = 2.3 in all the tests considered in Section 5).
Note that the only unknown in (14a) is the center ˜xi. If we assume the conservation of volume fraction of the form
1
∆xi
Z xi+1/2 xi−1/2
Φi(x) dx = ¯φi(t),
it can be determined uniquely as
˜ xi = 1
2β ln
"
exp β 1 + σi − 2 ¯φi+1/2 /σi
1 − exp β 1 − σi− 2 ¯φi+1/2 /σi
#
. (14b)
It should be mentioned that using (14) as a building block, we may re- construct any physical field ψ(x, t) for x ∈ [xi−1/2, xi+1/2] where it has a jump
∆ψi = ψmax− ψmin by,
Ψi(x) = ψmin+ Φi(x)∆ψi. (15) Here lower and upper bounds of the jump ψmin and ψmax can be deter- mined from the neighboring cells such as ψmin = min(ψi−1/2L , ψi+1/2R ) and ψmax = max(ψLi−1/2, ψi+1/2R ), where ψLi−1/2 denotes the value computed from the reconstruction over the left-side cell [xi−3/2, xi−1/2] and ψi+1/2R from the right-side cell [xi+1/2, xi+3/2].
Reconstruction (15) can be applied to various physical variables, such as primitive variables, characteristic variables and flux functions. As long as the THINC reconstruction (15) is built up, we can get directly the values at the cell boundaries as a function of time t through either a point mapping or a trajectory integration for a wave propagation (transport) equation. In the latter instance, we consider the following wave equation as an example,
∂ψ
∂t + u∂ψ
∂x = 0,
where u is the characteristic velocity field. We split u into right-moving and left-moving parts, i.e., u = u++ u− where u+ = max(u, 0) and u− =
min(u, 0). The solution at the right-side of cell boundary xi−1/2 reads ψi−1/2R (t) = Ψi
xi−1/2− u−i−1/2t
= ψmin+ Φi
xi−1/2− u−i−1/2t
∆ψi
= ψmin+1 2
"
1 + σitanh β u−i−1/2t
∆xi
− ˜xi
!!#
∆ψi,
(16a)
while that at the left-side of cell boundary xi+1/2 reads ψi+1/2L (t) = Ψi
xi+1/2− u+i+1/2t
= ψmin+ Φi
xi+1/2− u+i+1/2t
∆ψi
= ψmin+ 1 2
"
1 + σitanh β 1 − u+i+1/2t
∆xi
− ˜xi
!!#
∆ψi,
(16b)
which can be used to any approximate Riemann solver in a full- or semi- discrete form, see the formulas given in [2] for that purpose also.
For the mass conservation law,
∂ψ
∂t +∂f
∂x = 0,
where f = uψ is the flux function, we have the total leftward flux across cell boundary xi−1/2 over t as,
fˆi−1/2− (t) = Z t
0
f (xi−1/2, t) dt =
Z xi−1/2 xi−1/2−u−i−1/2t
Ψi(x) dx
= ψ0u−i−1/2t − σi∆xi
2β ln
cosh
β
˜ xi+u
− i−1/2t
∆xi
cosh (β ˜xi)
∆ψi.
(17a)
Analogously, the rightward flux across xi+1/2 is fˆi+1/2+ (t) =
Z t 0
f (xi+1/2, t) dt =
Z xi+1/2 xi+1/2−u+i+1/2t
Ψi(x) dx
= ψ0u+i+1/2t − σi∆xi
2β ln
cosh
β
1 − ˜xi− u
+ i+1/2t
∆xi
cosh (β (1 − ˜xi))
∆ψi.
(17b) where ψ0 = (ψmin + ψmax)/2. So, a one-step finite volume formulation to transport mass can be written as,
ψ¯i(t) = ¯ψi(0) − 1
∆xi
ˆfi+1/2(t) − ˆfi−1/2(t) ,
which is used in the existing THNIC method for incompressible flows.
To extend the above THINC reconstruction scheme to more than one space dimension, here we take a simple dimension-by-dimension approach as proposed in Xiao et al. [60]; this one-dimensional version of the method works quite well with the semi-discrete wave propagation method described in Section 2 and also the interface-sharpening reconstruction scheme for com- pressible two-phase flow described next.
4. Interface sharpening algorithm
In each time step, our interface-sharpening algorithm for compressible two-phase flow problems consists of the following steps:
(1) Reconstruct a piecewise polynomial function, denoted by ˜q(xi, tn), for all xi based on standard MUSCL/WENO reconstruction procedure from the cell average Qn at time tn to more than first order.
(2) Modify ˜q(xi, tn) for interface cells using a variant of THINC scheme from Qn to sharpen the resolution of interfaces.
(3) Solve Riemann problems with interpolated initial data from ˜q(xi, tn) obtained in steps 1 and 2 for spatial discretization.
(4) Employ a semi-discrete method in wave propagation form to update Qn from the current time to the next Qn+1 over a time step ∆t.
Note that if step 2 is omitted in the above algorithm, it is simply the stan- dard semi-discrete method proposed by Ketcheson et al. [15, 17] for general hyperbolic systems with new applications to compressible two-phase flow, see Section 2. In this instance, the implementation of the MUSCL/WENO reconstruction scheme as suggested in [13, 20, 34] for step 1 is enough for that matter. Our goal here is to describe step 2 of the reconstruction procedure that is essential to our interface-sharpening algorithm.
4.1. Homogeneous-equilibrium-consistent reconstruction scheme
We begin step 2 by employing the one-dimensional THINC scheme de- scribed in Section 3 for the reconstruction of the volume fraction function q(N +4) = α1, in the five-equation model (4). Suppose that with Qn given at a time tn, we have obtained the interpolated volume fraction for the ith cell at the edges (α1)Ri−1/2 and (α1)Li+1/2 based on the relations (16a) and (16b) at t = 0, respectively.
With that, to construct the sub-grid structure of the partial density, we follow the approach proposed by So et al. [51] (see also Appendix B and [54]) in that the phasic densities ρ1 and ρ2 are assumed to remain constants within the cells during the reconstruction step, yielding the definition of the cell- edge states for the partial densities, i.e., the first and second component of q; q(k)= αkρk for k = 1, 2, as
(αkρk)Ri−1/2 = (αkρk)i+ (ρk)i(αk)Ri−1/2− (αk)i , (18a) (αkρk)Li+1/2 = (αkρk)i+ (ρk)i(αk)Li+1/2− (αk)i . (18b) Here (αk)i (ρk)i, and (αkρk)i are the cell-average variables obtained from Qni in an interface cell Ci.
Now to find the reconstructed states for the total momentum q(j+2) = ρuj
for j = 1, 2, . . . , N, the velocity field ~u = (u1, u2, . . . , uN) is assumed to remain unchanged within the cells also; this is true in a region for the contact discontinuity (the case we are interested in this work). Thus, we have
(ρ~u)Ri−1/2 = (ρ~u)i+ ~ui ρRi−1/2− ρi , (18c) (ρ~u)Li+1/2 = (ρ~u)i+ ~ui ρLi+1/2− ρi , (18d) where ρRi−1/2 =P2
k=1(αkρk)Ri−1/2 and ρLi+1/2 = P2
k=1(αkρk)Li+1/2 denotes the reconstructed mixture density at the left- and right-edge of the grid cell i, respectively.
Finally, to reconstruct the total energy q(N +3) = E, we assume further the pressure equilibrium within interface cells, yielding
Ei−1/2R = Ei+ Ki ρRi−1/2− ρi +
2
X
k=1
(ρkek)ih
(αk)Ri−1/2− (αk)ii
, (18e)
Ei+1/2L = Ei+ Ki ρLi+1/2− ρi +
2
X
k=1
(ρkek)ih
(αk)Li+1/2− (αk)ii
, (18f)
where Ki = ~ui · ~ui/2 denotes the specific kinetic energy in cell Ci. Note that since the above sub-grid reconstruction procedure (18) makes use of the basic homogeneous-equilibrium assumptions of the underlying contin- uum model (1), it should be pertinent to call it a homogeneous-equilibrium- consistent reconstruction scheme.
4.2. Interface only problem
To see how our algorithm works for sharpening interfaces, it is instruc- tive to consider an interface only problem where the initial data consists of uniform pressure p = p0, constant velocity ~u = ~u0, and constant phasic den- sities ρk = ρk0 for k = 1, 2 throughout the domain, while there are jumps on the other variables such as partial densities and volume fractions across some interfaces. Without loss of generality, we consider a one-dimensional problem with a positive velocity u1 = u0 > 0. In this case, suppose that a 3-wave HLLC approximate solver (cf. [9, 44, 55]) is used in step 3 of the algorithm for solving the Riemann problem (6) (similar result follows if a Roe solver [1] is employed instead), it is easy to find the fluctuations for each cell Ci, i = 1, 2, . . . , M1, as
A−1∆Qi+1/2 = 0,
A+1∆Qi−1/2 = u0 qi−1/2R − qi−1/2L , A1∆Qi = u0 qi+1/2L − qi−1/2R .
Inserting the above expression for fluctuations into (5b), we have the spatial discretization of this interface only problem:
L1(Qi(t)) = − 1
∆x1u0 qi+1/2L − qLi−1/2 .
Now if the Euler’s method (8a) is employed for the time integration in (5a), together with the above spatial discretization term the cell average Qni is updated by
Qn+1i = Qni − ∆t
∆x1u0 qLi+1/2− qi−1/2L , or equivalently by
α1ρ1 α2ρ2
ρu1
E α1
n+1
i
=
α1ρ1 α2ρ2
ρu1
E α1
n
i
− ∆t
∆x1u0
(α1ρ1)Li+1/2− (α1ρ1)Li−1/2 (α2ρ2)Li+1/2− (α2ρ2)Li−1/2
u0
ρLi+1/2− ρLi−1/2 Ei+1/2L − Ei−1/2L (α1)Li+1/2− (α1)Li−1/2
. (19)
Here we have assumed a consistent approximation of the true solution of the problem from the earlier time steps so that the equilibrium conditions pm = p0, um1 = u0, and ρmk = ρk0, k = 1, 2, are fulfilled for 0 ≤ m ≤ n.
With that, if we substitute the first two components of (19) into the third one, we arrive at readily the expected state of the particle velocity
(u1)n+1i = (u1)ni = u0.
When we continue applying this result to the fourth component of (19), after simple algebraic manipulations, we get the update of the total internal energy
2
X
k=1
(αkρkek)n+1i =
2
X
k=1
(αkρkek)ni − ∆t
∆x1u0
h(αkρkek)Li+1/2− (αkρkek)Li−1/2i ,
or alternatively
αk
p − p∞,k
Γk
+ ρke∞,k
n+1 i
=
αk
p − p∞,k
Γk
+ ρke∞,k
n i
− ∆t
∆x1u0·
"
αk
p − p∞,k
Γk
+ ρke∞,k
L i+1/2
−
αk
p − p∞,k
Γk
+ ρke∞,k
L i−1/2
# ,
when it is written componentwise for k = 1, 2 and with the Mie-Gr¨uneisen equation of state (2). Note that based on the results for αkρk and α1
from (19), we find the phasic density ρn+1k retains its expected value ρnk for k = 1, 2. Using this fact together with pn= p0 and the consistent reconstruc- tion of ρkek at each cell edges, it is not difficult to show the fulfillment of pressure equilibrium pn+1= p0 from the above equation for internal energy.
Having these results in mind, it should be easy to comprehend that the numerical resolution of volume fractions gives the basic characterization of the discontinuous interfacial profile for both density and total internal energy in which the THINC reconstruction scheme proposed here may be advanta- geous for some practical problems with immiscible interfaces; this will be justified numerically in Section 5.
5. Numerical results
We now present sample numerical results obtained using our semi-discrete wave propagation method with and without THINC-based solution recon- struction described in the previous sections for compressible two-phase flows in one and two space dimensions. In all the cases considered here, we have used a second-order MUSCL method with the HLLC approximate Riemann solver in step 1, a second-order SSP method in step 4, and the Courant num- ber ν = 1/2 (cf. [20]), while performing the runs. To demonstrate more proof of competitiveness of the approach, we have also included results obtained using the anti-diffusion interface sharpening technique (see Appendix B for a brief overview of the method).
5.1. One-dimensional case
Example 5.1. We begin by considering a simple interface only problem in one dimension where the exact solution of the problem consists of a square liquid column evolving in gas with uniform equilibrium pressure p = p0 = 105 Pa and constant particle velocity u1 = u0 = 100 m/s in a shock tube of one meter. In this test, inside the region x1 ∈ [0.4, 0.6]m the fluid is nearly liquid that contains a uniformly distributed small amount of the gas volume fraction α2 = 10−8, and it is nearly gas that contains the liquid volume fraction α1 = 10−8 elsewhere; the density for the pure liquid and gas phase in the domain are ρ1 = 1000 kg/m3 and ρ2 = 1kg/m3, respectively. We use the stiffened gas equation of state to model the thermodynamic behavior of liquid and gas where the material-dependent functions appeared in (2) are Γk = γk− 1, p∞,k = −γkBk, and e∞,k = 0 with the parameter values taken
in turn to be γ1 = 4.4, B1 = 6 × 108Pa, and γ2 = 1.4, B2 = 0 for the liquid and gas phases.
We carried out the computation using the proposed semi-discrete method with and without THINC reconstruction. Periodic boundary condition is used on the left and right boundaries during the computations. Figure 1 shows numerical results at time t = 10ms using a 100 mesh; for compari- son purposes results obtained using a wave-propagation based anti-diffusion method are included also. From the figure, it is easy to see that, with each of the methods, we have retained the correct pressure equilibrium and parti- cle velocity, without introducing any spurious oscillations near the interfaces.
Comparing with the exact solution, we observe an improved resolution of the interfaces for the density and volume fraction when the interface-sharpening method (with either the THINC- or anti-diffusion-based) is employed in the simulation.
Example 5.2. Our next example is a two-phase impact problem pro- posed by Saurel and Abgrall [39] (see also [43]). Initially, under the atmo- spheric condition (i.e., with uniform pressure p0 = 105 Pa and temperature T0 = 300K throughout the domain), there is a rightward going copper plate with the speed u1 = 1500 m/s interacting with a solid explosive at rest on the right of the plate, see Fig. 2 for illustration. In this problem, as in the previous test, there is a uniformly distributed small amount of volume frac- tion α2 = 10−8 and α1 = 10−8 inside the copper and explosive, respectively.
To model the material properties of the copper and (solid) explosive, we use the Cochran-Chan equation of state that in (2) we set the same Γk as in the stiffened gas case, but with p∞,k, e∞,k defined by
p∞,k(ρk) = B1k ρ0k
ρk
−E1k
− B2k ρ0k
ρk
−E2k
,
e∞,k(ρk) = −B1k
ρ0k(1 − E1k)
"
ρ0k
ρk
1−E1k
− 1
# + B2k
ρ0k(1 − E2k)
"
ρ0k ρk
1−E2k
− 1
#
− CvkT0.
(20)
Here γk, B1k, B2k, E1k, E2k, Cvk, and ρ0k are material-dependent quantities, see Table 1 for a typical set of numerical values for copper and explosive considered here (cf. [26, 29, 59]).
It is known that the exact solution of this impact problem is composed
0 0.2 0.4 0.6 0.8 1 0
0.2 0.4 0.6 0.8 1
ρ(Mg/m3 )
0 0.2 0.4 0.6 0.8 1
0 0.05 0.1 0.15
0.2 with THINC
no THINC with antiD exact
u1(km/s)
0 0.2 0.4 0.6 0.8 1
0 0.5 1 1.5 2
p(105 Pa)
x1 (m)
0 0.2 0.4 0.6 0.8 1
0 0.2 0.4 0.6 0.8 1
α1
x1 (m)
Figure 1: Numerical results for a passive advection of a square liquid column at time t= 10ms. The solid line is the exact solution and the points show the computed solution with 100 mesh points obtained using each of the methods with and without THINC re- construction; in each graph the anti-diffusion results (marked by triangles) are included also as for comparison.
of a leftward-going shock wave to the copper, a rightward-going shock waves to the inert explosive, and a material interface lying in between that sep- arates these two different materials (cf. [39]). We run this problem using the same numerical methods as performed in the previous examples for the single interface only problem but with a 200 grid. Figure 2 shows the re- sulting solutions at time t = 85µs for the density, velocity, pressure, and the copper volume fraction. By comparing the computed solutions with the fine grid one obtained using the THINC-based interface-sharpening method but with 5000 meshes, we see good agreement of the solution behaviors in all the methods under concerned with the correct shock speeds and free of spurious oscillations in the pressure near the interface. With the THINC reconstruction, we again observe some improvement on the structure of the density and volume fraction near the interface. Notice that there is a slight
Table 1: Material quantities for solid copper (k = 1) and explosive (k = 2) in the Cochran- Chan equation of state (20).
k ρ0k(kg/m3) B1k(GPa) B2k(GPa) E1k E2k γk Cvk (J/kg·K)
1 8900 145.67 147.75 2.99 1.99 3 393
2 1840 12.87 13.42 4.1 3.1 1.93 1087
overshoot on the phasic density α1ρ1 on the left of the interface; this error diminishes however when the the mesh is refined. This overshoot does not appear in the present anti-diffusion results.
Example 5.3. We continue our test by considering the interaction of a shock wave in molybdenum with an encapsulated MORB (Mid-Ocean Ridge Basalt) liquid. As in [30, 31, 43], we use (2) to model the thermodynamic behavior of MORB and molybdenum with Γk, p∞,k, and e∞,k defined as follows:
Γk(ρk) = Γ0k
ρ0k
ρk
ηk
,
p∞,k(ρk) = p0+
ρ0kc20k 1 −ρρ0k
k
h1 − ζk
1 −ρρ0k
k
i2, e∞,k(ρk) = e0 + 1
2ρ0k
1 − ρ0k
ρk
(p0+ p∞,k(ρk)) .
(21)
Here Γ0k = γk− 1 represents a reference Mie-Gr¨uneisen coefficient at ρk = ρ0k, ηk ∈ [0, 1] is a dimensionless parameter, c0k denotes the zero-pressure isentropic speed of sound, and ζk is a dimensionless parameter. With the zero initial reference state for p0 and e0, typical set of material quantities for the MORB (k = 1) and molybdenum (k = 2) are given in Table 2.
Table 2: Material quantities for MORB (k = 1) and molybdenum (k = 2) in the Mie- Gr¨uneisen equation of state (2) with (21).
k ρ0k(kg/m3) c0k(m/s) ζk γk ηk
1 2660 2100 1.68 1.18 1
2 9961 4770 1.43 2.56 1
0 0.2 0.4 0.6 0.8 1 0
2 4 6 8 10
0 0.2 0.4 0.6 0.8 10
0.5 1 1.5 2 2.5
α1ρ1(Mg/m3 ) α2ρ2(Mg/m3 )
0 0.2 0.4 0.6 0.8 1
0 0.5 1 1.5
with THINC no THINC with antiD fine grid
u1(km/s)
0 0.2 0.4 0.6 0.8 1
0 2 4 6 8 10 12
p(GPa)
x1 (m)
0 0.2 0.4 0.6 0.8 1
0 0.2 0.4 0.6 0.8 1
α1
x1 (m)
Figure 2: Numerical results for a two-phase (solid explosive-copper) impact problem at time t = 85µs. The solid line is the fine grid solution computed by 5000 meshes, and the points show the solutions with 200 meshes.
In this problem, the initial condition is composed of a stationary (MORB- molybdenum) interface at x1 = 0.6 m and a rightward going Mach 1.163 shock wave in molybdenum at x1 = 0.4 m traveling from left to right in a shock tube of unit length. The material on the right of the interface is a MORB liquid with the data
ρ1, ρ2, u1, p, α1)R= (ρ01, ρ02, 0, 0, 1 − 10−6 ,
and the material on the left of the interface, (i.e., on the middle and the preshock state), is molybdenum with data
ρ1, ρ2, u1, p, α1)M = (ρ01, ρ02, 0, 0, 10−6 . The state behind the shock in the molybdenum is
ρ1, ρ2, u1, p, α1)L= (ρ01, 11042 kg/m3, 543 m/s, 3 × 1010 Pa, 10−6 .
We note that this gives us one example, in which the (MORB-molybdenum) interface is accelerated by a shock wave coming from the heavy-fluid to the light-fluid region, and it is known that the resulting wave pattern after the interaction would consist of a transmitted shock wave, an interface, and a reflected rarefaction wave (cf. [10]).
Numerical results for this problem are shown in Fig. 3, where the snap shot of density, velocity, pressure, and volume fraction for MORB are shown at time t = 120µs using a 200 grid. From the figure, we again observe sensible improvement of the interface structure obtained using our THINC- based method; this result is comparable with the anti-diffusion one in all region, except in a narrow region on the left of the interface as in the previous case where a slight overshoot on α1ρ1 is present. Despite this, our THINC scheme works in a satisfactory manner for the other physical variables near the interfaces. A two-dimensional version of this problem will be considered in Example 5.8.
Example 5.4. To end this subsection, we are interested in a well-studied UNDEX (UNDerwater EXplosion) problem that involves the growth and collapse of an underwater gas bubble in spherically symmetric geometry.
Similar to the test case considered in [25] (see [58] also), the initial condition we take is composed of 300g of spherical TNT charge at a depth of 94.1m in water with ambient pressure of 1×107dyn/cm2. In this work, the constitutive law for both the explosive (k = 1) and water (k = 2) are described by the JWL (Jones-Wilkins-Lee) equation of state where in (2) we have
Γk = γk− 1,
p∞,k(ρk) = A1kexp −ρ0kR1k
ρk
+ A2kexp −ρ0kR2k
ρk
, e∞,k(ρk) = A1k
ρ0kR1k
exp −ρ0kR1k
ρk
+ A2k
ρ0kR2k
exp −ρ0kR2k
ρk
,
(22)
see Table 3 for a sample set of material quantities in cgs units for TNT and water (cf. [58]). With that, inside the explosive sphere of radius 3.5287cm (for a 300g TNT) the state variables are
ρ1, ρ2, u1, p, α1) = (ρ01, ρ02, 0, 83812408875.7788 dyn, 1 − 10−8 , and outside the sphere the state variables are
ρ1, ρ2, u1, p, α1) = (ρ01, ρ02, 0, 1 × 107 dyn, 10−8 .
0 0.2 0.4 0.6 0.8 1 0
2 4
0 0.2 0.4 0.6 0.8 10
10 20
α1ρ1(Mg/m3 ) α2ρ2(Mg/m3 )
0 0.2 0.4 0.6 0.8 1
0 0.2 0.4 0.6 0.8 1
with THINC no THINC with antiD fine grid
u1(km/s)
0 0.2 0.4 0.6 0.8 1
0 5 10 15 20 25 30 35
p(GPa)
x1 (m)
0 0.2 0.4 0.6 0.8 1
0 0.2 0.4 0.6 0.8 1
α1
x1 (m)
Figure 3: Numerical results for a shock wave in molybdenum interacting with an encap- sulated MORB liquid at time 120µs. The graphs are displayed in the same manner as in Fig. 2.
It should be mentioned that unlike the work done in [25, 58] where a hy- brid barotropic and non-barotropic equation of state (i.e., the Tait equation of state for water, and the JWL for TNT) is used in their numerical simu- lation of the problem with ALE-type methods, here with the five-equation model (4) and the definition of mixture pressure (3) it is convenient to con- sider a non-barotropic equation of state for water (such as the JWL) as well in our numerical method for simulation. It has been demonstrated in [58] and the results shown below that the influence on the equation of state for water is not high for this spherical UNDEX problem; the compressibility effect of water counts however.
We carry out the computation using the THINC-based interface-sharpening method proposed here together with the source term treatment described in Section (2.3) for spherical symmetry, where the solid wall is used on the left boundary, and the non-reflection boundary is used on the right. As in the
Table 3: Material quantities in cgs units for gaseous explosive (k = 1) and water (k = 2) in the JWL equation of state (22) for a spherical underwater explosion test.
k ρ0k(g/cm3) A1k(dyn/cm2) A2k(dyn/cm2) R1k R2k γk
1 1.63 3.712 × 1012 3.23 × 1010 4.15 0.95 1.3 2 1.00381 15.82 × 1012 −4.668 × 1010 8.94 1.45 2.172
results present in [25], Figs. 4(a)-(d) show the time sequence of the den- sity and pressure through four different solution phases, namely, the initial phase, the shock-interface wave interaction phase, the incompressible phase, and bubble collapse and rebound phase, respectively. First, note that we have used an adaptive domain size, i.e., 20, 60, 160, and 5000cm for the so- lution phases (a)-(d) individually, so that the leading shock wave remains in the computational domain; the mesh size is kept as a constant ∆x1 = 0.05cm throughout the runs. Comparing our results with the ones shown in [25], we see a qualitative agreement on the basic structure of the solution at all the phases.
To give a quantitative assessment of the solutions for this spherical UN- DEX problem, Fig. 5 shows the time history of the bubble radius and the bubble-water interface pressure obtained using both the THINC- and anti- diffusion-based methods (cf. [11, 25, 58]). From the figure, it is easy to estimate the maximum bubble radius (denoted by rmax) by checking the maximum of the discrete set of the bubble radius. Moreover, we may also es- timate the bubble period (denoted by Tb) by measuring either the minimum radius or the peak interface after the maximum volume expansion. Table 4 summarizes the findings together with those results reported in the literature as for comparison.
Note that the incompressible solution shown in the table is computed based on solving a variant of the Rayleigh-Plesset equation with the gas pressure on the bubble surface following the isentropic JWL equation of state, see [36, 58] for the details. From the table, it is clear that the incompressible solution does not give a good prediction on the solution which indicates that the water compressibility is of fundamental importance to this problem.
Comparing our results with those appeared in the literature (cf. [52, 25, 58]) and the anti-diffusion ones, we observe good agreement of results.
It should be mentioned that this spherical UNDEX problem is known to
