Theoretical Study of Supercavitation Bubble Formation Based on Gillespie’s Algorithm

Abstract

Understanding the creation and development of a supercavitation bubble is essential for the design of supercavitational underwater vehicles and applications. The pressure field of the supercavitation bubble is one of the most significant factors in these processes, and it should be taken into account in the analysis. The underwater vessel is surrounded by a supercavitation bubble which is, in fact, an inhomogeneous fluid containing cavities (also described as microbubbles). The distribution of the cavities in the supercavitation volume dictates the pressure field and thus determines the stresses and forces that act on the vessel and affect its motion and stability. In this research, we suggest a new approach to studying the bubbles’ formation and learning about the cavities’ distribution in the low-pressure volume that envelops the underwater vehicle. We used Logvinovich’s principle to describe a two-dimensional ring of fluid that is created at the front edge of the supercavitation body and moves downstream along the vessel. To describe the distribution of the cavities we used Gillespie’s algorithm, which is usually used to describe biological and chemical systems. The algorithm succeeded in describing the random movement of the cavities in the cross-section under various conditions and also in describing their distribution and effects on the macroscopic system. A few factors of the physical characteristics of the fluid and the flow conditions were examined (the initial bubble supply, and the rate coefficients of creation and collapse). The results led to the conclusion that with an examination of those factors and using Gillespie’s algorithm, predictions of the distribution and thus the development of supercavitation could be achieved. The main finding of the analysis was that asymmetric development of the bubbles took place, in spite of the symmetry of the physical problem, as observed in high-resolution experiments.

1. Introduction

Cavitation is a unique phenomenon in hydrodynamics, where vapor appears within a homogenous liquid medium. It occurs in different situations commonly characterized by the pressure decreasing below the vapor pressure of the liquid. Consequently, the liquid starts to evaporate locally, generating volumes of vapor, referred to as “cavitation bubbles”. The vapor’s structures are often unstable, and when they reach a zone of increased pressure, they may violently collapse. Under certain conditions of velocity and geometry, a cavitation bubble formed over a moving underwater vehicle may grow and extend, enveloping the entire body. Such situation is called supercavitation and may be used purposefully to reduce friction and drag, allowing movement underwater at extremely high speeds. To prevent a sudden collapse of the supercavitation bubble and ensure its completeness, mechanisms for injecting gas or liquid from the body may be applied. By means of these injection mechanisms, supercavitation bubbles can be created under different conditions than those where bubbles occur naturally. This phenomenon is called “artificial supercavitation” and it is of interest to researchers as the basis for developing underwater vehicles which can implement the phenomenon, as well as research tools allowing future investigations of supercavitation under conditions that are conducive to research [1,2,3].

Understanding the creation and development of a supercavitation bubble is essential to the design of supercavitational underwater vehicles and applications. The pressure field of the supercavitation bubble is one of the most significant factors in these processes, and it should be taken into account for accurate analyses and predictions. This subject has not been studied yet; the common assumption of a constant pressure and a homogenous vapor regime in the supercavitation bubble may be inaccurate and misleading. An underwater vehicle is surrounded by the supercavitation bubble, which is, in fact, an inhomogeneous fluid containing microscopic cavitation bubbles commonly referred as cavities. The distribution of the microbubbles in the supercavitation volume dictates the density and the pressure field, and thus determines the stresses and forces that act on the vessel and affect its motion and stability [2]. In spite of its significant meaning and important role in the creation, development and stability of bubbles, the long-term study of the topic is far from being complete. There are still many unsolved questions regarding the subject, and the description of the pressure field remains unknown [4,5,6,7,8,9,10,11]. One of the reasons for the unfortunate lack of progress in the topic is the extreme experimental challenges. Measuring pressure inside a cavitational flow requires the insertion of disturbances in the flow, which will interrupt the natural creation of cavitation [3,12,13,14,15]. Inevitably, cavitation will appear on the gauges and will not only bring to errors in to the measurements but will also change the pressure field. Another reason for the difficulty of research in this field is the macroscopic approach required for describing the flow outside the bubble. Numerical attempts were proposed to investigate this phenomenon with computational fluid dynamics (CFD) programs by solving the continuity equation that used the equation of state to distinguish between the liquid’s density and the vapor’s density. The disadvantage of this method is that it cannot distinguish between the liquid’s vapor and the noncondensable gases used in ventilated supercavitating vehicles [16]. Recent studies have succeeded in improving the results by modeling the flow with the multiphase Navier–Stokes equations, which separate the flow to three components (the liquid, the vapor and the noncondensable gas) [17,18]. Other recent attempts in the past years have succeeded only in determining the flow field under certain conditions of symmetry, uniformity and steadiness, and with no treatment of the boundaries, difficulties and effects [19,20,21,22,23]. Most numerical studies have presented the flow regime and an approximation of the bubbles’ geometry, overlooking the accurate distribution of pressure and structure of the inner field [23,24,25].

The modeling of cavitational flow is based on fluid mechanics and continuum theory and does not include stochastic descriptions, which are commonly used for describing phase transitions and nonequilibrium thermodynamics. The pressure inside the supercavitation bubble is usually assumed to be constant and equal to the vapor pressure of the liquid [1,2,3]. This inaccuracy does not harm macroscopic modeling, but it has left researchers and engineers with an inability to comprehend the regime inside the bubble. This has caused a problem regarding the gas supply required in artificial supercavitation as well as difficulties in designing vessels which are planned to change their velocity and direction in the water, while their pressure field remains unknown.

In this research, we suggested a new approach in order to study bubble formation and to learn about the cavities’ distribution in the low-pressure volume that envelops the underwater vehicle. We used Logvinovich’s principle to describe a two-dimensional ring of fluid created at the front edge of the supercavitation body and moving downstream along the vessel [26,27]. To describe the distribution of the cavities, we used Gillespie’s algorithm, which is usually used to describe biological and chemical systems [28,29].

2. The Physical Model

We considered supercavitation bubbles developing along an axisymmetric cylindrical object in a uniform flow of water. When the pressure decreases below the equilibrium vapor pressure of the liquid, the water starts to evaporate, developing a supercavitation bubble over the body. As the flow velocity increases, the pressure decreases further downstream, and the bubble grows and can envelop the entire body. Logvinovich’s principle allows us to analyze the supercavitation bubble as a series of rings [26,27]. Every ring represents a cross-section of the bubble and its expansion. We examined a single ring, which described a cross-section of the flow around the cylindrical object that existed at a low pressure. This pressure was a critical pressure that enabled phase transition so cavitation could occur. The ring was divided to microvolumes of fluid. Each volume consisted of fluid particles which could cavitate or diffuse. A particle that cavitated could last but could also collapse.

2.1. Cavitation

Cavitational flows are described by the nondimensional cavitation number $\sigma$ (defined by Equation (1)).

$$\sigma =\frac{p_{\infty }-p_v\left(T_{\infty }\right)}{\frac{1}{2}\rho U_{\infty }{}^2}$$

(1)

where $p_{\infty },U_{\infty }$ are the pressure and velocity, respectively, at some reference point far from the cavitation, $p_v$ is the vapor pressure of the liquid, $\rho$ is the density of the liquid and $T_{\infty }$ is the reference temperature of the liquid. The cavitation number shows what may happen in a given flow when either the overall pressure is decreased or the flow velocity is increased, so that the pressure at some point in the flow approaches the vapor pressure of the liquid at the reference temperature. It indicates the potential of cavitation to occur locally in a flowing liquid.

Nucleation in flowing liquids will occur as the cavitation number is reduced. The particular value of the cavitation number $\sigma$ at which nucleation will first occur is called the incipient cavitation number and is denoted as ${\sigma }_i$. Further reductions in $\sigma$ below ${\sigma }_i$ cause an increase in the number and extent of vapor bubbles.

For cavitation and bubble creation, the deposition energy in a volume must reach a certain value, namely the critical deposition energy of Gibbs [30]

$$W_{CR}=\frac{16\pi S^3}{3{\left(\Delta p_C\right)}^2}$$

(2)

where $S$ is the surface tension of the liquid. For a bubble in equilibrium, the pressure difference between the gas phase and the liquid phase balances the surface tension force (Equation (3))

$$\Delta p_c=\frac{2S}{R_c}$$

(3)

where $R_c$ is the critical radius, which must be the intermolecular distance of ${10}^{-4}\mathsf{\mu }\mathrm{m}$ in the case of cavitation. It is of interest to substitute a typical surface tension of $0.05\mathrm{kg}/{\mathrm{s}}^2$ [3].

The probability that the stochastic nature of the thermal motions of the molecules would lead to a local energy perturbation of magnitude $W_{CR}$ related to typical kinetic energy of the molecules, $k_{B}T$, is described by the Gibbs number (Equation (4)), where $k_B$ is Boltzman’s constant and T is the temperature, which was suggested by Lineard and Karini [31] to be the critical temperature of the liquid.

$$G_b=W_{CR}/k_{B}T$$

(4)

Thus the probability is proportional to $e^{-G_b}$ and the nucleation rate $J$, which is defined as the number of nucleation events occurring in a unit of volume per unit of time, as described in Equation (5).

$$J=J_0{e}^{-G_b}$$

(5)

where $J_0$ is a factor of proportionality. Various functional forms have been suggested for $J_0$. A typical form is that given by Blander and Katz [32], as described in Equation (6).

$$J_0=N{\left(\frac{2S}{\pi m}\right)}^{\frac{1}{2}}$$

(6)

where $N$ is the numeric density of the liquid $[{\mathrm{molecules}/\mathrm{m}}^3]$ and $m$ is the mass of a molecule. Though $J_0$ may be a function of temperature, the effect of an error in $J_0$ is small compared with the effect on the exponent, $G_b$, in Equation (7). After nucleation, an immediate process of mass diffusion starts and the bubble expands. The diffusion process is slow and typical diffusion times are long, in the order of seconds. This was studied by Plesset and Hsieh [33], who gave the following expression for the time required for doubling the radius of the nucleus:

$$\tau =\frac{9}{4}\frac{R_0{}^2{\rho }_g}{{\epsilon }^2{D}_b{C}_{\infty }}$$

(7)

where $\epsilon$ is the relative amplitude of the pressure oscillations, $D_b$ is the diffusion coefficient of the mass and $C_{\infty }$ is the concentration far from the nucleus. This time is longer in comparison with the time needed for the bubble to collapse, which typically takes milliseconds. According to Blander and Katz [32], the nucleation rate is about ${10}^{-6}-{10}^{-8}{\mathsf{\mu }\mathrm{m}}^{-3}{\mathrm{s}}^{-1}$ and its value changes by about three or four orders of magnitude per degree Celsius. To describe the production of the cavities (microbubbles) in the cross-section, we will use several rates.

2.2. Diffusion

The microbubbles in the cross-section at temperatures above absolute zero have a thermal motion and can diffuse in the water. To calculate the diffusion coefficient for gas in liquid, we used the Arnold equation (Equation (8)) [34], where the empirical constant 0.01 is accurate at 20 °C

$$D=\frac{0.01{\left(\frac{1}{M_g}+\frac{1}{M_l}\right)}^{1/2}}{{\mu }_l{}^{1/2}{\left(V_g{}^{1/3}+V_l{}^{1/3}\right)}^2}$$

(8)

where $M_g,M_l$ are the molecular weight of the gas and of the liquid, respectively; $V_g,V_l$ are the molar volume of the gas and liquid, respectively; and ${\mu }_l$ is the absolute viscosity of the liquid.

The correct values of the diffusion coefficient for cavitation bubbles in water at a normal temperature can vary from $0.5\times {10}^{-5}{\mathrm{cm}}^2/\mathrm{s}$ to $2.5\times {10}^{-5}{\mathrm{cm}}^2/\mathrm{s}$ [1,34,35].

The cavities can diffuse in the cross-section radially from the cavitator toward the supercavitation bubble’s surface (the edge of the cross-section); this is known as expansion, and can also diffuse azimuthally, circling the cavitator in both directions (clockwise and counterclockwise).

To describe the diffusion of the microbubbles in the cross-section, we used a diffusion coefficient of $2\times {10}^3{\mathsf{\mu }\mathrm{m}}^2/\mathrm{s}$ for both directions (radial and azimuthal).

2.3. Collapse

Cavitation bubbles are usually unstable and are subjected to fluctuations in pressure that can interrupt the equilibrium (as described in Equation (3)) in which the bubble exists. Rayleigh time (see Equation (9)) [1] is the characteristic time of the bubble’s collapse in a homogenous liquid. This characteristic time was used to evaluate the rate of collapse of cavitation bubbles. It should be noted that, in the derivation of the Rayleigh time, surface tension is ignored; if were not so, the collapse would be slightly accelerated.

$$\tau =\sqrt{\frac{3}{2}\frac{\rho }{p_{\infty }-p_v}}{\displaystyle \underset{0}{\overset{R_0}{\int }}\frac{dR}{\sqrt{\frac{R_0{}^3}{R^3}-1}}\cong 0.915R_0\sqrt{\frac{\rho }{p_{\infty }-p_v}}}$$

(9)

The factor 0.915 is the approximate value of $\sqrt{\frac{\pi }{6}}\frac{\Gamma (5/6)}{\Gamma (4/3)}$, where $\Gamma$ is the factorial gamma function. This characteristic time was derived from the Rayleigh–Plesset equation and allowed us to determine the temporal evolution of the bubbles’ radius and, consequently, the pressure field in the liquid. In experimental studies, the characteristic time of bubble creation is about ${10}^2$ times larger. Hence, the average rate of bubble collapse is ${10}^2$ times larger than the creation rate, which consists of nucleation and expansion [1,2,3]. The collapse rate does not depend directly on the bubble’s size but only on the pressure difference, which, of course, depends on the radius of the bubble. The characteristic rate of microbubbles’ collapse (degradation) could be of the order from ${10}^{-8}{\mathrm{s}}^{-1}$ to ${10}^{-3}{\mathrm{s}}^{-1}$. This depends on the temperature, density, pressure difference and surface tension.

These three processes (cavitation/creation, diffusion and collapse) determine the distribution of concentration in the bubble together with the pressure field inside the bubble. They determine the bubble’s formation and stability and affect its dynamics. To simulate the change in the cross-section and the creation and distribution of microbubbles, we used Gillespie’s algorithm.

3. The Mathematical Model and the Description of the Algorithm

To describe the location of the processes occurring in the cross-section of the flow around the cylindrical object, we referred to the cross-section as a series of rings, divided into cells. Each cell can be referred to as a site, in which all the three physical processes mentioned in Section 2 can occur. To model the sites, we divided the cross-section into $h$ rings. Every ring consisted of $K$ cells. The thickness of the ring was $H$ and thus this was the length of each cell. The width of a cell depends on the ring within which it is located. The volume of the cell was calculated by multiplying the length, thickness and the average width, as described in Equation (10).

$$V_i=H\cdot H\cdot (H+H\cdot (1+\frac{\left(i+1\right)\pi }{K}))/2$$

(10)

where $i$ is the index of the ring. The area of contact between adjacent cells was used to describe the diffusion of the microbubbles between the cells of the cross-section. If the diffusion was between cells in the same ring, the contact area was $H\cdot H$, and if the diffusion was between cells from different rings, the area between the adjacent cells was calculated by multiplying the thickness by the average width of the cell from which the diffusion was directed, as described in Equation (11).

$$A_i=H\cdot H\cdot (2+\frac{i\pi }{K}+\frac{(i+1)\pi }{K})/2$$

(11)

The cells’ dimensions are described in Figure 1.

In order to simulate the cavitation, collapse, and diffusion of the microbubbles in the cross-section, we used Gillespie’s algorithm [28,29]. A scheme of the algorithm is presented in Figure 2. We used four random numbers. The first random number was the time of occurrence, which indicated the time at which an event occurred. The second random number indicated what the event was, and, if the event was diffusion, what the direction was. In the case of diffusion, the direction could be azimuthal or radial. If the diffusion was in the azimuthal direction of the cross-section, the third random number indicated whether the diffusion was clockwise or counterclockwise; and if the diffusion was in the radial direction, the fourth random number indicated whether the diffusion was toward the cavitator at the center of the cross-section or in the opposite direction toward the edge of the supercavitation bubble.

In order to calculate the probability of each event, an array of all the possibilities was used so that each cell was checked with regard to the second random number in order to find the event that occurred. Alpha, namely the array of possibilities, was built from $4\cdot K\cdot h$ cells. Each of the four possible events (cavitation, collapse, azimuthal diffusion and radial diffusion) could occur in every cell of the $K$ cells in each ring of the $h$ rings. Each cell represented the location of the cell, so that the first ring of cells were located at the first $4\cdot K$ cells, the second ring of cells were located at Cells $4K+1$ to $8K$, and so on. For Cell No. 1 in Ring No. 1, the possible events (for example) were at the following array locations: $1,K+1,2K+1,3K+1$. For Cell No. 2 in Ring No. 4, the possible occasions (for example) were at the following locations: $12K+2,13K+2,14K+2,15K+2$. According to the location in the array, an identification of the event was made, and also of the cell in which the event occurred. All events depended on the cell’s location, constants that characterized the process (k1, k2, D) and the dimensions of the cell.

3.1. Cell Location

The cross-section was limited to the inner radius and outer radius. The inner radius represents the cavitator’s radius. In the first ring, which was adjacent to the cavitator, radial diffusion toward the center of the cross-section (toward the cavitator) was not possible. Moreover, in the outer ring (the last ring), radial diffusion directed outside the cross-section was not possible. The supercavitation bubble’s geometry including the outer radius of the cross-section was determined by the flow regime outside the bubble. Thus, the outer radius was limited.

At the first ring next to the cavitator’s edge, the rapid production of microbubbles can occur due to the injection of gas from the supercavitation body in the case of artificial supercavitation. This was reflected as an initial condition of the number of bubbles at the beginning of the simulation. This was examined and is discussed in the results. The values that were examined are shown in

Table 1. Initial conditions and supply of bubbles.

Initial conditions (No. of microbubbles)

0

1

5

10

50

100

.

3.2. Characteristics of the Process

Every process had a rate which affected the probability of its occurrence. In every cell of the cross-section, cavitation occurred uniformly at a constant rate $k_2$, and microbubbles are produced in the system. The microbubbles collapse at the rate $k_1$ and diffuse with the radial diffusion coefficient $D_r$ and the azimuthal diffusion coefficient $D_{\theta }$. The constants characterizing the processes of cavitation, collapse and diffusion were chosen according to the physical model described in Section 2 and are shown in

Table 2. Characteristic constants of the processes.

Rate of cavitation $k_1$	${10}^{-6}{\mathsf{\mu }\mathrm{m}}^{-3}{\mathrm{s}}^{-1}$	${10}^{-8}{\mathsf{\mu }\mathrm{m}}^{-3}{\mathrm{s}}^{-1}$
Rate of collapse $k_2$	${10}^{-3}{\mathrm{s}}^{-1}$	${10}^{-4}{\mathrm{s}}^{-1}$	${10}^{-8}{\mathrm{s}}^{-1}$
Radial diffusion coefficient $D_r$	$2\cdot {10}^3{\mathsf{\mu }\mathrm{m}}^2{\mathrm{s}}^{-1}$
Azimuthal diffusion coefficient $D_{\theta }$	$2\cdot {10}^3{\mathsf{\mu }\mathrm{m}}^2{\mathrm{s}}^{-1}$

.

3.3. Dimensions of the Cell

The dimensions of the cells were different from one another due to the division of the cross-section described at the beginning of this section. The cell volume affects the production of microbubbles (cavitation), and the area of contact between adjacent cells affects the process of diffusion. To simplify the problem, the volume of each cell was calculated as described in Equation (12), in which the width of the cell was taken as the average between its upper arc and its lower arc

$$V_i=H\cdot H\cdot \frac{1}{2}\left(H\left(1+\frac{2\pi \left(i-1\right)}{K}\right)+H\left(1+\frac{2\pi i}{K}\right)\right)$$

(12)

where $V_i$ is the volume of a cell in Ring $i$. The area was calculated in the same way for the contact area between adjacent cells in different rings, as described in Equation (13), and the contact area between adjacent cells in the same ring was constant and equalled $H\cdot H$.

$$A_i=H\cdot \frac{1}{2}\left(H\left(1+\frac{2\pi \left(i-1\right)}{K}\right)+H\left(1+\frac{2\pi i}{K}\right)\right)$$

(13)

The error for not determining the direction of the diffusion from the cell and not taking the appropriate arc in the calculation was of the order of $O\left(\frac{H}{HK}\right)$. Since $K\gg 1$, the error was negligible, as was the error of calculating the cell as a quadrangle and not as a sector of a circle that is of the same order.

These three elements that affect the processes were taken into account in the description of alpha for each event.

The summation of the microbubbles in each cell was carried out by the matrix Bnum. The matrix rows represent the different rings of the cross-section, and the matrix columns represent the cells in every ring.

In each time step of the simulation, only one event can occur, and in every cell, there could be an unlimited number of bubbles.

4. Simulation and Results

To analyze the development of supercavitation, we chose to examine the system from a number of aspects. We ran several simulations. Each simulation focused on another condition of the system. The results were described in plots of the number of microbubbles in every cell in the cross-section, shown as circles, for which the sizes were proportional to the number of microbubbles in the cells’ volumes. At the center of the cross-section, the cavitator was displayed as a black full circle. The simulations were applied to examine the results, the flow conditions, the physical model and also the mathematical model. They also examined the initial conditions and compared some of them in order to study the effect of every aspect of the development of supercavitation.

4.1. Simulation 1—Display of the General Results

Before studying the system and examining the processes, the first simulation was chosen to consider how to display the results, in order to make sure that their representation was optimal and did not diverge from the conclusions. Simulation 1 took 10 s and focused on displaying the results in the supercavitation bubble’s cross-section. To describe the distribution of the microbubbles in the cross-section, we used several ways. All descriptions were based on the coefficients and conditions shown in

Table 3. Simulation 1: Constants and conditions.

$k_1$	${10}^{-8}\mathsf{\mu }{\mathrm{m}}^{-3}{\mathrm{s}}^{-1}$
$k_2$	${10}^{-6}{\mathrm{s}}^{-1}$
$D_r$	$2\cdot {10}^3\mathsf{\mu }{\mathrm{m}}^2{\mathrm{s}}^{-1}$
$D_{\theta }$	$2\cdot {10}^3\mathsf{\mu }{\mathrm{m}}^2{\mathrm{s}}^{-1}$
Initial conditions (No. of microbubbles)	$0$

. The first description displays the number of microbubbles in every cell (see Figure 3a). In this description, it is easy to notice the change in the number of microbubbles from one ring to another; on the other hand, it is hard to learn about the distribution in the cross-section. If one wishes to study the expansion and the spreading of the gas in the supercavitation bubble, it would be better to use a description in which the number of microbubbles related to the cell volume is displayed and not the number of cavities, since the volume of the cells changes from one ring to another (see Figure 3b) Here, the number of microbubbles was normalized with the ratio between the cell’s volume and the cells’ volume in the first ring, so the first rings in both figures were identical and the changes would be seen when the cell volumes changed, as would the spreading of the microbubbles in the cells. Otherwise, the display of the number of the microbubbles could be misleading in the sense of their concentration in the cross-section, etc. Contrarily, if the goal were to study the processes in the cross-section, the actual number of microbubbles would be required in order to examine the change and the growth. In these circumstances, the first display will be more suitable. In both displays, we can see the nonuniform spread of the microbubbles. The cross-section was inhomogeneous and there was not even inner symmetry in the microbubbles’ distribution.

Another way to present the bubbles’ distribution in the cross-section was to neglect the change in the cell’s area and volume from one ring to another. This would be legitimate according to Equations (12) and (13) only with a very large $K$ and a small number of rings, which represent a long and narrow supercavitation bubble at a high velocity. In this case, there was no difference between the two types of display mentioned above. Figure 4 shows the results for identical cells under the same conditions in Table 3. As it can be seen from the figure, the nonuniformity in the microbubbles’ distribution did not change whether the cells’ dimensions were identical or not. For conditions in which this presentation would be more suitable, one would probably talk about fluid particles and not microbubbles. Thus, the constants and coefficients will be different than those presented in Table 1 and Table 2. Moreover, the size of the cells would be smaller. To ease the calculations and analyses, we chose this description of the cross-section with identical cells in our examinations in some of the next subsections.

4.2. Simulation 2—Changes in Time

One of the most significant issues of the development of a supercavitation bubble is the question whether there is a time at which the system reaches to a steady state and whether the nonuniform distribution of the microbubbles in the cross-section mentioned in the previous subsection disappears with time. Simulation 2 was run for different times: 5 s, 10 s, 30 s, 1 min, 1 min and 30 s, and 2 min (see Figure 5). In this simulation, we used the same conditions as in Table 3.

To analyze the results, we added two counters to the simulation: the bubble creation counter and the bubble collapse counter. Figure 6 presents the difference between the two over time. The results confirmed that there was no convergence in time. The maximum of the difference between the production of bubbles and those which degraded could take place a long time after the start of their development. For example, in case of a 60 s simulation, the maximum occurred at approximately 39 s. These results were caused due to the nature of the processes and the nonlinearity in the system. However, there was a constant increase in the production of bubbles over time with the current characteristic coefficient. The increased production was much larger than the increase in degradation, and so there was an increase in the addition of bubbles in the cross-section over time (see Figure 5 and Figure 6).

4.3. Simulation 3—Changes in the Constants and Coefficients

The rate of production and degradation of microbubbles affects the propensities and the probability of events occurring. This affects the number of microbubbles in the cross-section and their distribution. To study how the change in the coefficient affected the number of microbubbles in the cross-section, we compared the number of bubbles created with those collapsed. Simulation 3 used the constants from Table 2 and ran for 30 s without extra bubble supply and an initial condition of zero. To ease the calculations, we used identical cells. The first simulation was for a production rate coefficient of ${10}^{-8}{\mathsf{\mu }\mathrm{m}}^{-3}/\mathrm{s}$ and three different degradation rates, and the second simulation was for a production rate coefficient of ${10}^{-6}{\mathsf{\mu }\mathrm{m}}^{-3}/\mathrm{s}$ and three different degradation rates. The reason for the two simulations was the scale needed for graphing the results. There was a significant difference in the number of bubbles produced for the two rates of cavitation (see Figure 7). As for the change in the degradation rates, we could clearly see the increase in the collapse of bubbles as the value of the degradation rate increased. From Figure 7b,d, it appears that only one order of growth in the degradation rate could lead to a number of collapsed bubbles that was seven times larger.

In both simulations, the number of microbubbles increased over time. The number of cavitation bubbles was much higher than that of collapsed bubbles, and so the total number increased significantly. The number of cavitation bubbles increased together with the number of collapsed bubbles. The difference in the increase was due to the rate coefficients. The yellow curves in both simulations in the phase figures (Figure 8b,c) also showed an increase in the number of collapsed bubbles but, as compared with the results of different rate coefficients, it seemed to tend to zero (this was only due to the graph’s proportions).

4.4. Simulation 4—Changes in the Initial Conditions

Simulation 4 studied the change in the initial conditions of the system and how these affected the number of microbubbles in the cross-section. We compared the number of bubbles created with those collapsed. Simulation 4 used the coefficients presented in

Table 4. Simulation 4: constants and coefficients.

$k_1$	${10}^{-8}{\mathsf{\mu }\mathrm{m}}^{-3}{\mathrm{s}}^{-1}$	$D_r$	$2\cdot {10}^3{\mathsf{\mu }\mathrm{m}}^2{\mathrm{s}}^{-1}$
$k_2$	${10}^{-3}{\mathrm{s}}^{-1}$	$D_{\theta }$	$2\cdot {10}^3{\mathsf{\mu }\mathrm{m}}^2{\mathrm{s}}^{-1}$

. The simulation ran for 30 s with changes in the initial conditions presented in Table 1. To ease the calculations, we used identical cells.

Figure 8 presents the number of bubbles in the cross-section over time. The cavitation bubbles produced in the cross-section and the collapsed bubbles which degraded are both presented in Figure 9. Due to the values of the rate coefficients, the number of cavitation bubbles was much bigger and it increased significantly, as opposed to those that collapsed. Regarding the initial conditions, it is better to look at Figure 8, where the change in the production of bubbles versus bubble degradation is presented.

According to Figure 9, the initial number of bubbles in the first ring did not change the process and did not affect the number of bubbles produced in the cross-section as well as the number of bubbles that collapsed in the cross-section. The only thing that changed was the total number of bubbles in the cross-section. In addition, we cannot say anything about the change in the curves’ inclination due to the change in the initial conditions. The number of cavitations and collapsed bubbles increased together without any apparent influence from the initial number of bubbles in the system.

In Figure 10 (where we use a normalized display for different cell sizes), it can be seen that the initial condition of the number of bubbles in the first ring did not affect the bubbles’ distribution in the 30 s simulation of their development. However, the fact that at the beginning of the movement, there was a large number of bubbles close to the cavitator enveloping the body with the gas medium encouraged a reduction in drag and improved the supercavitation vessel’s movement and stability. We can see that although there was no constant gas supply but only an initial supply of bubbles close to the cavitator, the number of bubbles in the first ring did not become smaller over time and the gas medium along the cavitator was maintained in the cross-section when this moved along the body with the flow stream. This describes a situation in which there is a microbubble (or gas) supply from the edge of the cavitator only at a constant rate, which suits the rate at which the cross-sections at the front edge are exchanged and an individual cross-section of the supercavitation bubble moves along the body.

4.5. Error Approximations and Implications

The simulations of stochastic processes can create two kinds of errors: errors that occur from the random numbers in each simulation, and errors resulting from the geometric division of the cross-section.

The first kind of error may relate to a change in the values in a local cell in a certain run resulting from the nature of the stochastic process, which used random numbers. Nevertheless, observations of the overall result of a large number of runs (>100) showed that the trend and rate of growth were found to converge to within +/−1%. Simulation 2 was run 150 times and the difference in the maximum and minimum number of bubbles was up to 10% in a specific time. However, the resulting overall rate of growth of the supercavitation bubble implied an uncertainty of less than 1% (0.88%). Simulation 3 was run 100 times and the outcomes showed practically the same rate of growth in the number of microbubbles with time (within less than a 0.7% difference). Simulation 4 was also run 100 times, presenting the same trend over time, even if a difference in the number of bubbles was found locally between various runs of the simulation. One can conclude that even if there are slight local deviations from run to run, the effect on the overall result is minimal regarding the rate of growth and the indication of asymmetrical trends.

The second kind of errors can be caused by the choice of the number of cells in the cross-section $K$, which determines the accuracy of the simulation and its results. For a larger number of cells in a ring, the size of the cells will be smaller by the same proportion according to Equation (14):

$$HK=\pi d$$

(14)

The angle of every cell can change with a change in $K$ according to Equation (15)

$$\Delta \alpha =\Delta \left(\frac{2\pi }{K}\right)\sim \frac{2\pi }{K^2}\Delta K$$

(15)

indicating an order of $O\left(\frac{1}{K^2}\right)$. As the number of cells was large, this error was very small. Even the shift in the angle around the cross-section would be negligible, reaching no more than the order of $O\left(\frac{1}{K}\right)$, derived by multiplying the change in the angle by the number of cells). The change in the cells’ area would be of the order of $O\left(\frac{1}{K^3}\right)$ and the volume would be of the order of $O\left(\frac{1}{K^4}\right)$. The change in the distance between neighboring cells would be of the order of $O\left(\frac{1}{K^2}\right)$, as would the change in the distance between the closest cells in different rings and different angles, which has no area of contact. As mentioned in Section 3.3, the effect on diffusion would be also of the order of the change in the angle. As a result, the effect of the number of cells on the symmetry is negligible; hence the errors caused by the choice of $K$ could not affect the solutions and lead to a different trend in the distribution of microbubbles in the cross-section.

Finally, the simulations presented clear trends and an obvious lack of symmetry, which were be affected by the method’s errors or calculations.

5. Discussion

This discussion aimed to present how the novel approach used in the present analysis can contribute to better predicting the trends observed in experimental research.

Previous studies which examined axisymmetric supercavitational flows and presented the flow regime and approximations of the bubbles’ shape and geometry, did not account for the accurate distribution of pressure and the inner structure of the supercavitation bubble, and hence resulted in a symmetric supercavitation bubble [5,6,8]. Moreover, past numerical studies [17,18,19] applying multiphase CFD simulations have shown symmetric results for supercavitational flows when the angle of attack of the body was zero.

In contrast to theoretical studies, experimental studies have shown deviations from symmetrical behavior. Brennen [1] summarized some representative examples of experimental observations of supercavitational flows, including large-scale cavitation, sheet cavitation, attached cavitation, cloud and vortex cavitation, and more. Observations of bubble cavitation and a bubbly wake prior to the transition to the fully developed cavity were presented by Brennen [36] as well. Examples illustrating the attached tails formed behind a traveling cavitation bubble were given by Ceccio and Brennen [37], as well as by Brenen [1]. All those observations have shown the asymmetrical pattern of the cavitational structure. The authors of the present study revealed such behavior as well in an experimental investigation of supercavitation bubbles in a free-surface flow [15]. It should be emphasized that these asymmetrical features resulted from systems that, in most cases, exhibited clear symmetry (the supercavitational bodies, the cavitators, the motion and the flow field), and that the asymmetric shape and structure of supercavitational bubbles were common experimental findings.

Furthermore, even in cases where macroscopic observations of supercavitation flows via video recording at 30 pictures per second (pps) showed an apparent symmetry, a recording with a higher resolution (at 50,000 pps) revealed asymmetric behavior (see Ahn et al. [12,13]). Similarly, in a previous investigation of supercavitation over cylindrical bodies in a duct flow [24,25], these authors found, both numerically and experimentally, a symmetrical bubble shape. However, the inner structure deviated from symmetry.

The new analytical approach suggested and presented in the present study, which took the inner field of the supercavitation bubble into account and simulated the stochastic processes of the microbubbles’ distribution, and thus the development of the supercavitation bubble, and predicted an asymmetric and uneven structure. Hence, this study contributes to a better evaluation of the geometry of the bubble, revealing the trends observed in empirical investigations.

6. Conclusions

The problem of the development of a supercavitation bubble consisting of microbubbles, which are created, diffuse, and could also collapse, proved to be a complicated problem which is not easy to solve, and the process is hard to predict. The creation, diffusion, and collapse, which are random processes, cannot lead to a uniform distribution in the supercavitation’s cross-section despite the system’s symmetry. This explains the instability of the bubble and the challenges in the navigation of vehicles using such technology. With time, the number of microbubbles in the cross-section increases and the medium fills with gas. In this case, the supercavitation becomes more stable and easier to control. The initial number of microbubbles produced in the system (describing cases of gas supply in artificial supercavitation) affects the total number of bubbles in the system and also the bubbles in the first ring during the whole process. This can help reduce drag and is a very useful tool. For supercavitational vehicles, these results can be helpful for evaluating the injection requirements (rates and duration). The most significant influences on the bubbles’ distribution in the cross-section are the rate coefficients of cavitation and collapse. These are factors of the physical characteristics of the fluid and the flow conditions. Changing the rates will change the distribution, and only a comprehensive examination of these will enable to predict the distribution of bubbles and thus the development of supercavitation. Moreover, by taking the inner structure of the bubble into account, an accurate analysis of the stability could be carried out, which would be able to improve the design of supercavitation vessels and applications. Previous experimental studies [1,12,14,15,36,37] revealed asymmetric images of the inner ring of supercavitation bubbles that can also be explained by the results of the present research. The previous numerical studies, which have presented the flow regime and an approximation of the bubbles’ geometry, by overlooking the accurate pressure distribution, could make use of the new approach of this study and present a better evaluation of the inner field of the bubble, and thus the flow field and the geometry of the supercavitation bubble at the macroscopic scale. Future work that combines the previous numerical methods describing the flow regime, together with evaluating the pressure distribution of the inner field of the bubble using the method presented in the present study, would be valuable and effective, and would contribute to higher accuracy in the descriptions of the shape, formation and stability of a supercavitation bubble.