Chair of Separation Science and Technology, TU Kaiserslautern, PO Box 3049, Germany
Received date: June 28, 2016; Accepted da te: July 07, 2016; Published date: July 17, 2016
Citation: Hlawitschka MW, Drefenstedt S, Bart HJ (2016) Local Analysis of CO2 Chemisorption in a Rectangular Bubble Column Using a Multiphase Euler-Euler CFD Code. J Chem Eng Process Technol 7:300. doi:10.4172/2157-7048.1000300
Copyright: © 2016 Hlawitschka MW, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Visit for more related articles at Journal of Chemical Engineering & Process Technology
Hydrodynamics and reaction process phenomena in bubble column reactors are complex and intrinsic interlinked and not fully understood. To obtain a better understanding, for the first time, a multiphase Euler-Euler solver is developed and validated to experimental data for the absorption of carbon dioxide in water offering the possibility of detailed temporal local and spatial analysis. The open source CFD tool OpenFOAM® with the multiphase solver multiphase Euler Foam is extended: (i) an absorption model was implemented, enabling absorption from different bubble size or phases to the liquid (ii) a chemical reaction model accounts for any number of reactions. The extensions were validated based on literature data and own experimental results at higher pH-value. The simulation results show a satisfactory agreement to experimental and simulative data from literature as well as to own experiments concerning bubble velocity, concentrations and pH profiles
CFD simulation; Bubble column; Chemisorption; Absorption; Reactive flow; Euler-Euler; Multiphase CFD
Reactive bubble columns are widely used in chemical, petrochemical, biochemical industries . Different intrinsic interlinked phenomena such as hydrodynamics, absorption, reactions, coalescence and break-up of bubbles are involved and make predictions and scale-up difficult. Due to the complexity of the hydrodynamics and its mutual dependency on chemical species transport (reactions, interfacial mass transfer), the industrial approach to handle such systems is very often dependent on a ‘rule of thumb’ basis. The problem is either to neglect local flow phenomena near stirrers and dead zones (integral approach) or not to account for the poly-dispersity of the bubble swarm (pseudo-homogeneity) . State of the art models, such as axial dispersion models, are still based on an integral description of the time-dependent hydrodynamics (one axial dispersion coefficient accounts for all non-idealities), mass transfer (integral mass transfer coefficient constant over the column height) and an integral bubble size (such as the constant mean Sauter diameter, d32). Hence, these models are developed for a specific reactor geometry and chemical system and thus limited in predictability.
Computational fluid dynamics (CFD) simulations can exactly predict local hydrodynamics considering any changes in apparatus geometry. In addition, population balance equations became a frequently applied tool to describe the changing bubble size distribution due to break-up and coalescence, while the bottleneck is still the predictability of coalescence and break-up kernels. However, recent research indicated a good agreement when coupling simulations with mass transfer models.
Concerning reactive mass transfer,  used a one-dimensional two-phase simulation combined with population balance simulation to describe the chemisorption of carbon dioxide into an aqueous solution of sodium hydroxide. The simulations reveal a high sensitivity against small variations in hydrodynamic, mass transfer, and kinetic parameters. By Ref.  developed a model for the simulation of hydrodynamics and an irreversible chemical reaction in a gas-lift reactor. A solver capable of calculating compressible two-phase bubbly flows with chemisorption has been introduced by Ref.  and the simulations were compared to experimental data of Ref. . While a good fit was found for the pHprofile, the used geometrical design of the column was different to the experimental data  developed a transient Eulerian-Eulerian model for hydrodynamics and mass transfer in rectangular bubble columns but chemical reactions were not included. By Ref.  introduced a hybrid approach for coupling turbulent mixing and chemical reaction. The mixing of phases is estimated by averaging flow and concentration profiles from preliminary CFD flow field calculations and a numerical tracer experiment. The validation was with data of carbon dioxide absorption into an alkali solution. Darmana et al.  used an Eulerian-Lagrangian approach to model hydrodynamics, mass transfer and chemical reactions for twophase flows. The chemisorption process of carbon dioxide in an aqueous solution of sodium hydroxide and the resulting reversible reactions were simulated in a bubble column reactor. A Lagrangian model enables an exact consideration of dispersed bubbles. However, it is not suitable for higher dispersed phase fractions and in general, this type of approach leads to a higher computational effort compared to an Eulerian one.
In this study reaction algorithm are combined with a multiphase Euler-Euler CFD code based on the open source toolbox OpenFOAM®. To our knowledge, it is the first coupling of a CFD multi-phase model with reaction kinetic models to simulate the reactive mass transfer, while the coupling with population balances is presented elsewhere . The resulting solver is capable of computing one liquid bulk phase and any amount of gaseous phase fractions in 3-D simulations. The model will be described in the next section, followed by a validation against literature test cases and own experimental data derived in a rectangular bubble column.
The multi-fluid model is based on general conservation equations such as the conservation of volume and conservation of momentum for each phase. The continuity equation for the liquid phase l is given by:
Here represents the velocity of the continuous (liquid) phase and is the mass transfer from the liquid phase l to the dispersed (gaseous) phase g and is the mass transfer vice versa. Sl defines the source term applied to the liquid phase entry and αl is the liquid holdup. The first term on the left hand side describes the accumulation, the second term the convection and the third term a possible interface sharpening, while the interface compression velocity is given by :
The conservation of momentum of the liquid phase is represented by:
The viscosity of the liquid phase is given by μl, the gravitation vector is named, and the forcesanddescribe the drag force, lift force and virtual mass force. The term describes the momentum influence due to mass transfer. The lift force is not relevant and not considered but the drag term is given by:
The dispersed phase fraction given by αd and the corresponding velocity is described by . The drag coefficient CD can be described by a various number of correlations. For example, the first one used in this work is Ishii and Zuber :
The second one used is a variant of Tomiyama et al. :
The virtual mass force includes the virtual mass effect occurring when a phase accelerates relative to another phase:
The turbulence modelling is based on large eddy simulation  with a simulation of large-scale turbulence structures and a finetuned modelling of the smaller scales. In this work the latter is by the Smagorinsky model , where the turbulent viscosity is given by:
Where Δ represents the grid size and is the rate-of-strain tensor. In OpenFOAM 2.3.1 the Smagorinsky constant Cs is defined as follows:
The model constant Cs=0.1 is taken from literature  and in ordered to ensure comparability the constants Ck=0.03742 and Ce=1.048 were used throughout.
Mass transfer and species transport
The chemisorption of carbon dioxide in aqueous sodium hydroxide solution starts with the absorption process of carbon dioxide from the dispersed phase into the surrounding continuous phase:
From a computational point, the presence of gas phase in a mesh cell leads to an increase of the corresponding species in the continuous phase until the solubility limit is reached.
The absorption rates appear in the source terms of the species transport equations, which describe the transport of the species inside the numerical domain:
Only Ns-1 transport equations are solved, where Ns corresponds to the number of species observed in the process. The mass fractions of the species must fulfill the following constraint:
The mass transfer is based on the two-film theory with a linear gradient of a species j in the liquid (continuous) phase and the dispersed phase. Figure 1 shows a schematic concentration profile and the mass fractions in the bulk phases and at the interface.and represent the mass fractions in the bulk of the dispersed (gaseous) phase and the surrounding liquid phase, respectively. The mass fractions and describe the corresponding mass fractions at the interface. Based on the assumption of a liquid side mass transfer resistance, the mass transfer is given by:
is the density of the continuous phase, Vcell is the cell volume and E is the enhancement factor.
The interface area corresponds to the summation of the interface areas of each single bubble in a numerical cell. Assuming spherical bubbles of constant size, it corresponds to:
The enhancement factor E describes the influence of chemical reaction to the absorption. Fleischer et al.  analyzed the dependence of the enhancement factor E on the pH value. A correlation is used to approximate the results of Fleischer in dependence of the hydroxide mass fraction as depicted in Figure 2.
Up to a pH-value of 10, the enhancement factor E is close to unity (see horizontal dotted line), but in the range of pH 10 to 14 rises rapidly and describes an enhancement of the physical mass transfer. Albeit its simple structure, the average deviation of the correlation of the findings of Fleischer is about 4% and never higher than 10%.
The mass transfer coefficient Kj is calculated based on the Sherwood Shj number and the diffusion coefficient Dj.
corresponds to the equilibrium mass fraction of species j in the liquid phase. It is defined by the Henry constant H.
With the above equations the mass transfer from the dispersed to the continuous phase is then:
Reactions describe the transformation of one set of chemical substances into another. The developed solver enables a description of reactions in the liquid phase. The general stoichiometry equation of a reversible reaction for J species taking part in a reaction is:
Xj represents the summation formula of species andare the stoichiometric coefficients of the reaction m for the educt side (left) and product side (right). For a chemical system with M reactions the production rate of a species is:
The production rates are used as source terms in the species transport equation. The reaction velocity m ω of the m-th reaction is dependent on the rate coefficients and for the forward chemical reaction and the backward chemical reaction as well as on the concentrations cj of the participating species j.
The concentration cj of the species j are derived from the mass fractions Yj, the density ρl and the molar mass Wj.
The temperature dependent rate coefficients of the forward reaction is described by the Arrhenius equation..
Am is the pre-exponential factor. The activation temperature TA can be interpreted as thermal energy required to start the reaction. The rate constant of the backward reaction is commonly derived from the equilibrium constant GG, m K :
For the rate constants and the equilibrium constant KGG,m different empirical models available (s. chapter 4.3). Besides the existing models in OpenFOAM 2.3.1 an additional model for the rate constant was implemented:
Where Am, Bm, Cm, Dm and Em are reaction specific constants.
The model is embedded in the standard multiphaseEulerFoam solver in OpenFOAM 2.3.1 and enables the simulation of reactive mass transfer. Declarations of the different phases i as well as hydrodynamic calculation of these are done within the original multiphaseEulerFoam environment. Figure 3 shows the interaction between the phases i, the chemical species j and the physical phenomena where one continuous phase with an infinite number of gaseous phases or different bubble sizes can be coupled. The composition of the liquid phase is being described by the OpenFOAM class basic Multi Component Mixture, where scalar fields Yj represent the mass fractions of the chemical species. It has no influence on the hydrodynamic behavior of the phases thus neglecting the last term in Eq. (4). Any modification induced by irreversible, consecutive, parallel or reverse reactions may be defined. A change in composition of the gaseous phases is not specified as it is always a pure gas. Figure 4 shows the computational sequence of the developed solver for a single time step. Newly introduced extensions to the original solver are highlighted in grey.
The new model is validated against literature data in three steps: First, the predicted hydrodynamics are compared with the results of Deen et al. . The work of Ref.  is used to verify the capabilities of the absorption model and then the chemistry model. The latter is further validated through own experiments. All simulations are performed in rectangular columns as shown in Figure 5. The continuous aqueous phase is not being exchanged in the semi-batch reactors. Our experiments reveal that temperature changes occur in a range of 1-2 K, which is in accordance to literature . This small change is neglected and all simulations are performed for 300 K and 1 bar. The utilized discretization schemes and boundary conditions are identical in all simulations. They are shown in Tables 1 and 2. The virtual mass coefficient CVM is set to be 0.5.
|ddt Schemes||fixed Value|
|grad Schemes||Gauss linear|
|div Schemes||Gauss upwind|
|Laplacian Schemes||Gauss linear corrected|
Table 1: Utilized discretization schemes for the occurring terms in differential equations.
|Phase fraction α||fixed Value||inlet Outlet||zero Gradient|
|Mass fraction of
chemical species Y
|zero Gradient||inlet Outlet||zero Gradient|
|Pressure p||fixed Flux Pressure||fixed Value||fixed Flux Pressure|
|Temperature T||zero Gradient||inlet Outlet||zero Gradient|
|Velocity||fixed Value||pressure Inlet-Outlet Velocity||fixed Value|
Table 2: Utilized boundary conditions for different physical variables.
Validation of hydrodynamics
Hydrodynamics play a crucial role when it comes to modelling of reactive mass transfer as it determines the local concentration changes in the liquid and gas phase. An initial simulation of a rectangular semi-batch bubble column (0.15 m × 0.15 m × 0.60 m) is based on the publication of Ref.  with a numerical mesh constructed by 30 × 30 × 120 (width × depth × height) hexahedrons. The superficial gas velocity is set to 4.9 mm/s. Gaseous carbon dioxide enters through a square pitch of 30 mm side length positioned at the center of the base area. The bubble diameter is set to the constant value of d=4 mm. At the start of the simulation the filling level is at a height of 0.45 m. During the simulation, the time step is adjusted to match a Courant number of 0.3. Turbulence is represented by the Smagorinsky model with a parameter of CS=0.1. Drag forces are calculated according to the model of Ishii-Zuber and mass transfer is neglected.
The Figure 6 shows the time averaged axial velocity of the continuous water phase in m/s plotted over the normalized side length x/X of the reactor. Velocities are measured at the centerline of the reactor at a height of 0.25 m. The simulative values are averaged over the entire duration of the calculation and the results are in good agreement with literature data. On the left hand side velocities are underestimated in the simulations compared to the experiment, whereas on the right side values are over estimated. Simulated values of the literature and own results deviate especially on the left side of the column, in the area of 0.1-0.4 x/X. The maximum mean velocity calculated by Ref.  is 0.186 m/s whereas our simulation produces a value of 0.176 m/s, which accounts to a deviation of 5%. Fluid phase velocity is being fixated by the boundary conditions to 0 m/s at the reactor walls. Concerning the symmetry of the flow profiles, the performed CFD simulation shows a high symmetry, while literature values slightly deviate.
Validation of physical absorption
The performed simulation is based on the publication of Ref.  and the setup is identical to the previous simulations. In addition to hydrodynamics the physical absorption of carbon dioxide is calculated using a Sherwood correlation for moving spheres :
Here initially the terminal rise velocity of the CO2-bubbles and in preceding simulation a value of =0.231m/s in accordance with the literature  is used. For the later the Sherwood number is being calculated in each cell for every time step according to the relative velocity between the gaseous and liquid phases.
In both simulations a rate of 2E-9 m²/s is used for the diffusion of CO2 in water. The solubility of CO2 in pure water is approximated with the formula from Versteeg and van Swaaij :
Darmana et al.  estimates the evolution of mass fraction of diluted gas according to:
Ag,tot is the total surface area of all gas bubbles inside the reactor. Vl,tot is the volume of the bulk phase. The assumptions of constant bubble size and quantity are made. The enhancement factor in Eqn. 21 was set to 1. Figure 7 compares simulations of the physical absorption process at global constant and a local (variable) Sh number with literature  and an absorption estimation according Eq. (32). The normalized concentration of diluted carbon dioxide in the batch liquid phase is shown in dependence of simulation time t in seconds. Values are taken at a height of 0.225 m in the center line of the reactor. The own CFD simulations as well as the simulation by Darmana et al.  show a final deviation of 2% compared to the absorption estimation. Also, the simulation with a constant Sh number (Sh=437) based on a relative gas velocity of =0.231 m/s does not show a huge deviation to the simulation using a variable Sh number (Sh=75-850) based on the local bubble velocity. A detailed study of the Sh number revealed, that the Sh ranges from 437 ± 5 in the area, with high gas fraction and only shows a higher deviation in areas with low gas fraction. For other geometrical designs and disperser designs, the significance of a local Sh number calculation may be more distinct.
Comparison to literature data: The chemisorption of carbon dioxide in aqueous sodium hydroxide solution is simulated with respect to the work of Ref.  at reactor measures of 0.2 × 0.03 × 1.5 m3. Initially the level of the liquid batch phase has a height of 1 m and a pH of 12.5. All bubbles have a constant diameter of 5.5 mm. The mesh grid consists of 27 × 4 × 200 (width × depth × height) regular hexahedrons. Gas is being introduced with a superficial velocity of 7 mm/s and the central inlet on the bottom has the dimensions 35 × 15 mm2. As above, the Smagorinsky model with a coefficient of 0.1 is used for modelling turbulence. For both phases the drag coefficient is being calculated with the formula
According to Ref.  a constant Sherwood number of 562 is chosen to model mass transport. The diffusion rate of carbon dioxide in water is calculated with the help of the formula from Versteeg and Swaalj .
In deviation to literature the enhancement factor E is initially constant and set to be 1. Two reversible chemical reactions take place in the reactor:
CO2 +OH− HCO3− (35)
The production rates result to:
The derivation of the reaction rate constants k follows literature. The constant of the first forward reaction is estimated with the relation of Pohorecki and Moniuk .
Eigen et al.  concluded that the reaction rate constants of processes with proton transfer have a magnitude of 1010-1011. The rate of the second reaction is set to 106 for reasons of CPU time. The constant for the backward reaction is again calculated with the help of the equilibrium constant.
For this purpose the approximation of Hikita et al.  is used:
Previous to the actual simulation of chemisorption, 120 s of pure hydrodynamic movement had been simulated. This prelude enabled the bubble plume to fully establish itself. After activating the absorption and reaction models another 250 s of simulation time were calculated.
In Figure 8 shows the evolution of species concentrations during the conducted simulation of chemisorption and the results of Ref. . The concentrations in kmol/m3 are hereby plotted against time t in seconds. During the first 80 s the production of CO32- is favored. Then the equilibrium changes towards the production of HCO3- due to a drop of the pH-value. After about 100 s all of the hydroxide-molecules are consumed and CO2 starts to accumulate. The simulation is in very good agreement with the values from literature. A time delay of approximately 10 s can be observed between the two calculations. It is being assumed that this is a result of the constant diameter model used in this work as well as the neglect of the enhancement factor E.
A shift of equilibrium of the first reaction happens at about 80 s of simulation time leading to the formation of bicarbonate. After approximately 200 s the concentration of bicarbonate reaches its maximum level as the accumulation of carbon dioxide inside the reactor begins. Both events can be observed as shifts in the evolution of the pH-value. In Figure 9 shows the temporal evolution of the pHvalue observed in both the experiment and simulation performed by Ref.  as well as own simulation results. All three curves decrease over time from their initial value of 12.5 to about 6.9 pH after 250 s. And all of them show the shifts mentioned above. But the simulative evolutions are delayed compared to the experiment, which is more serious for the second bump. However, the results of Ref.  are in slightly better agreement with the experiment than our own simulations. Yet both simulations capture the general trend of the experimental data. Darmana et al.  explain the delay with an underprediction of the mass transfer due to the chosen model. The time delay between the calculation in literature and the own results can be explained in the same manner as above when analyzing the evolutions of the species concentration.
The simulation is improved when considering a variable enhancement factor E using the above introduced Eqn. 18. The result of the pH-evolution with variable enhancement factor (—) is compared in Figure 10 to the data by Darmana et al.  (- - -) and our results with constant enhancement factor (····). The general behavior of the pH-value is similar to the previous simulations in Figure 9. The pH decreases over time having a value of about 6.9 after 250 s. Also the two mentioned pH value steps, which refer to changes in chemical equilibrium, can be observed. In contrast to the previous simulation the first step is now predicted slightly prematurely and the second step is again too late. However, it can be concluded that the use of a variable enhancement factor according Eqn. 18 greatly improves the agreement with the experimental data.
Comparison of pH-value to own experimental results: The capab ility of the solver to predict evolutions of the pH-value is further tested in comparison to own experiments. Again the chemisorption of carbon dioxide in aqueous sodium hydroxide solution is taken as reference case. The reactor consists of paraglas and measures 0.18 × 0.03 × 1.5 m3. Initially the level of the liquid phase has a height of 1 m and the solution has a pH of 13. The gas is introduced into the reactor through 3 × 7 needles in rectangular position at the center of the bottom plate. The superficial velocity is held constant at 8 mm/s, resulting in a bubble size of approximated 5.5 mm. Every 30 s the pH value is determined with the help of the portable measuring device SP80PI VWR sympHony. Before the measurements, the column is gassed with nitrogen and to start chemisorption, the gas feed is instantaneously shifted by a valve to carbon dioxide. Initially both the liquid and gaseous phase are at room temperature. During the experiment a rise in temperature from 294.5 K to 296.6 K is being measured.
A simulation is performed accordingly to the first case (chapter 4.3.1) with adapted boundary conditions as inlet velocity and mesh geometry. The mesh consists of 27 × 4 × 200 (width × depth ×height) regular hexahedrons resulting in a total cell number of 21.600 cells. As in the first case and corresponding to the injection of N2 in the experiments, pure hydrodynamic movement is simulated till 120 s as a starting point for the CO2 injection.
Figure 11 shows the temporal evolution of the pH-value observed in the performed experiment and simulation. The pH values decrease over time from their initial value of 13 to approximately 10 pH after 250 s. Compared to the first presented case, the higher starting pHvalue leads to a change of the temporal reaction progress. While OH- and CO32- were totally converted after 250 s in the first case, in this second case, this reaction is in the middle of the process.
While the simulation slightly underestimates the experiment after 120 s, the deviation is never higher than 3%. Hence, for this period of reaction, the simulation results are in very good agreement with the experimental ones. The overall computational time on 10 cores for a simulation of 250 s is 72 h considering a mesh size of 21600 cells.
A new CFD multiphase solver enabling reactive mass transfer has been developed. The solver enables for the first time a simulation of different bubble sizes or different gases in a single computational cell using Euler-Euler framework. The solver has been applied to study reactive mass transfer in bubble column. The results were compared to experimental and simulative data from literature  as well as own experiments. At no point, model parameters were fitted to match results with literature. The correct prediction of the hydrodynamic behavior of the phases was shown in comparison to the results of Ref. . It could be shown, that the evolution of concentration changes due to absorption and chemical reaction is in very good agreement with the simulative data from Ref.  for an initial pH value of 12.5, who used Euler-Lagrangian framework. A small time delay between own results and literature is being attributed to differences in the models for bubble size and the enhancement factor E. It could be shown, that accounting a variable enhancement factor in the simulation improves significantly the agreement with experimental results. The first time comparison to a second a second case experimental setup, showed a temporal delay in the overall reaction process, which is related to the higher pH-starting value. In conclusion, a very good agreement between simulation and experiment can be observed without any further adjustment, with exception of the initial starting values that must fit to the experimental setup.
For an independent CFD based column design, the code requires further validation, for example, at higher gas fractions. The use of different gas phases in a gas mixture is already possible but not sufficiently tested. Momentarily, the assumption of spherical bubbles of constant size is being used and requires an adaption to account for changes in the interfacial area due to (elliptical) shape variations of the bubbles. The introduction of an adequate PBM bubble model by combining the developed solver with the work of Ref.  enables more sophisticated simulations for heterogeneous flow conditions. With that, bubble coalescence and break-up as well as different and changing bubble shapes can be modelled. This will affect the absorption via the bubble surface area and thus reactive mass transfer. The composition of the gas phases should be specified just as it is being done for the liquid phase to account for a back diffusion of dissolved gas into the bubbles. This would allow simulating the mixing of gases or the absorption of vaporized liquid into a gas phase. Also, the gasliquid interface requires additional investigations. For example, the coalescence of bubbles induces a locally higher mass transfer, which is not accounted in the simulations till now. In order to correctly predict bubble movement, the lift force needs to be considered. Several models have already been developed and tested in literature . Additionally, the code can be extended to allow the simulation of other industrially relevant combinations of systems: solid-liquid, solid-gas and liquidliquid, as the multiphase approach is not restricted to gas fractions. An example to this are liquid-liquid systems, such as (reactive) extraction, when using OpenFOAM with respect to this .
This work was supported by the German Research Foundation (DFG), Priority Program SPP1740 “Reactive Bubbly Flows” (DFG HL-67/1-1).
|dc||phase material time derivative for the continuous phase|
|dd||phase material time derivative for the dispersed phase|
|gravitational acceleration||kg/m s²|
|k||mass transfer coefficient||m/s|
|k||reaction rate constant||varies|
|mass transfer||kg/m³ s|
|Cα||interface sharpening constant|
|H||dimensionless Henry constant|
|K||drag exchange coefficient||kg/m³ s|
|S||source term||kg/m³ s|
|mean rate of strain tensor|
|α||volumetric phase fraction||m³/m³|
|µ||dynamic viscosity||kg/m s|
|µsgs||turbulent viscosity||kg/m s|
|ω||reaction velocity||kmol/m³ s|