Numerical investigation of the dissolution mechanism of a freely oscillating CO 2 gas bubble by the method of lines

The dissolution process of a CO2 gas bubble undergoing adiabatic free oscillation is investigated by coupling bubble dynamics with convective diffusion mass transport in the liquid phase. The minimum equilibrium radius necessary to the adiabatic free oscillation and the related damped frequency are determined. By formulating the governing equations using stretched Lagrangean spatial coordinates, the model system can be effectively solved by the method of lines. It was found that the gas dissolution significantly enhances the damping of the free oscillations. In addition, a periodic detachment of the concentration boundary layer takes place at the bubble wall, which induces short periods of gas desorption within a cycle of the oscillation. We point out that this causes a retarding effect on the dissolution process.

The importance of research on bubble dynamics coupled with mass transport originates mainly from the chemical and bioprocessing industry, since gas and fluid phases are present simultaneously in many operations.The efficiency of such processes can be optimized by increasing the interfacial density using bubbly mixtures of the reagents.Bubbles can be classified according to their content into two basic types.Vapor filled bubbles are mainly formed in boiling fluids in heat exchangers or rectification columns, while gas filled bubbles are mostly generated by external injection in technologies used for gas purification [1], wastewater treatment [2] or synthesis of chemical materials [3].
External forcing of bubble oscillation can produce enhanced evaporation of liquid and release of dissolved gases, thus it expedites to the phenomena of boiling and cavitation.For instance, the ultrasonic excitation of bubbles, as a specific procedure, is widely applied in sonochemistry [3], metal production [4] and ultrasonic medicinal treatment [5].However, the less frequently studied free oscillation of bubbles, which also affect the mass and heat transport processes, is also frequently encountered in many operational conditions (e.g. in cases of bubble detachment, sudden decrease in ambient pressure etc).
The most typical mode of bubble oscillation is the volumetric one, which can be satisfactorily modeled by assuming spherical symmetry.This problem has been studied thoroughly in the literature.The well-known Rayleigh-Plesset equation [6] describes this system as a specific nonlinear oscillator.Minnaert [7] was the first one who, by neglecting the surface tension effect, determined the natural frequency of the bubble oscillation in the small amplitude linear range.Prosperetti [8] extended Minnaert's analysis to the nonlinear range by including the surface tension and dissipation effects.Lauterborn [9] provided detailed numerical results for the free oscillation with large amplitudes.Chang and Chen [10], using the Hamiltonian function proposed by Ma and Wang [11], provided complete evaluation to the various solutions of free oscillation problem.Vokurka [12] has released a comprehensive study on the evaluation methods of the physical parameters measured for free oscillation of a bubble.Recently, Hegedűs et al. [13] revealed in his work the effect of ambient temperature and pressure on the damping of free oscillation by both finite difference and spectral numerical methods.
The mutual effects of the mass transfer rate of dissolved gases and the forced oscillation of bubbles have also been extensively studied.Epstein and Plesset [14] prepared an analytical derivation for the passive dissolution of a quiescent gas bubble.In the case of forced oscillation, Blake [15] established that, beyond a certain threshold amplitude, desorption process takes place, causing the gradual expansion of the bubble.This phenomenon termed as rectified mass diffusion is closely investigated by Eller and Flynn, [16].They applied a nonlinear timeaveraging method in their calculation of the mass transport by neglecting the high frequency oscillating terms.In addition, they assumed that the thickness of the diffusive boundary layer is negligible as compared to the bubble radius.Fyrillas and Szeri [17] relaxed the limiting assumptions included in Eller's derivation by splitting the convective diffusion problem to an oscillatory and a smooth part, which were solved by singular perturbation method.Recently, Naji Meidani and Hasan [18] investigated both numerically and experimentally the rectified diffusion for an air gas bubble excited by ultrasound field and obtained qualitative agreement in the results.
Surprisingly, the interaction of gas diffusion and the unforced, free oscillation of bubbles is not described by itself in the literature, even though the investigation of the free oscillation is an essential case, by which the fundamental characteristics of the system dynamics can be revealed.This has motivated the authors to investigate the free oscillation of a gas filled bubble in this paper.
The reduction of CO 2 emission as the greenhouse gas emitted in the largest volume by human activities has become a primary focus in the last decades.Generally, in CO 2 gas purification technologies using chemical or physical absorption, the regeneration of solvent is carried out on high temperature, which is quite energy intensive and solvent consuming due to the evaporation loss.The ultrasonic excitation of the solution as a pioneering degassing technology increases the desorption rate of the dissolved gases and, therefore, enhances the developing of gas bubbles, [19].Taking advantage of this, significant energy and solvent mass can be saved.The combined environmental and industrial interest has motivated the authors to choose the system of CO 2 gas filled bubble in an aquatic environment as the specific example while investigating the free oscillation of a diffusing gas bubble.However, the model and the method outlined below can be easily adapted to other gas-liquid combinations.
This paper is organized as follows.First, we present and substantiate the coupled model of diffusion and bubble dynamics, leading to a coupled system of ordinary and partial differential equations.Next, the methods by which this mathematical problem can be treated numerically are discussed.Finally, the results obtained from the solution are analyzed.The symbolic notations we use are assembled in the Nomenclature, for brevity only those quantities are defined in the text, which deserve special consideration in this problem.

Model assumptions
The basic assumptions of our model setup include the following simplifications: 4 The concentration equilibrium is continuously maintained through the bubble wall.
5 No equilibrium reactions of the dissolved gas are considered in the liquid phase.
The material properties of the investigated system are considered at standard conditions (p = 101325 kPa, T = 29815 K).
Tab. 1 lists all the values corresponding to the specific choice of materials that occur in the model Two additional, less evident assumptions are discussed in more detail below in this section.
6. Initial state.The free oscillation of a bubble can be induced either by changing the internal energy of the bubble or by temporarily changing the ambient pressure, [12].In practice, it can occur for example when a wave of depression propagates in the liquid phase, as the result of which the bubble expands to the initial bubble radius on the expense of doing work on its surroundings and the interface.In this paper the excitation is achieved by setting the initial bubble radius R 0 larger than the 'equilibrium' value R eq .It is also assumed that initially the liquid contains no dissolved gas at all.Accordingly, under 'equilibrium' of the bubble, above, we have meant only full mechanical and thermal equilibrium with its liquid surroundings, without material balance.The quantity R eq defined this way thus depends exclusively on the total amount of gas and thus remains, as a useful reference, constant in time.
7.Adiabaticity of the bubble.Our work exclusively focuses on the size range of bubbles in which their thermodynamical behavior can be considered adiabatic.In case of a forced linear oscillation with small amplitudes, Prosperetti established the connection between the bubble oscillation and the thermal processes.He introduced the following two non-dimensional parameters, [27]: The physical meaning of G 1 can be interpreted as the square of the ratio of the thickness of the thermal boundary layer and the sound wavelength in the gas.Similarly, G 2 represents the square of the ratio of the bubble radius and the thermal penetration depth.The angular frequency (ω) in the former expressions corresponds to the driving frequency, which in the case of free oscillation has to be replaced by the oscillation frequency.The latter quantity can be reasonably well estimated by the damped frequency of the system, for which Prosperetti derived the approximate analytical formula, [8]. Here is the natural frequency of the bubble oscillation, ω 0 is its appropriately non-dimensionalized value, b is the dimensionless damping coefficient, C 0 is associated with the steady-state value of the amplitude of the 1/2 order subharmonic component, and h 3 is a specific function of the polytropic exponent k and the bubble's Weber number; all the necessary formulae for these quantities are given in [8].The product of G 1 and G 2 from (Eq. 1) yields which is equal to the square of the ratio of the sound wavelength in the gas and the bubble size and from which the uniformity of the gas pressure can be inferred.Eq. ( 4) can be combined with the practical criteria of Zhang and Li [28], which states that the pressure field is uniform in both phases if the respective sound wavelengths are at least ten times exceed the equilibrium bubble radius.This yields the following formulas: for the gas phase, and for the liquid phase.These relations determine the respective upper limits 0.39 and 12.1 for G 1 G 2 in the case of the CO 2 gas and water combination used.If, in addition to Eq. 6, the criteria are also fulfilled, than the thermodynamical properties of the bubble change adiabatically, according to Zhang and Li [28].These restrictions establish 2.8× 10 −4 m as a lower bound for the equilibrium bubble radius in the investigated case.Accordingly, we have decided to choose the equilibrium radius as Numerical investigation of the dissolution mechanism 65 in our calculations, which satisfies the adiabatic condition with a wide margin.For a bubble of this size, Prosperetti's formula (Eq.2) yields the damped angular frequency (or, equivalently, 3.087 kHz).

Governing equations
The oscillating bubble motion is described by the Rayleigh-Plesset equation [6], which can be cast to the nondimensionalized form [29] with the initial condition (10) by changing to the new dimensionless variables The spatial variable is normalized by the equilibrium bubble radius R eq that was specified in the previous section, c.f. Eq. ( 7), while the new time scale τ can be chosen freely.Actually, for τ, we use the dissolution time of the investigated CO 2 bubble without oscillation, the value of which is determined in the next paragraph.The constants (12) in Eq. 9 can be interpreted as characteristic time scales, which provide additional information about the magnitude of the effects influencing the oscillation of the bubble with R eq equilibrium radius.Considering the actual time scales of the system, which are tabulated in Tab. 2 for the R eq value chosen in Eq. 7, it can be established that the bubble motion is mainly controlled by the pressure terms in the initial state.However, due to the dissolution process, the mean bubble radius gradually decreases, which, in turn, causes some variation in the magnitude of the influencing factors (Eq.12).The dissolution mechanism of a stagnant gas bubble is well known, it was modeled by Epstein and Plesset, [14].In their derivation, the convective transport term was neglected, since the absorption is controlled exclusively by mass diffusion.In addition, the gas temperature is considered constant, due to the two order of magnitude smaller timescale of heat transport compared to that of the CO 2 diffusion in the water as expressed by the Lewis number Le = α l /D ≈ 80 (13) These approximations lead to the following dimensionless parametrical differential equation [14]:

Tab. 2. Characteristic times
where the non-dimensional constants and the scaled time variable are defined as Using the ODE15s solver in M R , we have obtained for the total dissolution time of a CO 2 filled bubble with the diameter specified in Eq. 7.
According to the model assumptions 2 and 7, the state of the gas in the bubble is described by the polytropic relation, where the bubble's mass normalized by its initial mass, is introduced as a new dynamical variable.The dispersion of the dissolved gas in the liquid phase is modeled by the convective-diffusive transport equation, which can be re-written in non-dimensional form as, [17]: where is the Péclet number of the system.Here the concentration is given in mass fraction (kg dissolved gas per kg solution), therefore it is already a dimensionless quantity.In accordance with assumptions 6 and 4, the initial and boundary conditions for (Eq.19) can be specified as follows: Note that initially the concentration profile is discontinuous at the bubble wall.The concentration equilibrium described by Henry's law is sustained, according to assumption 4, at the bubble wall, making the inner boundary concentration proportional to the gas pressure.The variation of the bubble mass (Eq.18) is determined by the diffusion mass flux across the bubble interface according to Fick's law: Proper normalization according to (Eq.12) and (Eq.25) yields the non-dimensional initial value problem: The obtained mathematical model contains two nonlinear, stiff ordinary differential equations (ODE's), Eq. 9 and (Eq.25), and a parabolic partial differential equation (PDE), Eq. ( 19), for which the initial and boundary conditions are set in Eq. 10,Eq.26,Eq.20,Eq.16 and Eq.23.The differential equations are coupled by the algebraic equation (Eq.17).Hereupon, this set of equations will be referred as the Eulerian model equations.From the physical point of view, it can be observed that, the mass transport affects indirectly the bubble dynamics through Fick's diffusion formula.However, in the reverse direction, the bubble motion controls the mass transport process not only through Fickian diffusion, but directly by the bubble radius, as well.
The moving bubble wall in boundary condition (Eq.23) (Stefan problem) can raise difficulties in the solution.To eliminate this problem, Fyrillas [17] introduced the new spatial coordinate If one neglects the mass flow rate contribution of the dissolved gas -which is in accordance with assumption 3 -, σ becomes a Lagrangean coordinate as it easily follows from the principle of continuity, [17].By using the new variable, one obtains the subsequent form for the mass transport equation the initial and boundary conditions of which are modified as follows: Similarly, the mass flux equation ( 19) must also be reformulated as By replacing (Eq.13), Eq. 21-Eq.23 and (Eq.25) from the set of the Eulerian model equations by (28-33) above, we obtain the governing equations of our Lagrangean model.

Applied numerical method
The method of lines is a fundamental numerical tool for solving PDE's and coupled ODE-PDE problems.The essential element of this method is that the spatial derivatives are discretized by applying various order finite differential schemes, which turns the PDE into an N dimensional coupled ODE system (N refers to the number of grid points).In case of an ODE-PDE problem, the original ODE's have to be included the system.Thereafter it is solved as an initial value problem by multistep numerical integrators.The adequate modeling of the concentration boundary layer requires a non-uniform grid.The contradicting demands can be resolved by applying yet another transformation to the spatial coordinate.In our case, the following stretching function, suggested by Gupta et al. [30], is applied: (33) This provides a transformation between an equidistant grid in ξ space and a grid that follows a geometrical series in the σ space.With the parameter setting where ∆σ in and ∆σ out Numerical investigation of the dissolution mechanism 67 correspond to the values of the grid spacing at the inner and outer boundaries and N is the number of grid points.The transformation (Eq.33) maps the intervals on each other.Note that for a practical solution, the outer boundary σ * has to be made finite.These settings have proved to be appropriate in the course of the model calculations.The transformation of Eqs.(Eq.128) and (Eq.32) to the stretched coordinates involves using the identities, [31]: Having re-written in the new coordinates, second order, central and upwind finite differential schemes are used to discretize these equations, respectively, thus retaining the numerical accuracy of the difference scheme.The final form of the discretized equations can be written in the following manner: Equations (Eq.37 and (Eq.38 together with (Eq. 9 form the system of ODE's to be solved.The initial values to these are given in (Eq. 10, (Eq.26 and the latter one realizes (Eq.21) in the numerical model.The boundary conditions appear in (Eq.38) through c 0 and c N , the concentrations at the bubble wall and at the outer boundary; these are determined by the algebraic constraints: Finally, the numerical model are closed by the algebraic polytropic equation (Eq.17).
The obtained system of coupled ODEs is numerically integrated by using again the ODE15s solver of M R , the absolute and relative error tolerances are set to 10 −10 and 10 −9 , respectively.
A detailed sensitivity analysis has been carried out to identify the effects of the grid parameters (σ * , ∆σ in ∆σ out , N).The settings σ * = 0.5 for the outer boundary and the value 10 −7 for the maximum grid stretching ratio ∆σ in ∆σ out have proved to be appropriate to resolve the concentration boundary layer in the course of the model calculations.The limiting factor is the number of grid points N. Quantitative information on this can be obtained by analyzing the truncation error ε of the numerical method by Richardson extrapolation, which was performed in the first point σ 1 of the non-uniform grid at t .Fig. 1 presents the concentration values and the relative truncation errors as the function of the number N of the grid points.The quantitative results of this examination are summarized in Tab. 3. Since the change of the concentration value calculated between cases c and d diminishes, we have concluded that our final calculations can be safely performed on the numerical grid consisting of N = 500 points.

Results and discussion
The numerical computations for the free oscillation following the initialization of a diffusing CO 2 bubble have been carried out for four cases of initial perturbation, the initial values corresponding to (Eq.10) are listed in Tab. 4. In order to demonstrate the effects of diffusive bubble oscillation we use two different types of references for comparison: (a) the oscillation of CO 2 bubbles without diffusion (i.e.constant mass bubbles) to show the effect of gas diffusion on bubble oscillation and (b) the diffusion from a non-oscillating ('stagnant') CO 2 bubble to show the reverse effect, that of the oscillation upon gas diffusion Reference time series of type (a) has been generated for each case of initial perturbation listed in Table 4 by solving the Rayleigh-Plesset equation (Eq.6) with adiabatic transitions (Eq.17), but without coupling the dissolution process, i.e. by keeping the initial mass of the bubble constant.Two reference solutions of type (b) have been obtained: an adiabatic one by solving the model system with the initial conditionR 0 = 1 and an isothermal one by solving the Epstein-Plesset equation (Eq.14).very close towards the non-oscillatory decay curves.(We note that the non-oscillatory reference solutions of type (b) -the adiabatic one and the isothermal one predicted by the Epstein-Plesset theory -are indistinguishable on the scale of Fig. 2) The frequency spectra obtained by FFT (Fast Fourier Transform) from each time series of reference (a) can be inspected in Fig. 3.It can be noted that in the cases with small amplitudes (cases I and II), the dominating frequency well agrees the damped frequency (8) calculated by Prosperetti's relation (Eq.2).This match provides an a posteriori verification of the choice Eq. 6 of the bubble radius and of assumption 7 as a whole.However, the increase in the initial perturbation of bubble radius brings about the intensification of subharmonic components, which indicates the nonlinear character of the oscillation.In addition, the initial frequency peaks undergo a significant broadening: analysis with short time FFT (not presented here) has revealed that this is actually a consequence of a gradual change of the eigenfrequencies during the relaxation.In the cases with coupled absorption, only a slight increase in the eigenfrequencies has been found (c.f.Fig. 3, which can be attributed to the bubble shrinkage due to the dissolution.The intensification of the damping process, discussed above, occurs to the subharmonic components as well.

Tab. 4. The investigated cases
The strength of the damping can be well characterized by the logarithmic decrement, which is defined as the logarithm of the ratio of two consecutive amplitudes: In case of an ordinary linear oscillatory system this quantity is constant: independent both of the time elapsed and of the initial state.In Fig. 4 the magnitude of the logarithmic decrement is presented for all the investigated cases of the oscillating CO 2 bubble, both with and without diffusion [comparison of type (a)].The figure clearly demonstrates the nonlinear character of the bubble oscillation, since, on the one hand, larger initial perturbations correspond to higher logarithmic decrement values and, on the other hand, the logarithmic decrements decrease in time.As for the effect of diffusion on the oscillation it can be established that the damping rate is more intense during the first one third part of the investigated time period.However, it can be also stated, that the damping rate decreases significantly later in the course of the dissolution In Fig. 5, the logarithmic decre- ment values plotted against the amplitude completely collapse to one curve each for the diffusive and reference cases.Hence it can be stated, that the temporal dependency presented in Fig. 4 is indirect, it occurs exclusively via the changing amplitude, not uncommon in non-linear oscillations.
To obtain an inside view of the mass transport process occurring through the bubble wall, the developing concentration boundary layer should be investigated in detail.In Fig. 6, the concentration profiles in six subsequent phases of the oscillation are presented (numbered from 1 to 6) during the 2nd oscillation period.The upper inset illustrates the variation of the bubble radius, while the panel on the right side shows the development of the mass flow rate through the bubble wall.In the compressive half period (before phase 3) the saturation concentration at the bubble wall increases by several times with respect to the initial value due to the increased gas pressure within the bubble according to Henry's law, Eq. 23.As Fig. 6 illustrates in case of large oscillations (Case IV), this jump in the wall concentration may be in the order of magnitudes relative to the initial value.This coincides with the development of a strong negative radial concentration gradient, yielding an increase of similar magnitude in the mass flow rate outward of the bubble [negative values, according to (Eq.24), since the order of the reduction of the bubble surface area during the same half period is much lower.In the expansive half period (after phase 3) the concentration maximum detaches from the bubble wall and this starts a damped spherical concentration wave, propagating outward and helping part of the dissolved gas transported deeper into the liquid phase.During the same time, while the bubble is expanding, a strong positive radial concentration gradient builds up inside of the concentration maximum causing an inwardly directed mass flow rate back to the bubble (positive values) The resultant of the two opposite phenomena causes a slight retardation of the absorption process altogether.Fig.7 illustrates the temporal development of the bubble mass in certain cases of the initial perturbation.It is clearly shown that, the more powerful the oscillation, the more the absorption of the gas bubble is retarded, that is the mass of the CO 2 bubble decreases to a lesser extent.Fig. 8 represents the mass ratio of the bubble with respect to the mass corresponding to the adiabatic dissolution of the bubble without initial perturbation [a comparison of type (b)] As the effect of the oscillation on the bubble mass, a relative increase can be noticed at the beginning up to a local maximum, the value of which grows with the initial perturbation magnitude.

Summary
In this study we have set up and numerically realized a model that allows the simultaneous investigation of both the oscilla-tions and the dissolution process of CO 2 filled gas bubbles in an aqueous medium.
We determined the minimum size of a CO 2 bubble for adiabatic free oscillations under standard conditions using Prosperetti's formula [27].We assumed zero dissolved gas concentration at the beginning and permanent concentration equilibrium controlled by Henry's law at the bubble wall.Our assumptions lead to a model that includes the advection-diffusion equation (PDE) coupled by the Rayleigh-Plesset equation (ODE), the diffusive flux equation (ODE) and the equation for the adiabatic change of state.Furthermore, following Fyrillas' suggestion [17], the mathematical model was transformed from Eulerian to Lagrangean coordinates to eliminate the spatial dependency in the boundary condition at the bubble interface.A nonuniform grid was generated by a suitable coordinate transformation in order to resolve precisely the concentration boundary layer and to retain the numerical precision.The obtained numerical model was solved by the method of lines.According to the error analysis, the use of 500 grid points has proved to be satisfactory for the required accuracy.
Our investigation comprises four cases with different strength of initial perturbation.The model calculations have indicated the followings.Most importantly, the diffusion from the bubble has a significant effect on its free oscillation: the oscillation undergoes much faster damping if it is accompanied with absorption.This enhanced damping process rapidly approaches the dissolution curve of a quiescent bubble.Conversely, the free oscillations also affect somewhat the diffusion process: it can be stated that strong oscillations have a retarding effect on the CO 2 absorption rate.Both the damping rate and the absorption delay increase with the oscillation amplitude.These findings demonstrate the mutual effect of free oscillations and absorption in case of an industrially relevant gas bubble-water system.The parameterization of this relationship may be proved useful in improved modeling, design or control of certain treatment processes in industrial applications.
Our coupled model is capable of reflecting several phenomena in this complex system.From the point of view of bubble dynamics it clearly reveals the nonlinear character of the oscillation: the existence of subharmonics on the one hand, and the fact that both the frequency and the damping rate depend on the amplitude on the other hand (The latter two quantities can be combined into an amplitude-dependent complex eigenfrequency.)Furthermore, the temporal change of the damping rate characterized by the logarithmic decrement is controlled by the oscillation amplitude.Concerning the diffusion/absorption process the following features are revealed: During the compressive phase of the oscillation the CO 2 absorption resulted from the increased gas pressure develops a concentration boundary layer while during the expansive phase a reverse mass flow towards the bubble occurs due to the dropping pressure which in turn yields the detachment of the boundary layer from the bubble wall.

Case R 0 Fig. 2
Fig. 2 illustrates the obtained oscillation patterns in the two extreme investigated cases: I and IV.The envelopes of oscillations are shown by the curves of oscillation extrema and the dark grey region in between represents the range of oscillations.A

Fig. 2 .Fig. 3 .
Fig. 2. Temporal change of the bubble radius in various perturbation cases.The oscillations for soluble and insoluble gas content are depicted by dark grey and bright grey regions, respectively.The Epstein-Plesset solutions are denoted by white dashed lines

Fig. 6 .
Fig. 6.Left: concentration profiles at certain phases of oscillation, compression and expansion are depicted by white and black symbols, respectively.The inset and the panel on the right show the simultaneous variation of bubble radius and mass flow rate.

Fig. 5 .
Fig. 5. Logarithmic decrement in the function of normalized amplitude for the varying and constant mass cases (dashed lines indicate the initial radii)