Proceedings of ASME-JSME-KSME Joint Fluids Engineering Conference 2011 AJK2011-FED July 24-29, 2011, Hamamatsu, Shizuoka, Japan
AJK2011-04009
AN ADAPTIVE MOVING-MESH RELAXATION SCHEME FOR COMPRESSIBLE
TWO-PHASE BAROTROPIC FLOW WITH CAVITATION
Keh-Ming Shyue
Department of Mathematics National Taiwan University
Taipei, Taiwan 10617 Email: [email protected]
ABSTRACT
We describe a simple relaxation scheme for the efficient nu-merical resolution of compressible two-phase barotropic flow with and without cavitation on moving meshes. The algo-rithm uses a curvilinear-coordinate formulation of the relaxation model proposed by Saurelet al. (J. Comput. Phys. 228 (2009) 1678-1712) as the basis, and employs a wave-propagation based relaxed scheme to solve the model system on a mapped grid that is constructed by a conventional mesh-redistribution procedure for better solution adaptation. Sample numerical results in both one and two space dimensions are present that show the feasibil-ity of the proposed method for practical problems.
NOMENCLATURE
αk Volume fraction for the fluid phase k. Yk=ρk/ρ Mass fraction for the fluid phase k. ρk Density for the fluid phase k.
pk Pressure for the fluid phase k.
ck Speed of sound for the fluid phase k. ρ=∑2
k=1αkρk Total density. p=∑2
k=1αkpk Total pressure.
ui Particle velocity in the xi-direction.
Uj=∑Ni=1d uiJji Contravariant velocity in theξj-direction.
c Mixture speed of sound. cf Frozen speed of sound.
µ Relaxation parameter.
Qni j Approximate solution of the cell average in Ci jat time tn. ∆ξi Mesh size in theξi-direction.
∆t Time step from the current time tnto the next tn+1.
κi j= J(
C
i j) Jacobian of grid mapping for cell Ci j.Mw Total number of waves.
INTRODUCTION
Cavitation is commonly defined as a phenomenon in a liquid-flowing system when the pressure of the liquid falls suf-ficiently low in some region of the flow so that vapor bubbles are formed. The study of the dynamics of cavitation is an ac-tive research in many fields of science and engineering. Typi-cal examples in relation to various features and characteristics of cavitating flows can be found in [1–3], for example.
To compute cavitating flow numerically, one popular ap-proach among them is to use a two-phase barotropic model (cf. [4]) in that, if we ignore the physical effects such as mass transfer, surface tension, and viscosity, the Eulerian formulation of the basic conservation laws in Nd≥ 1 space dimension takes
the form ∂ ∂t αα12ρρ12 ρui + Nd
∑
j=1 ∂ ∂xj α1ρ1 uj α2ρ2uj ρuiuj+ p(ρ)δi j = 0 (1)for i= 1, 2, . . . , Nd. To close the system, the phasic pressure
pk(ρ) for k = 1, 2 is assumed to be a one-to-one function of the
density (this should be true at least locally), and so we may use the saturation conditionα1+α2= 1 directly, yielding a nonlinear
algebraic equation of the form g(p) =ρα1ρ1 1(p) +ρα2ρ2 2(p) − 1 = 0 (2)
to be solved for the pressure p, where the quantitiesα1ρ1and
α2ρ2are known a priori. With that, it is easy to find the
remain-ing flow variables such asρ1,ρ2,α1, andα2.
In this case, it is known that combining (1) with (2) gives a hyperbolic model that is viable for a class of homogeneous two-phase barotropic flow problems with and without cavita-tion. However, due to the non-monotonic behavior of the mixture sound speed (denoted by c) versus the volume fraction, 1/ρc2=
α1/ρ1c21+α2/ρ2c22, in the two-phase coexistent regions, it poses
a major difficulty to attain a suitable stability condition when the model is discretized by a diffuse-interface method explicitly.
To overcome this numerical difficulty, we are interested in a relaxation approach proposed by Saurel et al. [5] in that in addition to (1) a transport equation with a stiff relaxation source term is included in the model for the volume fraction, such asα1,
of the form ∂α1 ∂t + Nd
∑
j=1 uj ∂α1 ∂xj = µ (p1(ρ1) − p2(ρ2)) . (3)In contrast with the aforementioned conventional model that makes use of the saturation condition (2), here the equilibrium pressure p is obtained by taking the limit of infinite relaxation µ→∞to the solution of (3), yielding p= p1(ρ1) = p2(ρ2), and
so an algebraic equation for the relaxed volume fractionα1,
ˆ g(α1) = p1 α 1ρ1 α1 − p2 α 2ρ2 1−α1 = 0. (4)
It is important to note that since this relaxation model will be solved by a fractional-step method in the zero relaxation limit µ→ 0, it possesses a nice monotonic behavior of the frozen speed
of sound versus the mass fractions, c2f= Y1c21+Y2c22, and so is an
easier one to use as compared to the above conventional model for numerical approximation.
It is worthwhile to mention that the single-phase barotropic flow model devised in [6, 7] works well for isentropic cavitat-ing problems, but is not suitable for general non-cavitatcavitat-ing two-phase problems. This is unlike the fluid-mixture model pro-posed by the author [8] which works quite well for the two-phase barotropic scenario, but has experienced numerical difficulties for problems with cavitation.
Our goal in this work is to employ a state-of-the-art adap-tive moving-mesh method (cf. [9]) for the efficient numerical
resolution of compressible two-phase barotropic flow with and without cavitation in general non-rectangular domains. This is a fundamental step in our further development of the method to more complicated cavitating flows of practical importance (cf. [5, 10, 11]).
MODEL EQUATIONS IN CURVILINEAR COORDINATES We begin our discussion by introducing a coordinate map-ping from the physical domain(x1, x2, x3) in three dimensions
Nd= 3 to the computational domain (ξ1,ξ2,ξ3) via the relations
dx1= a1dξ1+ a2dξ2+ a3dξ3,
dx2= b1dξ1+ b2dξ2+ b3dξ3,
dx3= c1dξ1+ c2dξ2+ c3dξ3,
(5)
where ai, bi, cifor i= 1, 2, 3 are the metric terms of the mapping.
Then under this mapping, the relaxation model described above can be transformed into the new coordinate system as
∂ ∂t(α1ρ1) + 1 J Nd
∑
j=1 ∂ ∂ξj(α1ρ1Uj) = 0, ∂ ∂t(α2ρ2) + 1 J Nd∑
j=1 ∂ ∂ξj(α2ρ2Uj) = 0, ∂ ∂t(ρui) + 1 J Nd∑
j=1 ∂ ∂ξj (ρuiUj+ pJji) = 0, i= 1, 2, . . . , Nd, ∂α1 ∂t + 1 J Nd∑
j=1 Uj ∂α1 ∂ξj = µ (p1(ρ1) − p2(ρ2)) , (6)that is fundamental in our method on adaptive moving meshes. Here the quantities Ji jfor i, j = 1, 2, 3 are as a result of the
coor-dinate change that satisfies the following expressions:
J11 J12J13 J21 J22J23 J31 J32J33 = b2c3− b3c2a3c2− a2c3a2b3− a3b2 b3c1− b1c3a1c3− a3c1a3b1− a1b3 b1c2− b2c1a2c1− a1c2a1b2− a2b1 , (7) and the quantity J= det |∂(x1, x2, x3)/∂(ξ1,ξ2,ξ3)| is the
Jaco-bian of the mapping which can be computed by
J= 3
∑
i=1 aiJ1i= 3∑
i=1 biJ2i= 3∑
i=1 ciJ3i. (8)Note that during the initialization step, all the coordinate trans-formation variables such as ai, bi, ci, J1i, J2i, J3ifor i= 1, 2, 3,
and J would be determined and remained fixed at all time when a mapped grid is constructed by a chosen numerical grid genera-tor (cf. [12, 13]).
It is easy to see that (5) would be a two-dimensional coor-dinate mapping from(x1, x2) to (ξ1,ξ2) for any spatial location
x3in the physical domain, if we have a simplified data set where
the quantities a3, b3, c1, and c2 are all zero, and c3is equal to
one. In this instance, if we set Nd= 2 in (6) with the
coordi-nate transformation variables defined as in (7) and (8), we would have the same relaxation model in a two-dimensional curvilinear coordinate when a mapping of the form
dx1= a1dξ1+ a2dξ2,
dx2= b1dξ1+ b2dξ2,
(9)
is used in the derivation (cf. [14, 15]). Thus, without causing any confusion, we may simply use the symbol Nd as in the
Carte-sian case, to represent the number of spatial dimension in the curvilinear-coordinate formulation of equations.
For convenience, we write the relaxation model described above into a more compact expression by
∂q ∂t + 1 J Nd
∑
j=1 ∂ ∂ξjfj(q) + Bj(q) ∂q ∂ξj = µψ(q), (10)with q, fj, Bj, andψdefined in turn by
q= α1ρ1,α2ρ2,ρu1, . . . ,ρuNd,α1 T , fj= α1ρ1Uj,α2ρ2Uj,ρu1Uj+ pJj1, . . . ,ρuNdUj+ pJj,Nd, 0 T , Bj= diag (0, . . . , 0,Uj) , ψ= (0, . . . , 0, p1(ρ1) − p2(ρ2))T.
Note that in the Cartesian coordinates case where the coordinate mapping quantities a1, b2, c3are all equal to one, while the
re-maining ones are all zeros, (10) reduces to
∂q ∂t + Nd
∑
j=1 ∂ ∂xj ˘ fj(q) + ˘Bj(q) ∂q ∂xj = µψ(q), (11) with ˘ fj= α1ρ1uj,α2ρ2uj,ρu1uj+ pδ1 j, . . . ,ρuNduj+ pδ3 j, 0 T , ˘ Bj= diag (0, . . . , 0, uj) .Then it is easy to check that fjand Bjare related to ˘fjand ˘Bjvia
fj=∑Ni=1d f˘iJjiand Bj=∑Ni=1d B˘iJji, respectively.
With these notations, by assuming the proper smoothness of the solutions, the quasi-linear form of our model (10) can be written as ∂q ∂t + 1 J Nd
∑
j=1 (Aj(q) + Bj(q)) ∂ q ∂ξj = µψ(q), (12)where Aj=∂fj/∂q=∑iN=1d A˘iJjiis the Jacobian matrix of fjwith
˘
Ai=∂f˘i/∂q for i= 1, 2, . . . , Nd. If we assume further that the
thermodynamic description of the materials of interest is lim-ited by the stability requirement, it is a straightforward matter to show that any linear combination of the matrices ˘Ai+ ˘Bi for
i= 1, 2, . . . , Ndis diagonalizable with real eigenvalues and a
com-plete set of linearly independent right eigenvectors (cf. [16]). Hence, we may conclude that this relaxation model is hyperbolic. Regarding discontinuous solutions of the system, such as shock waves or contact discontinuities, we find the usual form of the Rankine-Hugoniot jump conditions across the waves (cf. [17]).
NUMERICAL METHODS ON MAPPED GRIDS
To set the groundwork for the later development of an adap-tive moving mesh method, we describe a finite volume method in wave-propagation form (cf. [12, 18, 19]) for the numerical ap-proximation of our relaxation model (without the source terms)
∂q ∂t + 1 J Nd
∑
j=1 ∂ ∂ξjfj(q) + Bj(q) ∂q ∂ξj = 0 (13)on a mapped grid. The method is based on solving one-dimensional Riemann 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.
To review the basic idea of the method, we consider the two-dimensional Nd= 2 quadrilateral grid case as illustrated in Fig. 1,
for example. In a finite volume method, the approximate value of the cell average of the solution q over the(i, j)th grid cell at a
time tncan be written as
Qni j≈ 1
M
(C
i j) Z Ci j q(x1, x2,tn) dx1dx2 = 1 κi j∆ξ1∆ξ2 Z ˆ Ci j q(ξ1,ξ2,tn) dξ1dξ2,where
C
i j and ˆC
i j denote the regions occupied by the grid cell(i, j) in the physical and computational domains, respectively,
i− 1 i− 1 i i j j j+ 1 j+ 1
C
i j ˆC
i j ξ1 ξ2 mapping ∆ξ1 ∆ξ2 logical domain physical domain ←− x1= x1(ξ1,ξ2) x2= x2(ξ1,ξ2) x1 x2 ~ni−1 2, jFigure 1. A SAMPLE QUADRILATERAL GRID IN OUR
TWO-DIMENSIONAL NUMEIRCAL METHOD.
In this setup, a fully discrete version of the first-order wave propagation method for the equations (13) is a Godunov-type scheme on a quadrilateral grid that takes the form
Qni j+1= Qni j− 1 κi j ∆t ∆ξ1
A
1+∆Qi−1 2, j+A
− 1∆Qi+1 2, j − 1 κi j ∆t ∆ξ2A
2+∆Qi, j−1 2+A
− 2∆Qi, j+1 2 . (14) HereA
1+∆Qi−1 2, j,A
− 1∆Qi+1 2, j,A
+ 2∆Qi, j−1 2, andA
− 2∆Qi, j+1 2 arethe right-, left-, up-, and down-moving fluctuations, respectively, that are entering into the grid cell. To determine these fluctua-tions, we need to solve the one-dimensional Riemann problems normal to the cell edges.
Computing Fluctuations
Considering the fluctuations
A
1±∆Qi−12, j arising from the
edge (i −1
2, j) between cells (i − 1, j) and (i, j), for example.
This amounts to solve a Cauchy problem in theξ1-direction that
consists of ∂q ∂t + 1 J ∂f1(q) ∂ξ1 + B1(q) 1 J ∂q ∂ξ1 = 0,
as for the equations and the piecewise constant data
q(ξ1,ξ2,tn) = ( Qni−1, j if ξ1< (ξ1)i−1 2 Qn i j if ξ1> (ξ1)i−1 2,
as for the initial condition at a time tn.
Let~ni−1
2, j= (ˆb2, − ˆa2)i−12, jand~ti−12, j= ( ˆa2, ˆb2)i−12, jbe the
unit normal and tangential vectors to the cell edge(i −1 2, j) in
the physical grid, where ˆai= ai/Siand ˆbi= bi/Siare the scaled
version of the metric elements aiand biin (9) with Si=
q
a2
i+ b2i
for i= 1, 2. Then to compute
A
1±∆Qi−12, j, as in [12, 18], we
perform the following steps:
(1) Transform the data Qn
i−1, j and Qni, jinto the new data ˘QL=
R
i−1 2, jQ n i−1, jand ˘QR=R
i−1 2, jQ n i, j. HereR
i−1 2, jis a rotation matrix defined byR
i−1 2, j= 1 0 0 0 0 0 1 0 0 0 0 0(ˆb2)i−1 2, j−( ˆa2)i−12, j0 0 0( ˆa2)i−1 2, j (ˆb2)i− 1 2, j 0 0 0 0 0 1 .(2) Solve Riemann problem in the “x1” direction for
∂q ∂t + ∂ ∂x1 ˘ f1(q) + ˘B1(q) ∂q ∂x1 = 0 (15)
with ˘f1 and ˘B1 defined by (11) and the Riemann data ˘QL
and ˘QR. When an approximate Riemann solver is used for
the numerical resolution, this would result in propagating discontinuities that are moving with speeds ˘λ1,m
i−1 2, j
and the jumps ˘
W
1,mi−12, j across each of them for m= 1, 2, . . . , Mw, see [8, 16] for an example.
(3) Define scaled speeds λ1,m
i−1 2, j = (S2)i−1 2, j ˘ λ1,m i−1 2, j and rotate jumps back to the Cartesian coordinates by
W
1,mi−1 2, j =
R
T i−1 2, j ˘W
1,m i−1 2, j for m= 1, 2, . . . , Mw.(4) Determine the left- and right-moving fluctuations as
A
1±∆Qi−1 2, j= Mw∑
m=1 λ1,m i−1 2, j ±W
1,m i−1 2, j .As usual, the notations for the quantitiesλ±are set byλ+=
max(λ, 0) andλ−= min (λ, 0).
In a similar manner, we may determine the up- and down-moving fluctuations at the edge(i, j −12) in the form
A
2±∆Qi, j−1 2 = Mw∑
m=1 λ2,m i, j−1 2 ±W
2,m i, j−1 2that is as the result of solving
∂q ∂t + 1 J ∂f2(q) ∂ξ2 + B2(q) 1 J ∂q ∂ξ2 = 0
High-Resolution Correction
To achieve high resolution (i.e., second-order accurate on smooth solutions, and sharp and monotone profiles on discontin-uous solutions), the speeds and the limited version of the jumps are used to construct the piecewise linear correction terms as be-fore (cf. [18]), and are added to (14) in flux difference form as
Qni j+1:= Qn+1 i j − 1 κi j ∆t ∆ξ1 e
F
1 i+1 2, j − eF
1 i−1 2, j − 1 κi j ∆t ∆ξ2 eF
2 i, j+1 2 − eF
2 i, j−1 2 .Here at the edge(i −1
2, j) the correction flux takes the form
e
F
1 i−12, j= 1 2 Mw∑
m=1 λ1i−,m12, j 1− ∆ t κi−1 2, j∆ξ1 λ1i−,m12, j ! fW
1,m i−12, j, whereκi−12, j= (κi−1, j+κi, j)/2. The quantity f
W
1,m is a lim-ited value of
W
1,m obtained by comparingW
1,m with the cor-respondingW
1,mfrom the neighboring Riemann problem to the left (ifλ1,m> 0) or to the right (ifλ1,m< 0) for m = 1, 2, . . . , Mw.
The correction flux e
F
2i, j−1 2
at the edge(i, j −12) can be defined in
a similar manner.
In addition to that, a transverse propagation of wave is also included in the method as a part of the high-resolution correction terms, see [18] for the details. With the transverse wave propa-gation, the method is typically stable as long as the time step∆t satisfies a variant of the CFL (Courant-Friedrichs-Lewy) condi-tion of the form
ν=∆t max i, j,m λ1i−,m1 2, j κip, j∆ξ1 , λ2i, j−,m1 2 κi, jp∆ξ2 ≤ 1, (16) where ip= i ifλ1i−,m1 2, j > 0 and i − 1 ifλ1,m i−1 2, j < 0; jpis defined analogously.
RELAXATION SOLVER ON MOVING MESHES
One simple way to extend the above mapped grid method to a moving grid version is to take an interpolation-based approach proposed by Tang and Tang [20]. In the current case with the relaxation model (10), in each time step, our method consists of the following steps:
(1) (Mesh redistribution step) Solve an elliptic-type mesh-redistribution equation for the new location of the mesh ver-tices, and then interpolate the numerical solution on the re-sulting grid conservatively. This step should be done in an iterative manner until convergence.
(2) (Zero relaxation step) Solve the homogeneous part of the relaxation model (10) on a new grid obtained in step 1. (3) (Pressure relaxation step) Solve the model system with only
the source terms in the infinite relaxation limit, yielding the relaxed volume fraction and so the equilibrium pressure. Since the numerical method for step 2 of the algorithm has been described before, we now focus our discussion on steps 1 and 3 below.
Mesh Redistribution Scheme
In the interpolation-based moving mesh method considered here, it is a common approach to redistribute the mesh vertices xi for i= 1, 2, . . . , Nd based on the solution of the elliptic-type
equation of the form
Nd
∑
j=1 ∂ ∂ξj D(q)∂ξj∂xi = 0 (17)with prescribed boundary conditions on the domain. Note that (17) are the Euler-Lagrange equations associated with a vari-ational problem for equidistribution adaptive meshes (cf. [21]). Here D is called a monitor function that is chosen to measure the degree of regularity of the underlying solution. In this work, it is taken as
D(q) = (1 −β)|∇q| +βk∇qk1 (1 −β)k∇qk∞+βk∇qk1
,
whereβ∈ [0, 1] is a free parameter, see [22] for more details.
Consider the two-dimensional case Nd= 2 as an example,
we may discretize (17) using a standard 5-point stencil finite-difference formula 1 ∆ξ2 1 Dk i−1 2, j Zik−1, j−Dk i−1 2, j+ D k i+1 2, j Zi jk+ Dki+1 2, j Zki+1, j + 1 ∆ξ2 2 Dk i, j−1 2 Zki, j−1−Dk i, j−1 2 + Dk i, j+1 2 Zki j+ Dk i, j+1 2 Zik, j+1= 0,
where Zkrepresents the kth iterate of the computed mesh vertices for x1 and x2. The resulting linear system can be solved by a
Figure 2. A SAMPLE NEW MESH (GRAPHED AS SOLID LINES) AF-TER MESH REDISTRIBUTION OF INITIAL MESH (DASHED LINES). THE NUMERICAL SOLUTIONS IN SHADED REGION OF THE NEW MESH, FOR EXAMPLE, SHOULD BE UPDATED CONSERVATIVELY.
We note that after each mesh redistribution iterate k, the nu-merical solutions on this mesh should be interpolated conserva-tively,i.e.,
∑
∀sM
Csk+1 Qks+1=∑
∀sM
Csk Qks.This can be done quite easily by using either a finite-volume or a geometric approach (cf. [9]).
It should be mentioned also that rather than using the above iterative procedure for mesh redistribution, for some applica-tions such as problems with free-surface or moving boundary, it may be more efficient to use a non-iterative Lagrangian approach (cf. [23]) instead over a prescribed time step∆t; this can be done at the least before the time when grid tangling occurs (cf. [24]).
Pressure Relaxation
Finally, in step 3, using the solution obtained in step 2 as the initial condition, we solve the model system with the source term of the form
∂q
∂t = µψ(q), which, in the infinite relaxation limit, leads to
ˆ g αn1+1= p1 (α1ρ1)n+1 αn+1 1 ! − p2 (α2ρ2)n+1 1−αn1+1 ! = 0
for the relaxed volume fractionαn1+1which can be solved by a standard iterative root-finding solver such as the secant method.
NUMERICAL RESULTS
We now present numerical results obtained using our mov-ing mesh method for two-phase barotropic flow problems with and without cavitation. For simplicity, we assume that the con-stitutive law for each of the fluid phases satisfies the Tait equation of state of the form
pk(ρ) = (p0k+
B
k)ρ
ρ0k
γk
−
B
k, (18)where p0k and ρ0k are the pressure and density at a reference
state,γk is the ratio of specific heats, and
B
k is a pressure-likeconstant for k= 1, 2. The set of parameters we take
through-out the tests are, for phase 1 (the liquid phase),(γ,
B
,ρ0, p0)1=7, 3000 bar, 103kg/m3, 1 bar, and for phase 2 (the gas phase),
(γ,
B
,ρ0, p0)2= 1.4, 0, 1 kg/m3, 1 bar
.
One-Dimensional Water-Vapor Cavitation
Our first example is a water-vapor cavitation problem in that inside a shock tube of one-meter length with closed ends the fluid is a homogeneous air-water mixture at the standard atmospheric condition(ρ1,ρ2,α1) = 103kg/m3, 1 kg/m3, 1 −ε
, withε=
10−2. Initially, inside the tube, there is a jump on the velocity
at x1= 1/2 m with speed u1= −100 m/s on the left and u1=
100 m/s on the right of the tube. In addition to that, since it is a closed tube, there are jumps also on the velocity at the both ends. With this condition, as times go on, two rarefaction waves are formed, causing the decrease of the pressure and the forma-tion of the cavitaforma-tion zone inside the tube. In the meantime, there are inward-moving shock waves propagating from the bound-aries, yielding the collapse of the region of cavitation due to the shock-cavitation interaction. Figure 3 shows contours of the vol-ume fraction and pressure in the x1-t plane up to time 0.8ms, see
Fig. 4 for the meshes over time. We observe reasonable resolu-tion of the results as compared to the ones shown in [6].
Two-Dimensional Underwater Explosion
Our second example concerns with a model two-phase un-derwater explosion problem (cf. [25]). In this test, we take a rectangular domain(x1, x2) ∈ [−2, 2] × [−1.5, 1]m2, and
em-ploy the initial condition that consists of a stationary horizon-tal free surface at the x2= 0 axis and a circular gas bubble
in water with the center (x0
1, x02) = (0, −0.3)m and of the
ra-dius r0= 0.12m. Here above the free surface, the fluid is the
gas phase at the standard atmospheric condition(ρ1,ρ2,α1) =
1 kg/m3, 103kg/m3, 1 −ε. with ε= 10−6. Below the free surface, in region inside the gas bubble the fluid is modeled as a perfect gas also with the state variables (ρ1,ρ2,α1) =
719.686 kg/m3, 103kg/m3, 1 −ε, and in region outside the gas
0.2 0.4 0.6 0.8 0 2 4 6 8x 10 −4 x1 ti m e Volume fraction cavitation 0.2 0.4 0.6 0.8 0 2 4 6 8x 10 −4 x1 ti m e Pressure
Figure 3. MOVING MESH RESULTS FOR A ONE-DIMENSIONAL CAV-ITATION PROBLEM. CONTOUR PLOTS OF THE VOLUME FRACTION AND PRESSURE ARE SHOWN IN THEx1-tPLANE UP TO0.8ms.
0 0.2 0.4 0.6 0.8 1 0 2 4 6 8x 10 −4 x1 ti m e
Figure 4. MOVING MESHES FOR THE RESULTS SHOWN IN FIG. 3.
1 kg/m3, 103kg/m3,ε, Here we takeε= 10−6. The boundary conditions considered here are solid wall on the left, right, and bottom sides, and non-reflecting on the top side of the domain.
Note that due to the pressure difference between the fluids at r= r0, 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 sepa-rates the gas and water. Soon after this shock wave is diffracted through the nearby flat air-water interface, it is known in the lit-erature (cf. [25]) 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. 5, where
we have performed the computation using with a 200× 125 grid.
From the density plot, we observe clearly the basic feature of the solution structure as mentioned above, and from the pressure plot, we see the smooth variation of the solution near the inter-face, without introducing any spurious oscillations. Interestingly, these barotropic flow results are qualitatively similar to the non-barotropic results as shown earlier in [24], for example. In Fig. 6, we present the mesh system for the run shown in Fig. 5; the dy-namical movement of the grids is clearly seen.
t=0.24ms Density Pressure air water t= 0.4ms t= 0.8ms t= 1.2ms
Figure 5. MOVING MESH RESULTS FOR A MODEL UNDERWATER EXPLOSION PROBLEM. DENSITY AND PRESSURE CONTOURS ARE SHOWN AT FOUR DIFFERENT TIMESt = 0.2,0.4,0.8, and1.2ms OBTAINED USING A200× 125GRID.
−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 6. MOVING MESHES FOR THE RESULTS SHOWN IN FIG. 5.
Two-Dimensional Supersonic Flow Over Bluntbody Our final example is a steady-state calculation of a super-sonic flow over a bluntbody. Initially, the fluid is a homoge-neous air-water mixture at the standard atmospheric condition
(ρ1,ρ2,α1) = 103kg/m3, 1 kg/m3, 1 −ε
, withε= 10−2in the whole domain. Supersonic air-water mixture is flowing into the domain with speed u1= 600m/s from the left to right. Numerical
com-0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Volume fraction 0 0.5 1 1.5 2 2.5 x 108 Pressure
Figure 7. STEADY-STATE RESULTS FOR A TWO-DIMENSIONAL SU-PERSONIC FLOW OVER A BLUNTBODY. PSEUDO COLORS OF THE VOLUME FRACTION AND PRESSURE ARE SHOWN.
putation. A cavitation zone near the top and bottom edge of the bluntbody is clearly seen. Note that there is a smooth transition in the pressure across the cavitation boundary without any spuri-ous oscillations. Results obtained using a moving mesh version of the code will be reported elsewhere in that a suitable boundary redistribution scheme needs to be devised for the mesh adapata-tion near a general non-rectangular boundary, see [26, 27] for a possible approach.
SUMMARY
We have presented a simple interpolation-based adaptive moving mesh method for the efficient numerical computation of compressible two-phase barotropic flow with and without cav-itation. The algorithm uses a curvilinear-coordinate formula-tion of the relaxaformula-tion model as the basis, and employs a wave-propagation based relaxed scheme to solve the model system on a mapped grid that is constructed by a conventional mesh redistri-bution procedure based on equidistriredistri-bution principle. Numerical results for water-vapor cavitation in one dimension, underwater explosion in two dimensions, and steady state two-dimensional supersonic flow over bluntbody are present. Ongoing works are to extend this approach further to problems with phase transitions for general non-barotropic multiphase flow.
ACKNOWLEDGMENT
The author was supported in part by National Science Coun-cil of Taiwan Grants #96-2115-M-002-008-MY3 and 99-2115-M-002-005-MY2.
REFERENCES
[1] Franc, J.-P., and Michel, J.-M., 2005. Fundamentals of Cavitation. Dordrecht: Springer Science+Business Media. [2] Brennen, C. E., 2005. Fundamentals of Multiphase Flow.
Cambridge University Press.
[3] Young, F. R., 1989. Cavitation. McGraw-Hill.
[4] Venkateswaran, S., Lindau, J. W., Kunz, R. F., and Merkle, C. L., 2002. “Computation of multiphase mixture flows with compressibility effects”. J. Comput. Phys., 180, pp. 54–77.
[5] Saurel, R., Petitpas, F., and Berry, R. A., 2009. “Simple and efficient relaxation methods for interfaces separating com-pressible fluids, cavitating flows and shocks in multiphase mixtures”. J. Comput. Phys., 228, pp. 1678–1712.
[6] Liu, T. G., Khoo, B. C., and Xie, W. F., 2004. “Isentropic one-fluid modelling of unsteady cavitating flow”. J. Com-put. Phys., 201, pp. 80–108.
[7] D.P. Schmidt, C.J. Rutland, M. C., 1999. “A fully com-pressible, two-dimensional model of small, high speed, cavitating nozzles”. Atomization Sprays, 9, pp. 255–276. [8] Shyue, K.-M., 2004. “A fluid-mixture type algorithm for
barotropic two-fluid flow problems”. J. Comput. Phys.,
200, pp. 718–748.
[9] Tang, T., 2005. “Moving mesh methods for computational fluid dynamics”. In Contemporary Mathematics, Vol. 383, Z.-C. Shi, Z. Chen, T. Tang, and D. Yu, eds., American Mathematical Society, Rhode Island, pp. 620–625. [10] Petitpas, F., Massoni, J., Saurel, R., Lapebie, E., and
Mu-nier, L., 2009. “Diffuse interface model for high speed cav-itating underwater systems”. Int. J. Multiphase Flow, 35, pp. 747–759.
[11] Saurel, R., Petitpas, F., and Abgrall, R., 2008. “Modelling phase transition in metastable liquids: application to cavi-tating and flashing flows”. J. Fluid. Mech., 607, pp. 313– 350.
[12] Calhoun, D. A., Helzel, C., and LeVeque, R. J., 2008. “Logically rectangular grids and finite volume methods for PDEs in circular and spherical domains”. SIAM Review,
50(4), pp. 723–752. Webpage to accompany this paper:
http://www.amath.washington.edu/∼rjl/pubs/circles/index.html.
[13] Thompson, J. F., Soni, B. K., and Weatherill, N. P., 1999. Handbook of Grid Generation. CRC Press.
[14] Hoffmann, K. A., and Chiang, S. T., 2000. Computational Fluid Dynamics, 4 ed. Engineering Education System, Wi-chita, Kansas, USA.
[15] Wesseling, P., 2001. Principles of Computational Fluid Dy-namics. Springer-Verlag.
[16] Allaire, G., Clerc, S., and Kokh, S., 2002. “A five-equation model for the simulation of interface between compressible fluids”. J. Comput. Phys., 181, pp. 577–616.
[17] Courant, R., and Friedrichs, K. O., 1948. Supersonic Flow and Shock waves. Wiley-Interscience, New York.
[18] LeVeque, R. J., 2002. Finite Volume Methods for Hyper-bolic Problems. Cambridge University Press.
[19] Shyue, K.-M., 2010. “A high-resolution mapped grid al-gorithm for compressible multiphase flow problems”. J. Comput. Phys., 229, pp. 8780–8801.
[20] Tang, H., and Tang, T., 2003. “Adaptive mesh methods for one- and two-dimensional hyperbolic conservation laws”. SIAM J. Numer. Anal., 41, pp. 487–515.
[21] Ceniceros, H. D., and Hou, T. Y., 2001. “An efficient dy-namically adaptive mesh for potentially singular solutions”. J. Comput. Phys., 172, pp. 609–639.
[22] van Dam, A., and Zegeling, P. A., 2010. “Balanced monitoring of flow phenomena in moving mesh methods”. Comm. Comput. Phys., 7, pp. 138–170.
[23] Addessio, F. L., Baumgardner, J. R., Dukowicz, J. K., John-son, N. L., Kashiwa, B. A., Rauenzahn, R. M., and Zemach, C., 1990. “CAVEAT: A computer code for fluid dynam-ics problems with large distortion and internal slip”. LA-10613-MS-REV. 1.
[24] Shyue, K.-M., 2009. “A simple unified coordinates method for compressible homogeneous two-phase flows”. In Proc. Symp. Appl. Math., Vol. 67, E. Tadmor, J.-G. Liu, and A. Tzavaras, eds., American Mathematical Society, pp. 949–958.
[25] Grove, J., and Menikoff, R., 1990. “The anomalous reflec-tion of a shock wave at a material interface”. J. Fluid Mech.,
219, pp. 313–336.
[26] Li, R., Tang, T., and Zhang, P., 2001. “Moving mesh meh-ods in multiple dimensions based on harmonic maps”. J. Comput. Phys., 170, pp. 562–588.
[27] Tang, H. Z., 2006. “A moving mesh method for the Euler calculations using a directional monitor function”. Comm. Comput. Phys., 1, pp. 656–676.
