A wave-propagation based volume tracking method
for compressible multicomponent flow in two space dimensions
Keh-Ming Shyue
*Department of Mathematics, National Taiwan University, Taipei 106, Taiwan Received 18 June 2005; received in revised form 25 October 2005; accepted 27 October 2005
Abstract
We present a simple volume-of-fluid approach to interface tracking for inviscid compressible multicomponent flow problems in two space dimensions. The algorithm uses a uniform Cartesian grid with some grid cells subdivided by tracked interfaces, approximately aligned with the material interfaces in the flow field. A standard volume-moving procedure that consists of two basic steps: (1) the update of a discrete set of volume fractions from the current time to the next and (2) the reconstruction of interfaces from the resulting volume fractions, is employed to find the new location of the tracked interfaces in piecewise linear form at the end of a time step. As in the previous work by LeVeque and the author to front tracking based on a surface-moving procedure (R.J. LeVeque, K.-M. Shyue, Two-dimensional front tracking based on high-resolution wave propagation methods, J. Comput. Phys. 123 (1996) 354–368), a conservative high-resolution wave propagation method is applied on the resulting slightly non-uniform grid to update all the cell values, while the stability of the method is maintained by using a large time step idea even in the presence of small cells and the use of a time step with respect to the uniform grid cells. We validate our algorithm by performing the simulation of a Mach 1.22 shock wave in air over a circular R22 gas bubble, where sensible agreement of some key flow features of the computed solutions are observed when direct comparison of our results are made with the existing experimental and numerical ones that appear in the lit- erature. Other computations are also presented that show the feasibility of the algorithm together with a mixture type of the model equations developed by the author (K.-M. Shyue, A fluid-mixture type algorithm for compressible multicom- ponent flow with Mie–Gru¨neisen equation of state, J. Comput. Phys. 171 (2001) 678–707) for practical multicomponent problems with general compressible materials characterized by a Mie–Gru¨neisen form of the equation of state.
2005 Elsevier Inc. All rights reserved.
AMS: 65M06; 76L05; 76M20; 76T05
Keywords: Volume tracking; Wave propagation method; Multicomponent flows; Mie–Gru¨neisen equation of state; Impact problem
0021-9991/$ - see front matter 2005 Elsevier Inc. All rights reserved.
doi:10.1016/j.jcp.2005.10.030
* Tel.: +8862 3366 2866; fax: +8862 2391 4439.
E-mail address:shyue@math.ntu.edu.tw.
Journal of Computational Physics xxx (2005) xxx–xxx
www.elsevier.com/locate/jcp
1. Introduction
Our goal in this paper is to describe a simple volume tracking approach for the efficient numerical resolu- tion of material interfaces arising from inviscid compressible multicomponent problems in two space dimen- sions. We use a Eulerian formulation of the equations in which the principal motion of the single-component flow is governed by the basic conservation laws for the mass, momenta in the x- and y-direction, and energy as
o ot
q qu qv qE 0 BB B@
1 CC CAþ o
ox
qu qu2þ p
quv qEuþ pu 0
BB B@
1 CC CAþ o
oy
qv quv qv2þ p qEvþ pv 0
BB B@
1 CC
CA¼ 0; ð1Þ
respectively. Here q denotes the density, u is the particle velocity in the x-direction, v is the particle velocity in the y-direction, p is the pressure, and E is the specific total energy. To complete the model, the constitutive law for the thermodynamic behavior of the fluid component of interests is assumed to satisfy a Mie–Gru¨neisen equation of state of the form
pðq; eÞ ¼ prefðqÞ þ CðqÞ qe qe½ refðqÞ; ð2Þ
where e represents the specific internal energy, C = (1/q)(op/oe)jqis the Gru¨neisen coefficient, and pref, erefare the properly chosen states of the pressure and internal energy along some reference curve (e.g., along an isentrope, a single shock Hugoniot, or other empirically fit curves) in order to match the experimental data of the material being examined, see [54] and references cited therein for more details. As usual, we have E = e + (u2+ v2)/2.
As in the previous work by LeVeque and the author[31]to general front tracking for nonlinear hyperbolic systems of conservation laws in two dimensions, our approach to volume tracking is based on a fully discret- ized version of the high-resolution finite volume method in wave-propagation form. The algorithm uses an underlying uniform Cartesian grid with some rectangles subdivided by the tracked interfaces into two parts, where discontinuities in the material interfaces (or slip lines for that matter) are expected. Rather than employ- ing a marker-and-cell representation of the interface and a surface-moving procedure to advance the interface location as before, which can be done for some simple cases (cf.[51,52]) but is very troublesome (at least in view of the code development aspect) in dealing with situations such as the splitting or merging of interfaces and the changes of complex multidimensional interfacial topology (cf.[9]), we adopt a volume-fraction based interface-moving procedure for an easy treatment of the aforementioned issues instead, see [19,38] for the other possible alternates. That is to say, we introduce a scalar field Y in the algorithm that is defined to be the fraction of each cell that is on one particular side of the interface. In the current interest of two-dimen- sional compressible flow problems, this volume-fraction function Y would be updated in each time step by solving an evolution equation of the form
oY ot þ uoY
oxþ voY
oy ¼ 0 ð3Þ
with a standard wave-propagation method. Here u and v are the particle velocities governed by the full Euler equations(1) which is unlike the case for many incompressible flow problems where a divergence-free con- straint should be fulfilled on them for the fundamental agreement of the mathematical model being used.
When this is done, the new interface can then be reconstructed in a cell-by-cell fashion from the information on Y, yielding a set of disconnected line segments (this is as opposed to a connected line segment case obtained by using a typical surface-moving procedure[6,31]) in the cells containing parts of the interface, see Section3 for more exposition.
An important property of the wave-propagation form of the method is that reasonable time steps can be taken even if some of the subcells created by the tracked interfaces are orders of magnitude smaller than the underlying Cartesian cells. Uniform time steps are used throughout the computation, with the time step cho- sen so that the usual CFL (Courant–Friedrichs–Levy) condition is satisfied relative to the size of the regular grid cells. The method would typically remain stable in this case because near the small subcells the numerical
domain of dependence is extended automatically based on a large time step technique developed by LeVeque for general hyperbolic systems of partial differential equations (cf. [24–27]). High-resolution results can be achieved in a straightforward manner in most part of the domain, except possibly for cells near the tracked interfaces. This is discussed further in Section2.
It should be mentioned that as in our surface-based front tracking algorithm[31], the basic philosophy of the proposed volume tracking method is that the subdivision of cells is not assumed to give the definitive loca- tion of the true interface, but is viewed rather as an approximate location yielding a refined grid that is better able to represent the discontinuous solution than the Cartesian grid alone. As a result, there is minimal smear- ing of the interface in the solution. To prevent the occurrence of spurious oscillation in the pressure when such a scenario occurs in the multicomponent context, the algorithm uses a mixture type of the model system described in[54]as the basis, and a standard finite volume method based on the wave-propagation formula- tion as for the consistent and accurate approximation of the system. Numerical results presented in Section5 show that this gives an efficient way for the improvement of the resolution of the interfaces, and the correct convergence behavior of the weak solution for shock waves.
It is without doubt that there are many instances to have the desire of keeping the interface sharp, and con- sequently the advantage of a tracking method over a fixed grid capturing method. One good example among them is in the simulation of compressible two-phase flows with very different fluid component separated by interfaces, and under extreme flow condition, see[18,32,37,65,66,69]for an example. In this case, despite some simple situations (cf.[54]), numerical experiences indicate that it would be very difficult in general to devise a pressure-oscillation-free interface-capturing method to deal with the occurrence of more than one fluid com- ponent within a grid cell. In fact, it is one of the main motivations of the present attempt to develop a robust volume tracking method that is viable for a wide variety of two-phase flow problems with real materials.
There have been quite a few efforts over the years to develop efficient volume-of-fluid interface tracking methods for various practical applications. In regards to incompressible flow problems governed by the Navier–Stokes equations, for example, representative works in this direction are such as Hirt and Nichols for free boundary problems [16], Welch and Wilson for phase transition problems[67], and Gueyffier et al.
[13], Lafaurie et al. [22], Rudman [47,48], and Sussman and Puckett [59] for two-phase flow problems. In regards to the current compressible flow problems governed by the inviscid Euler equations, there are works done by Chern and Colella[5]for nonlinear hyperbolic systems of conservation laws (see[1]also for a version of the method with the local adaptive mesh refinement included), and Miller and Puckett[37]for problems in multiple condensed matters.
The approach we take here is in essence similar to that of Chern and Colella [5]in which they also use a volume-of-fluid approach to keep track of the interface location, and then apply high-resolution conser- vative finite volume methods on the resulting grid. The key difference is that they use a more traditional flux difference form, together with a flux redistribution algorithm to maintain stability in the presence of small cells. It has been demonstrated before (cf. [31,52]) that our wave-propagation form can give good results when it is combined with a surface-based front tracking algorithm in one and two dimensions.
One major objective of this paper is to show that this method can still be very effective when it is used in conjunction with a volume-based interface tracking method for two-dimensional compressible hydrody- namics problems.
The format of this paper is outlined as follows. In Section2, we review the basic idea of a finite volume method in wave-propagation form on a slightly non-uniform Cartesian grid with a moving interface. In Sec- tion3, we describe the standard piecewise linear interface construction method based on a given set of discrete volume fractions. The volume tracking algorithm for inviscid compressible 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 Section5.
2. Wave-propagation methods
To set the groundwork for the later development of a volume-based interface tracking algorithm, we begin our discussion by reviewing a finite volume method in wave-propagation form for the numerical approxima- tion of general hyperbolic systems of conservation laws in two space dimensions,
oq
otþofðqÞ
ox þogðqÞ
oy ¼ 0. ð4Þ
Here q2 Rmdenotes the vector of m conservative quantities, and f(q), g(q) denote the flux vectors, see(1)for an example. For all practical purposes as concerned below, it is essential to describe the method on a slightly non-uniform grid that is composed of a set of regular cells with a fixed mesh spacing Dx and Dy in the x- and y-direction, respectively, and some irregular cells which are formed by inserting a piecewise linear representa- tion of the tracked interfaces into the grid, yielding the subdivision of some regular cells into two subcells, see Fig. 1a for an illustration.
In a finite volume formulation of the method, the cell average of each grid cell is defined by integrating the solution over the cell and dividing by the area of the cell. Hence the approximate value of the cell average of the solution over any given cell S at time tncan be written as
QnS 1 MðSÞ
Z
S
qðx; y; tnÞ dx dy; ð5Þ
where MðSÞ is the measure (area) of this cell. The methods we use are based on solving one-dimensional Rie- mann problems at each cell edge, and the waves (i.e., discontinuities moving at constant speeds) arising from the Riemann problem are employed to update the cell averages in the cells neighboring each edge (cf.
[26,27,31]).
2.1. First-order method
As a first example to explain the basic idea of our first-order method, we consider the cell edge between irregular cells C and D as illustrated inFig. 1a. We solve the one-dimensional Riemann problem normal to this face, which in this case will be
Fig. 1. First-order wave-propagation method. (a) The typical grid created by our volume tracking algorithm. (b) A wave propagating normally from the cell edge between cells C and D updates the solution in both cells D and F. (c) A wave propagating normally from the cell edge between cells C and E updates four neighboring cell averages. (d) A wave propagating both normally and tangentially from the cell edge between cells C and D updates the solution in cells D, F, and H. (e) A wave propagating both normally and tangentially from the cell edge between cells C and E updates five neighboring cell averages. Here Wprepresents the region swept out by the p-wave, and Wpq
represents the region swept out by the p-wave and q-wave in the normal and tangential to the cell edge, respectively, over a time step Dt for p,q = 1, 2, . . . , m.
oq
otþofðqÞ ox ¼ 0;
with the initial data given by QnC and QnD. If RoeÕs approximate solver is chosen to solve the Riemann problem, the above nonlinear equation would be replaced by a linear system of the form
oq
otþ A Q nC; QnD oq ox¼ 0;
where AðQnC; QnDÞ is a constant matrix that depends on the initial data and is a local linearization of the Jaco- bian matrix of f(q) about an average state (cf.[29,46,62]). The result is that the jump QnD QnC is decomposed into eigenvectors of the matrix A,
QnD QnC¼Xm
p¼1
aprp;
where rpis the pth eigenvector of A and Arp= kprpwith kpthe corresponding eigenvalue. The scalar apgives the wave strength.
Clearly, if the wave speed kp> 0 the pth wave propagates towards the right, while if kp< 0 it propagates towards the left.Fig. 1b shows a typical p-wave with kp> 0, where part of cell D and part of cell F are affected by this wave over a time step Dt from the current time tnto the next tn+1. In the simplest case of the wave- propagation method, i.e., a two-dimensional extension of the first-order GodunovÕs method, these cell aver- ages are updated by
Qnþ1S :¼ Qnþ1S MðWp\ SÞ
MðSÞ aprp ð6Þ
for S = D, F. Here Wprepresents the region swept out by the wave and aprpis the jump across the wave. Note that we have employed a standard initialization procedure that assigns each cell average to its values at the previous time step, Qnþ1S :¼ QnS for all S.
Now, if we are concerned with the cell edge between cells C and E (this is one segment of the tracked inter- face), a Riemann problem would be solved in the direction normal to the interface with data QnC and QnE. Note that the Euler equations (1)are isotropic and easily rotated to this frame (cf. [26]). The p-wave indicated in Fig. 1c overlaps four cells and using (6) in a similar manner gives modification of the cell average in each of these cells by the jump across this wave, weighted by the fraction of the cell area covered by the wave.
To improve the stability of this first-order method, it is a customary approach by introducing transverse propagation of waves based on the solutions of Riemann problems in the orthogonal direction (cf.
[26,27,31]). Again we use the cell edge between cells C and D inFig. 1a, as an example. Then for the equation we take the tangential portion of(4), i.e., the portion of equation in the y-direction,
oq
otþogðqÞ oy ¼ 0;
and for the initial data if we are working with the p-wave we take QL¼ QnCþXp1
q¼1
aqrq; QR¼ QLþ aprp.
With the Roe solver, the solution to this Riemann problem yields a splitting of QR QL= aprpinto eigenvec- tors of the Roe matrix B(QL, QR) of the form
aprp¼Xm
q¼1
bqwq
where B is a suitable approximation of the Jacobian matrix of g(q), and B wq= lqwq. The wave shown in Fig. 1b is thus split into m waves propagating in the y-direction with speed given by the corresponding
eigenvalue lq, and a typical one is shown in Fig. 1d for the case lq> 0. Here to update the cell averages af- fected by this wave, we employ an analogous formula as in(6),
Qnþ1S :¼ Qnþ1S M W pq\ S
MðSÞ bqwq ð7Þ
for S = D, F and H, where Wpqrepresents the region swept out by the p-wave and q-wave over a time step Dt, and bqwqis the jump across the wave as usual.
We note that, using the same idea as described above, there is no difficulty to split any p-wave into m sub- waves Wpqfor q = 1, 2, . . . , m, moving oblique to the cell edge rather than normal to it. Each of these waves can then be propagated in a similar manner, computing its intersection with the neighboring cells and updat- ing the cell averages, seeFig. 1e for another illustration of a subwave, lq> 0, from the cell edge between cells C and E.
It is important to mention that the method remains fully conservative when the transverse propagation of waves is included in the algorithm. Furthermore, because waves are allowed to propagate through more than one grid cell and are not confined to stay within neighboring cells, this method is typically stable as long as the time step Dt satisfies the usual CFL (Courant–Friedrichs–Lewy) condition relative to the regular cell sizes: Dx and Dy,
m¼ Dtmax
p;q ðkp;lqÞ
minðDx; DyÞ 61; ð8Þ
even when there are very small irregular cells near the tracked interfaces. Nevertheless, the method is still only first-order accurate, see[26,28]for more details.
2.2. High-resolution method
To achieve high resolution of the method (i.e., second-order accurate on smooth solutions, and sharp and monotone profiles on discontinuous solutions), as usual, we use the solution of the Riemann problem in the direction normal to each cell edge as a basis for the introduction of correction waves in a piecewise-linear form with zero mean value. Having that, we then propagate each of these waves by the corresponding speed over a time step Dt to update the cell averages it overlaps. To demonstrate this,Fig. 2a shows a typical structure of a correction wave at time tnfor a case with kp> 0 from the cell edge between regular cells G and H, i.e., the graph of the three-dimensional piecewise-linear function
Zpðx; yÞ ¼ rpðx xGÞ if ðx; yÞ 2 ½xG Dx=2; xGþ Dx=2 ½yaGH; ybGH;
0 if otherwise,
ð9Þ where rp= aprp/Dx is the slope of Zp, xGis the x-coordinate of the geometric center of cell G, and yaGH; ybGHare the lower and upper ends of the y-coordinate of this cell edge, respectively. Note that, when kp< 0, Zpcan be made quite easily also by simply changing xGin(9)to xH. The new location of Zpmoving in the x-direction over a time step Dt is Zpðx kpDt; yÞ, seeFig. 2b. The cell averages G and H are then modified by
Qnþ1G :¼ Qnþ1G 1 2
jkpjDt
Dx Dx jkpjDt rp; Qnþ1H :¼ Qnþ1H þ1
2 jkpjDt
Dx Dx jkpjDt rp;
ð10Þ
where the correction terms are the volume of the portion of Zpðx kpDt; yÞ (weighted by the cell size) that overlaps the grid cell. In practice, the strength of each wave should be limited by using a slope-limiter (cf.
[10,29]), and so each ap in the above rp is replaced by a limited value eap obtained by comparing ap with the corresponding ap from the neighboring Riemann problem to the left (if kp> 0) or to the right (if kp< 0). Conservation is clearly maintained in this step with any choice of slopes, because the above two cor- rection terms sum to zero.
To extend the above approach to the more delicate case where the cell edge is adjacent to one or two irreg- ular cells, we may first define a general form of Zpin which the normal direction to the cell edge is at an arbi- trary angle h to the x-axis. Without loss of generality, we consider the one between cells C and E shown in Fig. 1a, as an example. Then, for the p-wave, one simple choice of the slope is rp= aprp/Dd, where Dd is some measure of the normal distance between cells C and E, such as Dd ¼ nE nC with nS ¼ xScos hþ ySsin h rep- resenting the n-coordinate of the geometric center of cell S in a local (n, g)-coordinate system (i.e., n¼ x cos h þ y sin h, g ¼ x sin h þ y cos h), for S = C, E (cf. [26]). With these definitions, for a case with kp> 0, we form the correction wave as
Zpðn; gÞ ¼ rpn n0
if ðn; gÞ 2 ½n0 Dd=2; n0þ Dd=2 ½gaCE;gbCE;
0 if otherwise,
(
ð11Þ where n0is the n-coordinate of the center of the rectangular region:½nCE;nCE Dd ½gaCE;gbCE, seeFig. 2c for an illustration. Having that, the new location of Zppropagating in the normal direction to the cell edge over a time step Dt is Zpðn kpDt;gÞ, seeFig. 2d. Now, the cell averages affected by this correction wave can be up- dated by
Qnþ1S :¼ Qnþ1S þM Z p\ S
MðSÞ ; ð12Þ
where MðZp\ SÞ is the measure (volume) of the portion of Zpðn kpDt;gÞ that partly coincides the cell S. It is true that the unlimited slope employed in(11)would be a sensible one to use only if there is no large jump in the solution across the wave. As a general practice, it should be good also to limit the strength of each wave ap
using a slope-limiter in an upwind manner based on the sign of kp. But if we now look at the case shown in Fig. 2c, it is not obvious right away how to obtain a limited value ~ap, because we do not know exactly what is the most appropriate neighboring apto be used in the slope limiter. There are some discussions in[26]on how such corrections might be applied on irregular cells, and the algorithm would perhaps be improved by imple- menting this. However, since there is expected to be a large jump at the tracked interface, this is precisely
Fig. 2. High-resolution wave propagation method. (a) A correction wave at tnfor a case with kp> 0 from the cell edge between cells G and H. (b) Propagation of the correction wave in (a) with the speed kp> 0 over a time step Dt updates the solution in cells G and H. (c) A correction wave at tnfor a case with kp> 0 from the tracked interface between cells C and E. (d) Propagation of the correction wave in (c) with the speed kp> 0 over a time step Dt updates the solution in neighboring cells. SeeFig. 1a for the basic labeling of the grid cell.
where the limiters are expected to minimized the effect of these second-order corrections and to reduce the method to Godunov. It seems reasonable to drop these corrections altogether with considerable simplification of the algorithm.
To end this section, it should be mentioned that the wave-propagation method described here is not limited to solve hyperbolic systems in conservation form, but works equally well to solve hyperbolic systems which are not in conservation form, see[29]for the details.
3. Interface reconstruction scheme
Denote Yijas being the approximate value of the cell average of the solution(3)over the rectangular region Cij= {(x, y) | xi6x 6 xi+ Dx, yj6y 6 yj+ Dy} of grid cell (i, j). Then, when Yij2 (0, 1), we expect there is an interface (represents as a linear segment) lying somewhere within Cij. Here we are concerned with a popular piecewise linear interface construction (PLIC) method (cf.[42,45]) for the reconstruction of interfaces from the discrete set of volume fractions Y on uniform Cartesian cells which will be one of the fundamental elements in our volume tracking algorithm, see Section4. Note that in a standard PLIC method each of the linear seg- ments would be made by first computing an approximate normal vector to the interface within the grid cell, and then determining its location in such a way that the cell volume fraction is preserved in this step. In the following, we review the aforementioned PLIC steps in more details.
3.1. Compute interface normal
As a preliminary, we consider a gradient method of Parker and Youngs[40,70]in which, when Yij2 (0, 1), the gradient of the volume fraction at the grid cell (i, j), denoted by $Yij, is computed explicitly by a central- difference like formula of the form
rYij¼ oY ox;oY
oy
ij
YR YL
2Dx ;YT YB
2Dy
ð13Þ with the quantities YR, YL, YT, and YB defined as
YR ¼14ðYiþ1;j1þ 2Yiþ1;jþ Yiþ1;jþ1Þ; YL¼14ðYi1;j1þ 2Yi1;jþ Yi1;jþ1Þ;
YT ¼14ðYi1;jþ1þ 2Yi;jþ1þ Yiþ1;jþ1Þ; YB¼14ðYi1;j1þ 2Yi;j1þ Yiþ1;j1Þ.
Having that, the unit normal of the interface at cell (i, j) may be set to ~Nij¼ rYij=jrYijj, and so in this cell
~Nij points in the direction of the most rapidly decrease of Y. While the use of(13)is a very simple one for the determination of the interface normal, it is reported in[42], for instance, that this method exhibits only overall first-order accuracy on the location of the interface. To achieve higher-order accuracy, such as the second order, one popular approach is to use the least-squares method of Puckett[42,43]which is described next.
Let Eij be the sum squared error of volume fractions for the cell (i, j) of the form
Eij¼ Xiþ1
m¼i1
Xjþ1
n¼j1
Ymn ~Ymn
2
; ð14Þ
where Ymnis the original volume fraction for cell (m, n), and ~Ymnis the approximate volume fraction for cell (m, n) based on an extended reconstructed interface at cell (i, j). The aim of the PuckettÕs method is to obtain a local minimum of Eij by adjusting $Yij(and so the location of the interface) iteratively in the center of a 3· 3 mesh block until convergence. Note that, during each iteration, i.e., for any given $Yij, we have to perform an interface-finding procedure, see Section3.2, so that the condition ~Yij¼ Yij is fulfilled in the center cell. Here the numerical method used to solve the minimization problem is a combination of the golden section search and the successive parabolic interpolation (cf.[3]and the NETLIBroutineFMIN). To ensure convergence of the correct result by the method, it is important to choose a suitable interval to work with in which Eij attains a minimum. In practice, this is typically done by a trial-and-error procedure (cf. [15,42,50]) where an initial guess of $Yijis picked up first, say based on(13)or more favorably the YoungsÕ second method (see below),
and then it is expanded slowly about this initial guess to an interval until Eij at the endpoints of the interval has a bigger value than the one given by the initial guess.
Note that in PuckettÕs method $Yijis determined implicitly and also rather costly as a result of the least- squares problem(14). It is interesting to mention that a much simpler method of this type is YoungsÕ second method[71]where $Yijis obtained by solving a weighted least-squares problem of the form
Exij ¼ Xiþ1
m¼i1
Xjþ1
n¼j1
xmnðYmn ~YmnÞ2; ð15Þ
where in this work xmnis taken as the inverse distance from the (m, n)th cell center to the (i, j)th cell center, and Y~mnis defined to be a Taylor series expansion of the volume fraction from the (m, n)th cell center to the (i, j)th cell center, i.e.,
Y~mn¼ Yijþ rYij ðxmþ1=2 xiþ1=2; ynþ1=2 yjþ1=2Þ.
With that, the solution of the minimization problem to(15), denoted by z = ($Yij)T, can be shown to satisfy the linear system of the form
ATAz¼ ATb; ð16Þ
where
A¼
Dx=Ds Dy=Ds
1 0
Dx=Ds Dy=Ds
0 1
0 1
Dx=Ds Dy=Ds
1 0
Dx=Ds Dy=Ds 0
BB BB BB BB BB BB BB
@
1 CC CC CC CC CC CC CC A
; b¼
ðYi1;j1 YijÞ=Ds ðYi1;j YijÞ=Dx ðYi1;jþ1 YijÞ=Ds
ðYi;j1 YijÞ=Dy ðYi;jþ1 YijÞ=Dy ðYiþ1;j1 YijÞ=Ds
ðYiþ1;j YijÞ=Dx ðYiþ1;jþ1 YijÞ=Ds 0
BB BB BB BB BB BB BB
@
1 CC CC CC CC CC CC CC A
;
and Ds¼
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðDxÞ2þ ðDyÞ2 q
. In general, one may want to use either a normal equation-based or a QR-based algorithm for the solution of this least-squares problem(15), see[63]. But after a simple algebraic calculation, we find that the matrix ATA is simply a diagonal matrix of the form
ATA¼ 2 diag 2 Dx Ds
2
þ 1; 2 Dy Ds
2
þ 1
!
;
and so we may find the solution of(16)very easily, and so set ~Nij¼ rYij=jrYijj as the normal of the inter- face for cell (i, j).
3.2. Determine interface location
Having gotten ~Nij, we now describe a numerical procedure (cf.[45]) for an efficient reconstruction of the location of the interface over the rectangular region Cij. For convenience of the later description, let us assume thatðx; yÞ 2 Cijis a reference point known a priori on the interface and ~Nij¼ ðaij;bijÞ with a2ijþ b2ij¼ 1. Then, the linear-segment representation of the desired location of the interface in cell (i, j) in parametric form may be expressed by LijðsÞ : ðxðsÞ; yðsÞÞ ¼ ðx bijs; yþ aijsÞ, where s is a non-negative parameter value. Clearly, with the proper choice of the interval on s, LijðsÞ would form a line, subdividing Cijinto two parts. Denote Cij and Cþij to be the region occupied by subcells in the back and in the front of Lij, respectively. Our aim here is to determine ðx; yÞ and so Lij in such a way that MðCijÞ=MðCijÞ ¼ Yij. Surely, the latter condition ensures the basic conservation of the volume fraction if we set Y = 1 for Cij and Y = 0 for Cþij. Note that this setting is consistent with our definition of the direction of ~Nij, i.e., a unit vector pointing in the steepest descent direction of Y. Recall that MðSÞ is the measure (area) of the region S.
To accomplish this, we define eðkÞij ¼ MðC;kij Þ=MðCijÞ Yijas being the error of the cell volume fraction (i, j) at the iteration step k, where MðC;kij Þ is the measure of C;kij using a reference point at (x(k), y(k)). With that, by choosing a suitable initial guess: (x(0), y(0)) and (x(1), y(1)) such that eð0Þij eð1Þij <0, it is not difficult to apply the standard secant method for finding the zero of eðkÞij , yielding the convergence of the solution ðxðkÞ; yðkÞÞ ! ðx; yÞ in a superlinear rate. Note that, given Lij, the basic geometric information of the splitted grid cell, such as the area and center of the cell, can be found quite easily based on algorithms in computa- tional geometry[39].
Finally,Fig. 3shows results for a typical situation in a volume-of-fluid algorithm for interface reconstruc- tion from volume fractions, where we have used both the methods of Parker and Youngs, and Puckett on a very coarse 5· 5 grid for illustrations, see [42] for an extensive study of the accuracy of these schemes and many others. Note that, in the graphs, we have set Y = 1 for grid cells in the shaded region, and Y = 0 for grid cells in the non-shadded region. As mentioned earlier, the tracked interface is approximated by a piece- wise linear segment with knots at the points where this segment intersects with the grid lines. We use the same data structure as described in[31]for the realization of the tracked interfaces in our computer program.
4. Volume tracking algorithm
In each time step, our volume-based interface tracking method for compressible fluid dynamics problems in two dimensions consists of the following steps:
(1) Take a time step on the current non-uniform grid using a first-order version of the wave-propagation method to update the cell averages of volume fractions governed by(3)at the next time step.
(2) Find the new interface location based on the volume fractions obtained in step 1 using an interface reconstruction scheme. Some cells will be subdivided, and the values in each subcell must be initialized.
(3) Take the same time interval as in step 1, but use the high-resolution wave-propagation method to update the cell averages of our model system such as(1) and (2)for single component flow or a fluid–mixture model proposed in[54]for multicomponent flow on the new non-uniform grid created in step 2.
Note that steps 1 and 2 constitute the basic volume-moving procedure for the modeling of the propagation of interfaces, and step 3 is simply the construction of approximate solutions to the physical problems of inter- ests on a carefully chosen non-uniform grid. Since many of the key ingredients of the algorithm have been discussed before, see Sections2 and 3, our goal here is to provide some more details with application to com- pressible flow problems.
4.1. Volume fraction update
Without loss of generality, suppose thatFig. 3c is the grid system we are working on at time tnwith the approximate state for the volume fraction of (3) defined as YnS¼ 0 and YnS¼ 1 when each of the grid cells
Fig. 3. An example of the reconstruction of the interface location from volume fractions. (a) The original cell-averaged volume fraction functions for an elliptical-shape interface. (b) Result obtained by using the method of Parker and Youngs. (c) Result obtained by using the method of Puckett. Here in the graphs we set Y = 1 and Y = 0 for grid cells in the shaded and non-shadded region, respectively.
S is in the non-shaded and shaded regions of the graph, respectively. In the meantime, the state vector for our (single or multiple component) physical model, denoted by QnS, is set by(5)to represent the approximate value of the cell average of the solution over any given cell S. The algorithm uses a less stringent time step which is determined by the CFL condition(8)based on the uniform grid sizes and the appropriate estimate of the wave speeds in QnS.
Notice that in this setup the value of Ynis only having jump with magnitude 1 across some polygonal edges of irregular cells, and is uniform elsewhere. Thus, in step 1, it should be reasonable to simply employ the wave- propagation method in first-order form (preferably the one with the transverse propagation of waves, see Section2.1), to find the new volume fraction Yn+1on this non-uniform grid at the next time step. Here the basic solution of the one-dimensional Riemann problem in the direction normal or tangential to the cell edge can be computed quite easily by solving the variable-coefficient linear transport equation(3)with the velocity field came from Qn for the coefficient, and Yn for the initial data (cf. [29]), yielding analogously the wave structure illustrated in Fig. 1for m = 1 over a given time step Dt.
With the above comments in mind,Fig. 4a gives numerical result of a one-step calculation for the volume fraction on the uniform rectangular grid, where the computation is performed by assuming a constant velocity field (u, v) = (1, 1) throughout the domain, and a time step Dt = 0.06. We note that in case the grid cell is an irregular one, the result present in the figure is computed by an area-weighted average of the original data on the subcells; this is necessary in the algorithm because the reconstruction of the new interface location to be done next depends on that. It is not difficult to show that our result is in good agreement as in comparison with the exact solution of this simple model problem.
4.2. New interface reconstruction
When the new volume fraction Yn+1becomes available on the uniform grid, the approximate location of the tracked interface at the next time step can be obtained in principle by applying any interface reconstruction scheme described in Section3to the data set. If we now employ the PLIC method of Puckett to the case shown inFig. 4a, for instance, we would have the disconnected set of linear segments as displayed inFig. 4b for the representation of the new location of the old tracked interface over the time step Dt = 0.06. Surely, when com- paring to the exact location of the interface, one may argue that our volume-based result does not seems to be superior to the result gotten by using such as a more typical surface-based approach for the same problem (not shown here, see [58] for more details). But as far as demonstrating the essential feature of the method is concerned, this simple example serves well for our purpose here (some sample results that demonstrate the feasibility of our method to practical problems will be present in Section5). Nevertheless, we want to empha- size again that in our algorithm we do not assume the tracked interface to be the definitive location of the true
Fig. 4. An example of the volume-moving procedure in our interface tracking algorithm with a uniform velocity field (u,v) = (1,1) throughout the domain. (a) The updated volume fraction on the uniform rectangular grid based on the solution of(3)with the initial condition shown inFig. 3c over a time step Dt = 0.06. Here j1= 5.7· 103and j2= 1.3· 103. (b) New interface location obtained by using the PLIC method of Puckett from the volume fractions appeared in (a), where the dotted-line represents the initial location of the tracked interface. Note that the volume fraction at the next time step is set to be unity and zero for grid cells in the shaded and non- shadded region, respectively.
interface, but is only considered as an approximate location yielding a refined grid that is better able to rep- resent the solution for the interface than the uniform grid alone.
It is clear that when a new interface is created, some cells will be subdivided, and the values in each subcell must be assigned. Here the method we use is implemented by taking a data structure that keeps the old and new grid systems separate, and as discussed in[31] doing so would be more efficient than the use of a data structure for an intermediate grid that contains both the old and new subdivisions. In this manner, when the subcell is in front and back of the (directed) new tracked interface within the grid cell, we set the volume fraction at the next time step to be zero and unity, respectively. This yields the new volume fraction that will be used in the next run of our algorithm.
In addition to that, in this work we also need to initialize the state variables for the compressible model of interests at the current time so that they can be updated in the subsequent step. For that, we should distinguish the situation on whether the new tracked interface is inserted into a regular cell or not. Firstly, suppose the new interface splits a regular grid cell A into two subcells B and C. In this instance, one simple way is to assign the cell averages QnB and QnC based on the value QnA as
QnB:¼ QnA; QnC :¼ QnA; ð17Þ
see[5,30]for some discussions on the other possibilities. Secondly, we have an irregular grid cell that is sub- divided by the old interface into subcells A and B, and also by the new interface into the subcells C and D.
Now the cell averages QnC and QnD may be initialized by the weighted averages of the values QnA and QnB as QnS :¼ ½M A \ Sð ÞQnAþ M B \ Sð ÞQnB=MðSÞ for S ¼ C; D. ð18Þ
4.3. Physical solution update
Our final step of the algorithm is to find the cell averages of the solution of our (single or multicomponent) compressible model system on the new non-uniform grid that contains some of the cell edges as the approx- imate location of the new interface at the next time step. As before (cf.[31]), this is accomplished in the same manner as described in Section2 for the solution update of hyperbolic conservation laws on a slightly non- uniform grid. That is to say, with a properly chosen Riemann solver (cf.[52–54]), we first solve the associated one-dimensional Riemann problems in the directions normal and tangential to each of the cell edges which belongs to the original grid system when entered into the algorithm for computation. Having done so, due to the fact that we have chosen the data structure and the initializations(17) and (18)so that there are no waves generated across any polygonal edges of the irregular cell to the new grid, we may then propagate each resulting wave independently of all others to update the cell averages of the aforementioned new grid in cells neighboring to the cell edge over a given time step.
It should be remarked that as in the work of Chern and Colella[5]for interface tracking, or Berger et al.[2]
and Colella et al.[7] for handling irregular geometries in a Cartesian grid method, it is quite common to employ a post-processing step of the solution that reassigns the cell average of each irregular cells to a new one using an interpolation scheme based on the cell averages over an extended region neighboring to the grid cell under concerned. Here, in our algorithm, we take the same one-sided interpolation scheme as discussed in [5]for that matter, observing good behavior of the numerical solutions with that for many practical problems, see Section5for an example.
To end this section, we want to mention that the basic idea of the two-dimensional volume tracking algo- rithm described here can be extended in a straightforward manner to three space dimensions. That is, in steps 1 and 3 of the algorithm, we may now apply a variant of the three-dimensional wave-propagation method proposed by Langseth and LeVeque[23]on a slightly non-uniform Cartesian grid to update the volume frac- tions and the physical solutions, respectively. Clearly, one major endeavor to realize this is to devise a robust numerical procedure so that the measure (volume) of the cubical region swept out by the wave obtained from solving the one-dimensional Riemann problem over Dt and the neighboring hexahedral cells can be computed accurately and efficiently. It is fortunate however that in the field of computational geometry (cf. [8,11,39]) there have been quite a few tools developed for that purpose already, and so we may simply take and imple- ment one of the popular approaches from there in our computer program. As to the interface-reconstruction
step of the algorithm, following the work appeared in[20,21], for instance, it can be done in a similar fashion without too much difficulty. It is important to note that to overcome the outstanding small cell problem when there is a tracked interface cutting through the underlying uniform grid, we are interested in the large time step method as usual. A detailed description of this three-dimensional algorithm (cf. [57]) will be reported else- where in the future.
5. Numerical results
We now present some sample numerical results obtained using our algorithm described in Section4for track- ing material interfaces arising from compressible multicomponent problems in two dimensions. Without stated otherwise, we have carried out all the tests in this section using the Courant number m = 0.9 defined by(8), the YoungsÕ second method in the interface reconstruction scheme (for the sake of its simplicity and reasonably good prediction of the interface normal), and theMINMODlimiter in the high-resolution version of the finite- volume wave-propagation method. For comparison purposes, in many cases considered below, we have also included non-tracking results when the interface capturing version of the method is employed to the problems.
Example 5.1. To begin with, we consider a well-studied shock-bubble interaction problem that involves the collision of a shock wave in air with a circular gas bubble (cf. [14,33,44]). To setup the problem, we take a shock tube of size: (x, y)2 [0, 445] · [0, 89] mm2, and use an initial condition that consists of a planar leftward- moving Mach 1.22 shock wave in air located at x = 275 mm traveling from right to left, and a stationary gas bubble with center (x0, y0) = (225, 44.5) mm and of radius r0= 25 mm in front of the shock which contains an R22 gas inside. Here both the air and R22 are modeled as perfect gases so that we have C = c 1 (c > 1 denotes the usual ratio of specific heats) and pref= eref= 0 in the Mie–Gru¨neisen equation of state(2). For this problem, inside the R22 gas bubble the state variables are
ðq; u; v; p; c; Y Þ ¼ ð3:863 kg=m3;0; 0; 1:01325 105Pa; 1:249; 1Þ;
and outside the bubble they are
ðq; u; v; p; c; Y Þ ¼ ð1:225 kg=m3;0; 0; 1:01325 105Pa; 1:4; 0Þ and
ðq; u; v; p; c; Y Þ ¼ ð1:686 kg=m3;113:5 m=s; 0; 1:59 105Pa; 1:4; 0Þ
in the preshock and postshock regions, respectively. We note that this gives us one example that the (air–R22) material interface is accelerated by a shock wave coming from the light-fluid to the heavy-fluid region, and the early stage of the resulting wave pattern after the interaction would consist of a transmitted shock wave, an interface, and a reflected shock wave.
To solve this two-fluid problem numerically, our algorithm (tracking or non-tracking version) uses a mix- ture type of the model equations proposed in [52] as the basis for approximation. Results of a sample run obtained using each of the tracking and non-tracking algorithms with an underlying uniform 3560· 356 Cartesian grid are shown inFig. 5where a sequence of the schlieren-type images of the density in some short distances around the gas bubble is presented at nine different times t = 55, 115, 135, 187, 247, 318, 342, 417, and 1020 ls (measured relative to the time where the incident shock first hits the upstream bubble wall). From the figure, it is easy to observe the sensible improvement of the non-tracking (i.e., interface capturing) result near the interface when the volume tracking algorithm is employed in the computation. It is also interesting to see that, as far as the qualitative structure of the solution is concerned, our tracking results are in good agree- ment with the existing experimental and numerical solutions appeared in[44, Fig. 7]. The approximate loca- tion of the tracked interfaces for the volume tracking run is displayed in Fig. 6, noticing the complicated topological changes of the interfaces during the computations. We note that it would be a notorious task to do if we want to track the interface of this degree of complexity by using an alternative surface-based inter- face tracking method like the one discussed in [31], for example.
As in[44], to get a quantitative assessment of several prominent flow features, we make a diagnosis of the space–time locations of the incident shock wave, the upstream bubble wall, the downstream bubble wall, the
Fig. 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 tracking and non-tracking version of the methods is shown at nine different times t = 55, 115, 135, 187, 247, 318, 342, 417, and 1020 ls with an underlying uniform 3560· 356 Cartesian grid.
refracted shock wave, and the transmitted shock wave at some selected times, and report our findings inFig. 7 only up to time t = 250 ls, seeFig. 5for an illustration of the sample positions of the marked symbols plotted in the figure and also Fig. 12 of[44]for comparison. Here we obtain the position of the incident shock, marked by symbol ‘‘·’’, by checking the jumps in the pressure to see whether there are greater than a prescribed tol- erance, at the horizontal sections of 5 mm above the bottom boundary. For the positions of the refracted shock (marked by symbol ‘‘n’’ when time t 6 202 ls), and the two transmitted shocks (marked by symbols
‘‘
’’ and ‘‘n’’ when time t > 202 ls), we do the same thing as in the incident shock case, but the pressure jumps at the axis of symmetry are checked instead. Now, as to the positions of the upstream and downstream bubble walls which are marked by symbols ’’+’’ and ‘‘e’’, respectively, there are two cases to be considered. That is, in the volume tracking run, we record the right-most and left-most of location of the tracking interface at the axis of symmetry, while in the non-tracking run we check the volume fractions to see if there are close to the midpoint value 1/2.Fig. 5 (continued)
With these trajectories, it is a common practice to perform a linear least-squares fit of each set of points in Fig. 7separately, and take the slope of the respective line as one measure of the wave speed of interests. A comparison of the various computed velocities obtained using our tracking and non-tracking algorithms with those reported in the literature[14,44]is summarized inTable 1, where Vsrepresents the speed of the incident shock wave during the time t2 [0, 250] ls, Vui(Vuf) is the speed of the initial (final) upstream bubble wall when time t2 [0, 400] ls (t 2 [400, 1000] ls), Vdi(Vdf) is the speed of the initial (final) downstream bubble wall when time t2 [200, 400] ls (t 2 [400, 1000] ls), VRis the speed of the refracted shock when time t2 [0, 202] ls, and
time=55µs
air R22
time=115µs time=135µs
time=187µs time=247µs time=200µs
time=342µs time=417µs time=1020µs
Fig. 6. Approximate locations of the interfaces for the volume tracking run shown inFig. 5, where the dashed line in each subplot is the initial location of the air–R22 interface.
0 20 40 60
0 50 100 150 200 250
x (mm)
t (µs)
Tracking
0 20 40 60
0 50 100 150 200 250
x (mm)
t (µs)
Capturing
Fig. 7. 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 ‘‘e’’), the refracted shock (marked by symbol ‘‘n’’), and the transmitted shock (marked by symbols ‘‘’’ and ‘‘n’’) for the runs shown inFig. 5. These trajectories can be used to estimate the speed such as Vs, Vu, Vd, VR, VTin the order mentioned above.
VT is the speed of the transmitted shock when time t2 [202, 250] ls. From the table, with reference to the experimental velocities, it is clear that our algorithm (either tracking or non-tracking) gives a better prediction of the speeds for shock waves such as Vs, VR, and VT than the numerical results of Quirk and Karni. In the other instance to the prediction of the speeds for bubble walls, we still observe a better agreement of the down- stream case: Vdiand Vdf, but a slightly less accurate ones for the upstream case: Vuiand Vuf. It should be men- tioned that because the needle-shape like bubble wall in the downstream is not being tracked by the algorithm, seeFigs. 5 and 6, we do not find a good match of the speed Vdffor that. Nevertheless, due to the sensitivity of the interface structure on the grid resolution used in the computation as well as the uncertainty on the exper- imental error, it should be fair to say that we have come up with a reasonable set of data for the bubble veloc- ities in this application.
Finally, to see how our algorithms works for the pressure, inFig. 8we record the time history of the runs shown inFig. 5at six stations, xp= 3, 11, 27, 43, 67, and 99 mm, downstream of the R22 gas bubble along the axis of symmetry. As compared this result with the one shown in Fig. 14 of[44], for example, we again notice
Table 1
A comparison of the computed velocities obtained using our tracking and non-tracking algorithms forExample 5.1with those reported in the literature[14,44]; see the text for the definition of the notations used in the table
Velocity (m/s) Vs VR VT Vui Vuf Vdi Vdf
Experiment[14] 415 240 540 73 90 78 78
Quirk and Karni[44] 420 254 560 74 90 116 82
Our result (tracking) 411 243 538 64 87 82 60
Our result (capturing) 411 244 534 65 86 98 76
100 200 300
1 2 3 4 5
xp = 3mm
p (bar)
200 250 300 350 1
1.5 2 2.5 3
xp = 11mm
200 300 400
1 1.5 2 2.5
xp = 27mm
250 300 350 400 1
1.2 1.4 1.6 1.8 2
xp = 43mm
time (µs)
p (bar)
300 400 500
1 1.2 1.4 1.6 1.8 2
xp = 67mm
time (µs) 350 400 450 500 550
1 1.2 1.4 1.6 1.8 2
xp = 99mm
time (µs)
Fig. 8. Plot of the time history of the pressure for the runs shown inFig. 5at six stations, xp= 3, 11, 27, 43, 67, and 99 mm, downstream of the R22 gas bubble along the axis of symmetry. The dotted points are the results obtained using the volume tracking method, and the solid lines are the results with interface capturing.
the sensible agreement of the solution behavior to the problem. Note that in carrying out the test we have used the solid wall boundary condition on the top and bottom, and the non-reflecting boundary condition on the left and right sides.
Example 5.2. Our next example of tracking material interfaces is concerned with a model underwater explosion problem (cf.[12,55]). In this test, we take a rectangular domain (x, y)2 [2, 2] · [1.5, 1] m2, and consider the initial condition that is composed of a horizontal air–water interface at the y = 0 axis and a circular gas bubble in water with the center (x0, y0) = (0,0.3) m and of the radius r0= 0.12 m. Here above the air–water interface, the fluid is a perfect gas at the standard atmospheric condition,
ðq; p; c; Y Þ ¼ ð1:225 kg=m3;1:01325 105Pa; 1:4; 1Þ.
Below the air–water interface, in region inside the gas bubble the fluid is a perfect gas also with the state variables
ðq; p; c; Y Þ ¼ ð1250 kg=m3;109Pa; 1:4; 1Þ;
and in region outside the gas bubble the fluid is water modeled by the simplified form of the Mie–Gru¨neisen equation of state(2)as
p¼ ðc 1Þqe cB; ð19Þ
where B is a pressure-like constant, and the state variables for that are ðq; p; c; B; Y Þ ¼ ð103kg=m3;1:01325 105Pa; 4:4; 6 108Pa; 0Þ.
Note that, initially, all the fluid components of interests is in a stationary and equilibrium state, but due to the pressure difference between the fluids, breaking of the underwater bubble occurs instantaneously, yielding a circularly outward-going shock wave in water, an inward-going rarefaction wave in gas, and an interface lying in between that separates the gas and the water. Soon after this shock wave is diffracted through the nearby air–water interface, it is known (cf.[12,55]) 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 con- tinue rising upward, causing the subsequent deformation of the horizontal air-water interface.
Fig. 9shows the schlieren-type images of the density and pressure at four different times t = 0.2, 0.4, 0.8, and 1.2 ms obtained using our algorithms (tracking and non-tracking) with a 400· 250 grid. From the density plot, we clearly observe the improvement on the use of the volume tracking method over the non-tracking 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.
As far as the global structure of the solution is concerned, our result agrees quite well with the result shown in [12], for instance, where a more sophisticated surface tracking method is used in the computation. In Fig. 10, we present the location of the tracked interfaces for the volume tracking run shown in9. Note that here only the circular bubble interface is tracked initially, and the horizontal air–water interface is tracked in the algorithm only in time when its location deviates from the initial equilibrium position y = 0 axis which is aligned with the underlying grid. The cross-section of the density and pressure for the same run along line x = 0 is drawn in Fig. 11, giving some quantitative information about the differences of the tracking and non-tracking results at the selected times. The boundary conditions employed in the current computations are the solid wall boundary on the bottom, and the non-reflecting boundary on the remaining sides.
Example 5.3. To show how our algorithm works on tracking material interfaces with complex equation of state, we consider the simulation of the interaction of a shock wave in molybdenum with a block of encapsulated MORB (Mid-Ocean Ridge Basalt) liquid (cf.[37,54]). Here the computational region we take is a unit square domain, and we use an initial condition that is composed of a planar rightward-moving Mach 1.163 shock wave in molybdenum at x = 0.3 m traveling from left to right, and a rectangular region [0.4, 0.7]· [0, 0.5] m2of MORB liquid in front of the shock. For this problem, inside the region of the MORB liquid, we have the state variables
ðq; u; v; p; Y Þ ¼ ð2260 kg=m3;0; 0; 0; 1Þ;
Fig. 9. Numerical results for the simulation of a model underwater explosion problem. Schlieren-type images are shown for (a) the density, and (b) the pressure at four different times t = 0.2, 0.4, 0.8, and 1.2 ms obtained using both the tracking and capturing version of the methods with a 400· 250 grid.
time=0.2ms
air
water
time=0.4ms
air
water
time=0.8ms
air
water
time=1.2ms
air
water
Fig. 10. Approximate locations of the tracked interface for the run shown inFig. 9, where the dashed line in each subplot is the initial location of the horizontal and circular bubble interfaces at time t = 0.
and outside the MORB, we have the state variables in the preshock region ðq; u; v; p; Y Þ ¼ ð9961 kg=m3;0; 0; 0; 0Þ
and the state variables in the postshock region
ðq; u; v; p; Y Þ ¼ ð11; 042 kg=m3;543 m=s; 0; 3 1010Pa; 0Þ.
As in[34,37,54], to model the thermodynamic behavior of solid media such as MORB and molybdenum, we use(2)together with the explicit expressions for C, pref, and erefas follows:
CðqÞ ¼ C0
q0 q
a
;
prefðqÞ ¼ p0þ
q0c201qq0 1 s 1 qq0
h i2;
erefðqÞ ¼ e0þ 1
2q0 1q0 q
p0þ prefðqÞ
ð Þ.
ð20Þ
We note that in above C0= c0 1 represents a reference Mie–Gru¨neisen coefficient at q = q0, a2 [0, 1] is a dimensionless parameter, c0 denotes the zero-pressure isentropic speed of sound, and s is a dimensionless parameter which is related to the pressure derivative of the isentropic bulk modulus KS= q(op/oq)jS by (oKS/op)jS= 4s 1 (cf.[49]). With the zero initial reference state for p0and e0, typical sets of material quan- tities for the molybdenum are
ðq0; c0; s;C0;aÞ ¼ ð9961 kg=m3;4770 m=s; 1:43; 2:56; 1Þ;
and for the MORB liquid are
ðq0; c0; s;C0;aÞ ¼ ð2660 kg=m3;2100 m=s; 1:68; 1:18; 1Þ.
Fig. 11. Cross-sectional plots of the results for the run shown inFig. 9along x = 0, where the solid lines are the results obtained using the interface capturing method, and the dotted points are the tracking results.
