An Eulerian interface sharpening algorithm for compressible two-phase flow: The algebraic THINC
approach
Keh-Ming Shyue∗,a, Feng Xiaob
aDepartment of Mathematics, National Taiwan University, Taipei 106, 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 algo- rithm, in regions near the interfaces where two different fluid components are present within a cell, the THINC (Tangent of Hyperbola for INterface Cap- turing) reconstruction scheme is used as a basis for the reconstruction of a sub-grid discontinuity of volume fractions at each cell edge, and it is comple- mented by a homogeneous-equilibrium-consistent reconstruction 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 re- constructions are then used as the initial data for Riemann problems, and the resulting fluctuations form the basis for the spatial discretization. Time in- tegration of the algorithm is done by employing a strong stability-preserving Runge-Kutta method. Numerical results are shown for sample problems with the Mie-Gr¨uneisen equation of state for characterizing the materials of inter- ests in both one and two space dimensions that demonstrate the feasibility
∗Corresponding author
Email addresses: shyue@ntu.edu.tw(Keh-Ming Shyue ), xiao@es.titech.ac.jp (Feng Xiao)
of the proposed method for interface-sharpening of compressible two-phase flow.
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, 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 (cf. [20, 23]). 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.
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, 11, 47, 48, 49]). 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. 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 [41] 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 applied for the interpolation of state variables to high order (cf. [17, 31, 45]); 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. [13, 15] 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 [14, 18].
It should be mentioned that our THINC-based interface-sharpening ap- proach is simpler than the ones proposed recently by Shukla et al. [32] and So et al. [41] (see also [38, 39]) in the sense that there is no need to introduce additional interface-compression and anti-diffusion source terms to the com- pressible two-phase flow model, respectively, for the purpose of sharpening interfaces numerically. In addition, our method is based on fixed underly- ing grids which is different from time-varying grid approaches such as the Lagrangian or ALE (Arbitrary Lagrangian-Eulerian) [4, 5, 16, 19, 21, 40], interface tracking [26, 29, 37], and adaptive moving grid method [42] for sharpening interfaces. Numerical results presented in Section 5 show the feasibility of the method for practical compressible two-phase flow 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 solving the quasi-conservative five-equation model. In Section 3, we describe the THINC scheme for the sharp reconstruction of interfaces based on a given set of discrete 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.
2. Semi-discrete wave propagation method
We begin our discussion by reviewing a semi-discrete finite volume method in wave-propagation form [13, 15] 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, . . . , 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 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 xi−1/2, respectively.
If an approximate solver is used for the numerical resolution of the above Riemann problem (cf. [1, 43]), 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 similiar 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 strong stability-preserving (SSP) 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 a time interval ∆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. [31]).
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, 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 ob- tained by solving the one-dimensional Riemann problems in the direction nor- mal to the x1-axis (cf. (6)), while the fluctuations A+2∆Qi,j−1/2, A−2∆Qi,j+1/2
and A2∆Qij are 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 (this is the case considered here), 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 [50] for the details. We will not discuss the fully multidimensional spatial discretization in wave-propagation form here which is beyond the scope of this paper (see [17, 31] and references cited therein for more discussions on how this can be done), but will devote our attention to the devise of a new reconstruction procedure for the purpose of sharpening interfaces from compressible two-phase flows.
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 here. 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.
(10)
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, (11) 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 safe measure to prevent from incurring spurious oscillations in the method (cf. [2]).
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
, (12a)
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 (12a) 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
#
. (12b)
It should be mentioned that using (12) 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 = ψimax− ψmini by,
Ψi(x) = ψmini + Φi(x)∆ψi. (13) Here lower and upper bounds of the jump ψimin and ψmaxi can be deter- mined from the neighboring cells such as ψimin = min(ψi−1/2L , ψi+1/2R ) and ψmaxi = max(ψi−1/2L , ψRi+1/2), where ψi−1/2L 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 (13) can be applied to various physical variables, such as primitive variables, characteristic variables and flux functions. As long as the THINC reconstruction (13) is built up, we can get directly the values at the cell boundaries as a function of time t through either a point map- ping or a trajectory integration for a wave propagation (transport) equation.
Considering the latter instance, we take 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
ψRi−1/2(t) = Ψi
xi−1/2− u−i−1/2t
= ψimin+ Φi
xi−1/2− u−i−1/2t
∆ψi
= ψimin+ 1 2
"
1 + σitanh β u−i−1/2t
∆xi
− ˜xi
!!#
∆ψi,
(14a)
while that at the left-side of cell boundary xi+1/2 reads ψLi+1/2(t) = Ψi
xi+1/2− u+i+1/2t
= ψmini + Φi
xi+1/2− u+i+1/2t
∆ψi
= ψmini + 1 2
"
1 + σitanh β 1 −u+i+1/2t
∆xi
− ˜xi
!!#
∆ψi,
(14b)
which can be used to any approximate Riemann solver in a full- or semi- discrete form.
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
= ψiaveu−i−1/2t − σi∆xi
2β ln
cosh
β
˜ xi+u
− i−1/2t
∆xi
cosh (β ˜xi)
∆ψi.
(15a)
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
= ψiaveu+i+1/2t − σi∆xi
2β ln
cosh
β
1 − ˜xi− u
+ i+1/2t
∆xi
cosh (β (1 − ˜xi))
∆ψi.
(15b) where ψiave = (ψimin + ψimax)/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. [47]; 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 or WENO reconstruction procedure from the cell average Qn at time tn to more than first order in space.
(2) Modify ˜q(xi, tn) for interface cells using a variant of THINC scheme from Qn to sharp 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. [13, 15] 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 [12, 17, 27] 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 fifth component of q, namely, the volume fraction function q(5) = α1, in the five-equation model (4). Sup- pose 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 (14a) and (14b) 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. [41] in that the phasic densities ρ1 and ρ2 are assumed to remain constants within the cells during the recon- struction 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 , (16a) (αkρk)Li+1/2 = (αkρk)i+ (ρk)i(αk)Li+1/2− (αk)i . (16b) Here (αk)i, (ρk)i, and (αkρk)i are the cell-average variables obtained from Qni in an interface cell Ci. For simplicity, we have surpressed the superscript n in (16a), and will continue to do so in the remaining subequations of (16).
Now to find the reconstructed states for the total momentum q(3) = ρ~u, the velocity field ~u 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 , (16c) (ρ~u)Li+1/2 = (ρ~u)i+ ~ui ρLi+1/2− ρi , (16d) 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(4) = 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
, (16e)
Ei+1/2L = Ei+ Ki ρLi+1/2− ρi +
2
X
k=1
(ρkek)ih
(αk)Li+1/2− (αk)ii
, (16f)
where Ki = ~ui · ~ui/2 denotes the specific kinetic energy in cell Ci. Note that since the above sub-grid reconstruction procedure (16) 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 better to consider an interface only problem where the initial data consists of uniform pressure p = p0 and constant velocity ~u = ~u0 throughout the domain, while there are jumps on the other variables such as volume fractions and material- dependent quantities 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, 36, 43]) is used in step 3 of the algorithm for solving the Riemann problem (6), it is easy to find the fluctuations for each cell Ci 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
. (17)
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 and um = u0 are fulfilled for 0 ≤ m ≤ n.
With that, if we substitute the first two components of (17) 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 (17), after simple algebraic manipulations, we get
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 expressing componentwise for each phase k = 1, 2, and employing the Mie-Gr¨uneisen equation of state (2) in it. Note that based on the results for αkρk and α1 from (17), we find the phasic density ρn+1k retains its expected value ρnk for k = 1, 2. Using this fact together with pn= p0, it is not difficult to show the maintain of pressure equilibrium pn+1 = p0 from the proposed method.
Having the above results in mind, it should be easy to comprehend that the numerical resolution of volume fractions gives the basic characterization of the discontinuous profile for both density and total internal energy in which the THINC reconstruction scheme plays an important role towards sharpening interfaces; this will be justified numerically in Section 5. We note that the same result shown in this subsection will follow also if a Roe solver [1] is employed in the method instead (cf. [33, 34]).
5. Numerical results
We now present sample numerical results obtained using our semi-discrete wave propagation method with and without interface-sharpening 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 in step 1, and a second-order SSP method in step 4, while performing the runs.
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 = 102 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 = 103kg/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, 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 method with and without interface sharpening. 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 and the Courant number ν = 1/2 (cf. [17]). From the figure, it is easy to see that, with each of the methods, we have retained the correct pressure equilibrium and particle 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 interface-sharpening is employed in the simulation.
Example 5.2. Our next example is a two-phase impact problem pro- posed by Saurel and Abgrall [30] (see also [35]). 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 where in (2) we set the same Γk as in
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 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 shows the computed solution with 100 mesh points obtained using each of the methods with and without interface- sharpening.
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.
(18)
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. [20, 23, 46]).
It is known that the exact solution of this impact problem is composed of
Table 1: Material quantities for copper (k = 1) and explosive (k = 2) in Cochran-Chan equation of state (18).
k ρ0k(kg/m3) B1k(GPa) B2k(GPa) E1k E2k γk Cvk
1 8900 145.67 147.75 2.99 1.99 3 393J/kg·K
2 1840 12.87 13.42 4.1 3.1 1.93 1087J/kg·K
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 separates these two different materials (cf. [30]). We run this problem using exactly the same method as performed in the previous examples for the single interface only problem but with a 200 grid. Figure 2 shows the resulting solution at time t = 85µs for the partial densities, velocity, pressure, and the copper volume fraction. By comparing the computed solutions with the fine grid one obtained using the interface-sharpening method but with 5000 meshes, we see good agreement of the solutions with the correct shock speeds and free of spurious oscillations in the pressure near the interface. With our interface- sharpening method, we again observe some improvement of the the structure of the solution for the density and volume fraction near the interface. Notice that there is a slight overshoot on the partial density α1ρ1 (shown in blue dots) behind the interface; this error diminishes however when the the mesh is refined.
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 [24, 25, 35], 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)) .
(19)
Here Γ0k = γk− 1 represents a reference Mie-Gr¨uneisen coefficient at ρk =
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 exact initial
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 solution with 200 meshes; the plot on the top-left corner is for the partial densities where the blue dots are for α1ρ1, abd the green dots are for α2ρ2.
ρ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 sets of material quantities for the MORB (k = 1) and molybdenum (k = 2) are given in Table 2.
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 − ε) ,
and the material on the left of the interface, (i.e., on the middle and the
Table 2: Material quantities for MORB (k = 1) and molybdenum (k = 2) in (19).
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
preshock state), is molybdenum with data
(ρ1, ρ2, u1, p, α1)M = (ρ01, ρ02, 0, 0, ε) . 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, ε , where ε = 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. [3, 10]).
Numerical results for this problem are shown in Fig. 3 at time t = 120µs for the states α1ρ1, α2ρ2, u1, p, and α1. We again observe sensible improve- ment of the interface structure, without introducing spurious oscillation in the pressure, while maintaining the correct behavior of the solution for shock and rarefaction waves. A two-dimensional version of this problem will be considered in Example 5.6.
5.2. Two-dimensional case
Example 5.4. To show how our method works in two dimensions, we begin by considering a simple interface only problem in a unit square do- main where the exact solution consists of a passive evolution of a square liquid column in gas of size (x1, x2) ∈ [0.3, 0.7] × [0.3, 0.7]m2 with uniform equilibrium pressure p = 105 Pa and constant particle velocity (u1, u2) = (102 m/s, 102 m/s). In this test, inside the square region (x1, x2) ∈ [0.3, 0.7]×
[0.3, 0.7]m2 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 = 103kg/m3 and ρ2 = 1kg/m3,
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 exact
u1(km/s)
0 0.2 0.4 0.6 0.8 1
0 5 10 15 20 25 30
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.
respectively. We use the linearized Mie-Gr¨uneisen equation of state to model the liquid and gas where the material-dependent functions appeared in (2) are
Γk = γk− 1, p∞,k(ρk) = c20k(ρk− ρ0k) , e∞,k = 0,
with the parameter values taken in turn to be γ1 = 4.4, c01 = 1624.8m/s, ρ01 = 103kg/m3, and γ2 = 1.4, c02 = 0m/s, ρ02 = 1kg/m3. Here periodic boundary condition is used in all sides.
Figure 4 shows contour plots of the density obtained using the method with and without interface-sharpening at time t = 20ms (this is the time the water column travelled over two periodic distance of the domain) using a 100 × 100 grid. Excellent interface-sharpening result is observed, whereas severely diffused interface is seen using the standard interface-capturing method.
It is not difficult to show that both the pressure and the velocity retains their equilibrium states, see Section 4.2 as an example.
With THINC No THINC
Figure 4: Numerical results for a passive evolution of a square liquid column in gas.
Contours of the density are shown at time t = 20ms obtained using the method with and without interface-sharpening on a 100 × 100 grid.
Example 5.5. We are next concerned with a well-known shock-bubble interaction problem that involves the collision of a shock wave in air with a circular R22 gas bubble (cf. [8, 16, 22, 28, 37] and cited references therein).
Initially, there is a planar leftward-moving Mach 1.22 shock wave in air lo- cated at x1 = 275mm traveling towards a stationary R22 gas bubble with center (x10, x20) = (225, 44.5)mm and of radius r0 = 25mm lying in front of it. Here both the air and R22 are modeled as perfect gases so that in (2) we have
Γk= γk− 1, p∞,k = 0, e∞,k = 0,
where γ1 = 1.249, ρ01 = 3.863 kg/m3 and γ2 = 1.4, ρ02 = 1.225 kg/m3 are used in turn for R22 (k = 1) and air (k = 2). Inside the R22 gas bubble, the state variables are
(ρ1, ρ2, u1, u2, p, α1) = (ρ01, ρ02, 0, 0, 1.01325 × 105 Pa, 1 − ε), while outside the bubble they are
(ρ1, ρ2, u1, u2, p, α1) = (ρ01, ρ02, 0, 0, 1.01325 × 105 Pa, ε) and
(ρ1, ρ2, u1, u2, p, α1) = (ρ01, 1.686 kg/m3, −113.5 m/s, 0, 1.59 × 105 Pa, ε) in the preshock and postshock regions, respectively, where ε = 10−8. The computational domain considered here is a shock tube of size: (x1, x2) ∈
[0, 445] × [0, 89]mm2, where the boundary conditions are solid wall on the top and bottom, and non-reflecting on the left and right.
As before, we carry out the computations using each of the methods with and without interface sharpening. Figure 5 shows a sequence of the schlieren-type images of the density at eight different times t = 55, 115, 187, 247, 318, 342, 417, and 1020µs (measured relative to the time where the incident shock first hits the upstream bubble wall) with a 3560 × 356 grid. Here for clarity, only the results in some short distances around the gas bubble are presented. From the figure, it is easy to observe the sensible improvement on the sharpness of the interface when the interface-sharpening version of the method is employed in the runs.
To get a quantitative assessment of several prominent features of the flow, following [28], in Fig. 6 we make a diagnosis plot of the space-time locations of the incident shock wave, the upstream bubble wall, the downstream bubble wall, the refracted shock wave, and the transmitted shock wave at some selected times up to t = 250 µs, see [37] for the details on how this graph is constructed. With these trajectories, we perform a linear least-squares fit to each set of points separately, and take the slope of the respective line as one measure of the wave speed of interests. Table 3 provides a comparison of the various computed velocities between our results with those appeared in the literature. Due to the sensitivity of the interface structure on the grid resolution and numerical methods used in the computation as well as the uncertainty on the experimental error, it should be fair to say that we have come up with a reasonable set of data for the bubble velocities in this application.
Here Vs represents the speed of the incident shock wave during the time t ∈ [0, 250]µs, Vui(Vuf) is the speed of the initial (final) upstream bubble wall when time t ∈ [0, 400]µs (t ∈ [400, 1000]µs), Vdi(Vdf) is the speed of the initial (final) downstream bubble wall when time t ∈ [200, 400]µs (t ∈ [400, 1000]µs), VRis the speed of the refracted shock when time t ∈ [0, 202]µs, and VT is the speed of the transmitted shock when time t ∈ [202, 250]µs.
Example 5.6. Finally, we consider a test problem proposed by Miller and Puckett [25] in which a shock wave in molybdenum is interacting with a region of encapsulated MORB liquid. Similar to the initial condition used in Example 5.3, at x1 = 0.3 m, there is a planarly rightward-moving Mach 1.163 shock wave in molybdenum traveling from left to right that is about to collide with a rectangular region [0.4, 0.7] × [0, 0.5] m2 which contains a MORB liquid inside. As before, we use the Mie-Gr¨uneisen equation of
With THINC No THINC
Figure 5: Numerical results for a planar Mach 1.22 shock wave in air interacting with a circular R22 gas bubble. A sequence of the schlieren-type images of the density obtained using each of the methods with and without THINC interface-sharpening is shown (from top to bottom) at eight different times t = 55, 115, 187, 247, 318, 342, 417, and 1020µs
With THINC No THINC
Figure 5: (Continued)
Table 3: A comparison of the computed velocities obtained using our algorithm for Ex- ample 5.5 with those reported in the literature; see the text for the definition of the notations used in the table.
Velocity (m/s) Vs VR VT Vui Vuf Vdi Vdf
Experiment [8] 415 240 540 73 90 78 78
Quirk and Karni [28] 420 254 560 74 90 116 82 Kokh and Lagoutiere [16] 411 243 525 65 86 86 64 Ullah et al. [44] 410 246 535 65 86 76 60
Shyue [37] 411 243 538 64 87 82 60
Our results (no THINC) 410 244 536 65 86 98 76 Our results (with THINC) 410 244 538 65 86 87 64
0 20 40 60 80
0 50 100 150 200
250 With THINC
x (mm)
t(µs)
0 20 40 60 80
0 50 100 150 200
250 No THINC
x (mm)
t(µs)
Figure 6: Space-time locations of the incident shock wave (marked by symbol ”×”), the upstream bubble wall (marked by symbol ”+”), the downstream bubble wall (marked by symbol ”⋄”), the refracted shock (marked by symbol ”△”), i and the transmitted shock (marked by symbols ”∗” and ”△”) for the runs shown in Fig. 5. These trajectories can be used to estimate the speed such as Vs, Vu, Vd, VR, VT in the order mentioned above.
