## 1 Introduction

The dynamics of dispersions of small particles or bubbles in a fluid are important to many engineering and medical applications. In medicine, ultrasounds, generated via small cavitating bubbles, are employed during cataract removal [kelman1967phaco], to stop internal bleeding [vaezy1997liver, vaezy1999hemostasis], and in other procedures like tumor necrosis [bailey2003physical]. Focused shockwaves can cavitate bubbles that ablate kidney stones during lithotripsy treatment [coleman1987acoustic, pishchalnikov2003cavitation]. Their interaction with biological tissue or manufactured soft materials also attracts the medical [brennen2015cavitation, pan2018phenomenology, oguri2018cavitation, dollet2019bubble, spratt2021characterizing] and material science communities [gaudron2015bubble, barajas2017effects, turangan2017numerical]. Bubble cavitation is also responsible for damage and noise in hydraulic pipe systems [weyler1971investigation, streeter1983transient], hydro turbines [escaler2006detection, kumar2010study, luo2016review], and propellers [sharma1990cavitation, ji2012numerical]. At the same time, soots are critical to combustion [kazakov1998dynamic, balthasar2003stochastic, pedel2014large, mueller2009joint] and aerosols are used in many industrial processes [sibra2017simulation, laurent2001multi, hussain2015new]. In nature, cavitation is used as part of the hunting strategies of some marine animals, including humpback whales [leighton2004trapped, leighton2007acoustical, bryngelson2020simulation], mantis shrimps [patek2004deadly], and snapping shrimps [bauer2004remarkable, koukouvinis2017unveiling].

While the dynamics of these particles can be simulated directly for a specific (sampled) dispersion by tracking each particle, distribution statistics are typically sought in applications. In flows with large spatial gradients, a large ensemble of such simulations (Monte Carlo) is required to gather these statistics [zhao2007analysis, rosner2003multivariate]. The poor scaling of Monte Carlo makes such simulations expensive, and particle tracking also interferes with efficient parallelization. By instead phase-averaging the equations of motion [zhang1994ensemble], a two-way coupled set of Eulerian equations that are more suitable to parallelization and GPU processing is obtained. However, the averaged equations involve solving the generalized population balance equation (PBE) [ramkrishna2000population]. The PBE evolves the dispersed phase number density function (NDF) as a function of its dynamic variables [bryngelson2020qbmmlib]

. For example, the relevant variables for bubbles dynamics are the bubble radii and their radial velocities. However, further treatment is still required. The PBE is a partial differential equation in the dynamic variables, separate from the spatial and temporal variables of the flow equations, making this approach intractable for large simulations.

Quadrature-based moment methods (QBMMs) are a low-cost approach to approximately solving a PBE. Introduced in [mcgraw1997description], QBMMs have seen rapid improvement [marchisio2013computational]

. In brief, they prescribe a finite moment set and invert it to an optimal set of quadrature nodes and weights in the dynamic system phase space. The success of QBMMs has led to the creation of open-source libraries for them

[bryngelson2020qbmmlib, passalacqua2018open]. In the case of multiple dynamic variables, conditional QBMMs like conditional-QMOM (CQMOM) [yuan2011conditional] and conditional hyperbolic-MOM (CHyQMOM) [fox2018conditional, patel2019three] are preferred. These methods can efficiently solve many problems but suffer from a combinatorial explosion of their computational cost when higher accuracy is needed. This problem stems from the need to evolve all moments up to a higher order to increase accuracy. Worse still, these methods can exhibit numerical instabilities when third- or higher-order moments are evolved [marchisio2013computational].To circumvent these stability issues, this work employs neural networks to enhance the predictive abilities of standard -by--node (-node) CHyQMOM, which only requires access to first and second-order moments. This approach avoids both the numerical instabilities and high computational costs of evolving higher-order moments. The method follows the recent success of deep neural networks for improving multiphase flow models [ma2015using, ma2016using, charalampopoulos2021machine, bryngelson2020gaussian]. We expand on a previous such effort that used neural networks to close strictly-Gaussian moment transport equations [bryngelson2020gaussian]. Here, we instead seek data-informed corrections to a CHyQMOM method [fox2018conditional, patel2019three]. By doing this, one has control over the resulting quadrature nodes and weights. This makes correcting moment approximations straightforward and consolidates the two neural networks of [bryngelson2020gaussian] to one. This allows for computation of even out-of-training-set moments, in contrast to data-informed moment methods that use low-order moments to learn specific high-order ones [huang2021machinea, huang2021machineb, huang2021machinec]. Extension from [bryngelson2020gaussian] includes non-uniform and long-time pressure forcings, making the trained model appropriate for CFD solvers.

We emphasize that incorporating a trained neural network into a numerical method comes with its constraints. The neural networks train on a finite set of pressure profiles thought appropriate for a class of physical problems. However, one can not ensure that this fully generalizes, something which must be verified. Here, pressure profiles are sampled from a distribution that can represent many practical bubbly flow problems. Yet, the method is not guaranteed to generalize well for drastically different external forcings.

Section 2 formulates a model problem that serves as the basis for the proposed extension of CHyQMOM, namely the dynamics of a population of cavitating bubbles whose statistics are significantly non-Gaussian. In section 3, the new hybrid method is described. Section 4 shows that this approach can improve low-order moment predictions while extrapolating out of the moment space to compute required high-order moment predictions. This section also investigates the utility of additional quadrature points whose locations are selected by the RNN. It also compares the computational costs of the present approach and the classical CHyQMOM. Section 5 summarizes our conclusions.

## 2 Problem Formulation

### 2.1 Ensemble-averaged flow equations

This work focuses on the fluid-coupled dynamics of a dispersion of small, spherical bubbles transported in a compressible carrier fluid. The mixture phase-averaged evolution equations for the continuous fluid are

(1) |

where , and

being the mixture density, velocity vector, pressure, and total energy, respectively. The system of equations is complemented by appropriate initial and boundary or radiation conditions specific to each individual problem. The void fraction of the bubbles is

and a dilute assumption is made. The bubbles are defined by their instantaneous bubble radii , its time derivative . The bubbles are assumed monodisperse and so have the same equilibrium radius .The mixture pressure is deduced from the ensemble phase-averaging method [zhang1994ensemble, bryngelson2019quantitative] as

(2) |

where and are the bubble wall and liquid pressures, respectively [ando2011numerical]. Liquid pressure follows from the stiffened-gas equation of state [menikoff1989riemann], though this model can be substituted for another if required. The usual coefficients for water are used [bryngelson2019quantitative].

The overbars in 2 denote raw moments of the bubble dispersion as

(3) |

where is the number density function of the bubbles. This paper focused on a new, improved method for computing these moments, which will be introduced in section 3.

The void fraction transports as [bryngelson2019quantitative]:

(4) |

The moments required to close the governing flow equations are thus

(5) |

### 2.2 Bubble model

To close the governing equations of the previous subsection, a model for the bubble dynamics, in terms of the dynamical variables and , is required. We use a Rayleigh–Plesset equation for this:

(6) |

which is dimensionless via the reference bubble size , liquid pressure , and liquid density . In (6), is the ratio between the fluid and bubble pressures and is a Reynolds number

(7) |

where is the liquid kinematic viscosity. For the cases considered here .

This model assumes the bubbles remain spherical and compress via a polytropic process with coefficient . While this model can be generalized to include heat exchange and liquid compressibility, these effects are not critical to our study and thus omitted here. Based on this model, the bubble wall pressure simplifies the last moment of as

(8) |

We also define a dimensionless time , where is the natural frequency of the bubbles. To simplify the notation, will be used in place of hereon.

### 2.3 Population balance formulation

A number density function describes the statistics of the bubbles. The generalized population balance equation is

(9) |

assuming the bubbles do not coalesce or break up, though these effects can be included via empirically modeled terms if desired. QBMMs solve (9) by representing as a set of raw moments [fox2018conditional, bryngelson2020qbmmlib]. Through an appropriate inversion procedure, these methods can transform these moments into quadrature nodes and weights in phase space. This allows for the approximation of via a weighted sum of Dirac delta functions. Hence, the following quadrature rule can approximate any raw moment

(10) |

where is the -th quadrature point locations for the -th internal coordinate.

## 3 Hybrid quadrature moment method formulation

We now present the hybrid, data-informed method for predicting the moments of cavitating bubbles. Subsection 3.1 presents -node CHyQMOM predictions [bryngelson2020qbmmlib]. Subsection 3.2 details the hybrid neural-network model that improves the predictions.

### 3.1 CHyQMOM

For two-dynamic-variable cases, conditioned moment methods are computationally preferable to traditional QMOM [yuan2011conditional]. We use CHyQMOM because it can close a second-order moment system with fewer carried moments than CQMOM [fox2018conditional]. For a -node quadrature rule, it uses the first- and second-order moments

(11) |

The CHyQMOM inversion process for obtaining the nodes and weights is presented in appendix A.

Taking the time-derivative of each of the (11) moments and applying (6) results in

(12) |

which are called the moment transport equations. The quadrature rule (10) approximates unclosed moments in (12).

While this scheme is computationally cheap, it is challenging to extend to include additional quadrature points without high computational cost or numerical instabilities. Thus, truncation errors can affect approximation of the right-hand-side of (12) and the extrapolation out of the low-order moment space to the moments of of (5). We next present an augmented method that treats these issues without introducing numerical instability or high computational cost.

### 3.2 Data-informed corrections

We improve the CHyQMOM moment inversion procedure by adding a correction term to the -node quadrature rule and introducing additional quadrature nodes. The unaugmented CHyQMOM quadrature rule is denoted via .

For these corrections, an LSTM RNN is employed. The LSTM incorporates non-Markovian memory effects into the reduced-order model. This approach is known to be capable of improving predictions of reduced-order models [vlachas2018data, bryngelson2020gaussian].

The corrections serve as input predictions for the first- and second-order moments as well as the pressure . They are modelled as

(13) |

where the vector

denotes hyperparameters and optimized parameters of the neural network as obtained during training. More detail on the implementation is in section

4, subsection 4.1. The chosen values for the hyperparameters are included in appendix B. The history of the reduced-order model states is(14) |

The hybrid quadrature rule follows as

(15) |

The neural network loss function is designed to ensure the target high-order moments

can be accurately computed and that the low-order moments evolve accurately. Hence, the right-hand-side of (12) is included in the loss function as(16) |

where,

(17) |

and

(18) |

The first term in (16) minimizes, in the sense, the right-hand-sides of (12) (given ). The second term in (16) minimizes prediction error for both and , while it also penalizes the network when the weights don’t sum up to (under the assumption that ). The last term in (16) penalizes negative-valued weights.

The discretized moment transport equation (12) and the quadrature rule 10 compute the time-derivatives required in (16) as

(19) |

Once trained, the scheme

results in a new quadrature rule that evaluates the right-hand-side of (12). The moment transport equations (12) then evolve via an adaptive th-order Runge–Kutta time-stepper. Algorithm 1 describes this procedure.

Note that the closure terms need to be evaluated at times , , and . The neural network does not make predictions at , so the equations are instead integrated in time by instead of .

## 4 Results

### 4.1 Pressure signals

The capabilities of the data-enhanced CHyQMOM method to predict the statistics of bubble populations are explored. The pressure term excites the bubbles causing oscillations. The representation of used here should be general enough to include pressure profiles seen in actual fluid flows. In a generic framework, let have a finite Fourier-series expansion

(20) |

where corresponds to nondimensional time, nondimenionalized by the natural oscillation frequency of the bubbles , are the included dimensionless frequencies, and and are the corresponding amplitude and phase. It is stressed that corresponds to the equilibrium pressure of the bubbles (for which and ).

Most cavitating flow applications do not contain pressure frequencies higher than the natural oscillation frequency of the bubbles [brennen2014cavitation] (with a notable exception of some high-intensity focused ultrasound treatments). We operate under this constraint hereon, though higher frequencies could be included if desired. On the other hand, very low frequencies are uninteresting because they cause the bubbles to evolve quasi-statically. Hence, without significant loss of generality, the dimensionless frequencies of are in the interval . The phases of the waveforms that make up

are independently sampled from a uniform distribution

in(21) |

Applications dictate the possible observed pressure amplitudes. For example, significantly low pressures are not relevant for most applications. To set an empirical threshold approximating this condition, the pressures must not cause the used Monte Carlo simulation configuration to become numerically unstable. The solver itself uses an adaptive 3rd-order Runge–Kutta scheme with minimum time-step and relative error tolerance of . Thus, we design a pressure distribution from which all samples are numerically stable and physically realistic. Algorithm 2 details this process.

Previous experimental works can also be used to justify that the forcing constraint in algorithm 2

avoids abrupt cavitation. This is estimated by the cavitation number

(22) |

where is the vapor pressure of the liquid at reference temperature and is the reference velocity [brennen2014cavitation]. If the liquid cannot withstand negative pressures then vapor bubbles appear when the liquid pressure is . Thus, nucleation occurs when

(23) |

Without loss of generality, we can choose to simplify the following calculations. For flows around axisymmetric headforms, with Reynolds number in the range of , if , the formed nuclei grow explosively up to a certain bubble size [ceccio1991observations]. This phenomenon renders numerical simulations for flows with considerably more expensive compared to cases in order to achieve the same numerical accuracy. For the pressure profiles presented here, the case is avoided when using algorithm 2.

Figure 1 shows example pressure profiles that are used to test the fidelity of the hybrid moment inversion method. Herein, the end of this time window, , is used to assess model fidelity. This enables the bubble dynamics to evolve from a specific initial state to one more representative of those found in actual flows.

### 4.2 LSTM RNN training procedure

We simulate samples of individual bubbles for each realization of

. The bubbles are initialized via samples from normal distributions with variances

and for and respectively. Each case is evolved until , which in this dimensionless system corresponds to natural periods of bubble oscillations. The individual bubble dynamics are then averaged to obtain the Monte Carlo reference statistics for each realization.For the numerical investigation, samples of from (20) are used. From those, are randomly selected for training, with the remaining cases used during testing. The number of samples used for training is chosen so that it is large enough to avoid over-fitting but small enough to still allow for the sampling of qualitatively different pressure profiles during testing.

Figure 2 shows and the for one realization of at different time instances. The same figure displays the CHyQMOM nodes as estimated by both the standard -node CHyQMOM scheme and the - and -node hybrid schemes.

During traning, the LSTM memory size is set to time-instances, with each apart. The Adam method [kingma2014adam] trained each neural network for epochs, minimizing the loss function (16). A table with the values of the hyperparameters of the neural networks is presented in appendix B. A 4th-order Runge–Kutta adaptive time stepper evolves the predictions of the hybrid scheme. The time integration scheme is adaptive, but the neural network predictions are uniform, so the neural network corrections are limited to the associated fixed time step .

### 4.3 Low-order moment evolution and error quantification

The model-form relative error is computed via a discrete error

(24) |

where (MC) indicates Monte Carlo surrogate truth data and . The are uniformly spaced times in the interval . Results regarding the low-order moments are presented in figure 3.

Figure 3 (i) shows for the first- and second-order moments for the -node schemes. Rows (ii–iv) shows the evolution of select moments for and row (v) shows the corresponding . All results correspond to randomly selected testing samples (a)–(d) as labeled. We observe a smaller for the hybrid scheme than standard CHyQMOM for all cases considered. The largest errors for both approaches appear for moment (row (iv)), which exhibits large and intermittent fluctuations when the bubbles collapse. The CHyQMOM method deviates most from the MC surrogate-truth data during intervals of high compression (small ), with hybrid CHyQMOM tracking the reference solution more accurately. Thus, the observed increase in accuracy varies significantly from case-to-case and moment-to-moment, from times smaller errors to only improvements.

### 4.4 High-order moment extrapolation

The next quantity of interest is the -error in predicting the target higher-order moments . Figure 4 presents these results for the same testing samples presented in figure 3. For all moments (ii)–(iv), hybrid CHyQMOM significantly improves the predictions of . This improvement is associated with the more accurate evolution of the low-order moments and these targeted moments in the neural network training procedure.

To better study the typical reduction in error for the -node hybrid CHyQMOM scheme, we compute the percent improvement of the error as

(25) |

Figure 5 shows a histogram of calculated for example low-order moments and target moments . For all testing samples, the -node hybrid scheme improves the accuracy of the standard CHyQMOM method.

Further, for both low-order moments () and target moments (), the error is reduced by more than for more than of the sampled cases. The variation in error improvement is due to the amplitude range of the sampled and how closely the time-series of corresponds to one of the training samples. Thus, results can improve by including more training samples.

### 4.5 Additional quadrature nodes

Another potential route to method improvement is to increase the number of quadrature nodes. While the number of quadrature nodes can change, the evolved moments remain . The algorithm for this is included in appendix A.

To quantify the effect of this change, is computed, which is the median error among the test samples. We then define

(26) |

to quantify the decrease in the -error when using higher-node-count hybrid CHyQMOM compared to the standard -node CHyQMOM.

Figure 6 shows that the accuracy of the hybrid predictions is improved as the number of nodes increases.
However, the error improvement gains diminish once nodes are reached.
Further, including additional nodes to the quadrature rule increases the computational time needed to perform a single time-step evolution for the system.
The computational cost of -node hybrid CHyQMOM per time-step is times the cost of CHyQMOM.
For , , and nodes, the hybrid method costs per time-step are , , and times that of -node CHyQMOM^{1}^{1}1
These simulations were performed using [bryngelson2020qbmmlib] on a single core of a Intel Core i CPU..
Hence, diminishing improvements are observed as the number of nodes increases to more than , as the simulations require significantly more computations per time-step for comparable accuracy.

## 5 Conclusions

A data-informed conditional hyperbolic quadrature method for statistical moments was presented. The method was applied to the statistics of a population of spherical bubbles oscillating in response to time-varying randomized forcing. The forcing is designed to resemble any possible function with a banded frequency spectrum from to the natural frequency of the bubbles. Results showed that the method reduces closure errors when compared against a standard -node CHyQMOM scheme. The hybrid method reduced errors more significantly for the extrapolated higher-order moments that close the phase-averaged flow equations. This improvement was achieved via recurrent neural networks that include time history during training. This result is significant because higher-order QBMM schemes are generally numerically unstable for this problem, so another approach is required to improve accuracy. Thus, while the presented hybrid scheme is about a factor of more expensive than CHyQMOM, its numerical cost should be viable for many applications.

All authors jointly conceived the problem formulation. A.C and S.B. designed the study. A.C. developed the ML code, conducted the numerical experiments, and drafted the manuscript. S.B. contributed the bubble dynamics code and further developments of . All authors reviewed and edited the manuscript.

We declare we have no competing interests.

The US Office of Naval Research supported this work under grant numbers N0014-17-1-2676 and N0014-18-1-2625. Computations were performed via the Extreme Science and Engineering Discovery Environment (XSEDE) under allocation CTS120005, supported by National Science Foundation grant number ACI-1548562.

The authors thank Rodney Fox and Alberto Passalacqua for valuable discussion of quadrature moment methods.

## References

## Appendix A CHyQMOM inversion algorithm

This is the inversion algorithm for the -node CHyQMOM scheme. Given the first and second-order moments it computes the nodes and corresponding weights for , in phase-space. In this work we assume . To tail the algorithm to our hybrid scheme of arbitrary number of quadrature nodes, the algorithm adds some fictitious extra nodes to the scheme with zero-valued weights to match the desired number of nodes of the hybrid scheme.

## Appendix B Neural network hyperparameters

Hyperparameter | Value |
---|---|

Epochs | |

Learning rate | |

Batch size | |

Activation function | tanh |

Recurrent activation function | hard sigmoid |

Dropout coefficient | |

Recurrent dropout coeffificent | |

LSTM is stateful | True |

Kernel initializer | Zeros |

Recurrent initializer | Zeros |

Bias initializer | Zeros |