Investigating the Advantages and Limitations of Modeling Physical Mass Transfer of CO 2 on Flat Plate by One Fluid Formulation in OpenFOAM

One fluid formulation is an approach used for modeling and analysis of mass transfer between two immiscible phases. In this study we implement and analyze the advantages and limitations of this approach for CO 2 physical mass transfer into MEA. The domain is a flat plate and gas liquid flow is counter current. The analysis was carried for operating parameters like liquid phase Reynolds number, MEA mass fraction and the angle of inclination of flat plate. The results clearly show that the model effectively captures the deviation in liquid side mass transfer coefficient due to the surface instabilities and liquid properties which are generally neglected by standard correlations. Also the model shows that the standard Higbie correlation is preferable at low Reynolds number at any angle of inclination. The grid independent studies show that a size of 6.25 µm is required in the interface region for effectively using this approach. The computational resource time at this resolution was found as the only limitation for using this approach and we suggest a procedure to overcome this limitation. The present simulation results can help CFD researchers investigating immiscible gas-liquid mass using OpenFOAM.

1 Introduction CO 2 is the main greenhouse gas that is causing an increase in average surface temperature of earth and need to be captured at the source to reduce its impact on the environment. The main sources of CO 2 emission are power plants, industrial processes and domestic consumption of fuels. Nearly 40 % of this total CO 2 emission originates from industries [1]. Hence, technologies required to capture the CO 2 gas at the source need to be developed. Several technologies have been developed in this regard like membrane separations (utilizing a membrane to separate the CO 2 from flue gases) [2], adsorption (adsorbing the CO 2 on the surface of the chemically activated solids) [3] and using solvents like ionic liquids [4], amines blends [5] (with and without reactions). Among all of them the absorption of CO 2 by MEA (mono ethyl amine) solvent is one of the most preferred technique [6]. CO 2 absorption into MEA is majorly carried using packed bed reactors. Various kinds of reactors like bubble column, stirred tank, packed bed reactors, fluidized bed reactor, membrane reactors and hybrid reactors combining above reactor types [7] are available for effective utilization of absorption phenomena. Among them, the availability of high interfacial area for mass transfer at low pressure drop makes the packed bed reactors preferable, commercially tested at industrial scale [8], for gas liquid absorption process.
The packed bed reactors facilitate high interfacial areas using either structured packings or random packings. Structured packings are found to give higher mass transfer rates [9][10][11][12]. The efficiency of a structured packing depends majorly upon micro level interactions between gas liquid immiscible phases next to packing surface. Hence understanding the microlevel mass transfer helps in developing better design and effective utilization of packed bed reactors [13,14].
Computational fluid dynamics (CFD) is an effective tool in understanding the dynamics of multi-phase flows in complex geometries. CFD reduces the overall number of experiments and allow access to velocity, temperature and other scalar and vector fields of interest at any location of the geometry, which otherwise is very difficult to be accessed by experimental techniques [15]. A detailed overview of utilizing CFD as an effective tool in modeling absorption in packed bed can be found the recent article published by Haroun and Raynal [16].
The micro level analysis of structured packing essentially involves modeling of multiphase fluid dynamics as immiscible phase flow where the fluids are not interpenetrating. The immiscible two phase flow is generally modeled using Volume of Fluid (VOF) approach and results are comparable to experimental values derived using particle image velocimetry techniques [17]. The mass transfer across these immiscible fluids is modeled majorly in two different ways. One approach involves the modeling of mass transfer across the immiscible phases by explicitly adding source terms to the governing equations of VOF approach. Sebastia-Saez et al. [15,18,19] modeled physical and reactive mass transfer on flat and texture plates by this approach. The geometry in their work was same as that of Hoffmann et al. [17]. A similar approach was used by Xu et al. [20] for modeling the mass transfer of propane gas into toluene liquid.
Another approach for modeling the mass transfer across immiscible phases is familiarly known as "one fluid approach" or "continuous specie transfer (CST)" approach. This approach was proposed and validated by Haroun et al. [21]. The model has an advantage of avoiding any kind of assumptions in terms of input parameters like mass transfer coefficient values derived from standard correlations. This model was also used in detailed numerical simulation of (direct numerical simulation, DNS) structured packing sheet also to derive mass transfer and liquid hold up in structured packing element by same group. The group of Haroun used JADIM multiphase software developed by IMFT [22]. Nieves-Remacha et al. [23] implemented this one fluid formulation for simulation of mass transfer in an industrial advanced flow reactor. A similar approach was developed by Marschall et al. [24] and was implemented in OpenFOAM [25,26] for mass transfer in gas liquid bubble flow. In Nieves-Remacha [27] dissertation thesis, the two formulations of Haroun et al. [21] and Marschall et al. [24] were compared and found to be producing same steady state results for various simple test cases. Recently, Wang et al. [28] used this approach for simulating the gas liquid mass transfer in wetted wall column. They found the CFD results in reasonable agreement with the experimental values.
The aim of this work is to discuss the advantages and limitations of utilizing one fluid formulation for modeling the physical mass transfer of CO 2 in OpenFOAM CFD software. The one fluid formulation was implemented as an additional transport equation in the existing code of OpenFOAM CFD software and a separate procedure for reducing the overall computational time is proposed. The physical absorption of CO 2 was studied using N 2 O analogy. A comparison of liquid side mass transfer coefficient value derived from the simulations with the standard theoretical correlation proposed by Higbie [29] is carried to study the aim. The investigations were carried for operating parameters such as flow rate, concentration of MEA and angle of inclination of plate. The simulation domain is two dimensional (2D) and the gas-liquid flow is counter current.

Two phase flow modeling
The simulations are carried using "interFoam" module available as part of open source software OpenFOAM. The interFoam is based on color function volume of fluid method and is used in modeling incompressible two phase flows. The implementation involves description of each phase on fixed mesh using Eulerian approach, Eq. (1) and modeling interface using a color function. The color function used in OpenFOAM is volume fraction of secondary phase (α) and its value is derived by solving an additional transport, Eq. (2). The advection of volume fraction is modeled using multidimensional universal limiter with explicit solution (MULES) algorithm. .
By solving these equations , the values of velocity vector, and pressure and phase fraction values at each grid cell is derived at each time step. The cells fully occupied with secondary phase will have a value of α = 1 and the cells fully occupied by primary phase will have value of 0 < α < 1 depending upon the flow conditions. The VOF approach involves solving of a single momentum equation for both the phases with volume averaged properties of at the interface defined by: density, Eq. (3) and viscosity, Eq. (4): Where the subscripts "L" and "G" stands for liquid and gas phase respectively. On similar lines, a single continuity Eq. (5) is solved for mass conservation: The last term in Eq. (1), F st represent surface tension force and is modeled by using formulation proposed by Brackbill et al. [30], Eq. (6): Where σ st is the surface tension coefficient; κ = −∇.  n represents the surface curvature where  n is unit interface normal vector. This normal vector is given by Eq. (7): This normal vector is corrected at the wall boundaries to consider the effect of wall adhesion by relating with the liquid contact angle by Eq. (8): Where  n w and  t w are the normal and tangential vectors to the wall and θ w is the contact angle between liquid and wall.
In interFoam formulation an additional interface compression term is added to the Eq. (2) in order to reduce the numerical diffusion at the interface and to have a sharp interface between the two immiscible phases. The modified Eq. (2) is given by Eq. (9): where  v r is the relative velocity between the phases.
The modeling of this additional term and also the impact of various input values to this additional term can be found in detail in Deshpande et al. [31].

Mass transfer modeling
The one fluid formulation developed by Haroun et al. [21] was implemented for modeling the mass transfer between immiscible phases. This formulation involves modeling the effect of thermodynamic properties like diffusion and solubility within a scalar transport equation instead of using values derived from correlations. The modified additional equation for the specie transport using this approach is given by Eq. (10): The one fluid formulation solves a single transport, Eq. (10) for specie concentration across the two phases by treating the concentration as a function of phase fraction. The VOF involves the same approach for defining the flow properties like density and viscosity. The effect of concentration jump at the interface is included using the additional mass flux term Φ. This additional mass flux is due to the solubility of gaseous specie into the liquid and thus can be calculated using Henry's law. The flux term Φ modeled to include the solubility is given by the Eq. (11): where D stands for the effective diffusivity of the specie in the two-phase mixture. The D is calculated as a harmonic average of diffusion of specie in each phase and is given by Eq. (12): where the D L and D G represent the diffusion coefficient of specie in liquid and gas phase respectively. The H in Eq. (11) represent the dimensionless Henry's constant for the specie under consideration. It is defined as the ratio of concentration of specie in gas ( C G ) to the concentration of specie in liquid ( C L ) which is H = C G / C L . Note that the additional flux term Φ will have numerical values only at the interface and will be zero within the individual phases. This is because the ∇α L term in Eq. (11) is essentially zero in grid cells fully occupied by individual phase fraction.

Input parameters
The physical mass transfer of CO 2 was modeled using surrogate N 2 O gas. This kind of approach was used in many studies involving physical mass transfer of CO 2 and is known as N 2 O analogy. The N 2 O analogy avoids the effect of reaction, between CO 2 and MEA, on absorption. An excellent review of using N 2 O analogy in the context of CO 2 capturing analysis can be found in Monteiro and Svendsen [32]. The simulations were carried at isothermal conditions at a temperature of 298 K. The CO 2 loading on MEA solvent for current investigation was fixed at 20 %. It was chosen so that the current studies can be used as a preliminary work for reactive mass transfer. In reactive mass transfer, the value of CO 2 concentration dissolved in MEA will highly influence the equilibrium partial pressure on the gas side of the interface and for CO 2 loading less than 30 % this influence was found negligible by Aronu et al. [33] and Wang et al. [28]. The solvent properties are chosen based on the concentration of MEA (% wt) which in our study vary between 10 to 40 %.
In the current simulation, pure N 2 O gas was used on gas side. The gas density values were calculated using ideal gas law and viscosity values were taken from open literature for pure N 2 O. The N 2 O diffusivity within gas medium is taken equal to that of kinematic viscosity. The diffusivity value of N 2 O in the solvent is calculated from correlation proposed by Versteeg and van Swaaij [34] as shown in Eq. (13).
where the µ water and µ soln are viscosities of pure water and carbonated solvents respectively. Properties like viscosity and density of carbonated solvents were taken from Weiland et al. [35] and surface tension values were taken from Fu et al. [36]. The Henry constant values were calculated from the Eq. (14)- (16) proposed by Penttilä et al. [37].
Where R is the universal gas constant with value 8.314 Pa.m 3 .mole -1 .K -1 and T is temperature in K.
The contact angle value for solvent was taken as fixed at 40° as given in [28].
The ranges of values used in our simulations are as tabulated in Table 1. Units of all the variables used in the current simulations have SI units.

Geometry, mesh and boundary conditions
The current study was carried on 2D (two dimensional) geometry. The N 2 O gas is the primary phase and the MEA solvent is the secondary phase and the flow was modeled by interFOAM. The geometry under simulation is shown in Fig. 1. The liquid enters the computational domain through a 1 mm inlet and leaves through a 1 mm outlet. The gas enters counter currently through a 4 mm inlet and leaves through a 4 mm outlet.
To avoid the gas phase entering the computational domain from the gas outlet, which is possible when the gas phase near the liquid inlet attains a higher relative velocity (because of high liquid velocity compared to gas) into  Fig. 1 Overall geometry under simulation the domain, the gas chamber size in the current simulation was considered 10 % higher (55 mm) when compared to liquid plate height (50 mm). Generic boundary conditions available in OpenFOAM were used in this simulation and are listed in Table 2. A comprehensive work for setting up a two dimensional investigation in OpenFOAM using generic boundary conditions can be found in [38], for further reading.
In their investigation of counter current gas liquid mass transfer Xu et al. [20] found that to capture the mass transfer effect with reasonable accuracy a mesh size of 0.07 h, where h is the film thickness, is required. Since the average film thickness in our study is 0.4 mm the minimum size of mesh required was 0.028 mm. But through our preliminary investigations to find the grid independency, we found that a mesh size of 0.00625 mm (approx. 0.015 h) is required in order to capture the mass transfer effectively. The requirement of higher mesh refinement can be due to the one fluid formulation used in our study which involves modeling of the very low diffusivity of gas in the liquid solvent. The lower diffusivity confines the gas concentration layer to regions near interface. A similar observation of high mesh refinement was made by Cooke [39]. In the remaining gas side geometry, a uniform mesh size of 0.0125 mm in the horizontal direction was used. Further refinement in this region wasn't affecting the final solution of concentration field, which can be due to the higher diffusivity of gas. In order to reduce the computational time, the size of 0.0125 mm was used in the gas side of domain. In the direction of height, a size of 0.1 mm was used. A snapshot of refined mesh is given in Fig. 2. The final mesh size was approximately 0.3 million hexahedra cells. The details of grid independency study can be found in the Appendix 1.

Solver setting and analysis of results
Higher order vanLeer scheme available in OpenFOAM was used to solve the continuity, momentum and concentration equations. The tolerance of 10 -14 was used as convergence criteria. All the flow equations were solved transiently until steady state flow field was achieved. An adjustable time step was used with a maximum Courant number of 0.9. In order to reduce the simulation time, the flow field was solved initially on a coarse mesh with 0.028 mm in the film thickness region. Later this flow field was mapped on refined mesh with 0.00625 mm and further simulated transiently until steady state flow field was achieved. The concentration equation was then solved independently on the resulting flow field with a fixed time step of 10 -5 s, until steady state was achieved. This procedure has resulted in reducing the overall computational time from several days to hours without compromising the "quality of the model", that is the final two phase flow field will be same as the flow field obtained by refined mesh. Note that this suggested procedure is suitable for problems involving steady state solutions and may not be suitable for obtaining transient solutions which are sensitive to initial conditions. Despite the suggested procedure to reduce the computation time, the time to obtain result for one flow condition that is from solving flow field on coarse mesh to solving concentration Eq. (10) was not less than 30 hours in total on a Dell workstation (DELL PRECISION T7810 (v2)) with 2 Intel Xeon E5-2630 v3 processors (16 core).
On obtaining the concentration values at steady state, the liquid side mass transfer coefficient is calculated from the simulation results using the Eq. (18) where F is the mass flux at the gas-liquid interface per unit area and ∆C is the difference between the concentration at the interface and bulk liquid. The results are compared with theoretical correlation proposed by Higbie [29] as Eq. (19) k D Where τ is the time of exposure calculated using an expression in Eq. (20) reported by Haroun et al. [40]: Where L is the length of the plate and v i is the interfacial velocity.
The interfacial velocity is calculated using Nusselt [41] theory which is valid for laminar flows. The Nusselt expression (Eq. (21)) predicts the liquid velocity within the liquid film region as a function of distance from the plate surface. The expression also includes the effect of inclination of the plate with respect to the horizontal surface.
6 Results and discussion 6.1 Comparing the effect of flow rate The counter current gas flow reduces and flattens the liquid surface velocity and the maximum liquid velocity can occur at a different position than that of interface. This kind of flow filed instabilities at micro scale near the interface can enhance the mass transfer by multitude, even in laminar flow conditions [42]. Hence a study of influence of liquid flow rate is of prime importance in understanding the mass transfer under counter current flow conditions. In this study, the effect of liquid flow rate on absorption of gas was studied at gas phase Reynolds number of Re G = 200 and for liquid phase Reynolds number of 83 < Re L < 249. This Re G was the minimum value used in Yu et al. [42] experiments of counter current gas liquid absorption. Also, it was found in the Yu et al. [42] experiments that the gas phase velocity has little influence on the interface liquid velocity and hence the effect of varying gas phase velocity wasn't studied in our investigations. The liquid and gas Reynolds numbers are defined as: The liquid side mass transfer coefficient, k L ( with units mole / Pa.s.m 2 ) is calculated using Eq. (18) where the flux value F is calculated based on overall mass consumption at steady state. The mass conservation inherently states that the amount of N 2 O dissolved through gas liquid interface is equal to the amount of N 2 O removed by solvent at liquid outlet. This amount of N 2 O removed at outlet when divided by interfacial area gives the mass flux value, F, in mole / m 2 -s. The resulting liquid side mass transfer coefficient is then compared to Higbie [29] correlation for various Re L and results are presented in Fig. 3.
The results agree very well with Higbie penetration theory at low Re L and deviate at higher Re L . The deviation at higher Re L can be due to the influence of surface instabilities whose effects aren't included in the Higbie [29] model as mentioned in Yu et al. [42]. The surface instabilities will lead to deviation in the velocity from standard Nusselt [41] laminar profile. The deviations can be observed from the velocity profile of the liquid phase at the liquid outlet as presented in the Fig. 4(a)-(d).
Also, it can be due to the small size of geometry and parameters used from experimental data published in open literature. Since the current model is promising by accurately predicting the mass transfer coefficient values at low Re L , it may be used by coupling with studies involving DNS to develop new correlations for higher Re L and this will be part of our future work.

Comparing the influence of MEA concentration
Since at Re L of 125 the Higbie [29] model was accurately approximating the liquid side mass transfer coefficient for solvent with MEA mass fraction of 30 %, as explained in the previous section, the current study was carried at Re L of 125. The influence of MEA mass fraction in the liquid solvent on the liquid side mass transfer coefficient has Fig. 3 Liquid side mass transfer coefficient is compared to Higbie [29] correlation for varying Re L been studied for four values of: 10 %, 20 %, 30 % and 40 % respectively. The results are as shown in Fig. 5.
For MEA mass fraction less than 30 % the Higbie [29] model over estimates and for mass fraction of 40 % underestimates the liquid side mass transfer coefficient. Also for MEA mass fraction of 30 % the Higbie [29] model accurately estimates the liquid side mass transfer coefficient. This kind of behavior was observed by Sebastia-Saez et al. [18]. Such variations occurs presumably due to the variation of liquid properties like kinematic viscosity and density of the solvent based on the MEA mass fraction. Hence it is safe to say that, using this approach a better estimate of physical mass transfer can be obtained than theoretical models. Since the MEA mass fraction is an important parameter influencing the chemical absorption rate [43] the current approach is suggested for estimating the enhancement factor and studies involving the design of CO 2 reactors.

Comparing the influence of angle of inclination of plate to the horizontal plane
In the two previous sections the film was flowing downward on a vertical plate. But in practice corrugated sheets are used for designing packed beds. These packed beds are generally made of metal sheets of steel and have corrugated textures with an angle of 45° or 60° between the corrugations respectively [13]. In this study, the influence of the plate inclination on the liquid side mass transfer coefficient was also investigated for inclinations of 45° and 60° at Re L of 83 for solvent with MEA mass fraction of 30 %. Fig. 6 shows simulation domain of an inclined plate in our simulation. The results are shown in Fig. 7. The liquid side mass transfer coefficient is exactly approximated by Higbie [29] model. Hence it is concluded that the current model allows the modeling of physical mass transfer in counter current flow on corrugated sheets, which are used to increase the available surface area per unit volume of the packed bed, accurately.

Conclusions
In this study we implement and analyze the advantages and limitations of one fluid formulation approach for CO 2 physical mass transfer into MEA. The domain considered is a flat plate and gas liquid flow is counter current. The analysis was carried for operating parameters like liquid phase Reynolds number in the range of 83 < Re L < 249, MEA mass fraction in the range of 0.1 to 0.4 and for the angle of inclination of flat plate varying between 45° to 90° respectively.
The results clearly show that the model effectively captures the deviation in liquid side mass transfer coefficient due to the surface instabilities which are significant at higher Reynolds numbers. The effect of liquid properties which vary with the mass fraction of MEA in the solvent are also predicted with greater accuracy. These effects are generally neglected in the standard correlations. The model also shows that the standard Higbie correlation is well suitable for modeling the flows at low Reynolds number.
The grid independent studies show that a size of 6.25 µm is required in the interface region for effectively using this approach. The requirement of this high mesh resolution results in high computational resource time. In order to reduce the overall computational time we adopted a sequential procedure of solving the concentration equation independently on the desired flow field. Also the time required for obtaining the flow field was further reduced by first solving the flow field on a coarse mesh and mapping the result on refined a mesh. This procedure was found to reduce the overall computational time from days to hours. However it should be noted that the suggested procedure is suitable for problems involving steady state solutions and may not be suitable for obtaining transient solutions which are highly dependent on initial conditions.
In conclusion, it can be said that the CFD modeling of mass transfer by one fluid formulation is proven to be a promising approach and an alternative to experimentation. The higher accuracy of this approach is due to the consideration of the effect of thermodynamic properties like diffusion and solubility on mass transfer coefficient instead of using derived values from correlations. Hence future CFD investigations of micro structure impact on CO 2 mass transfer can be carried using this approach and our suggested procedure can be adopted for reducing the simulation time effectively.

Acknowledgement
This study was financially supported by the Anadolu University Research Fund (Project No: 1603F123).   v r relative velocity between the phases (m/s)

Appendix 1
The grid independence study was aimed to find out the minimum number of nodes after which an increase in the number of nodes wouldn't affect the concentration measured at the outlet significantly. The study was carried for MEA weight percentage of 30 %, at Re L of 125 and Re G of 200. The resulting graph, Appendix 1, of liquid outlet concentration as a function of number of nodes clearly shows that a number of 160 nodes (6.25 µm) in the liquid inlet region of 1 mm are sufficient to capture the outlet concentration.