行政院國家科學委員會專題研究計畫 成果報告
具毛細管效應之可壓縮性流體數值計算法
研究成果報告(精簡版)
計 畫 類 別 ： 個別型
計 畫 編 號 ： NSC 95-2115-M-002-010-
執 行 期 間 ： 95 年 08 月 01 日至 96 年 10 月 31 日
執 行 單 位 ： 國立臺灣大學數學系暨研究所
計 畫 主 持 人 ： 薛克民
計畫參與人員： 大學生-兼任助理：林保同、張紘睿、李念達
處 理 方 式 ： 本計畫可公開查詢
中 華 民 國 97 年 04 月 02 日
A wave-propagation based generalized Lagrangian method for
compressible multicomponent flow problems
Keh-Ming Shyue
Department of Mathematics, National Taiwan University, Taipei 106, Taiwan
E-mail address: [email protected]
Abstract
A simple Lagrange-like moving grid approach is developed for the efficient numerical resolution of inviscid compressible multicomponent flow problems with a stiffened gas equation of state in more than one space dimension. The algorithm uses a variant of the fluid-mixture formulation of equations (cf. K.-M. Shyue, An efficient shock-capturing algorithm for compressible multicomponent problems, J. Comput. Phys. 142 (1998) 208–242) that are written in the generalized curvilinear coordinate as a basis to the modeling of the numerically induced mixing between two different fluid components within a computationally logical rectangular grid cell. In the algorithm, the temporal evolution of the physical grid system is chosen to follow approximately to the underlying flow fields, and a set of geometric conservation laws is included to describe the basic motion of the grid metrics associated with the coordinate transformation between the physical and logical grid. A standard flux-based wave decomposition method devised by Bale et al. (D. Bale, R.J. LeVeque, S. Mitran, J.A. Rossmanith, A wave propagation method for conservation laws and balance laws with spatially varying flux functions, SIAM J. Sci. Comput. 24 (2002) 955–978) is employed to solve the proposed multicomponent model with the dimensional-splitting technique incorporated in the method for multidimensional problems. Several numerical results are presented in one and two space dimensions that show the feasibility of the algorithm to the improvement of numerical resolution of interfaces when it is applied to solve a sample class of problems of practical importance.
AMS: 65M06, 76L05, 76M20, 76T05
Keywords: Generalized Lagrangian method; Flux-based wave decomposition; Compressible multicompo-nent flow; Stiffened gas equation of state; Spatially varying fluxes
1
Introduction
Our goal of this paper is to present a simple moving grid approach for the efficient numerical resolution of multicomponent flow problems with general inviscid compressible materials in more than one space dimension. We consider a simplified two-fluid problem as an example for the basic method development, where the flow regime of interest is assumed to be homogeneous with no jumps in the pressure and velocity (the normal component of it) across a material interface that separates regions of two different fluid components in an Nd≥ 1 spatial domain. The algorithm uses the Euler equations of gas dynamics
as a model system for the principal motion of each fluid component in that the conservation form of the equations in an Nd-dimensional Cartesian coordinate ~X = (x1, x2, . . . , xNd) and time t can be written as
∂ ∂t ρ ρui E + Nd X j=1 ∂ ∂xj ρuj ρuiuj+ pδij Euj+ puj = 0 for i = 1, 2, . . . , Nd, (1)
where ρ, ui, p, E, and δij denote the density, the particle velocity in the xi-direction, the pressure, the
total energy, and the Kronecker delta function, respectively. To complete the model, a linearized version of the Mie-Gr¨uneisen (i.e., the linearly density-dependent stiffened gas) equation of state of the form
p(ρ, e) = (γ − 1) ρe + (ρ − ρ0) B (2)
is assumed for the fundamental characterization of the thermodynamic behavior of the fluid component of concerned (cf. [24, 31, 33, 34]). Here e, γ, ρ0, and B are in turn the specific internal energy, the
adiabatic constant (γ > 1), the reference density, and the reference speed of sound. With that, the reference pressure is defined by p0= ρ0B/γ, and we have E = ρe + ρPNj=1d u2j/2 as usual.
To solve this homogeneous two-fluid problem numerically, one popular approach among them is to take a uniform Cartesian grid and employ a fluid-mixture type interface-capturing method proposed previously by the author (cf. [44, 45, 46, 48, 49]). In this approach, to model grid cells that contain more than one fluid component, the method may be based on a volume-fraction formulation of the model system that is formed by combining the Euler equations (1) for the basic conservative variables of the fluid mixtures and an additional set of the effective equations
∂ ∂t ρB γ − 1 + Nd X j=1 ∂ ∂xj ρB γ − 1uj = 0, ∂α ∂t + Nd X j=1 uj ∂α ∂xj = 0, (3)
for the mixtures of the material-dependent variables ρB/(γ − 1) and α, respectively. Here α ∈ [0, 1] denotes the volume fraction of the fluid component within a grid cell. It should be mentioned that the above effective equations are introduced in the method primarily for a direct computation of the pressure from the equation of state,
p = " E − PNd j=1(ρuj)2 2ρ + ρB γ − 1 − 2 X ι=1 αι ρ0ιBι γι− 1 # 2 X ι=1 αι γι− 1, (4)
and are derived so as to ensure a consistent modeling of the energy equation near the interfaces where two or more fluid components are present in a grid cell, and also the fulfillment of the mass equation in the other single component regions, see [48] for the details. Note that in (4) we have set α1 = α and
α2 = 1 − α for the volume fractions occupied by the fluid-component 1 and 2, in a respective manner,
and the mixtures of γ and ρ0B = γp0 are computed by
γ = 1 + 1 2 X ι=1 αι γι− 1 and ρ0B = 2 X ι=1 αι ρ0ιBι γι− 1 2 X ι=1 αι γι− 1 .
A standard high-resolution wave propagation method (cf. [26]) can then be used to find approximate solutions of the model system for multicomponent problems, yielding results that are free of spurious oscillations in the pressure near a somewhat smeared interface, see Section 6 for an example.
In this work, motivated by the success of a unified coordinate approach advocated by Hui and cowork-ers for single component flow problems with interfaces and complex geometries (cf. [16, 17, 18, 19]), we are interested in a novel moving grid method that is a generalization of the aforementioned fluid-mixture method and is also an easy extension of the unified coordinate approach for multicomponent flows. Here as a first endeavor, we focus our attention on the improved numerical resolution of interfaces, but leave the problems with complex geometries in a sequel. In this instance, a simple Lagrange-like condition of the form
∂ ~X
is assumed as the basis for the temporal evolution of the physical grid coordinate ~X, where ~u = (u1, u2, . . . , uN d) is the vector of the particle velocity, and h ∈ [0, 1] is a freely chosen parameter
that takes on the value of zero when we have a fixed Eulerian grid, and on the value of unity when we have a Lagrangian grid, see [5, 9, 28, 29, 30, 41] and the references cited therein for some background in-formation of a Lagrangian-type method. Note that, in practice, we have chosen a fixed constant h = 0.99 in almost all the numerical tests carried out in Section 6, but do not consider the more general space-and time-dependent case (cf. [17, 19, 20]) in the present study.
With the mesh-moving condition (5), the algorithm employs a variant of the model equations (1) and (3) that are written in a generalized curvilinear coordinate as a basis to the modeling of the numeri-cally induced mixing between two different fluid components within a computational (logical) rectangular grid cell. In addition to that, a set of geometric conservation laws is included to describe the basic motion of the grid metrics associated with the coordinate transformation between the physical and computational grid. This will be described further in Section 2.
We use a finite volume method based on an f -waves decomposition viewpoint developed by Bale et al. [2] to find approximate solutions of our proposed multicomponent model in generalized coordinates with the dimensional-splitting technique incorporated in the method for multidimensional problems. The method is a variant of the Godunov scheme, and has been widely used in many applications including the single-component fluid of ideal gases on general quadrilateral grids [26], and hyperbolic systems on curved manifolds [40] and on spheres [39]. The current use of the method extends the previous one from a single-component to multicomponent problems and also from the equations to be solved in the Cartesian coordinate in the physical space to the generalized coordinate in the logical space. It is an efficient and yet accurate scheme without introducing any spurious oscillations in the pressure near an interface as illustrated by numerical results presented in this paper. The basic idea of this method will be reviewed briefly in Section 4.
It should be mentioned that due to the Lagrangian nature of the above moving grid condition in problems when there is a large interfacial deformation or collapse of the underlying velocity field into a narrow region (cf.[46, 50]), for example, the grid system that is generated by the algorithm during the computation would be suffering from the usual grid tangling and so eventually the break down of the method in some instances. Despite this, there are some problems of sufficient interest that the current approach is worthwhile, particularly since it is relatively easy to implement the method in more than one space dimension and can provide a better resolution of interfaces than standard Eulerian methods. Moreover, the methodology we have proposed here can be extended quite straightforwardly to problems involving more than two fluid components and more complicated equations of state. Without going into the details for that, our goal is to establish the basic solution strategy and validate its use via some sample numerical experimentations.
This paper is organized as follows. In Section 2, we describe our mathematical model for the present homogeneous two-fluid problems in a generalized curvilinear coordinate. In Section 3, we discuss the hyperbolicity of the proposed model equations. In Section 4, we give a brief review of the dimensional-splitting version of the high-resolution flux-based wave decomposition method for multi-dimensional prob-lems. In Section 5, we explain the Riemann problem and a shock-only Riemann solver that is fundamental in our algorithm for numerical approximation. Numerical results of some sample examples in one and two space dimensions are presented in Section 6.
2
Mathematical models in generalized coordinates
We begin our discussion by describing a general coordinate transformation of the aforementioned fluid-mixture type model system, i.e., the combined equations (1) and (3), from the Cartesian coordinate
−1 0 1 −1.5 −1 −0.5 0 −1 −0.5 0 0.5 −2 −1.5 −1 −0.5 x1 x2 ξ1 ξ2 Ω mapping−→ Ωˆ ξ1= ξ1(x1, x2) ξ2= ξ2(x1, x2) logical domain physical domain
Figure 1: An example of a general non-rectangular domain Ω in two dimensions on the left that is mapped to a logical domain ˆΩ on the right via the transformation (6).
~
X = (x1, x2, . . . , xNd) in a physical domain Ω to the generalized coordinate ~Ξ = (ξ1, ξ2, . . . , ξNd) in a
computationally logical domain ˆΩ via the relations
τ = t, ξj= ξj( ~X, t) for j = 1, 2, . . . , Nd (6)
that is fundamental in our moving grid algorithm for constructing approximate solutions of the present homogeneous two-fluid problems with the stiffened gas equation of state (2), see Fig. 1 for an illustration of the spatial domain of interest when Nd= 2. In this instance, using the chain rule of partial differentiation,
we have ∂ ∂t = ∂ ∂τ + Nd X i=1 ∂ξi ∂t ∂ ∂ξi , ∂ ∂xj = Nd X i=1 ∂ξi ∂xj ∂ ∂ξi for j = 1, 2, . . . , Nd. (7)
Here ∂ξi/∂t and ∂ξi/∂xjappearing in the above equations for i, j = 1, 2, . . . , Nd, are the so-called metrics,
and have Nd(Nd+1) of them in total. For computational purposes, it is often desirable to find expressions
that relate the above metrics in the logical domain to the derivatives in the physical domain. We can do this quite straightforwardly (cf. [1, 15, 54]) if the existence of the inverse transformation of (6) is assumed, i.e.,
t = τ, xj= xj(~Ξ, t) for j = 1, 2, . . . , Nd, (8)
and the elementary differential rule from multi variable calculus as
∂(τ, ~Ξ) ∂(t, ~X) = ∂(t, ~X) ∂(τ, ~Ξ) −1 (9)
is used in the derivation. Note that in matrix form ∂(τ, ~Ξ)/∂(t, ~X) and ∂(t, ~X)/∂(τ, ~Ξ) are given by
1 0 . . . 0 ∂tξ1 ∂x1ξ1 . . . ∂xNdξ1 ∂tξ2 ∂x1ξ2 . . . ∂xNdξ2 .. . ... ... ... ∂tξNd ∂x1ξNd . . . ∂xNdξNd and 1 0 . . . 0 ∂τx1 ∂ξ1x1 . . . ∂ξNdx1 ∂τx2 ∂ξ1x2 . . . ∂ξNdx2 .. . ... ... ... ∂τxNd ∂ξ1xNd . . . ∂ξNdxNd ,
respectively, where to simplify the presentation we have used the notation ∂az = ∂z/∂a for any pair of
As an example, we consider the three-dimensional case Nd= 3, and so from (9) we would have 1 0 0 0 ∂tξ1 ∂x1ξ1 ∂x2ξ1 ∂x3ξ1 ∂tξ2 ∂x1ξ2 ∂x2ξ2 ∂x3ξ2 ∂tξ3 ∂x1ξ3 ∂x2ξ3 ∂x3ξ3 = 1 0 0 0 ∂τx1 ∂ξ1x1 ∂ξ2x1 ∂ξ3x1 ∂τx2 ∂ξ1x2 ∂ξ2x2 ∂ξ3x2 ∂τx3 ∂ξ1x3 ∂ξ2x3 ∂ξ3x3 −1 = 1 J J 0 0 0 J01 J11 J21 J31 J02 J12 J22 J32 J03 J13 J23 J33 ,
where J = |∂(x1, x2, x3)/∂(ξ1, ξ2, ξ3)| is the determinant of the matrix ∂(x1, x2, x3)/∂(ξ1, ξ2, ξ3), and Jij
for i, j = 1, 2, 3, are the entries satisfying the following expressions:
J11= |∂(x2, x3)/∂(ξ2, ξ3)| , J21= |∂(x1, x3)/∂(ξ3, ξ2)| , J31= |∂(x1, x2)/∂(ξ2, ξ3)| ,
J12= |∂(x2, x3)/∂(ξ3, ξ1)| , J22= |∂(x1, x3)/∂(ξ1, ξ3)| , J32= |∂(x1, x2)/∂(ξ3, ξ1)| ,
J13= |∂(x2, x3)/∂(ξ1, ξ2)| , J23= |∂(x1, x3)/∂(ξ2, ξ1)| , J33= |∂(x1, x2)/∂(ξ1, ξ2)| .
In addition to that, we also have J0j= −PNi=1d Jij∂τxi, j = 1, 2, 3. Now the metrics in (7) which consist
of twelve spatial- and temporal-varying quantities in a three-dimensional logical domain may take the form
Dξj= ∂tξj, ∇X~ξj = (∂tξj, ∂x1ξj, ∂x2ξj, ∂x3ξj) =
1
J (J0j, J1j, J2j, J3j) (10) for j = 1, 2, 3, that are conveniently dependent on the derivatives in the physical domain.
With these definitions, after some simple algebraic manipulations (cf. [15]), the transformed version of our two-fluid model in the logical domain ˆΩ can be written as
∂ ∂τ (ρJ) + Nd X j=1 ∂ ∂ξj (ρJUj) = 0, ∂ ∂τ (ρJui) + Nd X j=1 ∂ ∂ξj J ρuiUj+ p∂ξj ∂xi = 0 for i = 1, 2, . . . , Nd, ∂ ∂τ (JE) + Nd X j=1 ∂ ∂ξj J EUj+ pUj− p ∂ξj ∂t = 0, ∂ ∂τ ρB γ − 1J + Nd X j=1 ∂ ∂ξj ρB γ − 1JUj = 0, ∂α ∂τ + Nd X j=1 Uj∂α ∂ξj = 0, (11)
where Uj= ∂tξj+PNi=1d ui∂xiξjdenotes the contravariant velocity in the ξj-direction for j = 1, 2, . . . , Nd.
It is clear that in this system the first Nd+ 2 are simply the compressible Euler equations in strong
conservation-law form, while the remaining are the effective equations that are introduced to ensure the correct mixing of the problem-dependent material variables near the interfaces. It is easy to observe also that the flux functions of the conservation laws in the transformed model depend not only on the conserved variables but also on the spatially-varying metrics which is unlike the case in the original model where the fluxes are dependent on the conserved variables only. We should have this in mind when we want to develop a method that discretizes the model for numerical approximation.
It is important to note that in the general ∂ ~X/∂t 6= 0 case of (5), to complete the model, we have to devise a way for the numerical resolution of the metrics that are presence in (11) throughout the computations. Among a wide variety of methodologies proposed in the literature (cf. [1, 9, 21, 30, 32, 53]),
in this work, we are interested in an approach that is based on the compatibility condition of the mixed second order partial derivatives ∂τ∂ξjxiand ∂ξj∂τxito yield the geometric conservation laws of the form
∂ ∂τ ∂xi ∂ξj + ∂ ∂ξj −∂xi ∂τ = 0 for i, j = 1, 2, . . . , Nd.
This approach has been used quite successfully in a unified coordinate method (cf.[16, 17, 19]) for solving general fluid flow problems with stationary or moving complex geometries and also in a Lagrangian scheme for gas dynamics (cf. [9]). It is a simple way to be coped with our current implementation of a dimension-by-dimension splitting method described in Section 4 for numerical approximation of multi-dimensional problems. Note that, to avoid any possible confusion that may cause in the above equations, we introduce a new set of notation χji = ∂ξjxi as for the unknowns (N
2
d of them in total for i, j = 1, 2, . . . , Nd) and
when combining the mesh-moving relation (5) for ∂τxi in the problem formulation, we therefore have
∂χji
∂τ + ∂ ∂ξj
(−hui) = 0 for i, j = 1, 2, . . . , Nd. (12)
It is certain that, with a prescribed region of the physical domain for the computation (rectangular or not), the initial grid system and also the initial condition of the above equations can be obtained by a chosen technique from numerical grid generators (cf. [4, 54]) when it is necessary. Once the solutions of (12) are known during a time step, it is an easy manner to compute the various determinants that are required in (10) for the metrics.
In summary, together with the stiffened gas equation of state (2) and the Lagrange-like moving grid condition (5), the model equations we propose to solve compressible homogeneous two-fluid problems consist of the physical part (11) and the geometrical part (12); this gives us totally N2
d+ Nd+ 4 equations
to be solved in an Ndspatial (logical) domain. Surely, when h approaches to zero, this model reduces to
the physical part (11) only. With a model system formulated in this way, there is no difficulty to compute the pressure from the equation of state (4) in the entire computational region at all time.
For the ease of the later reference, this two-fluid model is written in the following form
∂q ∂τ + Nd X j=1 fj ∂ ∂ξj , q = 0, (13)
where in the two-dimensional case Nd= 2, for example, the state vector q and the vector-value functions
f1, f2are defined by
q =
ρJ, ρJu1, ρJu2, EJ, ρB
γ − 1J, α, χ11, χ12, χ21, χ22 T , (14a) f1= ∂ ∂ξ1(ρJU1) , ∂ ∂ξ1(ρJu1U1+ pχ22) , ∂ ∂ξ1(ρJu2U1− pχ21) , ∂ ∂ξ1(EJU1+ p(u1χ22− u2χ21)), ∂ ∂ξ1 ρB γ − 1JU1 , U1∂α ∂ξ1 , ∂ ∂ξ1 (−hu1) , ∂ ∂ξ1 (−hu2) , 0, 0 T , (14b) and f2= ∂ ∂ξ2 (ρJU2) , ∂ ∂ξ2 (ρJu1U2− pχ12) , ∂ ∂ξ2 (ρJu2U2+ pχ11) , ∂ ∂ξ2 (EJU2+ p(u2χ11− u1χ12)), ∂ ∂ξ2 ρB γ − 1JU2 , U2∂α ∂ξ2 , 0, 0, ∂ ∂ξ2 (−hu1) , ∂ ∂ξ2 (−hu2) T , (14c)
in a respective manner, with
U1= (1 − h)(u1χ22− u2χ21)/J,
U2= (1 − h)(u2χ11− u1χ12)/J,
(15)
and J = χ11χ22− χ21χ22. Note that when h 6= 1, the proposed system is not written in the full
conservation form, but is rather a quasi-conservative system of equations. In the special Lagrangian case where h = 1, however, we have a conservative system where (14b) is simplified to
f1= 0, ∂pχ22 ∂ξ1 , −∂pχ21 ∂ξ1 , ∂ ∂ξ1 p(u1χ22− u2χ21), 0, 0, − ∂u1 ∂ξ1 , −∂u2 ∂ξ1 , 0, 0 T , (16a) and (14c) is simplified to f2= 0, −∂pχ12 ∂ξ2 , ∂pχ11 ∂ξ2 , ∂ ∂ξ2 p(u2χ11− u1χ12), 0, 0, 0, 0, − ∂u1 ∂ξ2 , −∂u2 ∂ξ2 T , (16b)
see [9] for an alternative and yet interesting form of the Lagrangian system when the additional con-straints: ∂χ21 ∂ξ1 = ∂χ11 ∂ξ2 and ∂χ22 ∂ξ1 =∂χ12 ∂ξ2 ,
are being imposed in the mathematical formulation.
In addition, we use the notations ˘fι to denote the flux functions for the conservation laws of (13) in
the ξι-direction, yielding
˘ f1=
ρJU1, ρJu1U1+ pχ22, ρJu2U1− pχ21, EJU1+ p(u1χ22− u2χ21),
ρB γ − 1JU1,
0, −hu1, −hu2, 0, 0
T (17a)
for the fluxes in the ξ1-direction, and
˘ f2=
ρJU2, ρJu1U2− pχ12, ρJu2U2+ pχ11, EJU2+ p(u2χ11− u1χ12),
ρB γ − 1JU2,
0, 0, 0, −hu1, −hu2
T (17b)
for the fluxes in the ξ2-direction when Nd= 2.
3
Hyperbolicity of the model equations
To examine the hyperbolicity of (13), it is a common practice to inquire the characteristic structure of the quasi-linear form of the equations
∂q ∂τ + Nd X j=1 Aj(q) ∂q ∂ξj = 0, (18)
when the proper smoothness of the solutions is assumed. Without loss of generality, we consider the two-dimensional case Nd= 2 as an example, and so have the state vector q in (18) defined by (14a) and
the matrices Aj for j = 1, 2 that are derived from (14b) and (14c), respectively, to be of the form
with Aj,m, m = 1, 2, . . . , 10, denoting the mth row vector for matrix Aj. Here, for matrix A1, we have A1,1= 0, (1 − h)χ22 J , − (1 − h)χ21
J , 0, 0, 0, −ρU1χ22, ρU1χ21, −ρU2χ22, ρU2χ21 , A1,2 = Kχ22 J − u1U1, u1(1 − h − Γ)χ22 J + U1, − u1(1 − h)χ21 J − u2Γχ22 J , Γχ22 J , Γχ22 J , ϕχ22, − ρu1U1+ (p + ρ0B)χ22 J χ22, ρu1U1+ (p + ρ0B)χ22 J χ21, − ρu1U2− (p + ρ0B)χ12 J χ22, ρu1U2− (p + ρ0B)χ12 J χ21− ρ0B , A1,3 = −Kχ21 J − u2U1, u2(1 − h)χ22 J + u1Γχ21 J , − u2(1 − h − Γ)χ21 J + U1, − Γχ21 J , − Γχ21 J , −ϕχ21, − ρu2U1−(p + ρ0B)χ21 J χ22, ρu2U1−(p + ρ0B)χ21 J χ21, − ρu2U2+(p + ρ0B)χ11 J χ22+ ρ0B, ρu2U2+(p + ρ0B)χ11 J χ21 , A1,4 = (K − M )V1, M χ22 J − u1V1Γ, − M χ21 J − u2V1Γ, (1 − h + Γ)V1, ΓV1, ϕV1, −(ρM + ρ0B)V1χ22, (ρM + ρ0B)V1χ21, −(ρM + ρ0B)V2χ22+ ρ0Bu2, (ρM + ρ0B)V2χ21− ρ0Bu1 , A1,5 = −BU1 Γ , B(1 − h)χ22 ΓJ , − B(1 − h)χ21 ΓJ , 0, U1, 0, − ρBU1χ22 Γ , ρBU1χ21 Γ , −ρBU2χ22 Γ , ρBU2χ21 Γ , A1,6 = [0, 0, 0, 0, 0, U1, 0, 0, 0, 0] , A1,7 = u1h ρJ , − h ρJ, 0, 0, 0, 0, 0, 0, 0, 0 , A1,8 = u2h ρJ , 0, − h ρJ, 0, 0, 0, 0, 0, 0, 0 , A1,9 = 0, A1,10 = 0,
and for matrix A2, we have
A2,1 =
0, −(1 − h)χ12 J ,
(1 − h)χ11
J , 0, 0, 0, ρU1χ12, −ρU1χ11, ρU2χ12, −ρU2χ11 , A2,2 = −Kχ12 J − u1U2, − u1(1 − h − Γ)χ12 J + U2, u1(1 − h)χ11 J + u2Γχ12 J , − Γχ12 J , − Γχ12 J , −ϕχ12, ρu1U1+(p + ρ0B)χ22 J χ12, − ρu1U1+(p + ρ0B)χ22 J χ11+ ρ0B, ρu1U2− (p + ρ0B)χ12 J χ12, − ρu1U2− (p + ρ0B)χ12 J χ11 ,
A2,3= Kχ11 J − u2U2, − u2(1 − h)χ12 J − u1Γχ11 J , u2(1 − h − Γ)χ11 J + U2, Γχ11 J , Γχ11 J , ϕχ11, ρu2U1−(p + ρ0B)χ21 J χ12− ρ0B, − ρu2U1−(p + ρ0B)χ21 J χ11, ρu2U2+(p + ρ0B)χ11 J χ12, − ρu2U2+(p + ρ0B)χ11 J χ11 , A2,4= (K − M )V2, − M χ12 J − u1V2Γ, M χ11 J − u2V2Γ, (1 − h + Γ)V2, ΓV2, ϕV2, (ρM + ρ0B)V1χ12− ρ0Bu2, −(ρM + ρ0B)V1χ11+ ρ0Bu1, (ρM + ρ0B)V2χ12, − (ρM + ρ0B)V2χ11 , A2,5= −BU2 Γ , − B(1 − h)χ12 ΓJ , B(1 − h)χ11 ΓJ , 0, U2, 0, ρBU1χ12 Γ , − ρBU1χ11 Γ , ρBU2χ12 Γ , − ρBU2χ11 Γ , A2,6= [0, 0, 0, 0, 0, U2, 0, 0, 0, 0] , A2,7= 0, A2,8= 0, A2,9= u1h ρJ , − h ρJ, 0, 0, 0, 0, 0, 0, 0, 0 , A2,10= u2h ρJ , 0, − h ρJ, 0, 0, 0, 0, 0, 0, 0 .
Note that in the above matrix entries we have used the following notations for abbreviation: Γ = γ − 1, K = Γ(u2
1+ u22)/2, M = ((1 − h)E + p)/ρ, V1 = ~u · ∇X~ξ1 = (u1χ22− u2χ21)/J, V2 = ~u · ∇X~ξ2 =
(u2χ11− u1χ12)/J, and ϕ = −Γ[(p + ρ01B1)/(γ1− 1) − (p + ρ02B2)/(γ2− 1)].
Before proceeding further, we recall that, for any given Nd, a nonlinear system of partial differential
equations of the form (13) is said to be (strongly) hyperbolic if any linear combination of the matrices Aj for j = 1, 2, . . . , Nd appearing in (18) has (i) real eigenvalues and (ii) a complete set of linearly
independent right eigenvectors. It is said to be weakly hyperbolic if (i) is satisfied but (ii) is not fulfilled due to the presence of a Jordan block (cf. [13, 23]).
Now let ~nj = ∇X~ξj/|∇X~ξj| =
n(1)j , n(2)j be the unit normal vector to the ξj-direction, and let
~tj = t(1)j , t (2) j =−n(2)j , n (1) j
be the respective unit transverse vector. Then the velocity vector ~u in the Cartesian coordinate can be projected into the local ~nj and ~tj coordinate as ~u = υj~nj+ ωj~tj, where
υj= ~u · ~nj is the component in the ~nj-direction and ωj= ~u · ~tj is the component in the ~tj-direction.
With these definitions, the eigenvalues λj of matrix Aj for j = 1, 2 can be computed in a standard
way (cf. [52]), and the results for the general moving grid case, h 6= 0, are λj,0= 0 (algebraic multiplicity 4), λj,1= Uj− c ∇X~ξj (algebraic multiplicity 1), λj,2= Uj (algebraic multiplicity 4), λj,3= Uj+ c ∇X~ξj (algebraic multiplicity 1), (19)
where c =pγ(p + p0)/ρ is the speed of sound of the fluid. Note that in this work we assume that the
thermodynamic description of the materials of interest is limited by the stability requirement that the speed of sound c belongs to a set of real numbers, and so all the eigenvalues λj that appear in (19) are
In the current instance, with the help of a symbolic algebra software such as Maple, the corresponding right eigenvectors for matrix A1 are: for λ1,1,
r1,1 = " 1, u1− n(1)1 c (−) 1 , u2− n(2)1 c (−) 1 , H − υ1c(−)1 − phβ1(−) ρ , B Γ, 0, hn(1)1 c (−) 1 ρJλ1,1 , hn (2) 1 c (−) 1 ρJλ1,1 , 0, 0 #T , for λ1,2, r1,2= 1, u1, u2, K Γ, 0, 0, 0, 0, 0, 0 T , r1,3= " 0, U1t(1)1 , U1t(2)1 , ω1U1, 0, 0, hn (2) 1 ρJ , − hn(1)1 ρJ , 0, 0 #T , r1,4= [0, 0, 0, −1, 1, 0, 0, 0, 0, 0]T, r1,5= 0, 0, 0, ϕJ Γ , 0, 1, 0 , 0, 0, 0 T , for λ1,3, r1,6= " 1, u1+ n(1)1 c (+) 1 , u2+ n(2)1 c (+) 1 , H + υ1c(+)1 + phβ1(+) ρ , B Γ, 0, − hn(1)1 c(+)1 ρJλ1,3 , − hn(2)1 c(+)1 ρJλ1,3 , 0, 0 #T , and for λ1,0, r1,7= [0, 0, 0, 0, 0, 0, χ21, χ22, 0, 0]T, r1,8= 1, u1, u2, E ρ, 0, 0, χ11 ρJ, χ12 ρJ , 0, 0 T .
Moreover, the right eigenvectors for matrix A2 are: for λ2,1,
r2,1= " 1, u1− n(1)2 c (−) 2 , u2− n(2)2 c (−) 2 , H − υ2c(−)2 − phβ2(−) ρ , B Γ, 0, 0, 0, hn(1)2 c (−) 2 ρJλ2,1 , hn(2)2 c (−) 2 ρJλ2,1 #T , for λ2,2, r2,2= 1, u1, u2, K Γ, 0, 0, 0, 0, 0, 0 T , r2,3= " 0, U2t(1)2 , U2t(2)2 , ω2U2, 0, 0, 0, 0, hn (2) 2 ρJ , − hn(1)2 ρJ #T , r2,4= [0, 0, 0, −1, 1, 0, 0, 0, 0, 0]T, r2,5= 0, 0, 0, ϕJ Γ , 0, 1, 0 , 0, 0, 0 T , for λ2,3, r2,6= " 1, u1+ n(1)2 c (+) 2 , u2+ n(2)2 c (+) 2 , H + υ2c(+)2 + phβ2(+) ρ , B Γ, 0, 0, 0, − hn(1)2 c(+)2 ρJλ2,3 , − hn(2)2 c(+)2 ρJλ2,3 #T ,
and for λ2,0, r2,7= [0, 0, 0, 0, 0, 0, 0, 0, χ11, χ12]T, r2,8= 1, u1, u2, E ρ, 0, 0, 0, 0, χ21 ρJ, χ22 ρJ T . Here we have c(±)j = c(1±hβ (±) j ), β (−) j = c|∇X~ξj|/(λj,1+ch|∇X~ξj|), and β (+) j = c|∇X~ξj|/(λj,3−ch|∇X~ξj|) for j = 1, 2.
From the above, it is easy to observe that in either case the geometric multiplicity of the zero eigenvalue λ1,0 or λ2,0 is 2, but is not 4, yielding two missing eigenvectors for matrices A1 and A2 (this is not
surprising however because there are two null vectors in each of these matrices). In addition, in case the contravariant velocity Uj in the ξj-direction is happened to be zero which would occur in particular
when h = 1 or ~u = 0 as seen from (15), we loss another eigenvector for matrix Aj because the two
previously linearly independent eigenvectors rj,3 and rj,7 becomes parallel to each other. In view of this
fact that, we conclude that our two-fluid model in the generalized Lagrangian coordinate (13) is only weakly hyperbolic in two dimensions; this can be shown to be true in any Nd space dimension, with the
exception of Nd= 1 (cf. [16, 17]). A more detailed discussion about the mathematical significance of the
loss of strong hyperbolicity for the Lagrangian gas dynamics system in two dimensions is given in [9], for example. Surely, in the Eulerian case, h = 0, where the basic motion of the flow is governed by (11), our model system can be shown to be of (strongly) hyperbolic type that satisfies all the requirement mentioned above, irrespective of the number of space dimension.
4
Wave propagation methods
To find approximate solutions of our homogeneous two-fluid flow model (13) in a generalized coordinate that contains spatially varying metrics in the flux functions of the physical conservation laws, we use a flux-based wave decomposition method developed by Bale et al. [2, 40] with the dimensional-splitting technique incorporated in the method for multidimensional problems. This method is a variant of the standard wave-propagation scheme [26, 48] in that we solve one-dimensional Riemann problems at each cell interface as usual, but rather than using the resulting jumps of the state variables moving at constant speeds to update the solutions in neighboring grid cells, we use the resulting jumps of fluxes for that instead, which has shown (cf. [2]) to be a much more accurate approach than the former one for a class of hyperbolic problems with spatially varying fluxes (this is the type of problems concerned here).
To explain the basic idea of the method, we consider the two-dimensional case Nd= 2 as an example,
and so in a simple dimensional-splitting approach the equations to be solved ∂q ∂τ + f1 ∂ ∂ξ1 , q + f2 ∂ ∂ξ2 , q = 0
with q, f1 and f2 defined by (14a), (14b), and (14c), respectively, are split into a sequence of
one-dimensional problems as ξ1-sweeps: ∂q ∂τ + f1 ∂ ∂ξ1 , q = 0, (20a) ξ2-sweeps: ∂q ∂τ + f2 ∂ ∂ξ2 , q = 0. (20b)
Assuming a uniform rectangular grid with a fixed mesh spacing ∆ξ1 in the ξ1-direction and ∆ξ2 in the
i − 1 i − 1 i i j j j + 1 j + 1 Cij ˆ Cij ξ1 ξ2 mapping ∆ξ1 ∆ξ2 logical domain physical domain ←− x1= x1(ξ1, ξ2) x2= x2(ξ1, ξ2) x1 x2
Figure 2: A sample grid system that is formed at a time in our two-dimensional generalized Lagrangian scheme. The numerical solution on the rectangular grid cell Cij in the logical domain gives distinctively
the result on the mapped quadrilateral grid cell ˆCij in the physical domain for all the grid cell (i, j).
standard finite-volume formulation in which the approximate value Qn
ij of the cell average of the solution
over the (i, j)th grid cell at time τn can be written as
Qnij≈ 1 ∆ξ1∆ξ2 Z Cij q(ξ1, ξ2, τn) dξ1dξ2,
where Cij denotes the rectangular region occupied by the grid cell (i, j). Note that in the algorithm the
numerical solution on the rectangular grid cell Cij in the logical domain gives distinctively the result on
the mapped quadrilateral grid cell ˆCij in the physical domain for all the grid cell (i, j), see Fig. 2. The
time step from the current time τn to the next τn+1is denoted by ∆τ .
In this numerical setup, a dimensional-splitting (or called Godunov-splitting) version of the first-order wave-propagation method in two dimensions can be written as
Q∗ij = Qnij− ∆τ ∆ξ1 h A+1∆Q n i−1/2,j+ A − 1∆Q n i+1/2,j i , (21a) Qn+1ij = Q∗ ij− ∆τ ∆ξ2 h A+2∆Q ∗ i,j−1/2+ A − 2∆Q ∗ i,j+1/2 i . (21b)
Here in the ξ1-sweeps we start with cell average Qnij at time τn and solve (20a) along each row of cells
Cij with j fixed, updating Qnij to Q∗ij by the use of (21a) with the fluctuations
(A+1∆Q)ni−1/2,j = X m:(λ1,m)ni−1/2,j>0 (Z1,m)ni−1/2,j and (A−1∆Q)ni+1/2,j = X m:(λ1,m)ni+1/2,j<0 (Z1,m)ni+1/2,j,
where (λ1,m)nι−1/2,j and (Z1,m)nι−1/2,j are in turn the wave speed and the f -waves (the flux-based wave
decomposition) for the mth family of the solutions obtained from solving the one-dimensional Riemann problems in the direction normal to the cell interface between Cι−1,j and Cιj with Qnι−1,j and Qnιj as
initial data, see Section 5 for the details. Then in the ξ2-sweeps we can use the Q∗ij values as data for
solving (20b) along each column of cells Cijwith i fixed, which gives us the solution of the next time step
Qn+1ij from (21b) with the fluctuations
(A+2∆Q) ∗ i,j−1/2= X m:(λ2,m)∗i,j−1/2>0 (Z2,m)∗i,j−1/2
and (A− 2∆Q) ∗ i,j+1/2 = X m:(λ2,m)∗i,j+1/2<0 (Z2,m)∗i,j+1/2.
It is clear that this method belongs to a class of upwind schemes, and is stable when the typical CFL (Courant-Friedrichs-Lewy) condition:
ν = ∆τ maxm(λ1,m, λ2,m) min (∆ξ1, ∆ξ2)
≤ 1, (22)
is satisfied (cf. [11, 25, 26]). Moreover, it is not difficult to show that the method is quasi-conservative in the sense that when applying the method to (13) not only the conservation laws but also the transport equations are approximated in a consistent manner by the method.
To extend this splitting method to a high-resolution version (i.e., second-order accurate on smooth solutions, and sharp and monotone profiles on discontinuous solutions), it is a common practice to modify (21a) and (21b), in a respective manner, as
Q∗ ij= Qnij− ∆τ ∆ξ1 h A+1∆Q n i−1/2,j+ A − 1∆Q n i+1/2,j i − ∆τ ∆ξ1 ˜F1n i+1/2,j− ˜F1n i−1/2,j and Qn+1ij = Q ∗ ij− ∆τ ∆ξ2 h A+2∆Q ∗ i,j−1/2+ A − 2∆Q ∗ i,j+1/2 i − ∆τ ∆ξ2 ˜F2∗ i,j+1/2− ˜F2∗ i,j−1/2 ,
where the add in correction terms such as ( ˜F1)ni−1/2,j and ( ˜F2)∗i,j−1/2, for example, written in view of the
f -waves approach with mw wave family in total may have the form
( ˜F1)ni−1/2,j= 1 2 mw X m=1 sign (λ1,m) 1 − ∆τ ∆ξ1 |λ1,m| ˜ Z1,m n i−1/2,j and ( ˜F2)∗i,j−1/2= 1 2 mw X m=1 sign (λ2,m) 1 − ∆τ ∆ξ2 |λ2,m| ˜ Z2,m ∗ i,j−1/2 ,
see [2, 26] for more expositions. Note that ˜Zι,m is a limited value of Zι,mobtained by comparing Zι,m
with the corresponding Zι,mfrom the neighboring Riemann problem to the left (if λι,m > 0) or to the
right (if λι,m< 0) for ι = 1, 2.
To end this section, we want to mention that the extension of this wave propagation method from two to three space dimensions can be done in a straightforward manner when the solutions of one-dimensional Riemann problems in each dimensional-sweep can be computed readily (cf. [17]).
5
Riemann problem and approximate solutions
To determine the aforementioned wave speeds and fluctuations in the wave-propagation methods, we need to solve one-dimensional Riemann problems in direction normal to each cell interfaces. If we consider the case between cells Ci−1,jand Cijas illustrated in Fig. 2, for example, it would be a standard one when we
define the normal Riemann problem to this face as a Cauchy problem that consists of the equations (20a) with the piecewise constant initial data
q(ξ1, 0) = Qn i−1,j if ξ1< (ξ1)i−1/2, Qn ij if ξ1> (ξ1)i−1/2. (23)
In this instance, from (20a) it is easy to observe that all the flux functions for the conservation laws,i.e., ˘
f1 given in (17a), of the Riemann problem bears some relations to the time-independent part of the
geometric variables χ21 = ∂ξ2x1 and χ22 = ∂ξ2x2, while there is no relation for any of them to the
remaining time-dependent geometric variables χ11 = ∂ξ1x1 and χ12 = ∂ξ1x2. With this in mind, it
should be sensible (cf. [2, 17, 18, 19, 36, 37]) to put the original Riemann problem into a new generalized form as ∂qi−1,j ∂τ + f1 ∂ ∂ξ1 , qi−1,j = 0 if ξ1< (ξ1)i−1/2, ∂qij ∂τ + f1 ∂ ∂ξ1 , qij = 0 if ξ1> (ξ1)i−1/2, (24)
for the equations and (23) for the initial data as before. Here in the above we have used a cell-average approach where q and f1(∂ξ1, q) are discretized to yield a new function qij = q|(χ21,χ22)ij and f1(∂ξ1, qij) =
f1(∂ξ1, q)|(χ21,χ22)ij, respectively, that holds throughout the grid cell Cij.
We note that as to the normal-direction (~n1)i−1/2,j = (∇X~ξ1/|∇X~ξ1|)i−1/2,j to the cell interface
(ξ1)i−1/2 which should be a constant vector with unit length in this Riemann problem, we use the
relation
(~n1)i−1/2,j = ( ¯χ22, − ¯χ21) / ¯S (25)
for that, where ¯z = (zi−1,j+ zij)/2 is the average value of zi−1,j and zij for z = χ21, χ22, and ¯S =
p ¯χ2
21+ ¯χ222 is the Euclidean distance of the vector ( ¯χ22, − ¯χ21) in two dimensions. With that, the
transverse-direction to this cell interface is taken to be (~t1)i−1/2,j= ( ¯χ21, ¯χ22)/ ¯S.
Recall that, for the single component case with the ideal gas law, it has been demonstrated by Hui and coworkers [17, 19] that with a reasonable set of the initial data solutions of this generalized Riemann problem exist that consist of the classical self-similar solutions such as shock waves, rarefaction waves, and contact discontinuities, for gas dynamics, and an additional stationary discontinuity at the cell interface (to be called a cell-edge discontinuity below) where there are only jumps on the geometric variables but are uniform on the physical flow variables across it, see Fig. 3 for an illustration of the basic solution structure. As a consequence, in their work, they are able to take the usual procedures as if there is no such a cell-edge discontinuity to the construction of the exact solution for the Riemann problem, see [18] also for more discussion when the shallow water equations are involved in the problem formulation.
Due to the possible loss of strong hyperbolicity of the proposed model, see Section 3, which would lead to numerical difficulties in using a linearized Riemann solver based on eigenvector decomposition of state variables (cf. [38, 55]), in the current work with the homogeneous two-fluid model (13) and the stiffened gas equation of state (2), we are interested in a popular shock-only (or called two-shock) approach (cf. [6]) that ignores the possibility of rarefaction waves and simply build a solution in which each pair of the states is connected along the Hugoniot locus for a shock. To simplify the notation in the following discussion, we denote zL= zi−1,jn and zR= zijn as being the Riemann data for a state variable
z on the left and right of the cell interface, respectively. We then write the rotated velocity components of the Riemann data in the normal- and transverse-direction to the cell interface as υι = ~uι· (~n1)i−1/2,j
and ωι = ~uι· ~t1i−1/2,j, in an respective manner, for ι = L or R.
As in the standard Riemann problem on general quadrilateral grids (cf. [3, 26]), the key step in this shock-only approach is to find the midstate (υm, pm) in the υ-p phase plane so that it can connect to
(υL, pL) by a 1-shock, and to (υR, pR) by a 3-shock. For that, this amounts to solving the following
nonlinear equation in an iterative manner for the pressure pm:
ψ(pm) = υmR(pm) − υmL(pm) = 0. (26)
ξ1 τ Qn i−1,j Qnij 3-shock contact 1-rarefaction 1-shock cell-edge
∂τqi−1,j+ f1(∂ξ1, qi−1,j) = 0 (ξ1)i−12 ∂τqij+ f1(∂ξ1, qij) = 0
qmL− q+mL qmR
Figure 3: Typical solution structure of the generalized Riemann problem for our homogeneous two-fluid model discussed in Section 5. The key step in obtaining this solution is to find the pressure and the velocity (the normal component of it) in the regions mL and mR where both of them are continuous across the
contact and the cell-edge discontinuities. Note that there are only jumps in the geometric variables across the cell-edge discontinuity that separates the constant states q−
mLand q +
mL. In a shock-only approximate
Riemann solver, the rarefaction wave is replaced by an entropy-violating shock.
curves, respectively,
υmL(p) = υL−p − pL
ML(p), υmR(p) = υR+
p − pR
MR(p), (27)
with Mι denoting the Lagrangian shock speed, for ι = L or R. In the current application of the pressure
law (2), we may compute Mι quite easily by evaluating the formula
Mι2(p) = Cι2 1 + γι+ 1 2γι p + p0ι pι+ p0ι − 1 (28)
where Cι= ριcι is the Lagrangian sound speed, and ρmι is the midstate density on the ι side,
ρmι(p) = ρ−1ι − p − pι M2 ι(p) −1 (29)
Note that (28) and (29) are as a result derived from the Rankine-Hugoniot jump conditions across the shock waves (cf. [8, 17, 19]).
When applying a standard root-finding approach such as the secant method to (26), we have a 2-step iteration scheme as follows,
p(k+1)m = p(k)m − p (k) m − p(k−1)m υ (k) mL− υ (k−1) mL + υ (k) mR− υ (k−1) mR h υ(k)mR− υmL(k)i, (30) where υmι(k) = υmι h p(k)m i
, for ι = L or R, and k = 1, 2, . . . (until convergence). With a suitable choice of the starting values p(0)m and p(1)m, method (30) typically converges to the exact solution pm at a super
linear rate [22]. For gas dynamics, it is a common practice to set p(0)m and p(1)m by
p(0)m = pRCL+ pLCR− (υR− υL)CLCR CL+ CR , p(1) m = pRML(0)+ pLMR(0)− (υR− υL)ML(0)M (0) R ML(0)+ MR(0) , (31)
where Mι(0)= Mι
h p(0)m
i
. Having that, we may assign υmL(0) and υ(0)mRby
υ(0)mL= υL− p(0)m − pL CL , υmR(0) = υR+ p(0)m − pR CR ,
and υmL(1) and υmR(1) according to (27). After a satisfactory convergence of the scheme, υm can then be
calculated based on the formula:
υm=
pL− pR+ υLML(pm) + υRMR(pm)
ML(pm) + MR(pm)
.
Figure 3 illustrates a typical solution structure of the present two-fluid generalized Riemann problem. Clearly, in a two-shock approximate solver, we replace the leftward-going 1-rarefaction wave by a 1-shock, and so the solution would consist of four discontinuities moving at constant speeds. Here the propagation speed of each discontinuity is determined by
(λ1,0)i−1/2,j = 0, (λ1,1)i−1/2,j = (1 − h)υm− ML(pm) ρmL(pm) ∇X~ξ1 i−1/2,j, (λ1,2)i−1/2,j = (1 − h)υm ∇X~ξ1 i−1/2,j, (λ1,3)i−1/2,j = (1 − h)υm+ MR(pm) ρmR(pm) ∇X~ξ1 i−1/2,j, (32) where ∇X~ξ1
i−1/2,j is a scale factor for the wave speed which can be taken as an average value based
on ∇X~ξ1 Land ∇X~ξ1
R, for example. We note that the zero th- and the second-wave families, i.e., the
case with speeds λ1,0 and λ1,2, respectively, would correspond to the linearly degenerate field such as the
cell-edge and contact discontinuities, and the remaining wave families, i.e., the case with speeds λ1,1 and
λ1,3, would correspond to the genuinely nonlinear field such as the shock waves and rarefactions (cf. [51]).
To define the jumps across each of these discontinuities, in a spatially varying fluxes case as the one considered here, it is known (cf. [2, 40]) to have the advantage of using a variant of the flux-based wave decomposition approach in that the flux difference ˘f1(QR) − ˘f1(QL) (see (17a) for the definition of ˘f1)
is decomposed into a sum of the difference between the fluxes to the left and right of the discontinuity,
˘ f1(QR) − ˘f1(QL) = 3 X m=0 (Z1,m)i−1/2,j,
where in the wave configuration shown in Fig. 3, for example, we have
(Z1,0)i−1/2,j = ˘f1 qmL+ − ˘f1 qmL− ,
(Z1,1)i−1/2,j = ˘f1 qmL− − ˘f1(QL) ,
(Z1,2)i−1/2,j = ˘f1(qmR) − ˘f1 q+mL ,
(Z1,3)i−1/2,j = ˘f1(QR) − ˘f1(qmR) ,
(33)
with a modification on the sixth-component of the second-wave family (i.e., for the contact discontinuity) as
Z1,2(6)
i−1/2,j = (λ1,2)i−1/2,j(αR− αL)
that takes account of the effect due to the non-conservative transport equation for the volume fraction when h 6= 1. Note that in (33) ˘f1 qmL− is calculated using the data: ρmL, υm, ωL, pm, γL, BL, p0L,
(χ21)L, (χ22)L, ˘f1 qmL+ is calculated using the data: ρmL, υm, ωL, pm, γL, BL, p0L, (χ21)R, (χ22)R, and
˘
f1(qmR) is calculated using the data: ρmR, υm, ωR, pm, γR, BR, p0R, (χ21)R, (χ22)R. As usual, wave
propagation methods are based on using these propagating discontinuities to update the cell averages in the cells neighboring to each interface, see Section 4.
As a final remark, when a more general equation of state (cf. [7, 46]) is included in the problem formulation, due to the complexity that is involved to find the midstate (υm, pm), this shock-only solver
may not be an effective one to be used for the numerical resolution of the Riemann problem. The development of a simplier HLL-type (cf. [14, 55]) or Roe-type (cf. [36, 37, 38]) solver for the proposed generalized Riemann problem would become necessary; this is the subject to be discussed in a sequel paper elsewhere.
6
Numerical results
We now present some sample numerical results obtained using our algorithm described in Sections 4 and 5 for single- and two-component inviscid compressible flow problems in one and two space dimen-sions. Without stated otherwise, we have carried out all the tests using the Courant number ν = 0.5 defined by (22), the grid-movement parameter h = 0.99, and the minmod limiter in the high-resolution version of the finite-volume method based on f -wave formulation. The material-dependent parame-ters in the stiffened gas equation of state (2) are set to be (γ, ρ0, B) = (1.4, 1.2kg/m3, 0) and
(4.4, 103kg/m3
, 2.64 × 106(m/s)2) for the gas- and liquid-phase (i.e., for the air and water),
respec-tively. For comparison purposes, we have also included Eulerian results when h = 0 is employed in the method to the problems.
6.1
One-dimensional case
Example 6.1.1. To begin with, we consider a two-component version of the Lax Riemann problem where a variant of this problem in the single component case is one of the popular tests for numerical validation of a Lagrangian-type scheme (cf. [9, 43]). Initially, on the left when x1∈ [0, 0.5), we have the
perfect gas with the state variables
(ρ, u1, p, γ, B, α)L= (0.445, 0.698, 3.528, 1.4, 0, 1) ,
while on the right when x1∈ [0.5, 1], we have a different perfect gas with the state variables
(ρ, u1, p, γ, B, α)R= (0.5, 0, 0.571, 1.2, 0, 0) .
The exact solution of this problem consists of a leftward-going rarefaction wave, a rightward-going contact discontinuity, and a shock wave.
Figure 4 shows numerical results for the density, velocity, and pressure at time t = 0.14 with a 100 mesh points. From the figure, it is easy to observe a sharper resolution of the contact discontinuity, when the density profile is in comparison with the one obtained using the Eulerian h = 0 version of the method. In addition, the shock and rarefaction waves are in good agreement with the exact solution as well. To demonstrate the time variation of the physical grid coordinates generated by our algorithm, in Fig. 5 we display a sample of them at eleven different times t = i × 0.014, for i = 0, 1, . . . , 10. The dynamic movement of the grid system is clearly seen. It should be noted that each little dashed lines displayed in that graph gives a cell-center location of the grid at a specific output time. Here the non-reflecting outflow boundaries are used on the left and right computational domain.
Example 6.1.2. To see how our algorithm works for complex wave interactions, we are concerned with a model single-component (gas-phase only) problem studied by Woodward and Colella [56, 57]. In
0 0.5 1 0 0.5 1 1.5 2 ρ Exact h=0 h=0.99 0 0.5 1 0 0.5 1 1.5 2 0 0.5 1 0 1 2 3 4 u1 p x1 x1 x1
Figure 4: Numerical results for a two-component Lax Riemann problem at time t = 0.14. The solid line is the exact solution, and the dotted and triangular points are the computed solutions obtained using two different grid-movement parameters h = 0 and h = 0.99, respectively, with a 100 mesh points.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.05 0.1 0.15 x1 ti m e
Figure 5: Physical grid coordinates for the generalized Lagrangian run shown in Fig. 4 at eleven different times t = i × 0.014, for i = 0, 1, . . . , 10. Each little dashed line displayed in the graph gives a cell-center location of the grid system at a specific output time.
this problem, the initial condition consists of three constant states with data
ρ u1 p L = 1 0 103 , ρ u1 p M = 1 0 10−2 , ρ u1 p R = 1 0 102 ,
where L is the state used for x1∈ [0, 0.1), M is the state used for x1∈ [0.1, 0.9), and R is the state used
for x1∈ [0.9, 1]. Here there are two solid walls at the left and right boundaries.
In this setup, it is known that, after breaking the membranes at x = 0.1 and 0.9, a shock wave, contact discontinuity, and rarefaction wave develop at each discontinuity individually. As time progresses, the shock waves move toward each other and then collide, yielding a new contact discontinuity from the collision. Further collisions then occur, see [27, 56] for a complete wave pattern of this problem in the space-time plane.
To solve this problem numerically, we use a 200 mesh points, and perform the computation in the same manner as in Example 6.1.1. In Fig. 6, the density, velocity, and pressure are presented at three different times t = 0.016, 0.032, and 0.038. From there, we again observe a sensible improvement of the contact discontinuities, especially, the emergent one after the head-on collision of the initial shock waves, when they are in comparison with the Eulerian results. Besides, we find good resolution of the
shock wave and rarefaction as we compare our numerical results with the fine grid solutions obtained using the Eulerian grid approach with a 2000 mesh points. The physical grid coordinates generated by our generalized Lagrangian algorithm are plotted in Fig. 7 at eleven different times t = i × 0.0038, for 0, 1, . . . , 10, noticing interesting numerical temporal behavior of the grid movement that follows closely to the main feature of the underlying flow field.
6.2
Two-dimensional case
Example 6.2.1. As a first example, we consider a single-component two-dimensional Riemann problem, i.e., Configuration 4 studied by Schulz-Rinne et al. [42]. In this case, the initial condition is composed of four shock waves with the data in the four quadrants given by
(ρ, u1, u2, p)1= (1.1, 0, 0, 1.1) , (ρ, u1, u2, p)2= (0.5065, 0.8939, 0, 0.35) ,
(ρ, u1, u2, p)3= (1.1, 0.8939, 0.8939, 1.1) , (ρ, u1, u2, p)4= (0.5065, 0, 0.8939, 0.35) .
Here, the fluid component is a perfect gas throughout the whole unit square domain, and the boundary conditions are non-reflecting on all sides. Note that this problem has been used as a benchmark test to verify a cell-by-cell adaptive mesh Lagrangian scheme proposed by Morrell [35].
Numerical results for a sample run of this problem with two different grid-movement parameters h = 0 and 0.99 are shown in Fig. 8, where contour plots of the density and pressure, and the physical grid system are presented at time t = 0.2 with a 200 × 200 grid. From the figure, it is interesting to see that the collisions between the initial shock waves creates an oval shape region bounded by the incident and reflected shock waves. When we make a comparison of our results with the one appeared in Fig. 6 of [42], for instance, (which was done using a state-of-the-art second order Eulerian method with a 400 × 400 grids), we find notably a better resolution of our generalized result than the Eulerian one for the slip lines that are situated in the shock-waves bounded oval region, and similar solution behaviors for the location and structure of the reflected shock waves. To give a quantitative assessment of these solutions, Fig. 9 plots the cross section of the results for the same run along the diagonal axis of the computational grid, observing good agreement of the density and pressure profiles in most places, and some deviation in the density near the center of the oval region where the slip lines from the upper-left and lower-right corners collide with each other. We note that at those points of discrepancy the topology of the physical grid cells in the generalized Lagrangian case differ a lot from the Eulerian case, see Fig. 8, and so the finite-volume cell averages displayed there do not give exactly the same point-wise physical-state correspondence. Recall that J is the Jacobian of the grid metrics.
Example 6.2.2. We are next concerned with a two-component (gas-liquid) radially symmetric prob-lem that the computed solution in two dimensions can be compared to the one-dimensional results for numerical validation. We use the following set of data for experiments in which, inside a circle of radius r0= 0.2m, the fluid is a gas-phase with
(ρ, u1, u2, p, α)r≤r0 =
1250 kg/m3, 0, 0, 109Pa, 1,
while outside the circle, the fluid is a liquid-phase with
(ρ, u1, u2, p, α)r>r0 =
103kg/m3, 0, 0, 105Pa, 0,
where r =p(x1− x01)2+ (x2− x02)2 and (x01, x02) = (0, 0). We note that due to the pressure difference
between the fluids at r = r0, breaking of the membrane occurs instantaneously, yielding an outward-going
shock wave in liquid, an inward-going rarefaction wave in gas, and a outward-going contact discontinuity lying in between that separates the gas and liquid. We note that because of the symmetry of the solution,
0 0.5 1 0 2 4 6 8 ρ 0 0.5 1 −10 0 10 20 0 0.5 1 0 100 200 300 400 u1 p t = 0.016 t = 0.016 t = 0.016 0 0.5 1 0 5 10 15 20 ρ Fine grid h=0 h=0.99 0 0.5 1 −5 0 5 10 15 0 0.5 1 0 200 400 600 u1 p t = 0.032 t = 0.032 t = 0.032 0 0.5 1 0 2 4 6 8 ρ 0 0.5 1 −5 0 5 10 15 0 0.5 1 0 200 400 600 u1 p x1 x1 x1 t = 0.038 t = 0.038 t = 0.038
Figure 6: Numerical results for the Woodward-Colella problem at three different times t = 0.016, 0.032, and 0.038. Here the solid line is the fine grid solution obtained using the Eulerian grid approach with a 2000 mesh points, and the dotted and triangular points are the computed solutions obtained using h = 0 and h = 0.99, respectively, with a 200 mesh points.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.01 0.02 0.03 0.04 x1 ti m e
Figure 7: Physical grid coordinates for the generalized Lagrangian run shown in Fig. 6 at eleven different times t = i × 0.0038, for 0, 1, . . . , 10. Here the grid system is displayed in the same manner as in Fig. 5.
0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8
Density Pressure Physical grid a) h = 0.99 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8
Density Pressure Physical grid b) h = 0
Figure 8: Numerical results for a two-dimensional Riemann problem, a 4-shocks initial condition case. Contour plots of the density and pressure (30 lines range from 0.033 to 1.998 and 0.044 to 2.685, respec-tively), and the physical grid coordinates are shown at time t = 0.2 obtained using both h = 0 and 0.99 with a 200 × 200 grid. Note that the computed velocity field is superimposed into the pressure contours. For clarity, solution coarsening factors of 10 and 4 in each x1- and x2-direction are used to graph the
velocity vectors and the physical grid system, respectively.
0 0.5 1 0 0.5 1 1 1.2 1.4 1.6 1.8 0 0.5 1 0 0.5 1 1 1.2 1.4 1.6 1.8 2 0 0.5 1 0 0.5 1 0.6 0.7 0.8 0.9 1 1.1 x1 x1 x1 x2 x2 x2 ρ p J
Figure 9: Cross-sectional plots of the results for the shown shown in Fig. 8 along the diagonal axis of the computational grid. The dotted and triangular points are the numerical solutions obtained using h = 0 and 0.99, respectively.
for simplicity, we only take a quarter of the unit square domain, i.e., (x1, x2) ∈ [0, 0.5] × [0, 0.5] m2, and
apply the line of symmetry boundary conditions to the bottom and the left sides during the computations. Fig. 10 shows contours of the density and pressure, and the physical grid system at time at time t = 120µs with a 100 × 100 grid, where the graphs are displayed in the same manner as in Fig. 8. Good resolution of the wave pattern (i.e., both the shock and interface remain circular and appear to be very well located) are easily seen from the contour plots. The scatter plots shown in Fig. 11 provide the validation of our two-dimensional results as in comparison with the “true” solution obtained from solving the one-dimensional model with appropriate source terms for the radial symmetry, using the Eulerian high-resolution method with a 1200 mesh points (cf. [44]). Here ¯u =pu2
1+ u22denotes the radial velocity.
From the figure, it is clear that our results agree quite well with the ”true” physical solutions, and are free of spurious fluctuations in the pressure near the gas-liquid interface. We also observe some improvement in the density across the interface when the current generalized Lagrangian grid approach is in use as opposed to the standard Eulerian grid approach.
0 0.2 0.4 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0 0.1 0.2 0.3 0.4 0.5
Density Pressure Physical grid
air water a) h = 0.99 0 0.2 0.4 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0 0.1 0.2 0.3 0.4 0.5
Density Pressure Physical grid
air
water b) h = 0
Figure 10: Numerical results for a two-component radially symmetric problem at time t = 120µs with a 100 × 100 grid. Here the solutions are plotted in the same manner as in Fig. 8 with the contours ranging from 20.8 to 1270.8 kg/m3 and 16.7 to 1016MPa, for the density and pressure, respectively.
Example 6.2.3. As an example to show how our algorithm works on interaction between shock wave and material interface, we are interested in a model underwater explosion problem (cf. [12, 47, 50]). In this test, we take a rectangular domain (x1, x2) ∈ [−2, 2] × [−1.5, 1]m2, and consider the initial condition
that is composed of a stationary horizontal air-water interface at the x2= 0 axis and a circular gas bubble
in water with the center (x0
0 0.1 0.2 0.3 0.4 0.5 0.2 0.4 0.6 0.8 1 1.2 1−d h=0 h=0.99 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 1.5 0 0.1 0.2 0.3 0.4 0.5 0.5 1 1.5 2 r(m) r(m) ρ (M g / m 3 ) ¯u (k m / s) p (G P a ) J air water
Figure 11: Scatter plots of the results for the run shown in Fig. 10. The solid line is the “true” solution obtained from solving the one-dimensional model with appropriate source terms for the radial symmetry using the high-resolution method and 1200 mesh points. The dotted and triangular points are the two-dimensional results obtained using h = 0 and 0.99, respectively. The dashed line is the approximate location of the air-water interface.
interface, the fluid is a perfect gas at the standard atmospheric condition,
(ρ, u1, u2, p, α) =
1.2 kg/m3, 0, 0, 105Pa, 1.
Below the air-water interface, in region inside the gas bubble the fluid is modeled as a perfect gas also with the state variables
(ρ, u1, u2, p, α) =
1250 kg/m3, 0, 0, 109 Pa, 1,
and in region outside the gas bubble the fluid is water with the state variables
(ρ, u1, u2, p, α) =
103 kg/m3, 0, 0, 105 Pa, 0.
There are three solid walls at the left, right, and bottom boundaries in the current problem formulation. As in Example 6.2.2, breaking of this underwater bubble would results in an outward-going shock wave in water, an inward-going rarefaction wave in gas, and a material interface lying in between that separates the gas and water. Soon after this shock wave is diffracted through the nearby flat air-water interface, it is known in the literature (cf. [12]) that the topology of the underwater bubble will undergo a change from the original circular-shape to an oval-like shape. As time evolves, this gas bubble would continue rising upward, causing the subsequent deformation of the horizontal air-water interface.
Contour plots of the density and pressure at four different times t = 0.2, 0.4, 0.8, and 1.2ms are presented in Fig. 12, where we have performed the computation using h = 0 and 0.95 with a 400 × 250 grid. From the density plot, we clearly observe the improvement on the use of the generalized Lagrangian method over the Eulerian method to the sharpness of the solution structure near the interfaces. Moreover, from the pressure plot, we see the smooth variation of the solution near the interface, without introducing any spurious oscillations. In Fig. 13, we present the physical grid system for the generalized Lagrangian
run shown in Fig. 12 where the dynamical movement of the grids that follow closely to the main feature of the flow is observed. The cross-section of the density and pressure for the same run along line x1= 0
is drawn in Fig. 14, giving some quantitative information about the differences of the Eulerian and generalized Lagrangian results at the selected times. The time history of the computed grid Jacobian J for the generalized Lagrangian run is plotted in Fig. 15, where Jave, Jmin, and Jmax denote the average, minimum, and maximum values of J, respectively.
To end this section, it should be mentioned that all the simulations done here were run on an Al-phaServer DS20E under the Tru64 Unix operating system. As the whole, it is observed that, the CPU time that requires for carrying out the generalized Lagrangian computations depends strongly on the uni-formity of the grid structure. Consider the case for the underwater explosion problem in Example 6.2.3, for example, it took in turn 695.2253 and 1528.683 seconds for the Eulerian and the Lagrangian run, while for the radially symmetric problem in Example 6.2.2 it took 14.6341 and 17.7257 seconds, respectively. We note that the computer programs (written in Fortran 77) used to obtain the results above are available from the author.
7
Conclusion
We have presented a simple moving grid approach for the numerical simulation of compressible homoge-neous two-fluid flow problems with a stiffened gas equation of state in more than one space dimension. The algorithm uses a Lagrangian-like condition for the temporal evolution of the underlying grid system, and employs a volume-fraction based model equations in generalized curvilinear coordinates as a basis for the principal motion of the fluid mixture. A set of geometric conservation laws is included in the mathematical formulation of the problem also for an easy computation of the grid metrics that come as a result of the coordinate transformation between the physical and logical grid. A standard high-resolution method based on the f-waves viewpoint is used to solve the proposed model equations with the dimensional-splitting technique included in the method for multidimensional problems. Numerical results presented in this paper demonstrate clearly the feasibility of the approach for a reasonable class of two-fluid problems. Ongoing work is to study more general grid-movement strategy that is required for a class of solid-solid or fluid-solid impact problems (cf. [10]) and problems with complex geometries (cf. [3]).
Acknowledgement
The author would like to thank Professors W.-H. Hui and J.-J. Hu for valuable conversations and help in the early stage of the work presented here. This research was supported in part by National Science Council of Republic of China Grants NSC 93-2115-M-002-010, 94-2115-M-002-016, and 95-2115-M-002-010.
References
[1] D. A. Anderson, J. C. Tannehill, and R. H. Pletcher. Computational Fluid Mechanics and Heat Transfer. McGraw-Hill, New York, 1984.
[2] D. Bale, R. J. LeVeque, S. Mitran, and J. A. Rossmanith. A wave propagation method for conserva-tion laws and balance laws with spatially varying flux funcconserva-tions. SIAM J. Sci. Comput., 24:955–978, 2002.
a) Density air water h = 0 h = 0.9 t = 0.2ms t = 0.4ms t = 0.8ms t = 1.2ms
Figure 12: Numerical results for the simulation of a model underwater explosion problem. Contour plots are shown for (a) the density (30 lines range from 18.6 to 1137.2 kg/m3), and (b) the pressure (30 lines range from 6.43 to 392.18MPa) at four different times t = 0.2, 0.4, 0.8, and 1.2ms obtained using both h = 0 and 0.95 with a 400 × 250 grid.
b) Pressure air water h = 0 h = 0.95 t = 0.2ms t = 0.4ms t = 0.8ms t = 1.2ms Figure 12: (continued)
−1 0 1 −1 −0.5 0 0.5 −1 0 1 −1 −0.5 0 0.5 t = 0.2ms t = 0.4ms −1 0 1 −1 −0.5 0 0.5 −1 0 1 −1 −0.5 0 0.5 t = 0.8ms t = 1.2ms
Figure 13: Physical grid system for the generalized Lagrangian run with h = 0.95 shown in Fig. 12. For clarity, a grid coarsening factor of 5 in each x1- and x2-direction is used to make each of the graphs.
−1 0 1 0 0.5 1 −1 0 1 0 0.5 1 −1 0 1 0 0.5 1 −1 0 1 0 0.5 1 −1 0 1 0 0.1 0.2 0.3 0.4 −1 0 1 0 0.1 0.2 0.3 0.4 −1 0 1 0 0.1 0.2 0.3 0.4 −1 0 1 0 0.1 0.2 0.3 0.4 t = 0.2ms t = 0.4ms t = 0.8ms t = 1.2ms ρ (M g / m 3 ) p (G P a ) x2(m) x2(m) x2(m) x2(m)
Figure 14: Cross-sectional plots of the results for the run shown in Fig. 12 along x1= 0, where the solid
0 0.2 0.4 0.6 0.8 1 1.2 10−1 100 101 102 J ave J min J max time (×10−3ms) lo g J
Figure 15: Time history of the grid Jacobian J for the generalized Lagrangian run shown in Fig. 12, where Jave, Jmin, and Jmax denote the average, minimum, and maximum values of J, respectively.
[3] J. W. Banks, D. W. Schwendeman, A. K. Kapila, and W. D. Henshaw. A high-resolution Godunov method for compressible multi-material flow on overlapping grids. J. Comput. Phys., Article in Press, 2006.
[4] D. A. Calhoun, C. Helzel, and R. J. LeVeque. Logically rectangular grids and finite volume methods for PDEs in circular and spherical domains. preprint, 2006.
[5] J. Cheng and C.-W. Shu. A high order ENO conservative Lagrangian type scheme for the compress-ible Euler equations. J. Comput. Phys., 227:1567–1596, 2007.
[6] P. Colella. Glimm’s method for gas dynamics. SIAM J. Sci. Stat. Comput., 3(1):76–110, 1982.
[7] P. Colella and H. M. Glaz. Efficient solution algorithms for the Riemann problem for real gases. J. Comput. Phys., 59:264–289, 1985.
[8] R. Courant and K. O. Friedrichs. Supersonic Flow and Shock waves. Wiley-Interscience, New York, 1948.
[9] B. Despres and C. Mazeran. Lagrangian gas dynamics in two dimensions and Lagrangian systems. Arch. Rat. Mech. Anal., 178:327–372, 2005.
[10] G. R. Gisler, R. P. Weaver, C. L. Mader, and M. L. Gittings. Two- and three-dimensional asteroid impact simulations. IEEE Computing in Science and Engineering, May/June:46–55, 2004.
[11] E. Godlewski and P.-A. Raviart. Numerical Approximation of Hyperbolic Systems of Conservation Laws. Springer-Verlag, 1996. Applied Mathematical Science 118.
[12] J. Grove and R. Menikoff. The anomalous reflection of a shock wave at a material interface. J. Fluid Mech., 219:313–336, 1990.
[13] B. Gustafsson, H.-O. Kreiss, and J. Oliger. Time Dependent Problems and Difference Methods. John Wiley & Sons, Inc., New York, 1995.
[14] A. Harten, P. D. Lax, and B. van Leer. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Review, 25:35–61, 1983.
[15] K. A. Hoffmann and S. T. Chiang. Computational Fluid Dynamics. Engineering Education System, Wichita, Kansas, USA, 4 edition, 2000.
[16] W. H. Hui. The unified coordinate system in computational fluid dynamics. Comm. Comput. Phys., 2:577–610, 2007.
[17] W. H. Hui and S. Koudriakov. A unified coordinate system for solving the three-dimensional Euler equtions. J. Comput. Phys., 172:235–260, 2001.
[18] W. H. Hui and S. Koudriakov. Computation of the shallow water equations using the unified coordinates. SIAM J. Sci. Comput., 22(5):1615–1654, 2002.
[19] W. H. Hui, P. Y. Li, and Z. W. Li. A unified coordinate system for solving the two-dimensional Euler equtions. J. Comput. Phys., 153:596–637, 1999.
[20] P. Jia, S. Jiang, and G. P. Zhao. Two-dimensional compressible milti-material flow calculations in a unified coordinate system. Computers and Fluids, 35:168–188, 2006.
[21] C. Jin and K. Xu. A unified moving grid gas-kinetic method in Eulerian space for viscous flow computation. J. Comput. Phys., 222:155–175, 2007.
[22] D. Kincaid and W. Cheney. Numerical Analysis. Brooks/Cole, 1990.
[23] H.-O. Kreiss and J. Lorenz. Initial-Boundary Value Problems and Navier-Stokes Equations. Acad-emic Press, San Diego, 1989.
[24] A. G. Kulikovskii, N. V. Pogorelov, and A. Yu. Semenov. Mathematical aspects of Numerical Solution of Hyperbolic Systems. Chapman & Hall/CRC, Boca Raton, 2001.
[25] R. J. LeVeque. Numerical Methods for Conservation Laws. Birkh¨auser-Verlag, 2 edition, 1992.
[26] R. J. LeVeque. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, 2002.
[27] R. J. LeVeque and K.-M. Shyue. One-dimensional front tracking based on high resolution wave propagation methods. SIAM J. Sci. Comput., 16:348–377, 1995.
[28] M.-S. Liou. An extended Lagrangian method. J. Comput. Phys., 118:294–309, 1995.
[29] H. Luo, J. D. Baum, and R. L¨ohner. On the computation of multi-material flows using ALE formu-lation. J. Comput. Phys., 194:304–328, 2004.
[30] P.-H. Maire, R. Abgrall, J. Breil, and J. Ovadia. A cell-centered Lagrangian scheme for two-dimensional compressible flow problems. SIAM J. Sci. Comput., 29:1781–1824, 2007.
[31] S. P. Marsh. LASL Shock Hugoniot Data. University of California Press, Berkeley, 1980.
[32] D. J. Mavriplis and Z. Yang. Construction of the discrete geometric conservation law for high-order time-accurate simulations on dynamic meshes. J. Comput. Phys., 213:557–573, 2006.
[33] H. J. Melosh. Impact Cratering: A Geologic Process. Oxford University Press, 1989.
[34] R. Menikoff and B. Plohr. The Riemann problem for fluid flow of real materials. Rev. Mod. Phys., 61:75–130, 1989.
