Conditional moment methods for polydisperse cavitating flows

A dilute suspension of dynamically evolving bubbles flow in a compressible liquid is considered. For simplicity, we assume no slip between the bubbles and the surrounding liquid and that the gas density is much smaller than the liquid density. Under these assumptions, the mixture-averaged form of the compressible flow equations are

where $\rho$, $\bm{u}$, $p$, and $E$ are the density, velocity vector, pressure, and total energy. Terms associated with the bubbles modify these quantities and transport them in space according to ensemble phase-averaging, discussed next.

The ensemble-averaged equations follow from Zhang and Prosperetti [8] and Bryngelson et al. [10]. The disperse phase has a void fraction $\alpha$ and is assumed to be a dilute ($\alpha\ll 1$) population of spherical bubbles. The bubbles are represented statistically via random variables $R$, $\dot{R}$, and $R_{o}$ corresponding to the instantaneous bubble radius, time derivative, and equilibrium bubble radius. Section 2.3 presents this bubble dynamics model in detail. The mixture-averaged pressure field is computed as

where $p_{bw}$ are the associated bubble wall pressure and $p_{\ell}(\bm{x},t)$ is the liquid pressure according to the stiffened-gas equation of state [37]:

The coefficients of (3) represent water, with specific heat ratio $\gamma_{\ell}=1+1/\Gamma_{\ell}=7.15$ and stiffness $\Pi_{\infty,\ell}=$356\text{\,}\mathrm{MPa}$$ [38].

The bubble number density per unit volume $n(\bm{x},t)$ is conserved in the absence of coalescence or breakup:

For the spherical bubbles considered here, $n$ is related to the void fraction $\alpha$ via

and thus the void fraction $\alpha(\bm{x},t)$ transports as

where the right-hand-side represents the change in void fraction due to bubble growth and collapse.

The over-barred terms appearing in (2), (5), and (6),

denote average quantities of the bubble dispersion. In particular, they are raw moments $\mu_{lmn}$ with respect to a bubble number density function $f(R,\dot{R},R_{o})$,

which are computed via the methods of section 3.

Even for spherical models, dozens of available models employ differing assumptions about the bubble contents and simplifications of the physics. These models range from a system of two or more ODEs up to a set of PDEs (in a radial coordinate and time) describing the balance of mass, momentum, and energy of each bubble. PDE-based models are deemed intractable (at present), and we focus on ODE-based modeling. We choose one of the simplest possible models to demonstrate our methodology but note that our framework can be extended to arbitrarily complex ODE-based models.

We assume the spherical bubbles are filled with noncondensible gas and that, insofar as their dynamics are concerned, the liquid may be assumed incompressible. We assume that the gas undergoes a polytropic process during compression and expansion. Under these assumptions, the bubble radius is governed by a Rayleigh–Plesslet-like equation

which is dimensionless via the reference bubble size $R_{o}^{\ast}$ and ambient liquid pressure $p_{0}$, and density $\rho_{0}$. The polytropic index is $\gamma$; in examples below we use $\gamma=1.4$. In (9), $C_{p}\equiv p_{0}/p_{\ell}$ is the forcing pressure ratio and Reynolds and Weber numbers correspond to viscous and surface tension effects as

where $S$ is the air–water surface tension coefficient and $\nu_{0}$ is the liquid kinematic viscosity. Under these conditions, the bubble wall pressure of (2) simplifies to $p_{bw}=\left(R_{o}/R\right)^{3\gamma}$ and the last moment of (7) reduces to $\mkern 1.5mu\overline{\mkern-1.5muR^{3}(R_{o}/R)^{3\gamma}\mkern-1.5mu}\mkern 1.5mu$.

Following the usual procedure, a finite set of raw moments $\vec{\boldsymbol{\mu}}$ represent $f$ per (8) [14]. The specific moments that make up $\vec{\boldsymbol{\mu}}$ depend on the inversion algorithm and are defined in the appendices A and B. These moments transport on the grid and evolve according to the bubble dynamics of section 2.3 as

where

and $\Omega=\Omega_{R}\times\Omega_{\dot{R}}\times\Omega_{R_{o}}=(0,\infty)\times(-\infty,\infty)\times(0,\infty)$. The integrand of (13) is closed via the bubble dynamics model (9) for $\ddot{R}$, and the integral is approximated via quadrature, which is discussed next.

Since $R_{o}$ is not a dynamic variable, the number density function is split as

The raw moments (8) are then

with $N_{R_{o}}$ time-independent weights $w_{i}$ and nodes $\widehat{R}_{o,i}$. Simpson’s rule computes these nodes and weights and the accuracy in the approximation of (16) depends on $N_{R_{o}}$. Other numerical integration methods are suitable and will be discussed in section 5.3. The $R_{o}$-conditioned moments are

The moment indices comprising the moment set of (18) are determined by the conditional quadrature moment method used to invert those moments for a set of quadrature points and weights (and the number of points desired) [33, 35]. In particular, $\mu_{lm}(\widehat{R}_{o,i})$ is traded for quadrature points $\{\widehat{R}_{j},\widehat{\dot{R}}_{k}\}(\widehat{R}_{o,i})$ and weights $\widehat{w}_{j,k}(\widehat{R}_{o,i})$ for each $i=1,\dots,N_{R_{o}}$ (with $j=1,\dots,N_{R}$; $k=1,\dots,N_{\dot{R}}$). The total moments (16) are then approximated by substituting (18) into (16) as

The CHyQMOM and CQMOM algorithms are described in appendices A and B. Their implementation is verified by ensuring that the error in closing a linear set of moment transport equations is comparable to a finite-precision round-off [39].

We consider the ability of the QBMMs to represent the statistics of the bubble dynamics. For this example, we take a monodisperse population with the same equilibrium radius $R_{o}=1$ (and thus $N_{R_{o}}=1$). $R$ and $\dot{R}$ are initialized with independent log-normal and normal distributions with variances $\sigma_{R}^{2}$ and $\sigma_{\dot{R}}^{2}$, respectively. The bubbles are forced by a step-change in pressure at $t=0^{+}$ from equilibrium to a value $C_{p}$ (which is varied). This represents a stringent test of model fidelity as high values of $C_{p}$ will induce strongly nonlinear bubble dynamics.

Figure 2 shows the evolution of the carried moment set for $C_{p}=0.3$. To assess the accuracy, the CHyQMOM results are compared against Monte-Carlo solutions that used $10^{4}$ samples to ensure that the sampling error was at least 10-times smaller than the QBMM model-form error. Thus the Monte-Carlo solutions can be regarded as exact solutions for these comparisons. The radial moments $\mu_{10}$ and $\mu_{20}$ show how the bubbles oscillate in response to the change in external pressure, eventually reaching a statistical equilibrium [40]. These oscillations are damped rapidly for smaller bubbles due to stronger viscous effects (smaller Re). Thus, larger bubbles entail larger model-form errors as they build over time. This is most prominent for the $\dot{R}$ moments, such as $\mu_{02}$, and smaller bubbles ((b) $R_{o}^{\ast}=$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$$) as the velocity statistics change more quickly than the radial ones.

A QBMM model-form relative $L_{2}$ error is

where Monte-Carlo (MC) simulations serve as surrogate truth data and $t_{i}$ are $N_{t}=10^{3}$ uniformly spaced times in the time interval $t\in[0,T]$. Figure 3 shows $\varepsilon$ for select first- and second-order moments $\mu$. We see that $\varepsilon$ is smaller for CHyQMOM and CQMOM than Gaussian closure, though the difference is modest and more strongly associated with $C_{p}$. The larger errors for smaller $C_{p}$ are associated with stronger bubble dynamics and the formation of non-Gaussian statistics, like skewness and kurtosis, that these closures do not represent [41]. One can represent higher-order statistics by carrying higher-order moments and thus inverting for a higher-order quadrature rule. However, this is computationally cumbersome and numerically unstable in the cases tested here.

Model accuracy was similar for all of the closure methods, but the computational time to solution is not. Figure 4 shows the total relative simulation time for the cases of figure 3 under the same adaptive time-step tolerance. CHyQMOM simulations are about 10-times cheaper than those using CQMOM or Gaussian closures, even though the accuracy is the same. This is attributed to larger time step sizes (given the same error constraint) and a smaller carried moment set than CQMOM. Due to its low cost for the same accuracy, the fully-coupled simulation algorithm employs CHyQMOM and is discussed next.

We initialize the flow as follows. The number density function $f(R,\dot{R},R_{o})$ is independently distributed with log-normal, normal, and log-normal shapes in the $R$, $\dot{R}$, and $R_{o}$ directions with expected values $\mathbb{E}[R]=\mathbb{E}[R_{o}]=1$ and $\mathbb{E}[\dot{R}]=0$ and shape parameters $\sigma_{\cdot}$.

The NDF $f$ initializes the moment set $\vec{\boldsymbol{\mu}}$ via integration. The primitive variables $\vec{\bm{q}}_{p}\equiv\{\rho,\bm{u},p,\alpha,\vec{\boldsymbol{\mu}}\}$ are initialized as appropriate. The bubble number density $n$ follows from (5), $\vec{\bm{q}}_{p}$, and $\mu_{300}$. Lastly, the known primitive variables are used to compute the conservative ones as $\vec{\bm{q}}_{c}\equiv\{\rho,\rho\bm{u},E,\alpha,n\vec{\boldsymbol{\mu}}\}$.

A fifth-order-accurate WENO [43] scheme reconstructs the primitive variables $\vec{\bm{q}}_{p}$ and the HLLC approximate Riemann solver [44] computes the fluxes. Algorithm 2 describes this process in detail. High-order WENO reconstructions do not guarantee that the reconstructed moments are realizable, though the moment sets remained invertible in the subsequent simulations.

The conservative variables are integrated in time using third-order-accurate SSP–RK3 time integration [45]. Once the spatial derivatives have been approximated, (21) becomes a semi-discrete system of ordinary differential equations in time. We treat the temporal derivative using a Runge–Kutta time-marching scheme for the state variables. To achieve high-order accuracy and avoid spurious oscillations, we use the third-order-accurate total variation diminishing scheme of Gottlieb and Shu [46]:

where superscripts $(1)$ and $(2)$ indicate intermediate time-step stages, $\bm{L}$ represents the portion of (21) that does not include the time derivative, and $\tilde{n}$ is the time-step index.

We assess the statistics of bubble dynamics in an idealized model problem consisting of an acoustically excited dilute bubble screen. The bubble screen parameterization matches that of Bryngelson et al. [10], with a length of $L_{s}=$5\text{\,}\mathrm{mm}$$, initial void fraction $\alpha_{o}=10^{-4}$, median bubble equilibrium size $R_{o}^{\ast}=$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$$ and log-normal variance $\sigma_{R_{o}}$. The domain is $L_{d}=5L_{s}$ long, and its boundaries are non-reflective via characteristic-based boundary conditions [47]. The one-way (positive $\hat{x}$-direction) sound wave $p_{\infty}(t)$ is generated via source terms in the governing equations (21) according to Bryngelson et al. [10]. Its form is a single period of a sinusoid with peak amplitude $0.3p_{0}$ and frequency $300\text{\,}\mathrm{kHz}$.

We start by considering a screen with fixed dynamic coordinate distributions $\sigma_{R}=\sigma_{\dot{R}}$, but varying distributions of equilibrium sizes $\sigma_{R_{o}}$. Polydispersity in $R_{o}$ is integrated via Simpson’s rule $61$ quadrature points for all cases. Figure 6 shows the bubble screen pressure for these cases as they evolve in time. For larger $\sigma_{R_{o}}$ (or broader distributions or bubble equilibrium sizes), the pressure is less oscillatory in time. A similar observation was made by Bryngelson et al. [10] for cases with no $R$ or $\dot{R}$ distributions. In the case of figure 6, we instead observe high-frequency oscillations in addition to the long-wavelength behaviors associated with the impinging pressure wave $p_{\infty}$. The origin of these oscillations is discussed next.

Figure 7 shows the dynamics associated with a bubble screen in varying degrees of statistical disequilibrium, represented via different $\sigma_{R}$ and $\sigma_{\dot{R}}$. Figure 7 (a) fixes $\sigma_{\dot{R}}$ and varies $\sigma_{R}$. We observe the shorter-wavelength oscillatory behavior, observed in figure 6, becoming more prominent for larger $\sigma_{R}$. These wavelengths are commensurate with the mean bubble natural frequencies, which superimpose the longer wavelength acoustics associated with the impinging $p_{\infty}$ wave. Figure 7 (b) shows a smoother pressure profile for larger $\sigma_{\dot{R}}$. Phase-cancellation between the larger waves associated with broader $\sigma_{\dot{R}}$ distributions and those of the $\sigma_{R}$ distributions may account for this behavior. Notably, these behaviors are qualitatively similar to those associated with varying $R_{o}$ distribution widths. Thus, parameterizing an $R_{o}$ distribution based on single-probe pressure measurements is insufficient.

We quantify the moment closure error, $\varepsilon_{c}$, via the mismatch in bubble screen pressure $p(t,x=0)$ due to truncated $R_{o}$ integration. Otherwise, the definition of $\varepsilon_{c}$ has the same form as $\varepsilon$ of (20), though the truth values are changed from Monte-Carlo simulations (which are unavailable) to high-resolution $N_{R_{o}}=401$ simulations as

This choice is made for two reasons. First, extending the QBMM method to additional $R$ and $\dot{R}$ quadrature points (or moments, equivalently) was found numerically unstable for most bubble cavitation problems. One approach to addressing this specific problem is introducing a recurrent neural network to correct the quadrature points and weights [48], though we do not discuss it further here. Second, we will see that the closure errors are most strongly associated with $N_{R_{o}}$, and thus focus on the influence of this parameter.

Figure 8 shows the $N_{R_{o}}$ closure errors associated with variations in all three PDF directions (panels a–c). For varying $\sigma_{R}$ (a) and $\sigma_{R_{o}}$ (c) we see that increasing variance $\sigma_{\cdot}$ results in larger closure errors $\varepsilon_{c}$. This appears to be associated with the larger pressure oscillations $p(x=0,t)$ observed for such cases. For figure 8 (b), we see the reverse trend, with larger $\sigma_{\dot{R}}$ corresponding to smaller closure errors, though this effect is small. Indeed, this effect matches that of the $R$- and $R_{o}$-direction effects, where larger $\sigma_{\dot{R}}$ results in smoother pressure histories. We also investigated the effect of Re, which results in smoother pressure profiles for smaller Re, and found the same trend.

Compelled by the closure error associated with $N_{R_{o}}$ dominating the total simulation error, we also implemented Gaussian–Hermite and Gauss–Legendre quadrature methods in the $R_{o}$-direction in an attempt to reduce the error. Figure 9 shows these closures errors $\varepsilon_{c}$ for increasing $N_{R_{o}}$. We see that these alternative quadrature rules do not significantly change the closure error, despite being optimal for a given $N_{R_{o}}$. This is because the integrand becomes highly oscillatory without viscous damping, as bubbles of different equilibrium radii $R_{o}$ become out of phase [40]. Long-time simulations without significant viscous effects, broadly polydisperse bubble populations (large $\sigma_{R_{o}}$) will require treatment of these developing oscillations. For example, a time-stationary analysis could be used to approximate the $R_{o}$-moments, so long as the bubble dynamics are fast compared to the fluid flow [40]. We will investigate this in future work.