Mathematical Modeling of Oxygen Transfer Using a Bubble Generator at a High Reynolds Number: A Partial Differential Equation Approach for Air-to-Water Transfer

Abstract

This paper presents the mathematical modeling of the oxygen transfer process using partial differential equations (PDEs). This process is crucial in various environmental and engineering applications, such as wastewater treatment, aeration systems, and natural water bodies, in order to maintain water quality. The authors solved the typical PDE for describing the change in oxygen concentration over time and present the developed model of the differential equation with the term “source”, indicating that the model could be used to optimize oxygen transfer in various environmental and engineering applications, contributing to improved water quality and system efficiency.

1. Introduction

The urgent need to address the challenges of the global water sector continues to be the driving force behind digital transformation. Optimizing water use is a growing issue for industries and institutions around the globe. Most industries use data-driven applications and artificial intelligence for wastewater treatment and oxygenation. In the current era of digital transformation, the adoption of modern technologies has already revolutionized several industries, and the water treatment and oxygenation sectors are no exception.

As environmental concerns grow, the need to oxygenate water using green energy is linked to the challenges associated with managing wastewater in a more efficient and sustainable manner. The benefits arising from the exploitation of renewable energy systems can be divided into three categories: saving energy, generating new jobs, and reducing environmental pollution. The most important benefit of renewable energy systems is the reduction of environmental pollution.

The need for water oxygenation is highlighted by one of the problems faced by various stationary water sources, namely eutrophication, a critical global environmental problem [1,2,3], which plays a key role in decreasing the oxygen available for aquatic life.

To develop a mathematical model for expanding suitable areas for fish survival in a polluted river, the dynamics of both pollution concentration and dissolved oxygen (DO) levels in the water must be determined. The authors of [4] proved that high pollutant concentrations can be reduced by increasing the concentration of dissolved oxygen.

Blue-green algae are also called cyanobacteria. Their main morphological characteristic is that they have gas bubbles (pulsating) that play a role in buoyancy. Thus, they migrate according to favorable conditions in the water mass. They decant hard or not at all. They form large, macroscopic colonies on the surface of the water, like dense felt, apparently green. If the water is stagnant or the flow speed is low, this mat prevents gas exchange, the oxygen concentration drops, and the ecosystem is destroyed. Colonies are strong because these species produce a mucilaginous mass that holds them together, plus they attract other smaller colonies through electrostatic charging. When the blue algae die due to old age or a lack of oxygen (they produce it during the day but consume it at night), they quickly putrefy, produce gases, smell of gases, have a marshy appearance, and oxygenation of the waters becomes mandatory.

Aquaculture, characterized by a high stocking density of fish with a high input of materials and energy, has induced serious negative effects on the aquatic environment, such as eutrophication and oxygen depletion [5,6,7]. Low oxygen in water can kill fish and other organisms present in the water [8]. Recently, research has been conducted for the application of a solar-powered water oxygenation system for pond aquaculture in areas with limited resources [9,10].

Researchers have already initiated the development and testing of a study generated by modified technology, which will be of significant importance in the operation of a pond as an economic alternative for improving the aeration system [11]. Other hybrid sustainable energy systems for aquaculture farms using photovoltaic systems are presented in papers [12,13].

In previous research [14,15,16,17,18,19], the performance of microbubble generators built by special technologies was demonstrated. Theoretical developments and numerical simulations are presented in [20], where the authors demonstrate that aeration installations using mobile fine bubble generators are more efficient than classic ones with fixed fine bubble generators. Other experimental results can be found in [21], where it was found that the oxygen transfer rate in water increases when fine bubble generators are used. Also, the generation of smaller and more numerous bubbles, leading to the improvement of the water aeration process, was demonstrated in [22].

A fundamental transport process in fluid mechanics is diffusion. Diffusion differs from advection in that it is random in nature (it does not necessarily follow a fluid particle). Diffusion is a process of spontaneously initiated mixing (due to the principle of “Brownian motion”) of the particles of two substances, gradually acting until their complete and irreversible homogeneity. The two substances can be gaseous, liquid, or solid substances. Diffusion has two primary properties: it is random in nature, and transport occurs from regions of high concentration to regions of low concentration [23].

Gas transfer at the gas–liquid interface is described by the two-film theory, established by “Lewis and Whitman” [24,25]. In this model, it is considered that on each side of the interface there is a thin film. For the gas molecules to dissolve in the liquid, they must cross the two films by molecular diffusion.

According to Henry’s law, the concentration of dissolved gas in the liquid phase is in equilibrium with the concentration in the gas phase. Fick’s second law applies in the case of a non-stationary diffusion C = f(x,y,z,τ), when there are no driving forces to keep the concentration gradient constant. The law connects the variation in concentration in time and in space [26].

Differential equations with partial derivatives (or equations with partial derivatives—PDEs) are inevitably encountered when modeling real phenomena, the complexity of which involves the determination of functions with two or more independent variables. A more realistic description of the phenomena taking place in space requires at least two independent variables, and tracking their development may require the introduction of a new independent variable, time.

2. Materials and Methods

2.1. The Fine Bubble Generator with the Plate with Orifices in Rotational Motion

Fine bubble generators are widely used in various applications, such as wastewater treatment, aquaculture, and chemical processes. These devices create small bubbles that increase the surface area for gas transfer, enhancing the efficiency of processes such as oxygenation or aeration. When combined with a rotating motion, the distribution and dynamics of the bubbles can be significantly influenced, leading to improved performance.

In this sense, a fine bubble generator whose perforated plate is in rotational motion was designed and built. The orifices made on the plate were executed on a micromachining machine, with a positioning error of up to 10 μm. The orifices were arranged in a single row, with 18 orifices on one side of the fine bubble generator and 19 orifices on the other side. The orientation of the air emission orifices is opposite to the direction of rotation of the generator of rectangular shape (Figure 1a). The constructive solution of the fine bubble generator ensures the oxygenation of the volume of water above the orifice plate.

Considering the orientation of the rotation movement, the area A-B with orifices and the area C-D with orifices are found “downstream” of the direction of the rotation movement (Figure 1b). Each orifice will emit air bubbles at different tangential speeds.

In the xoy plane, the air bubble exits tangentially at the radius of rotation R.

The location of the orifices is crucial for ensuring uniform bubble distribution and effective mixing. The orifices are typically placed at a specific radial distance from the center of rotation. This distance decides the initial velocity of the emitted bubbles due to the rotational motion. The velocity distribution of bubbles emitted from the orifices is influenced by the rotational motion of the generator and the characteristics of the emitted air jet. The resultant velocity of a bubble immediately after emission is a combination of the tangential velocity due to rotation and the initial jet velocity.

The precise positioning of orifices and the resultant velocity distribution of the emitted bubbles are key factors that influence the efficiency of the bubble generation process. Proper design ensures uniform bubble distribution and optimal gas transfer rates.

To ensure the model is well understood and its limitations are clear, the following assumptions are made:

-

A high Reynolds number indicates turbulent flow, which enhances mixing and mass transfer rates.

-

The bubbles have a constant size, neglecting variations due to changes in hydrostatic pressure or coalescence/breakup dynamics.

-

The bubbles are uniformly distributed within the fluid.

-

Homogeneity in fluid properties (density, viscosity, diffusivity).

2.2. Solving the Differential Equation of Diffusion with the Term “Source”

The diversity of equations and systems of differential equations is extraordinary, and boundary and initial conditions (whose role is particularly important) only add new cases and types of problems to be solved. In the field of EDP (especially nonlinear equations), the mathematical theory is not sufficiently elaborated. For these reasons, the numerical approach cannot “keep up with theory”, recommending extreme caution in confirming the results obtained in cases not sufficiently covered by theory.

Equations with partial derivatives, known as partial differential equations (PDEs), can be classified according to several criteria based on their characteristics and the nature of their solutions. The primary classification criteria include linearity, order, type, and the number of independent variables. Here is a detailed overview of these criteria [27,28]:

-

According to the order of the partial derivatives, they are classified into first-order, second-order, and n-order equations;

-

According to the character of linearity, they are divided into linear, quasi-linear, and nonlinear equations;

-

According to the type of influence of the integration domain on the solution at a point, they are divided into elliptic, parabolic, and hyperbolic equations;

-

According to the general form of a linear differential equation with partial derivatives with two variables [29]

$$A{u}_{xx}+B{u}_{xy}+C{u}_{yy}+D{u}_x+E{u}_y+Fu=G,$$

(1)

where the coefficients are functions of the independent variables x and y (for example, A = A(x,y), B = B(x,y)), and

$$u_{xx}=\frac{{\partial }^{2}u}{\partial x^2},$$

(2)

The quantity B² − 4AC refers to the discrimination of the equation.

The differential equation of diffusion in the non-stationary regime in the vertical direction (oz) has the form

$$\frac{\partial C}{\partial \tau }=D\frac{{\partial }^{2}C}{\partial z^2},$$

(3)

so, it is a parabolic equation.

By examining linearity, order, type, number of independent variables, and the associated boundary and initial conditions, one can decide the most suitable approach to address a given PDE.

3. Results

Parabolic-type problems involve a variational and complementary formulation appearing in the mathematical models of several applications in engineering, economics, biology, and different branches of physics. These types of problems present a series of analytical and numerical difficulties related, for example, to the evolution in time and movement on the border [30,31,32]. There are a considerable number of analytical and numerical methods used to obtain the solution to a parabolic type of problem.

Since the finite difference method in the default variant is much more stable than the other variants, it will be used to solve the diffusion equation. In this scheme, the derivative with respect to space is studied at level j + 1.

To express the first-order partial derivative with respect to time, it is preferable to use the direct difference; therefore

$${\left.\frac{\partial C}{\partial \tau }\right|}_{\left(z_i,{\tau }_j\right)}\approx \frac{C_{z_{i,}{\tau }_{j+1}}-C_{z_{i,}{\tau }_j}}{\Delta \tau },$$

(4)

The spatial second derivative at z_i, τ_j can be approximated using a central difference:

$${\left.\frac{\partial C}{\partial \tau }\right|}_{\left(z_i,{\tau }_j\right)}\approx \frac{C_{z_{i+1,}{\tau }_j}-2C_{z_{i,}{\tau }_j}+C_{z_{i-1},{\tau }_j}}{{\left(\Delta z\right)}^2},$$

(5)

Using these approximations, the finite difference form of the equation at z_i, τ_j becomes

$$\frac{C_{z_{i,}{\tau }_{j+1}}-C_{z_{i,}{\tau }_j}}{\Delta \tau }=D\cdot \frac{C_{z_{i+1,}{\tau }_j}-2C_{z_{i,}{\tau }_j}+C_{z_{i-1},{\tau }_j}}{{\left(\Delta z\right)}^2},$$

(6)

Rewriting this to emphasize that the left-hand side is calculated at (z_i,τ_j)

$$C_{z_{i,}{\tau }_{j+1}}\approx C_{z_{i,}{\tau }_j}+\frac{D\Delta \tau }{{\left(\Delta z\right)}^2}\left[C_{z_{i+1,}{\tau }_j}-2C_{z_{i,}{\tau }_j}+C_{z_{i-1},{\tau }_j}\right],$$

(7)

To solve for the unknown $C_{z_{i,}{\tau }_{j+1}}$ using a finite difference method, one needs to discretize the partial differential equation appropriately (see Figure 2).

The solution C(z,τ) is computed using the values at the previous time step. This iterative process allows the solution to propagate forward in time, and we can visualize this propagation in a space–time plane.

To specify the dividing points of the intervals for the variables t and z, one can define them in a systematic manner. For this, the time step is considered Δτ and the spatial step Δz.

The time variable τ is defined as ${\tau }_j=j\cdot \Delta \tau ; j=0,1,2,\dots ,N_{\tau }$, where τ_j represents the time points and N_τ is the total number of time intervals.

The spatial variable z is defined as $z_k=k\cdot \Delta z; k=0,1,2,\dots ,N_z$, where τ_j represents the spatial points and N_z is the total number of spatial intervals.

In Equation (7), it is noted that

$$\lambda =\frac{D\cdot \Delta \tau }{{\left(\Delta z\right)}^2},$$

(8)

and it results in

$$C_{z_{i,}{\tau }_{j+1}}=\lambda C_{z_{i-1},{\tau }_j}+\left(1-2\lambda \right)C_{z_{i,}{\tau }_j}+\lambda C_{z_{i+1,}{\tau }_j},$$

(9)

Boundary conditions are set as

$$C_1^{{\tau }_{j+1}}=C_0^0,$$

(10)

Equation (10) allows the calculation of the concentration in a certain node of the grid, i.e., in the node (i, j + 1), which has the abscissa z, at the moment τ + Δτ depending on the concentration values in three adjacent nodes, $C_{z_{i-1},{\tau }_j}$, $C_{z_{i,}{\tau }_j}$, and $C_{z_{i+1,}{\tau }_j}$.

In numerical simulations, a computational grid (or mesh) is used to discretize the domain for solving partial differential equations. For a default scheme using the finite difference method (FDM), the computational grid defines the points at which the solution is approximated.

For the tank within the experimental stand, the computational grid is presented in Figure 3.

The computational grid for the default finite difference scheme involves discretizing the spatial and temporal domains into a grid of points where the differential equations are approximated. The solution progresses iteratively over time, using boundary and initial conditions to ensure the accuracy and stability of the numerical simulation.

Considering the initial concentration of oxygen dissolved in water C₀ = 3.12 mg/L, the concentration at saturation C_s = 8.3 mg/dm³, and the volumetric mass transfer coefficient ak_L = 0.09 s⁻¹ [33,34,35,36], the source (S—the source term represents the rate of oxygen addition or consumption) of dissolved oxygen introduced into the water at each time step can be calculated as

$$S=\frac{\partial C}{\partial \tau }=a{k}_L(C_s-C_0)=0.09\cdot (8.3-3.12)=0.4662 \mathrm{mg}/{\mathrm{dm}}^3,$$

(11)

For “j + 1” time steps, Equation (9) becomes

$$C_{z_{i,}{\tau }_{j+1}}=\lambda C_{z_{i-1},{\tau }_j}+\left(1-2\lambda \right)C_{z_{i,}{\tau }_j}+\lambda C_{z_{i+1,}{\tau }_j}-\left(j+1\right)\cdot S,$$

(12)

To solve Equation (12), the following conditions are set:

$$C_{i,0}=3.12 \mathrm{mg}/{\mathrm{dm}}^3;C_{0,j}=3.12 \mathrm{mg}/{\mathrm{dm}}^3;C_{1,j+1}=3.12 \mathrm{mg}/{\mathrm{dm}}^3,$$

(13)

$$\lambda =\frac{D\cdot \Delta \tau }{{\left(\Delta z\right)}^2}=\frac{0.26\cdot {10}^{-4}\cdot 15}{{(1/50)}^2}=0.975,$$

(14)

$$\begin{array}{c}\begin{array}{l}For i=1, j=0:C_{1,1}=\lambda C_{0,1}+\left(1-2\lambda \right)C_{1,1}+\lambda C_{2,1}-S\to C_{2,1}=3.59 \mathrm{mg}/{\mathrm{dm}}^3\\ For i=2, j=0:C_{2,1}=\lambda C_{1,1}+\left(1-2\lambda \right)C_{2,1}+\lambda C_{3,1}-S\to C_{3,1}=4.05 \mathrm{mg}/{\mathrm{dm}}^3\\ For i=1, j=1:C_{1,2}=\lambda C_{0,2}+\left(1-2\lambda \right)C_{1,2}+\lambda C_{2,2}-2\cdot S\to C_{2,2}=4.07 \mathrm{mg}/{\mathrm{dm}}^3\\ For i=2, j=1:C_{2,2}=\lambda C_{1,2}+\left(1-2\lambda \right)C_{2,2}+\lambda C_{3,2}-2\cdot S\to C_{3,2}=5.97 \mathrm{mg}/{\mathrm{dm}}^3\\ For i=1, j=2:C_{1,3}=\lambda C_{0,3}+\left(1-2\lambda \right)C_{1,3}+\lambda C_{2,3}-3\cdot S\to C_{2,3}=4.55 \mathrm{mg}/{\mathrm{dm}}^3\\ For i=2, j=2:C_{2,3}=\lambda C_{1,3}+\left(1-2\lambda \right)C_{2,3}+\lambda C_{3,3}-3\cdot S\to C_{3,3}=7.41 \mathrm{mg}/{\mathrm{dm}}^3\end{array},\end{array}$$

(15)

By iterating this equation over time steps for each node of the grid, it results in a system of linear equations that can be expressed in matrix form; finally, the variation in dissolved oxygen concentration in water (see Figure 4), in space and time, is obtained.

The initial steep slope of the curve suggests a rapid increase in dissolved oxygen concentration, likely due to high diffusion rates when the concentration gradient between air and water is the greatest. As time progresses, the rate of increase in dissolved oxygen slows down, indicating a decrease in the concentration gradient.

The rotating fine bubble generator was used to oxygenate mains water with the characteristics from

Table 1. The experimental conditions.

Characteristic	Notation	Value
Initial dissolved oxygen concentration	C0	3.12 mg/dm3
Dissolved oxygen concentration at saturation	Cs	8.3 mg/dm3
Water temperature	t	24 °C

:

Considering these conditions, the variation curve of the concentration of dissolved oxygen in water was determined.

In the present model, the size of bubbles has been assumed to be constant and independent of hydrostatic pressure. Future models can be improved by incorporating this dependence, leading to more accurate representations of the physical processes involved. However, for the scope of this study, the assumption of constant bubble size provides a reasonable approximation that simplifies the analysis and focuses on the primary aspects of diffusion.

4. Discussion

The finite difference method provides a straightforward way to discretize and solve the convection–diffusion equation, enabling the computation of the concentration C(z,τ) at successive time steps. The propagation of the solution in the space–time plane is visualized by iteratively updating the concentration values at each grid point based on the discretized equation. This approach allows the tracking of how oxygen concentration evolves over time across the spatial domain.

Specifying the dividing points of the intervals of the variables t and z is essential for clarity and reproducibility in any analysis or research. These dividing points can significantly impact the interpretation and outcome of this study.

Considering the diffusion equation in the non-stationary regime, Equation (3), one can introduce the following dimensionless variables:

Dimensionless time: $\tau =\frac{t}{t_0}$

Dimensionless vertical position: $\varsigma =\frac{z}{L}$

Dimensionless concentration: $\theta \left(\varsigma ,\tau \right)=\frac{C\left(z,t\right)}{C_0}$

Where t₀ is a characteristic time scale, L is the characteristic length scale, and C₀ is a reference concentration.

Rewriting the diffusion equation using these dimensionless variables, one can obtain

$$\frac{\theta \left(\varsigma ,\tau \right)}{\partial \tau }=\frac{t_{0}D}{L^2}\frac{{\partial }^2\theta \left(\varsigma ,\tau \right)}{\partial {\varsigma }^2},$$

(16)

If the dimensionless diffusion coefficient is defined as

$$Fo=\frac{t_{0}D}{L^2},$$

(17)

where Fo is the Fourier number, a dimensionless number that characterizes heat conduction (or diffusion) relative to the storage, the dimensionless diffusion equation becomes

$$\frac{\theta \left(\varsigma ,\tau \right)}{\partial \tau }=Fo\frac{{\partial }^2\theta \left(\varsigma ,\tau \right)}{\partial {\varsigma }^2},$$

(18)

Considering the dimensionless boundary $\theta \left(0,\tau \right)=0; \theta \left(1,\tau \right)=0$ and initial conditions $\theta \left(\varsigma ,0\right)=\frac{f\left(z\right)}{C_0}$, the analytical solution in dimensionless form follows the same steps as the dimensional solution but now uses dimensionless variables.

This approach enhances the generality and applicability of the model, making it easier to compare with other studies and apply it to different systems.

To adapt the model for situations with laminar mixing or a stagnant liquid, it is necessary to modify the assumptions and approach to reflect the lower Reynolds number (Re) regime, where mixing is less intense and diffusion processes become more dominant. In laminar flow, the Reynolds number is low, and in a stagnant liquid, the flow is essentially zero, and molecular diffusion is the primary transport mechanism. So, the model described in the manuscript, with appropriate boundary and initial conditions that must be defined to reflect the specific setup, can be used.

5. Conclusions

The contributions made in this paper include the study of methods for solving differential equations with partial derivatives that are the basis of the description of the oxygen transfer process from air to water; the approach and presentation of the mathematical model for solving the equation that determines the transfer rate of oxygen to water; the classification of methods for solving second-order differential equations; the finite difference method; and solving the differential equation with the term “source”.

The mathematical modeling of the oxygen transfer process from air to water involves setting up a differential equation that captures the key dynamics of the system. Equation (12) represents the change in oxygen concentration over time due to molecular diffusion from the area with high concentrations to the area with low concentrations. The concentration of dissolved oxygen in the water will increase over time if air is blown into the water. Since the differential diffusion equation is a parabolic equation, it was chosen to be solved by finite differences in the default version.

The finite difference method provides a practical and efficient approach to solving the PDEs governing the oxygen transfer process from air to water. By discretizing the spatial and temporal domains and iteratively solving the resulting system of equations, one can model and analyze the dynamics of oxygen concentration under various conditions and configurations.

The approach and presentation of the mathematical model for solving the differential equation of the parabolic type, of the diffusion in the non-stationary regime, in the vertical direction, by the finite difference method, the implicit scheme is a success of the authors and brings benefits to the field.