Regressing bubble cluster dynamics as a disordered many-body system

Abstract

The coherent dynamics of bubble clusters are of fundamental and industrial importance, and are elusive due to the complex interactions of disordered bubble oscillations. Here we introduce and demonstrate a method for decomposition of the Lagrangian time series of bubble dynamics data by combining theory and principal component analysis. The decomposition extracts coherent features of bubble oscillations based on their energy, in a way similar to proper orthogonal decomposition of Eulerian flow field data. This method is applied to a dataset of spherical clusters under harmonic excitation at different amplitudes, with various nuclei density and polydispersity parameters. Results indicate that the underlying correlated mode of oscillations is isolated in a single dominant feature in cavitating regimes, independent of the nuclei's parameters. A systematic data analysis procedure further suggests that this feature is globally controlled by the dynamic cloud interaction parameter of Maeda and Colonius (J. Fluid Mech., vol. 862, 2019, pp. 1105–1134) that quantifies the mean-field interactions, regardless of initial polydispersity or nonlinearity. The method provides a simplified and comprehensive representation of complex bubble dynamics as well as a new path to reduced-order modelling of cavitation and nucleation.

1. Introduction

Cavitation bubble clusters nucleate when the liquid pressure rapidly falls below a certain threshold. These clusters coherently oscillate and violently collapse to cause extreme energy concentration that leads to various critical consequences, and use in applications as diverse as injectors and pumps (e.g. Plesset & Ellis 1955; Blander & Katz 1975; Mørch 1980; Chen & Heister 1994; Hashimoto et al. 1997; Prosperetti 2017), hydraulic machines (Arndt 1981), underwater propulsion and hydrofoils (e.g. Kubota, Kato & Yamaguchi 1992; Merkle, Feng & Buelow 1998; Kunz et al. 2000; Schnerr & Sauer 2001; Ganesh, Mäkiharju & Ceccio 2016; Gnanaskandan & Mahesh 2016; Venning, Pearce & Brandner 2022), medical ultrasound (e.g. Ikeda et al. 2006; Krimmel et al. 2010; Maxwell et al. 2011; Pishchalnikov, Williams & McAteer 2011; Maeda et al. 2015; Movahed et al. 2016; Maeda 2020), surface cleaning (Verhaagen & Rivas 2016; Yamashita & Ando 2019), chemical synthesis (Suslick et al. 1999; Cairós & Mettin 2017), and bio-inspired devices (Tang & Staack 2019). Characterizing the dynamics is a challenge due to the complex interactions of bubbles involving disorders and stochasticity. Nuclei are typically micro-sized and polydisperse and randomly distributed. Their rapid, nonlinear oscillations are, except in controlled experiments (Bremond et al. 2006), practically not measurable. Molecular and hydrodynamics simulations can provide detailed insights into nucleation (Angélil et al. 2014; Gallo, Magaletti & Casciola 2021), while their time- and spatial scales have not reached those of practical cluster oscillations. Analyses have been made on the interaction dynamics in various regimes (Brennen 2014), yet no common knowledge has been established if the many-body coherence globally exists and if so scaling is possible, beyond the consensus that polydispersity induces strong disorders.

For the past decades, the Rayleigh–Plesset (R–P) equation and its variations have been actively explored to investigate the dynamics of single bubbles (Plesset 1949; Plesset & Prosperetti 1977). Relatively few studies addressed the theory of clusters. By using the mean field approach to interacting bubbles modelled by the R–P equation, d'Agostino & Brennen (1989) derived a non-dimensional parameter that dictates the linear coherent oscillations of monodisperse clusters, the so-called ‘cloud interaction parameter’. Zeravcic, Lohse & Van Saarloos (2011) used the coupled R–P equations and identified disorders represented by the Anderson localization of acoustic energy in polydisperse, lattice-like clouds under weak excitation. We have recently extended the interaction parameter to the non-equilibrium, cavitating clusters under strong excitation by considering the effective interaction at excited states (Maeda et al. 2018; Maeda & Colonius 2019; Maeda & Maxwell 2021). To recall, we scale the mean kinetic energy of liquid induced by $N ({\gg }1)$ bubbles as (Maeda & Colonius 2019)

(1.1)\begin{equation} \langle K\rangle \sim \langle K_s \rangle (1+B_d),\end{equation}

where $K_s$ is the energy of a single bubble: $K_s=2{\rm \pi} \rho R^3_{b,c}\dot {R}^2_{b,c}$; $B_d$ is the parameter controlling the effective contribution of hydrodynamic inter-bubble interaction: $B_d = N{\langle \overline {R(t)} \rangle }/R_C$, and $R_{b,c}$ and $R_C$ denote the characteristic (reference) bubble radius and the cluster radius; $\rho _l$ is the liquid density; $\langle {\cdot } \rangle$ and $\overline {({\cdot })}$ denote time average during a period in which bubble dynamics are statistically stationary and the mean value about the bubbles ($i=1,2,\ldots,N$), respectively. The scaling can be simply derived from the coupled R–P equation for the correlated (synchronized) limit of monodisperse bubbles ($R_1=R_2\cdots =R_{b,c}$). Although realistic correlations are imperfect due to polydispersity and nonlienarity, $B_d$ was found to control well both the coherent dynamics of polydisperse cavitating clusters, and their acoustic emission in numerical simulations and experiments (Maeda & Colonius 2019; Maeda & Maxwell 2021). Overall, previous studies indicate that the coherence can depend on both polydispersity and nonliearity in a non-separable manner, posing perplexing questions about the universality of scaling. The theoretical characterization of the nonlinear dynamics of disordered many-body systems is, in general, not a simple task. Meanwhile, greater computing power has enabled learning physics by analysing big data. Principal component analysis (PCA) is a powerful method for unsupervised learning which has seen recent success in characterizing the coherent physics of many-body and high-dimensional systems in fields ranging from quantum information to fluid dynamics (Milano & Koumoutsakos 2002; Holmes et al. 2012; Lloyd, Mohseni & Rebentrost 2014; Taira et al. 2017).

In this study, we introduce and demonstrate a method for unsupervised data decomposition to study the coherent bubble cluster dynamics by combining theory and PCA. PCA extracts dominant states and dynamical features, such as coherent quantum states and turbulent structures, and their amplitudes as the eigenfunction (feature) and the eigenvalue (variance) of the co-variance matrix of physical data. When applied to spatio-temporal data of dynamical systems, PCA is often denoted as proper orthogonal decomposition (POD). Those data can be properly weighted prior to PCA such that the variance becomes consistent with the norm induced by an energetic inner product of state variables (e.g. kinetic energy) (Lall, Marsden & Glavaški 1999; Rowley 2005). This weighting allows a physical interpretation that resulting features associated with a large variance are energetically dominant coherent structures. Proper weighting of Lagrangain bubble dynamics data is non-trivial since the linear variance of extracted features needs to account for the nonlinear interaction energy. For meaningful analysis, we introduce strategic pre-processing of the data prior to PCA such that the PC-variance becomes theoretically consistent with the energy modelled by the coupled R–P equation. Analysing simulation datasets of clusters, we show that the PCA can systematically extract not only coherent but also incoherent features whose magnitudes are respectively measured by the PC-variance and the entropy. We discover that the coherence is lost by disorders induced by polydispersity and nonlinearity, while under strong excitation, the underlying correlations are globally isolated in a single coherent feature whose variance (energy) is scaled by $B_d$, regardless of the disorders.

The remainder of this paper is as follows. In § 2, we describe the method. The first PC variance, spectral entropy and the coherence measure are introduced as quantifiable measures to characterize the coherent bubble dynamics from extracted features. In § 3, we verify and demonstrate the method using a numerical dataset of bubble clouds with different density and polydispersity parameters under various amplitudes of harmonic excitation. The deviation of the PC-variance from the physical energy is quantified for two weighting methods. The amplitude dependencies of the measures are quantified. The PC-spectra and their correlations with the coherent energy are analysed. Moreover, the extracted coherent dynamics is related to $B_d$ and its universality is discussed for cavitating clouds. In § 4, the physical significance of the method is discussed. In § 5, we state conclusions.

2. Methods

2.1. Principal component analysis of bubble dynamics data

For clusters modelled by the coupled R–P equations, observable dynamical variables are the bubbles’ radial velocities and radii. Consider a data matrix $\boldsymbol {Q}$ containing the $N_t$ snapshots of the radial velocities with a constant temporal interval,

(2.1)\begin{equation} \boldsymbol{Q}=[\boldsymbol{q}_1,\boldsymbol{q}_2,\dots,\boldsymbol{q}_{Nt}], \end{equation}

where $\boldsymbol {q}_k$ denotes the vector containing the radial velocities of the $N$ bubbles at time $t_k$ (in the $k$th snapshot),

(2.2)\begin{equation} \boldsymbol{q}_k=[\dot{R}_1(t_k),\dot{R}_2(t_k),\dots,\dot{R}_{N}(t_k)]^{\rm T}. \end{equation}

For later convenience, we also define the vector $\boldsymbol {r}_k$, containing the radii of the bubbles at the same instances:

(2.3)\begin{equation} \boldsymbol{r}_k=[R_1(t_k),R_2(t_k),\dots,R_N(t_k)]^{\rm T}. \end{equation}

PCA can be performed on $\boldsymbol {Q}$ by using the singular value decomposition (SVD) (e.g. Abdi & Williams 2010; Jolliffe & Cadima 2016):

(2.4)\begin{equation} \boldsymbol{Q}=\boldsymbol{U}\boldsymbol{\varSigma}\boldsymbol{X}^*. \end{equation}

The $i$th principal component (feature) is stored in the $i$th score matrix,

(2.5)\begin{equation} \boldsymbol{\varPi}_{i}=\boldsymbol{U}\boldsymbol{\varSigma}_{i}, \end{equation}

where $\boldsymbol {\varSigma }_{i}$ contains the $i$th largest singular value and zeros elsewhere.

2.2. Weighted principal component analysis

Although this procedure for PCA is simple and straightforward, the physical meaning of the extracted features are obscure since PCA itself is not informed on the underlying dynamics of the system. A related issue of PCA for fluid flow data has been addressed in the context of POD. In POD of the snapshots of Eulerian fluid flow data, the state vectors consisting of the data matrix are often weighted such that the corresponding weighted inner product of the state becomes consistent with the kinetic energy of the original system. The dominant features (POD modes) can then be interpreted as energetically dominant coherent flow structures. Moreover, in our previous studies (Maeda & Colonius 2019; Maeda & Maxwell 2021), the scaling of bubble cloud dynamics by $B$ was successfully demonstrated based on the total kinetic energy of liquid induced by interacting bubbles. These considerations motivate us to relate the energy to PCA for addressing the physics of bubble cloud.

Inspired by the POD, we consider weighting the present bubble dynamics data prior to PCA. To find an appropriate weight, we revisit the potential theory behind (1.1). The kinetic energy of the fluid induced by the spherical bubble oscillations is explicitly expressed as

(2.6) \begin{equation} K = 2{\rm \pi}\rho_l\sum^N_{i=1} \Big[R_i^3\dot{R}_i^2+\sum^N_{j\ne i}\frac{R_i^2R_j^2 \dot{R}_i \dot{R}_j}{r_{ij}}\Big]+(H.O.T),\end{equation}

where $R_i$, $\dot {R}_i$ and $r_{ij}$ are the radius and the radial velocity of bubble $i$, and distance between the centres of bubble $i$ and $j$, respectively. The second term in the brackets represents the contribution of the long-range interactions. At each instant, $K$ can be expressed as a weighted inner product of $\boldsymbol {q}$:

(2.7)\begin{equation} K=\boldsymbol{q}^T\boldsymbol{T}\boldsymbol{q}=(\boldsymbol{W}\boldsymbol{q})^{\rm T}(\boldsymbol{W}\boldsymbol{q}), \end{equation}

where

(2.8)\begin{equation} T_{ij}(\boldsymbol{r}_k) = \begin{cases} 2{\rm \pi}\rho R^3_i(t_k) & (i=j),\\[2pt] 2{\rm \pi}\rho \dfrac{R^2_i(t_k)R^2_j(t_k)}{r_{ij}} & (i\ne j). \end{cases} \end{equation}

The weight matrix, $\boldsymbol {W}$, can be obtained through the Cholesky factorization of $\boldsymbol {T}(\boldsymbol {r})$:

(2.9)\begin{equation} \boldsymbol{T}(\boldsymbol{r})=\boldsymbol{W}\boldsymbol{W}^*. \end{equation}

This expression suggests the weighted data, $\boldsymbol {W}\boldsymbol {Q}$, is appropriate for use in PCA. However, $\boldsymbol {T}(\boldsymbol {r})$ and $\boldsymbol {W}$ are time dependent since ${\boldsymbol {r}}$ can change in time during large-amplitude oscillations. The choice of $\boldsymbol {r}_k$ in the time series to define $\boldsymbol {W}$ is unclear. A straightforward choice is to use the temporal mean of the radius for each bubble, $\langle \boldsymbol {r}\rangle : \langle \boldsymbol {r}\rangle =\sum ^{Nt}_{k=1}{\boldsymbol {r}_k}/Nt$, but $T(\langle \boldsymbol {r}\rangle )$ is not a favourable approximation for $T(\boldsymbol {r})$ for bubbles with large-amplitude oscillations.

We address this obstacle by variable transformation as a means of pre-processing the data. The schematic is shown in figure 1. First, we transform the variables from $(\boldsymbol {q}, \boldsymbol {r})$ to $(\boldsymbol {\xi }, \boldsymbol {\eta })$, where $\boldsymbol {\xi }=[\xi _1, \xi _2,\ldots,\xi _N]^{\rm T}$ and $\boldsymbol {\eta }=[\eta _1, \eta _2,\ldots,\eta _N]^{\rm T}$, and $(\xi _i,\eta _i)=(q_ir_i, r_i)=(\dot {R}_iR_i,R_i)$ for $i\in [1,N]$. Here, $\xi _i$ is nothing but the velocity potential evaluated at the surface of bubble $i$. Supplemental discussions and justifications for this transformation are provided in Appendix A. The instantaneous energy of the system is expressed by a weighted inner product as

(2.10)\begin{equation} K=\boldsymbol{\xi}^T\boldsymbol{P}(\boldsymbol{\eta})\boldsymbol{\xi}=(\boldsymbol{\varOmega}\boldsymbol{\eta})^{\rm T} (\boldsymbol{\varOmega}\boldsymbol{\eta}), \end{equation}

and $\boldsymbol {\varOmega }$ is the weight matrix satisfying

(2.11)\begin{equation} \boldsymbol{P}(\boldsymbol{r})=\boldsymbol{\varOmega}\boldsymbol{\varOmega}^*, \end{equation}

where

(2.12)\begin{equation} P_{ij}(\eta) = \begin{cases} 2{\rm \pi}\rho \eta_i, & i=j,\\[2pt] 2{\rm \pi}\rho \dfrac{\eta_i \eta_j}{r_{ij}}, & i\ne j. \end{cases} \end{equation}

Second, we partially replace $\boldsymbol {\eta }$ with $\langle \boldsymbol {\eta }\rangle$ to approximate the system. Using $(\hat {\boldsymbol {\xi }},\hat {\boldsymbol {\eta }})$, the instantaneous energy of the approximate system is expressed as

(2.13)\begin{equation} \hat{K}=\hat{\boldsymbol{\xi}}^T{\boldsymbol{P}}(\hat{\boldsymbol{\eta}})\hat{\boldsymbol{\xi}}.\end{equation}

The temporal mean of the energy of the original system can then be approximated as

(2.14)\begin{equation} \langle K\rangle \approx\langle\hat{K}\rangle =\langle\hat{\boldsymbol{\xi}}^T{\boldsymbol{P}(\hat{\boldsymbol{\eta}})}\hat{\boldsymbol{\xi}}\rangle \approx\langle\boldsymbol{\xi}^T\boldsymbol{P}(\langle\boldsymbol{\eta} \rangle)\boldsymbol{\xi}\rangle.\end{equation}

This is a critical result in the present context of PCA, since the energy is now related to the weighted inner product with the constant weight, $\boldsymbol {P}(\langle \boldsymbol {\eta } \rangle )$. Quantitative verification of this approximation is addressed through numerical experiments in the following section. Using the new set of variables, the PC decomposition is performed as

(2.15)\begin{equation} \boldsymbol{\varOmega}\boldsymbol{\varXi}=\boldsymbol{U}_{\xi}\boldsymbol{\varSigma}_{\xi}\boldsymbol{X}_{\xi}^*.\end{equation}

The $i$th PC is stored in $\boldsymbol {\varPi }_{\xi,i}$: $\boldsymbol {\varPi }_{\xi,i}=\boldsymbol {U}_{\xi }\boldsymbol {\varSigma }_{\xi,i}$. We denote the $i$th largest singular value of $\boldsymbol {\varSigma }_{\xi }$ as $\sigma _i$. The degree of coherence can be measured by the normalized variance of the first PC:

(2.16)\begin{equation} \hat{\sigma}_1^2=\sigma_1^2/\mathrm{tr}(\boldsymbol{\varSigma}^2), \end{equation}

which represents the ratio of the energy occupied by the first PC to the total energy of the system. In the following sections, to distinguish the present method from the standard PCA, we denote (2.15) as the weighted principal component analysis of transformed data (WPCA-TD).

Figure 1. Schematic of the feature extraction by PCA from Lagrangian bubble dynamics data, $\boldsymbol {Q}$, after pre-processing. The variance of the resulting features is consistent with the interaction energy predicted by the coupled Rayleigh–Plesset equation.

2.3. Spectral entropy and coherence measure

We introduce two additional key measures to characterize the coherent physics of bubble clusters through WPCA-TD. First, to quantify the degree of incoherent bubble oscillations, we define the spectral entropy of the weighted co-variance matrix, $\boldsymbol {\varOmega }\boldsymbol {\varXi }$, as

(2.17)\begin{equation} \hat{S}_{vN}=-\sum\frac{\hat{\sigma}^2_k(\ln\hat{\sigma}^2_k)}{\ln(N)}, \end{equation}

where $\ln _2(N)$ is the normalization factor. The spectral entropy of data is a discrete analogue of the Shannon entropy based on the spectrum of data and has been used to quantify the randomness of data in various applications (Kullback 1997; Alter, Brown & Botstein 2000; Hu et al. 2005; De Domenico & Biamonte 2016). Aubry, Guyonnet & Lima (1991) introduced a similar definition of entropy for spatio-temporal signals of canonical fluid flows from their POD eigenvalues and used the entropy to characterize flow instabilities. In the present context, when bubbles are in perfect correlation, we expect to excite only the first PC capturing the entire energy of the system: $\hat {\sigma }^2_1=1$ and $\hat {\sigma }^2_k=0$ for $k: k\in [2,N]$, yielding $\hat {S}_{vN}=0$. In contrast, if the energy is equi-partitioned into all PCs, $\hat {\sigma }^2_k=1/N$ for all $k: k\in [1,N]$ and $\hat {S}_{vN}=1$. Second, we define the coherence measure $C$:

(2.18)\begin{equation} C=\frac{B_1}{B_d}, \end{equation}

where

(2.19)\begin{equation} B_1=\frac{\sigma^2_1}{\sigma'^2_1}-1. \end{equation}

Here, $\sigma '^2_1$ is the first PC-variance excluding the contribution of interactions, obtained from the WPCA-TD of the same data using a diagonal weight matrix $\boldsymbol {\varOmega }'$ whose diagonal entries are those of $\boldsymbol {\varOmega }$. Further details of $\boldsymbol {\varOmega }'$ are provided in Appendix B. In the limit of perfect correlation (also see (1.1)),

(2.20)\begin{equation} \langle K\rangle \approx \sigma^2_1 \approx (1+B_d)\sigma'^2_1\end{equation}

and

(2.21)\begin{equation} B_d\approx B_1. \end{equation}

Therefore, $C\approx 1$. This condition is typically realized in monodisperse clusters under weak (linear) oscillations (d'Agostino & Brennen 1989). For real clusters under strong excitation, bubbles are not perfectly correlated and the energy can be distributed in broad features. In this regime, approximation (2.20) is not necessarily expected to hold and $C$ can take values far from unity. Phenomenologicaly speaking, $C$ quantifies from data the degree to which the most coherent mode of oscillation is represented by the cluster's mean-field interaction. A schematic of the process to obtain $C$ from data is provided in Appendix B.

3. Numerical experiments

3.1. Datasets

To verify and demonstrate the WPCA-TD, we use datasets of spherical clusters excited by 40 cycles of harmonic pressure excitation. Each cluster contains $O(10-10^3)$ bubbles with their initial radii following log-normal distributions with a reference radius of $R_{ref}=O(10) \mathrm {\mu }{\rm m}$; $\ln (R_0/R_{ref})=N(0,s_d)$, where $s_d$ is the lognormal standard deviation as the measure of polydispersity (Maeda & Colonius 2019). We address $s_d=[0.1,0.3,0.5,0.7]$. In real bubble clouds, nuclei are expected to be polydisperse. Here, $s_d= 0.7$ may be a representative estimation based on previous studies (Katz 1978; Ando, Liu & Ohl 2012; Maeda & Colonius 2019) and is used unless noted. The smaller values of $s_d$ are used in some cases to quantify the effect of polydispersity. The bubbles are randomly distributed in the spherical region with a specified cluster radius. Similar parameters of clusters were previously simulated to compare with experiments (Maeda & Colonius 2019; Maeda & Maxwell 2021). The density of bubbles is characterized by $B_0$, the value of $B_d$ at rest. The far-field pressure is given as $p_{\infty }(t)=p_0[1+A\sin (2{\rm \pi} ft)]$. The frequency of excitation is $f=\ $500 kHz unless noted, near the adiabatic resonant frequency of the reference bubble. The amplitude of excitation is defined by $A$, relative to the ambient pressure at $p_0=1.0$ atm. For each set of parameters, we compute an ensemble average by taking a mean of the results from 20 bubble clouds with distinct spatial placements of bubbles in the clouds. For data generation, we use mesh-free, coupled Keller–Miksis equations by modifying previous methods. Details of this method are provided in Appendix D.

3.2. Visualization of representative data

In figure 2, we show the evolution of representative quantities of a bubble cloud in the dataset with $(B_0,s_d,N)=(0.5,0.1,100)$ excited at $A=20$, during a stationary state. Figure 2(a) shows the radius of a representative bubble. The plot presents familiar features of cavitation bubbles forced by continuous, strong excitation, including fast events of collapse/rebound and slow growth/decay between them. During the slow phase, small-amplitude oscillations at excitation amplitude are evident. For the same bubble, figure 2(b) shows the evolution of the velocity potential obtained from the raw data, and those from the first and second dominant features extracted through the WPCA-TD. Compared with the radius, the raw data of the velocity potential looks much more symmetric about zero, even around the collapse events. The potential of the first feature presents a sinusoidal profile at the excitation frequency and with a constant amplitude which is close to the peak amplitude of the original data. The collapse events are not captured in this feature. The potential of the second feature is symmetric but has a much lower amplitude compared with the first feature. There is no clear similarity between the original data and the second feature, unlike that between the original data and the first-PC. Figure 2(c) shows the void fraction. The fraction oscillates with an amplitude of approximately $5.0\times 10^{-4}$ near $3.0\times 10^{-3}$ at the excitation frequency. Slow, small amplitude of variations are also observed. Figure 2(d) shows the kinetic energy of fluid induced by bubble oscillations. The energy oscillates approximately $0.75\ \mathrm {\mu }{\rm J}$ at doubled the excitation frequency. Although the frequency is as expected, the peak amplitude largely fluctuates in the window as well as the waveform is not symmetric, unlike the void fraction. The fluctuation and asymmetry indicated can be associated with the incoherent oscillations to the kinetic energy. For instance, if one bubble is expanding and another bubble is collapsing out of phase, their net contribution may not appear in the void fraction due to mutual cancellation but can appear in the energy. Overall, the averaged quantities, void fraction and the energy, are much smoother than the individual bubble dynamics. This can be trivially explained by the incoherence of violent collapse events among bubbles and coherence of the linear response against fundamental frequencies. Meanwhile, the quantitative nature of the coherent response is not predictable or easily analysable from these plots due to strong nonlinearity, especially under inter-bubble interactions with disordering factors including randomness of bubble position and polydispersity.

Figure 2. Evolution of representative quantities of a bubble cloud with $(B_0,s_d,N)=(0.5,0.1,100)$ excited at $A=20$, during a stationary state. (a) Radius of a representative bubble. (b) For the same bubble, $\xi$ obtained from the raw data, and those from the first and second dominant features extracted using WPCA-TD. (c) Void fraction. (d) Kinetic energy of the fluid induced by the bubble cloud.

Figure 3(a–c) shows the projected side-views of the three-dimensional bubble cloud of figure 2. The size of spheres can be interpreted as the mean energy of oscillations of bubbles at those locations in each feature. The overall distribution of the spheres of the first feature resembles that of the raw data, while the energy of the bubbles in the second feature is much smaller than that of the first feature. The plots therefore visually confirm that the first feature represents the most energetic mode of oscillations in the original data. This result also agrees with the observation of figure 3(b).

Figure 3. Projected side-views of the three-dimensional bubble cloud of figure 2. The size of the spheres denotes the root-mean-square amplitude of the velocity potential evaluated at the bubble surface at corresponding locations during the stationary state of oscillations, for (a) the raw data, (b) the first principal feature and (c) the second principal feature.

3.3. Error analysis

To show the effectiveness of the pre-processing, in figure 4, we plot the relative errors of the time-averaged kinetic energy of bubble clusters, one approximated using the original variable ($\langle \boldsymbol {q}^T\boldsymbol {T}(\langle \boldsymbol {r}\rangle )\boldsymbol {q} \rangle$) and the other using the transformed variables ($\langle \boldsymbol {\xi }^T\boldsymbol {P}(\langle \boldsymbol {\eta }\rangle )\boldsymbol {\xi } \rangle$), against the excitation amplitude. At $A<10^{-1}$, the error is nearly zero for both approximations. At $A>10^{-1}$, at which bubble dynamics become nonlinear, the error grows with $A$ for the former, while it remains small for the latter. This result confirms the improved approximation by the pre-processing.

Figure 4. Relative error of the mean kinetic energy of the fluid induced by clusters for the two approximations, $(\times ){:}\langle K\rangle \approx \langle \boldsymbol {q}^T\boldsymbol {T}(\langle \boldsymbol {q} \rangle )\boldsymbol {q}\rangle$ and $(\circ ){:}\langle K\rangle \approx \langle \boldsymbol {\xi }^T\boldsymbol {P}(\langle \boldsymbol {\eta }\rangle )\boldsymbol {\xi }\rangle$, against the excitation pressure amplitude with various polydispersities and values of $B_0$. Black, $(s_d,B_0)=(0.1,0.5)$; red, $(s_d,B_0)=(0.7,5.0)$.

3.4. Amplitude dependence of the key measures

The coherent dynamics critically depend on the excitation amplitude. Figure 5 shows the dependence of $\hat {\sigma }^2_1$, $\hat {S}_{vN}$ and $C$ against $A$ for various density and polidispersity parameters of clouds. This dependence is best highlighted in the result of a sparse, weakly polydisperse cluster ($B_0=0.5, s_d=0.1$) in figure 5(a). The relative importance of the first PC, $\hat {\sigma }_1^2$, decays nearly monotonically from 0.9 to 0.2 through three distinct regimes. For $A\lessapprox 0.2$, $\hat {\sigma }^2_1\approx 1$, meaning the entire energy is captured by the first PC. Then $\hat {\sigma }^2_1$ rapidly decays to 0.4 and stays nearly constant up to $A\approx 3$. At $A>3$, $\hat {\sigma }^2_1$ decays again and stays nearly constant around 0.2. The decay indicates the decrease in the coherence with increasing $A$, and can be explained by the excitation of the nonlinear oscillations and cavitation triggered at $A\approx 1$ and above. Here, $\hat {S}_{vN}$ has a profile vertically mirrored to $\hat {\sigma }_1^2$; $\hat {S}_{vN}$ increases from approximately 0.1 to 0.5 through the three regimes, indicating more partitioning of the energy into multiple PCs and increase of incoherence, by increasing $A$. The mirrored profiles of $\hat {\sigma }_1^2$ and $\hat {S}_{vN}$ suggest that these parameters are complementary.

Figure 5. The first PC-variance ($\hat {\sigma }_1^2$), the normalized von Neumann entropy ($\hat {S}_{vN}$) and the coherence measure ($C$), against the excitation amplitude with various initial density and polydispersity parameters, with $(B_0,s_d)= (a)$ $(0.5,0.1)$, (b) $(0.5,0.7)$, (c) $(5.0,0.1)$ and (d) $(5.0,0.7)$.

Remarkably, $C$ draws a square-well like profile, where $C\approx 1$ for both low ($A<0.2$) and large ($A>5$) amplitudes reaching a minimum value ($C\approx 0.5$) at $A\approx 1$. This counter-intuitive result suggests that the energy of the first PC is scaled by the mean-field parameter $B_d$ regardless of the increase of incoherence for large $A$, and the scaling is lost only for the intermediate range at $1< A<5$. Comparisons with the result of dense, weakly polydisperse clusters $(B_0,s_d)=(5.0,0.1)$ shown in figure 5(b) highlights the effect of the density of bubbles. In figure 5(b), $\hat {\sigma }^2_1$ takes values near unity at small $A$, and decays at $A\approx 1.0$ to approximately 0.4 and stays nearly constant at $A>5.0$. At $A>10$, $\hat {\sigma }^2_1$ slightly grows against $A$. Here, $\hat {S}_{vN}$ has a profile vertically mirrored to $\hat {\sigma }_1^2$. The features of $\hat {\sigma }^2_1$ and $\hat {S}_{vN}$ are similar to those observed in figure 5(a), except that the transition occurs at a larger range of $A$. The sudden decay and growth of $\hat {\sigma }^2_1$ and $\hat {S}_{vN}$ can likewise be associated with the linear-to-nonlinear transition of bubble dynamics which results in incoherence. Meanwhile, the positive shift of the transition range of $A$ indicates that the dense bubble clouds tend to behave more coherently than sparse bubble clouds, agreeing with the previous theory (d'Agostino & Brennen 1989; Maeda & Colonius 2019). Here, $C$ is greater than 0.5 at almost all values of $A$, and is relatively more insensitive against $A$ compared with in figure 5(a), indicating that the cluster dynamics is moderately controlled by $B_d$. Figures 3(c) and 3(d) respectively correspond to sparse and dense clusters with strongly polydisperse nuclei; $(B_0,s_d)=(0.5,0.7)$ and $(5.0,0.7)$. These plots show clear differences from those of weekly polydisperse clusters in figures 5(a) and 5(b). For the sparse clusters (figure 5c), $\hat {\sigma }^2_1$ mildly decays from 0.8 to 0.7 at $0.1< A<0.2$, sharply decays to the minimum of 0.2 at $A\approx 0.4$ and then mildly grows to 0.3 at $A=1.2$. The sharp decay resembles that in figure 5(b), while the slope is milder. Additionally, $\hat {S}_{vN}$ has a profile mirrored to $\hat {\sigma }^2_1$, but it is less symmetric than in figures 5(a) and 5(b). Here, $\hat {S}_{vN}$ is nearly constant around 0.2 at $0.1< A<2.0$ and draws a concave curve with its maximum of approximately 0.5 at $A\approx 5.0$ followed by a smooth decay to 0.3 at $A=20$. The $C$ almost constantly increases from 0.2 to $O(1)$ throughout the plot, indicating that the coherent dynamics is controlled by $B_d$, only at large $A$, unlike the weakly polydisperse case of figure 5(a). The overall trend of the plots in figure 5(b) is similar to that in figure 5(a), although the changes of variables against $A$ are milder in figure 5(b).

Overall, the mutual trends of variables at $A>5$ are similar between figures 5(a) and 5(c), and between figures 5(b) and 5(d), indicating the decreasing influence of the initial polydispersity in the nonlinear regime. At $A< O(1)$, the polydisperse clouds tend to have smaller values of $\hat {\sigma }^2_1$, and larger values of $\hat {S}_{vN}$ and $C$. This can be explained by the enhancement of incoherent dynamics induced by polydispersity.

3.5. Principal component spectrum

To gain deeper insights to the meaning of $C$, in figures 6(a)–6(d), we show the PC-variances obtained at $A=2\times 10^{-2}$, $1.2$ and $20$, for the first 10 PCs, obtained from the sparse, weakly polydisperse clouds blueplotted in figure 5(a). The insets show the evolution of the square-root of the normalized total energy ($\sqrt {K}$), and those of the first and the second PC-variances ($\sqrt {K}_1$ and $\sqrt {K}_2$), during the four periods of excitation in statistically stationary states. As expected, in the linear regime (figure 6a), the first PC occupies nearly the entire energy and $\sqrt {K}_1$ evolves at the fundamental frequency. In the transition regime (figure 6b), the first PC occupies 40 % of the energy and the rest is partitioned into the sub-dominant PCs with a smooth decay. Here, $\sqrt {K}$ evolves more chaotic than the linear regime, as expected due to the nonlinear response of bubbles. Both $\sqrt {K}_1$ and $\sqrt {K}_2$ evolve with similar quasi-periodic profiles. We interpret that both PCs represent the coherent part of the energy.

Figure 6. (a–d) PC spectral profiles at $A=2\times 10^{-2}$, 1.2, $8.6$ and 20. Insets show the evolution of the square root of the normalized total energy ($\sqrt {K}$, black) and those of the first (blue) and the second (red) PC-variances ($\sqrt {K_1}$ and $\sqrt {K_2}$), during the four periods of harmonic excitation, with $\hat {t}$ being a non-dimensional time $\hat {t}=tf$. The $y$-axis of each inset is normalized by the maximum value of $\sqrt {K}$.

Interestingly, in the nonlinear regime (figure 6c), energy partitioning is non-smooth; the first PC occupies 40 % of the total energy similar to the transition regime, but the rest of the energy is broadly distributed into the other PCs with much smaller amplitudes. The evolution of $\sqrt {K}$ is non-periodic with noisy, fine structures of spikes. These spikes are expected due to the incoherent collapse events. The evolution of $\sqrt {K_1}$ is, in contrast, highly periodic and somewhat resembles that of figure 6(b). The evolution of $\sqrt {K_2}$ is more chaotic and less smooth than $\sqrt {K_1}$. The difference between the first-PC and the rest of the PCs suggest that, in this regime, only the first PC captures a major coherent feature and the rest of the PCs represent more the incoherent dynamics as broadband noise. Figure 6(d) shows the result of the nonlinear regime with a stronger excitation amplitude ($A=20$). Overall, both the spectrum and the evolution of ${K}$ look similar to that of figure 6(c), other than that the amplitude of the first PC-variance is increased to 60 %. The resemblance of figures 6(c) and 6(d) indicates that the dynamical features identified from these two plots are common in the nonlinear regimes.

The resemblance of the evolution of the first PCs in the linear and the nonlinear regimes can explain the recovery of $C$ at $A>5$ in figure 6(a). Although the overall dynamics are much more chaotic in the nonlinear regime, the contribution of the coherent interactions to the system's energy effectively appears only in the first PC in both regimes and therefore the relative contribution of the interaction to the first PC is commonly scaled by $B_d$. This result also implies that the underlying coherence in the nonlinear regime represents the perfect correlation (synchronized oscillations) like that of linear, monodisperse clouds.

3.6. Evaluation of the coherence measure

To assess the variation of the coherence measure dependent on the three regimes, in figures 7(a)–7(c), we plot $B_1$ against $B_d$ for $O(10^3)$ clusters with various values of $s_d$, $N$ and $R_C$, at the three distinct excitation amplitudes ($A=2\times 10^{-2}$, $1.2$ and $20$). Appendix D summarizes the parameters used. With the weak excitation (figure 7a), bubble oscillations are in a linear regime and $B_d\approx B_0$. The data points are collapsed on the line of $C=1$ for $s_d=0.1$, while data points are scattered for the other values of $s_d$. This result is expected as $B_d$ was originally defined to scale the coherence of the monodisperse, perfectly correlated bubbles in the linear regime (d'Agostino & Brennen 1989; Zeravcic et al. 2011). With the intermediate excitation (figure 7b), the bubble oscillations are nonlinear. The data points are scattered from $C=1$, regardless of the value of $s_d$. In this regime, the result indicates that the coherence is lost due to the nonlinear dynamics, regardless of polydispersity. Surprisingly, with the strong excitation (figure 7c), the data points are collapsed on the line of $C=1$, meaning that the variance of the first PC is scaled by $B_d$, regardless of the parameters. This collapse is not observed for the second PC, regardless of the parameters (Appendix E). The results suggest that the scaling of $B_1$ is universal for cavitating clusters, which are typically excited at $O(1)$ MPa and above. Physically speaking, the scaling indicates that in cavitating clusters, the energy is partitioned into a single coherent mode of correlated (synchronized) oscillations and incoherent modes generalizing the aforementioned interpretation of figures 6(c) and 6(d). The coherent energy is controlled by the mean field originally derived for monodisperse, near-equilibrium bubbles. It is suggested that this partitioning and scaling are universal regardless of nuclei's polydispersity and the degree of nonlinearity. This finding explains the successful use of $B_d$ in characterizing and controlling seemingly disordered clusters in our previous studies (Maeda & Colonius 2019; Maeda & Maxwell 2021). The isolation of coherence implies that the details of the microscopic scale as well as of the many-body interactions represented by higher-order PCs could be modelled as fast variables, which force the macroscopic (relevant) scales in the form of noise. A rigorous way to corroborate this hypothesis might be the Mori–Zwanzig projection operator formalism (Mori 1965) by assuming that the fluctuations of fast and slow variables are uncorrelated. Although the choice of relevant slow variables may be unclear, the present approach is promising since the WPCA-TD can be seen as a prompter for more refined models of the evolution of slow variables.

Figure 7. $B_1$ against $B_d$ for $O(10^3)$ clusters with various values of $s_d$, $N$, $R_C$ and other physical parameters (Appendix D) in the (a) linear ($A=2\times 10^{-2}$), (b) transition ($A=1.2$) and (c) nonlinear ($A=20$) regimes.

4. Significance of the WPCA-TD to inform on cavitation physics

In general, PCA can be considered as a mathematical technique for data decomposition and the physical meaning of extracted features can be left to one's interpretation. To the authors’ knowledge, studies of complex fluid flows employing data decomposition and feature extraction techniques often rely on one's intuition to discuss the physical meaning of features. In fact, in § 3.5, the analyses of the PC spectra required rigorous interpretations based on previous knowledge on cavitation physics. Meanwhile, the correlation between $B_1$ and $B_d$ identified in the last section shows that the present WPCA-TD can directly extract $B_d$ from data without additional signal processing like spectral filtering, which are often required for Fourier-based analysis, despite the presence of noisy and incoherent features. To generalize, WPCA-TD is not only a tool for data decomposition and feature extraction, but also can be used to directly inform on the non-dimensional number that controls the coherent physics of interests. This informative aspect of the WPCA-TD can provide a meaningful shortcut to address the physics of cloud cavitation.

To perform WPCA-TD, bubble dynamics data can be obtained from experimental measurements as well as other numerical approaches (e.g. Kameda & Matsumoto 1996; Maeda & Colonius 2018; Pishchalnikov et al. 2019). The cluster's shape can be arbitrary. Although a bubble's translation and deformation are neglected in the present numerical experiments, WPCA-TD can incorporate dynamical variables controlling these effects (e.g. Ilinskii, Hamilton & Zabolotskaya 2007; Murakami, Gaudron & Johnsen 2020). Physically meaningful data decomposition requires the fine temporal resolution of individual bubble dynamics. In practice, such information would be difficult to obtain in experiments except for a small number of bubbles in a highly controlled environment. We thus emphasize that WPCA-TD would primarily be useful to process fine temporal resolution of numerical data.

5. Conclusion

In conclusion, to corroborate the coherent dynamics of bubble clusters, we introduced and demonstrated WPCA-TD, a method of PCA to comprehensively decompose the time series of Lagrangian bubble dynamics data into coherent dynamical features, in way similar to the modal decomposition of Eulerian flow field data. The data are pre-processed such that the PC-variance of the features becomes consistent with the hydrodynamic potential energy induced by bubble oscillations that is predicted by the R–P equation. By analysing simulation datasets of clusters under harmonic excitation, we demonstrated that the coherent energy and the degree of incoherence are respectively quantified by the variance and the spectral entropy. The coherence was lowered by disorders induced by the nuclei's polydispersity and nonlinear response of bubbles, as expected. Meanwhile, in cavitating regimes, underlying, correlated modes of oscillations were isolated in a single dominant feature. The variance of this feature was found to be controlled by the previously identified mean-field parameter, $B_d$, regardless of the disordering factors, indicating that the underlying coherent dynamics may be universal in cavitating clusters. These results suggest that the method can provide a simplified and comprehensive representation of complex bubble dynamics. Analogous to the use of POD for reduced-order modelling (ROM) of various flows, the dynamical features extracted by PCA may be used for ROM, without directly solving many-body interactions. Such a model may be promising for controlling cavitation and nucleation without tracking individual nuclei.