Numerical modelling of shock-bubble interactions using a pressure-based algorithm without Riemann solvers

TL;DR

This paper evaluates the accuracy of fully-coupled pressure-based algorithms without Riemann solvers in simulating shock-bubble interactions across different flow configurations, highlighting their effectiveness and limitations.

Contribution

It demonstrates the predictive capabilities of pressure-based algorithms without Riemann solvers for shock-bubble interactions in gas-gas and liquid-gas flows, with insights on mesh resolution effects.

Findings

Mesh resolution has minor influence in gas-gas flows.

Mesh resolution critically affects liquid-gas bubble behavior.

Pressure-based algorithms accurately model shock-bubble interactions.

Abstract

The interaction of a shock wave with a bubble features in many engineering and emerging technological applications, and has been used widely to test new numerical methods for compressible interfacial flows. Recently, density-based algorithms with pressure-correction methods as well as fully-coupled pressure-based algorithms have been established as promising alternatives to classical density-based algorithms based on Riemann solvers. The current paper investigates the predictive accuracy of fully-coupled pressure-based algorithms without Riemann solvers in modelling the interaction of shock waves with one-dimensional and two-dimensional bubbles in gas-gas and liquid-gas flows. For a gas bubble suspended in another gas, the mesh resolution and the applied advection schemes are found to only have a minor influence on the bubble shape and position, as well as the behaviour of the dominant shock waves and rarefaction fans. For a gas bubble suspended in a liquid, however, the mesh resolution has a critical influence on the shape, the position and the post-shock evolution of the bubble, as well as the pressure and temperature distribution.

Equations

$$\frac{\partial \rho }{\partial t}+\frac{\partial \rho u_i}{\partial x_i}=0,$$

$$\frac{\partial \rho u_j}{\partial t}+\frac{\partial \rho u_i{u}_j}{\partial x_i}=-\frac{\partial p}{\partial x_j},$$

$$\frac{\partial \rho h}{\partial t}+\frac{\partial \rho u_{i}h}{\partial x_i}=\frac{\partial p}{\partial t},$$

$$\rho =\frac{p+{\gamma }_0{\Pi }_0}{R_{0}T},$$

$$a=\sqrt{{\gamma }_0\frac{p+{\Pi }_0}{\rho }}=\sqrt{({\gamma }_0-1)c_{p}T}$$

$$c_p=c_{p,0}\frac{p+{\Pi }_0}{p+{\gamma }_0{\Pi }_0}.$$

$$\psi (\mathbf{x})=\{\begin{array}{ll}0& \text{if}\mathbf{x}\in {\Omega }_{\text{a}}\\ 1& \text{if}\mathbf{x}\in {\Omega }_{\text{b}},\end{array}$$

$$\frac{\partial \psi }{\partial t}+u_i\frac{\partial \psi }{\partial x_i}=0.$$

$${\iiint }_{V_P}{\frac{\partial \varphi }{\partial t}|}_{P}dV\approx \frac{{\varphi }_P-{\varphi }_P^{(t-\Delta t)}}{\Delta t}{V}_P,$$

$${\iiint }_{V_P}{\frac{\partial \varphi }{\partial t}|}_{P}dV\approx \frac{3{\varphi }_P^{(t)}-4{\varphi }_P^{(t-\Delta t)}+{\varphi }_P^{(t-2\Delta t)}}{2\Delta t}{V}_P,$$

$${\iiint }_{V_P}{\frac{\partial \rho u_i\varphi }{\partial x_i}|}_P\text{d}V\approx \sum _f{\tilde{\rho }}_f{\vartheta }_f{\tilde{\varphi }}_f{A}_f,$$

$${\tilde{\varphi }}_f={\varphi }_U+\frac{{\xi }_f}{2}({\varphi }_D-{\varphi }_U),$$

$${\vartheta }_f={\overline{u}}_{i,f}{n}_{i,f}-{\hat{d}}_f\left[{\frac{\partial p}{\partial x_i}|}_f-\frac{{\rho }_f^{\ast }}{2}\left({\frac{1}{{\rho }_P}\frac{\partial p}{\partial x_i}|}_P+{\frac{1}{{\rho }_Q}\frac{\partial p}{\partial x_i}|}_Q\right)\right]n_{i,f}+{\hat{d}}_f\frac{{\rho }_f^{\ast (t-\Delta t_1)}}{\Delta t}\left({\vartheta }_f^{(t-\Delta t)}-{\overline{u}}_{i,f}^{(t-\Delta t)}{n}_{i,f}\right),$$

$${\frac{\partial p}{\partial x_i}|}_f{n}_{i,f}\approx \frac{p_Q-p_P}{\Delta x}.$$

$$\frac{{\rho }_P^{(n+1)}-{\rho }_P^{(t-\Delta t)}}{\Delta t}{V}_P+\sum _f{\tilde{\rho }}_f^{(n)}{\vartheta }_f^{(n+1)}+{\tilde{\rho }}_f^{(n+1)}{\vartheta }_f^{(n)}-{\tilde{\rho }}_f^{(n)}{\vartheta }_f^{(n)}{A}_f=0,$$

$$\begin{array}{rl}{\vartheta }_f^{(n+1)}& ={\overline{u}}_{i,f}^{(n+1)}{n}_{i,f}-{\hat{d}}_f\left[\frac{p_Q^{(n+1)}-p_P^{(n+1)}}{\Delta x}-\frac{{\rho }_f^{\ast }}{2}\left({\frac{1}{{\rho }_P}\frac{\partial p}{\partial x_i}|}_P^{(n)}+{\frac{1}{{\rho }_Q}\frac{\partial p}{\partial x_i}|}_Q^{(n)}\right)n_{i,f}\right]\\ & +{\hat{d}}_f\frac{{\rho }_f^{\ast (t-\Delta t)}}{\Delta t}\left({\vartheta }_f^{(t-\Delta t)}-{\overline{u}}_{i,f}^{(t-\Delta t)}{n}_{i,f}\right)\end{array}$$

$${\rho }^{(n+1)}=\frac{p^{(n+1)}+{\gamma }_0{\Pi }_0}{R_{0}T}.$$

$$\begin{array}{rl}& \frac{{\rho }_P^{(n)}{u}_{j,P}^{(n+1)}+{\rho }_P^{(n+1)}{u}_{j,P}^{(n)}-{\rho }_P^{(n)}{u}_{j,P}^{(n)}-{\rho }_P^{(t-\Delta t)}{u}_{j,P}^{(t-\Delta t)}}{\Delta t}{V}_P\\ & +\sum _f\left({\tilde{\rho }}_f^{(n)}{\vartheta }_f^{(n)}{\tilde{u}}_{j,f}^{(n+1)}+{\tilde{\rho }}_f^{(n)}{\vartheta }_f^{(n+1)}{\tilde{u}}_{j,f}^{(n)}+{\tilde{\rho }}_f^{(n+1)}{\vartheta }_f^{(n)}{\tilde{u}}_{j,f}^{(n)}-2{\tilde{\rho }}_f^{(n)}{\vartheta }_f^{(n)}{\tilde{u}}_{j,f}^{(n)}\right)A_f=-\sum _f{\overline{p}}_f^{(n+1)}{n}_{j,f}{A}_f\end{array}$$

$$\begin{array}{rl}& \frac{{\rho }_P^{(n)}{h}_P^{(n+1)}+{\rho }_P^{(n+1)}{h}_P^{(n)}-{\rho }_P^{(n)}{h}_P^{(n)}-{\rho }_P^{(t-\Delta t)}{h}_P^{(t-\Delta t)}}{\Delta t}{V}_P\\ & +\sum _f\left({\tilde{\rho }}_f^{(n)}{\vartheta }_f^{(n)}{\tilde{h}}_f^{(n+1)}+{\tilde{\rho }}_f^{(n)}{\vartheta }_f^{(n+1)}{\tilde{h}}_f^{(n)}+{\tilde{\rho }}_f^{(n+1)}{\vartheta }_f^{(n)}{\tilde{h}}_f^{(n)}-2{\tilde{\rho }}_f^{(n)}{\vartheta }_f^{(n)}{\tilde{h}}_f^{(n)}\right)A_f=\frac{p_P^{(n+1)}-p_P^{(t-\Delta t)}}{\Delta t}{V}_P.\end{array}$$

$$\frac{\partial \psi }{\partial t}+\frac{\partial u_i\psi }{\partial x_i}-\psi \frac{\partial u_i}{\partial x_i}=0.$$

$$\frac{{\psi }_P-{\psi }_P^{(t-\Delta t_{\psi })}}{\Delta t_{\psi }}{V}_P+\sum _f\frac{{\psi }_f+{\psi }_f^{(t-\Delta t_{\psi })}}{2}{\vartheta }_f{A}_f-\frac{{\psi }_P+{\psi }_P^{(t-\Delta t_{\psi })}}{2}\sum _f{\vartheta }_f{A}_f=0,$$

$$\rho =(1-\psi ){\rho }_{\text{a}}+\psi {\rho }_{\text{b}},$$

$$\frac{1}{\gamma -1}=\frac{1-\psi }{{\gamma }_{0,\text{a}}-1}+\frac{\psi }{{\gamma }_{0,\text{b}}-1}.$$

$$c_p=\frac{(1-\psi ){\rho }_{\text{a}}{c}_{p,\text{a}}+\psi {\rho }_{\text{b}}{c}_{p,\text{b}}}{\rho },$$

$${\tilde{\rho }}_f={\rho }_U^{\star }+\frac{{\xi }_f}{2}\left({\rho }_D^{\star }-{\rho }_U^{\star }\right).$$

$${\rho }_U^{\star }={\rho }_{\text{a},U}+{\psi }_P\left({\rho }_{\text{b},U}-{\rho }_{\text{a},U}\right)$$

$${\rho }_D^{\star }={\rho }_{\text{a},D}+{\psi }_P\left({\rho }_{\text{b},D}-{\rho }_{\text{a},D}\right).$$

$${\rho }_P^{(t-\Delta t)}={\rho }_{\text{a},P}^{(t-\Delta t)}+{\psi }_P\left({\rho }_{\text{b},P}^{(t-\Delta t)}-{\rho }_{\text{a},P}^{(t-\Delta t)}\right)$$

$${\tilde{h}}_f=\frac{1}{{\tilde{\rho }}_f}\left[{\rho }_U^{\star }{h}_U^{\star }+\frac{{\xi }_f}{2}\left({\rho }_D^{\star }{h}_D^{\star }-{\rho }_U^{\star }{h}_U^{\star }\right)\right],$$

$$h_U^{\star }$$

Full text

\geometry

textheight=25.7cm, textwidth=17.0cm

1 Introduction

The interaction of a shock wave with a bubble is a process of broad academic and engineering interest, with applications in combustion and detonation [Michael and Nikiforakis, 2019], medical applications, e.g. shock-wave lithotripsy [Johnsen, 2007, Pan et al., 2018], in geophysics [Delale, 2013] and in microfluidics [Ando et al., 2012, Ohl and Ohl, 2016], featuring a rich variety of fluid dynamic and thermodynamic phenomena, such as compression and expansion waves, strong local heating, cavitation, evaporation and condensation, as well as the production of vorticity and turbulence [Ranjan et al., 2011, Delale, 2013]. Shock-bubble interaction in a gas-gas flow is also used to study the interaction of shocks with gas inhomogeneities, whereas the shock-induced collapse of gas bubbles in liquids is of direct relevance to cavitating flows [Delale, 2013, Ohl and Ohl, 2016, Fuster, 2019]. As a result, the fluid dynamics of shock-bubble interactions have been studied extensively, both experimentally [Haas and Sturtevant, 1987, Layes et al., 2003, 2005, Ranjan et al., 2007, Zhai et al., 2011] and computationally [Quirk and Karni, 1996, Bagabir and Drikakis, 2001, Johnsen, 2007, Niederhaus et al., 2008b, a, Johnsen and Colonius, 2009, Zhai et al., 2011, Hejazialhosseini et al., 2013, Xiang and Wang, 2017, Pan et al., 2018, Yoo and Sung, 2018, Michael and Nikiforakis, 2019]. In addition, shock-bubble interactions have also been widely used as a canonical reference system to test and scrutinise new numerical schemes, see e.g. [Saurel and Abgrall, 1999, Allaire et al., 2002, Hu and Khoo, 2004, Nourgaliev et al., 2006, Chang and Liou, 2007, Terashima and Tryggvason, 2009, Kokh and Lagoutière, 2010, Shukla et al., 2010, Shukla, 2014, Wong and Lele, 2017, Denner et al., 2018].

Computational fluid dynamics has assumed an increasingly prominent role for the study and analysis of compressible interfacial flows, and especially for the study of shock-bubble interactions, over the past decades, as a result of rapidly advancing developments of the relevant numerical algorithms as well as the substantial computational resources routinely available nowadays. The numerical modelling of compressible interfacial flows thereby requires a consistent numerical treatment of the fluid interface that retains the main features of the solution, in particular the propagation of pressure waves [Abgrall and Saurel, 2003, Coralic and Colonius, 2014, Denner et al., 2018]. However, the typically sharp change in Mach number at the fluid interface and the associated change in dominant physical mechanisms lead to distinct, and often contrasting, numerical requirements, which complicate the accurate and robust modelling of compressible fluid phenomena, such as the interaction of shock waves with gas bubbles. For instance, while the numerical algorithm has to ensure that the density is independent of the pressure and recovers a divergence-free velocity field in the incompressible limit [Chorin and Marsden, 1993, Hauke and Hughes, 1998], the accurate prediction of shock waves requires a conservative discretisation of the governing conservation laws [Hou and Floch, 1994].

Contemporary numerical methods for compressible interfacial flows are typically predicated on density-based algorithms, where the governing conservation equations are solved for the density, momentum and total energy of the flow [Baer and Nunziato, 1986, Allaire et al., 2002, Murrone and Guillard, 2005, Coralic and Colonius, 2014]. In these models, an exact or approximate Riemann solver is usually applied to evaluate the fluxes, with HLLC-type solvers [Toro et al., 1994] having gained particular popularity for interfacial flows [Shyue, 2006, Tokareva and Toro, 2010, Tian et al., 2011, Coralic and Colonius, 2014]. The Ghost-Fluid method (GFM) [Fedkiw et al., 1999a] has established itself as a promising alternative to solving a Riemann problem [Fedkiw et al., 1999b, Terashima and Tryggvason, 2009, Bo and Grove, 2014], with recent extensions to improve the stability of simulations with strong shock-interface interactions and compressible gas-liquid flows [Liu et al., 2003, Wang et al., 2006, Liu and Hu, 2017]. While density-based algorithms are naturally suited for compressible flows, they are poorly suited for low-Mach number flows [Chorin, 1967, Karimian and Schneider, 1994, Wesseling, 2001, Cordier et al., 2012], where the coupling of pressure and density vanishes. Traditionally, tailored pre-conditioning techniques have been applied to extend density-based methods to low-Mach number flows [Turkel et al., 1993, Turkel, 2006], which are however computationally very expensive for transient problems. This has been motivating recent work on combining density-based methods with segregated pressure-correction algorithms [Xiao, 2004, Moguen et al., 2012, Fuster and Popinet, 2018, Moguen et al., 2019] and hybrid density/pressure-based algorithms [van der Heul et al., 2003, Park and Munz, 2005], in which the continuity equation is solved for density but the energy equation is reformulated as an equation for pressure. An additional difficulty for interfacial flows associated with density-based methods is that the pressure field has to be reconstructed based on the applied thermodynamic closure model, which has proved to be a considerable difficulty in interfacial cells where two bulk phases coexist [Abgrall and Karni, 2001, Allaire et al., 2002, Murrone and Guillard, 2005].

Pressure-based algorithms for compressible flows, in which the continuity equation is solved for pressure, are less prominent than their density-based counterparts. Deriving stable and efficient numerical schemes for the transonic regime, and formulating consistent shock-capturing schemes, is known to be difficult for pressure-based algorithms [Wesseling, 2001]. However, because pressure plays an important role in all Mach number regimes, i.e. the pressure-velocity coupling dominates at low Mach numbers and the pressure-density coupling dominantes at high Mach numbers [Van Doormaal et al., 1987, Moukalled et al., 2016], pressure-based algorithms potentially offer a distinct advantage for applications in all Mach number regimes or in which the Mach number varies strongly, such as interfacial flows. Although, starting with the seminal work of Harlow and Amsden [1971b], a variety of pressure-based algorithms for compressible single-phase flows has been proposed, notably [Van Doormaal et al., 1987, Demirdžić et al., 1993, Karimian and Schneider, 1994, Xiao et al., 2017], it was only recently that Denner et al. [2018] proposed a conservative pressure-based algorithm for compressible interfacial flows at all speeds. This algorithm was proposed in conjunction with a new interface discretisation, the acoustically-conservative interface discretisation (ACID), that retains the acoustic features of the flow, which facilitates a rational definition of fluid properties in interfacial cells and which does not require Riemann solvers to compute the fluxes through the fluid interface. Denner et al. [2018] showed that such an algorithm yields unique definitions of the speed of sound and the Rankine-Hugoniot conditions in the interface region, and demonstrated a reliable prediction of acoustic and shock waves even in interfacial flows with acoustic impedance matching and shock impedance matching.

To model shock-bubble interactions, density-based algorithms in conjunction with Riemann-type solvers have been applied in most published studies to date [Quirk and Karni, 1996, Bagabir and Drikakis, 2001, Nourgaliev et al., 2006, Johnsen, 2007, Niederhaus et al., 2008b, Hejazialhosseini et al., 2013, Xiang and Wang, 2017]. While exact Riemann solvers are prohibitively time consuming, approximate Riemann solvers require an a priori approximation of the characteristic wave speeds, which ensues a substantial complexity of the numerical algorithms [Saurel and Pantano, 2018]. Thereby a strong dependence of the solution on the spatial resolution and the applied numerical schemes has been generally observed. For instance, a distinct feature of the shock-bubble interaction in gas-gas systems observed in numerical simulations of shock-bubble interactions, are instabilities (cf. Richtmyer-Meshkov instability [Brouillette, 2002]) forming on the interface as the shock wave passes. Noticeably, the shape and evolution of these interface instabilities depend strongly on the applied numerical methods [Johnsen and Colonius, 2006, Denner et al., 2018, Saurel and Pantano, 2018], e.g. the interface treatment or the advection schemes. In addition, these interface instabilities feature ever smaller lengthscales with an increasing spatial resolution of the simulation [Wong and Lele, 2017, Denner et al., 2018], although this is to be expected if surface tension, viscous stresses and heat conduction are neglected, a common assumption, supported by experimental observations [Layes et al., 2003, Zhai et al., 2011], with reference to the small timescales considered in typical numerical simulations and the associated marginal influence of these effects, since there is no physical means to regulate or dissipate the small-scale flow features. However, the influence of the spatial resolution of the computational mesh and of the choice of discretisation schemes on the predictive quality of the modelling of shock-bubble interactions have not yet been studied comprehensively, especially in the context of pressure-based algorithms without Riemann solvers. It is, moreover, as yet unclear what influence such differences in interface instabilities have on the development and evolution of shock waves and rarefaction fans during and after the shock-bubble interaction, which is important for many of the associated engineering applications, especially in medical and bioengineering applications.

This article investigates the modelling of shock-bubble interactions using the pressure-based algorithm proposed by Denner et al. [2018], where the fluxes are evaluated with the ACID method and no Riemann solvers are applied. The aim is to identify the minimum spatial resolution requirements for a converged solution with respect to the primary flow quantities, as well as the influence of the discretisation scheme and of interface instabilities on the predictive accuracy of the main flow features, for shock-bubble interactions in both gas-gas and liquid-gas flows. As test-cases the shock interaction with a one-dimensional helium-bubble in air, a one-dimensional air-bubble in water, a two-dimensional R22-bubble in air and a two-dimensional air-bubble in water are considered. While the shock-bubble interaction in gas-gas flows is not very sensitive to the employed discretisation schemes or the resolution of the computational mesh, the presented results demonstrate a very strong dependency of the primary flow quantities, especially temperature, on the spatial resolution of the computational mesh during the interaction of a shock wave with an air bubble suspended in water.

The governing equations are introduced in Section 2 and the numerical framework is presented in Section 3. The results of this study are presented and discussed in Section 4, and conclusions are drawn in Section 5.

2 Governing equations

The conservation laws governing fluid flow at all speeds, assuming viscous stresses and heat conduction are neglected, are the Euler equations, consisting of the conservation of mass

[TABLE]

the conversation of momentum

[TABLE]

and the conservation of energy

[TABLE]

where $t$ is time, $\boldsymbol{u}$ is the velocity vector, $p$ is pressure, $\rho$ is the density and $h=c_{p}\,T+\boldsymbol{u}^{2}/2$ is the specific total enthalpy, with $c_{p}$ the specific isobaric heat capacity and $T$ the temperature. Gravity and surface tension are neglected in this study.

The stiffened-gas model [Harlow and Amsden, 1971a, Saurel et al., 2007] is applied to define the thermodynamic properties of the fluid and close the governing conservation laws. The density-pressure relationship is defined by the stiffened-gas equation of state (EOS)

[TABLE]

where $\Pi_{0}$ is a fluid-dependent pressure constant, $R_{0}=c_{p,0}-c_{v,0}$ is the specific heat capacity and $\gamma_{0}=c_{p,0}/c_{v,0}$ is the heat capacity ratio, with the reference specific isobaric heat capacity $c_{p,0}$ and the reference specific isochoric heat capacity $c_{v,0}$. The speed of sound is given as

[TABLE]

and the specific isobaric heat capacity is [Denner et al., 2018]

[TABLE]

For $\Pi_{0}=0$, the stiffened-gas EOS reverts to the ideal-gas EOS, and the fluid is calorically perfect with $c_{p}=c_{p,0}$.

The Volume-of-Fluid (VOF) method [Hirt and Nichols, 1981] is adopted to capture the fluid interface between two immiscible bulk phases. To this end, the VOF method applies a colour function field $\psi$, defined as

[TABLE]

where $\Omega_{\textup{a}}$ and $\Omega_{\textup{b}}$ are the subdomains occupied by fluid a and b, respectively, and $\Omega=\Omega_{\textup{a}}\cup\Omega_{\textup{b}}$ is the computational domain. The interface is located in every cell where $0<\psi<1$. Because the interface is a material front propagating with the flow [Denner et al., 2018], the colour function $\psi$ is advected with the underlying fluid velocity by the advection equation

[TABLE]

3 Numerical framework

The numerical framework is based on the fully-coupled pressure-based algorithm of Denner et al. [2018], which is predicated on a finite-volume discretisation with collocated variable arrangement.

3.1 Temporal and spatial discretisation

The First-Order Backward Euler scheme (BDF1) and the Second-Order Backward Euler scheme (BDF2) are used to discretise the transient terms of the governing flow equations. The BDF1 scheme applied to the integrated transient term of a general flow variable $\phi$ is given for cell $P$ as

[TABLE]

and the BDF2 scheme is defined as

[TABLE]

where $\Delta t$ is the time-step, superscript $(t-\Delta t$) denotes values of the previous time-level, superscript $(t-2\Delta t)$ denotes values of the previous-previous time-level and $V_{P}$ is the volume of mesh cell $P$. As previously suggested by Denner et al. [2018], for consistency all transient terms of the governing equations (1)-(3) are discretised with the same scheme.

Applying the divergence theorem, assuming the surface of the control volume has a finite number of flat faces $f$ and applying the midpoint rule, the discretised advection terms of Eqs. (1)-(3) are given as

[TABLE]

where $\vartheta_{f}=\boldsymbol{u}_{f}\cdot\boldsymbol{n}_{f}$ is the advecting velocity of face $f$ (see Section 3.2), $\boldsymbol{n}_{f}$ is the outward-pointing unit normal vector of face $f$ and $A_{f}$ is the area of face $f$. The advected variable $\phi$ at face $f$ is interpolated from the adjacent cell-centred values using a total variation diminishing (TVD) interpolation scheme [Denner and van Wachem, 2015], with which the face value is given as

[TABLE]

where subscripts $U$ and $D$ denote the upwind and downwind cells, respectively, and $\xi_{f}$ is the flux limiter. In this study, the first-order upwind scheme ($\xi_{f}=0$), the Minmod scheme and the Superbee scheme [Roe, 1986] are considered.

3.2 Advecting velocity

The momentum-weighted interpolation (MWI) is applied to define an advecting velocity $\vartheta_{f}=\boldsymbol{u}_{f}\cdot\boldsymbol{n}_{f}$ at cell faces, which is used in the discretised advection terms of the governing equations. Following the unified formulation of the MWI proposed by Bartholomew et al. [2018], the advecting velocity $\vartheta_{f}$ at face $f$ is given as

[TABLE]

where subscript $Q$ denotes the neighbour cell of $P$ adjacent to face $f$, the interpolated face velocities $\overline{\boldsymbol{u}}_{f}$ and $\overline{\boldsymbol{u}}_{f}^{(t-\Delta t)}$ are obtained by linear interpolation, and the pressure gradient normal to face $f$ is discretised as

[TABLE]

The face density $\rho_{f}^{\ast}$ is interpolated by a harmonic average [Bartholomew et al., 2018] and the coefficient $\hat{{d}}_{f}$ follows from the coefficients associated with the advection term (and shear stress term if viscosity is considered) of the momentum equations, as detailed in [Bartholomew et al., 2018].

MWI provides a robust pressure-velocity coupling for incompressible and low Mach number flows by applying a low-pass filter acting on the third derivative of pressure [Bartholomew et al., 2018], thus avoiding pressure-velocity decoupling due to the collocated variable arrangement. As a result of the additional terms required to ensure a robust pressure-velocity coupling, the MWI introduces an unphysical dissipation of kinetic energy, which however diminishes with $\Delta x^{3}$ and is independent of $\Delta t$ [Bartholomew et al., 2018]. The transient term of Eq. (13) ensures a time-step independent contribution of the MWI in conjunction with the coefficient $\hat{{d}}_{f}$ of the pressure term [Bartholomew et al., 2018] and including the transient term is important for a correct temporal evolution of pressure waves [Moguen et al., 2015, Bartholomew et al., 2018].

3.3 Discretised governing conservation laws

The discretised continuity equation (1) for cell $P$, applying the BDF1 scheme for clarity of presentation, is given as

[TABLE]

where the superscript $(n)$ denotes known values of the most recent available solution and superscript $(n+1)$ denotes quantities which are solved implicitly. The advection term is linearised by a Newton linearisation to facilitate a smooth transition from low to high Mach number regions [Karimian and Schneider, 1994, Kunz et al., 1999, Xiao et al., 2017]. Following previous studies [Denner et al., 2018, Denner, 2018], the semi-implicit formulation of the advecting velocity is given as

[TABLE]

and the pressure-implicit formulation of the density is given as

[TABLE]

The discretised momentum equations (2) of cell $P$ are given, with both the transient term and the advection term linearised by a Newton linearisation [Denner et al., 2018, Denner, 2018], as

[TABLE]

with $\rho^{(n+1)}_{P}$ given by Eq. (17) and $\vartheta^{(n+1)}_{f}$ given by Eq. (16). Similarly, the discretised energy equation (3) of cell $P$ is given as

[TABLE]

3.4 Interface advection

The VOF advection equation (8) is discretised using a compressive VOF method [Denner and van Wachem, 2014, Denner et al., 2018]. Following Denner et al. [2018], Eq. (8) is reformulated as

[TABLE]

Using the Crank-Nicolson scheme for the discretisation of the transient term, the semi-discretised form of Eq. (20) becomes

[TABLE]

where $\Delta t_{\psi}$ is the time-step applied to advect the colour function $\psi$. The advection of the colour function is discretised using the same advecting velocity $\vartheta_{f}$ as for all advection terms of the governing equations. The face value $\psi_{f}$ is interpolated using the CICSAM scheme [Ubbink and Issa, 1999], taking into account the orientation of the interface and the available flux volume. Excellent volume conservation has previously been demonstrated for this compressive VOF method for incompressible [Denner and van Wachem, 2014] and compressible flows [Denner et al., 2018].

3.5 Coupling of the bulk phases

The discretised governing equations presented in Section 3.3 are extended to interfacial flows using the acoustically-conservative interface discretisation (ACID) [Denner et al., 2018]. The ACID method assumes that, for the purpose of discretising the governing conservation laws for a given cell, all cells in its finite-volume stencil are assigned the same colour function value, i.e. the colour function is assumed to be constant in the entire finite-volume stencil. The relevant thermodynamic properties that are discontinuous at the interface, i.e. density and enthalpy, are evaluated based on the constant colour function field in the applied finite-volume stencil. This recovers the contact discontinuity associated with the interface [Anderson, 2003, Denner et al., 2018] and enables the application of the fully-conservative discretisation scheme presented in Section 3.3, identical to the one applied for single-phase flows.

3.5.1 Fluid properties

A hydrodynamic and thermodynamic consistent definition of the fluid properties for interfacial flows requires special consideration. The density of the fluid is defined based on the colour function $\psi$ and the densities of the bulk phases as

[TABLE]

where the partial densities $\rho_{\textup{a}}$ and $\rho_{\textup{b}}$ of the bulk phases are given by Eq. (4). This linear interpolation of the density is required for the discrete conservation of mass, momentum and energy and is equivalent to an isobaric closure assumption for compressible interfacial flows [Allaire et al., 2002, Shyue, 2006]. The heat capacity ratio also follows from the isobaric closure assumption as

[TABLE]

The specific isobaric heat capacity is defined by a mass-weighted interpolation [Denner et al., 2018], which is essential for the conservation of the total energy, given as

[TABLE]

where the partial densities $\rho_{\textup{a}}$ and $\rho_{\textup{b}}$ are given by Eq. (4), density $\rho$ is given by Eq. (22), and the partial specific isobaric heat capacities $c_{p,\textup{a}}$ and $c_{p,\textup{b}}$ are given by Eq. (6). As shown by Denner et al. [2018], the speed of sound is defined throughout the domain based on Eq. (5) as $a=\sqrt{(\gamma-1)\,c_{p}\,T}$, and the material-dependent pressure constant of the stiffened-gas model is given as $\Pi=[(\gamma-1)\,\rho\,c_{p}\,T/\gamma]-p$, with the density $\rho$ given by Eq. (22), the specific isobaric heat capacity $c_{p}$ given by Eq. (24), and $(\gamma-1)$ as well as $\gamma$ given by Eq. (23).

3.5.2 Density treatment

Under the assumption that the colour function $\psi$ is constant throughout the finite-volume stencil of cell $P$, the density interpolated to face $f$ from the adjacent cell centre is given as

[TABLE]

The density $\rho_{U}$ at the upwind cell $U$ and $\rho_{D}$ at the downwind cell $D$ are given based on the colour function value of cell $P$ by Eq. (22), so that

[TABLE]

and

[TABLE]

The density at previous time-levels is evaluated in a similar fashion based on the colour function value of cell $P$, with [Denner et al., 2018]

[TABLE]

and likewise for $\rho_{P}^{(t-2\Delta t)}$, if required.

3.5.3 Enthalpy treatment

The specific total enthalpy at face $f$ is given, again assuming the colour function $\psi$ is constant throughout the finite-volume stencil of cell $P$, as [Denner et al., 2018]

[TABLE]

with $\tilde{\rho}_{f}$ given by Eq. (25), where the specific total enthalpy of the upwind and downwind cells are given as

[TABLE]

respectively, $\rho_{U}^{\star}$ is given by Eq. (26) and $\rho_{D}^{\star}$ is given by Eq. (27). The specific isobaric heat capacities $c_{p,U}^{\star}$ and $c_{p,D}^{\star}$ are defined by Eq. (24) with $\psi_{P}$ as

[TABLE]

and

[TABLE]

Since the specific total enthalpy is a primary solution variable, a deferred correction approach as proposed by Denner et al. [2018] is applied to enforce Eq. (29).

The specific total enthalpy at the previous time-levels follow analogously as [Denner et al., 2018]

[TABLE]

with

[TABLE]

and likewise for $h_{P}^{(t-2\Delta t)}$ and $c_{p,P}^{\star,(t-2\Delta t)}$, if required.

3.6 Solution procedure

The discretised governing equations presented in Section 3.3 are solved simultaneously in a single linear system of equations [Denner et al., 2018, Denner, 2018], $\boldsymbol{A}\boldsymbol{\chi}=\boldsymbol{b}$, with $\boldsymbol{A}$ being the coefficient matrix of size $5N\times 5N$, $\boldsymbol{\chi}\equiv(\boldsymbol{u},p,h)^{T}$ is the solution vector of length $5N$ of the primary solution variables and $\boldsymbol{b}$ is the right-hand side vector containing all known contributions, where $N$ is the number of mesh cells of the three-dimensional computational mesh. The solution procedure performs nonlinear iterations in which the linear system of governing equations is solved using the Block-Jacobi preconditioner and the BiCGSTAB solver of the software library PETSc [Balay et al., 2017], as described in detail by Denner [2018].

4 Results

The presented results focus on the spatial resolution requirements and discretisation necessary for the accurate prediction of shock-bubble interactions using a pressure-based algorithm. As already comprehensively demonstrated by Denner et al. [2018], the applied numerical algorithm captures shock waves and rarefaction fans accurately in single-phase flows and interfacial flows, with a robust convergence under mesh refinement and a precise prediction of the Rankine-Hugoniot relations also in interfacial cells.

4.1 One-dimensional helium-bubble in air

The interaction of a shock wave travelling in air with Mach number $M_{\textup{s}}=1.1$ and interacting with a helium bubble in a one-dimensional domain with a length of $1\,\textup{m}$ is considered. The initial post-shock conditions (I) are

[TABLE]

and the pre-shock conditions (II) are

[TABLE]

Air is taken to have a heat capacity ratio of $\gamma_{\textup{0,Air}}=1.4$ and a specific gas constant of $R_{\textup{0,Air}}=288.0\,\textup{J}\,\textup{kg}^{-1}\,\textup{K}^{-1}$, and helium is assumed to have a heat capacity ratio of $\gamma_{\textup{0,He}}=1.648$ and a specific gas constant of $R_{\textup{0,He}}=1581.2\,\textup{J}\,\textup{kg}^{-1}\,\textup{K}^{-1}$. The shock is initially located at $x_{0}=0.3\,\textup{m}$, the helium bubble occupies the interval $0.5\leq x\leq 0.7$ and the applied time-step corresponds to a Courant number of $\textup{Co}=a_{\textup{II,He}}\Delta t/\Delta x=0.44$.

The results, shown in Fig. 1, are obtained on equidistant meshes with three mesh resolutions resolving the one-dimensional domain with $200$, $1000$ and $5000$ cells, which corresponds to a spatial resolution of $40$, $200$ and $1000$ cells for the initial length of the helium bubble, respectively. The distribution of pressure, velocity and temperature within the one-dimensional domain depends considerably on the applied mesh resolution. While this may be a problem for the accurate local prediction of pressure, velocity and temperature, as well as the related thermodynamic or chemical processes, the position of the primary shock wave as it travels through the bubble is, apart from the sharpness of the discontinuity, unaffected by the mesh resolution. In addition, the colour function $\psi$ and the density $\rho$ (both not shown) are in very good agreement on the different meshes. All quantities obtained with the finest mesh resolution, $1000$ cells for the initial length of the helium bubble, are in excellent agreement with the corresponding exact Riemann solution.

4.2 One-dimensional air-bubble in water

The interaction of a shock wave travelling in water with Mach number $M_{\textup{s}}=1.1$ and interacting with an air bubble in a one-dimensional domain with a length of $2\,\textup{m}$ is considered. The initial post-shock conditions (I) are

[TABLE]

and the pre-shock conditions (II) are

[TABLE]

Water is taken to have a heat capacity ratio of $\gamma_{\textup{0,Water}}=4.1$, a pressure constant of $\Pi_{\textup{0,Water}}=4.4\times 10^{8}\,\textup{Pa}$ and a specific gas constant of $R_{\textup{0,Water}}=6000\,\textup{J}\,\textup{kg}^{-1}\,\textup{K}^{-1}$, and air is taken to have a heat capacity ratio of $\gamma_{\textup{0,Air}}=1.4$, a pressure constant of $\Pi_{\textup{0,Air}}=0\,\textup{Pa}$ and a specific gas constant of $R_{\textup{0,Air}}=288.0\,\textup{J}\,\textup{kg}^{-1}\,\textup{K}^{-1}$. The shock is initially located at $x_{0}=1.1\,\textup{m}$, the air bubble occupies the interval $1.3\leq x\leq 1.5$ and the applied time-step corresponds to $\textup{Co}=a_{\textup{II,Water}}\Delta t/\Delta x=0.45$.

The results are obtained on equidistant meshes with three mesh resolutions, resolving the one-dimensional domain with $400$, $2000$ and $10000$ cells, which corresponds to a spatial resolution of $40$, $200$ and $1000$ cells for the initial length of the air bubble, respectively, as considered in the previous section for the helium bubble in air. As observed in Fig. 2, the density is principally in good agreement on all three meshes, whereas the temperature distribution inside the bubble appears to be especially sensitive to the mesh resolution, with visible differences between the results obtained on the coarsest mesh compared to the results obtained on the two meshes with higher resolution. Furthermore, the pressure upstream of the water-air interface after the shock wave has passed exhibits a considerable dependency on the mesh resolution, as seen in Fig. 3. The pressure is significantly underpredicted compared to the exact Riemann solution at both time-instances shown in Fig. 3, with an underprediction of approximately $1.5\times 10^{5}\,\textup{Pa}$ on the coarsest mesh and approximately $0.6\times 10^{5}\,\textup{Pa}$ on the finest mesh. Despite these differences in pressure and temperature, the position of the primary shock wave as it travels through the bubble is, apart from the sharpness of the discontinuity, still in very good agreement on the meshes corresponding to $200$ and $1000$ cells for the initial length of the air bubble. On the coarsest considered mesh, corresponding to $40$ cells for the initial length of the air bubble, however, the position of the shock wave has an offset in the downstream direction, which may be attributed to an overprediction of the speed of sound as the shock wave passes the interface.

4.3 Two-dimensional R22-bubble in air

The interaction of a shock wave with $M_{\textup{s}}=1.22$ in air with a circular R22 bubble is simulated, a shock-bubble interaction which has previously been studied experimentally [Haas and Sturtevant, 1987] and numerically [Quirk and Karni, 1996, Niederhaus et al., 2008b]. The computational setup is schematically illustrated in Fig. 4. The shock wave is initially situated at $x=0.17\,\textup{m}$ and travels from left to right at speed $u_{\textup{s}}$. The shock wave separates the post-shock region (I) and the pre-shock region (II), which are initialised with

[TABLE]

Air is taken to have a heat capacity ratio of $\gamma_{\textup{0,Air}}=1.4$ and a specific gas constant of $R_{\textup{0,Air}}=288.0\,\textup{J}\,\textup{kg}^{-1}\,\textup{K}^{-1}$, and R22 is assumed to have a heat capacity ratio of $\gamma_{\textup{0,R22}}=1.249$ and a specific gas constant of [Denner et al., 2018] $R_{\textup{0,R22}}=90.885\,\textup{J}\,\textup{kg}^{-1}\,\textup{K}^{-1}$. The applied computational mesh is equidistant and Cartesian, and the applied time-step corresponds to a Courant number of $\textup{Co}=a_{\textup{Air,II}}\Delta t/\Delta x=0.38$.

Figures 5 and 6 show the contours of the density gradient and the pressure distribution at different dimensionless times $\tau=t\,a_{\textup{R22,II}}/d_{0}$ for equidistant Cartesian meshes with a mesh resolution of $200$, $300$ and $500$ cells per initial bubble diameter $d_{0}$. While the overall shape as well as the position of the bubble predicted on the different meshes are largely the same, interface instabilities with smaller structures develop as the mesh resolution increases. As mentioned in the introduction, this is to be expected, yet a coherent and sufficiently accurate description of the magnitude and frequency with which these instabilities occur in reality is presently not available. These interface instabilities generate acoustic waves, as seen in Figs. 5 and 6, which however do not affect the position, shape and strength of the dominant flow structures, i.e. shock waves and rarefaction fans. In general it is noticeable in Figs. 5 and 6, that the overall impact of the mesh resolution on the observed flow features is minor, apart from the interface instabilities developing as a result of the passing shock wave and the resolution of the shock waves and rarefaction fans. This observation is supported by the density profiles along the $x$-axis (direction of travel of the primary shock wave) shown in Fig. 7, which exhibit very little differences for the three considered mesh resolutions.

Considering different TVD differencing schemes for the discretisation of advected variables (see Section 3.1) leads to very similar observations as mesh refinement; applying a more compressive differencing scheme facilitates and increases the generation of interface instabilities. Figure 8 shows the contours of the density gradient and the pressure distribution at dimensionless time $\tau=1.15$ on a mesh with a mesh resolution of $500$ cells per initial bubble diameter $d_{0}$, using (in order of increasing compression) the first-order upwind scheme, the Minmod scheme and the Superbee scheme. The interface advection is unaffected by this choice and identical for all these cases, discretised as described in Section 3.4. Applying different TVD advection schemes only influences the development of interface instabilities, while the position and overall shape of the R22 bubble is largely the same. The strong instabilities observed at the interface when the Superbee scheme is applied, and the ensuing acoustic waves, can be clearly observed in Fig. 8, yet it is also apparent that the different resolution of the discontinuities and the interface instabilities developing with the Minmod and Superbee schemes have very little influence on the position and strength of the dominant shock waves and rarefaction fans. The density profiles along the $x$-axis in Fig. 9 support this observation.

4.4 Two-dimensional air-bubble in water

The interaction of a shock wave with $M_{\textup{s}}=1.72$ in water with a circular air bubble is simulated, as considered previously in other studies [Nourgaliev et al., 2006, Shukla, 2014, Haimovich and Frankel, 2017, Goncalves et al., 2019]. The computational setup is schematically illustrated in Fig. 10. The shock is initially situated at $x=6.6\times 10^{-3}\,\textup{m}$ and travels from left to right at speed $u_{\textup{s}}$. The shock wave separates the post-shock region (I) and the pre-shock region (II), which are initialised with

[TABLE]

Water is taken to have a heat capacity ratio of $\gamma_{\textup{0,Water}}=4.1$, a pressure constant of $\Pi_{\textup{0,Water}}=4.4\times 10^{8}\,\textup{Pa}$ and a specific gas constant of $R_{\textup{0,Water}}=6000\,\textup{J}\,\textup{kg}^{-1}\,\textup{K}^{-1}$, and air is taken to have a heat capacity ratio of $\gamma_{\textup{0,Air}}=1.4$, a pressure constant of $\Pi_{\textup{0,Air}}=0\,\textup{Pa}$ and a specific gas constant of $R_{\textup{0,Air}}=288.0\,\textup{J}\,\textup{kg}^{-1}\,\textup{K}^{-1}$. The applied computational mesh is equidistant and Cartesian, and the applied time-step corresponds to a Courant number of $\textup{Co}=a_{\textup{Water,II}}\Delta t/\Delta x=0.11$.

Figures 11-13 show the contours of the pressure and temperature distribution of the water-air system at different times $t$ for equidistant Cartesian meshes with a mesh resolution of $200$, $400$ and $600$ cells per initial bubble diameter $d_{0}$. Contrary to the rather similar solutions obtained on different meshes for the gas-gas shock-bubble interaction in Section 4.3, the evolution of the shock-bubble interaction of the air bubble in water is strongly dependent on the mesh resolution. In particular the temperature distribution inside the air bubble exhibits distinct differences on the considered meshes, with generally higher temperatures predicted when the mesh resolution is increased. These differences are especially pronounced when the primary shock wave travels through the bubble, e.g. at $t=3.0\,\mu\textup{s}$, where the higher temperature appears to influence the position of the shock wave considerably, as seen in Fig. 14c. Despite these differences in temperature distribution and position of the shock wave at $t=3.0\,\mu\textup{s}$, which are much less pronounced in the pressure and density fields shown in Fig. 14, the results obtained on the meshes with $400$ and $600$ cells per initial bubble diameter $d_{0}$ are in reasonably good agreement at the later stages of the shock-bubble interaction, as seen in Figs. 15 and 16. In fact, similar observations with respect to the mesh dependency for the same shock-bubble interaction were recently reported by Shukla [2014] and Goncalves et al. [2019] using density-based methods. The mesh with $200$ cells per initial bubble diameter$d_{0}$, on the other hand, yields significantly different results compared to the meshes with higher resolution, which affects the position of the shock wave as well as the shape of the bubble, as evident by comparing Fig. 11c with Fig. 12c, and by the density profiles in Fig. 14b, 15b and 16b.

5 Conclusions

In the current paper, the numerical modelling of shock-bubble interactions using the pressure-based algorithm proposed by Denner et al. [2018], where the fluxes are evaluated with the ACID method and no Riemann solvers are applied, has been investigated. While shock-bubble interactions in gas-gas flows are largely of academic interest, the interaction of shock waves with a bubble suspended in a liquid is encountered in many different engineering and emerging technological applications, especially in microfluidics and medical applications. Of course, an accurate prediction of shock-bubble interactions is, therefore, a prerequisite for numerical methods to be utilised in research and process development pertaining to these applications.

The presented results have demonstrated a very strong dependency of the interaction of a shock wave with an air bubble suspended in water on the spatial resolution of the computational mesh. In particular, the temperature field has been found to exhibit large differences between different mesh resolutions, which also contributes to differences in the propagation of the primary shock wave. To this end, for the shock interaction with the air bubble in water, the presented results suggest a mesh resolution of at least $400-600$ cells per initial diameter to yield reasonably converged results. However, even though the position of the shock wave as well as the pressure, temperature and density fields have been found to yield agreement for such a mesh resolution at later stages of the shock-bubble interaction, i.e. after the shock-wave has passed the bubble, these quantities still exhibit considerable differences while the shock wave passes through the bubble. Considering the rapid and significant increase in pressure, temperature and density as the bubble is compressed when the shock wave passes, the accuracy of the ideal-gas model also warrants further study, since the evolution of the bubble collapse has been shown by the presented results to be strongly dependent on the quality and accuracy of the prediction of pressure and temperature.

Shock-bubble interactions also provide a convenient canonical reference system to test and scrutinise new numerical schemes for the simulation of compressible interfacial flows; shock-bubble interaction can be found in most publications that propose a new numerical scheme for compressible interfacial flows. These tests mostly focus on gas-gas flows, e.g. the R22 bubble in air also considered in this study, while the computationally more expensive and challenging shock-bubble interaction in liquid-gas flows is frequently neglected. However, the presented results show that the shock-bubble interaction in a gas-gas flow is not sensitive to the employed numerical methods and the spatial resolution of the computational mesh, contrary to the shock-bubble interaction in a liquid-gas flow. This puts the informative value of validating numerical schemes using gas-gas shock-bubble interactions into question and strongly suggests that the shock-bubble interaction in liquid-gas flows are generally better suited to scrutinise and compare numerical methods for compressible interfacial flows.