1. Introduction
A gas bubble is formed when an atomically or molecularly dissolved gas becomes supersaturated in a liquid solvent as a result of the reduction in imposed gas pressure, change in liquid temperature, or change in solute or solvent character Rosner et al., 1972, [
1]. The study of gas bubbles is of major interest due to their appearance in many real-world problems. One of the important applications of bubble hydrodynamics is in chemical process industries, for example in the production of foamed plastics Elshereef et al., 2010, [
2]. When a gas-generating substance such as a blowing agent is mixed with a high-pressure molten polymer, the resulting product turns out to be thermoplastic Arefmanesh et al., 1992, [
3]. In this process, gas bubbles emerge and have a considerable effect on product quality. Therefore, it is necessary to understand the behavior of bubbles under different process parameter conditions. High-density foamed thermoplastics, otherwise called cellular plastics, are used in household furniture, transportation, and building products; on the other hand, low-density thermoplastics are frequently used in rigid packing Lee et al., 1996, [
4]. The formation and growth of bubbles due to de-gassing or reduction in pressure in a supersaturated gas–liquid solution is observed in a broader spectrum of industrial and natural processes. For example, a very well-known process in which de-gassing is observed are carbonated beverages, such as beer, soda, and champagne (Bisperink et al., 1994, [
5]; Jones et al., 1999, [
6]; Barker et al., 2002, [
7]; Liger-Belair 2005, [
8]; Lee et al., 2011, [
9]; Enríquez et al., 2013, [
10] Enríquez et al., 2014, [
11]). The study of bubble dynamics is vital in production industries, where molten polymers, metals, and glasses are of major interest Amon and Denson, 1984, [
12] and a bubble prediction theory is important in the exsolution of gases during oil extraction Pooladi-Darvish et al., 1999, [
13].
Several mathematical models have been developed to predict the bubble size evolution in various industrial processes. For instance, Epstein and Plesset, 1950, [
14] derived an approximate analytical solution, by neglecting inertia, to an unbounded single bubble growth/dissolution in a gas–liquid solution due to pure mass transfer (diffusion) for supersaturated and undersaturated conditions. Epstein’s formulation suggests that the bubble grows as the square root of time, i.e.,
, where
is the radius of the bubble. However, their formulation lacks in explaining the hydrodynamic effects on bubble growth, including inertia, surface tension, etc. Barlow and Langlois, 1962, [
15] were the first to combine diffusion with hydrodynamics, wherein they introduced a very complicated integro-differential equation based on a thin shell assumption. The formulation of Barlow et al. is complicated and computationally time-consuming to solve for larger bubble growth rates. Rosner and Epstein, 1972, [
1] assumed a parabolic concentration profile in a thin boundary layer to generate an approximate solution of the diffusion equation. This work has been adopted by many researchers including Elshereef et al., 2010, [
2], Patel, 1980, [
16] and Han and Yoo, 1981, [
17] formulation does not account for the change in gas pressure inside the bubble with time. Patel, 1980, [
16] developed two coupled ordinary differential equations (ODEs) for predicting the unbounded growth of a single bubble in a Newtonian liquid; however, he neglected the effect of inertia in his formulation. Later, Amon and Denson, 1984, [
12] introduced a cell-based model that incorporates the effect of available gas from the surrounding bubbles. Amon and Denson’s formulation is developed based on a cell model assumption, where they have considered the foam as a summation of an equal microscopic unit of spherical cells with a constant mass in it and every cell has a spherical gas bubble that grows by diffusion of gas from the microscopic unit.
As discussed earlier, Barlow et al., 1962, [
15] and Patel, 1980, [
16] developed models for pure Newtonian fluid cases, hence neglected the effect of the elastic nature of the polymer. To fill this gap, Han and Yoo, 1981, [
17] and Ramesh et al., 1991, [
18] introduced a model that includes the effect of the elasticity of the fluid (polymer). Elshereef et al., 2010, [
2] compared two popular bubble growth models. The first model is known as the Patel model or single bubble growth model, which is developed on assumption that a single bubble grows in a pool of liquid with infinite availability of gas, and the second model is called a cell model or Amon and Denson model, which is developed by incorporating the finiteness of gas availability and considering the proximity of gas bubbles. The main motivation of the Elshereef et al., 2010, [
2] investigation was to compare these two models in terms of numerical implementations and accuracy in bubble growth prediction. In this regard, they compared the models with Han and Yoo’s experimental findings. In recent years, Soto et al., 2019, 2020, [
19,
20] investigated experimentally carbon dioxide (CO
2) and nitrogen (N
2) bubble growth in water solutions with and without confinements. Their finding suggests that after the initial period of diffusion-driven bubble growth, the mass transfer is further accelerated due to density-driven convective flow.
Although researchers have performed ample work in understanding the hydrodynamics of the bubbles in different processes, clear insight into the diffusion process coupling with hydrodynamics and an explanation of process flow parameters’ effects on hydrodynamics are lacking. The current work emphasizes solving the diffusion process numerically and closely studying how different liquid parameters such as liquid viscosity, surface tension, diffusion coefficient, system pressure, and solubility of gas affect the hydrodynamics of bubble growth. Though the current numerical framework developed in this work is for Newtonian liquids, the authors aim to explore how the current model compares with the different Newtonian liquid models of Elshereef et al., 2010, [
2] and the viscoelastic experimental data of Han and Yoo, 1981, [
17].
3. Non-Dimensionalization and Solution Procedure
3.1. The Dimensionless Problem
The equations and their initial and boundary conditions are non-dimensionalized as follows. The velocity scale is taken as
, which is related to the initial driving pressure difference, the length scale is the initial bubble radius
, the pressure scale is the initial gas pressure
, and the concentration scale is the equilibrium concentration
. In this case, the time scale is naturally
. The dimensionless variables become
There are five non-dimensional groups appearing in the problem, three familiar groups: the Reynolds number
, the capillary number
and the Péclet number
. More explicitly:
Here, the Reynolds number () compares the inertial force due to bubble growth in the liquid region with the liquid viscosity. The capillary number () weighs between viscous forces from the liquid to the surface tension forces at the interface of the bubble and the liquid and the Péclet number describes the ratio between the convection mass transfer to the diffusive mass transfer of gas from the liquid into the bubble.
The additional two new non-dimensional parameters are denoted by
and
, the former being the ratio of the initial pressure to the pressure difference, and the latter reflects the initial level of gas solubility:
Finally, a sixth additional parameter in the problem is the dimensionless atmospheric-to-gas pressure ratio .
3.2. Domain Mapping
The interface of the bubble changes with time, which makes the numerical procedure for solving the concentration distribution in the liquid more complicated and time-consuming. We implement an implicit finite difference in space and integrate the resulting equations with respect to time. One obvious but costly approach is to track the interface of the bubble with time and re-mesh the computational domain at each time step.
Alternatively, we recast the concentration Equation (13) in terms of Lagrangian coordinates
, such that at all time intervals the interface is fixed. Therefore, after non-dimensionalization and coordinate transformation the Equations (7), (11) and (13) takes the form:
The rescaled initial and boundary conditions are deduced from (3), (12), and (14) to:
3.3. Numerical Implementation
Equation (20) is a non-linear, second-order ODE that describes the bubble growth. If the pressure in the bubble is constant, Equation (20) can be solved for the bubble growth and its interface velocity with the use of any readily available numerical time integration solver, such as ode45 in MathWorks MATLAB version R2019b. However, the difficulty arises when the pressure inside the bubble varies with time, and it then needs to be coupled with the scalar diffusion equation to solve for the concentration gradient at the interface. Additionally, the scalar diffusion Equation (22) contains a highly non-linear convective term in terms of bubble radius and interface velocity. This combination makes the equations stiffer and involves solving Equations (20)–(22) simultaneously. Therefore, solving the highly stiff equations with ode45 takes a tremendous amount of time. Instead of ode45, a variable order of accuracy solver, ode15s, is used to integrate the equations. Here, ode15s uses first to fifth orders, changing the order as required. This solver takes much less time compared to the ode45 solver without compromising accuracy.
To solve these two equations simultaneously, the second-order non-linear hydrodynamic Equation (20) primarily needs to be converted into the system of first-order ODEs by letting
. Therefore, the system of first-order ODEs is given as
This way, when Equation (25) is integrated, one can obtain , which is bubble interface velocity, and similarly Equation (26) is integrated to obtain , which is the bubble radius. Since Equation (26) includes partial derivates in time and space, one can approximate either time or space using the finite difference methods. For convenience, the space partial derivative is approximated with a finite difference, up to second-order accuracy.
Let
be the node position and
be the total number of nodes (see
Figure 2) in the gas–liquid solution, starting from the interface
= 0 to infinity. The central difference scheme is adopted for the derivates. Therefore, the finite difference approximation for the first- and second-order derivatives with central difference schemes are written as
The discretized form of the scalar diffusion equation using Equations (27) and (28) takes the form
The discretized form of diffusion Equation (29) needs to be solved at
−2 (
) nodes, starting from
= 2 to
=
− 1. Whereas at the interface, i.e., at
= 1, the boundary condition (24a) can be written in terms of ODE as
The final node serves as a boundary and the value of concentration is known from the boundary condition (24b), therefore at
=
,
Similarly, the concentration gradient at the interface in Equation (21) is discretized using a forward finite difference scheme and is given as
and substituting Equation (32) in (21) results in
To be consistent with the notation used for the hydrodynamic ODEs (25) to (26), Equations (29) to (33) are rewritten in terms of y as follows:
For the nodes between 1 and
(
) is written as
at the interface node
,
and at the final boundary node
,
Finally, the pressure equation takes the form:
Therefore, the total ( +3) equations starting from (25) to (37) are the final system of ODEs that has to be solved simultaneously subjected to the initial and boundary conditions (23) to (24).
3.4. Grid Independence Test
For the numerical simulations, the infinite spatial domain is assumed to be 10 times the maximum radius of the bubble. Furthermore, the maximum radius of the bubble is anticipated to be 250 µm. This suggests that the physical infinity of the domain is 250
10 = 2500 µm and in terms of
it is 2250. (Note that
). The grid independence test seeks to minimize discretization error by making the numerical solution independent of the grid spacing.
Figure 3 shows that the solution converges with increasing number of nodes. When the domain is discretized from 100 to 1000 nodes a 15% of maximum error is observed in the bubble radius and the error reduced to 2% as the number of nodes increased from 1000 to 3000; Therefore, to achieve accurate results in the numerical simulations, the domain is equally discretized with 3000 nodes.
5. Concluding Remarks
The hydrodynamics of a single bubble in the pool of Newtonian liquid that expands due to mass transfer was investigated in the current work. This study directly relates to foaming processes, carbonated beverages, and any other problem in which the bubble grows due to mass transfer.
Rigorous non-dimensional formulations were derived to incorporate interfacial, viscosity, diffusivity, and solubility effect on bubble growth. Especially the inertia of the liquid was included in the formulation, along with full scalar advection–diffusion processes. A strong numerical approach to the highly non-linear stiff coupled equations was discussed. The moving interface of the bubble was tackled by mapping the domain to the new coordinate (x).
The results obtained with the present formulation and numerical solution to the advection–diffusion equation was compared with the Elshereef et al., 2010, [
2] models. The present numerical model predicts accurate bubble growth in comparison to Elshereef et al., 2010, [
2] models. These results were validated by comparing with the Han and Yoo, 1981, [
17] experimental data set.
To our knowledge, the influence and behavior of the concentration of the gas in the liquid has not been reported in the literature. In this work, a clear insight is provided on the concentration profiles of gas in the liquid and a boundary layer variation around the bubble. A simple numerical investigation was conducted to compare the variation in the approximated diffusion equation results against the present numerical results. We showed that that the gas concentration profile in the liquid deviates from the traditional concentration profile.
With the validated numerical model, a comprehensive parametric study was performed on the bubble growth. The results show that the rate of bubble growth depends primarily on the viscosity of the liquid, initial pressure difference, diffusion, and solubility. The effect of surface tension on the overall bubble growth process is limited.
We showed that the higher viscosity of the liquid lowers the bubble growth rate, and vice versa. The initial pressure difference between the bubble and the system has a significant effect on the overall bubble growth process. The higher the initial pressure difference, the greater is the bubble growth. With a lower initial pressure difference, the bubble growth is limited.
The investigation shows that the effect of diffusion and solubility of the gas in the liquid play an important role in the overall bubble growth process. Higher magnitude of these parameters leads to a higher bubble growth rate, and vice versa. It is concluded that these parameters have a similar effect on bubble growth.