1   Introduction

Multiphase flows are prevalent in numerous applications in the scientific and engineering fields. Theoretical and experimental methods are important tools for investigating multiphase flow problems and have achieved considerable success. Various theoretical and experimental research results have been obtained for multiphase flow problems. Petalas and Aziz (2000) presented a new mechanical model applicable to all pipe geometries and fluid properties and proposed new empirical correlations for the friction coefficients at the liquid/wall and liquid/gas interfaces in stratified flows. Nazeer et al. (2021) conducted a theoretical study of the flow of a rheological fluid containing two types of nanoparticles in a steep flow channel and compared the flow dynamics of Casson multiphase flow and Newtonian fluid suspension flow. Marrone et al. (2011b) investigated the wave patterns generated by a high-speed elongated ship with a sharp stem and verified the reliability of the study by comparing its results with the results of a series of experiments on the high-speed bow-breaking problem based on the Athena model ship conducted at the Institute of Marine Engineering-National Research Council Rome (CNR-INSEAN) Towed Pool Laboratory. Cui et al. (2016) conducted a series of small-charge underwater explosion experiments to investigate the behavior of bubbles under gravity and various boundary conditions. Zhang et al. (2015) generated bubbles with an enhanced buoyancy effect by discharging them in a low-pressure tank to analyze the dynamics of large bubbles affected by different buoyancy effects. Despite the availability of various theoretical and experimental methods, the continuous advancements in computational science and technology have led to the proliferation of computational fluid dynamics (CFD) solvers. Mesh-based CFD solvers have been extensively utilized to tackle a variety of multiphase flow problems. However, Eulerian-mesh-based methods track and capture these complex free interfaces only by employing special numerical techniques when dealing with free surfaces or multifluid interfaces. The level-set (Hu et al., 2009) or volume-of-fluid (Li et al., 2019a) methods to capture interfaces in Eulerian-mesh-based methods have reached a high level of maturity and have been extensively applied to various scientific and engineering problems (Gibou et al., 2018). However, their precision and conservation are constrained by the resolution of the mesh utilized. Lagrangian-mesh-based methods can capture interfaces by directly tracking the deformation of the mesh (Li et al., 2019b). However, when the interface deformation is large, such as the undulation or rupture of the free surface or the formation of a jet of cavitated bubbles, this technique becomes less feasible. Hence, CFD researchers have sought alternative techniques for simulating multiphase flow. Because of the inherent limitations of grid-based methods in simulating such problems and the rapid development of grid-free methods in recent years, an increasing number of researchers are utilizing grid-free methods, which are not subject to grid constraints, to solve a variety of fluid dynamics problems (Gotoh and Khayyer, 2018; Khayyer and Gotoh, 2013; Liu and Liu, 2010; Shadloo et al., 2016; Sun et al., 2018; Zhang et al., 2017). Among them, the meshfree particle method is the most popular because it is not limited by the amount of boundary deformation when simulating large deformation problems owing to its natural Lagrangian characteristics and the advantages of gradually complete particle approximation theory. The smooth particle hydrodynamics (SPH) method is a representative Lagrangian particle technique, which demonstrates distinct advantages in addressing issues related to free surfaces or multiphase interfaces (Hu and Adams, 2006; Luo et al., 2016; Monaghan and Rafiee, 2013; Sun et al., 2019b; Zhang et al., 2015b).

The Lagrangian meshfree particle method known as SPH was initially introduced by Lucy (1977) and Gingold and Monaghan (1977) in 1977 to simulate 3D astrophysical phenomena. The SPH method has been developed and applied to address the fluid dynamics problems (Monaghan, 1994; Adami et al., 2012; Morris et al., 1997) and large deformation problems in solid mechanics (Bui et al., 2008; Nonoyama et al., 2015). In recent years, SPH has evolved to encompass numerous variants for simulating various hydrodynamic problems, including free-surface flows (Zhang and Liu, 2018; Kazemi et al., 2020; Antuono et al., 2010; Xie et al., 2021), multiphase flows (Colagrossi and Landrini, 2003; Hammani et al., 2020; Gong et al., 2016; Sun et al., 2021a), and fluid–structure interactions (Zhang et al., 2017; Khayyer et al., 2021b; Yilmaz et al., 2021; Sun et al., 2018; Khayyer et al., 2021a). Although SPH has achieved significant success in resolving incompressible (Hu and Adams, 2007; Lind et al., 2012) or weakly compressible (Colagrossi and Landrini, 2003; Monaghan, 2012) multiphase flows, only a few SPH models are available for simulating strongly compressible multiphase flows. The study of strongly compressible multiphase flows has always played a crucial role in hydrodynamics, as all fluids exhibit some level of compressibility in practical engineering applications. Therefore, the simulation of strongly compressible multiphase flows is practical and relevant. In particular, underwater explosions and cavitation bubbles have attracted increasing attention from CFD scientists. The pressure wave caused by the collapse of cavitation bubbles has been widely utilized in marine exploration (Zhang et al., 2019; Li et al., 2020b). In ocean and hydraulic engineering applications, the impact pressure from high-speed water jets and the radiation pressure waves resulting from bubble collapse can cause significant damage or erosion to marine structures (Lauer et al., 2012). Given that underwater explosions or cavitation bubbles often occur under axisymmetric conditions, the axisymmetric SPH model has been developed to enable rapid numerical simulation and accurate resolution of these issues.

Because of the substantial computational expense associated with the large number of particles required for a full 3D simulation, most SPH simulations are limited to the 2D Cartesian coordinate system. Nonetheless, SPH has made significant advancements in computational acceleration, such as through General Processing Unit (GPU) (Crespo et al., 2015; Mokos et al., 2015; Chen and Wan, 2019) and Message Passing Interface (MPI) (Ferrari et al., 2009; Oger et al., 2016). However, achieving 3D simulations of practical engineering problems still poses significant challenges. Considering this, to avoid the need for full 3D simulation in certain problems exhibiting axisymmetric characteristics, SPH simulations are conducted in a cylindrical coordinate system to reduce computational expenses. Axisymmetric SPH models have been widely used in the field of modeling and solving strongly compressible multiphase flow problems. Omang et al. (2006) constructed new kernel functions based on the fundamental interpolation theory of SPH in both cylindrical and spherical coordinate systems and verified the effective capture of excitations by this method through excitation tube experiments. Li et al. (2020a) proposed a simple, intuitive, and physically meaningful method to derive the SPH formula in the column coordinate system and conducted a simulation study of the classic surface tension test and the bubble rise problem under different conditions. Huang et al. (2022) utilized the axisymmetric SPH model to solve various inlet problems with different axisymmetric structures (including spheres and cones with different inclination angles). Wang et al. (2022) used the axisymmetric Riemann SPH method to simulate high-pressure bubbles near different boundaries. Ming et al. (2014) used the axisymmetric SPH method to analyze the effects of charge parameters on underwater contact explosions and proposed improvements and developments based on the axisymmetric SPH methods. The precise numerical resolution of the continuity equation is important for the SPH method. In recent years, many researchers have conducted in-depth studies. Sun et al. (2019a) proposed the δ-plus-SPH model by employing particle shifting technology (PST) and a diffusion term in the continuity equation. Shi and Huang (2022) introduced a modified numerical diffusion term and a special shift treatment near the multiphase interface to the original δ-plusSPH model. However, their work did not provide a satisfactory solution for the two incompressibility conditions corresponding to the invariant density and divergence-free velocity conditions. Recently, Khayyer et al. (2023) proposed two new schemes, namely, the velocity-divergence error mitigation scheme and the volume-conserving shift scheme. The use of these two schemes improved the resolution of the continuity equation and satisfied the two incompressibility conditions more effectively. Khayyer et al. (2023) further proposed improvements and developments based on the δ-plus-SPH model. However, an urgent problem in simulating strongly compressible multiphase flows is the significant difference in volume change before and after the motion of strongly compressible fluids, which presents an urgent challenge to be addressed in the simulation of strongly compressible fluids.

Recently, the axisymmetric SPH model has been gradually developed for strongly compressible multiphase flow. Nevertheless, only a few studies have addressed the substantial volume change issue in strongly compressible multiphase flows, such as those encountered in underwater explosion scenarios. In strongly compressible multiphase flow problems, the physical process involves one fluid compressing another fluid. During this process, the volume of some fluid particles in the expansion area becomes large, whereas the volume of some fluid particles in the compression area becomes small. When using the SPH method for simulations, the volume of fluid particles will change accordingly, leading to significant changes in the distance between particles. Therefore, to maintain the number of neighboring fluid particles within a certain range, some CFD researchers use variable smoothing length technology (Benz, 1990) to expand or compress the particle support domain based on the density of the particles themselves. Nevertheless, in certain extreme scenarios, particles may undergo excessive expansion or compression, and simply updating the smoothing length may not be sufficient to maintain a sufficient number of neighboring particles. To address this challenge, Sun et al. (2021b) introduced a volume adaptive scheme (VAS) to adjust the particle volumes and maintain a uniform volume distribution within the context of strongly compressible multiphase flow scenarios. However, uneven particle distributions near the interface of multiphase flow are detected by the VAS. In the present study, an improved volume adaptive scheme (IVAS) is developed to further enhance the uniformity of particle distribution.

Recently, the multiresolution technique has been rapidly developed for both grid-based and particle-based methods. To enhance computational efficiency and ensure computational accuracy, the adaptive mesh refinement technique (Freret et al., 2022; Berger and Oliger, 1984) for Finite Element Method (FEM) and the adaptive particle refinement (APR) technique for SPH have been developed (Yang et al., 2023a; Liang et al., 2023). The SPH method utilizes APR through different techniques, such as particle splitting and coalescing (Alimi et al., 2003; Marsh et al., 2011). Kitsionas and Whitworth (2002) implemented the particle splitting technique in astrophysics; subsequently, they made some developments and improvements to the particle splitting technique (Kitsionas and Whitworth, 2002; 2007). Feldman and Bonet (2007) significantly advanced the development of APR technology, which involves the splitting of parent particles into multiple subparticles. In addition, the three conservation laws are satisfied by the mass, volume, density, velocity, and pressure of the subparticles. This approach was further enhanced by Reyes López et al. (2013), Omidvar et al. (2012), and Vacondio et al. (2012). In the study conducted by Vacondio et al. (2012), one coarse particle is divided into seven fine particles, and two fine particles are combined into one coarse particle. However, in the process of splitting and merging, the thinning and coarsening speeds do not match. The particle coarsening technology created by Vacondio et al. (2012) was enhanced by Barcarolo et al. (2014). In their method, a coarse particle is not eliminated when it splits into four fine particles but is excluded from the SPH calculation. The fine particles are eliminated, and the coarse particles reenter the SPH calculation after exiting the refining zone. However, pressure or velocity oscillations could occur close to the edge of the refining zone. APR was developed by Lyu et al. (2022) for use in 3D water entry simulations. In the present study, the IVAS technique is combined with the APR technique to further improve the computational accuracy for modeling the strongly compressible multiphase flow problem.

The excessive variation of fluid particle volume, the pressure oscillation at the interface of multiphase flow, and the capture of pressure shock waves are key points and difficulties in simulating strongly compressible multiphase flow using the SPH method. This work presents an improved adaptive multiresolution SPH model for strongly compressible multiphase flow. The IVAS technique is developed based on VAS to regulate the particle volume and obtain more uniform particle distributions near the interface of multiphase flows. On this basis, the IVAS technique combined with the APR technique (IVAS-APR) is developed to further improve computational accuracy and efficiency. In the present IVAS-APR technique, first, the IVAS is used to regulate the particle volume, and the variation in particle volume is controlled within a narrow range. Then, the APR is applied, a high resolution is used for the local region that covers the front and back of the pressure shock wave and near the interface of the multiphase flow, and the other computational domains are modeled using coarse particles.

The remainder of this paper is organized as follows: In Section 2, the axisymmetric SPH method and numerical treatments are introduced. In Section 3, the adaptive multiresolution scheme (AMRS) is described in detail. In Section 4, four numerical examples are investigated to verify the effectiveness of the proposed method. Finally, in Section 5, the conclusions are drawn.

2   Strongly compressible multiphase axisymmetric SPH model

2.1   Governing equations

Lagrangian meshfree techniques include the SPH method. SPH uses a large number of particles to describe the spatial distribution of a physical field. The weights of nearby particles are used to approximate the physical properties of each particle. The governing equations of the SPH method are derived based on the basic principles of fluid mechanics and typically include mass, momentum, and energy conservation equations. The same applies to the basic principles of the governing equations of the axisymmetric SPH model, but they need to be expressed in the polar coordinate system (Fang et al., 2022) as follows:

where the terms pressure, density, velocity, gravitational acceleration, and internal energy per unit mass are represented by the symbols $p, \rho, u, g$, and $e$, respectively. The coordinates in the axial, circumferential, and radial directions are represented by the symbols $r, \theta$, and z, respectively. Following the simplified treatment of $\partial f / \partial \theta=0(f=u, p, e)$ and $u^{\theta}=0$, the fluid problem exhibits rotational symmetry when it is axisymmetric, indicating that the properties of the fluid are uniform along the circumferential direction and do not change with the angle. The governing equations for the axisymmetric issue are reduced as follows:

2.2   Discretization equations

By applying SPH approximations (Liu and Zhang, 2019), the discretization equations can be expressed as follows:

where the terms pressure, density, velocity, gravitational acceleration, and internal energy per unit mass are represented by the symbols $p, \rho, u, g$, and $e$, respectively. The coordinates in the axial, circumferential, and radial directions are represented by the symbols $r, \theta$, and $z$, respectively. The subscripts $i$ and $j$ are the coordinates of particles, and W is the kernel function with a smoothing length h. In this study, the B-spline kernel function is adopted (Monaghan and Lattanzio, 1985).

Based on Brookshaw's model (Brookshaw, 2002), the axisymmetric adaptive multiresolution SPH method is developed in this study. In 2D fluid problems, the initial volume of an SPH particle is $V_{0}=(\Delta x)^{2}$, where $\Delta x$ is the initial distance between particles. Then, according to the density of the fluid itself, the particle mass is calculated as $\rho(\Delta x)^{2}$. However, in axisymmetric SPH, the 3D fluid problem with axisymmetric characteristics is simplified into a 2D fluid problem. Each SPH particle should still represent a 3D fluid torus with a mass of $2 \pi r \rho(\Delta x)^{2}$. Consequently, the torus cross-section of the fluid is considered the particle volume, which is calculated as follows, m stands for mass:

2.2.1   Conservative corrected SPH differential operators

The integration of the function is approximated by particles, resulting in the following discrete particle approximation equation, $(f=u, p, e)$:

where V is the volume of the particle. The term $W_{i j}$ is the condensed form of $W\left(\left|\boldsymbol{r}_{i}-\boldsymbol{r}_{j}\right|; h_{i j}\right)$, and $h_{i j}=\left(h_{i}+h_{j}\right) / 2$ is utilized to ensure particle interaction symmetry and, thus, the conservation property of the numerical scheme. The initial smoothing length to particle spacing $\Delta x_{i}$ ratio in the 2D simulation is $h_{i} / \Delta x_{i}=2.0$. This work applies the SPH approach to the strongly compressible multiphase flow problem with abrupt changes in particle spacing. To increase accuracy and efficiency, the smoothing length should be adjusted over time.

To enhance the accuracy of the kernel gradient, kernel gradient correction is utilized in this study (Randles and Libersky, 1996). Then, the approximations of $\nabla f$ and $\nabla \cdot \boldsymbol{f}$($f=u, p, e$) can be expressed as follows:

where the velocity divergence is discretized using the corrected divergence operator. However, the antisymmetry of the particle pair i and j disappears, making it impossible to ensure proper linear momentum conservation when using the gradient operator to discretize the pressure gradient.

In this study, the method proposed by Sun et al. (2021b) is employed to discretize the pressure gradient. The antisymmetry $\nabla_{i} W_{i j}=-\nabla_{j} W_{i j}$ arises from the use of $h_{i j}=\left(h_{i}+h_{j}\right) / 2$. The conventional gradient operator can be rewritten as follows:

Alternatively, another type of variable gradient can be created by substituting the modified kernel gradient $\nabla W_{i j}^{C}$ for $\nabla W_{i j}$, i.e.,

The aforementioned equation can be rewritten as follows (Sun et al., 2021b):

which ensures the precise conservation of total linear momentum.

2.2.2   Discretization equations with diffusive terms

In this work, the pressure and velocity fields for all conservation laws (i.e., mass, momentum, and energy) are stabilized using the full diffusion term formula created by Antuono et al. (2010). The governing equations of the axisymmetric SPH model are expressed as follows:

where $\rho_{i j}=\left(\rho_{i}+\rho_{j}\right) / 2, \boldsymbol{u}_{i j}=\boldsymbol{u}_{j}-\boldsymbol{u}_{i}, \boldsymbol{r}_{i j}=\boldsymbol{r}_{j}-\boldsymbol{r}_{i}$, and $c_{i j}= \left(c_{i}+c_{j}\right) / 2$. In addition, when particles $i$ and $j$ have the same phase, $\varphi_{i j}=1$; otherwise, $\varphi_{i j}=0$.

Using the artificial density diffusion method proposed by Marrone et al. (2011a) to control the cumulative error of the evolution of mass density, $D_{i j}$ can be expressed as follows:

where the fluid particles only interact with other nearby particles in the same phase, and the density gradient $\langle\nabla \rho\rangle^{L}$ is computed using the adjusted spatial gradient in Eq. (6). Similarly, $\varepsilon_{i j}$ is expressed as follows:

The two parameters $\delta$ and $\kappa$ do not need to be modified to account for varied problems because they are only capable of changing within a certain range, specifically δ = 0.1 and κ = 0.1. However, the values of the two parameters α and β need to be set in artificial viscosity, as their values directly impact the quality of the simulation results. The results obtained with α = 1 and β = 2, which are relatively large values, are the most favorable for simulating explosion pressure waves and forceful impacts (Monaghan, 2005). When the artificial viscosity coefficient has a large value, the calculation and simulation processes tend to be stable because, for high-speed, high-pressure, and instantaneous problems, such as explosions or impacts, the simulation itself is less stable and often prone to calculation collapse and numerical divergence. However, for weakly compressible or incompressible fluids, smaller artificial viscosity factors, such as α = 0.02 and β = 0, are typically used (Marrone et al., 2011). When simulating such problems, the artificial viscous force cannot contribute a large proportion to the control equation, which would result in significant errors in the calculation results and significant numerical dissipation.

To prevent the penetration of particles near the multiphase flow interface and maintain the crispness of the interface, an interface sharpness force $F_{i}^{\text {sharpness }}$ is typically introduced to the momentum equation. Where particle j are neighboring particles of a phase different from particle $i$, $F_{i}^{\text {sharpness }}$ can be written as follows:

In the current work, ε is set as 0.08.

2.3   Equation of state

The equation of state, denoted as $p=f_{\text {state }}(\rho, e)$, which describes the relationship between the pressure, density, and energy of the fluid, must be applied to close the governing equations. The ideal gas equation of state (Toro, 2013) for the fluid simulated in Subsections 4.1, 4.2, and 4.3 is expressed as follows:

where the value of the adiabatic index, denoted by γ, is 1.4. The sound speed is calculated using the formula $c= \sqrt{\gamma p / \rho}$, as per Toro (Toro, 2013).

The energy and physical compressibility of water must be considered when solving strongly compressible multiphase flow problems, such as underwater explosions discussed in Subsection 4.4. Therefore, water is subjected to the Mie–Gruneisen equation of state (Steinberg, 1987). The pressure of water during compression is determined using the following formula:

The pressure of water when it expands is determined using the following formula:

In the aforementioned equations, $a$ is the volume correction coefficient, $c_{0}$ is the reference sound, $\rho_{0}$ is the initial density of water, $\gamma_{0}$ is the Gruneisen coefficient, and $S_{1}, S_{2}$, and $S_{3}$ are the linear Hugoniot coefficients. $\eta=\rho / \rho_{0}$ with $\rho_{0}$ being the initial mass density of the water. The variable $\mu=\eta-1=\rho / \rho_{0}-1$, where $\rho_{0}$ is equal to 1 000 kg/m³. The states of compression and expansion of water are denoted by $\mu>0$ and $\mu < 0$, respectively. For the explosive gas generated by the TNT charge following the underwater explosion instance discussed in Subsection 4.4, the widely used Jones–Wilkins–Lee (JWL) equation of state (Dobratz, 1981) is utilized as follows:

Herein, $\eta=\rho / \rho_{0}$ with $\rho_{0}$ being the initial mass density of the explosive gas and $\rho_{0}$ = 1 630 kg/m³. The fitting coefficients $A, B, R_{1}, R_{2}$, and $\omega$ are 3.712 × 10¹¹, 3.21 × 10⁹, 4.15, 0.95, and 0.3, respectively. In addition, e = 4.29 × 10⁶ J/kg.

2.4   Numerical treatments for SPH

To increase computational accuracy and stability, PST and Extended-SPH (X-SPH) technologies are coupled in this work. PST is applied to the interior of the fluid, whereas X-SPH is applied close to the multiphase flow interface. In addition, computational stability is enhanced using adaptive smoothing lengths.

This study employs the PST developed by Xu et al. (2009), which can be expressed as follows:

where $\hat{x}_{i}, \delta \boldsymbol{x}_{i}, u_{\text {max }}$, and $C$ are the average distance between neighboring particles j and i, adjustment distance, current maximum velocity, and constant coefficient, respectively. N is the number of neighboring particles that surround particle i. The unit displacement vector between particles i and j is denoted by n_ij. The proposed adaptive multiresolution SPH model only considers the adjustment of the tangent direction for the particles close to the multiphase flow interface. Numerous studies describe in detail the implementation of tangential shifting. Khayyer et al. (2019) used an optimized PST to maintain the regularity of phase interface and free-surface particles. Lyu and Sun (2022) used transition PST to prevent disordered particle distribution and tensile instability caused by negative pressure. To ensure the smoothness of the multiphase flow interface, the PST of particles near the interface needs to be modified. Only the tangential modified displacement vector is retained, the normal vector is set as 0, and the expression is written as follows:

in which $\boldsymbol{x}_{i\text{new}}$ is the new coordinate, $\mathbf{I}$ is an identity matrix, and $\boldsymbol{n}_{i}$ is the normal vector along the interface of particle i.

The variables of density, velocity, and pressure are modified by the Taylor series, which should be written as follows:

where $\varphi_{i\text{new }}$ is the new speed, density, and pressure.

Monaghan (1989) was the first to propose X-SPH technology as a means of enhancing particle dispersion quality. X-SPH technology was utilized by Balsara (1995) to adjust the velocity and enhance the computational stability. In this research, the particle velocity near the multiphase interface is smoothed using X-SPH, thereby enhancing the stability of the multiphase flow interface calculation, as follows:

where $\varepsilon_{\mathrm{X} \text {-SPH }}$ is the coefficient with a value of 0.001, and $\delta \boldsymbol{u}_{i}$ is the adjusted velocity vector. The updated velocity vector $\boldsymbol{u}_{i\text{new}}$ is written as follows:

The new velocity vector is represented by $\boldsymbol{u}_{i\text{new}}$.

The proposed adaptive multiresolution SPH model utilized adaptive smoothing lengths in this study. The particle density can be utilized to adjust the smoothing length of the particles (Benz, 1990) as follows:

where $h_{0}$ is the initial value of the smoothing length and dim is the dimension of the problem. Benz (1990) derived Eq. (26), and by taking the derivative of both sides, the formula can be transformed into the following expression:

2.5   Boundary condition

In axisymmetric SPH algorithms, dealing with symmetry axes has been consistently challenging because of the absence of neighboring particles in the r < 0 region (García-Senz et al., 2009). Moreover, the simulation exhibits a nonphysical penetration phenomenon because the particles near the axial symmetry axis have smaller support domains than the internal particles. However, when r = 0, the 1/r term might appear singular, which would cause the simulation to end. A dynamic mirror border is applied to the symmetry axis to solve this issue (Joshi et al., 2021). Every sub-time step generates a similar mirror particle on the opposite side of the symmetry axis for every particle in the 0 < r < kh region, the k generation represents a constant, usually taking the value 2 or 3. The mirror particle conveys information, which can be written as follows:

where the SPH fluid real particle is denoted by the subscript s and the associated dynamic mirror particle that is close to the symmetry axis is denoted by the subscript w.

2.6   Time stepping approach

The Navier – Stokes equation is solved by employing the leapfrog integral approach. To ensure computational stability, the time step must be constrained to meet specific requirements.

Given that the changes in internal energy within SPH simulations necessitate the careful selection of the time step size, the magnitude of time step $\Delta t$ is obtained using the following equation:

Monaghan et al. (1992) proposed a method to estimate the time step $\Delta t_{\mathrm{cfl}}$ in accordance with the CFL criteria. This formula can be expressed as follows:

The time step estimation formula with viscosity term was introduced by Morris et al. (1997), which can be written as follows:

The final time step $\Delta t$ should be assigned the minimum value, which is derived as follows:

3   AMRS

The axisymmetric adaptive multiresolution SPH model of strongly compressible multiphase flow used in this study is inspired by the research of Sun et al. (2021a). The VAS proposed by Sun et al. (2021a) suffers from particle aggregation near the excitation front when simulating shock waves. In addition, they only consider the problem of excessive differences in particle volume changes in strongly compressible multiphase flows. In the simulation of strongly compressible multiphase flows, the multiphase interface and the pressure shock wave are the core regions of the simulation, which should be computed with higher resolution. An IVAS has been developed to control the particle volume and achieve a more uniform particle distribution near the interface of multiphase flows. In addition, we propose the AMRS combined with IVAS-APR to further enhance the computational accuracy of the pressure shock wave and crucial regions of the interface. In the current AMRS, first, the change in particle volume is controlled within a small range by adjusting the particle volume using IVAS. Then, the APR is employed, where coarse particles mimic other computational domains, and a high resolution is applied to the local region covering the interface of the multiphase flow and the front and rear of the pressure shock wave.

3.1   IVAS

During the simulation of strongly compressible multiphase flow, the particle volume may vary significantly from the initial value, resulting in a notable alteration in the particle distance. By adjusting the particle volume using the VAS developed by Sun et al. (2021a), all fluid particles in the computational domain can maintain a uniform volume distribution. However, particle aggregation near the front of the shock wave due to VAS is a cause for concern. This study proposes IVAS to solve this problem.

3.1.1   VAS

In the VAS, fluid particles are adaptively divided or combined based on the volume values of each particle. In particular, as shown on Figure 1(a), a particle is divided into four daughter particles that are dispersed at the four vertices of a square when its volume exceeds 1.6V₀, where V₀ is the original volume of the particle. The physical details of the newly created daughter particles are expressed as follows:

Figure  1  Sketch of particle splitting and merging (Fang et al., 2022)

Download: Full-Size Img

where the mother particle is denoted by the subscript m and the daughter particle is denoted by the subscript $d$. $F_{\text {state }}^{-1}$ state stands for inverse state function. Two particles combine to form a single mother particle when their respective volumes are less than $(2 / 3) V_{0}$, and their separation is less than $\sqrt{(2 / 3) V_{0}}$, as shown on Figure 1(b). The physical details of the mother particle created during the merging process are expressed as follows:

The use of VAS can help maintain a homogeneous particle distribution by preventing particle aggregation in the compression region and particle sparsity in the expansion region. Most significantly, the use of VAS guarantees the conservation of mass, momentum, and energy, which is a critical aspect of maintaining high numerical accuracy in the simulation of situations involving extremely compressible multiphase flows. However, in the compression region, VAS causes particles to gather near the front of the shock wave.

3.1.2   IVAS

To address the problem of particle aggregation at the front of the shock wave, this study develops IVAS. In IVAS, the PST is utilized to make the particle distributions more uniform, and the calculations of particle splitting and merging are both conducted before the time integration of the SPH calculation of the discrete control equation. When using the PST for multiphase flow, the PST developed by Xu et al. (2009) should be modified. For particles at the interface of multiphase flows, they only move along the tangent direction of the interface, and the normal component of the particle shift vector is set as 0. Moreover, to deal with different particle sizes, the particle movement vector is further modified as follows (Long et al., 2017):

As shown in Figure 2, the particle merging process in VAS is placed after the SPH calculation, which allows more particles to participate in the SPH calculation. However, a large number of particles accumulate at the front edge of the compressed fluid region, resulting in an extremely uneven particle distribution, which makes the calculation highly unstable. The calculation of particle splitting and merging is placed before the SPH calculation, and the PST is utilized to make the particle distribution more uniform, which can prevent calculation instability caused by the aggregation of some particles that are not merged in time and enable the entire calculation domain to reach a unified elementary resolution.

Figure  2  Flowchart of the VAS and IVAS

Download: Full-Size Img

3.2   APR technique

This study employs the APR technique, which was estab lished by Yang et al. (2023a), after the IVAS computation. The specific physical variables define the high-resolution region. A coarse particle is divided into four smaller particles at the four vertices of the square when it enters a highresolution refining zone of space. Refining zones can be moved over time, locked in place, or dynamically expanded and contracted in response to the needs of the computing task. The graphic depicting particle splitting is shown in Figure 3. In contrast to IVAS particle splitting when a coarse particle in APR is refined into four daughter particles, the parent particle remains in the system but is no longer involved in the computation of the Navier–Stokes governing equations or other particles and is referred to as inert particles in this study. Conversely, particles that participate in the SPH calculation are referred to as active particles. As shown in Figure 4, a transition region must be included to ensure high accuracy and stability in the interaction between fine and coarse particles. The fine particles produced by splitting become inert particles and do not participate in the calculation of the Navier–Stokes governing equations, whereas the coarse particles remain active in the transition region and are continuously used in the computation of the physical parameters of other particles. Coarse particles become inert as they pass through the designated high-resolution refining zone, whereas fine particles become active as they pass through the designated high-resolution refining zone.

Figure  3  Schematic of the particle splitting pattern (Yang et al., 2023)

Download: Full-Size Img

Figure  4  Schematic of APR technology (Yang et al., 2023a)

Download: Full-Size Img

This study forms a high-resolution area in the computational domain using APR technology. However, the uniformity of particles in this area still needs to be ensured using PST. Moreover, the coupling of APR and PST is not a simple process. Recently, Yang et al. (2023b) proposed the coupling of APR and PST to solve the water entry and exit problems. In this study, an auxiliary virtual particle technology is applied to realize the coupling of APR and PST. As shown in Figure 5, a series of auxiliary virtual particles are generated with the same resolution as the refined particles near the outside boundary of the refined area. The density, velocity, pressure, and other physical quantities of the auxiliary virtual particles are obtained by interpolation of the fluid particles around the auxiliary virtual particles. The specific interpolation formula is written as follows:

Figure  5  Auxiliary virtual particle technology

Download: Full-Size Img

The auxiliary virtual particles do not participate in the calculation of the Navier–Stokes equation, only in PST, as the auxiliary particles replace the coarse particles that participate in the displacement correction calculation of refined particles.

Figure 6 illustrates the flowchart of the adaptive multiresolution SPH model of the axisymmetric strongly compressible multiphase flow in this study. In the present AMRS, first, IVAS is used to adjust the particle volume, ensuring that the change in particle volume is controlled within a small range. Subsequently, APR technology is applied to improve the calculation accuracy of the core area.

Figure  6  Flowchart of the AMRS

Download: Full-Size Img

4   Example tests

To verify the accuracy and effectiveness of the proposed adaptive multiresolution SPH model for modeling the axisymmetric strongly compressible multiphase flow, four numerical examples of typical strongly compressible multiphase flow problems are tested.

The first calculation example is the classic Sod shock tube problem, which verifies the accuracy and validity of the IVAS and ensures that the adaptive multiresolution SPH model for axisymmetric strongly compressible multiphase flow can achieve a high resolution at the multiphase flow interface. The second calculation example is the blast shock wave problem, which further validates the accuracy of the adaptive multiresolution SPH model for axisymmetric strongly compressible multiphase flow and the IVAS in the problem of strongly compressible multiphase flow with lower compressibility. The third calculation example is the Sedov point explosion problem, which verifies the effectiveness of the adaptive multiresolution SPH model for axisymmetric strongly compressible multiphase flow in capturing pressure shock waves with high precision and validates the accuracy and effectiveness of IVAS in simulating strongly compressible multiphase flow problems under extremely large compression ratios. The fourth calculation example is the underwater explosion impact problem, which further verifies that the adaptive multiresolution SPH model for axisymmetric strongly compressible multiphase flow simultaneously performs high-resolution calculations at the multiphase flow interface and pressure shock waves. The accuracy and effectiveness of the model also demonstrate its suitability for the numerical simulation of practical problems, such as underwater explosion shock wave propagation. The model can provide valuable reference and assistance in ocean engineering problems and serve as a foundation for further research on strongly compressible multiphase flows. Scientific researchers have contributed their own modest efforts to this issue.

4.1   Sod problem

To verify the effectiveness of the IVAS and the present AMRS SPH model for modeling strongly compressible multiphase flows, the Sod problem was investigated. The Sod test (Sod, 1978) is a typical example of a strongly compressible multiphase flow problem. Its basic idea involves a membrane separating an endlessly long tube. Two distinct gases are placed within this tube, i.e., low-density gas on one side and high-density gas on the other side. As illustrated in Figure 7, the computational domain of the 2D Sod (Toro, 2013) is typically characterized by a circular region. Here, F1 stands for Fluid 1, whose initial condition is defined as ($\boldsymbol{u}, \rho, p, e$) $=\left(0, 1 \mathrm{~kg} / \mathrm{m}^{3}, 1 \mathrm{~Pa}, 2.5 \mathrm{~J} / \mathrm{kg}\right)$, and F2 stands for Fluid 2, whose initial condition is ($\boldsymbol{u}, \rho, p$, $e)=\left(0, 0.125 \mathrm{~kg} / \mathrm{m}^{3}, 0.1 \mathrm{~Pa}, 2 \mathrm{~J} / \mathrm{kg}\right)$. Redefining the computational domain is necessary to verify the accuracy of the present axisymmetric AMRS SPH model. Figure 8 illustrates the dispersion of the fluid particles in the 2D axisymmetric plane. The r-axis has only one layer of fluid particles, and these fluid particles are copied up and down at each time step to obtain a total of 16 layers of duplicated particles. $\Delta x$ is the separation between successive layers, and the physical properties of the duplicated particles correspond to those of the fluid particles. Then, dynamic mirror virtual particles are generated at the left and right edges of the computational domain to prevent truncation errors in the kernel function. Each fluid particle interacts with replicated particles, mirror particles, and other fluid particles within the support domain.

Figure  7  Sketch of the computational model

Download: Full-Size Img

Figure  8  Sketch of the corresponding particle distribution in 2D axisymmetric plane

Download: Full-Size Img

Figure 9 illustrates the particle distribution at four different times during the entire simulation process. Notably, to enhance clarity in the particle distribution, the particle distribution on the r-axis is magnified in the region of z ∈ [− 0.3 m, 0.3 m ]. Under high-pressure and high-density conditions, Fluid 1 gradually moves to the left and constantly compresses Fluid 2. Fluid 2 is also forced to move to the left under the influence of Fluid 1. At the same time, the interface between the two fluids constantly moves to the left. In this process, Fluid 1 expands, and its particle volume increases. By contrast, Fluid 2 is continuously compressed, and its particle volume becomes smaller. This situation results in a sparse particle distribution in the expansion region, a dense particle distribution in the compressed region, and an extremely uneven particle distribution throughout the entire problem domain. Using IVAS, particles in the expansion region are divided in time, whereas particles in the compressed region are merged in time. In addition, particles are adjusted to suitable positions by an improved multiphase PST to maintain a uniform particle distribution and further enhance calculation accuracy. The particles near the interface of the multiphase fluid are further refined, and the refining zone changes with the movement of the interface. Throughout the entire simulation process, the multiphase interface remains in the high-resolution calculation region, thereby improving the calculation accuracy of particles near the interface and effectively reducing the nonphysical shocks caused by discontinuous pressure at the interface during the simulation process. At the same time, this approach can effectively reduce the computational costs.

Figure  9  Particle distributions at different times

Download: Full-Size Img

Figure 10 shows the cloud map of velocity, pressure, and density distributions at t = 0.2 s. Notably, the rarefaction wave propagates to the left and the contact wave moves to the right. Figure 11 displays the velocity, pressure, and density distributions on the z = 0 plane at a resolution of $\Delta x$ = 0.0025 m, comparing them with the reference solution (Toro, 2013) and Riemann SPH results (Fang et al., 2022). The results indicate that the numerical results of this study are consistent with the reference solution (Toro, 2013), and the nonphysical oscillations are significantly reduced compared with the Riemann SPH results, especially in the regions of r ∈ (0.64 m, 0.8 m) and r ∈ (0.72 m, 0.86 m), because the particles at the multiphase interface are always in the high-resolution region for the AMRS SPH model.

Figure  10  Density, pressure, and velocity distributions at t = 0.2 s

Download: Full-Size Img

Figure  11  Pressure, density, and velocity distributions on the z = 0 plane with $\Delta x$ = 0.002 5 m

Download: Full-Size Img

In addition, the present SPH results are compared with the 2D Moving-least-squares weighted essentially nonoscillatory SPH (MWSPH) (Avesani et al., 2014) and Riemann SPH (Fang et al., 2022) results at a resolution of $\Delta x$ = 0.01 m, as shown in Figure 12. Notably, the current AMRS SPH, Riemann SPH (Fang et al., 2022), and MWSPH (Avesani et al., 2014) can handle the discontinuity problem effectively. However, pressure and density produce some small oscillations in the present SPH results. The numerical dissipation is small in the region r ∈ (0.6 m, 0.8 m), which is consistent with the reference solution (Toro, 2013).

Figure  12  Pressure and density distributions on the z = 0 plane with $\Delta x$ = 0.01 m

Download: Full-Size Img

Figure 13 presents the history of the number of particles in Fluids 1 and 2 over time at a resolution of $\Delta x$ = 0.002 5 m. As the particles expand, the volume of Fluid 1 becomes larger, and the number of particles increases from 12 932 to 16 226. Conversely, as the particles are compressed, the volume of Fluid 2 becomes smaller, and the number of particles decreases from 25 132 to 21 533. Figure 14 shows the time histories of Fluid 1, Fluid 2, and total fluid particles. The use of IVAS resulted in a 25.47% increase in the number of particles in Fluid 1 and a 14.32% decrease in the number of particles in Fluid 2 at the end of the simulation. However, the total number of particles remains nearly constant throughout the simulation process, indicating that the average resolution of particles, except in the high-resolution region, changes within a limited range under the influence of IVAS. Thus, IVAS improves the numerical accuracy in strongly compressible multiphase flow simulation.

Figure  13  Time histories of the number of particles

Download: Full-Size Img

Figure  14  Time histories of the number of particles

Download: Full-Size Img

This case demonstrates that IVAS effectively handles substantial volume changes of particles, and the current AMRS is effective and robust in forming refined high-resolution regions around the multiphase interface to improve computational accuracy. The current SPH model is well-suited for simulating strongly compressible multiphase flows.

4.2   Blast shock wave problem

To further verify the effectiveness of the present AMRS SPH model for simulating strongly compressible multiphase flow with the low compression process, the blast shock wave test is utilized. Avesani et al. (2014) introduced the blast shock wave problem, which serves as the second benchmark test. Two different fluids of the same density are placed in an airtight space, one of which is in a state of high pressure, and the other is in a state of normal pressure. The schematic diagram of the computational domain of the 2D axisymmetric plane, as well as the related particle distribution, is identical to that of the 2D Sod problem, as shown in Figures 7 and 8. However, the initial conditions of the two fluids in the blast shock wave problem differ from those of the 2D Sod problem, x and y represent the horizontal and vertical coordinates of the particle respectively, as follows:

Figure 15 shows the particle distributions at three different times. Under high-pressure conditions, Fluid 1 gradually moves to the left, whereas Fluid 2 is compressed and also moves to the left. However, in contrast to the Sod problem, the movement of the multiphase interface is slow, and the degree of fluid compression is low. Similarly, the particles near the multiphase interface are refined, and the local region near the interface is set as a high-resolution region.

Figure  15  Particle distributions at different times obtained by the present SPH method

Download: Full-Size Img

Figure 16 illustrates the velocity, pressure, and density distributions derived by the present axisymmetric AMRS SPH model. The velocity and pressure distributions are similar to those in the Sod problem, with the contact wave moving to the right and the rarefaction wave propagating to the left. The difference is that a high-density region forms near the multiphase interface on the left, whereas a low-density region forms on the right.

Figure  16  Density, pressure, and velocity distributions at t = 0.2 s

Download: Full-Size Img

Figure 17 shows the velocity, pressure, and density distributions at a resolution of $\Delta x$ = 0.002 5 m, which are compared with those obtained using the reference (Avesani et al., 2014) and Riemann SPH (Fang et al., 2022) solutions. The results show that the calculated results obtained by the present SPH are consistent with those obtained using the reference (Avesani et al., 2014) and Riemann SPH (Fang et al., 2022) solutions and play a significant role in reducing nonphysical vibration.

Figure  17  Pressure, density, and velocity distribution with $\Delta x$ =0.002 5 m

Download: Full-Size Img

This case demonstrates that the present AMRS SPH model is effective and suitable for modeling strongly compressible multiphase flow with the low compression process.

4.3   Sedov point explosion problem

To further verify the effectiveness of the present AMRS SPH model in addressing high compressibility and accurately capturing pressure shock waves in explosion problems, the Sedov point explosion problem (Sedov, 2018) is employed. The sketch of the Sedov point explosion problem is shown in Figure 18. The Sedov point explosion problem can be explained as follows: Under the initial pressure disturbance, cylindrical explosion waves propagate in a self-similar manner with a large impact force and low density. The extreme discontinuity and approaching vacuum conditions coexist in this problem, making the SPH simulation extremely challenging. The initial particle distribution is shown in Figure 19. With r = 0.025 m, Fluid 1 represents a cylindrical explosion point, and its initial condition is set as $(\boldsymbol{u}, \rho, p, e)=\left(0, 1 \mathrm{~kg} / \mathrm{m}^{3}, \frac{(\gamma-1) e_{0}}{\pi r^{2}} \mathrm{~Pa}, 0\right)$, where $e_{0}=1 \mathrm{~J} / \mathrm{kg}$ is the specific internal energy and γ = 1.4 is the adiabatic constant. The remaining computational domain is occupied by Fluid 2, with the initial condition set as ($\boldsymbol{u}$, $\rho, p, e)=\left(0, 1 \mathrm{~kg} / \mathrm{m}^{3}, 10^{-5} \mathrm{~Pa}, 0\right)$ and a particle resolution of $\Delta x=r / 8 \mathrm{~m}$.

Figure  18  Sketch of the computational model

Download: Full-Size Img

Figure  19  Sketch of the corresponding particle distribution in the 2D axisymmetric plane

Download: Full-Size Img

Figure 20 depicts the particle distributions at four different moments, illustrating the propagation process of a pressure shock wave. Notably, to make the propagation process of the pressure shock wave clearer, the particle distribution on the r-axis is enlarged in the region of z ∈[− 0.06 m, 0.06 m]. The explosion shock wave spreads swiftly to the left during the initial pressure disturbance, and the adaptive multiresolution SPH accurately captures the pressure wave, resulting in a high-resolution region in front of and behind the wave. However, because of the high compression ratio and significant deformation, the overall particle distribution is not particularly uniform. Moreover, the particle spacing at the interface is large because of the secondary reflection of the pressure wave at the multiphase interface.

Figure  20  Particle distributions and pressure waves at different times

Download: Full-Size Img

Figure 21 shows the velocity, pressure, and density distributions at t = 0.05 s. As shown in Figure 22, the velocity, pressure, and density distributions on the z = 0 plane obtained by the present axisymmetric AMRS SPH model are compared with the results obtained by Sedov et al. (2018), Sigalotti et al. (2006), Fu and Ji (2019), and Fang et al. (2022). The results obtained by the present AMRS SPH model have a certain error compared with the results obtained using the reference solution (Sedov, 2018), especially in terms of pressure, where the interface pressure oscillations are obvious. The reason for this finding is that, for the strong compressible multiphase flow with strong impact, a larger artificial viscosity coefficient is usually used to stabilize the numerical method and prevent particle penetration, resulting in numerical dissipation. In addition, this study does not use a method similar to Riemann SPH to deal with the interface pressure. We will conduct an in-depth discussion and analysis on this in subsequent research.

Figure  21  Pressure, density, and velocity distributions at t = 0.05 s

Download: Full-Size Img

Figure  22  Pressure, density, and velocity distributions on the z = 0 plane

Download: Full-Size Img

Figure 23 shows the time histories of the number of particles in Fluids 1 and 2 with time. The number of particles increases from 492 to 1 353 as Fluid 1 expands. However, in the region t ∈ [0, 0.01 s], the increase is modest because the high-resolution region of the pressure shock wave is transferred from Fluid 1 to Fluid 2. At the same time, the number of particles in Fluid 2 decreases from 3 608 to 2 173 because of compression. However, in the region t ∈ [0, 0.01 s], the high-resolution region of the pressure shock wave shifts from Fluid 1 to Fluid 2, resulting in a gradual reduction in the number of particles. Figure 24 shows the time histories of Fluid 1, Fluid 2, and the total number of particles. For the total number of particles, the variation range is relatively small throughout the simulation process. The maximum number of particles is 4 100, and the minimum number of particles is 3 526.

Figure  23  Time histories of the number of particles

Download: Full-Size Img

Figure  24  Time histories of the number of particles

Download: Full-Size Img

This scenario demonstrates that the current AMRS effectively forms refined high resolutions in the local region surrounding a pressure shock wave and that the current axisymmetric AMRS SPH accurately and effectively deals with significant discontinuity challenges.

4.4   Underwater explosion problem

To further verify the effectiveness and accuracy of the present AMRS SPH model in handling multiple specific high-resolution regions, the underwater explosion test is utilized. The local regions of the multiphase interface and pressure shock wave are set as high-resolution regions. Underwater explosions are widely used in practical engineering and have always been a focal point in the field of marine engineering research. However, experimental investigation and numerical modeling are extremely challenging because of the complexity of the nonlinear properties of underwater explosions, which include massive deformation and extreme discontinuities. The SPH approach has been proven to be a viable solution for this issue. However, current SPH simulations of underwater explosions are mostly conducted in a 2D framework because of the high computational cost of 3D models. Establishing a 2D axisymmetric underwater explosion model significantly reduces the computational expense of 3D simulation. The widely utilized empirical solution of Zamyshlyaev and Yakovlev (1973) is compared with the simulation results.

The computer model for underwater explosions is illustrated in Figure 25. The detonation radius of cylindrical TNT explosive is r = 0.018 m, and $\Delta x$ = 0.000 9 m is the particle resolution setting at the center of the cylindrical water tank, which has a height of H = 0.18 m and a radius of R = 0.45 m. Two monitoring stations are placed at 6r and 8r radians away from the explosive. In this test, the explosive gas created by the TNT explosion is solved using the JWL equation of state, whereas the water pressure is solved using the Mie–Gruneisen equation of state.

Figure  25  Sketch of the computational model of underwater explosions

Download: Full-Size Img

Figure 26 shows the particle distributions and the changes of the multiphase flow interface at four different times during underwater explosions. Notably, the high-pressure explosive gas rapidly spreads outward after the explosion, constantly squeezing the surrounding water particles. Because of the low compressibility of water, the high-pressure explosive gas encounters significant resistance when expanding outward, causing the interface to move slowly. Under the action of the adaptive multiresolution method, a high-resolution region is formed near the interface of multiphase flow, and this region changes continuously with the movement of the multiphase interface. Figure 27 shows the particle pressure distributions at four different times during underwater explosions, illustrating the propagation process of the pressure shock wave produced by the explosion. After the explosion, the pressure shock wave produced by the explosion propagates rapidly to the left. Under the action of the adaptive multiresolution SPH technology, the pressure wave is accurately and dynamically captured, forming a high-resolution region at the front of and behind the wave.

Figure  26  Particle distributions at different times

Download: Full-Size Img

Figure  27  Pressure waves at different times

Download: Full-Size Img

The pressure distributions of underwater explosions obtained by the present SPH method are shown in Figure 28. Notably, the pressure shock wave propagates rapidly to the left. However, because of rapid attenuation under the hindrance of water, the peak pressure decreases rapidly. The change in pressure at the detection point over time is shown in Figure 29. The results obtained by the present AMRS SPH model have a certain error compared with the results obtained using the empirical solution of Zamyshlyaev and Yakovlev (1973), especially in terms of pressure attenuation and increase. The reason for this finding is that, for the strong compressible multiphase flow with strong impact, a larger artificial viscosity coefficient is usually used to stabilize the numerical method and prevent particle penetration, resulting in numerical dissipation. In addition, this study does not use a method similar to Riemann SPH to deal with the interface pressure. We will conduct an in-depth discussion and analysis on this in subsequent research. The temporal histories of the particle counts of water and TNT are shown in Figure 30.

Figure  28  Pressure distributions at different times

Download: Full-Size Img

Figure  29  Pressure versus time of the measuring points

Download: Full-Size Img

Figure  30  Time histories of the number of particles

Download: Full-Size Img

This case further demonstrates that the present AMRS SPH effectively forms multiple local high-resolution regions for the multiphase interface and pressure shock wave to improve computational accuracy. The present axisymmetric AMRS SPH accurately and effectively deals with large discontinuity problems.

In this study, some classic strongly compressible multiphase flow problems are simulated based on the AMRS SPH model. Bubble dynamics has always been the focus of most scholars in multiphase flow. Plesset and Prosperetti (1977) provided a detailed introduction to the principle of bubble dynamics. Recently, Zhang et al. (2023) established a new oscillating bubble dynamics theory and proposed the bubble dynamics equation for the first time. The SPH has special advantages for simulating the bubble dynamics. The AMRS SPH model developed in this study is also suitable for simulating bubble dynamics, and a high-resolution region at the bubble interface can be formed to improve the computational accuracy of the interface. Moreover, the bubble dynamics will be investigated using the AMRS SPH model in our future work.

5   Conclusion

In this study, an axisymmetric adaptive multiresolution SPH model is developed to solve the strongly compressible multiphase flow problem. The difference in fluid volume for the strongly compressible multiphase flow problem is significant, which leads to an uneven particle distribution. In the axisymmetric AMRS SPH model, an IVAS is proposed to maintain the stability of particle volume and the uniformity of particle distribution. Moreover, an AMRS, the IVAS technique combined with the APR technique, for strongly compressible multiphase flow is developed, which realizes high-precision calculation by dynamic forming high-resolution local regions in the problem domain.

Four typical strongly compressible multiphase flow problems are investigated. Comparing the SPH results with the results from other resources, we can conclude that the IVAS is suitable for solving the problem of uneven particle distribution caused by excessive differences in particle volume. The AMRS, the IVAS technique combined with the APR technique, effectively and dynamically forms a high-resolution local region to improve the computational accuracy, which effectively reduces some nonphysical oscillations in simulating strongly compressible multiphase flow problems. The adaptive multiresolution SPH model effectively models strongly compressible multiphase flows.

However, the present AMRS SPH model still has certain flaws in dealing with strongly compressible multiphase flow problems with high compression ratios and strong impact properties. The present AMRS SPH model also has the problem of numerical dissipation, which causes certain errors in the results. The present AMRS SPH model will be combined with Riemann SPH to improve the computational accuracy and efficiency in future work.