A γ-based Model Interface Capturing Method to Near-field Underwater Explosion(UNDEX)Simulation

(China Ship Scientific Research Center,Wuxi 214082,China)

A γ-based Model Interface Capturing Method to Near-field Underwater Explosion(UNDEX)Simulation

YU Jun,PAN Jian-qiang,WANG Hai-kun,MAO Hai-bin

(China Ship Scientific Research Center,Wuxi 214082,China)

The Simulation of near-field Underwater Explosion(UNDEX)was heightened interest with the development of underwater explosion research.The primary concerns that arise in the simulation of near-field underwater explosion are multi-component flows evolvement and interface capturing. In this paper,a robust numerical method for simulation of near-field underwater explosion was proposed,which combine five-order WENO construction and HLLC approximate Riemann solver,and the multi-component flow interface capturing method is based on the γ-model.This scheme was tested by 1D shock tube problem and 2D axis-symmetric problem,which show good agreement with theory results or experiment data.Then this method was applied to simulate some phenomenon in nearfield underwater explosion such as propagation of shock wave and blast near a free surface,which all show approval results.This method could provide some enlightenment in the research of cavitation and water-jet effect in near-field underwater explosion for future.

compressible multi-component flows;Interface capturing; near-field underwater explosion(UNDEX);γ-based model; WENO reconstruction;HLLC Riemann solver

0 Introduction

An Underwater Explosion(UNDEX)usually contains complicated sequence of the events that include detonation waves,shock waves,fluid-fluid interactions of detonation-produced explosive gas and water,fluid-structure interactions and cavitations[1-4].These events occur over a wide range of space and time scales.Due to its complex physics,different simulation methods for near/far-field and early/later-time are used in practice.The early-time interaction in UNDEX is mainly shock wave propagation in water,and the later-time interaction in UNDEX is mainly explosion bubble evolution.The far-field UNDEX event often involves the global response of the structure such as whipping(i.e.ship hull girder vibration[1,3]).On the assumption of inviscid and incompressible flow,the fluid in far-field UNDEX is commonly assumed to be acoustic medium[5-6].But this assumption is not valid for the near-field UNDEX which contains shock-wave propagation and strong shock-bubble interaction.The UNDEX fluid flow involves multiple fluids with different physical and thermodynamic properties across different materialinterface.In order to capture the evolvement between the blasé gas and surround water,different EOSs are applied to each side across the material interface during the simulation.So the method of solving multi-component flows with shock wave and gas-water interface evolvement becomes much important.

In early algorithms for computing compressible multi-component flows,the interface position was usually captured by the mass fraction[7],volume of fluid(VOF)[8]and level set function[9],which evolved according to an advection equation coupled to the system equation(i.e. Euler equation).The system equation was solved by first-or second-order accurate reconstructions with a Roe or HLL solver[10-11].But spurious oscillations easily appeared at interfaces. This was identified by Abgrall[12]and he proposed a quasi-conservative method based on the mass fraction formulation for gases.Subsequently,this method was extended to other general equation of state and multiphase flows[13].The method was mainly used to stiffened equation of state with different ratio of specific heats in multi-component flows,so it was usually called γbased model.In this paper,our goal is to simulate multi-component flow problems to avoid spurious oscillations near shockwaves and interfaces.In order to achieve the goals,we have extended the former γ-model interface capturing methods by implementing a high-order accurate WENO reconstruction and HLLC solver[14].The system equations for multi-component flows are presented in Chapter 1 and numerical discretization schemes are stated in Chapter 2.In Chapter 3,the method is validated using two-dimensional problem.Finally,we applied our method to simulate several typical near-field underwater explosion problems.

1 Equation of motion

In underwater explosion,multi-component flows are a subset of multiphase flows where different components,characterized by their respective ratio of specific heats,are immiscible; diffusive effects,surface tension and cavitation are neglected.Therefore,the governing equations employed are usually the Euler equations.

For the explosive gas flow,gases are assumed ideal and the inviscid Euler equation is employed.We can write the governing system for a 2D arrangement in a consistent form as

where the respective expressions of Q，F，G and S for flows are given as

where ρ is the fluid density,u and v are the flow velocities in z and r directions,respectively, p is the pressure and E is the total energy per unit volume and given bythe internal energy per unit mass.n is a system parameter which takes on a value of 1 or 2.Ifn=1,Eq.(1)is for a 2D planar flows problem;if n=2,it is for a 2D axis-symmetric flow problem.For 1D planar flow problem,v=0&n=1.For 1D spherical-symmetric problem,u=0& n=2.To close the system,we use the stiffened equation of state[15]

For perfect gases,γ is the ratio of specific heats and p∞=0;for water,γ=5.5,p∞=0[16].This equation has been used widely to simulation multi-component flows with shock waves propagation.

The interface between two flows is showed by the discontinuity in the properties,γ and p∞. of the different fluid components.Since material interfaces are convected by the flow,the properties must obey the advection equation[17].

Eqs.(1)and(4)form the system governing equations.Euler Eq.(1)is conservative,while the advection Eq.(4)is not conservative.So the system governing equations are quasi-conservative.

2 Discretization scheme

2.1 Finite volume formulation in Spatial discretization

Consider the Euler Eq.(1)in conservation form

where the flow gird sizes are given by

Integrating(5)over the control volume Ii，j,we obtain the following semi-discrete relations[18-19]：

2.2 Finite volume reconstruction

In the FV formulation,the left and right states of Riemann problem are reconstructed from the cell averages in a piecewise constant fashion.ENO reconstruction[20-21]chooses the optimal stencil to build the cell average value of a function,and the function is interpolated on either side of the cell edges.This method provides high-order accuracy and essentially non-oscillatory(ENO)behavior,while WENO reconstruction[22-23]provides a convex combination of all the candidate stencils,and constitutes an improvement on ENO schemes on many levels[24].In this paper,we choose five-order WENO reconstruction scheme.

But Multi-dimensional ENO and WENO schemes can not easily be applied to two-dimensional FV reconstruction,because two one-dimensional reconstructions should be needed per grid point[17].In this situation,a Gaussian quadrature rule is used in multiple dimensions[25].

2.3 The HLLC approximate Riemann solver

As said in Section 2.1,in FV formulation,the physic fluxes are built on the approximate Riemann solver,which may use high-order reconstruction values.The approximate Riemann solver proposed by Harten Lax and van Leer(HLL)[11]requires estimates for the fastest signal velocities emerging from the initial discontinuity at the interface.In 1992,the HLLC Riemann solver was proposed,it was a modification of the HLL scheme,whereby the missing contact and shear waves in the Euler equations are restored[26].Meanwhile,the HLLC resolves dis continuities sharply,and isolates shockwaves and contacts exactly[27].In the condition of proper left and right states around the contact boundary,the HLLC solvers could preserve positivity.At the same time,WENO reconstruction is a convex function of the candidate stencils[22].So that the HLLC solver that combining the left and right states could keep this method positivity.

As the HLLC scheme is a modification of HLL scheme described in previous section,the HLLC scheme added the slowest and fastest signal speed sLand sRwith a middlewave speed s*is shown in Fig.1.The HLLC approximate Riemann solver is given as follows[14]：

Fig.1 HLLC approximate Riemann solver

where D=L or R.The wave speeds are given by

And intermediate wave speed is computed according to[27]：

2.4 Time discretization

As mentioned before,the space discretization with five-order WENO reconstruction,in order to retain uniformly high order of time accuracy the solution is advanced in time by means of TVD Runge-Kutta methods[28].Usually,the following third order three-stage method is used：

Meanwhile,the explicit scheme considered above requires the computation of a time step△t to be used in Eq.(15),such that stability of the numerical method is ensured.One way of choosing method is

2.4 Algorithm

We introduce the following algorithm to compute compressible multi-component flows, based on the system Eqs.(1)and(4).

(1)Build the average conservative variableand interface function

(2)Using five-order WENO scheme to reconstruct the left and right state around the cell centre

(4)Use the HLLC solver to computer the numerical flux

(5)Use the adaptation to the HLLC solver to computer the right-hand side of the advection Eq.(4).

3 Algorithm validation

3.1 One-dimensional shock-tube problems

Shock-tube test is a typical Riemann discontinuous problem for testing the schemes on shock wave positing and interface capturing.The Sod[29]and Lax[30]problems have been computed using the present scheme and produced good agreement with the analytical solution. But in this paper,we introduced a more difficult shock-tube problem with higher pressure ratio[31].We consider a stiff two-fluid one-dimensional Riemann problem with initial condition

where the specific heat ratio of air and helium are γ1=1.4 and γ2=1.6,respectively.In this problem,the right gas(Helium)is highly compressed by large pressure ratios.The domain has 800 points and WENO reconstruction is used,along with the HLLC solver.Fig.2 shows density,velocity,pressure and volume fraction profiles at t=0.015.The computed solution agrees well with the exact solution and compares favorably to previous findings[32].The wave speeds are correct and there is only a little oscillation near the interface.

Fig.2 Solution of the one-dimensional Riemann problem computed on a uniform grid with N=800.Solid lines depict the exact solution whereas(□)denote numerical results

3.2 Two-dimensional axis-symmetric problem

This test case is the collapse of an spherical helium cavity in air tube by a Mach 1.25 shock.This problem has been studied in experiment by Haas[33]due to its importance in a wide range of physical and practical phenomena.The rectangular computational domain for this problem is Ω=[0，0.325]×[0，0.044 5],and the two-dimensional symmetrical axis is y=0.0.A spherical helium bubble of 45 mm radius is placed in air with its center at(0.175，0),and the shock wave is initiated at x=0.2 at t=0. The schematic of domain is shown in Fig.3.In calculating only up half part of the domain is needed.The rectangular computational domain Ω is discretized using a 650×88 uniform grid,with reflecting boundary conditions on the top and along the centerline.

Fig.3 Schematic of the shock and helium bubble in air

The test of the initial condition iswhere the specific heat ratio of air and helium are γ1=1.4 and γ2=1.67,respectively.The flow field at different moment is shown in Fig.4.The computed solution shows good agreement between the simulation and experiment results.

Fig.4 Shock and air/helium bubble interaction.Top：shadow-photographs of Haas and Sturtevant(1987), bottom：Numerical Schlieren-type results at(a)t=20 μs,(b)t=223 μs,(c)t=350 μs,(d)t=600 μs

4 Algorithm in near-field underwater explosion(UNDEX)

4.1 Explosion shock wave loading

Underwater explosion is a chemical reaction in a substance which converts the original material into a gas at very high temperature and pressure,this process occurring with extreme rapidity and evolving a great deal of heat.Generally speaking,the majority of scattered wave shock wave will not play important role after the interaction of incident wave and structure in far-field of underwater explosion.But this is not valid in near-field of underwater explosion, because when the reflected wave from the structure enters the bubble,refraction wave generates at the fluid interface and plays an important role on the structure.So that in near-field underwater explosion,the propagation of the wave between the structure and the bubble must be considered.

Next,we consider the UNDEX of an 1 g TNT explosive surrounded by water at a depth of 0.6 m.In our model,the TNT charge is replaced by a gas bubble with equivalent mass and internal energy for the original High Explosive(HE)material.The cavity radius is r which initial value r0=3V0/4（ ）π1/3with the approximate initial volume V0for the HE material.For 1 g TNT,the equivalent initial value r0is about 5.27 mm.The initial conditions are shown in Fig.5.

The initial condition is[13,36]

Fig.5 Schematic of the underwater explosion near a rigid wall

Fig.6 Schematic of the experimental arrangement

where the water is modeled with γ=5.5,p∞=4.92e8 Pa,and the gas with γ=1.25,p∞=0.0 Pa. As shown in Fig.5,the domain is axis-symmetric,and only the right half part is modeled.The rectangular computational domain for this model is Ω=[0,2 m]×[0,0.6 m].The computed domain Ω is discretized using a 400×800 non-uniform grid,where the minimum space step is only 0.25 mm.Fig.6 shows a sketch of the experiment arrangement used in this study.The water tank had dimensions of 1 m×1 m×1 m and is made from glass plates and filled up with tap water.The target plate has dimensions of 0.4 m×0.4 m×0.01 m,and it is used as rigid wall.1 g TNT detonator charge is chose as explosive and it is placed on the normal liner through the centre of target plate with the distance of 0.16m.A pressure gauges is taped to the surface of the target plate.The type of gauges is PVDF piezoelectric pressure gauges.

Fig.7 shows the simulation solution and experiment results of the pressure-time curve on the center surface of the target plate during the time of 0.2 ms.The maximum pressure in experiment is 56 MPa and it is 60 MPa in simulation solution.The whole curve is few bigger in simulation than in the experiment.The reason may be the difference from detonate charge,flows character,and the boundary conditions.But the trend is the same and the difference is acceptable in engineering application.

Fig.8 shows the simulation solution of density and pressure contours in the flow field at various times.The rigid wall is denoted by the grey region on the left of the domain.In frame 1,the compressible wave expands rapidly outwards and rarefaction ware expands inwards.The rarefaction wave interacts with the explosion gas and creates a low pressure zone in the centre of the bubble.As shown in frame 2,shock wave reaches the rigid wave and reflects compressible wave.When the reflect wave overlaps the shock wave,some high pressure zone will appear especially round the centre of the target wall.In frame 3,the reflect wave from the wall has impacted on the bubble and reflected with strong rarefaction waves(Prandtl-Meyer rarefaction waves)while a transmitted shock is propagating into the bubble.Due to the much lower acoustic impedance of the gas medium,the transmitted shock is much weaker than the reflected Prandtl-Meyer rarefaction waves.The Prandtl-Meyer rarefaction waves are propagating in the opposite directions along the water-gas interface and cause the water pressure below the bubble surface to drop,as seen in frame 4.

Fig.7 Pressure-time curve of underwater explosion near a rigid wall.Solid line depicts the experiment result whereas(○)denote simulation solution

Fig.8 Density(top)and pressure(bottom)contours for shock wave propagation near a rigid wall at 0.05 ms,0.1 ms,0.15 ms and 0.2 ms

4.2 Underwater explosion near a free surface

The problem of underwater explosion near a free surface has been simulated by[31,34-35]as a typical phenomenon in multi-component flow simulation.The initial conditions are given as follows with no-dimension：a highly pressurized explosion gas cylinder of radius r=0.12 is located at the origin(0,-0.3)in water.The free surface is located at the straight-line y=0.The schematic of domain is shown in Fig.9.The initial condition is[31]：

Fig.9 Schematic of the underwaterexplosion near a free surface

For this two-dimensional planer problem,The rectangular computational domain is[-2, 2]×[-1.5,1.5],using a 600×450 uniform grid.The reflecting boundary condition is put at the bottom,while the other direction boundary conditions are non-reflecting condition.

Fig.10 Pressure(up)and interface function(down)contours for underwater explosion near a free surface at times t=0.06,0.25,0.4 and 0.5(from left to right)

The computing solution is listed in Fig.10.It depicts the flow fields at different times using Ref.[31]and this paper’s methods.In the paper[31],Shukla chose the five-equation modelas the system governing equations,and the advection variable was volume of fraction of the gas.This problem is computed by this paper’s method and the results is depicted in Fig.10.It is noted that the computed solutions are little numerical oscillations and agree well with previously published results.

5 Conclusions

We have chosen an accurate quasi-conservative scheme for simulating the compressible multi-component flows,using five-order WENO reconstruction and HLLC approximate Riemann solver.The scheme has been tested in 1D shock tube problem and 2D axis-symmetric problem.The numerical results are accurate with no numerical oscillation.These results are in good agreement with experimental results and previous numerical results obtained by more sophisticated numerical methods.Secondly,we used the method to simulate some phenomenon in near-field underwater explosion.The early shock wave loading could be captured properly. At last,we show the evolution of underwater explosion near a free surface.Multi-component flow simulation is only the first step to research near-field underwater explosion,where the physics contains complicate phenomenon such as cavitation and water-jet effect.Those factors will be considered in the next study.

[1]Cole R H.Underwater explosion[M].Dover publications,1948.

[2]Keil A H.The response of ship to underwater explosions[J].The Society of Naval Architects and Marine Engineers,1961.

[3]Liu J H.Theory and its applications of ship dynamic responses to non-contact underwater explosions[D].Wuxi：CSSRC doctor’s thesis,2002.

[4]Reid W D.The response of surface ships to underwater explosions[R].Technical Report,Department of Defense,1994.

[5]Dobratz B M,Crawford P C.LLNL explosives handbook：Properties of chemical explosives and explosive simulations[M]. Lawrence Livermore National Laboratory,University of California,Livermore,CA,31,1985.

[6]Geers T L.Transient response analysis of submerged structures[J].Finite Elements in Analysis and Design,1989,6：153-172.

[7]Naber J.A numerical solver for compressible two-fluid flow[R].Report MASE0505,Amsterdam：CWI,2005：1-48.

[8]Osher S,Fedkiw R P.Level set methods and dynamic implicit surface[M].Applied Mathematical Sciences,153,Springer, 2003.

[9]Fedkiw R P,Aslam T,Merriman B.A non-oscillatory Eulerian approach to interfaces in multi-material flows(The Ghost Fluid Method)[J].Journal of Computational Physics,1999,152：457-492.

[10]Roe P L.Approximate Riemann solvers,parameter vectors,and difference schemes[J].J Comput.Phys.,1981,43：357-372.

[11]Harten A,Lax P D,Van Leer B.On upstream differencing and Godunov-type schemes for hyperbolic conservation laws [J].SIAM Rev.,1983,25：35-61.

[12]Abgrall R.How to prevent pressure oscillations in multicomponent flow calculations：A quasi-conserbative approach[J].J Comput.Phys.,1993,125：150-160.

[13]Shyue K M.An efficient shock-capturing algorithm for compressible multicomponent probles[J].J Comput.Phys.,1998, 142：208-242.

[14]Toro E F.Riemann solvers and numerical methods for fluid dynamics[M].Springer,Berlin,2007.

[15]Harlow F.Amsden A.Fluid dynamics,Monograph LA-4700[R].Los Alamos National Laboratory,Los Alamos,NM,1971.

[16]Cocchi J P,Saurel R.CLoraud J.Treatment of interface problems with Godunov-type schemes[J].Shock Waces,1996,5：347-357.

[17]Johnsen E,Colonius T.Implementation of WENO schemes in compressible multicomponent flow problems[J].J Comput. Phys.,2006,219：715-732.

[18]Shu C W.High order ENO and WENO schemes for computational fluid dynamics,in High-Order Methods for Computational Physics[J].Edited by Barth T J and Deconinck H Lecture Notes in Computational Science and Engineering, (Springer-Verlag,Berlin/New York),1999,9：439-582.

[19]Ali Akturk.Two dimensional finite volume weighted essentially non-osciliatory Euler schemes with different flux algorithms [D].The Degree of Master’s thesis of Middle East Technical University,2005.

[20]Harten A,Engquist B,Osher S,Chakravarthy S.Uniformly high order essentially non-oscillatory schemes[J].J Comput. Phys.,1987,71：231-303.

[21]Shu C W,Osher S.Efficient implementation of essentially non-oscillatory shock capturing schemes[J].J Comput.Phys., 1988,77：439-471.

[22]Jiang G S,Shu C W.Efficient implementation of WENO schemes[J].J Comput.Phys.,1996,126：202-228.

[23]Liu X D,Osher S,Chan T.Weighted essentially non-oscillatory chemes[J].J Comput.Phys.,1994,115：200-212.

[24]Shu C W.High-order finite difference and finite volume WENO schemes and discontinuous Galerkin methods for CFD[J]. Int.J Comput.Fluid Dyn,2003,17：107-118.

[25]Titarev V A,Toro E F.Finite-volume WENO schemes for three-dimensional conservations laws[J].J Comput.Phys., 2004,201：238-260.

[26]Toro E F,Spruce M,Speares W.Restoration of the contact surface in the HLL-Riemann solver[R].Technical Report COA-9204,Department of Aerospace Science,Collegue of Aeronautics,Cranfield Institute of Technology,UK,1992.

[27]Batten P,Clarke N,Lambert C,Causon D M.On the choice of wavespeeds for the HLLC Riemann solver[J].SIAM J Sci. Comput.,1997,18：1553-1570.

[28]Shu C W.Total-variation-diminishing time discretizations.SIAM[J].J Scientific and Statistical Computing,1988,9：1073-1084.

[29]Sod G A.A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[J].J Comput. Phys.,1978,27：1-31.

[30]Lax P.Weak solutions of nonlinear hyperbolic equations and their numerical computations[J].Commun Puere Appl.Math., 1954,7：159-193.

[31]Shukla R K,Pantano C,Freund J B.An interface capturing method for the simulation of multi-phase compressible flows [J].J Comput.Phys.,2010,229：7411-7439.

[32]Abgrall R,Karni S.Computations of compressible multifluids[J].J Comput.Phys.,2001,169：594-623.

[33]Haas J F,Sturtevant B.Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities[J].J Fluid. Mech.,1987,181：41-76.

[34]Xie W F,Liu T G,Khoo B C.Application of a one-fluid model for large scale homogeneous unsteady cavitation：The modified Schmidt model[J].Comp.&Fluid.,2006,35：1177-1192.

[35]Liu.T G,Khoo B C,Yeo K S.The simulation of compressible multi-medium flowⅡ.Application to 2D underwater shock refraction[J].Computers&Fluids,2001,30：315-337.

[36]Liu G R,Liu M B.Smoothed particle hydrodynamics-a meshfree particle method[M].World Scientific Publishing Co,Singapore,2010.

基于γ界面捕捉模型的近场水下爆炸数值模拟

余 俊，潘建强，王海坤，毛海斌

（中国船舶科学研究中心，江苏 无锡 214082）

随着对水下爆炸现象研究的逐步深入，近场水下爆炸的数值模拟研究越来越受到重视。近场水下爆炸数值模拟当中的主要难点在于多相流运动以及界面捕捉方法的研究。文章提出一种较为稳定的近场水下爆炸数值模拟方法，其包含了五阶WENO重构以及HLLC近似黎曼求解器的格式，对于多相流界面的捕捉方法基于γ模型。文中的求解格式通过了一维激波管问题以及二维轴对称问题的验证，计算结果与理论解和实验结果吻合较好。在此基础上将该方法推广到近场水下爆炸过程中冲击波传播以及近水面爆炸现象的模拟，模拟结果较为满意。该文的研究方法能够为后续的近场水下爆炸过程中的空泡以及水射流的研究提供指导。

可压缩多相流；界面捕捉；近场水下爆炸；γ模型；WENO重构；HLLC求解器

O389 U661.4

A

余 俊（1984-），男，中国船舶科学研究中心工程师；

O389 U661.4

A

10.3969/j.issn.1007-7294.2015.06.003

1007-7294（2015）06-0641-13

潘建强（1971-），男，中国船舶科学研究中心研究员；

王海坤（1980-），男，中国船舶科学研究中心高级工程师；

毛海斌（1978-），男，中国船舶科学研究中心高级工程师。

Received date:2015-02-02

Biography:YU Jun(1984-),male,engineer,E-mail：feiyue617@163.com;PAN Jian-qing(1971-),male,researcher.