Application of CWENO Finite-Volume Interface Capturing Schemes for Modeling Cavitation Induced by an Underwater Explosion

1. Introduction

Cavitation, which refers to the formation of vapor bubbles within a liquid due to high velocities and strong acceleration at low pressure, has diverse applications in various scientific and engineering fields. It can be beneficial in marine design engineering for reducing skin friction drag and achieving high speed for submerged projectiles, as well as in medicine for treating kidney stones and cancerous bodies. However, cavitation can also have adverse effects such as material erosion, noise, and vibrations in marine propellers and turbomachinery, that can cause significant damage. Therefore, understanding the dynamics of cavitation is crucial, but it poses challenges in numerical simulations of compressible multiphase problems due to the complex physics involved in phase change, high gradients of flow variables, numerical discretisation choices, and validation with experimental results. Numerical studies often focus on attached cavitations, as they are the primary subject of investigation in most research [1,2,3]. However, hull/bulk cavitation, which occurs during underwater explosions, is different as it evolves rapidly and unsteadily, changing dimensions before collapsing violently due to the surge in pressure, resulting in erratic and rapidly changing cavitation patterns. This makes numerical simulations of hull/bulk cavitation more challenging.

In hull cavitation, vapour bubbles are created on the surface of a submerged rigid surface, such as a projectile, ship hull, or submarine, due to rarefaction caused by high-pressure shock waves from an explosion. On the other hand, bulk cavitation occurs near a free surface between two interacting fluids due to the explosive shock waves, leading to widespread cavitation throughout the liquid. Hull cavitation is often more challenging to study using numerical methods due to its rapid and unsteady nature, while bulk cavitation is characterised by the formation of large-scale vapour bubbles throughout the liquid. Both forms of cavitation can have detrimental effects on marine structures and machinery.

In the context of hull or bulk cavitation, where regions experience extreme variations in pressure, negative pressure values are observed, which are physically unrealistic and can cause numerical codes to fail. To address this issue, two common approaches are used in diffuse interface models for simulating UNDEX (Underwater Explosion) phenomena. The first approach is to incorporate a single-fluid model, which has been used by several researchers [4,5,6,7,8,9] in similar UNDEX simulations. Various methods in this category, such as the cut-off model, vacuum model, isentropic model, Schmidt model, and modified Schmidt model, have been used to avoid negative pressure values during cavitation. The cut-off model assumes that cavitation bubbles do not form until the pressure drops below a critical value, known as the cut-off pressure, and is generally accurate in predicting peak pressure at the cavitation boundary, but may not accurately predict cavitation collapse periods. The Schmidt model is accurate when the void fraction (ratio of vapour bubble density to total fluid density) is less than 10 ^-5 [4], but has limitations in calculating the speed of sound, making it suitable for high-speed cavitating flows. The modified Schmidt model is an extension of the original Schmidt model and addresses its limitations, thereby expanding its range of applicability to larger vapor-to-liquid density ratios. The vacuum cavitation model assumes zero mass inside the cavitation bubble but is limited in capturing the complex shapes and dynamics of cavitation bubbles in multi-dimensional flows. The isentropic model, which has been used in the five-equation diffuse interface models [5,10], is robust as it considers the mixture as consisting of isentropic vapour and liquid phases, but it requires more computational resources and the predetermination of a model parameter k [4]. The cut-off model is preferred in this case as it is relatively simple to implement and does not alter peak pressure accuracy at the cavitation boundary, which is important for the UNDEX test cases considered in the studies. The one-fluid cavitation model aligns with the no-phase transition assumption of the five-equation interface capturing models used in this study.

The accuracy of the one-fluid model in predicting the details of the cavitation process, including thermodynamic aspects, in compressible two-phase flows, is limited due to its inability to account for mass and heat transfer between the liquid and vapour phases. To overcome this limitation, researchers have extended the simplified diffuse-interface models such as the four-, five-, or six-equation models by incorporating source terms in the energy and mass conservation equations to represent mass and heat transfer. This approach, known as the "two-fluid cavitation models," aims to prevent non-physical pressure occurrences in the cavitation region. The source terms are implicitly expressed using an algebraic relaxation solver for phase transition modelling, as proposed by Pelanti and Shyue [11,12]. However, this method is more complex and computationally expensive, especially for the seven-equation diffuse interface models such as the full non-equilibrium models of Baer-Nunziato [13] or Saurel-Abgrall [14]. A computational efficient and simple strategy in UNDEX is to employ a reduced diffuse interface models with phase transition, such as the four- or five-equation models, which enforce mechanical, thermal, or both equilibrium. Jun et al [15] extended the phase transition model used by Chiapolino et al. [16] for the four-equation model to study cavitation in underwater explosion test cases. Other existing diffuse interface models with phase transitions include those proposed by Pelanti and Shyue [9,11], Martelot et al. [17], Zein et al.[18], and others [14]. Ma et al. [19] utilised the Kapila [20] variant of the five-equation two-phase flow model with thermo-chemical relaxation terms to investigate liquid-vapour phase transition in cavitating flows for UNDEX test cases.

The present study focuses on hull and bulk cavitation in compressible multiphase or multi-component flows, specifically in the context of underwater explosions. Two interface capturing models, namely the ones proposed by Allaire et al. [21] and Kapila et al. [20], will be used to investigate the cavitation phenomenon in this study. The performance of these models will be evaluated using four underwater explosion (UNDEX) test cases in the open-source finite-volume UCNS3D solver [22].

To solve the interface capturing models, high-order central-weighted essentially non-oscillatory (CWENO) schemes for multi-component species within the finite volume framework, as implemented by Tsoutsanis and Dumbser[23], will be employed. The CWENO schemes offer computational advantages such as lower cost and faster computation compared to classical WENO schemes, due to the use of low-order polynomials associated with directional stencils and their reduced size compared to the high-order polynomials related to the central stencil. These advantages will be fully utilised in capturing the dynamics of cavitation in the test cases considered herein.

The paper proceeds to Section 2, where it presents the governing equations that describe the diffuse interface models of Allaire et al. [21] and Kapila et al. [20], which consist of five equations. It includes details about the additional $K\nabla ·\mathbf{v}$ term that distinguishes these models, as well as the cut-off model used to prevent non-physical pressure values during cavitation. In Section 3, the implementation of the reconstruction process of the CWENO schemes is described, including the chosen fluxes and temporal discretisation, within the finite volume framework using the diffuse interface methods mentioned earlier. Section 4 is dedicated to the study of five test cases, where the proposed method is utilised to investigate cavitation occurrence in underwater explosions near different boundaries, such as a planar solid wall, a closed cylindrical container, and a free surface. The analysis focuses on parameters such as the influence of shock loading and cavitation collapse reloading on nearby structures like ship hulls or closed tanks. Additionally, the bubble dynamics of cavitating flows induced by underwater explosions near a free surface, as well as the impact of bubble separation distance on cavitation phenomena, are analysed. The numerical results obtained are compared with analytical, experimental, and reference solutions whenever possible. A comparison between the five-equation models of Allaire et al. [21] and Kapila et al. [20] is also presented. Finally, the last section provides concluding remarks and a summary drawn from this study.

2. The Five-equations Diffuse Interface Model

The five-equations model typically consists of two continuity equations, one momentum equation and one energy equation coupled with a non-conservative advection equation of the volume fraction equation. The common five-equations models are the Allaire et al. [21], Kapila et al. [20], Murrone et al.[24], and Wackers et al. [25] models, all these are simpler formulations or reductions to lesser unknowns of the full seven (7) Baer-Nunziato [13] equations, which are often too complex to solve numerically. In this paper, we shall only be considering the Allaire and Kapila five-equations models. The difference between Allaire’s and Kapila’s model is only in the non-conservative advection equation in the volume-fraction equation. The equations are presented below :

$$\frac{\partial \left({\alpha }_1{\rho }_1\right)}{\partial t}+\nabla ·\left({\alpha }_1{\rho }_1\mathbf{v}\right)=0,$$

(1)

$$\frac{\partial \left({\alpha }_2{\rho }_2\right)}{\partial t}+\nabla ·\left({\alpha }_2{\rho }_2\mathbf{v}\right)=0,$$

(2)

$$\frac{\partial \rho \mathbf{v}}{\partial t}+\nabla ·(\rho \mathbf{v}\mathbf{v}+p\mathbf{I})=0,$$

(3)

$$\frac{\partial \rho E}{\partial t}+\nabla ·(\rho E+p)\mathbf{v}=0,$$

(4)

$$\frac{\partial {\alpha }_1}{\partial t}+\mathbf{v}·\nabla {\alpha }_1=0(\mathrm{Allaire}\mathrm{et}\mathrm{al}.),$$

(5)

$$\frac{\partial {\alpha }_1}{\partial t}+\mathbf{v}·\nabla {\alpha }_1=K\nabla ·\mathbf{v}(\mathrm{Kapila}\mathrm{et}\mathrm{al}.)$$

(6)

where the subscript $1,2$ represents the fluid’s component, $\alpha$ is the volume fraction of each component, $\rho$ is the density, $\mathbf{v}=(u,v,w)$ is the velocity, p is the pressure and E is the total energy.

This class of five-equations models (Allaire et al. [21] and Kapila et al. [20] model) are often referred to as the mechanical equilibrium models. We assume that the momentum, energy, and mass transfer between the phases reach equilibrium due to the thermodynamic difference between each component. This reasonable flow assumption involves a free surface like the one treated in this paper. Note that this is a mechanical equilibrium of the fluid; each phase remains in thermal disequilibrium. All the models in this category can be expressed into the following general form of a non-linear system of PDE in multiple space dimensions:

$$\frac{\partial \mathbf{Q}}{\partial t}+\nabla ·\mathbf{F}\left(\mathbf{Q}\right)+\mathbf{H}\left(\mathbf{Q}\right)\nabla ·\mathbf{v}=\mathbf{S}\left(\mathbf{Q}\right)$$

(7)

where Q is the vector of evolution variables (conserved and not conserved), F is a flux function, $\mathbf{v}$ is the velocity field, and H and S are non-conservative quantities.

We then define Q, F and H in the form as presented in Eq.(8) ( Allaire et al. model [21]) and Eq. (9) (Kapila et al. model [21]) for a two-phase flow. The K∇· v term in Kapila’s model [20] describes the compressibility and thermodynamic properties of the mixture region of the two-phase fluid. It is derived from the asymptotic reduction of Baer-Nunziato’s seven-equation model[13,26,27]. From Eq.(5) and Eq.(6), $K\nabla ·\mathbf{v}$ term differentiates Kapila’s model et al.[20] from Allaire’s et al. [21].

The Allaire’s et al. 5-equation model presented in a vector form:

$$Q=\left[\begin{array}{c}{\alpha }_1{\rho }_1\\ {\alpha }_2{\rho }_2\\ \rho \mathbf{v}\\ \rho E\\ {\alpha }_1\end{array}\right],F=\left[\begin{array}{c}{\alpha }_1{\rho }_1\mathbf{v}\\ {\alpha }_2{\rho }_2\mathbf{v}\\ \rho {\mathbf{v}}^2+p\\ \mathbf{v}(\rho E+p)\\ {\alpha }_1\mathbf{v}\end{array}\right]H=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ -{\alpha }_1\end{array}\right].$$

(8)

The Kapila’s et al. 5-equation model presented in a vector form:

$$Q=\left[\begin{array}{c}{\alpha }_1{\rho }_1\\ {\alpha }_2{\rho }_2\\ \rho \mathbf{v}\\ \rho E\\ {\alpha }_1\end{array}\right],F=\left[\begin{array}{c}{\alpha }_1{\rho }_1\mathbf{v}\\ {\alpha }_2{\rho }_2\mathbf{v}\\ \rho {\mathbf{v}}^2+p\\ \mathbf{v}(\rho E+p)\\ {\alpha }_1\mathbf{v}\end{array}\right]H=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ -{\alpha }_1-\mathbf{k}\end{array}\right].$$

(9)

For each fluid’s component i, the mixture is assumed to be in mechanical equilibrium as stated previously, the velocity and pressure of each component is $v=v_1\equiv v_2$ and $p=p_1\equiv p_2$ (i.e. we maintain a single velocity and single pressure for the mixture).

The total density, momentum, the kinetic and internal energy of the mixture is stated below:

$$\rho ={\left(\alpha \rho \right)}_1+{\left(\alpha \rho \right)}_2=\sum _{i=1}^N{\left(\alpha \rho \right)}_i:\mathrm{mixture}\mathrm{density}$$

(10)

$$\rho \mathbf{v}=({\left(\alpha \rho \right)}_1+{\left(\alpha \rho \right)}_2))\mathbf{v}=\sum _{i=1}^N{\left(\alpha \rho \right)}_i\mathbf{v}:\mathrm{mixture}\mathrm{momentum}$$

(11)

$$\rho e^{\mathrm{K}}=\frac{{\rho }_i{\mathbf{v}}^2}{2}:\mathrm{mixture}\mathrm{kinetic}\mathrm{energy}.$$

(12)

$$\rho e^{\mathrm{I}}={\left(\alpha \rho \right)}_1{e}_1+{\left(\alpha \rho \right)}_2{e}_2=\sum _{i=1}^N{\left(\alpha \rho e\right)}_i:\mathrm{mixture}\mathrm{internal}\mathrm{energy}.$$

(13)

So the total energy can be expressed as:

$$\rho E^{\mathrm{T}}=({\left(\alpha \rho \right)}_1{e}_1+{\left(\alpha \rho \right)}_2{e}_2+\frac{{\rho }_i{\mathbf{v}}^2}{2})$$

(14)

To determine the internal energy and put closure to the equations, it is convenient to use the stiffened EOS for both fluids components:

$$p_i=\left({\gamma }_i-1\right){\rho }_i{e}_i-{\gamma }_i{p}_{c,i},i=1,2$$

(15)

Using the stiffened EOS, the internal energy can now be expressed as:

$$\rho e_i=\frac{p+{\gamma }_i{p}_c,i}{{\gamma }_i-1}$$

(16)

and the total energy below:

$$\rho E=\sum _{i=1}^N{\alpha }_i(\frac{p+{\gamma }_i{p}_c,i}{{\gamma }_i-1}+\frac{{\rho }_i.v^2}{2})$$

(17)

From the five-equations, as shown in Eq. (5) and Eq. (6), the advection equation of the volume fraction is in a non-conservative form which will present some difficulties when solving our five-model equations. We simplify it into a conservative term using Johnsen and Colonius [28] approach.

First, let f represent a vector in three-dimension and is defined as, $f={\alpha }_1\mathbf{v}={\left[{\alpha }_{1}u,{\alpha }_{1}v,{\alpha }_{1}w\right]}^{\mathrm{T}}.$

Divergence of f then is:

$$\begin{array}{cc}\hfill \nabla ·f& =\frac{\partial }{\partial x}\left({\alpha }_{1}u\right)+\frac{\partial }{\partial y}\left({\alpha }_{1}v\right)+\frac{\partial }{\partial z}\left({\alpha }_{1}w\right)\hfill \\ \\ & ={\alpha }_1\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\right)+\left(u\frac{\partial {\alpha }_1}{\partial x}+v\frac{\partial {\alpha }_1}{\partial y}+w\frac{\partial {\alpha }_1}{\partial z}\right)\hfill \\ \\ & ={\alpha }_1\nabla ·\mathbf{v}+\mathbf{v}·\nabla {\alpha }_1\hfill \end{array}$$

And then,

$$\mathbf{v}·\nabla {\alpha }_1=\nabla ·\mathbf{f}-{\alpha }_1\nabla ·\mathbf{v}.$$

Substituting $\mathbf{v}·\nabla {\alpha }_1$ into the non-conservative advection equation of the volume fraction in Eq. (5) and Eq. (6) to obtain the quasi-conservative form, the equation for the advection of the volume fraction can be written as Eq. (18) and Eq. (19). For Allaire’s et al:

$$\frac{\partial {\alpha }_1}{\partial t}+\nabla ·f=\left({\alpha }_1\right)\nabla ·v$$

(18)

For Kapila’s et al.

$$\frac{\partial {\alpha }_1}{\partial t}+\nabla ·f=\left({\alpha }_1+K\right)\nabla ·v$$

(19)

The K function in Kapila’s et al. can be determined from the Eq. (20) where the mixture’s speed of sound is gotten using the Wood’s speed of sound as seen in Eq. (22)

$$K={\alpha }_1{\alpha }_2\left(\frac{1}{{\rho }_1{c}_1^2}-\frac{1}{{\rho }_2{c}_2^2}\right)\rho c^2$$

(20)

where $c_1$, $c_2$ and $\rho$ are given as :

$$c_1=\sqrt{\frac{{\gamma }_1}{{\rho }_1}\left(p+p_{c,1}\right)},c_2=\sqrt{\frac{{\gamma }_2}{{\rho }_2}\left(p+p_{c,2}\right)}$$

(21)

$$\frac{1}{\rho c^2}=\frac{{\alpha }_1}{{\rho }_1{c}_1^2}+\frac{{\alpha }_2}{{\rho }_2{c}_2^2}=\sum _{i=1}^N\frac{{\alpha }_i}{{\rho }_i{c}_i^2}$$

(22)

The “mixture-mixture” speed of sound derived from the mixture stiffened-gas SG-EoS is used for the Allaire et al. five -equations [21] and is given as follows:

$$c=\sqrt{\frac{\gamma (p+p^{\infty })}{\rho }}$$

(23)

2.1. Cut-off Methods

The stiffened-gas equation of state (SG-EOS) introduced by Harlow and Amsden [29] has been commonly used by various authors [28,30,31,32,33] to achieve thermodynamic closure for multiphase compressible flow, including flows involving cavitation. However, in cases where the pressure difference between the two-phase medium is considerably high, such as in underwater explosions, the use of stiffened gas EOS alone may not be sufficient. This is because cavitation induced by the high-pressure surge is very fast, unsteady, and rapidly evolves into various dimensions before collapsing violently, resulting in non-physical negative pressures when using only the stiffened gas EOS. Similar findings have been reported by other authors in their numerical simulations.

To improve the accuracy of simulations for cavitating flows, modifications are needed. One common approach is to incorporate an additional single-fluid model to circumvent the issue, such as the cut-off model, isentropic model, or Schmidt model, which have been used by authors in previous studies. Among these models, the cut-off model is often preferred due to its simplicity of implementation and preservation of peak pressure accuracy during cavitation collapse, which is important in underwater explosion (UNDEX) test cases. Additionally, these models are pure phase models that do not consider phase exchange, as noted in previous studies [4,34]. Further details of these models can be found in the references cited. In the cut-off model, once the pressure in the cavitation region drops below a certain fixed saturation pressure p_sat, the non-physical pressure computed assumes a new value which is the value set for the p_sat.

$$p_i=\left\{\begin{array}{cc}SG-EOS,\hfill & \mathrm{if}{p}_i>p_{\mathrm{sat}}\hfill \\ p_{\mathrm{sat}},\hfill & \mathrm{if}{p}_i<p_{\mathrm{sat}}\hfill \end{array}\right.$$

The physical saturation pressure (${\mathrm{p}}_{\mathrm{sat}}$ ) is chosen to be a comparatively small value in relation to the significantly higher ambient pressure exerted by the fluid, and $p_i$ refers to the static pressure at cell i. To avoid the negative physical pressure, most authors [4,5,35,36] assume a value between 2000-5000 Pa.

The mixture’s internal energy can be expressed as :

$${\rho }_i{e}_i^{\mathrm{I}}sat=\left(\frac{p_{sat}+{\gamma }_i{p}_{c,i}}{{\gamma }_i-1}\right),$$

(24)

and the total energy is final given by:

$$\rho e_{\mathrm{sat}}=\sum _{i=1}^N{\alpha }_i(\frac{p_{sat}+{\gamma }_i{p}_{c,i}}{{\gamma }_i-1}+\frac{{\rho }_i.v^2}{2}).$$

(25)

3. Numerical model

3.1. Finite volume framework

Lets consider a 3D domain $\Omega$ consisting of conforming tetrahedral, hexahedral, prism, and pyramid, triangles, or quadrilateral cells each one of them indexed by a unique mono-index i, and the governing equations of the five-equations model written in vector form as follows:

$$\frac{\partial }{\partial t}\underset{{\Omega }_i}{\int }\mathbf{Q}d\Omega +\underset{\partial {\Omega }_i}{\oint }{\mathbf{F}}_n(\left(Q\right)dS=\underset{{\Omega }_i}{\int }\mathbf{G}d\Omega$$

(26)

where $\mathbf{Q}=\mathbf{Q}(\mathbf{x},t)$ is the vector of conserved variables and the volume fraction of one species, ${\mathbf{F}}_n$ is the non-linear flux in the direction normal to the cells interface as given below:

$$\underset{{\Omega }_i}{\overset{Allaire}{\int }}{\alpha }_1\nabla ·\mathbf{v}d\Omega \approx \underset{{\Omega }_i}{\int }{\alpha }_{1}d\Omega ·\underset{\partial {\Omega }_i}{\int }{\left(v_n\right)}^{Riem.}dS.$$

(27)

$$\underset{{\Omega }_i}{\overset{Kapila}{\int }}({\alpha }_1+K)\nabla ·\mathbf{v}d\Omega \approx \underset{{\Omega }_i}{\int }({\alpha }_1+K)d\Omega ·\underset{\partial {\Omega }_i}{\int }{\left(v_n\right)}^{Riem.}dS.$$

(28)

where $(n_x,n_y,n_z)$ is a unit normal vector across the boundary $\delta \Omega$, $v_n$ is the velocity normal to the bounded surface area, defined by $v_n=n_{x}u+n_{y}v+n_{z}w$. The source term $\mathbf{G}$ is with regards to the term ${\alpha }_1\nabla ·\mathbf{v}$ of Eq. (18). Following the approach of Johnsen and Colonius [28] the source term is numerically approximated as a surface integral, rather than a volume one while using the same velocity estimate as the one used for the evaluation of the fluxes. The numerical flux $F_n\left(Q\right)$ can be computed considering the states to the left and the right of the high-order approximations of the solutions for the cell interface i and j respectively and expressed as:

$$F_{{\mathbf{n}}_{\mathbf{ij}}}=F_{{\mathbf{n}}_{\mathbf{ij}}}\left({\mathbf{Q}}_{ij,L}^n({\mathbf{x}}_{ij,\alpha },t),{\mathbf{Q}}_{ij,R}^n({\mathbf{x}}_{ij,\alpha },t)\right)$$

(29)

Integrating Eq. (26) over the mesh element i using a high-order explicit finite-volume formulation, the following equation is obtained that incorporates the source term as previously defined:

$$\begin{array}{c}\hfill \frac{d{\mathbf{Q}}_i}{dt}=\frac{1}{\left|V_i\right|}\sum _{j=1}^{N_f}\sum _{\alpha =1}^{N_{qp}}\left(F_{{\mathbf{n}}_{\mathrm{ij}}}\left({\mathbf{Q}}_{ij,L}^n\left({\mathbf{x}}_{ij,\alpha },t\right),{\mathbf{Q}}_{ij,R}^n\left({\mathbf{x}}_{ij,\alpha },t\right)\right)-{\mathbf{a}}_{i,1}^n·v_n^{Riem}\left({\mathbf{x}}_{ij,\alpha },t\right)\right){\omega }_{\alpha }\left|S_{ij}\right|,\end{array}$$

(30)

where ${\mathbf{Q}}_i$ is the volume averaged vector of variables

$${\mathbf{Q}}_i=\frac{1}{\left|V_i\right|}{\int }_{V_i}\mathbf{Q}(x,y,z)dV,$$

(31)

and $F_{{\mathbf{n}}_{\mathbf{ij}}}$ is a numerical flux function in the direction normal to the cell interface between a considered cell i and one of its neighbouring cells j. $N_f$ is the number of faces per element, $N_{qp}$ is the number of quadrature points used for approximating the surface integrals, $|S_{ij}|$ is the surface area of the corresponding face, and ${\mathbf{Q}}_{ij,L}^n({\mathbf{x}}_{ij,\alpha },t)$ and ${\mathbf{Q}}_{ij,R}^n({\mathbf{x}}_{ij,\alpha },t)$ are the high-order approximations of the solutions for cell i and cell j respectively. $\alpha$ corresponds to different Gaussian integration points ${\mathbf{x}}_{\alpha }$ and weights ${\omega }_{\alpha }$ over each face. ${\mathbf{a}}_{i,1}^n$ corresponds to the volume-averaged volume fraction of cell i at time level n. A suitable Gauss-Legendre quadrature numerically approximates the volume, surface and line integrals.

The reconstruction process adopted in UCNS3D [22] follows the approaches of Tsoutsanis et al. [23,37,38,39,40,41,42,43], Titarev et al.[38] that have been previously applied to a wide-range of smooth and discontinuous flow compressible flow problems including laminar, transitional, and fully turbulent settings [44,45,46,47,48,49,50,51,52,53,54,55,56,57,58].

3.2. CWENO scheme

Tsoutsanis’s and Dumbser’s [23] CWENO and CWENOZ schemes for mixed-element unstructured meshes schemes were used to test compressible multiphase components [42,59,60]. The paper shows that the CWENO can produce lesser oscillations near the interface at a cheaper computational cost than the WENO schemes on the same grid. This is made possible by the fact that in the presence of discontinuity, CWENO, unlike the traditional WENO schemes, uses the lower-order polynomials from the directional stencils that might contain a smooth solution for its reconstruction, and in the presence of smooth use, the optimal polynomial in the desired order of accuracy. This slight but significant difference gives CWENO scheme advantages over the WENO. The CWENO will be discussed below briefly, and consequently, the same schemes will be used in capturing and resolving the material interface between the two-phase medium of UNDEX. More details of the CWENO schemes implemented by Tsoutsanis on an unstructured grid can be seen in the paper [23,32].

The optimal polynomial are defined as:

$$p_{opt}(\xi ,\eta ,\zeta )=\sum _{s=1}^{s_t}{\lambda }_s{p}_s(\xi ,\eta ,\zeta ),$$

(32)

where s is the stencil index, with $c=1$ being the central, $c=2,3,.$ being the directional, $s_t$ being the total number of stencils, and ${\lambda }_s$ being the linear coefficients for each stencil, whose sum is equal to 1. The $p_1$ polynomial is not computed directly, but computed by subtracting the lower-order polynomials from the optimum polynomial as follows:

$$p_1(\xi ,\eta ,\zeta )=\frac{1}{{\lambda }_1}\left(p_{opt}(\xi ,\eta ,\zeta )-\sum _{s=2}^{s_t}{\lambda }_s{p}_s(\xi ,\eta ,\zeta )\right).$$

(33)

The CWENO computes a reconstruction polynomial below as a non-linear combination of all the polynomials:

$$p{(\xi ,\eta ,\zeta )}^{cweno}=\sum _{s=1}^{s_t}{\omega }_s{p}_s(\xi ,\eta ,\zeta ),$$

(34)

where ${\omega }_s$ correspond to the non-linear weights assigned to each polynomial, and in regions with smooth data ${\omega }_s\approx {\lambda }_s$, hence obtaining the high-order approximation from the central stencil, and in regions of discontinuous solutions the reconstructed solution will be mostly influenced from the lower-order polynomials of the directional stencils using the Type 3 directional stencils. The reader is referred to the work of Tsoutsanis [61] for the stencil selection algorithms. ${\tilde{a}}_k$ are the reconstructed degrees of freedom; and the non-linear weight ${\omega }_s$ is defined as: the nonlinear coefficients ${\omega }_s$ are computed from the linear ones ${\omega }_s$ as:

$${\omega }_s=\frac{{\tilde{\omega }}_s}{{\displaystyle \sum _{s=1}^{s_t}}{\tilde{\omega }}_s}\mathrm{where}{\tilde{\omega }}_s=\frac{{\lambda }_m}{{(ϵ+{SI}_s)}^b}.$$

(35)

where $ϵ$ is taken as 10^-6 similarly to the WENO reconstructions, $S{I}_s$ denotes the smoothness indicator and $b=4$. The readers are reffered to [23] for a detailed description of the smoothness indicators.

3.3. Fluxes & Time advancement

For the numerical fluxes we employ the (Harten-Lax-van Leer-Contact) Riemann solver of Toro [62]. The solution is advanced in time by a third-order TVD Runge-Kutta method:

$$\begin{array}{cc}\hfill & {\mathbf{U}}_i^{(n+1/3)}={\mathbf{U}}_i^n+\Delta t{R}_i\left(\mathbf{U}\right),\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill & {\mathbf{U}}_i^{n+2/3}=\frac{3}{4}{\mathbf{U}}_i^n+\frac{1}{4}{\mathbf{U}}_i^{(n+1/3)}+\frac{1}{4}\Delta t{R}_i\left({\mathbf{U}}^{(n+1/3}\right),\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill & {\mathbf{U}}_i^{n+1}=\frac{1}{3}{\mathbf{U}}_i^n+\frac{2}{3}{\mathbf{U}}_i^{(n+2/3)}+\frac{2}{3}\Delta t{R}_i\left({\mathbf{U}}^{(n+2/3}\right),\hfill \end{array}$$

(38)

and the time step $\Delta t$ is selected according to:

$$\Delta t=\mathrm{CFL}\frac{1}{d}\underset{i}{min}\left(\frac{h_i}{S_i}\right),$$

(39)

where $S_i$ is an estimate of the maximum in absolute value of the propagation speed in the cell i, d is the number of dimensions, $h_i$ is a characteristic length of the element i and CFL refers to the Courant-Friedrichs-Lewy condition, and a $CFL=0.6$ is used for all tests unless otherwise specified. The reconstruction is performed with respect to the primitive variables since they were previously found [59] to be less oscillatory than the conserved variables.

All the volume/surface/line integrals are approximated by Gaussian quadrature rule suitable for the order of polynomial employed. All the schemes developed are implemented in the UCNS3D CFD code [22] which is written in object-oriented Fortran 2008, employing MPI message passing interface (MPI), and the Open Multi-Processing (OpenMP) application programming interface (API), and the reader is referred to [63] for more details on implementation and performance benchmarks.

4. Test Cases

Using the high-order CWENO schemes in conjunction with the 5-equation interface capturing schemes models [20,21], five underwater test problems were selected to study the cavitation effect of underwater shock near different boundaries.

1.

2D Underwater explosion near a planar wall.

2.

2D Underwater explosion near a free surface.

3.

3D Underwater explosions in a rigid cylinder.

4.

2D underwater explosion of two bubbles placed horizontally near a free surface.

5.

3D Underwater of two bubbles placed vertically in an enclosed region.

4.1. 2D underwater explosion near a free surface

To study the dynamics of bubbles at the interface between air and water in an underwater explosion, a two-dimensional model of an underwater explosion near the free surface of water is considered. Previous studies by various authors have used numerical [5,15,64,65,66,67,68,69] and experimental [69,70] methods to investigate this phenomenon. For the experimental setup, a spherical charge of silver azide with known initial energy is detonated in an open cylinder half-filled with water, and the process involves charge detonation, bubble pulsation, and jet formation [69].

For all of our UNDEX numerical studies in this paper, a pressurised bubble will be used as a substitute for explosives in underwater explosions. The pressurised gas bubble with a diameter of 0.24 m is placed at [0, -0.3] m in the two-phase medium, using a similar setup as seen in Shukla et al.’s work [68]. Our computational domain is a square area with $x\times y$ of dimensions ($[-2,2]\times [-1.5,2.5]$) m as shown in Figure 1, with the air-water interface positioned at the centre $[0,0]$.

Transparent boundary conditions were set up at the top, left, and right surfaces, while the bottom is set to reflecting. The simulation is carried out until $t=1.2$ ms as the final time, and the initial physical parameters of each fluid (density, velocities, pressure, specific ratio, and volume fraction) were set accordingly and presented below:

$$\left(\rho ,u,v,p,\gamma ,\alpha \right)=\left\{\begin{array}{cc}\left(1.225,0,0,{10}^5,2.0,0\right),\hfill & \mathrm{Air}\hfill \\ \left(1250,0,0,{10}^9,2.0,1\right),\hfill & \mathrm{Explosive}\hfill \\ \left(1000,0,0,{10}^5,7.15,0\right),\hfill & \mathrm{Water}.\hfill \end{array}\right.$$

(40)

We employ the high-order CWENO4 finite-volume method developed by Tsoutsanis et al. [23,32], on three grid resolutions (coarse, medium and fine) in our study to investigate the interaction between gas bubbles and the free surface of an air-water medium. Numerical Schlieren visualisations were generated to capture density and pressure gradients at various instants, as depicted in Figure 3 and Figure 4, respectively. The high-pressure air bubble undergoes an explosion, resulting in a shock wave that propagates outward towards the interface, reflects into the water upon hitting the air-water interface, and leaves behind cavitation bubbles, as evident from the pressure gradients in Figure 4. Subsequently, the bubble transforms from a spherical shape to an oval shape, as illustrated in the volume fraction contour in Figure 5.

We conducted a study to investigate the impact of grid resolution on the evolution of the interfacial region between the exploding gas bubble and the air-water free surface using the CWENO4 scheme, employing different grid sizes (400 × 400, 800 × 800, and 1600 × 1600). Changes in the CWENO4 scheme’s behaviour and its ability to capture the interfacial region between the gas bubble and the air-water free surface can be seen in Figure 6 as the grid size is refined. These include reductions in the smearing or numerical diffusion between the free surface and the expanding air bubble, resulting in a sharper description of the interface. The results are free from oscillations which also demonstrates the robustness of the numerical framework employed. However interface sharpening or anti-diffusion techniques could prove beneficial especially for up to 3rd-order methods at coarse grid resolutions since they are the ones that can benefit most from this strategy [35,59,68].

We have performed a quantitative investigation on the medium grid, for the obtained density and pressure profiles along the horizontal axis (x = 0) and compared them with published results that utilize sharpening techniques and anti-diffusive techniques in conjunction with high-order schemes [9,71]. The results, as shown in Figure 7 and Figure 8, are consistent with previously published findings[9,71,72].

Furthermore, we applied both the five-equation models proposed by Allaire et al. and Kapila et al. to investigate the influence of the additional term $K\nabla ·u$ in Kapila’s model, thereby differentiating between the two models. The density and pressure plots, depicted in Fig. Figure 9, demonstrate that both the five-equation diffuse-interface multi-component models of Allaire et al. and Kapila et al. yield comparable results, and are in good agreement with the published works of [73] and [71] . This shows that, for cavitation test cases such as the one under consideration in this paper, we can affirm that 5-equation model of Allaire et al. and Kapila et al. both accurately portrays cavitation formation and collapse. However, this accuracy does not extend to spherical bubble collapse in a liquid, as highlighted in the paper [26]. where the 5-equation model of Allaire et al. falls short in correctly representing the solution of spherical bubble collapse. In Figure 10, we observe that the maximum peak pressure on the finest mesh is approximately 8.5 MPa at x = -1.0 m, which is the location where the shock wave reflects upon hitting the expanding cavity from below. This finding is consistent with the observation made in [73] that traditional WENO schemes, although computationally more expensive than CWENO, fail to prevent the smearing of the air-water interfaces. Our results demonstrate that CWENO is capable of resolving material interfaces with reduced numerical smearing while maintaining their non-oscillatory nature, as reported in [59]. This makes CWENO schemes a favourable choice in terms of computational cost for simulating compressible multiphase or multicomponent flows involving cavitation.

Figure 7. UNDEX near a free surface; plot of the density distribution along the horizontal axis x = 0 at two instants, specifically t = 0.4 ms and t = 0.8 ms. The obtained results are compared to those reported in previous studies of Kokh (2010) [71] and Shyue (2006)[74].

Figure 7. UNDEX near a free surface; plot of the density distribution along the horizontal axis x = 0 at two instants, specifically t = 0.4 ms and t = 0.8 ms. The obtained results are compared to those reported in previous studies of Kokh (2010) [71] and Shyue (2006)[74].

Figure 8. UNDEX near a free surface; pressure plot along horizontal axis x = 0 at instants t = (0.4; 0.8)ms. Comparison are made with the results of [71,74] and [72]

Figure 8. UNDEX near a free surface; pressure plot along horizontal axis x = 0 at instants t = (0.4; 0.8)ms. Comparison are made with the results of [71,74] and [72]

Figure 9. UNDEX near a free surface; density and pressure plot along horizontal axis x = 0 at instants t = 0.2ms using Allaire et al. [21] and Kapila et al. [20] models.

Figure 9. UNDEX near a free surface; density and pressure plot along horizontal axis x = 0 at instants t = 0.2ms using Allaire et al. [21] and Kapila et al. [20] models.

Figure 10. UNDEX near a free surface; density and pressure plot along horizontal axis x = 0 at instants t = 1.2 ms using different grid resolutions.

Figure 10. UNDEX near a free surface; density and pressure plot along horizontal axis x = 0 at instants t = 1.2 ms using different grid resolutions.

4.2. 2D cylindrical underwater explosion near a planar rigid wall.

This test case is carried out to examine the impact of shock loading and cavitation loading or reloading, caused by an underwater explosion, on a rigid structure (e.g., a ship’s hull). Numerous studies, both numerical and experimental, have been conducted on this test case by various researchers [34,65,69,75,76]. The experiment involves detonating an explosive submerged in water underneath a plate.

The physical properties (density, velocities, pressure, specific ratio, and volume fraction) of each fluid involved in these test cases are presented below:

$$\left(\rho ,u,v,p,\gamma ,\alpha \right)=\left\{\begin{array}{cc}\left(1270,0,0,8.29\times {10}^8,2,1\right),\hfill & \mathrm{Explosive}\hfill \\ \left(1000,0,0,{10}^5,7.15,0\right),\hfill & \mathrm{Water}.\hfill \end{array}\right.$$

(41)

The rectangular computational domain in two dimensions (as shown in Figure 11) has dimensions of $\left(\right[-6,6]\times [-6,3\left]\right)$ m. Similar to the previous test case, the explosive is modelled as a high-pressure air bubble with a spherical gas bubble of 1 m diameter placed in water at the origin $[0,0]$. A rigid wall is located at $y=3$ m above the bubble, and a numerical sensor is placed at the centre point of the wall (denoted as P) to measure the load on the wall. The boundary conditions for the upper surface are set as reflecting, while the left, right, and bottom surfaces are set as non-reflecting. A coarse grid is used to discretise the domain, and the simulation is run for a final time of $t=100$ ms.

The 4th-order CWENO scheme is used with the 5-equation diffuse interface model for this test problem. To address the issue of non-physical negative pressure values (i.e., p < 0) that can occur in cavitation regions during the explosion, a cut-off model, described in Section 2.1, is introduced to maintain the pressure at 5000 Pa.

Contour plots of pressure obtained using the 4th-order CWENO scheme at different times during the explosion are shown in Figure 12. The plots shows that the primary loading occurs when the pressurised air bubble explodes, generating shock waves that hit the wall above it at t = 1.5 ms and immediately reflect into the water, causing a rapid drop in density. Rarefaction can be observed in Figure 12 at t = 2 ms, as the pressure decreases around the wall, forming a low-pressure region at t = 3 ms. This pressure drop intensifies, resulting in a cavitation region around the wall at t = 4 ms and t = 5.3 ms. The cavitation collapses and generates a reload impact on the wall, known as cavitation collapse reloading (secondary reloading).

The pressure history of the wall along the horizontal axis is shown in Figure 13. The first peak pressure created by the shockwaves is 6200 MPa for the Allaire et al. model and 6400 MPa for the Kapila et al. model, indicating good agreement with the results of other authors [65,75]. As the bubbles continue to expand, the cavitation reloading impact on the wall occurs at a second peak pressure of 720 MPa for Allaire et al.’s model and 760 MPa for Kapila et al.’s model. It should be noted that although the impact of the shock loading is higher than that of the cavitation reloading, the duration of the cavitation reloading is longer, resulting in significant loading on the wall.

4.3. 3D underwater explosion in an enclosed, rigid cylindrical wall.

This test case has been investigated numerically by various authors [9,34,75,78,79,80]. It serves as an ideal problem for studying shock and cavitation responses in enclosed structures such as a cylindrical tank. The initial conditions for this simulation are based on Xie et al.’s work [4] and are presented below:

$$\left(\rho ,u,v,p,\gamma ,\alpha \right)=\left\{\begin{array}{cc}\left(1770,0,0,2\times {10}^9,2,1\right),\hfill & \mathrm{Explosive}\hfill \\ \left(1000,0,0,{10}^5,7.15,0\right),\hfill & \mathrm{Water}.\hfill \end{array}\right.$$

(42)

The computational setup is illustrated in Figure 14, following that of [4,75].

The explosive, modeled as a gas bubble, has a diameter of 0.1 m and is placed within a rectangular cylindrical wall with dimensions of 0.2286 m in height and 0.0889 m in width. The water in the tank is partially aerated with 15% air. The test case can be simulated in two-dimensional axisymmetric form using cylindrical coordinates, or in three-dimensional form. In our study, the test is conducted in a three-dimensional setup.

Reflecting boundary conditions are imposed on all surfaces, and a hybrid hexa-dominant mesh is used to discretise the computational domain with approximately 0.5 million cells as shown in Figure 16. A CWENO4 scheme is utilized for numerical simulations, and the simulation is carried out until $t=120$ µs. The time evolution of the gas bubble’s pressure is shown in the pressure contour plot in Figure 17. It can be observed that as the bubble explodes and expands, shock waves are generated and impact the side walls of the cylinder at t = 15 µs. These waves are immediately reflected back into the water, striking the expanding bubble at t = 30 µs. This results in rapid changes in density and rarefaction. At t = 45 µs, the shape of the bubble transforms into an oval. Over time, a low-pressure region is observed around the bubble and close to the center of the lateral walls at t = 30 µs, leading to cavitation at t = 45 µs, 60 µs, and 70 µs, as evident in Figure 17. The cavitation pocket is most prominent on the wall at t = 60 µs. The violent collapse of cavitation generates reloading pressure on the walls of the cylindrical structure, as depicted in Figure 15. This phenomenon of cavitation-induced reloading can have a significant impact on enclosed structures, with effects comparable to those of the primary shock-wave loading.

The pressure history plot in Figure 15 reveals that the initial peak pressure obtained in our study using the diffuse interface methods proposed by Allaire et al. is 6600 MPa. This value is in close agreement with the experimental result of Wardlaw et al. [83], where P = 6300 MPa. Other studies by Shyue et al. [9] and Ying et al. [82], which utilized the Allaire et al. diffuse interface methods in conjunction with cut-off cavitation models, reported numerical values in the range of 7000-7200 MPa, respectively.

The second impact of cavitation collapse occurs at a time of 90 µs, with a pressure value of p = 2100 MPa, as evident in our results. This value is lower than the numerical results reported by Shyue et al. [9] and Ying et al. [82], but closer to the experimental value than the latter. Notably, the study conducted by Chen [84] using Smoothed Particle Hydrodynamics (SPH) models exhibited a delay in the cavitation collapse and reload, as the Schmidt model [4,77,85] was employed instead of the cut-off model, as utilized in the present study.

4.4. 2D underwater explosion of two bubbles placed horizontally near a free surface

When two bubbles are positioned close to one another, the local pressure field may vary due to the bubble-bubble interactions. The growth, collapse, and oscillation periods of the bubbles—i.e., the time it takes for a bubble to complete one cycle of expansion and collapse can be affected by these changes in pressure. In a two-dimensional setup Figure 18, two bubbles are positioned horizontally in close proximity to each other near a free surface as an attempt to investigate the influence of the separation distance between the bubbles on cavitation. The separation distance, denoted as $L_o$, is varied at four different values:$(0.3,0.5,0.7,0.9)$ m, between two identical bubbles located 0.3 m below the free surface. This configuration closely is similar to the one utilized by Yu et al.[86]. Additionally, a single bubble with similar computational parameters and dimensions is also simulated for comparative analysis. The domain is discretized using a quadrilateral coarse mesh comprising approximately 0.025 million cells, with all surfaces imposing outflow boundary conditions.

Figure 18. Computational domain setup for the two-dimensional underwater explosion of two bubbles placed horizontally near a free surface.

Figure 18. Computational domain setup for the two-dimensional underwater explosion of two bubbles placed horizontally near a free surface.

Figure 19. Underwater explosion of a single and two bubbles configurations placed horizontally near a free surface; plots of pressure history (at the centre location of the left wall ) for the underwater explosion in an enclosed cylindrical wall with CWENO4 scheme, different inter-bubble distances $(0.3,0.5,0.7,0.9)m$ are considered. [9,81,87]

Figure 19. Underwater explosion of a single and two bubbles configurations placed horizontally near a free surface; plots of pressure history (at the centre location of the left wall ) for the underwater explosion in an enclosed cylindrical wall with CWENO4 scheme, different inter-bubble distances $(0.3,0.5,0.7,0.9)m$ are considered. [9,81,87]

The physical characteristics of the fluid involved in these test cases are described below, which includes parameters such as density, velocity, pressure, specific ratio (or specific heat capacity ratio), and volume fraction:

$$\left(\rho ,u,v,p,\gamma ,\alpha \right)=\left\{\begin{array}{cc}\left(1270,0,0,8.29\times {10}^8,2,1\right),\hfill & \mathrm{Explosive}\hfill \\ \left(1000,0,0,{10}^5,7.15,0\right),\hfill & \mathrm{Water}.\hfill \end{array}\right.$$

(43)

The volume fraction contours presented in Figure 20 effectively illustrate the temporal evolution of the bubbles with the distances $L_o$ ranging from 0.3 m to 0.9 m. As the separation distance $L_o$, increases from 0.3 m to 0.9 m, the height of the free surface center consistently decreases at each time instant. In the same figure, it is observed that for smallest inter-bubble distances $L_o=0.3$ m, the bubbles expand and eventually coalesce or merge.

As observed in previous UNDEX (Underwater Explosion) test cases in this paper, the two bubbles undergo rapid expansion, generating a shock wave that propagates through the water and interacts with the surrounding bubbles. The pressure-time history at the center of the explosion, depicted in the pressure contour plot Figure 19, reveals that the peak pressure caused by the initial shock wave hitting the free surface for a single bubble is 170 MPa, while for the two-bubble case with the largest inter-bubble distance of $L_o=0.9$ m, it reaches 290 MPa, which is 1.7 times higher than that of the single bubble. Furthermore, as evident in Figure 19, the initial peak pressure which is as a result of the intensity of the shock wave on the free surface increases with decreasing inter-bubble distance, from the maximum value at $L_o=0.9$ m to the smallest value at $L_o=0.3$ m. The close proximity of the two bubbles allows for stronger interaction, leading to a higher shock-intensity on the free surface.

The occurrence of cavitation beneath the free surface during a two-bubble explosion is influenced by the size and spacing of the bubbles. In same plot Figure 19, as the inter-distance bubbles increase, the cavitation bubble collapse that results in the cavitation reloading on the free surface is further delayed as seen in the second pressure peak. The second peak pressure associated with cavitation reloading is also more pronounced for the two charged bubbles compared to a single bubble. This observation further supports the conception that two charged bubbles have a higher potential for intensifying structural damage due to cavitation reloading compared to a single bubble. The ratio of the initial peak pressure on the free surface caused by the shock wave to that of the cavitation reloading is about four times for the smaller inter-bubble $L_o=(0.3,0.5)$ m, and twice for the larger inter-bubble distances $L_o=(0.7,0.9)$ m. This shows that for larger inter-bubble distances, the effect of structural damage caused by cavitation is more comparable to that of the primary shock-wave loading than for smaller inter-bubble distances.

In Figure 21, the pressure gradient contour is depicted for that of a single charge and the two bubbles with varying inter-bubble distances. The contour illustrates a region of higher pressure near bubble(s), while the pressure diminishes further away as the bubbles expand. The magnitude of the pressure gradients varies based on the distance between the bubbles.

4.5. 3D Underwater of two bubbles placed vertically in an enclosed region.

To investigate the interaction between two vertically aligned bubbles in water, enclosed in a region, it is best to consider this test case in either an axis-symmetric or three-dimensional setup, taking into account gravitational effects. The bubbles used in this test case are of smaller and identical diameter (d = 5.3 mm). They are initially positioned one above the other, with the first bubble placed at the bottom of the enclosed region and the second bubble above it.

Our numerical setup (Figure 22) uses the similar computational dimensions as the Test case Section 4.3, but with smaller bubble radii. The distance between the first and second bubbles is denoted by $H_b$, and chosen to be 15, 30, and 60, all of which are given in millimetres. Reflecting boundary conditions are applied to all surfaces, and a quadrilateral mesh discretizes the domain with approximately 0.4 million cells. The initial physical properties of both fluid is given below:

$$\left(\rho ,u,v,p,\gamma ,\alpha \right)=\left\{\begin{array}{cc}\left(1.225,0,0,1.57\times {10}^7,2.0\right),\hfill & \mathrm{Air}\mathrm{bubble}\hfill \\ \left(998.2,0,0,{10}^5,7.15,0\right),\hfill & \mathrm{Water}.\hfill \end{array}\right.$$

(44)

The volume fraction contour in Figure 24 displays the shortest inter bubble distance $H_b$ = 15 mm, where the bubbles are in close proximity. As the explosion evolves with time, the water film between the interface of the two interacting gaseous bubbles thins out. This causes a surface tension force to pull the bubbles together until the thin film ruptures and the two bubbles merge together. The two bubbles combine to form a single, larger bubble after coalescing.

Figure 23. Underwater explosion of two bubbles vertically placed in an enclosed region; plots of pressure history (at the centre location of the left wall ) for the underwater explosion in an enclosed cylindrical wall with CWENO4 schemes for inter-bubble lengths of (15, 30, and 60) mm

Figure 23. Underwater explosion of two bubbles vertically placed in an enclosed region; plots of pressure history (at the centre location of the left wall ) for the underwater explosion in an enclosed cylindrical wall with CWENO4 schemes for inter-bubble lengths of (15, 30, and 60) mm

Figure 24. Underwater explosion of two bubbles vertically placed in an enclosed region; volume fraction contours for the 3D underwater explosion in an enclosed for $H_b=30$ mm, rigid cylindrical wall.Captured using CWENO4. The green-yellow areas shows the interface between water and explosive. The two bubbles eventually merge into a single larger bubble.

Figure 24. Underwater explosion of two bubbles vertically placed in an enclosed region; volume fraction contours for the 3D underwater explosion in an enclosed for $H_b=30$ mm, rigid cylindrical wall.Captured using CWENO4. The green-yellow areas shows the interface between water and explosive. The two bubbles eventually merge into a single larger bubble.

The pressure plot history at the centre of the left wall is displayed in Figure 23 for the three inter-bubble instances ($H_b$ = 15, 30, and 60 mm). The pressure gradient contours in Figure 25 shows the cavitation regions are more prominent for the shortest inter bubble distance $H_b=15$ mm for varying time frames t = (1.54, 6, 30, 60 and 90) µs. The cavitation regions in same Figure 25 can be observed to be majorly concentrated at the walls of the cylinder and the surrounding areas close to the bubbles.

5. Conclusions

This research extends the CWENO high-order finite-volume numerical framework employed for simulating multicomponent and multiphase flows, as detailed in [59], to specifically investigate cavitation phenomena in underwater explosions (UNDEX). The application of the CWENO schemes proved successful across a range of intricate two- and three-dimensional test cases related to underwater explosions, yielding noteworthy observations.

Notably, the CWENO schemes showcase a remarkable capability to provide high-resolution and sharply defined representations of shock-cavitation structures, as well as the infinitesimal interfacial region. This is achieved by resolving material interfaces with minimized numerical smearing, eliminating the necessity for additional interface sharpening [68] or compression [71]. This takes away the increase in computational cost associated with incorporating interface sharpening or compression into the diffuse interface models.

The compressibility effect introduced by the $K\nabla ·$ u term in the Kapila et al. five-equation models [20], which distinguishes it from the Allaire five-equation models [21] from our results shows that this additional term does not significantly impact cavitation. The similarity in results between the two models underscores their agreement with other published papers, affirming their reliability in simulating cavitation phenomena.

The results also show that the proximity of the bubbles played a crucial role in intensifying the shock wave impacts and the occurrence of cavitation around free surfaces and rigid boundaries. All these findings contribute valuable insights into dynamics of cavitation within the context of underwater explosions.