A Study of the Influence of Thermodynamic Parameters on the Dynamics of Cavitation Bubbles

1. Introduction

In this work, we consider the patterns of bubble behaviour within a cluster, taking into account the mutual influence of the effects accompanying the cavitation development process. Regardless of whether a dynamic effect at a local point in the fluid is a consequence of the collective influence of a large number of bubbles or is solely caused by the action of a single bubble (for example, the formation of a cumulative microjet or the emission of a spherical shock wave by a collapsing bubble), all bubbles participate directly or indirectly in its creation. The difference between the behaviour of a single bubble within a bubble cloud and the behaviour of an isolated bubble in an infinite volume of liquid is primarily due to the influence of its nearest neighbours on the bubble’s evolution. As they develop, these neighbours create their own pressure field in the vicinity of the bubble, which differs from the external pressure currently acting on the system. Therefore, the behaviour of an individual bubble within a cluster under variation of the external pressure is more complex than that represented in the single-bubble dynamics model.

Compared to the large volume of work dedicated to the study of single bubble dynamics, publications related to the study of the collective behaviour of bubbles and the evolution of an individual bubble within a dynamically developing ensemble are less prevalent. These studies can be broadly divided into three main groups.

The first group focuses on phenomena related to the dynamic effects of bubbles, for example, [1,2]. The authors show that hydrostatic pressure is the dominant factor in controlling cavitation activity, but due to the inertia of the cavitation cluster, dynamic effects are difficult to predict for clustered bubbles. A cavitation system possesses a specific bubble density, which can be manipulated by changing the hydrostatic pressure to enhance dynamic effects. This could serve as a promising approach for developing technologies, such as material synthesis or processing [3,4].

The core idea of the second group of studies focuses on the fact that bubbles can accumulate into clusters and thereby intensify the cavitation effect. However, many authors point to the ambiguity of such effects [5,6]. The behaviour of bubbles in clusters differs significantly from the dynamics of single bubbles. Their dynamics are mutually influenced by neighboring bubbles, leading to attraction or repulsion. These effects can be explained by the presence of the secondary Bjerknes force.

Some authors [7,8] indicate that the radial oscillations of a cavitation bubble are suppressed by larger neighbors when they are in close proximity, and the magnitude and direction of the secondary Bjerknes force change, leading to a transition between mutual attraction and repulsion of the bubbles. Authors in [9] assert that the suppression of bubble pulsations increases with increasing bubble density. If the bubble density exceeds a critical value, a repulsive region emerges, preventing their coalescence. For example, at a pressure of 3 MPa and a bubble density up to 1 · 10¹³ m⁻³, distinct zones of attraction appear, and as the hydrostatic pressure increases, the attraction zone gradually expands towards a higher density region. In fact, increasing hydrostatic pressure can increase the stable bubble density of the cluster, which will intensify the cavitation effect. It can be argued that while there are numerous publications devoted to experimental studies of single bubbles [10,11], there are significantly fewer studies on the mutual influence of double bubbles [12,13], and very few experimental reports focus on clusters, which constitute the third approach to research. Furthermore, the question of what causes the dynamic effect—whether it is the cumulative effect of the cluster or the influence of individual bubbles surrounding the cluster—remains unexplained. Modelling of dynamic effects [14,15,16] is more widely applied. Among the approaches presented in publications for describing cavitation processes, several well-established types of models can be distinguished.

1. Multiscale models [17,18,19,20]. This approach can describe a compressible multiscale cavitation region where the dynamics of large vapour cavities and microbubbles are modeled, taking into account the average compressibility. The cavitation cluster is treated as a liquid-vapor mixture continuum in the Eulerian sense. These models are developed by coupling an Eulerian medium and functional dependencies, expressed in terms of velocities and radii, which are tracked in the Lagrangian sense using various approaches. For a stable transition between the Eulerian and Lagrangian frameworks, the transition criterion can be defined based on a parameter that accounts for the non-uniformity of the domain discretization. In hybrid models, the onset of cavitation is first estimated in low-pressure regions, which limits predictions using Eulerian mass transfer models.

2. Models relating fluid characteristics to thermodynamic parameters [21,22,23,24]

This method, based on equations of state that link fluid characteristics to thermodynamic parameters, is limited only to very specific problems unrelated to either bubble growth dynamics or the dynamic effects of cavitation. In general, these models represent a simplified representation of real processes and are usually solved for an averaged substance whose properties reflect the weighted average of the volume fraction of the two media. Another approach involves modelling cavitation using real gas polynomial equations of state to accurately predict the properties of water in any state. Modelling complex dynamic effects relies on numerical methods which, although accurate, require significant computational resources and databases, and are not applicable for rapid prediction in practical scenarios.

3. Phase boundary modelling [25,26]

4. Modelling using transport equations [27,28,29,30]

Such models accurately predict the dynamics of the bubble cloud in the initial stage by considering bubble interaction. However, in the final stage, the inertia effect, represented as a second-order derivative, plays a more significant role and shows a higher collapse speed.

5. Predicting bubble dynamics using the Rayleigh–Plesset equation [31,32,33]

The Rayleigh–Plesset equation accurately describes bubble oscillations in incompressible liquids but neglects energy losses due to acoustic frequencies, as it was based on incompressible fluids. Therefore, this equation becomes inapplicable in cases where energy loss cannot be neglected. Nevertheless, this approach to describing cavitation processes is fundamental in the development of cavitation modelling, including the description of bubble cluster dynamics.

Still, the idea of modelling the impulse of a bubble cluster has not received significant attention in the literature, despite the fact that the mathematical interpretation of cavitation phenomena may be highly attractive for a large number of problematic cases. A critical question remains unanswered by existing models: how do neighbouring bubbles affect the dynamic effects within a cluster?

The principles of the theoretical approach to solving this problem are outlined in various works, e.g., in [34,35,36]. One simplified approach to modelling cavitation clusters involves describing the dynamics of an aggregate of bubbles using the cell model approximation, which allows the evolution of a monodisperse ensemble to be linked to the behaviour of a single bubble. Despite the small number of publications, the available information on the results of experimental and theoretical studies indicates the complexity and extraordinary nature of the processes being studied. These models are based on the fundamental assumption that at all stages of cloud evolution, the bubble sizes are much smaller than the distance between the nearest neighbors.

The results of theoretical calculations of pressure fields and flow velocities, even for two interacting bubbles, show poor agreement with experiments [37,38]. According to researchers, this is explained by the fact that the bubble dynamics equations are derived with the assumption of potential flow in its vicinity, although the real velocity field is neither symmetrical nor potential. Accounting for all the features of bubble cluster dynamics will likely remain an unresolved problem, but the development of simplified methods for assessing the behaviour of a cavitation cloud is a promising task.

2. Bubble Dynamics Within a Cluster using the Cell Model Approximation

Without detailing the complex spatial configurations of velocity and pressure fields within the inter-bubble regions of a dynamically evolving cluster, the general nature of the fluid flow can be conceptualized as follows. As the bubbles expand or contract within an elementary volume of the two-phase system, an unsteady liquid flow is induced in the space between the bubbles. Relative to any given point within this volume, the liquid flow exhibits inherent asymmetry, consequently leading to a non-symmetrical pressure field. This scenario is further complicated by the continuous temporal changes in the distribution of bubble sizes and their mutual distances throughout the volume. The microflows observed in the inter-bubble space are fundamentally vortical, which causes additional localized pulsations in the fluid’s velocity and pressure. When the amplitude and frequency of these pulsations are sufficiently large, a specific collective effect emerges, whereby the bubble cloud intensifies hydromechanical and heat and mass transfer processes.

The intricate pressure distribution surrounding each bubble complicates the analytical description of its dynamics. Under these conditions, the fundamental assumption of spherical symmetry for the velocity and pressure fields in the bubble’s vicinity becomes inaccurate, necessitating the inclusion of higher-order effects. In contrast to a single, isolated bubble, where developing reliable and consistent models presents no significant challenges, describing the behaviour of an individual bubble within a cloud, and the cluster as a whole, is fraught with difficulties of a fundamental nature.

In the review section of this article, attention is drawn to the dependence of bubble behaviour on their concentration. The use of cell models reflects this dependency by adopting two distinct approaches: one for dilute systems [39,40] and another for concentrated clusters [41,42].

These foundational principles of the cell model enable the utilization of established single-bubble dynamics equations as a first approximation for calculating the averaged ensemble parameters. Subsequent calculation steps then involve sequentially incorporating higher-order effects. Crucially, when a bubble is surrounded by other growing bubbles, the pressure fields generated by its nearest neighbours mutually overlap. Consequently, the character of the pressure distribution in the vicinity of each bubble, p_b=f(r), deviates significantly from the distribution described by the widely recognized single-bubble equations [43,44].

$$p_{\mathrm{b}}\left(\Delta ,\tau \right)=p_f+\left[p_b-p_f+{\displaystyle \frac{1}{2}}·{\rho }_b·u_R^2·\left(1-{\displaystyle \frac{R^3}{r^3}}\right)-{\displaystyle \frac{2\sigma }{R}}-{\displaystyle \frac{4\mu ·u_R}{R}}\right]·{\displaystyle \frac{R}{r}}$$

(1)

where Δ is the distance between the bubbles.

The dynamic evolution of each bubble is therefore not governed by the external pressure p_f, acting on the cluster as a whole, but rather by an average pressure within the inter-bubble space, denoted as p_b, which is characteristically greater than the external pressure p_b>p_f. In the approximation utilized by the cell model, the external boundary condition for the pressure of an individual bubble should be applied not at infinity, but at a distance corresponding to the average spacing between adjacent bubbles—specifically, at the radius of the spherical cell within which the bubble is centrally located (Figure 1).

Analysing the available literature, the fundamental principles of the cell model, as applied to the dynamics of vapour bubbles within a cluster, can be summarized. If the volumetric vapour content (void fraction) of the two-phase bubble medium, γ(τ), and the mean bubble radius, R(τ), are known, then the volume allocated per bubble is given by 4/3πR³γ(τ). This volume can be conceptually treated as an elementary spherical cell containing a single bubble at its center. The outer radius of this cell, Ψ, is defined as Ψ=R/γ^0.3, and the thickness of the liquid layer within the cell is Ψ=R+R(1-γ^0.3)/γ^0.3. It is postulated that the liquid motion within the boundaries of the cell, induced by the expansion or compression of the central bubble, is strictly radial. Crucially, there is no mass exchange across the external boundary of the cell with the rest of the cluster. Consequently, the mass of the substance contained within the cell is conserved (remains constant). Within the framework of the cell model, it is further postulated that the following conditions are satisfied at the outer boundary:

$${{\displaystyle \frac{dm}{dr}}|}_{\mathrm{r}=\psi }=0;{{\displaystyle \frac{dT}{dr}}|}_{r=\psi }=0.$$

(2)

This condition facilitates the incorporation of temperature variations within the cell’s liquid due to heat and mass transfer with the central bubble. The temperature at the external surface of the cell only begins to change once the thickness of the thermal layer surrounding the bubble, θ(τ), exceeds the radius of the cell [36].

The behaviour of each bubble is treated as that of an isolated single bubble confined within the volume of a closed cell. As the cell expands or contracts, it undergoes a force interaction with the surrounding portion of the two-phase system, which dictates the pressure exerted on the cell’s external surface. This pressure, p(Ψ)=p_k, is then employed as the external boundary condition during the derivation of the equation governing the fluid motion within the cell. The resulting pressure distribution in the liquid within the confines of the cell, according to the Lagrange–Cauchy integral [45], is given by the expression:

$$p\left(\theta \right)=p_{\mathrm{k}}-{\rho }_f·\left({\displaystyle \frac{\partial \varphi }{\partial \tau }}-{{\displaystyle \frac{\partial \varphi }{\partial \tau }}|}_{r=\psi }+{\displaystyle \frac{u^2}{2}}-{{\displaystyle \frac{u^2}{2}}|}_{r=\psi }\right).$$

(3)

By applying the boundary condition for the pressure at the bubble interface, p(R)=p_R, and subsequently defining the parameter values on the right-hand side of Equation (3) using the following formulas:

(4)

Consequently, by applying the condition for pressure at r=R (the bubble interface) and at r = Ψ (the cell boundary), we derive the equation for the radial motion of the fluid within the cell:

$${\displaystyle \frac{d{u}_R}{d\tau }}={\displaystyle \frac{p_R-p_k-0,5{\rho }_b{u}_R^2\left(3-4{\gamma }^{0,3}+{\gamma }^{1,3}\right)}{{\rho }_{f}R\left(1-{\gamma }^{0,3}\right)}}.$$

(5)

This resulting equation differs significantly from the standard equation of motion developed for a single bubble in an infinite liquid volume [46], due to the inclusion of the cell boundary effects.

$${\displaystyle \frac{d{u}_R}{d\tau }}={\displaystyle \frac{p_R-p_f-1,5{\rho }_b·u_R^2}{{\rho }_{b·}·R}}$$

(6)

This difference stems from two main factors: the presence of correction terms on the right-hand side of Equation (5) that account for the non-negligible volume occupied by the bubbles, and the replacement of the external pressure p_f with the currently unknown pressure at the cell boundary, p_k. As the vapor volume fraction γ→0, the cell radius Ψ tends towards infinity Ψ→∞, and the pressure at the cell’s outer boundary p_k converges to the external pressure p_k→ p_f . Consequently, Equation (5) simplifies to Equation (6) under dilute conditions. It is worth noting that for the initial growth stages of small bubbles or the final collapse phase, γ is small, rendering the correction terms in Equation (5) of secondary importance. However, the cell model approximation becomes inherently flawed for comparatively high values of the two-phase system’s vapour content (γ >0.5).

The principal distinction between Equation (5) and Equation (6) lies in their boundary conditions: the behaviour of the isolated bubble in Equation (6) is determined by the known, externally applied pressure p_f. Conversely, the behaviour of a bubble surrounded by others is dependent on the pressure p_k(τ) at the external cell boundary, a value that is not known a priori and is, in turn, determined by the intensity of expansion or compression across the entire bubble cluster. Finally, the precise magnitude of the pressure in the liquid phase p_R, as well as the quantitative influence of capillary and viscous forces on this pressure, remain open questions within this modelling framework.

3. Modified Model

We consider an alternative approach to modelling bubble behaviour within a cluster, one that is based on averaging the pressure solely within the liquid’s inter-bubble space. We assume that within a liquid volume, initially at temperature T_f0, there are uniformly distributed vapor bubbles of identical size, with an initial radius R₀. These bubbles are considered to be in thermodynamic equilibrium with the liquid. The initial vapour pressure inside each bubble is equal to the saturated vapor pressure at the liquid temperature, p_sat(T_f0). The initial pressure in the liquid is uniform throughout the entire volume and is determined by the condition: p_f0= p_sat(T_f0) - 2σ/R. The number of bubbles per unit mass of the system, k_b, is assumed to be known and remains constant throughout the system’s subsequent evolution, thereby excluding the effects of bubble coalescence and breakup. At time τ=0, the external pressure acting on the system, initially p_f = p_f0,, begins to change according to a known law, initiating the growth or compression of all bubbles. For a monodisperse system, the volumetric vapour content (void fraction) is related to the size and concentration of the bubbles and is determined by the expression:

$$\gamma ={\displaystyle \frac{4\pi R^3{\mathrm{k}}_b{\rho }_b}{3+4\pi R^3{k}_b{\rho }_b}}.$$

(7)

Within the vicinity of each bubble, specifically within the boundaries of the elementary cell of radius $\psi \left(\tau \right)=R\left(\tau \right)/{\gamma }^{0,3}$, a spherical pressure distribution, p_f(r,τ), is established. This pressure field can be determined using Equation (8):

$${\displaystyle \frac{\partial \varphi }{\partial \tau }}+{\displaystyle \frac{u^2}{2}}+{\displaystyle \frac{p}{{\rho }_b}}=f\left(\tau \right),$$

(8)

under the condition $p_f\left(R,\tau \right)=p_R\left(\tau \right)$.

When determining the function f(τ) on the right-hand side of Equation (8) using the boundary condition $p\left(\psi \right)=p_k$, the pressure distribution obtained takes the form of Equation (3). Conversely, if we utilize the boundary condition $p\left(R\right)=p_R$to determine f(τ), the pressure distribution inside the cell adopts the following form:

(9)

By utilizing the relationships defined in Equation (4) for the calculation of the quantities on the right-hand side of Equation (9), we then reduce this equation to the final form (Equation (10)):

$$p_f\left(r,\tau \right)=p_R-{\rho }_b\left(R{\displaystyle \frac{d{u}_R}{d\tau }}+2u_R^2\right){\displaystyle \frac{1-R^4}{r^4}}+0.5{\rho }_b{u}_R^2{\displaystyle \frac{1-R^4}{r^4}}.$$

(10)

By setting the radial coordinate equal to the cell radius, r=Ψ, Equation (10) can be re-expressed as follows:

$$p_f\left(r,\tau \right)=p_R-{\rho }_b\left(R{\displaystyle \frac{d{u}_R}{d\tau }}+2u_R^2\right)·\left(1-{\gamma }^{0.3}\right)+0.5{\rho }_b{u}_R^2·\left(1-{\gamma }^{1.3}\right).$$

(11)

Considering the right-hand side of Equation (5), we can derive the equation for calculating the pressure in the vicinity of the bubble:

$$p_b\left(r,\tau \right)=p_R-{\rho }_b\left({\displaystyle \frac{p_R-p_k+0.5{\rho }_b{u}_R^2\left(1-{\gamma }^{1.3}\right)}{1-{\gamma }^{0.3}}}\right)·\left(1-{\displaystyle \frac{R}{r}}\right)+0.5{\rho }_b{u}_R^2\left(1-{\displaystyle \frac{R^4}{r^4}}\right).$$

(12)

The expression for the averaged pressure inside the cell is as follows:

$$p_m\left(\tau \right)=p_k+{\displaystyle \frac{{\rho }^{0.3}\left[p_R-p_k+0.5{\rho }_b{u}_R^2\left(1-{\gamma }^{1.}\right)\right]}{1-{\gamma }^{0.3}}}·\left[{\displaystyle \frac{1.5\left(1-{\gamma }^{0.6}\right)}{1-{\gamma }^{0.6}}}-1\right]-0.5{\rho }_b{u}_R^2·\left({\displaystyle \frac{3\gamma ·\left(1-{\gamma }^{0.3}\right)}{1-{\gamma }^{0.3}}}-{\gamma }^{1.3}\right).$$

(13)

The following algorithm is employed to determine the pressure at the external cell boundary p_k, the pressure at the bubble boundary p_R, the mean pressure in the liquid p_m, as well as the bubble radius R and the velocity of its surface u_R at each stage of the two-phase system’s evolution.

Algorithm for dynamic parameter calculation. At the initial time τ=0, the external pressure acting on the system instantaneously changes from its initial equilibrium value p_f=p_f0 to some constant prescribed value p_f<p_f0. At this precise moment, the pressure within the inter-bubble space, including the external cell boundary, is equal to this external pressure. By substituting this value of the pressure into Equation (5) as the parameter p_k,, we determine the magnitude of the liquid acceleration at the bubble interface within the cell du_R/dτ. Subsequently, we calculate the new values for the system’s governing parameters, including R, u_R, γ, and p_R. Substituting these newly calculated quantities into the right-hand side of Equation (13) allows us to compute the value of the liquid pressure averaged over the cell volume, p_m. The pressure value p_m found in this manner is then utilized in the next calculation step as the cell boundary pressure p_k. This value of p_k is then inserted back into equation (5) to subsequently determine the new values for du_R/dτ, u_R and R. Other governing system equations are solved simultaneously with equation (5) to compute the new values of the remaining parameters. This procedure is repeated iteratively at each time step of the numerical solution for the bubble dynamics equations. If the external pressure p_f acting on the system does not change instantaneously but over a finite time interval Δτ according to some known function, for example:

$$p_f=p_{f0}-\left(p_{f0}-p_b\right)\sin \left(0.5\pi \tau /\Delta \tau \right)$$

(14)

the change in external pressure, dp_f, within the time interval 0<τ<Δτ is determined by the expression:

$$d{p}_f=-0.5\left(p_{f0}-p_b\right)·\cos \left(0.5\pi \tau /\Delta \tau \right)·\pi d\tau /\Delta \tau .$$

(15)

At each calculation step, the pressure p_f is computed according to the specified scheme, and the obtained value is then adjusted by the increment dp_f, substituting the time step interval as dτ. After the change in external pressure is completed τ> dτ, p_f remains constant p_f=const and dp_f=0. This iterative procedure allows for the calculation of the time-dependent pressure at the cell boundary, p_k, and the mean pressure within the inter-bubble space of the ensemble, p_m, for any variation in the external pressure.

When the vapor content of the bubble system changes, the liquid temperature must vary due to the release or absorption of heat associated with phase transitions occurring at the expanded interfacial surface. The mean temperature of the liquid phase in the bubble ensemble, T_m, can be identified with the mean temperature of the liquid within the elementary cell. In the vicinity of each bubble, within the confines of the non-stationary thermal layer, a temperature distribution T_f(r,τ) exists.

The mean liquid temperature within the cell, T_m, at any given time is determined by integrating the functional expression T_f(r,τ) over the volume of the liquid in the cell and subsequently dividing by that volume. As previously noted, the cell model inherently assumes that each elementary cell is thermally isolated. Therefore, the amount of heat transferred through the bubble surface due to phase transition and heat exchange, 4πR²(JQ+q) [6], uniquely determines the change in the heat content of the liquid in the cell and, consequently, its mean temperature. Thus, the energy balance can be formulated as:

$${\displaystyle \frac{d}{d\tau }}\left[{\displaystyle \frac{4}{3}}\pi {\rho }_{b}c\left({\psi }^3-R^3\right)·T_m\right]=4\pi R^2\left(JQ+q\right).$$

(16)

By substituting the expression for the void fraction γ from Equation (7) into Equation (16) and performing the necessary mathematical transformations, we arrive at a simplified equation for determining the mean temperature of the liquid in the inter-bubble space of the ensemble:

$${\displaystyle \frac{d{T}_m}{d\tau }}={\displaystyle \frac{{\mathrm{k}}_{\mathrm{b}}\left(JQ+q\right)}{c}}.$$

(17)

A new initial parameter introduced into this model is the concentration of bubbles per unit mass, k_b. The system of equations, transformed in this manner and incorporating the proposed iterative scheme for determining the pressure within the inter-bubble space, can be effectively utilized to predict the dynamic behaviour of a cluster of monodisperse bubbles when subjected to variations in external pressure.

4. Liquid Flow in the Inter-Bubble Space of a Cluster

Cell models offer a capability to estimate the average pressure and void fraction within the inter-bubble space of the system at any given moment. This, in principle, allows for the investigation of steady and unsteady two-phase bubble flows in pipes and channels, as well as tracking the evolution of an individual bubble within such a system. However, these models are fundamentally incapable of detailing the hydrodynamic and heat and mass transfer phenomena that occur at the micro-level within the inter-bubble space. The model achieves pressure or void fraction averaging over a volume that significantly exceeds the volume of a single bubble. Consequently, the cluster is effectively treated as a homogeneous medium, with bubbles represented merely as singularities.

To accurately describe the behaviour of bubbles at the micro-level and to model the velocity and pressure fields in the inter-bubble space, several determining factors influencing the hydrodynamic environment in the liquid phase during intense bubble surface motion must be considered.

1. The bubbles within the cluster vary in geometrical size, shape, velocity, and the acceleration of their radial surface motion. This disparity necessitates an individualized description of each bubble’s evolution.

2. The force interaction between dynamically developing bubbles results in their chaotic translational movement. This, in turn, causes continuous changes in the velocity and direction of the microflows over time, leading to pulsations of local liquid pressure.

3. The dynamic interaction between neighbouring bubbles causes their continuous deformation. This further induces oscillations on the surface and significantly complicates the hydrodynamic environment in the vicinity of each bubble.

In essence, the physical model that accounts for these complex factors is sufficiently straightforward and apparent. For instance, the non-uniformity of the bubbles within the ensemble can be addressed by specifying their initial size and spatial distribution. In estimating the force interaction between bubbles, the expression for the Bjerknes force, which arises between two oscillating bubbles [47,48], can be used as a first approximation. The nature of the surface deformation is managed by the condition of pressure equality at the interfacial surface, allowing the calculation of the bubble’s radius of curvature in any direction [49]. However, the major modeling difficulty lies in the excessive duration of numerical computations, even for a single time step. Consequently, a detailed analysis of the behaviour of even two interacting bubbles proves to be a highly complex and labor-intensive task [50].

To analyze the nature of microflows, pressure fields, and velocity fields in the liquid phase of a dynamically evolving bubble system, we consider a model as a first approximation based on the following key assumptions:

1) the influence of viscosity is manifested only at the interfacial surface;

2) bubbles retain a spherical shape during their evolution, this is justifiable for relatively small bubbles, where capillary pressure counteracts the external forces from other bubbles, provided that γ<0.55.

3) the force interaction between neighbouring bubbles, caused by the radial motion of their surfaces, is accounted for, which results in the mutual translational movement of bubbles within the ensemble;[47,48

4) the liquid flow in the inter-bubble space is potential.

We treat the expanding or compressing bubbles as discretely distributed dynamic and thermal sources (or sinks) within the liquid. Based on assumption (4) (potential flow), we utilize the superposition principle in accordance with potential theory: the resultant flow is the sum of flows from each source, and the magnitude of the scalar potential of the velocity field at an arbitrary point in the liquid is the sum of the potentials generated at that point by all sources.

Consider a cluster of vapour bubbles comprising k_b elements. We place the origin of the Cartesian coordinate system at some point in the center of the bubble system. If the center coordinates (x,y,z), radii R_k, and surface expansion velocity u_R, of each of the k_b bubbles are known at a specific instant in time, the total field potential at an arbitrary point in the liquid with chosen coordinates, according to the superposition rule of potentials, is determined by the formula:

$$\varphi \left(x,y,z\right)={\displaystyle \sum _{i=1}^{\mathrm{k}}-\frac{u_{Ri}{R}_i^2}{{\varpi }_i}},$$

(18)

where ${\varpi }_i=\sqrt{\Delta x^2-\Delta y^2-\Delta z^2}$.

The components of the resultant velocity vector at the point (x,y,z) are determined by the following expressions:

$$u_x={\displaystyle \sum _{i=1}^k\frac{u_{Ri}{R}_i^2\Delta x}{{\varpi }_i^3}};u_y={\displaystyle \sum _{i=1}^k\frac{u_{Ri}{R}_i^2\Delta y}{{\varpi }_i^3}};u_z={\displaystyle \sum _{i=1}^k\frac{u_{Ri}{R}_i^2\Delta z}{{\varpi }_i^3}},$$

(19)

The spatial direction of the velocity vector is subsequently determined using the standard formulas of analytical geometry, based on the known values of its projections:

$${\displaystyle \frac{\partial \varphi }{\partial \tau }}={\displaystyle \sum _{i=1}^k-\left(\frac{d{u}_{Ri}}{d\tau }·\frac{R_i^2}{{\varpi }_i}+\frac{2R_i{u}_{Ri}^2}{{\varpi }_i}\right)},$$

(20)

where du_Ri/dτ is the acceleration of the radial surface motion of the individual bubble, as determined by Equation (6).

By substituting the value of the partial derivative of the potential with respect to time, $\partial \varphi /\partial \tau$, from Equation (20) and the velocity vector from Equation (19) into the Lagrange–Cauchy integral (Equation (8)), the magnitude of the pressure at a local point (x,y,z) of the microflow can be determined. The function f(τ)on the right-hand side of Equation (8) is determined by imposing the condition that at a sufficiently large distance from the bubble cloud, the pressure in the liquid is constant and equal to p_f,, and both the flow velocity and the potential are zero. Under these far-field conditions, the local pressure within the inter-bubble space inside the cloud is determined by the following equation:

$$p_b\left(x,y,z\right)=p_f-{\rho }_b\left({\displaystyle \frac{\partial \varphi }{\partial \tau }}+{\displaystyle \frac{{\stackrel{\rightharpoonup }{u}}^2}{2}}\right)$$

(21)

by substituting the values of $\stackrel{\rightharpoonup }{u}$ and $\partial \varphi /\partial \tau$from Equations (20) and (19), respectively, into its right-hand side.

Since the point (x,y,z) is arbitrarily chosen, this scheme allows for the estimation of the velocity and pressure magnitudes at any point within the inter-bubble space. Consequently, the instantaneous distribution of velocities and pressures throughout the two-phase system under consideration can be determined. It is essential to ensure that the point where $\stackrel{\rightharpoonup }{u}$ and p_b are determined does not fall within the region currently occupied by the bubble’s gas phase. Determining the pressure field in the liquid between the bubbles makes it possible to average the pressure over the volume of each individual cell (which are not necessarily identical). This averaged pressure is then used in the next time step to solve the single-bubble dynamics equation for each of the k_b bubbles, thereby determining the new values for R, u_R, du_R/dτ, and the new spatial coordinates x,y,z. Subsequently, the entire aforementioned procedure for calculating the velocity and pressure fields is performed again.

In this model, the correct calculation of the dynamics of each individual bubble within the finite-volume ensemble is of fundamental importance. As previously indicated, equation (5) must be used as the equation of motion to determine the value of the radial acceleration du_R/dτ on the right-hand side of equation (20).

The main question centers on the appropriate values for the parameters p_k and γ to be substituted into these equations. During the derivation of equation (5), it was established that the void fraction γ in this equation characterizes the cell radius Ψ=R/γ^0.3 surrounding the specific bubble, while p_k represents the liquid pressure averaged over the volume of that cell. This bubbles located on the periphery of the cloud are in constant contact with an unlimited volume of the liquid phase, where the pressure p_f=const. For these peripheral bubbles, the cell radius Ψ is assumed to be half the distance from their centers to the center of their nearest neighbours, which determines the corresponding value of γ. The pressure calculated via equation (13) is substituted as p_k into Equation (5). These peripheral bubbles essentially behave as isolated bubbles in an infinite liquid volume throughout the process. However, the pressure that these bubbles create in the outer layer of the cluster, calculated using Equation (13), is hypothesized to be transmitted as the cell boundary pressure p_k for all bubbles situated in the next adjacent spherical layer—the layer bordering the outer cloud layer. This pressure p_k, created by the bubbles in the second layer, is then passed on as the boundary condition for the bubbles in the third spherical layer, and so on. Consequently, the averaged pressure within a spherically shaped bubble cloud will vary from the periphery to the center, and the pressure gradient along the cloud’s radius depends on the intensity of expansion or compression of the bubbles within each spherical layer. This macro-level procedure for determining pressures within a finite ensemble allows for the calculation of the current parameters of each bubble as a function of its distance from the cloud center and its initial size. Subsequently, at each computational time step, the known values of the radial acceleration du_R/dτ and velocity u_R for every bubble enable the estimation of the velocity, pressure magnitudes, and direction of the microflow at any point in the inter-bubble space, utilizing the field equations (18)–(21). In summary, the presented equations for both unlimited and finite numbers of vapour bubbles, combined with the single-bubble dynamics equations, provide a framework for predicting the behaviour of bubbles within a cluster and the cluster’s overall behaviour under varying external pressure, as well as for determining the velocity and pressure fields in the inter-bubble space and their temporal evolution under any operational parameters.

5. Discussion

Based on the modified model, an analysis of the behaviour of vapor bubbles in water was performed under different operating modes and various external pressure change conditions. Below, some of the prediction results are presented as examples to illustrate the model’s capabilities

5.1. Dynamics of Bubble Growth and Collapse in a Cluster

Within the framework of the modified model, the evolution of a monodisperse vapour bubble cloud was considered for known initial values of the liquid temperature T_f0 and initial pressure и p_f0. At the initial instant, all bubbles with radius R₀ and concentration k_b are in equilibrium with the liquid. At the initial time τ=0, the external pressure acting on the system begins to change. In analyzing various experimental methodologies related to the rapid change of pressure in a liquid volume [51,52,53], we established that, regardless of the method used to change the pressure in the system, the nature of the actual curve p_f=f(τ) is well approximated by a cosinusoidal dependence, taking the form:

(22)

where P represents the constant pressure established in the liquid far from the bubble as a result of the change in external pressure, and Δτ is the duration of this change. After the specified time interval Δτ, the external pressure assumes a new constant value P.

Figure 2 illustrates how the mean pressure in the liquid within the inter-bubble space of the cluster, p_m=p_k, changes following an instantaneous drop in external pressure from an initial equilibrium value of 60 kPа to 15 kPa. The pressure in the liquid rapidly increases, asymptotically approaching the saturated vapour pressure p_sat(T_f), which corresponds to approximately 30 kPa under the conditions presented in the figure. The figure also includes a graph showing the pressure pulsations within a cavitation cloud, based on the data presented in the study [54].

The presented data suggests that the pressure within the cluster undergoes an averaging process, but the dynamic effects originating from the bubbles at the cloud boundaries can be substantially greater than those observed inside the cluster. These internal effects within the cavitation cloud become more prominent as the bubble concentration increases. Figure 3 illustrates the change in the liquid pressure inside a monodisperse cloud of collapsing bubbles following an instantaneous increase in the external pressure from the equilibrium value of 70 kPa to 150 kPa. Due to the intense initial compression of the bubble ensemble and the associated stretching of the liquid, the pressure in the inter-bubble space initially drops sharply to nearly 70 kPa. However, the rapid compression leads to a significant increase in the vapour pressure inside the bubble, which in turn causes the pressure to rise in the liquid layers adjacent to the surface. This elevated pressure then propagates throughout the entire cell volume. Similar to the previous case of bubble growth, the pressure inside the cell is controlled by the vapour pressure in the bubble, and this control is more noticeable as the cell size decreases, meaning the higher the bubble concentration. Consequently, as shown in the figure, at high concentrations, the mean pressure in the liquid substantially exceeds the external pressure until the bubbles completely collapse.

The elevated pressure established within the cell acts to suppress the bubble’s expansion following maximum compression and its subsequent oscillation, which would be characteristic of the collapse of a single, isolated bubble. However, at comparatively low concentrations (curve 1 in Figure 3), the behaviour of the bubbles closely approaches that of a single bubble. During the collapse process, the pressure inside the cell exhibits only weak oscillation, tending toward the external pressure value, p_f=130 kPа. The collapse of bubbles at low concentrations is accompanied by oscillation of their size, as evident from Figure 3 (curve 1), whereas at high concentrations, the bubble radius decreases monotonically until complete collapse. Interestingly, despite the qualitative difference in the compression process, the lifetime of the bubbles in the cluster does not depend on the concentration and is practically identical to the duration of the single bubble collapse under the same conditions. Furthermore, for identical values of k_b, the initial bubble size R₀ has only a minor influence on their behaviour. A study of the growth process at various R₀ values, ranging from 5 µm to 100 µm, demonstrated that shortly after the onset of bubble expansion (the time depending on concentration), the magnitude of the current pressure in the inter-bubble space and the rate of vaporization are practically independent of R₀.

The magnitude of the pressure drop or the degree of liquid superheating exerts a similarly minor influence on the dynamics of bubble growth for identical values of T_f0 and k_b. For example, when the pressure in a monodisperse ensemble is reduced from an equilibrium value of 60 kPа to 20 kPа, the corresponding time dependencies for the averaged pressure p_m = f(τ) essentially coincide. This congruence also applies to the dependencies of R(τ), u_R(τ) and γ(τ). However, during the short initial growth phase, before the quasi-equilibrium pressure is established in the liquid phase, the intensity of bubble expansion is higher when the pressure differential is larger. The duration of the external pressure change, Δτ, also has only a weak effect on the characteristics of the dependencies p_m(τ), R(τ). In the inertial stage, the rate of vapor phase growth is higher when the rate of pressure reduction dp_f/dτ is higher. At the same time, intensive bubble growth facilitates a more rapid approach to the quasi-equilibrium value p_m. Regardless of the duration of the pressure drop, the same averaged pressure p_m is established inside the cluster after a very short period. This pressure is what governs the subsequent evolution of the bubbles, meaning that after only approximately 80 µs following the onset of the bubble ensemble expansion, the recorded parameter values for R, u_R or γ no longer allow one to discern the cluster’s development history or determine whether the change in external pressure has ceased or is still ongoing.

5.2. Velocity and Pressure Fields in the Inter-Bubble Space

The modified model allows for the individualization of each bubble, taking into account its current size, location within the cluster, kinematic and thermodynamic characteristics, and other properties. Information obtained through the cell model approximation is also incorporated into the description of the bubble cloud behaviour. Figure 4 is presented as an illustration, showing the instantaneous distribution of liquid microflow velocities inside the cluster at a specific moment in time. The microflows, which exhibit a complex nature, are caused by the expansion of the bubbles during the cavitation process. The arrows in the figure indicate the direction of the velocity at local points but not its magnitude. Clusters containing between 10 and 1500 bubbles were considered in the analysis. At any given time, the local values of the velocity vector $\stackrel{\rightharpoonup }{u}$ (magnitude and direction), the velocity field potential φ, and the pressure magnitude p are determined at various points in the liquid. To visualize the results, the projections of the velocity vector onto a specified plane passing through the volume occupied by the bubble cloud were calculated. It is characteristic that the velocity field exhibits a truly cellular structure. Directly at the surface of the bubbles, the flow is close to radial. However, as the distance from the surface increases, the microflow pattern becomes more complex, and at the boundaries of bubble interaction, the velocity vectors change direction sharply. This figure demonstrates that bubble growth occurs more intensely at the edges of the cloud, which is attributed to the non-uniform distribution of the averaged pressure p_k throughout the cloud volume: in the central part, p_k>>p_f, whereas at the periphery, p_k→p_f.

Figure 5 illustrates how the magnitudes of the local pressure p and local velocity φ change along the direction I-I (as defined in Figure 4). Significant velocity fluctuations are observed within the cluster along any direction, whereas outside the cluster, even near its boundaries, the velocity magnitude changes monotonically, maintaining a practically constant direction. Despite such abrupt velocity changes in the inter-bubble space, the pressure pulsations are, as evident from Figure 5, relatively weak. Nevertheless, the magnitude of the pressure gradient reaches up to 20 kPа/mm at local points, and the velocity gradients exceed 20 m/s per 1 mm. The figure clearly indicates that the pressure within the cloud volume is higher than at the periphery. This increased internal pressure leads to a situation where the bubble size and, consequently, the void fraction are higher in the peripheral zone of the cloud than in its interior.

The dashed lines in Figure 5 demarcate the region occupied by the cluster. The data obtained indicates that within the inter-bubble space of the cloud, even during the monotonic expansion of the bubbles, there are sharp jumps in pressures and velocities characteristic of turbulent flow. This nature of flow contributes to the intensification and stimulation of heat and mass transfer and hydromechanical processes within the liquid phase of the bubble system. When cavitation conditions are realized in a liquid, a multitude of vapour-gas bubbles grow and then simultaneously collapse, forming a cavitation cluster. The effectiveness of the cavitation impact is determined by the collective influence of all bubbles within the cluster. Cluster sizes typically range from 0.1 to 5 cm. The concentration of nuclei in the cluster depends on the liquid’s temperature, gas content, and the presence of impurities. The superposition of shock waves from each individual bubble forms a single destructive shock wave. However, this conventional concept does not account for the destruction of micro-objects inside the cluster itself. Furthermore, when modelling the behaviour of a cavitation cluster and developing analytical methods for cavitation reactors, temperature is often not even included as one of the model parameters [56,57]. The following section provides an assessment of the main characteristics of heat and mass transfer across the bubble surface.

5.3. Heat and Mass Transfer Across the Interfacial Surface During Bubble Collapse

As demonstrated earlier, the dynamic effects of individual bubbles at the cluster boundary during the collapse stage significantly exceed the average effects of the cluster itself. The most intense heat and mass transfer across the interfacial surface, and the most vigorous mass transfer processes (evaporation and condensation), occur during the stage of maximum bubble compression, as illustrated in Figure 6 and Figure 7. Although the amplitude values of the heat flux impulse JQ+q are extremely high during the maximum compression stage, the quantity of thermal energy transferred through the bubble wall during this time is very small. Due to the extremely short duration of this impulse, this part of the energy constitutes only a negligible fraction of the total thermal energy accumulated in the bubble after the transformation of kinetic energy into potential energy.

The intensity of mass transfer, J, also increases to an extreme degree during the stage of maximum compression, but this fractional change in acquired or lost vapours mass is negligibly small in the overall mass balance. For the remainder of the bubble’s development time, the inherently low rates of heat and mass transfer can still exert a relatively large influence on the bubble’s dynamics. Figure 8a,b illustrate how the vapours mass in the bubble, m_s, and the quantity of heat accumulated in the bubble, q, change during the oscillation process. Despite the strong fluctuations in the bubble size, both quantities, m_s and q, decrease practically monotonically, with the largest reduction in vapours mass and thermal energy occurring during the first period of oscillations. The absolute quantities of energy and mass transferred across the bubble surface per unit time are not very large. However, considering the minute dimensions of the bubble at the moment of maximum compression, the specific fluxes across a unit surface area are substantial. As shown in Figure 8 and Figure 9, the magnitude of the specific mass flux of vapours condensing into the liquid or evaporating into the bubble is estimated to be on the scale of 100 kg/m²s, and the specific energy flux is on the scale of 100 МW/m².

5. Conclusions

Based on the conducted analysis, a mathematical model for bubble cluster dynamics has been developed, which describes the evolution of a collection of vapour bubbles when thermodynamic equilibrium is disrupted in the vapour-liquid system. Computational experiments were performed to investigate the patterns of bubble evolution resulting from a sharp decrease or increase in external pressure. When the external pressure is released in a vapour-liquid bubble medium, the pressure in the liquid phase rapidly increases, asymptotically approaching the saturated vapour pressure corresponding to the liquid temperature. The rate of approach to saturation pressure is proportional to the bubble concentration. he degree of superheat, the rate of pressure release, and the initial size of the bubble nuclei only influence the initial stage of cavitation development.

The computational experiments provided quantitative information regarding the local values of pressure and velocity at any point in the liquid phase of the developing bubble ensemble. It was established that even during the monotonic expansion of bubbles in the thermal regime, the magnitudes of the velocity and pressure gradients reach extremely high values.

The dynamic action of cavitation is attributed to the following factors:

• – during bubble growth and compression, complex non-stationary microflows are formed inside the cluster due to the superposition of radial flows towards each bubble; these flows feature extremely high values of shear rates and shear stresses, the effect of which is present throughout the entire lifetime of the cluster;
• – calculations indicate that in the inter-bubble space, shear rates at local points on the boundaries of flow interaction can reach 10⁶ s^-1, which corresponds to a shear stress of approximately 1 kPа when water is the dispersed phase;
• – during collapse, each bubble in the cluster acts as the center of a spherical shock wave with a pressure amplitude that can reach 100 MPа.

The non-stationary fields of velocities and pressures in the liquid inside the cluster are extremely complex. During bubble compression, the volume-averaged pressure inside the cluster rapidly increases and can exceed 10 MPа in the final stage. The simultaneous collapse of bubbles causes the pressure within the cluster volume to drop sharply, and a high-amplitude resultant impulse is radiated into the surrounding medium.

A significant finding is that even during the stage of monotonic cluster expansion, when the bubble growth rate does not exceed 2 m/s, the shear rates in local liquid zones are quite high—on the order of 5·10⁴ s^-1. Thus, while a high level of thermal impact is only realized during the nanosecond-scale stage of maximum cluster compression, the effective influence of shear stresses is evident throughout the entire existence of the cluster.

The question of the total dynamic effect of the cluster remains a subject for discussion. If one assumes that dynamic effects are determined only by the collapse of the outer peripheral bubbles, theoretical models are likely to predict overestimated dynamic parameter values. This is because real processes exhibit effects of mutual compensation of counter-propagating impulses directed from different bubbles, a finding confirmed by certain experimental studies.

The developed modified mathematical model is designed for the analysis of pressure and temperature fields, as well as microflows, in the inter-bubble space of the ensemble during boiling and cavitation processes.