Accepted Manuscript Multiphase Phenomena in Underwater Explosion Nikita V. Petrov, Alexander A. Schmidt PII: DOI: Reference:
S0894-1777(14)00129-0 http://dx.doi.org/10.1016/j.expthermflusci.2014.05.008 ETF 8224
To appear in:
Experimental Thermal and Fluid Science
Received Date: Revised Date: Accepted Date:
26 March 2014 5 May 2014 6 May 2014
Please cite this article as: N.V. Petrov, A.A. Schmidt, Multiphase Phenomena in Underwater Explosion, Experimental Thermal and Fluid Science (2014), doi: http://dx.doi.org/10.1016/j.expthermflusci.2014.05.008
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Multiphase Phenomena in Underwater Explosion Nikita V. Petrov, Alexander A. Schmidt, Computational Physics Laboratory, Ioffe Physical Technical Institute, Saint Petersburg, 194021, Russia Keywords: numerical simulation, underwater explosion, blast wave, rarefaction wave, free surface, cavitation Abstract The paper is devoted to numerical investigations of processes accompanying an underwater explosion near the free surface. A model of evolution of liquid and gas separated by free surface is implemented. The model enables one to take into account cavitation in rarefaction waves propagating through the liquid. Thus, the study is directed at phenomena of propagation of compression and rarefaction waves in the liquid and in the gas, of interaction of the waves with each other and with moving free surface, and of cavitation due to pressure drop in the liquid. Governing equations are solved by a numerical method based on Godunov-type high resolution scheme. Location of the free surface is determined utilizing information on distribution of the liquid volume fraction taking into account compressibility of both media. Computations of underwater explosion near the free surface have demonstrated robustness of the proposed algorithm. media. The liquid (water the case under study) is considered a barotropic medium obeys the Tait constitutive equation, while gas is considered an ideal perfect medium. Besides, an attempt is made to introduce the stiffened gas model as the constitutive equations for both water and gas. An algorithm proposed for solution of the governing equations uses a high-resolution Godunov type numerical scheme  providing the second order accuracy for both spatial and temporal coordinates in domains of smooth gas dynamic function behaviour. Since there is no analytical solution of the Riemann problem for water, numerical procedure based on Newton iterations is implemented.
1. Introduction Although the problem of underwater explosion near free surface has been attracting attention of many researchers for a long time (see, in particular, ), many challenges still remain. Typical structure of flow, induced by energy release in water, is presented in Fig. 1, which shows Schlieren visualization of early stage of underwater detonation of 10 mg AgN3 . Compression and rarefaction waves propagating in liquid, shock wave transmitted in outer gas, and domain of cavitation behind the rarefaction wave are clearly seen in the figure. It should be noted that evolution of cavitation bubbles in changing liquid pressure is one more important process, accompanying the underwater explosion. Effects of ambient pressure and interphase mass transfer result in complex oscillating behaviour of the bubbles. For example, Fig. 2 demonstrates experimental data on time variation of radius of a bubble induced by laser energy release in glycerin . In this paper the main attention is paid to propagation of blast and rarefaction waves in liquid and gas, to their interaction with each other and with the free surface, and to inception and development of cavitation in domains of pressure drop. Analysis of these phenomena is important for both development of theory of heterogeneous media and a wide range of applications utilizing specific features of underwater explosions. One of goals of this study is development of an efficient, convenient, and flexible tool for investigation of such phenomena. Since due to cavitation the liquid becomes a two-phase bubbly medium simulation of the processes initiated by local energy release in such a medium is based on combined Euler-Lagrange approach treating the carrier phase (liquid) as continuum and the dispersed phase (cavitation bubbles) as a set of test particles. A modification of well known liquid volume fraction (VOF) method  is used to determine location of the interface between the bubbly liquid and the outer gas. To close the mathematical model it is important to choose adequate constitutive equations of the involved
Fig. 1 Schlieren visualization of underwater explosion .
Fig. 2. Time variation of radius of bubble induced by laser energy release in glycerin .
Nomenclature Greek letters
u v E t N H W U T r R V Na m s CHe cg q
component of velocity in x direction (m s-1) component of velocity in y direction (m s-1) total energy (J kg-1) time (s) number density of bubbles (m-3) bubble nucleation rate (m-3 s-1) the work of formation of a new bubble (J) energy flux (J s-1) temperature (K) radius of bubble (m) gas constant (m2 s-2 K-1) volume (m3) Avogadro constant (mol-1) molar mass (kg mol-1) entropy (J K-1) Henry constant (kg m-3 Pa-1) concentration of condensed gas (kg m-3) released energy per unit length (J m-1)
α β ρ Γ Σ σ μ η ε Subsripts
l v g b s cr eq fr
2.1. Euler stage of the algorithm In the framework of the proposed mathematical model a combined Euler-Lagrange approach is used for two-phase medium. Euler stage for the carrier phase is based on the Euler equations supplemented with terms accounting for interphase interaction. These equations can be written in the following form:
p − p∞
where p∞ = 489.115 MPa, γ=4.9 and p∞ = 0, γ = 1.4 for water and gas, respectively. "Stiffened gas" model is used in the outer gas and for both phases of the bubbly liquid in the vicinity of the liquid-gas interface. In the bulk of the liquid the barotropic Tait equation is used:
⎡⎛ ρ p = p a K ⎢⎜⎜ ⎢ ρ ⎣⎝ a
⎛ (1 − β )ρu ⎛ (1 − β )ρv ⎞ ⎞ ⎜ ⎜ ⎟ ⎟ 2 ( ) uv − 1 β ρ ⎜ (1 − β ) ρu + p ⎟ ⎜ ⎟ Fx = ⎜ ⎟ Fy = ⎜ ⎟ 2 ( ) ( ) uv u p − − + 1 β ρ 1 β ρ ⎜ ⎜ ⎟ ⎟ ⎜ (1 − β )(ρE + p )u ⎟ ⎜ (1 − β )(ρE + p )v ⎟ ⎝ ⎝ ⎠ ⎠
ρ (γ + 1)( p + p∞ ) + (γ − 1)( p0 + p∞ ) = ρ 0 (γ + 1)( p0 + p∞ ) + (γ − 1)( p + p∞ )
⎛ N b Γlv ⎞ ⎛ (1 − β )ρ ⎞ ⎜ ⎟ ⎜ ⎟ ( ) u − 1 β ρ N u Γ + Σ ⎜ ⎟ b lv x⎟ ⎜ z =⎜ = f (1 − β )ρv ⎟⎟ ⎜ N b Γlv v + Σ y ⎟ ⎜ ⎜ (1 − β )ρE ⎟ ⎜⎜ ⎟⎟ ⎝ ⎠ ⎝ N bU lv ⎠
liquid vapour gas bubbles saturation critical equilibrium front
with parameters depending on type of the medium. For rarefaction and for shock waves this equation can be reduced to:
2. Mathematical model
∂z ∂Fx ∂Fy + + =−f , ∂t ∂x ∂y
volume fraction of bubbly liquid (VOF) volume fraction of bubbles density (kg m-3) interphase mass flux (kg s-1) surface tension (N m-3) surface tension coefficient (N m-1) dynamic viscosity (Pa s) accommodation coefficient internal energy (J kg-1)
B ⎤ ⎞ ⎟⎟ − 1⎥ + p a ⎥ ⎠ ⎦
where K = 3045, B = 7.15, pa = 101325 Pa, ρa = 996.5 kg/m3 in the case of water. Application of barotropic EOS provides rapid convergence of the Riemann problem solution. Underwater explosion is simulated by high pressure pulse within a localized domain of liquid. Pulse amplitude is determined by the energy, released in this domain:
In set of equations (1) the terms responsible for interaction of the carrier phase with the dispersed one, Γlv , U lv , are calculated in the Lagrange stage of the algorithm. The terms accounting interface deformation, Σ x , Σ y , are calculated with the help of the modified VOF
⎛ 1 1⎞ pa K ε = pa ⎜ − ⎟ (1 − K ) + ρ ρ B ( − 1) ρa ⎝ a ⎠
⎡⎛ ρ ⎞ B −1 ⎤ ⎢⎜ ⎟ − 1⎥ . ⎢⎣⎝ ρ a ⎠ ⎥⎦
2.3. Modified VOF method Proposed technique extends the VOF method to the case of two compressible phases. In this case equation for the liquid volume fraction can be written as:
2.2. Equations of state (EOS) To close equations of the Euler stage two types of EOS are used. The first is the "stiffened gas" EOS  for both phases, and the second is the Tait EOS for water. The "stiffened gas" EOS utilizes approach of single equation
∂αρ l ∂αρ l u ∂αρ l v + + =0. (5) ∂t ∂x ∂y The proposed modified VOF algorithm is as follows: at each time step the governing equations are solved providing distribution of the liquid velocity and density. Extrapolation of these data into the outer gas domain enables one to solve the above equation in the vicinity of the liquid-gas interface. Obtained distribution of the liquid volume fraction provides location of liquid-gas interface. Surface tension on deformable interface is determined as a bulk force and can be written in form: Σ = σk n ,
grad α , k = div ( gradα ) . n= |grad α|
It is supposed that bubble appears on each site which is bigger than the current equilibrium vapor-gas bubble radius determined by: req =
Fig. 3. The spectrum of nucleation sites (dots present experimental data) .
where W =
3 ( p s (Tl ) − pl )
2.5. Bubble evolution As it easily can be seen from Fig. 3, initial distribution of the critical bubbles in liquid provides the bubble volume fraction which does not exceed 10-6. That is why variation of the bubble volume fraction due to dynamical processes and interphase mass transport should be necessarily taken into account. This analysis is performed at Lagrange stage of the algorithm. The bubbles are assumed to be in velocity equilibrium with the liquid. Lagrange stage is governed by a set of equations, including mass and energy conservation for the test bubble, as well as Rayleigh-Plesset equation for test bubble radius :
is the energy required for
2σ . p s − pl According to eq. (7) bubble nucleation can only occur at high stretching pressure (about 500 MPa at room temperature for water). But experiments show that nucleation starts at positive pressure. This fact can be explained by presence of impurities in liquid. The impurities serve as sites for heterogeneous nucleation of the bubbles. In the framework of the classical theory presence of nucleation sites in the liquid can be taken into account by modification of the energy required for production of the critical bubble (see, for example ). In this study heterogeneous bulk nucleation is considered to be predominant. The behavior of cavitating liquid is strongly affected by size distribution and by total number of the nucleation sites. Experimental data evidence that the size distribution function of the impurities can be written as : production of a bubble of critical radius rcr =
(V V* ) , N = N0 4 1 + (V V* )
Vapor and gas pressures in such a bubble are assumed to be equal to saturation and outer pressures, respectively.
2.4. Lagrange stage of algorithm Inception of cavitation is a result of pressure drop below the saturation pressure in rarefaction waves. In this case the key issue is nucleation of bubbles and their further evolution. Classical theory of nucleation gives the following expression for homogeneous nucleation rate in pure liquids : ⎧ W ⎫ ⎛ N ⎞ 2 2σ H = ρl ⎜ a ⎟ exp⎨− ⎬, π ⎝ m ⎠ ⎩ kTl ⎭
2σ . pv + pg − pl
d ⎡4 3 ⎤ π r ρv ⎥ = 4π r 2 Γlv dt ⎢⎣ 3 ⎦ d ⎡4 3 ⎤ π r ρ g ⎥ = 4π r 2 Γlg (10) dt ⎢⎣ 3 ⎦ d ⎛4 3 d ⎛4 3⎞ ⎞ 2 ⎜ π r ρbU b ⎟ = − pb ⎜ π r ⎟ − 4π r U lv dt ⎝ 3 dt 3 ⎠ ⎝ ⎠ 2
1 d 2 r 3 ⎛ dr ⎞ r 2 + ⎜ ⎟ = 2 ⎝ dt ⎠ ρl dt
2σ 4μl dr ⎞ ⎛ − ⎜ pb − pl − ⎟ r r dt ⎠ ⎝
Terms, accounting for interphase transport in bubbly liquid, can be written as follows (see, for example, ): - for liquid-vapor mass transfer (based on Hertz-Knudsen law):
here N0 is the total number of particles per unite volume and V* is the normalization parameter, V is the nucleation site volume (Fig. 3). Moreover, since usually there is a dissolved gas in the liquid, the bubbles, arising on these sites, contain both vapor and evolved gas.
Γlv = 4π r 2
η 2π Rv
⎛ p s (Tl ) pv ⎞ − ⎜ ⎟, ⎜ T ⎟ T l v ⎝ ⎠
the charge locates at distance of 3.3 cm below the free
- for liquid-gas mass transfer (based on Henry law):
L0 (c − C p ) , ( L0 − 1) r g He g
Γlg = − Dg
here L0 =
4π . 3β
Absence of gas diffusion on the free surface is assumed. In this case:
cg = CHe p outer .
Fig. 4. Comparison of data  with present simulations for cylindrical underwater explosion.
The interphase energy transport can be written as:
U lv = 4π r 2
(γ − 1)
Rv p s (Tl ) Tl − pv Tv 2π
surface. Photographs correspond to 14, 28, and 42 microseconds after the detonation. At conditions of the experiment cavitation starts behind the rarefaction wave propagating from free surface Results of numerical simulation based on the proposed algorithm are shown Fig. 5 (on the left). Satisfactory agreement of the predicted flow structure with the experimental one is seen. These tests demonstrate that the algorithm is adequate for simulation of flows characterized by compression and rarefaction waves, interacting with each other and with the free surface, as well as transmission of pressure pulse into the outer gas. Simulation gives the liquid pressure behind the rarefaction wave lower than the saturation pressure, this results in inception of cavitation.
The bubble volume fraction can obviously be calculated from:
β = ∑ N i π ri3 , i
here summation is performed over all test bubbles. 3. Numerical Scheme
The Euler equations (1) are solved by numerical method which is based on a high resolution Godunov type numerical scheme  with Van Leer limiter function . This method is explicit, monotonic, and shock-capturing one. It possesses the second order accuracy on the smooth solutions. Convergence of the solution is determined by the CFL condition. One of the main components of this method is solution of the Riemann problem. Such an approach was proposed in . Solution of the Riemann problem is generalized for the Tait equation of state (3) in liquid and for ”stiffened gas” equations (2) at interface. This solution was validated using acoustic impedances of the phases. Since the equations of the Lagrange stage, eq.(10), can be stiff ones, it is necessary to implement special methods of numerical solution, for example, the Adams method . 4. Results and Discussion
The proposed model of underwater explosion was validated by comparison of the computations with data presented in . For different values of the released energy, eq. (4), and sizes of domains of energy deposition Fig. 4 demonstrates dependences of pressure jump at compression wave front on reduced radius of the wave, rfr0 = rfr q , where q is released energy. One more benchmark is based on experimental study . In the simulations explosion was modelled by high pressure domain in the water. The pressure was chosen to provide the internal energy in this domain equal to the released energy. Fig. 5 (on the right) presents experimental Schlieren photographs of waves propagating in water as a result of explosion of a spherical charge of 10 mg AgN3. Center of
Fig. 5. Comparison of simulating results and experimental photographs obtained by Kleine .
Several characteristic problems relevant to underwater explosion are considered below. Evolution of the free surface and shape of a gas bubble at initial pressure 1 MPa is shown in Fig. 6. The process starts with propagation of compression and rarefaction
waves in the liquid and in the gas in bubble, respectively. In narrow domain between the bubble and the free surface, governed by eq. (5), the upward liquid velocity increases significantly, this effect and deformation of the interfaces are results of multiple interaction of the rarefaction and compression waves with the interfaces (see Fig. 7).
Fig. 6. High pressure gas bubble shape evolution.
(b) Fig. 8. Time variation of the bubble radius, vapour pressure, and vapour mass in the bubble, (a); the bubble volume fraction, the actual liquid pressure, and the bubble temperature, (b), for a typical cavitation bubble in rarefaction wave. Pressure drop rate is 8.55 1010Pa/s, initial equilibrium critical bubble radius is 6.8 1010m-3 pg = patm, pv = psat.
Fig. 7. Multiple interaction of compressible and rarefaction waves with gas bubble shape and free surface.
Fig. 8 presents simulation of evolution of a test cavitation bubble in rarefaction wave which is a result of interaction of the blast wave with the free surface. The blast wave is induced by explosion of 0.2 kg of TNT deposited on 0.2 m solid sphere locates at depth 0.5 m. The rate of pressure drop in the rarefaction wave can be estimated as 8.55 1010Pa/s. The numerical density of the bubbles is 6.8 1010 m-3. Actual ambient liquid pressure is governed by two processes: tension in the rarefaction wave and changing of the bubble volume fraction (because of inertia of the liquid rise of the bubble volume fraction leads to rise of the liquid pressure and vice versa). It is seen that initial growth of the bubble radius and volume fraction in the rarefaction wave, as well as rapid decrease of the vapour pressure and temperature in growing bubble are accompanied by permanent evaporation of the liquid into the bubble. Increase of the bubble mass leads to rise of vapour temperature and pressure while increase of the bubble volume fraction leads to rise of the ambient pressure and, as a result, in decrease of the bubble radius and volume fraction providing oscillating behavior of these values.
Fig. 9 presents flow structure induced by evolution of high pressure domain of liquid with initial diameter equal to 0.3 m, the distance of the domain center from the free surface is equal to 0.2 m, and initial pressure of the liquid in this domain is 1 MPa. The figure demonstrates flow structures at 0, 50, 95, and 105μs after the “detonation”. For 95 and 105μs profiles of the pressure and the bubble volume fraction along vertical centerline are shown. At these time moments two domains of cavitation behind the rarefaction waves can be distinguished (shown by white). The first one is in the rarefaction wave propagating to the center of high pressure domain in the liquid and the second locates behind the rarefaction wave produced in interaction of compression wave with the free surface.
induced cavitation bubbles in liquids, Appl. Physics B, 108, 345-351, (2012)  Hirt C., Amsden A., An arbitrary lagrangian-eulerian computing method for all speeds. J. of Comp. Phys. 14, 227 (1974)  Rodionov A.V. Improving the order of approximation of S.K. Godunov’s scheme. J Comp. Mathematics and Mathematical Physics 27, 1853–1860 (1987).  Harlow F., Amsden A., Fluid dynamics. LANL Monograph LA-4700 (1971)  Frenkel J. Kinetic theory of liquids. New York. Dover, (1955)  Alamgir Md., Lienhard J. H. Correlation of pressure undershoot during hot-water depressurization, Journal of Heat Transfer. 103, N 1, 52-55 (1981)  Plesset S.Milton, Prosperetti A., Bubble dynamics and cavitation, Ann. Rev. Fluid Mech. 9, 145-85 (1977)  Nigmatulin R.I., Dynamics of multiphase media, v. 1 & 2, Hemisphere, N.Y., (1990).  Van Leer B Towards the ultimate conservative difference scheme. IV. A new approach to numerical convection. Journal of Computational Physics, 23, 3, 276-299 (1977)  Godunov S. K. A finite difference method for numerical calculation of discontinuous solutions of the equation of fluid dynamics. Matematicheskii Sbornik 47, 271-306 (1959)  Oran E.S., Boris J.P. Numerical simulation of reactive flows. Elsevier Science Publ., 655 (1987).
Fig. 9. Evolution of pressure after underwater explosion near free surface taking into account inception of cavitation. (t = 0, 50, 95,105 μs) 5. Conclusions
In the present paper a combined Euler-Lagrange algorithm is proposed which allows to simulate multiphase flow structure induced by underwater explosion near the free surface. This method provides description of propagation of compression and rarefaction waves, their interaction with each other and with deformable free surface, inception of cavitation in the bulk of the liquid due to pressure drop in rarefaction waves. Evolution of cavitation bubbles in a rarefaction wave and, in particular, influence of accompanying processes on the ambient liquid pressure were considered. Computations revealed qualitative effect of the bubbles on the liquid pressure. The algorithm enables one to determine location of the liquid-gas interfaces (free surface and surface of the bubble produced by detonation). Modification of conventional VOF method was proposed for compressible media. The next step of the study is consideration of detonation processes and products. References
 Kedrinskii V.K., Hydrodynamics of Explosion: Experiments and Models (Shock Wave and High Pressure Phenomena), Springer; p. 362, (2005).  Kleine H., Tepper S, Takehara K., Etoh T.G., and Hiraki K. Cavitation induced by lowspeed underwater impact. 26th Int. Symp. ShockWaves, Proceedings 2, 895–900 (2009).  Koch S., Garen W., Hegedus F., Neu W., Reuter R., Teubner U. Time-resolved measurements of shock