The Effects of Viscosity and Capillarity on Nonequilibrium Distribution of Gas Bubbles in Swelling Liquid–Gas Solution

: A detailed statistical description of the evolution of supersaturated-by-gas solution at degassing has been presented on the basis of ﬁnding the time-dependent distribution in radii of overcritical gas bubbles. The inﬂuence of solution viscosity and capillarity via internal pressure in the bubbles on this distribution has been considered until the moment when the gas supersaturation drops due to depletion and stops nucleation of new overcritical gas bubbles. This study is based on our previous results for the nonstationary growth rates of overcritical bubbles depending on gas supersaturation, diffusivity and solubility in solution, solution viscosity, and surface tension on bubble surface. Other important factors are linked with the initial rate of homogeneous gas bubble nucleation and coupling between diffusivity and viscosity in the solution. Here, we numerically studied how all these factors affect the time-dependent distribution function of overcritical bubbles in their radii, maximal and mean bubble radii, and the time-dependent swelling ratio of a supersaturated-by-gas solution in a wide range of solution viscosities.


Introduction
The theory of swelling of a gas-supersaturated solution as a result of nucleation and further growth of overcritical gas bubbles is in demand in many problems of fundamental science and related technologies. Such a theory is needed for describing bubble nucleation and subsequent foaming as a general phenomenon [1][2][3], the processes of intensive foaming in food and beverages [4][5][6], embolism in tissues and blood vessels [7][8][9], the formation of porous materials from polymeric melts [10][11][12], and volcanic eruptions [13][14][15][16]. The foundations of this theory were previously formulated in Refs. [17][18][19][20] and extended to the multicomponent case in Refs. [21,22]. In these works, the importance of the high solubility of the gas and its strong supersaturation in the initial solution was underlined as the conditions for the implementation of the mode of self-similar nonstationary diffusion growth of overcritical bubbles [17][18][19][23][24][25] with a significant swelling of the entire solution. It was also shown there that, in the case of stationary diffusion growth of gas bubbles, the swelling of the solution is insignificant and stays within a few percent. Further studies [26][27][28][29] showed that a marked role in the growth of gas bubbles can be played by capillary and viscous effects, which, through the pressure inside a bubble, slow down the transition to a self-similar regime and may even prevent it at all. More detailed analytical study of the limiting growth regimes and a numerical calculation of the growth rate of overcritical bubbles in the presence of capillary and viscous effects in a wide range of solution viscosities developed in Ref. [29] demonstrated, that even in the absence of the self-similar diffusion regime of bubble growth, the nonstationary gas concentration decline around the overcritical bubbles can be significant and provide nonstationary growth rates of bubbles [30]. It has also been mentioned [30] that coupling between diffusivity and viscosity can enhance the damping effect on swelling with increasing solution viscosity. As a result of previous studies, we expect that the capillary and viscous effects may noticeably affect the general state of the supersaturated-by-gas solution and, in particular, its swelling at degassing. The first task of this paper is applying the analysis [29] to direct calculation of the distribution function of supercritical gas bubbles in a degassing solution at an arbitrary moment of the nucleation stage in the wide range of solution viscosities and nucleation rates with full allowance for the Laplace pressure in small overcritical bubbles. The novelty of this task lies in the fact that earlier, such calculations of the strong influence of the viscosity on the distribution function of the nucleated and growing gas bubbles with nonstationary diffusion were not carried out at all because of the lack of an analytical basis. Cooperation of the diffusivity and viscosity will be checked for the first time. The second task of this work is to study numerically the collective behavior of the ensemble of gas bubbles in degassing solution through the mean bubble radius and coefficient of swelling of the solution as functions of time and the parameters of the solution. Earlier, the collective characteristics of an ensemble of gas bubbles were considered only for the case of stationary diffusion in the absence of viscosity.

Concentration of Gas around a Single Growing Bubble
At nonstationary growth of a spherical gas bubble in a supersaturated-by-gas liquid solution, a concentric shell is formed around the growing bubble with a radially inhomogeneous distribution of the bulk gas concentration c(r, t), where r is the distance from the bubble center and t is the time. The concentration c(r, t) is described as a solution of the equation of nonstationary gas diffusion to the bubble with moving boundary, where D is the gas diffusivity in the liquid solution, R is the bubble radius, and . R(t) is the growth rate of the bubble radius. The boundary conditions for Equation (1) are where c 0 is bulk concentration of the dissolved gas in the liquid solution, c(R) is the equilibrium concentration of the dissolved gas at the bubble boundary, c ∞ = c(R)| R→∞ is the equilibrium concentration of the dissolved gas in liquid solution at the flat boundary with the gas phase, R * ≡ 2σ/P l , σ is the bubble surface tension, P l is the bulk pressure in the liquid solution, and η is the viscosity of the liquid solution. Here, we have used Henry's law of solubility at the bubble surface with c(R) = sn(R) (s is the gas solubility in the liquid solution) and uniform concentration of gas in the bubble n(R) = P g k B T = P l k B T 1 + 2σ P l R + 4η . R P l R (P g is the gas pressure in the bubble, k B is the Boltzmann constant, and T is the absolute temperature of solution, see details in Ref. [29]). In the case when contributions R * /R and 2η . RR * /σR cannot be neglected in Equation (3) (and this is just our case of interest with important capillary and viscous effects), there is an approximate solution of Equation (1) in the form [29,30]  where ρ = r/R(t), c ∞ = sP l /k B T, and ζ 0 ≡ (c 0 − c ∞ )/c ∞ is the bulk gas supersaturation in the liquid solution and R c = R * /ζ 0 is the critical radius of the bubble which is in unstable equilibrium with liquid solution at the bulk gas concentration c 0 , Function h(R) in Equation (5) is a correction function providing the balance of the number of gas molecules that have left the liquid solution and entered the growing bubble and is determined [29,30] as while the growth rate . R satisfies equation Thus, to determine h(R) and . R(R), it is necessary to solve jointly the system of Equations (6) and (7), then use Equation (5) to find the function Φ as a function of ρ and R, and then substitute found function Φ into Equation (4) to calculate the gas concentration c(ρ, R) in the solution around the single growing gas bubble.

Excluded Volume for a Single Growing Overcritical Bubble
Since the gas concentration c(r, t) in the solution near the bubble, according to solution (4), drops to the equilibrium value, the local supersaturation ζ(r, t) = (c(r, t) − c ∞ )/c ∞ of the gas also drops there. It is convenient to determine the local relative decline ϕ(r, t) of gas supersaturation from the bulk value ζ 0 by the relation Taking into account Equation (4), we can rewrite Equation (8) as In the shell around the growing overcritical bubble, where the decline ϕ(ρ, R) is maximal, the local nucleation rate I(ρ, R) of new overcritical bubbles is suppressed. At the same time, outside this shell, the nucleation rate increases to the initial level I 0 . According to Refs. [17][18][19][20]29], it is convenient to introduce the volume V ex (R), which is excluded for nucleation around the selected growing overcritical bubble with radius R(t) by the integral relation: According to Equation (10), the total number of bubbles formed per unit time at current profile of gas supersaturation in the solution volume V around a selected bubble is equal to the number of bubbles nucleated in the volume V − V ex (R) where the nucleation rate is assumed to be equal to I 0 at bulk gas supersaturation ζ 0 . The ratio q(R) ≡ V ex (R)/V R Colloids Interfaces 2023, 7, 39 4 of 14 of volume V ex (R) and volume V R = (4π/3)R 3 of the bubble with radius R becomes an important parameter. Figure 1 illustrates the physical meaning of the quantities V ex and q.
( ) 3 2 ex 0 1 ( , ) 4 1 According to Equation (10), the total number of bubbles formed per unit time at current profile of gas supersaturation in the solution volume V around a selected bubble is equal to the number of bubbles nucleated in the volume ex ( ) V V R − where the nucleation rate is assumed to be equal to 0 I at bulk gas supersaturation 0 ζ . The ratio of the bubble with radius R becomes an important parameter. Figure 1 illustrates the physical meaning of the quantities ex V and q. As is generally known [17][18][19][20], the ratio  As is generally known [17][18][19][20], the ratio I(ρ, R)/I 0 can be expressed through the decline ϕ(ρ, R) as I(ρ, R)/I 0 = e −Γϕ(ρ,R) , where Γ ≡ −ζ 0 (∂∆F c /∂ζ 0 ) and ∆F c is the activation barrier for formation of the critical gas bubble. With the help of Equation (10), we can represent ratio q as As we will see in the next section, this parameter plays an important role in finding the total excluded volume for an ensemble of growing bubbles.

Total Excluded Volume and Distribution of Overcritical Gas Bubbles in Radii
We assume that the dissolved gas at degassing only redistributes between the volume of the solution and the nucleated and growing bubbles. Let us now designate V tot ex (t) as the sum of individual excluded volumes around all bubbles at time t. If the initial volume of solution equals V L , then V L − V tot ex (t) is the solution volume, where initial nucleation rate I 0 remains preserved to time t.
where the Dirac delta-function δ(R − R(t − τ)) takes into account the spread in radii of the bubbles past the critical size at previous moments of time. Evidently, the total volume V g (t) of bubbles at time t is related to the distribution function f (R, t) as where R 0 is the radius of the smallest overcritical bubble. Because Equation (11) for the excluded volume ratio q(R) ≡ V ex (R)/V R can be applied to any overcritical bubble, we assume [29] that rate dV tot (which is a direct change in time of the total excluded volume) and, in view of the possibility of an accidental overlapping of diffuse layers of bubbles at any point in a liquid solution, is also proportional to the fraction 1 − V tot ex (t)/V L of initial liquid solution volume which is available for the nucleation of new bubbles. Thus, we may write As shown in Ref. [29], with the help of Equation (12), Equation (14) can be converted to equation For the nucleation stage, the first iteration solution of Equation (15) may be written as Substitution of Equation (16) into Equation (12) gives where θ is the Heaviside step function and t(R) is the time of overcritical bubble growth from the initial radius R 0 to radius R (the inverse function for R(t)), R(0) = R 0 . With the help of Equation (17) we can find the mean bubble radius R(t): As follows from Equation (16), the duration t 1 of the nucleation stage is determined by the condition where R(t 1 ) is the maximal bubble radius to the end of the nucleation stage.

Swelling of Solution
One of the important collective characteristics of the state of the degassing solution is its swelling over time. Assuming that the dissolved gas does not leave the solution through its external boundaries, we can determine the swelling coefficient κ as Substitution of Equations (13) and (17) into Equation (20) gives the following expression for swelling coefficient κ as a function of time t:

Results and Discussion
According to the tasks formulated in the Introduction, we are interested in direct calculation of the distribution function f (R, t) of supercritical gas bubbles in a degassing solution and finding the mean bubble radius R(t) and coefficient κ(t) of swelling of the solution at an arbitrary time moment of the nucleation stage in the wide range of solution viscosities and nucleation rates with full allowance for the Laplace pressure in small overcritical bubbles. As is seen from Equations (17) and (21), the quantities f (R, t) and κ(t) are determined by the dependence of the individual bubble radius R(t) on time (related to the bubble growth rate . R(R)) and the dependence of the excluded volume ratio q(R) (related to the local relative decline ϕ(r, R, t) of gas supersaturation ζ(r, R, t) around the bubble with radius R from the bulk supersaturation value ζ 0 ). To illustrate the predictions of the theory considered in Section 2 at different values of solution viscosity η and different nucleation rates I 0 , we carried out several numerical calculations. First, we substituted Equation (6) in Equation (7) (6) provided the correction function h(R). Then, with the help of Equation (9), we found the local relative decline ϕ(r, t) of gas supersaturation around a single bubble in variables ρ and R, and substituted the results in Equation (12) to calculate the excluded volume ratio q as a function of the bubble radius R. As a next step, we used the results for . R(R) and q(R) in Equations (17) and (18) to calculate the distribution function f (R, t) of overcritical bubbles and the dependence of the mean bubble radius R(t) on time. Equation (19) allows us to find the duration t 1 of the nucleation stage. Finally, with the help of Equation (20) we computed the dependence of the swelling coefficient κ on time t.

R(R) functions in Equation
The following values of parameters for the solution and the dissolved gas have been taken in calculation: that gave us Evidently, the parameters in (22) can have other values that are acceptable from a physical point of view. Choosing values (22), we did not pursue the goal of describing any particular system but only illustrating the role of viscosity in a wide range of its values. The accepted values of the parameters are close to those that can be appropriate for solutions of highly soluble gases in water.
Generally, the initial nucleation rate I 0 is the function of initial supersaturation ζ 0 [1,22] which is mostly characterized by derivative Γ ≡ −ζ 0 (∂∆F c /∂ζ 0 ) (mentioned in Section 2). Pre-exponential factor in the relation between I 0 and activation barrier ∆F c can give a second-order correction. Taking in mind the discussion in [31] concerning the relation between the activation barrier ∆F c for the classical theory of homogeneous bubble gas nucleation in gas-supersaturated liquid solution, the gas solubility, and supersaturation, we decided to take two values for the initial nucleation rate I 0 :I 0 = 10 5 m −3 s −1 and I 0 = 10 7 m −3 s −1 at the same values of ζ 0 and Γ fixed in Equation (22). The influence of the solution viscosity η has been checked at three values of viscosity: η = 1 Pa·s, η = 10 3 Pa·s, and η = 10 5 Pa·s.
Note that it is possible to study the degassing of various gases in a solvent, and this allows us to formally consider the case when the diffusion coefficient is fixed at different values of the viscosity of the solution. If we study the degassing of a given gas in solutions with different viscosities, then we should take into account the relationship between the diffusivity and viscosity. In such a case, we have used in our calculations the Einstein-Stokes formula which predicts that the product of the diffusion coefficient for gas molecules and viscosity of solution is constant at other fixed parameters.
The results of the computations are presented below. The dependence of an individual bubble radius R(t) on time t is shown in Figure 2 for three values of the solution viscosity. The dashed blue curve 0 depicts the limiting case at η = 0 and R >> R * which corresponds to self-similar growth of bubbles [17,19,29]. In this limiting case, all viscosity and capillary effects are absent. Solid curves 1-3 correspond at fixed diffusivity D = 10 −9 m 2 /s to η = 1 Pa·s, η = 10 3 Pa·s, and η = 10 5 Pa·s, respectively. Curve 2 is plotted for the solution with η = 10 3 Pa·s, but the product Dη was taken to be the same as for the solution with D = 10 −9 m 2 /s and η = 1 Pa·s. It is seen that increasing the viscosity at fixed diffusivity decreases the bubble growth rate and diminishes R(t). However, this effect is significant only for curves 3 and 2 . Evidently, nucleation rate does not affect R(t). corresponds to self-similar growth of bubbles [17,19,29]. In this limiting case, all viscosity and capillary effects are absent. Solid curves 1-3 correspond at fixed diffusivity is equal to unity for all bubble sizes. As the distance r from the bubble center increases, the value of ( )  Figure 3a,b illustrate how the behavior of the local relative decline ϕ(r, R(t)) of gas supersaturation changes with increasing viscosity at fixed diffusivity. If we neglect the capillarity and viscosity, the gas concentration near the bubble is equal to the equilibrium value c ∞ and, respectively, the local supersaturation of the gas near the bubble at r = R(t) is zero. Then, the local relative decline ϕ(r = R(t), R(t)) is equal to unity for all bubble sizes. As the distance r from the bubble center increases, the value of ϕ(r, R(t)) drops to zero. This pattern corresponds to the behavior of ϕ(r, R(t)) on surface 0 in Figure 3a, and we see that the decay of ϕ(r, R(t)) is quite fast for the values of the parameters taken from Equation (22) (nucleation rate does not affect ϕ(r, R(t))). It is also seen that rising viscosity leads to a slower decrease in the value of ϕ(r, R(t)) in the space around the bubble, especially for smaller bubble sizes. The presence of folds on surfaces 1 and 2 in the region of small droplet radii in Figure 3a is just due to viscous and capillary effects. Compared with the dependence of the radius on time shown in Figure 2, the influence of these effects on ϕ(r, R(t)), even at a fixed diffusion coefficient, is noticeably more significant.  As is seen in Figure 3b, viscosity has much a stronger influence on the value of ( )  Figure 3. (a) The 3D-surfaces of local relative decline ϕ(r, R(t)) of gas supersaturation at distance r around the bubble with radius R(t). Surface 0 plots the limiting case at η = 0, D = 10 −9 m 2 /s, and R >> R * . Surfaces 1 and 2 correspond to η = 1 Pa·s and η = 10 3 Pa·s at D = 10 −9 m 2 /s. (b) The 3D-surface of local relative decline ϕ(r, R(t)) of gas supersaturation at distance r around the bubble with radius R(t) at η = 10 5 Pa·s and D = 10 −9 m 2 /s. (c) The 3D-surfaces of local relative decline ϕ(r, R(t)) of gas supersaturation at distance r around the bubble with radius R(t) for solutions with η = 1 Pa·s, D = 10 −9 m 2 /s (surface 1) and η = 10 3 Pa·s, D = 10 −12 m 2 /s (surface 2 ).
As is seen in Figure 3b, viscosity has much a stronger influence on the value of ϕ(r, R(t)) at η = 10 5 Pa·s. However, we can observe here that ϕ(r, R(t)) still monotonically drops along the lines with fixed R(t) as in Figure 3a, and this fact is in accordance with the general behavior of the solution of the diffusion equation. The nonmonotonic behavior of function ϕ(r, R(t)) as a function of the bubble radius R(t) for small radii at a fixed distance r/R(t) from the bubble center, which is well observed in Figure 3b (and weakly noticeable in Figure 3a), is associated with the effect of viscosity on the bubble growth rate . R in the factor 1 − (R c /R) 1 + 2η . R/σ in the expression (9). We also compared surfaces of ϕ(r, R(t)) for solutions with fixed product Dη at η = 1 Pa·s, D = 10 −9 m 2 /s and η = 10 3 Pa·s at D = 10 −12 m 2 /s (i.e., in the case of coupling of diffusivity and viscosity). The result is shown in Figure 3c. As is seen, the functions ϕ(r, R(t)) in scaled variables r/R(t) and R(t)/R * virtually coincide in this case. This is not a random result. A study of the asymptotical forms of the growth rate for small and large R, carried out in [30], showed that they are directly proportional to the diffusion coefficient, and only the terms proportional to the product Dη are present in the corrections. Our calculations here additionally show that the whole curves for . (6), for correction function h(R) coincide at different viscosities and diffusion coefficients with a fixed product Dη.

R(R)/D and, in view of Equation
The next important characteristic of gas absorption by a bubble is the excluded volume ratio q. The dependence of the ratio q on bubble radius R(t) at fixed diffusivity D = 10 −9 m 2 /s and three values of viscosity η = 1 Pa·s, η = 10 3 Pa·s, and η = 10 5 Pa·s is shown in Figure 4 (curves 1, 2, and 3, respectively). Curve 2 is related to the solution with η = 10 3 Pa·s and D = 10 −12 m 2 /s. As is outlined by the blue color, this curve virtually coincides with the curve 1 for the solution with η = 1 Pa·s and D = 10 −9 m 2 /s (with the same product Dη). As follows from definition (11), this coincidence is a consequence of the coincidence of the surfaces ϕ(r, R(t)) for curves 1 and 2 . The value q 0 = 0.3399, which is shown by line 0, corresponds to self-similar growth of bubbles in the limiting case at η = 0 and R >> R * at the diffusivity D = 10 −9 m 2 /s with neglection of all viscosity and capillary corrections. This value is asymptotic for other curves. Increasing the viscosity at fixed diffusivity increases the bubble size to reach this asymptotic value.      Figure 5a illustrates the behavior of f (R, t) at η = 1 Pa·s (surface 1) and η = 10 3 Pa·s (surface 2) relative to the distribution function for the limiting case at η = 0 and R >> R * (surface 0) with neglection of viscosity and capillarity effects. We used the semi-transparency of the plotted surfaces. Therefore, it is clear that under the specified parameters, all distribution functions are similar and successively nested into each other. Figure 5b shows the behavior of f (R, t) at η = 10 5 Pa·s. We see that in all cases considered, in accordance with Formula (17), the presence of the Heaviside theta function leads to the appearance of a front, beyond which the distribution function is equal to zero. At any moment of time t there cannot be bubbles with radii greater than R(t). At the distribution front, the function f (R, t) reaches its maximum values for each bubble radius R. A more specific behavior of f (R, t) at small bubble radii in Figure 5b is related to the fact that the bubble growth rate . R(t) is very small for these sizes (see analytical results in [30]).   Of particular interest is the behavior of the distribution function of gas bubbles at higher nucleation rates and when the diffusion coefficient of gas molecules and the Of particular interest is the behavior of the distribution function of gas bubbles at higher nucleation rates and when the diffusion coefficient of gas molecules and the viscosity of the solution are coupled. The results of the corresponding calculations of f (R, t) for solution with η = 10 3 Pa·s and D = 10 −12 m 2 /s are shown in Figure 5c in relation to the distribution function for the case with neglection of viscosity and capillarity effects at D = 10 −12 m 2 /s. Note that the maximal value of time t in Figure 5a-c corresponds to the duration of the nucleation stage during solution degassing. According to Equation (19), the characteristic time t 1 of the nucleation stage at parameters (22) and nucleation rate I 0 = 10 5 m −3 s −1 is 0.30 s for a solution with at η = 1 Pa·s (the maximal bubble radius is 241R * ), 0.34 s for a solution with η = 10 3 Pa·s (the maximal bubble radius is 244R * ), and 0.66 s for a solution with η = 10 5 Pa·s (maximal bubble radius is 191R * ). At D = 10 −12 m 2 /s and η = 10 3 Pa·s, we find t 1 = 3.62 s and maximal bubble radius is R(t 1 ) = 24R * . In comparison with surface 2 in Figure 5a related to I 0 = 10 5 m −3 s −1 , we see in Figure 5c a sharp rise in the number of small bubbles with a significant increase in the degassing time. Figure 6 shows the dependence of the mean bubble radius R(t) on time at the nucleation stage, calculated by Formula (18). The numbering of the curves (with used values of physical parameters) in Figure 6 corresponds to the numbering of the curves in Figure 2 for the dependence of the bubble radius on time t. Obviously, the curves in Figure 2 can be interpreted as dependences of the maximal bubble size in the ensemble of bubbles.   (18). The numbering of the curves (with used values of physical parameters) in Figure 6 corresponds to the numbering of the curves in Figure  2 for the dependence of the bubble radius on time t . Obviously, the curves in Figure 2 can be interpreted as dependences of the maximal bubble size in the ensemble of bubbles. Comparing the curves in Figure 6 with the corresponding curves in Figure 2, we conclude that the average radius in an ensemble of bubbles always grows noticeably slower than the maximum bubble radius. An increase in viscosity slows down the growth of the mean radius ( ) R t , especially in the case of inverse proportionality between the diffusion coefficient and the viscosity of the solution. The maximal value of time t for curves in Figure 6 corresponds to the duration 1 t of the nucleation stage during solution degassing.
Finally, Figure 7 shows the time dependence of the swelling coefficient of the degassing solution during the nucleation stage. The numbering of the curves under different degassing conditions coincides with the numbering in Figure 6. It should be noted that at Comparing the curves in Figure 6 with the corresponding curves in Figure 2, we conclude that the average radius in an ensemble of bubbles always grows noticeably slower than the maximum bubble radius. An increase in viscosity slows down the growth of the mean radius R(t), especially in the case of inverse proportionality between the diffusion coefficient and the viscosity of the solution. The maximal value of time t for curves in Figure 6 corresponds to the duration t 1 of the nucleation stage during solution degassing.
Finally, Figure 7 shows the time dependence of the swelling coefficient of the degassing solution during the nucleation stage. The numbering of the curves under different degassing conditions coincides with the numbering in Figure 6. It should be noted that at all three values of viscosity (fixed and inversely proportional to the viscosity diffusivity), swelling of the solution at the nucleation stage rises approximately to close values, although at a significantly different rate.  The fact that the swelling coefficient of curve 2′ with inverse proportionality between the diffusion coefficient and viscosity reaches greater values by the end of the nucleation stage than the swelling coefficient of curve 2 at a fixed diffusion coefficient is due to the use of a higher nucleation rate in the calculations for the case 2′ . As a result, we have in the case of 2′ a huge number of small bubbles to the end of nucleation stage.

Conclusions
Our previous analysis [29,30] of individual growth rates of gas bubbles at degassing in viscous supersaturated-by-gas solutions under high initial gas supersaturation and gas solubility has been extended to study of cooperative effects for an ensemble of growing gas bubbles. As we have shown here, viscous effects may noticeably affect the distribution function of supercritical gas bubbles in a degassing solution on the nucleation stage, and we gave a description of the mechanism of a such influence. In particular, as a macroscopic effect at degassing, viscosity provides slowdown of the maximal and mean bubble growth as well as the rate of swelling of the total solution. With increasing solution viscosity by five decimal orders at fixed diffusivity, viscosity does provide the damping effect on the rate of growth of the mean bubble radius and the rate of swelling up to two times and, correspondingly, increases twice the time of the nucleation stage. Taking inverse coupling between diffusivity and viscosity into account allowed us to find much stronger dumping up to one decimal order for three decimal orders of viscosity. Nevertheless, we have shown that to the end of nucleation stage, the swelling coefficient will be several units. Thus, in spite of the fact that the self-similar theory of diffusional growth of bubble [17][18][19] (which requires 0 η = and R R * >> ) cannot be directly applied at large viscosities, many of its features are preserved.
The results obtained allow us to consider any specific system for which the physicochemical parameters, such as gas solubility, surface tension, diffusion and viscosity coefficients, nucleation rate, and work of formation are self-consistent. At present, such a set of values of parameters can only be obtained experimentally.
In this paper, we considered the case of instantaneous decompression of a gas-saturated solution. In the case of a finite decompression rate, the approach proposed here The fact that the swelling coefficient of curve 2 with inverse proportionality between the diffusion coefficient and viscosity reaches greater values by the end of the nucleation stage than the swelling coefficient of curve 2 at a fixed diffusion coefficient is due to the use of a higher nucleation rate in the calculations for the case 2 . As a result, we have in the case of 2 a huge number of small bubbles to the end of nucleation stage.

Conclusions
Our previous analysis [29,30] of individual growth rates of gas bubbles at degassing in viscous supersaturated-by-gas solutions under high initial gas supersaturation and gas solubility has been extended to study of cooperative effects for an ensemble of growing gas bubbles. As we have shown here, viscous effects may noticeably affect the distribution function of supercritical gas bubbles in a degassing solution on the nucleation stage, and we gave a description of the mechanism of a such influence. In particular, as a macroscopic effect at degassing, viscosity provides slowdown of the maximal and mean bubble growth as well as the rate of swelling of the total solution. With increasing solution viscosity by five decimal orders at fixed diffusivity, viscosity does provide the damping effect on the rate of growth of the mean bubble radius and the rate of swelling up to two times and, correspondingly, increases twice the time of the nucleation stage. Taking inverse coupling between diffusivity and viscosity into account allowed us to find much stronger dumping up to one decimal order for three decimal orders of viscosity. Nevertheless, we have shown that to the end of nucleation stage, the swelling coefficient will be several units. Thus, in spite of the fact that the self-similar theory of diffusional growth of bubble [17][18][19] (which requires η = 0 and R >> R * ) cannot be directly applied at large viscosities, many of its features are preserved.
The results obtained allow us to consider any specific system for which the physicochemical parameters, such as gas solubility, surface tension, diffusion and viscosity coefficients, nucleation rate, and work of formation are self-consistent. At present, such a set of values of parameters can only be obtained experimentally.
In this paper, we considered the case of instantaneous decompression of a gassaturated solution. In the case of a finite decompression rate, the approach proposed here requires modification, primarily related to the change in the growth dynamics of a single gas bubble. The corresponding problem of single bubble growth in a highly viscous liquid was considered in recently published paper by Chernov, Davydov, and Pil'nik [32].