An anti-diffusion based Eulerian interface-sharpening algorithm for compressible two-phase flow with cavitation

Proceedings of the 8th International Symposium on Cavitation CAV2012 - Abstract No. 198 August 14-16, 2012, Singapore

AN ANTI-DIFFUSION BASED EULERIAN INTERFACE-SHARPENING ALGORITHM

FOR COMPRESSIBLE TWO-PHASE FLOW WITH CAVITATION

Keh-Ming Shyue Department of Mathematics

National Taiwan University Taipei, Taiwan 10617 Email: [email protected]

SUMMARY

We describe a novel interface-sharpening approach for the efficient numerical simulation of compressible two-phase flow with cavitation. The algorithm uses a five-equation transport model developed by Allaireet al. (J. Comput. Phys. 181 (2002) 577-616) as the basis, and incorporates both the auxiliary anti-diffusive terms to the equations as a model for interface sharp-ening and an energy-preserving pressure cutoff procedure as a model for cavitation. A standard fractional-step method is em-ployed to solve the proposed model equations in two separate steps, yielding an easy implementation of the algorithm. Sam-ple numerical results for underwater explosions in two dimen-sions are shown to demonstrate the feasibility of the algorithm for practical problems.

INTRODUCTION

Consider an unsteady, inviscid, compressible homogeneous two-phase flow that is modeled by a five-equation transport equa-tions of the form

∂t(α1ρ1) + ∇ · (α1ρ1~u) = 0,

∂t(α2ρ2) + ∇ · (α2ρ2~u) = 0,

∂t(ρ~u) + ∇ · (ρ~u ⊗~u) + ∇p = 0,

∂tE+ ∇ · (E~u + p~u) = 0,

∂tα1+~u · ∇α1= 0,

(1)

for the principal motion of the state variables such as the partial densities, momentum, total energy, and volume fraction, respec-tively (cf. [1]). Here ρk and αk∈ [0, 1] denote in turn the kth

phasic density and volume fraction for k = 1, 2, ρ = α1ρ1+ α2ρ2

denotes the total density, ~u the vector of particle velocity, and p the mixture pressure. For simplicity, we assume that the consti-tutive law for the material of interest can be characterized by the stiffened gas equation of state

pk= (γk− 1)ρkek− γkp∞,k,

where ek is the specific internal energy, γk > 1 is the ratio of

specific heats, and p∞,k is a pressure-like constant. As usual,

E= ρe + ρ|~u|2/2 is the total energy with ρe = ∑2k=1αkρkek

rep-resenting the internal energy and |~u| being the Euclidean distance of the velocity vector. In this model, we have the saturation con-dition α1+ α2= 1 imposed for the volume fractions as well.

If we assume further the isobaric closure where p1= p2= p

in a region that contains more than one fluid component, we may derive easily the expression for the mixture pressure

p= ρ e − 2

∑

∑

and also the expression for the mixture acoustic impedance

ρ c2= 2

∑

∑

where c is the mixture speed of sound. Then it can be shown that this five-equation model is hyperbolic when each physically relevant value of the state variables of the flow are defined in the region of thermodynamic stability, see [1] for the detail.

It is clear that since the volume fraction plays an important role in this compressible two-phase flow modelling, the devise of a robust scheme that attains non-oscillatory, positivity-preserving resolution of volume fraction near the interface before and after wave interactions is highly desirable and necessary when solving this model system numerically. In this work, motivated by the recent success of the anti-diffusion Eulerian interface sharpening approach proposed by So, Hu, and Adams [2] for viscous incom-pressible two-phase flow, we are interested in a novel extension of the method from incompressible to compressible two-phase flow with and without cavitation. Our goal here is to describe the basic idea of this extension, and give a preliminary assessment of the method for some sample problems, see [3–5] for previous works in this direction.

It should be mentioned that for incompressible flow interface sharpening of some kind (cf. [6–13]) is a popular technique that is applied in an underlying advection scheme to compute sharp solution profile for the volume fraction equation arising from the modelling of numerical mixing between two different incom-pressible fluid components in a control volume. Since the result-ing volume fraction solution is used solely to the definition of the averaged material quantities that are present as coefficients in the Navier-Stokes equations, it is a common practice to decouple this volume fraction computation from the remaining part of the flow solver, yielding an efficient implementation of the method.

For compressible flow, however, because of the intrinsic cou-pling of genuinely nonlinear and linearly degenerate waves in the solution states, the development of an efficient interface sharp-ening method requires some thoughts, see [14] for an example based on artificial interface compression.

ANTI-DIFFUSION INTERFACE SHARPENING MODEL To begin with, we recall that the anti-diffusion interface sharpening model proposed by So, Hu, and Adams [2] for solving the volume-fraction transport equation

∂tα1+~u · ∇α1= 0 (4)

with given ~u and discontinuous initial α1takes the form

∂tα1+~u · ∇α1= −

1

µ∇ · (ε ∇α1) , (5)

where ε > 0 is the diffusion coefficient, and µ is a free parameter. To approximate (5) numerically, a fractional step method that consists of the following steps in each time iteration is employed:

ALGORITHM1

(1) (Advection step) Solve the homogeneous part of the model equation, i.e., (4), using a monotone scheme [15] over a time step ∆t.

(2) (Anti-diffusion step) Take the solution obtained in step 1 as the initial condition, and solve the model equation with only source term

∂τα1= −∇ · (ε∇α1)

using a simple explicit method, for example, over a time step ∆τ towards a “sharp layer”.

Here τ = t/µ is a scaled time variable. Note that numerical reg-ularization is required such as employing the MINMOD limiter to stabilize the computation of ∇α and so the flux ε∇α in step 2 of the method (cf. [16, 17]). Sample numerical results in two dimensions are present that show the feasibility of this interface-sharpening approach for passive advection problems and prob-lems arising from incompressible two-phase flow [14].

To extend this method for interfaces governed by the five-equation compressible two-phase model (1), as before (cf. [18]),

it is useful to begin with an interface-only problem where both the pressure and each component of the particle velocities are con-stant in the domain, while the other variables such as the density and the material quantities are having jumps across some inter-faces. Then we write (1) in the following non-conservative form,

∂t(α1ρ1) +~u · ∇ (α1ρ1) + α1ρ1∇ ·~u= 0, ∂t(α2ρ2) +~u · ∇ (α2ρ2) + α2ρ2∇ ·~u= 0,

∂t(ρ~u) +~u∇ · (ρ~u) + ρ~u · ∇~u + ∇p = 0,

 ∂t  ρ1 2|~u| 2  +~u∇ ·  ρ1 2|~u| 2  + ρ1 2|~u| 2 ∇ ·~u  +

[∂t(ρe) +~u · ∇ (ρe) + ρe∇ ·~u] + ∇ · (p~u) = 0,

∂tα1+~u · ∇α1= 0,

yielding easily the basic transport equations for this interface-only problem

∂t(α1ρ1) +~u · ∇ (α1ρ1) = 0,

∂t(α2ρ2) +~u · ∇ (α2ρ2) = 0,

~u (∂tρ +~u· ∇ρ) = 0,

|~u|2

2 (∂tρ +~u· ∇ρ) + [∂t(ρe) +~u · ∇ (ρe)] = 0, ∂tα1+~u · ∇α1= 0.

Having (5) in mind for the volume fraction, it should be sen-sible to assume that the anti-diffusion model for the phasic den-sities α1ρ1and α2ρ2take the form

∂t(α1ρ1) +~u · ∇ (α1ρ1) = − 1 µ∇ · (ε ∇α1ρ1) , (6) and ∂t(α2ρ2) +~u · ∇ (α2ρ2) = − 1 µ∇ · (ε ∇α2ρ2) , (7)

respectively. By summing up the above phasic-density equations, we have our model for the total density

∂tρ +~u· ∇ρ = −

1

µ∇ · (ε ∇ρ ) . (8)

With that, to ensure the velocity remains at a constant state across the interfaces the momentum equation should be modified by

~u (∂tρ +~u· ∇ρ) = −

1

µ~u ∇ · (ε∇ρ) . (9) Furthermore, to ensure the pressure retain in equilibrium, us-ing (2), (5), and (8), it is not difficult to show that the energy

equation should be modified also by

|~u|2

2 (∂tρ +~u· ∇ρ) + [∂t(ρe) +~u · ∇ (ρe)] = −1 µ  |~u|2 2 ∇ · (ε ∇ρ ) − ∆(ρ e)∇ · (ε ∇α1)  , (10) where ∆(ρe) = ρ2e2− ρ1e1.

Since in general we are interested in shock wave problems as well, we should apply the anti-diffusion terms described above only locally near the interfaces. For this reason, it is common to introduce an interface indicator, denoted by HI, to the model so

that it takes effect near the interfaces only, and has no effect on the other genuinely nonlinear shock and rarefaction waves.

In summary, with the stiffened gas equation of state (2), the anti-diffusion interface-sharpening model we propose to solve compressible two-phase flow takes the form

∂t(α1ρ1) + ∇ · (α1ρ1~u) = − 1 µHIDα1ρ1, ∂t(α2ρ2) + ∇ · (α2ρ2~u) = − 1 µHIDα2ρ2, ∂t(ρ~u) + ∇ · (ρ~u ⊗~u) + ∇p = −

1 µHIDρ~u , ∂tE+ ∇ · (E~u + p~u) = − 1 µHIDE, ∂tα1+~u · ∇α1= − 1 µDα1. (11)

Here, for convenience, we use the notation Dz= ∇ · (ε∇z) to

denote the differential term for the state variable z for z = α1,

α1ρ1, α2ρ2, ρ appearing in (5), (6), (7), and (8), respectively.

Thus, without causing any confusion, we may use the notation Dρ~u= ~uDρ to denote the differential term appearing in (9) and

DE= (|~u|2/2)Dρ+ ∆(ρe)Dα1 for the differential term appearing

in (10).

To find approximate solutions of our interface sharpening model (11), we use a similar fractional-step method as described

in ALGORITHM1 where in step 1 a standard wave propagation

method developed by LeVeque [19, 20] is employed to solve the homogeneous part of the model equations, and in step 2 the model equation with the source term only is solved using a standard ex-plicit finite difference method (cf. [21]) for the anti-diffusion heat equation towards a “sharp layer” (in practice, only 1 or 2 itera-tions are enough for the interface-sharpening purpose).

Here the local interface indicator HI is considered to be a

Heaviside function in the volume-fraction jump whose value is one when this jump in the nearby state exceeds some tolerance and zero otherwise. Alternatively, we may choose HI to be a

hyper-tangent function of the volume fraction, see [14], or use the other measure based on physical quantities such as entropy, pressure, and velocity on this matter. Note that we have taken the diffusion coefficient ε to be a function of the local velocity that varies both in space and time in the computation shown below, which is unlike the case considered in [14] where only a single diffusion constant is used in the entire domain at each time.

CAVITATION MODEL

Before proceeding further, we pause to discuss a simple pres-sure cutoff procedure that can be included in (11) as a model for cavitation, see [22–25] and references therein for other state-of-the-art cavitation models. In this case, when the pressure p computed by (2) is lower than a cavitation threshold such as the saturated vapor pressure psat, we assume it to be the value of this quantity, and so the cutoff phasic internal energy may be set by (ρkek)sat = (psat + γkp∞,k)/(γk− 1) in a cavitated region for

k= 1, 2. With that, it is quite common to implement this cav-itation model as a post-processing step in a numerical method, where after each time step, we apply the usual non-conservative energy correction E_{:= Esat =} 2

∑

see [26], or alternatively the energy-preserving volume-fraction correction proposed here,

α1:= (α1)sat =

ρ e − (ρ2e2)sat

(ρ1e1)sat − (ρ2e2)sat

. (13)

Here (α1)sat can be derived directly from the α-weighted

aver-age of ρe. Note that, in the former case, the other state variables such as α1ρ1, α2ρ2, ρ~u, and α1are kept as constants during the

correction, while in the latter case, all the conservative variables in (1) are kept as constants.

It should be noted that no matter what post-processing cor-rection, i.e., (12) or (13), is applied to the algorithm, the speed of sound in the cavitation region is computed consistently with (3). In addition, it is not difficult to see that the difference in the to-tal energy between these two different corrections is equal to a constant times the difference in the volume fraction, i.e.,

∆E = [(ρ2e2)sat − (ρ1e1)sat] ∆α1,

where ∆z = z − zsat for z = E and α1. Thus, even though there

may be no significant compression or expansion of the volume fraction during some stages of the cavitation inception and termi-nation where both ∆α1and ∆E are small, the energy-preserving

approach should be a more viable one to be used in this cutoff model as opposed to the non-conservative energy approach.

NUMERICAL EXAMPLES

We now present sample results for underwater explosions in two dimensions obtained using our interface sharpening method with cavitation. Additional results that further validate the pro-posed method will be reported elsewhere.

UNDERWATER EXPLOSIONS NEAR A FLAT WALL As a first example, we consider a model underwater explo-sion problem studied by Pishevar and Amirifar [27] near a solid

wall. In this test, the initial condition consists of a stationary cir-cular explosive charge of radius 1m that is placed at the origin of a rectangular region (x1, x2) ∈ [−15, 15] × [−3, 15]m2that is full

of water. Here, for simplicity, the explosive is modeled as a per-fect gas with the state variables ρ = 1, 270 kg/m3, p = 9, 000 bar, γ = 2.0, and p∞= 0, and the water is modeled as a stiffened gas

with ρ = 1, 000 kg/m3, p = 1 bar, γ = 4.4, and p∞= 6, 000 bar.

In carrying out the computation, the saturation pressure we use in the cutoff model for cavitation is psat = 0.3619 bar, and the boundary conditions are solid wall on the bottom, and non-reflecting on the remaining sides of the domain.

Due to the presure difference, the explosion of the gas bubble leads to an outward-going shock wave in water, an inward-going rarefaction wave in gas, and a contact discontinuity lying in be-tween that separates the gas and water. As times go along, the outward-going shock wave would be reflected from the flat wall that generates shock and bubble interaction, creating a reflected rarefaction wave that lower the pressure and the incept of cavita-tion in a region between the bubble and the wall.

Figure 1 shows a sequence of pressure contours at four dif-ferent times t = 1, 3, 4, 6ms obtained using our anti-diffusion interface sharpening algorithm using a 300 × 180 grid, where in the graph we have also included contour lines indicating the loca-tions of the gas bubble as well as the region of cavitation. From the figure, we observe the general features of the solutions as described above. It is easy to see that there is no spurious oscil-lations in the pressure near the bubble interface, showing that our anti-diffusion algorithm works in a satisfactory manner for this problem.

The time histories of the pressure and the velocity at the cen-ter of the wall are shown in Fig. 2, where a fine grid solution ob-tained using the method without the anti-diffusion is included for comparison. From the figure, we observe, on the one thing, good agreement between the results with and without anti-diffusion, and on the other thing, minor effect to the pressure and velocity field due to sharpening of interface at this plotted location at the least. We note that, qualitatively our results are in good agree-ment with the one appeared in [27] using an adaptive Arbitrary Lagrangian-Eulerian method and a barotropic model for water with cavitation, see [28, 29] for similar results but are obtained using different constitutive laws and numerical methods.

UNDERWATER EXPLOSION IN A RIGID CYLINDER Our second example is concerned with an explosion of a spherical gas bubble in a cylinder that contains water. Similar to the setup considered in [22, 28, 29], we take a cylinder that is of 0.0889m in diameter and 0.2286m height. Initially, a stationary spherical explosive charge of diameter 0.03m is placed at the cen-ter of the cylinder. Here, the explosive is modeled as a perfect gas again with the state variables ρ = 1, 770 kg/m3, p = 20, 000 bar, γ = 2.0, p∞= 0, and the water is modeled as a stiffened gas with

the same parameters as used in the previous example. The sat-uration pressure we use in this test for the cavitation cutoff is psat = 0.05 bar. Since the problem is symmetric with respect to the cylindrical axis, we only take the right-half of the cylinder during the computations. The boundary conditions we used are solid wall on the left and right, and non-reflecting on the top and

t = 1ms water shock wave interface t = 3ms reflected shock ace t = 4ms cavitation ace t = 6ms cavitation cavitation ace

Figure 1. ANTI-DIFFUSION RESULTS FOR AN UNDERWATER EXPLOSION PROBLEM NEAR A FLAT SOLID WALL. A SE-QUENCE OF PRESSURE CONTOURS IS SHOWN AT FOUR DIF-FERENT TIMES 1, 3, 4, 6MS USING A 300 × 180 GRID.

0 1 2 3 4 5 6 7 8 9 10 0 1000 2000 3000 4000 5000 6000 7000 With anti−difusion Without−anti−diffusion 0 1 2 3 4 5 6 7 8 9 10 −80 −60 −40 −20 0 20 With anti−difusion Without−anti−diffusion P res su re (at m ) V el o ci ty (m /s )

cavitation cutoff cavitation collapse

Time (ms)

Figure 2. THE TIME HISTORIES OF THE PRESSURE AND VE-LOCITY AT THE CENTER OF THE SOLID WALL. THE FILLED CIRCLES IN THE GRAPH ARE THE RESULTS AT THE SELECTED

TIMES SHOWN IN FIG. 1. RESULTS OBTAINED USING THE

METHOD WITHOUT ANTI-DIFFUSION USING A FINER 600 × 360 ARE INCLUDED ALSO FOR COMPARISON.

bottom sides.

For this problem, at the early time, the basic mechanism for the creation of the cavitation region is the same as in the previ-ous example. That is, due to the shock wave reflection from the cylinder wall and the resulting shock-bubble interactions, a lower pressure region is formed between the bubble and the cylinder wall. However, due to the geometric effect, complicated wave in-teraction occurs at a later time, yielding the collapse of the origi-nal cavitation, and the formation of the others in various parts of the cylinder.

In Fig. 3, we plot a sequence of the pressure contours at four different times t = 30, 45, 60, 120µs obtained using our anti-diffusion interface sharpening algorithm using a 140 × 720 grid, where we have made the graph for the solution in the whole cylin-drical region so that the basic features of the solution mentioned above can be seen more clearly. From the figure, we again ob-serve smooth behavior of the pressure field near the bubble inter-face computed by the algorithm.

The time histories of the pressure and the velocity at the cen-ter of the cylindrical wall are shown in Fig. 4, where a fine grid solution obtained using the method without the anti-diffusion is included for comparison. From the figure, we again observe good agreement between the results with and without anti-diffusion, and not so much effect to the pressure and velocity field due to sharpening of interface. Furthermore, we notice some qualitative agreement of solutions when our results are in comparison with those one appeared in the literature, see [22, 28, 29].

Note that to solve this problem numerically, geometric source terms have been included in the model (cf. [30]), and have been handled as one of the fractional steps before the anti-diffusion step of the original algorithm.

CONCLUSIONS

We have presented a simple interface-sharpening approach based on the anti-diffusion viewpoint for the numerical resolu-tion of compressible two-phase flow with cavitaresolu-tion. Numeri-cal results for underwater explosions show sensible behavior of the anti-diffusion solutions when comparisons are made with the finer grid solutions without anti-diffusion. Ongoing work is to extend this approach further to problems with phase transitions (cf. [31, 32]).

ACKNOWLEDGMENTS

The author was supported in part by the National Science Council of Taiwan Grants #96-2115-M-002-008-MY3 and 99-2115-M-002-005-MY2.

REFERENCES

[1] 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.

[2] So, K. K., Hu, X. Y., and Adams, N. A., 2011. “Anti-diffusion method for interface steepening in two-phase in-compressible flow”. J. Comput. Phys., 230, pp. 5155–5177.

t = 30µs shock reflected rarefaction reflected shock ace t = 45µs cavitation interface t = 60µs cavitation ace t = 120µs cavitation ace

Figure 3. ANTI-DIFFUSION RESULTS FOR AN UNDERWATER EXPLOSION PROBLEM IN THE SINGLE RIGID CYLINDER. A SE-QUENCE OF PRESSURE CONTOURS IS SHOWN AT FOUR DIF-FERENT TIMES t = 30, 45, 60, AND 120µS USING A 140 × 720 GRID.

[3] Despr´es, B., and Lagoutie`re, F., 2001. “Contact discontinu-ity capturing schemes for linear advection and compressible gas dynamics”. J. Sci. Comput., 16(4), pp. 479–524. [4] Kokh, S., and Lagouti`ere, F., 2010. “An anti-diffusive

nu-merical scheme for the simulation of interfaces between compressible fluids by means of a five-equation model”. J. Comput. Phys., 229, pp. 2773–2809.

[5] Xu, Z., and Shu, C.-W., 2005. “Anti-diffusive flux correc-tions for high order finite difference WENO schemes”. J. Comput. Phys., 205, pp. 458–485.

[6] Boris, J. P., and Book, D. L., 1973. “Flux-corrected trans-port I. SHASTA, a fluid transtrans-port algorithm that works”. J. Comput. Phys., 11, pp. 38–69.

0 20 40 60 80 100 120 0 2000 4000 6000 8000 With anti−difusion Without−anti−diffusion 0 20 40 60 80 100 120 −10 0 10 20 30 40 50 60 With anti−difusion Without−anti−diffusion P res su re (at m ) V el o ci ty (m /s ) Time (µs) cavitation cutoff

cavitation collapse and reload

Figure 4. THE TIME HISTORIES OF THE PRESSURE AND VE-LOCITY AT THE CENTER OF THE CYLINRICAL WALL. THE FILLED CIRCLES IN THE GRAPH ARE THE RESULTS AT THE SELECTED TIMES SHOWN IN FIG. 3. RESULTS OBTAINED US-ING THE METHOD WITHOUT ANTI-DIFFUSION USUS-ING A FINER 280 × 1440 GRID ARE INCLUDED ALSO.

interface calculation)”. In Proc. 5th Intl. Conf. on Numer. Meth. in Fluid Dynamics, A. I. van de Vooren and P. J. Zandbergen, eds., Springer-Verlag.

[8] Ubbink, O., and Issa, R. I., 1999. “A method for captur-ing sharp fluid interfaces on arbitrary meshes”. J. Comput. Phys., 153, pp. 26–50.

[9] LeVeque, R. J., 2002. Finite Volume Methods for Hyper-bolic Problems. Cambridge University Press.

[10] Xiao, F., Honma, Y., and Kono, T., 2005. “A simple alge-braic interface capturing scheme using hyperbolic tangent function”. Int. J. Numer. Mech. Fluids, 48, pp. 1023–1040. [11] Xiao, F., Li, S., and Chen, C., 2011. “Revisit to the THINC scheme: A simple algebraic VOF algorithm”. J. Comput. Phys., 230, pp. 7086–7092.

[12] Yokoi, K., 2007. “Efficient implementation of THINC scheme: A simple and practical smoothed VOF algorithm”. J. Comput. Phys., 226, pp. 1985–2002.

[13] Youngs, D. L., 1982. “Time-dependent multi-material flow with large fluid distortion”. In Numerical Methods for Fluid Dynamics, K. W. Morton and M. J. Baines, eds., Academic Press, pp. 273–285.

[14] Shukla, R. K., Pantano, C., and Freund, J. B., 2010. “An interface capturing method for the simulation of multi-phase compressible flows”. J. Comput. Phys., 229, pp. 7411– 7439.

[15] LeVeque, R. J., 1992. Numerical Methods for Conservation Laws, 2 ed. Birkh ¨auser-Verlag.

[16] Breuß, M., Brox, T., Sonar, T., and Weicker, J., 2005. “Sta-bilized nonlinear inverse diffusion for approximating hyper-bolic PDEs”. Proceedings scale space and PDE methods in computer vision, 3459, pp. 536–547.

[17] Breuß, M., and Welk, M., 2007. “Staircasing in semidis-crete stabilized inverse linear diffusion algorithms”. J. Com-put. Appl. Math., 206, pp. 520–533.

[18] Shyue, K.-M., 1998. “An efficient shock-capturing algo-rithm for compressible multicomponent problems”. J. Com-put. Phys., 142, pp. 208–242.

[19] LeVeque, R. J., 1988. “High resolution finite volume meth-ods on arbitrary grids via wave propagation”. J. Comput. Phys., 78, pp. 36–63.

[20] LeVeque, R. J., 1997. “Wave propagation algorithms for multi-dimensional hyperbolic systems”. J. Comput. Phys., 131, pp. 327–353.

[21] LeVeque, R. J., 2007. Finite Difference Methods for Ordi-nary and Partial Differential Equations: Steady-State and Time-Dependent Problems. SIAM, Philadelphia.

[22] 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.

[23] Saurel, R., and Abgrall, R., 1999. “A multiphase Godunov method for compressible multifluid and multiphase flows”. J. Comput. Phys., 150, pp. 425–467.

[24] Saurel, R., and LeMetayer, O., 2001. “A multiphase model for compressible flows with interfaces, shocks, detonation waves, and cavitation”. J. Fluid Mech., 431, pp. 239–271. [25] 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.

[26] Deiterding, R., Cirak, F., and Mauch, S. P., 2009. “Ef-ficient fluid structure interaction simulation of viscoplas-tic and fracturing thin-shells subjected to underwater shock loading”. In Intl. Workshop on Fluid Structure Interaction, S. Hartmann, A. Meister, M. Sch ¨afer, and S. Turek, eds., University Press GmbH, pp. 65–80.

[27] Pishevar, A. R., and Amirifar, R., 2010. “An adaptive ALE method for underwater explosion simulations including cav-itation”. Shock Waves, 20, pp. 425–439.

[28] Park, J., August, 2008. “A coupled runge kutta dis-continuous galerkin-direct ghost fluid (rkdg-dgf) method to near-field early-time underwater explosion (undex) simulations”. PhD thesis, Virginia Poly-technic Institute and State University. Available at http://www.openthesis.org/documents/Runge-Kutta-Discontinuous-Galerkin-Direct-591985.html.

[29] Wardlaw, Jr., A. B., and Luton, J. A., 2000. “Fluid-structure mechanisms for close-in explosions”. Shock and Vibration, 7, pp. 265–275.

[30] Shyue, K.-M., 2010. “A high-resolution mapped grid algo-rithm for compressible multiphase flow problems”. J. Com-put. Phys., 229, pp. 8780–8801.

[31] Saurel, R., Petitpas, F., and Abgrall, R., 2008. “Modelling phase transition in metastable liquids: application to cavitat-ing and flashcavitat-ing flows”. J. Fluid. Mech., 607, pp. 313–350. [32] Zein, A., Hantke, M., and Warnecke, G., 2010. “Modeling phase transition for compressible two-phase flows applied to metastable liquids”. J. Comput. Phys., 229, pp. 2964– 2998.