*Article* **Droplet Characteristics of Rotating Packed Bed in H2S Absorption: A Computational Fluid Dynamics Analysis**

#### **Zhihong Wang \*, Xuxiang Wu, Tao Yang \*, Shicheng Wang, Zhixi Liu and Xiaodong Dan**

School of Chemistry and Chemical Engineering, Southwest Petroleum University, Chengdu 610500, China; Wxx2alley@163.com (X.W.); 18200351814@163.com (S.W.); 18780273398@163.com (Z.L.); dongceswpu@163.com (X.D.)

**\*** Correspondence: wzhswpu@swpu.edu.cn (Z.W.); yangtao9264826@outlook.com (T.Y.); Tel.: +86-1390-8215-506 (Z.W.); +86-1354-1167-685 (T.Y.)

Received: 7 August 2019; Accepted: 8 October 2019; Published: 11 October 2019

**Abstract:** Rotating packed bed (RPB) has been demonstrated as a significant and emerging technology to be applied in natural gas desulfurization. However, droplet characteristics and principle in H2S selective absorption with *N*-methyldiethanolamine (MDEA) solution have seldom been fully investigated by experimental method. Therefore, a 3D Eulerian–Lagrangian approach has been established to investigate the droplet characteristics. The discrete phase model (DPM) is implemented to track the behavior of droplets, meanwhile the collision model and breakup model are employed to describe the coalescence and breakup of droplets. The simulation results indicate that rotating speed and radial position have a dominant impact on droplet velocity, average residence time and average diameter rather than initial droplet velocity. A short residence time of 0.039–0.085 s is credited in this study for faster mass transfer and reaction rate in RPB. The average droplet diameter decreases when the initial droplet velocity and rotating speed enhances. Restriction of minimum droplet diameter for it to be broken and an appropriate rotating speed have also been elaborated. Additional correlations on droplet velocity and diameter have been obtained mainly considering the rotating speed and radial position in RPB. This proposed formula leads to a much better understanding of droplet characteristics in RPB.

**Keywords:** rotating packed bed; natural gas desulfurization; droplet characteristic; Eulerian–Lagrangian approach; computational fluid dynamics

#### **1. Introduction**

Rotating packed bed (RPB) has been demonstrated as a novel and efficient equipment in process intensification of mass transfer and reaction [1], as displayed in Figure 1. Gas phase and liquid phase contact counter-currently in the packing by centrifugal force, and the height of mass transfer unit (HTU) in RPB can dramatically decrease by 1–3 orders of magnitude compared with that in a conventional column. Therefore, RPB can evidently enhance productivity while also remarkably reducing the device size and equipment investments. RPB has been maturely applied to the traditional chemical industry [2] in processes such as distillation, absorption and extraction. It has gradually begun to play a crucial role in gas purification, especially for H2S removal in natural gas.

**Figure 1.** Schematic diagram of rotating packed bed (RPB) (counter-currently).

Various investigations have been carried out to reveal the gas–liquid behavior in RPB by experimental method. Burns and Ramshaw [3] exhibited a visual liquid patterns, shown in Figure 2, in the way that a video camera was arranged rotating synchronously with the packing. Liquid patterns can be distinguished into pore flow, droplet flow and film flow in different stages and positions. At a lower rotating speed, the rivulet flow was dominant while droplet flow mainly existed at a higher rotating speed. Guo [4] found that liquid impinged and deformed intensively within the initial 7–10 mm from the inner diameter and that liquid performed two states in the packing: Films on the packing surface and films flying in the gap of packing. Moreover, Li [5] investigated the effect of a various operating conditions especially for the layer numbers on the liquid patterns and drew a conclusion that the droplet was the main pattern in packing. As for recent research, Sang [6] indicated a criterion to distinguish two typical liquid patterns: Ligament flow and droplet flow. The Computational fluid dynamics is called CFD for short, which is a combination of computer and numerical techniques. The nature of CFD is that some physical phenomena can be obtained by solving related partial differential equations in computer. Due to the lower cost and visualization of process and results of CFD technology, the application of CFD has become a powerful approach to depict the liquid patterns and behaviors in RPB. Shi [7] verified that the flow patterns were varied under different rotating speeds in CFD simulation. Ouyang emphasized that a higher viscosity of liquid led a liquid line predominantly [8]. In another simulation article, Ouyang et al. pointed out that droplet diameter increased with decreasing rotating speed, liquid initial velocity and number of layers in a Rotor–Stator reactor (RSR) [9]. Conclusively, Xie [10] figured out that the collision and merging of droplets occurred between droplets and packing and that the liquid flow tended to form droplets under the influence of surface tension (larger contact angles) and high rotating speeds (1000–1500 rpm).

**Figure 2.** Three types of liquid pattern in RPB with a video camera, reproduced with permission from [3]. Copyright Elsevier, 1996.

Aforementioned visual studies through experimental method and CFD method elaborate that typical flow pattern in RPB is liquid droplet. Although characteristics of droplets such as velocity, average residence time and average diameter exert great influence in design and improvement of RPB, there are still few studies on them. Consequently the present work investigated the influence of rotating speed, initial droplet velocity and initial droplet diameter on droplet velocity, average droplet residence time as well as average droplet diameter by CFD method. Through introducing the collision and breakup model, correlations of droplet velocity and average diameter were proposed. Additionally, a phenomenon called "end effect" in the end zone and a principle of processing intensification contributing to H2S selective absorption were also presented.

#### **2. Simulation**

#### *2.1. Physical Model and Grid Refinement of RPB*

The structure of RPB in this study is based on a pilot installation built in a natural gas purification plant, and Table 1 displays the main dimensions of RPB. It has a 24 mm height, 45 mm inner diameter and 160 mm outer diameter. Packing is regarded as a crucial part of RPB and it generally consists of uniform thickness wire, which has a 20 mm height, 48 mm inner diameter and 92 mm outer diameter.


**Table 1.** The main geometric dimensions of RPB.

In order to make the 3D model as consistent as possible with the actual structure, there are three parts in 3D model: End zone, packing zone and cavity zone (Figure 3a). The packing zone is built with 21 concentric layers, and the foursquare obstacles model is employed to simulate the wire behavior in RPB. Consequently, the side length of each foursquare is 0.8 mm and the center distances between two foursquares in radial and tangential directions are 2 mm and 5 mm, respectively. Four annulus distributors set uniformly in the inner radius of RPB are set as liquid inlets and an outer ring is the outlet.

A computational 3D model grid is built by ICEM 19.0. In order to describe more details, a grid refinement has been made, namely different numbers of cells are meshed in different parts. At the same time, the packing zone and the end zone are meshed by tetrahedral cell, and the cavity zone is meshed by mixed cell. The meshing result is shown in Figure 3b.

**Figure 3.** (**a**) 3D physical model diagram (1, gas inlet; 2, liquid inlet; 3, gas outlet; 4, foursquare obstacles; 5, liquid outlet; 6, end zone; 7, packing zone; 8, cavity zone). (**b**) Grid refinement.

#### *2.2. Mathematical Modelling*

In this study, based on an assumption that the second phase as a dispersed part occupies a low volume fraction (below 10%), the Eulerian–Lagrangian approach is used in which the natural gas is considered as a continuous phase while the MDEA solution is regarded as a discrete dispersed phase. Forces, like drag force, buoyancy force and virtual mass force, are taken into account to describe the particle force balance. A turbulence model is introduced to reflect the turbulent process when packing shears the MDEA solution intensely under high rotation speed. Models, such as Taylor analogy breakup (TAB) model and O'Rourke model, are carried out to predict the breakup and coalescence of droplets.

#### 2.2.1. Governing Equations

The basic equations coupled in CFD models for solving mass and momentum conservation equations are listed below for an incompressible, laminar flow, non-accelerating and unsteady state process:

**Mass conservation equation**

$$\frac{\partial \rho}{\partial t} + \nabla \cdot \left(\rho \overrightarrow{\boldsymbol{\nu}} \right) = 0,\tag{1}$$

#### **Momentum conservation equation**

$$\frac{\partial(\rho \overrightarrow{\boldsymbol{\upsilon}})}{\partial t} + \nabla \cdot \left(\rho \overrightarrow{\boldsymbol{\upsilon}} \overrightarrow{\boldsymbol{\upsilon}}\right) = -\nabla p + \rho \overrightarrow{\boldsymbol{\varsigma}} + \nabla \cdot \left(\mu\_{eff} \left(\overrightarrow{\boldsymbol{\upsilon}} \overrightarrow{\boldsymbol{\upsilon}} + \left(\overrightarrow{\boldsymbol{\upsilon}} \overrightarrow{\boldsymbol{\upsilon}}\right)^{T}\right)\right) + \overrightarrow{F},\tag{2}$$

where *p* is the static pressure, ρ is the fluid density, <sup>→</sup> *g* is the gravitational vector, μ*eff* is the molecular viscosity and <sup>→</sup> *F* is the external body force representing the interactions between the dispersed phase and the continuous phase. As for a flow in a rotating reference frame, like the rotation in RPB, the origin of the non-accelerating system is determined by a centrifugal acceleration <sup>→</sup> ω in a rotating system, hence the governing equations of flow in a rotating reference frame are as follows:

$$\frac{\partial \rho}{\partial t} + \nabla \cdot \left(\rho \overrightarrow{\boldsymbol{\upsilon}}\_{r}\right) = 0,\tag{3}$$

$$\frac{\partial \left(\rho \overrightarrow{\boldsymbol{\upsilon}}\_{r}\right)}{\partial t} + \nabla \cdot \left(\rho \overrightarrow{\boldsymbol{\upsilon}}\_{r} \overrightarrow{\boldsymbol{\upsilon}}\_{r}\right) = -\nabla p + \rho \left(2\overrightarrow{\boldsymbol{\omega}} \times \overrightarrow{\boldsymbol{\upsilon}}\_{r} + \overrightarrow{\boldsymbol{\omega}} \times \overrightarrow{\boldsymbol{u}}\_{r}\right) + \nabla \cdot \left(\mu\_{eff} \left(\nabla \overrightarrow{\boldsymbol{\upsilon}}\_{r} + \left(\nabla \overrightarrow{\boldsymbol{\upsilon}}\_{r}\right)^{T}\right)\right) + \overrightarrow{F}\_{r} \tag{4}$$

$$
\overrightarrow{\boldsymbol{v}}\_r = \overrightarrow{\boldsymbol{v}} - \overrightarrow{\boldsymbol{u}}\_r \tag{5}
$$

$$
\overrightarrow{\boldsymbol{\mu}}\_{r} = \overrightarrow{\boldsymbol{\omega}} \times \overrightarrow{\boldsymbol{r}}\_{r} \tag{6}
$$

where <sup>→</sup> *v <sup>r</sup>* is the relative velocity, <sup>→</sup> ω is the centrifugal acceleration, <sup>→</sup> <sup>ω</sup> <sup>×</sup> <sup>→</sup> *v <sup>r</sup>* is the Coriolis acceleration, <sup>→</sup> *ur* is the whirl velocity and <sup>→</sup> <sup>ω</sup> <sup>×</sup> <sup>→</sup> *ur* is the centripetal acceleration.

Although the Navier–Stokes equation can accurately describe the mass and momentum of the laminar flow, the turbulence caused by intensive movement can evidently influence the transport process. Therefore, the turbulence model determines additional variables as time-averaged or ensemble-averaged in the modified governing equations.

#### 2.2.2. Turbulence Model

The standard k–ε model is presented by molecular viscosity and dynamic viscosity due to the turbulence and it is a two-equation model. It is a semi-empirical model, and the solution of the k equation allows the turbulent kinetic energy to be determined while the result of the ε equation can be considered relying on phenomenological considerations and empiricism.

$$\frac{\partial}{\partial t}(\rho k) + \nabla \cdot \left(\rho k \overrightarrow{\boldsymbol{\nu}}\_{r}\right) = \nabla \cdot \left(\left(\mu + \frac{\mu\_{t}}{\sigma\_{k}}\right) \nabla \mathbf{k}\right) + \mathbf{G}\_{k} + \mathbf{G}\_{b} - \rho \varepsilon\_{r} \tag{7}$$

$$\frac{\partial}{\partial t}(\rho \varepsilon) + \nabla \cdot (\rho \varepsilon \stackrel{\scriptstyle \gamma}{\upsilon}\_{\mathcal{I}}) = \nabla \cdot \left( \left( \mu + \frac{\mu\_{\mathcal{I}}}{\sigma\_{\mathcal{I}}} \right) \nabla \varepsilon \right) + \mathfrak{c}\_{\mathcal{E}1} \frac{\varepsilon}{k} \left( \mathcal{G}\_{k} + \mathfrak{c}\_{\mathcal{E}3} \mathcal{G}\_{\mathcal{b}} \right) - \mathfrak{c}\_{\mathcal{E}2} \rho \frac{\varepsilon^{2}}{k} \tag{8}$$

$$G\_k = \mu\_t \nabla \overline{\boldsymbol{\upsilon}} \cdot \left(\nabla \overline{\boldsymbol{\upsilon}}^\circ + \left(\nabla \overline{\boldsymbol{\upsilon}}^\circ\right)^T\right) \tag{9}$$

$$\mathbf{G}\_b = -\frac{\mu\_t}{\rho Pr\_t} \mathbf{g} \nabla \rho\_\prime \tag{10}$$

$$
\mu\_t = \mathbb{C}\_{\mu} \rho \frac{k^2}{\varepsilon},
\tag{11}
$$

where *Gk* represents the influence of the mean velocity gradients, *Gb* represents the influence of the buoyancy force, μ*<sup>t</sup>* is the turbulent viscosity, *Prt* is the turbulent Prandtl number for energy, *k* is the turbulence kinetic energy and ε is the dissipation rate. *c*ε1, *c*ε2, *c*ε3, *C*μ, σ*k*, σε are the model constants and the default values are: *c*ε<sup>1</sup> = 1.44, *c*ε<sup>2</sup> = 1.92, *C*<sup>μ</sup> = 0.09, σ*<sup>k</sup>* = 1.0, σε = 1.3.

#### 2.2.3. Droplet Force Balance

In a Lagrangian reference frame, the force balance on the particle which drives the particle acceleration is due to the difference in velocity between fluid and particles as well as various influences of force, like drag force, gravity, virtual mass force and centrifugal force in RPB. Hence, a particle force balance can be written as:

$$\frac{du\_p}{dt} = F\_D(u - u\_p) + \frac{\left(\rho\_p - \rho\right)}{\rho\_p}g + F\_{x\_{\prime}} \tag{12}$$

**Drag force**

$$F\_D = \frac{18\mu}{\rho\_p d\_p^2} \frac{\mathbb{C}\_D \text{Re}}{24},\tag{13}$$

$$\text{Re} = \frac{\rho d\_p |u\_p - u|}{\mu},\tag{14}$$

where *FD u* − *up* is the drag force per unit particle mass, (ρ*p*−ρ) <sup>ρ</sup>*<sup>p</sup> g* is the gravity per unit particle mass, *Fx* is the additional acceleration term (force per unit particle mass), *CD* is the drag coefficient, Re is the relative Reynolds number, *up* is the particle velocity, *u* is the fluid velocity and ρ*<sup>p</sup>* is the particle density.

**Virtual mass force**

$$F\_m = \frac{1}{2} \frac{\rho}{\rho\_p} \frac{d}{dt} (u - u\_p)\_{\prime} \tag{15}$$

#### **Forces in rotating reference frames**

$$F\_{r, \mathbf{x}} = \left(1 - \frac{\rho}{\rho\_P}\right) \omega^2 \mathbf{x} + 2\omega \left(u\_{y, p} - \frac{\rho}{\rho\_P} u\_y\right) \tag{16}$$

$$F\_{r,y} = \left(1 - \frac{\rho}{\rho\_p}\right) \omega^2 y + 2\omega \left(u\_{x,p} - \frac{\rho}{\rho\_p} u\_x\right) \tag{17}$$

where *Fm* is the virtual mass force, *Fr*,*<sup>x</sup>* and *Fr*,*<sup>y</sup>* are the forces in rotating reference frames in *x* and y directions, respectively and ω is the rotating speed.

#### 2.2.4. Droplet Coalescence and Breakup Model

#### **The coalescence model**

The O'Rourke model a is model based on energy balance; the droplets will coalescence once the kinetic energy cannot overcome the surface energy of a new droplet [11]. We is defined as a ratio of inertial force to surface tension to quantify the collision and coalescence of droplets. Therefore, the critical offset *b* is a function of the collisional, shown in Figure 4.

$$\text{We} = \frac{\rho\_p \left| \overrightarrow{u}\_I - \overrightarrow{u}\_s \right|^2 \overline{d}}{\sigma}, \tag{18}$$

$$\text{We} \ge \frac{2.4(r\_l + r\_s)^2 f\left(\frac{r\_l}{r\_s}\right)}{b^2},\tag{19}$$

$$b\_{crit} = (r\_l + r\_s) \sqrt{\min\left\{ 1.0, \frac{2.4f\left(\frac{r\_l}{r\_s}\right)}{\text{We}} \right\}},\tag{20}$$

$$f\left(\frac{r\_l}{r\_s}\right) = \left(\frac{r\_l}{r\_s}\right)^3 - 2.4\left(\frac{r\_l}{r\_s}\right)^2 + 2.7\frac{r\_l}{r\_s},\tag{21}$$

where We is the Weber number, <sup>→</sup> *ul* and <sup>→</sup> *us* are the velocities of large droplets and small droplets, respectively, *d* is the arithmetic mean diameter of two droplets, *rl* and *rs* are the radii of large droplets and small droplets, respectively, *b* and *bcrit* are the actual collision parameter and the critical offset of collision, respectively. If *b* < *bcrit*, the result of collision is coalescence, and the new velocity based on conservation of momentum and kinetic energy is calculated below:

$$v' = \frac{m\_l v\_l + m\_s v\_s + m\_s (v\_l - v\_s)}{m\_l + m\_s} \left(\frac{b - b\_{crit}}{r\_l + r\_s - b\_{crit}}\right) \tag{22}$$

where *ml* and *ms* are the masses of large droplets and small droplets, respectively.

**Figure 4.** Schematic diagram of the coalescence process.

#### **Taylor Analogy Breakup (TAB) model**

The TAB model governs the oscillating and distorting droplet; when the droplet oscillations grow to a critical value, the large droplet will break up into small droplets. The droplet distortion governing equation is defined as follows:

$$F - kx - d\frac{d\mathbf{x}}{dt} = m\frac{d^2\mathbf{x}}{dt^2},\tag{23}$$

where *x* is the displacement of the droplet, *m* is the droplet mass, *F* is the aerodynamic force of the droplet, *kx* is the surface tension of the droplet and *ddx dt* is the viscous force.

$$\frac{F}{m} = C\_F \frac{\rho \left| \mu - \mu\_p \right|^2}{\rho\_P r\_p} \,, \tag{24}$$

$$\frac{d}{m} = \mathbb{C}\_d \frac{\rho \mu\_p}{\rho\_p r\_p^2},\tag{25}$$

$$\frac{k}{m} = \mathbb{C}\_k \frac{\sigma\_p}{\rho\_p r\_p^3} \tag{26}$$

$$\frac{d^2y}{dt^2} = \frac{\mathbb{C}\_F}{\mathbb{C}\_b} \frac{\rho \left| \overrightarrow{u} - \overrightarrow{u}\_p \right|^2}{\rho\_p r\_p^2} - \mathbb{C}\_k \frac{\sigma\_p}{\rho\_p r\_p^3} y - \mathbb{C}\_d \frac{\mu\_p}{\rho\_p r\_p^2} \frac{dy}{dt},\tag{27}$$

where *CF*, *Cd*, *Ck* are dimensionless constants, *CF* = 1/3, *Cd* = 5, *Ck* = 8. Setting y = x/ *Cbrp* , when *Cb* = 0.5 and y > 1, the non-dimensional equation for droplet breakup is shown in Equation (28).

$$\mathbf{y}(t) = \mathbf{W}\mathbf{e}\_{r} + \mathbf{e}^{-\frac{t}{t\_d}} \left\{ (y\_0 - \mathbf{W}\mathbf{e}\_r)\cos(\omega t) + \frac{1}{\omega} (\frac{dy\_0}{dt} + \frac{y\_0 - \mathbf{W}\mathbf{e}\_r}{t\_d}) \sin(\omega t) \right\},\tag{28}$$

$$\text{We}\_{r} = \frac{\prescript{C}{C\_F}}{\prescript{C\_b}{C\_k}} \frac{\rho \left| \overrightarrow{\boldsymbol{\mu}} - \overrightarrow{\boldsymbol{\mu}}\_p \right|^2 r\_p}{\sigma\_p},\tag{29}$$

$$t\_d = \frac{2\rho\_p r\_p^2}{\mathcal{C}\_d \mu\_p} \,\,\,\,\,\tag{30}$$

$$
\omega \omega^2 = \mathbb{C}\_k \frac{\sigma\_p}{\rho\_p r\_p^3} - \frac{1}{t\_d^2}. \tag{31}
$$

The size of droplet is determined by the energy conservation between the parent droplet and the child droplets which satisfies the Sauter distribution and the number of child droplets can be calculated by mass conservation.

$$r' = \frac{1 + \frac{1}{20}y^2 + \frac{1}{8}\frac{\rho\_p r\_p^3}{\sigma\_p} \left(\frac{dy}{dt}\right)^2}{r\_p} \tag{32}$$

#### *2.3. Fluid Properties*

MDEA aqueous solution is one of the conventional absorbents for H2S selective absorption. In this study, the gas mixture consists mainly of methane and the mass fraction of the aqueous solution is 35%. The properties of the gas and liquid used for the CFD simulations are shown in Table 2. The MDEA aqueous solution is assumed to operate at a constant temperature of 30 ◦C and 2 MPa, which are close to the real operation conditions for H2S selective absorption in RPB. Under these temperature and pressure conditions, the density and viscosity of the gas mixture are 15.99 kg/m3 and 1.203 <sup>×</sup> 10−<sup>5</sup> kg/(m·s), respectively; the density and viscosity of the aqueous solution are and 1027 kg/m3 and 1.203 <sup>×</sup> <sup>10</sup>−<sup>5</sup> kg/(m·s), respectively.

**Table 2.** Fluid properties.


#### *2.4. Solution Procedure*

The domain of the packing zone is defined as the moving reference frame (MFR) and the rest of RPB: End zone and cavity zone are relative to a stationary reference frame in the Fluent 19.0 software. Thus, the packing with a reflect boundary rotates in a rotating speed of 300 to 900 rpm and the wall condition with no-slip and escape boundary is set to simulate the casing of RPB. The gas inlet and outlet in RPB are set as the velocity inlet boundary and pressure outlet boundary conditions, respectively, while the liquid inlet and outlet are defined as discrete phase surface injection and trap boundary condition, respectively. The gas velocity is considered as a constant 0.5 m/s, and the droplets are injected into the RPB in this way that initial droplet velocity is specified by a range of 0.5 to 2.5 m/s, the initial droplet diameter is set to 1 to 5 mm and the number of droplets injected into the RPB is 400. Steady simulations are implemented with standard k–ε model and DPM model to investigate the droplet characteristics in RPB. The droplet discrete phase interacts with the continuous phase every five iterations and an unsteady particle tracking is applied to track the droplet trajectories individually with a time step size of 0.001 s for a maximum of 5000 steps. Drag force, gravity, virtual mass force and centrifugal force are considered in this solution and droplet coalescence and breakup models mentioned above are also introduced. The pressure–velocity coupling is resolved by SIMPLE scheme and spatial discretization methods for gradient and pressure are least squares cell-based and PRESTO, respectively. The second-order upwind scheme is employed for solving the momentum equations and turbulence equations. For the convergence absolute criteria, the residuals in mass conservation,

momentum conservation and turbulence are less than 1 <sup>×</sup> <sup>10</sup><sup>−</sup>5. In addition, the maximum number of iterations for steady calculation is 4000 to achieve the steady state in RPB.

#### *2.5. Grid Independence*

As is known, the number of grids has a great impact on the accuracy of the simulation result. Therefore, a grid independence study is performed to obtain a reasonable computational mesh. Fixing the droplet diameter of 2 mm, rotating speed of 300 rpm and initial droplet velocity of 0.5 m/s as the initial conditions and choosing the droplet average residence time as the only dependent variable, five different grids consisting of 1.4, 2.5, 3.3, 4.3 and 5 million unstructured cells have been implemented to analyze the effect of the cell on the droplet average residence time. As shown in Figure 5, the average residence time is almost same when the number of cells reaches 4.3 million, which may be considered as a responsible grid to predict droplet characteristics. Finally, the number of 4,349,428 cells is selected to investigate other simulation results while considering the computing resources and accuracy.

**Figure 5.** The grid independence checking.

#### **3. Results and Discussions**

Since droplet form in RPB is the dominant liquid pattern, its characteristics such as velocity, residence time and diameter in mass transfer and reactions will play a crucial role in novel research. Therefore, a CFD simulation on droplet characteristics will demonstrate a visual flow distribution as well as significant guide for H2S selective absorption into MDEA solution in RPB.

#### *3.1. Droplet Velocity in RPB*

Velocity distribution, as a common and crucial characteristic of droplets, has a direct effect on residence time and diameter, therefore the velocity value and vector under various rotating speeds and initial droplet velocities are shown in Figure 6. Similarly, the velocity is distributed uniformly and symmetrically along the radial and tangential directions. The liquid is injected and dispersed uniformly into packing zone from the liquid distributor and it rapidly obtains a synchronous tangential speed as it is packing, finally liquid droplets are spun out from the outer packing zone into the cavity randomly.

(**a**)

(**b**)

**Figure 6.** *Cont.*

(**c**)

(**d**)

**Figure 6.** *Cont.*

**Figure 6.** The velocity values and vectors under various initial droplet velocities and rotating speeds. (**a**) *u*<sup>0</sup> = 0.5 m/s, ω = 300 rpm; (**b**) *u*<sup>0</sup> = 0.5 m/s, ω = 600 rpm; (**c**) *u*<sup>0</sup> = 0.5 m/s, ω = 900 rpm; (**d**) *u*<sup>0</sup> = 1.5 m/s, ω = 600 rpm; (**e**) *u*<sup>0</sup> = 2.5 m/s, ω = 600 rpm.

In the end zone, a parabolic droplet motion is observed with the increase of rotating speed from 300–900 rpm, since the initial droplet velocity in the radial direction is lower than that in tangential direction. Further, a phenomena of back-mixing, called "end effect" in our review [12] on mass transfer process in RPB, was also confirmed to have a superiority in processing intensification on mass transfer [13]. In the packing zone, droplets move in a radial spiral and distribute irregularly. An analysis on different rotating speeds and initial droplet velocities illustrates that the radial velocity is mainly affected by initial droplet velocity while the tangential velocity is related to rotating speed and radial position. In the cavity zone, droplets move towards the outer wall of RPB with a constant angle which decreases with the increase of rotating speed, and it manifests that rotating speed has a major impact on resultant droplet velocity.

Under the sustaining shear from the packing, a large proportion of droplets are synchronous with packing and the droplet velocity increases linearly along the radial position. Some fluctuation of droplet velocity has also been observed because the collision, breakup and coalescence of droplets take place among droplet–droplet, droplet–gas and droplet–packing interactions.

#### 3.1.1. Effect of Initial Droplet Velocity on Droplet Velocity

Figure 7 show the effect of initial droplet velocity on droplet velocity in the packing zone under the condition that initial droplet diameter is 3 mm and rotating speed is 900 rpm. In this zone, the droplet velocity is seldom influenced by initial droplet velocity. The reason is that the initial droplet velocity makes less contribution to resultant velocity compared with the tangential velocity generated by centrifugal force. On the other hand, the inelastic collision on droplets results in the loss of kinetic energy once droplets are sheared by packing.

**Figure 7.** The effect of initial droplet velocity on droplet velocity in the packing zone.

#### 3.1.2. Effect of Rotating Speed on Droplet Velocity

Figure 8 establishes the effect of rotating speed on droplet velocity in the packing zone under the condition that initial droplet diameter is 3 mm and initial droplet velocity is 0.5 m/s. In the packing zone, the rotating speed is dominant for droplet velocity increasing and it illustrates that the droplet velocity increases along with rotating speed increasing. As a consequence, Equation (33) for predicting droplet velocity in RPB is proposed by considering the influence of initial droplet velocity, rotating speed and radial position; validation Equation (34), predicted from Sang [14], is implemented. Finally, the mathematic expression declares that rotating speed and radial position have the main impacts on droplet velocity.

$$
u\_p = 0.6088\omega^{0.9819} u\_0^{0.0008} R^{0.7702} \,\text{s}}\tag{33}$$

$$
\mu\_p = 1.013 \omega^{0.9994} \mu\_0^{0.008} R^{1.005} \mu^{-0.009} \sigma^{0.036}.\tag{34}
$$

**Figure 8.** The effect of rotating speed on droplet velocity in the packing zone.

#### *3.2. Average Residence Time Distribution in RPB*

Analysis on selective absorption with MDEA indicated that it was kinetically selective towards H2S while thermodynamically selective towards CO2 [15]. According to the theory of the diffusion–reaction process for H2S selective absorption, the reaction and mass transfer on H2S are both instantaneously fast, while those process are restrained in CO2 mass transfer into the liquid film. Therefore, average residence time of droplet is considered as a key parameter in the H2S selective absorption into MDEA process in RPB. A method to obtain the residence time distribution is applied by tracking droplets, and residence time distributions in RPB under various initial droplet velocities and rotating speeds are shown in Figure 9.

Droplets interact intensely with the rotating packing and they are driven by centrifugal force, gravity, gas viscous resistance and packing frictional resistance, resulting in a frequent path change. Most of the droplets are accelerated to a synchronous motion with packing and are ejected in a radial spiral. With the increase of initial droplet velocity and rotating speed, the turbulence intensity increases with the increase of random trajectory. Table 3 summarizes the number of tracked droplets and average residence time in RPB under different rotating speeds and initial droplet velocities.


**Table 3.** Average residence time under various initial droplet velocities and rotating speeds.

(**a**)

(**b**)

**Figure 9.** *Cont.*

(**c**)

(**d**)

**Figure 9.** *Cont.*

(**e**)

**Figure 9.** Residence time distribution in RPB under various initial droplet velocities and rotating speeds. (**a**) *u*<sup>0</sup> = 0.5 m/s, ω = 300 rpm; (**b**) *u*<sup>0</sup> = 0.5 m/s, ω = 600 rpm; (**c**) *u*<sup>0</sup> = 0.5 m/s, ω = 900 rpm; (**d**) *u*<sup>0</sup> = 1.5 m/s, ω = 600 rpm; (**e**) *u*<sup>0</sup> = 2.5 m/s, ω = 600 rpm.

3.2.1. Effect of Initial Droplet Velocity on Average Residence Time

At the same rotating speed of 900 rpm, the average residence time changes drastically from a low initial velocity to a high initial velocity. Because of the increase of initial velocity, the time to pass through the whole packing will shorten correspondingly. As it is displayed in Figure 10, there is a turning point when the initial velocity reaches 1.5 m/s and the average residence time is reduced gradually.

**Figure 10.** Effect of initial droplet velocity on average residence time.

#### 3.2.2. Effect of Rotating Speed on Average Residence Time

As shown in Figure 11, the droplet attains a tangential velocity when it enters the packing zone and a synchronous velocity will be achieved in a short time. For the sake of reducing residence time in RPB for H2S selective absorption, enhancing rotating speed seems to be an optimal way. Lower residence time in RPB compared with that in conventional column is ascribed to intensification process under centrifugal force. The higher the droplet velocity is, the lower the residence time obtained will be. Therefore, a commonsense calculation can be proposed to describe the residence time in RPB, and that is the ratio of radial distance to droplet velocity, calculated by Equation (35).

$$t = \frac{R}{u\_p} \tag{35}$$

In the process of natural gas desulfurization and purification, H2S selective absorption and CO2 partial removal must meet the requirements on commercial natural gas, and the energy consumption also needs to be reduced in MDEA solution regeneration process. Therefore, higher mass transfer efficiency and lower residence time in RPB raise various concerns in natural gas purification. In the quantitative description of Qian's work [16], it only took about 2.0 <sup>×</sup> <sup>10</sup>−<sup>9</sup> s for H2S to establish a steady concentration gradient while the process needed 1 s to complete for CO2. Thus, it can be highlighted from Figure 10, that the residence time can be only 0.039–0.085 s under a rotating speed of 900 rpm. Droplets go through the packing in a short residence time, especially with high initial droplet velocity and rotating speed.

**Figure 11.** Effect of rotating speed on average residence time.

#### *3.3. Droplet Diameter Distribution in RPB*

Droplet surface renewal frequency, which is defined as the percentage of surface updated per unit of time and effective mass transfer area, can be investigated by the droplet diameter distribution along radial position. Figure 12 shows the 3D droplet diameter distribution with an initial droplet velocity of 1.5 m/s, rotating speed of 600 rpm and initial droplet diameter of 3 mm. The droplet diameter changes a lot along the radial position, and collision and shear between droplet and packing result in sustaining breakup and coalescence. This also indicates that frequent surface renewal makes a prominent contribution to mass transfer.

**Figure 12.** A 3D droplet diameter distribution in RPB.

3.3.1. Effect of Initial Droplet Diameter on Droplet Diameter Distribution

Figure 13 reveals that the droplet diameter increases with an increasing initial droplet diameter. After a larger droplet is introduced from the inlet, it tends to be broken up by packing because of lower surface tension needed to be overcome, and this results in a frequent surface renewal and an intensifying mass transfer process. Guo [17] obtained a formula by fitting relevant experimental data, *d* = 0.7284 σ/ρω2*R* 0.5, which indicated that the droplet diameter decreases along the radial direction. Both Figures 13 and 14 show that the average droplet diameter is negatively correlated with the radius position, which is consistent with Guo's research. Therefore, it also manifests that the CFD model introduced in this study reasonably predicts the droplet diameter distribution.

(**a**)

**Figure 13.** *Cont.*

**Figure 13.** The droplet diameter distribution under different initial droplet diameters. (**a**) *d*<sup>0</sup> = 1 mm; (**b**) *d*<sup>0</sup> = 3 mm; (**c**) *d*<sup>0</sup> = 5 mm.

**Figure 14.** Effect of initial droplet diameter on average droplet diameter.

Specially, Figure 14 illuminates a particular phenomenon: There is rarely collision in the packing zone when the initial droplet diameter is 1 mm, compared with others. The reasons for the above phenomenon are as follows: (1) the packing size is larger than the initial droplet diameter; and (2) the inertial force under this condition is not enough to break the small droplet. Therefore, the study of the collision between 1-mm-diameter initial droplets and packing needs smaller packing size and larger inertial force. In a word, the effect of initial droplet diameter on droplet diameter distribution analyzed above can not only provide guidance for liquid distributor design, but also make a suitable proposal on choosing packing size.

#### 3.3.2. Effect of Initial Droplet Velocity and Rotating Speed on Droplet Diameter

Figures 15 and 16 illustrate that the average droplet diameter along the radial position decreases with the increase of initial droplet velocity from 0.5 to 2.5 m/s and with the increase of rotating speed from 300 rpm to 900 rpm. The droplet obtains the tangential velocity from the packing immediately, resulting in a violent collision and breakup among droplet and packing. Therefore, the droplet surface renewal frequency and effective mass transfer area have been enhanced remarkably.

(**a**)

(**b**)

**Figure 15.** *Cont.*

**Figure 15.** Droplet diameter under various initial droplet velocities and rotating speeds. (**a**) *u*<sup>0</sup> = 0.5 m/s, ω = 600 rpm; (**b**) *u*<sup>0</sup> = 2.5 m/s, ω = 600 rpm; (**c**) *u*<sup>0</sup> = 1.5 m/s, ω = 300 rpm; (**d**) *u*<sup>0</sup> = 1.5 m/s, ω = 900 rpm.

(**b**)

**Figure 16.** Effect of initial droplet velocity and rotating speed on average droplet diameter. (**a**) *d*<sup>0</sup> = 3 mm, ω = 600 rpm; (**b**) *d*<sup>0</sup> = 3 mm, *u*<sup>0</sup> = 1.5 m/s.

As shown in Figure 16a, under the same conditions of initial droplet diameter of 3 mm and rotating speed of 600 rpm, the average droplet diameter decreases drastically with the increase of initial droplet velocity from 0.5 to 2.5 m/s in the packing zone. In other words, the average droplet diameter with initial velocity of 2.5 m/s is less than the average droplet diameter with initial velocity of 1.5 m/s, which is less than the average droplet diameter with initial velocity of 0.5 m/s in same radial position. A higher initial droplet velocity means a greater initial kinetic energy, which is beneficial to the breakup of droplets when kinetic energy is greater than surface energy, so the phenomenon is obtained in Figure 16a. In addition, the initial droplet diameter of 3 mm will gradually break up into smaller droplets along the radial position and finally reaches a constant diameter of 1.25 mm at the outer packing zone. Thus, conclusions can be drawn that the minimum droplet diameter in RPB only depends on the rotating speed and that there exists a minimum rotating packing radial length for droplets to fully fragmentize.

From Figure 16b, it can be observed distinctly that average droplet diameter distribution changes regularly on initial droplet diameter and velocity of 3 mm and 1.5 m/s, respectively. The effect of rotating speed on droplet diameter indicates that the droplet diameter decreases with the increasing rotating speed. The centrifugal force is linear to radial position and is proportional to the second power of rotating speed. Therefore, a higher centrifugal force results in more drastic collisions and smaller droplet diameters. Further, a balance between shear force generated by rotating packing and surface tension on droplets is proposed and discussed under an appropriate rotating speed. When the rotating speed decreases to 300 rpm, the shear force cannot meet the surface tension, resulting in only a few droplet breakups at the end of packing. On the contrary (at 900 rpm), shear force is dominant instead of surface tension, so the droplet diameter changes gently.

Compared with the effect of initial droplet velocity on droplet diameter, a more obvious impact of rotating speed is obtained. Radial velocity and tangential velocity all contribute to droplet breakup, but the tangential velocity caused by rotating speed is dominant. Therefore, considering the effect of rotating speed, radial position and fluid density on droplet diameter, a fitting expression is obtained:

$$d\_p = 14.96 \left[ \frac{\rho}{(\rho\_p - \rho)\omega^2 R} \right]^{0.8805} \tag{36}$$

#### *3.4. Principle of Processing Intensification in RPB*

A comparison on particle diameter is conducted between RPB and conventional deposition process under gravity. Droplet force balance has been discussed in Equation (12) and the droplet diameter is calculated by two equations in Table 4, which were derived from Equations (12) and (13), considering mainly the centrifugal force in RPB (or gravity in deposition process) and drag force. Droplet diameter in different force fields is related to fluid density, droplet density, fluid velocity, droplet velocity and main acceleration. Generally speaking, the behavior in RPB can be considered as a processing intensification in force fields by replacing the term "*g*" with "ω2*R*". Thus, greater acceleration results in smaller droplet diameter and larger droplet surface area of mass transfer. According to the theory of mass transfer, the mass transfer rate is proportional to the surface area. Therefore, the mass transfer efficiency of RPB is 1–3 orders of magnitude greater than that in conventional equipment. In addition, some intensification in principle mentioned above will be implemented by coupling multiple fields like magnetic, electric and microwaves to explore their potential in RPB.

**Table 4.** Droplet diameter under different force fields.


On the other hand, an arrangement should be taken as soon as possible by increasing the relative velocity between droplet and packing to enhance the collision of droplet-packing. Thus, an impinging stream distributor in an impinging stream rotating packed bed (IS-RPB) [18] has been put forward to increase the droplet initial velocity from the liquid distributor. A novel RPB called split packing RPB [19,20] (SP-RPB) can enhance the relative velocity in the packing zone by operating two separate rotating packings in opposing directions, which promotes a large surface area and makes a great contribution to surface renewal and mass transfer in the process of H2S selective absorption into MDEA solution.

#### **4. Conclusions**

A 3D CFD Eulerian–Lagrangian approach has been built by introducing a droplet breakup and coalescence model to investigate the droplet characteristics of RPB in H2S selective absorption into MDEA solution. Droplet characteristics such as droplet velocity, average residence time and average diameter in RPB have been analyzed by diagrams and correlations, which are compared with available experimental data in the literature [14,17]. The results show that the velocity increases with increasing rotating speed and radial position, but the opposite conclusion is made on the average residence time. Specially, in the end zone, a phenomenon called "end effect" in mass transfer intensification has been observed and can be illuminated by droplet back-mixing. A correlation on droplet velocity has been deduced in Equation (33) mainly associated with rotating speed and radial position rather than initial droplet velocity. Under the condition of 900 rpm, a short average residence time 0.039–0.085 s in RPB has been recommended for H2S selective absorption into MDEA solution. This is because the reaction and mass transfer rate of H2S are both instantaneously fast compared with CO2, thus a short average residence time allows for efficient selective absorption between H2S and CO2. When the initial droplet velocity and rotating speed increase, the average droplet diameter decreases inordinately. However, the initial droplet diameter has a restriction (1 mm) to be captured and broken by the packing size under the simulation conditions. Furthermore, conclusions are made that the rotating speed determines the minimum droplet diameter and that a packing length in the radial direction is needed to meet droplet breakup completely. In addition, a balance between shear force and surface tension on the droplet indicates an appropriate rotating speed. A correlation (Equation (36)) on droplet diameter is obtained considering the effect of rotating speed, radial position and fluid density.

The simulation results indicate this CFD approach has the capability in describing droplet characteristics, and an investigation on principle in RPB has been conducted based on these results. Processing in RPB can be regarded as an intensification of a traditional separation device in gravity by replacing the term "*g*" with "ω2*R*" and some coupled field with centrifugal force like magnetic, electric and microwaves, which will inspire further potential in RPB. Finally, IS-RPB and SP-RPB as novel equipment have made a great contribution to surface renewal and mass transfer in the process of H2S selective absorption into MDEA solution through increasing the relative velocity and collision between droplets and packing.

**Author Contributions:** All authors contributed to the manuscript. Z.W. proposed the main idea; X.W. performed the CFD simulation; T.Y. wrote the manuscript; S.W. improved the manuscript; Z.L. provided many key suggestions; X.D. collected and analyzed data.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclatures**


#### **Latin symbols**

σ


#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **A Computational Fluid Dynamics Approach for the Modeling of Gas Separation in Membrane Modules**

#### **Salman Qadir, Arshad Hussain and Muhammad Ahsan \***

School of Chemical and Materials Engineering (SCME), National University of Sciences and Technology, (NUST), Islamabad 44000, Pakistan

**\*** Correspondence: ahsan@scme.nust.edu.pk; Tel.: +92-51-9085-5125

Received: 1 May 2019; Accepted: 28 June 2019; Published: 3 July 2019

**Abstract:** Natural gas demand has increased rapidly across the globe in the last decade, and it is set to play an important role in meeting future energy requirements. Natural gas is mainly produced from fossil fuel and is a side product of crude oil produced beneath the earth's crust. Materials hazardous to the environment, like CO2, H2S, and C2H4, are present in raw natural gas. Therefore, purification of the gaseous mixture is required for use in different industrial applications. A comprehensive computational fluid dynamics (CFD) model was proposed to perform the separation of natural gas from other gases using membrane modules. The CFD technique was utilized to estimate gas flow variations in membrane modules for gas separation. CFD was applied to different membrane modules to study gas transport through the membrane and flux, and to separate the binary gas mixtures. The different parameters of membrane modules, like feed and permeate pressure, module length, and membrane thickness, have been investigated successfully. CFD allows changing the specifications of membrane modules to better configure the simulation results. It was concluded that in a membrane module with increasing feed pressure, the pressure gradient also increased, which resulted in higher flux, higher permeation, and maximum purity of the permeate. Due to the high purity of the gaseous product in the permeate, the concentration polarization effect was determined to be negligible. The results obtained from the proposed CFD approach were verified by comparing with the values available in the literature.

**Keywords:** computational fluid dynamics; membrane module; gas separation; concentration polarization

#### **1. Introduction**

Natural gas is considered one of the significant fossil fuels. It is found in subsurface reservoirs and mostly produced as a byproduct of oil production. The demand for natural gas has seen a considerable rise in recent years [1,2]. As per the reports of the U.S. Energy Information Administration (EIA), global natural gas consumption is increasing 4% annually. It is estimated that 10% of the world power sector depends on natural gas [3]. By 2040, power production and industrial usage are expected to be 73% dependent on natural gas. Natural gas consists of several chemical species, including methane (CH4), ethane, propane, butane, water vapor, nitrogen, and acidic gases such as carbon dioxide (CO2) and hydrogen sulfide (H2S). The species comprising natural gas, like CO2, N2, water vapor, and H2S, are considered impurities [4]. The different concentrations of impurities in natural gas is in the range 4–50%, depending on the reservoir. The presence of these impurities can significantly affect pipelines, in terms of corrosion, and raise health and safety concerns [5]. Therefore, typical pipeline specifications usually mandate that the concentration of carbon dioxide in natural gas not exceed 2–5 volume percent, making it necessary to treat natural gas and remove the impurities before it is transported. In recent

years, the membrane separation process has been widely used because it is more economical compared to other methods [6].

Membrane gas separation can remove unnecessary species from a gas mixture. The membrane allows only the desired components of a mixture to pass through because of its selectivity [7]. It is a widely used process due to its reliability, separation performance, low maintenance, and easy operation [8]. Moreover, membrane technology is used for the separation of various fuel mixtures in different industries as a result of economic competitiveness and other current demanding situations related to competitive environments [9]. It has many applications in industrial sectors, hydrogen recovery from ammonia, hydrogen recovery in refineries, air separation for oxygen purification, sour gas treatment, and carbon dioxide removal from natural gas [10]. These processes are integrated with big industrial units to perform specific industrial operations. It is estimated that the use of membrane gas separation will increase substantially in 2020 [11]. It is an essential fact that the use of this technology will decrease the unit operation cost of gas separation and reduce the environmental hazards [12].

The different mathematical models for the separation of gaseous mixtures using a membrane were developed using altered assumptions [13]. Recently, it was found that the most commonly available commercial membrane modules for gas separation were hollow fiber and spiral wound membrane modules [14,15]. Many researchers have investigated flow behavior with different module arrangements for gas separation and desalination processes [16,17]. Alrehili's study [18] showed the different arrangements of fibers with parallel feed channels that made a hybrid membrane module. The simulation results of the hybrid module gave better membrane flux for both spiral wound and hollow fiber membrane module configurations. Ahsan [19] applied computational fluid dynamics (CFD) modeling to gas separation using a finite difference method in polymeric membranes. Saeed [20] described laminar flow behavior in spacers with narrow channels, and mass transfer coefficient calculated with different wire spacing. Karode [21] determined that pressure dropped with bounding surfaces in a rectangular channel. The effect of shear stresses on both sides of the membrane was also observed.

Other researchers developed a mathematical algorithm for flow behavior and membrane surfaces, and investigated the concentration polarization phenomena in gas separation processes [22]. The transfer of CO2 gas molecules through the membrane increased due to higher flux on the feed side, but rejected molecules of other gases that then accumulated on the membrane surface [23]. For this reason, concentration polarization occurs in membrane processes. Mourgues [24] showed the effect of concentration polarization on membrane separation processes for both counter-current and co-current patterns. The most important factors were analyzed based on feed pressure, permeability, and selectivity of the mixture.

Several researchers developed membrane processes to study the influence of concentration polarization on the feed side. Ahsan and Hussein [25], in their CFD model, studied energy transfer phenomena in the membrane using permeation flux. Coroneo [26] developed a three-dimensional single-tube membrane module to define flux, based on Sievert's law, considering both membrane sides. Recently, a non-isothermal model was developed to evaluate the effect of temperature on permeance [27]. In another study, Chen [28] considered co-current and counter-current flow patterns of the membrane process at different operating conditions using COMSOL Multiphysics 4.0a software.

Another study looked at plug flow and perfect mixing channel, commonly used for modeling membrane gas separation, and reported the fluid behavior in permeate [29]. Flat sheet membrane modules were widely used to evaluate membrane performance. The incompressible Navier–Stokes model was used to improve the fluid chamber while the solid stress–strain model was used to enhance the mechanical performance of the module [30].

In this study, the CFD technique was used to solve the model equations. The permeability of the membrane was measured by introducing the species of interest into the feed gas. The effect of gas flow profiles on gas separation in the membrane modules is reported. A three-dimensional (3D) model was

established using CFD simulation in COMSOL Multiphysics software. The geometry of the flat sheet and spiral wound membrane modules were determined, and co-current and counter-current models were used to find the flow profiles. The effect of molar flux on species transported en mass through the membrane was considered. The binary gas mixtures CO2/CH4 and CH4/C2H6 were used for separation in the flat sheet and spiral wound membrane modules, respectively. The investigated parameters, permeability, feed pressure, permeate pressure, and feed gas concentration, were compared with the data available in the literature.

#### **2. Numerical Methods**

The CFD technique (COMSOL Multiphysics®4.3, COMSOL, Inc., Burlington, MA, USA, 2012) is used in simulations looking at flow profiles in the membrane module. In this study, the flat sheet and spiral wound membrane modules were used for the simulations. The data used for the simulations were taken from the literature and are shown in Table 1.



The numerical simulations were performed using the COMSOL Multiphysics® package (COMSOL Multiphysics®4.3, COMSOL, Inc., Burlington, MA, USA, 2012). The spiral wound and flat sheet membranes were developed with 3D axisymmetric geometry. Fick's law of permeation was used for the main transport of diluted species through the membrane. The following assumptions were made [13].

#### *2.1. Assumptions*


#### *2.2. Mathematical Modeling of Mass Transport*

The model accounts for diffusion transport in a unit cell of the structure sketched in Figure 1. The unit cell is a small part of the membrane that is representative of the whole membrane. In this model, the initial flux in the membrane was studied, and corresponded to the largest difference in concentration between the two chambers on different sides of the membrane. The supporting structure in the membrane consisted of a mesh structure made up of a dense and rigid polymeric material. This system was used to classify the permeability properties of a membrane to certain species.

**Figure 1.** Membrane model unit.

#### 2.2.1. Convection and Diffusion Model

In the convection and diffusion model, the mass transport equation contained both phenomena. In chemical engineering, the convection and diffusion model is applicable for mass transport, which is described by Fick's Law. In the 19th century, Fick gave the simplest definition of diffusion for mass transport. The mass flux of any species is directly proportional to the concentration gradient due to diffusion. The rate of change of concentration at a point in space is proportional to the second derivative of concentration with space;

$$\frac{\partial \mathbf{C\_{a}}}{\partial \mathbf{t}} + \mathbf{v\_{x}} \frac{\partial \mathbf{C\_{a}}}{\partial \mathbf{x}} + \mathbf{v\_{y}} \frac{\partial \mathbf{C\_{a}}}{\partial \mathbf{y}} + \mathbf{v\_{z}} \frac{\partial \mathbf{C\_{a}}}{\partial \mathbf{z}} = \mathbf{D\_{AB}} \left[ \frac{\partial^{2} \mathbf{C\_{a}}}{\partial \mathbf{x}} + \frac{\partial^{2} \mathbf{C\_{a}}}{\partial \mathbf{y}} + \frac{\partial^{2} \mathbf{C\_{a}}}{\partial \mathbf{z}} \right] + \mathbf{R\_{a}} \tag{1}$$

Convective transport accounts for the bulk flow of fluid due to velocity v. This term v can be solved analytically or by solving the momentum equation with the mass balance equation. All these expressions include time (t), and spatial and velocity components are used for mass balance with convection. The vx, vy, and vz are velocity fields in three dimensions. Ra describes the rate of reaction and it should be equal to zero as no reaction is involved in gas separation in membrane modules.

#### 2.2.2. Diffusion Term

The chemical species is transferred from high to low regions due to diffusion, and becomes a mass transfer phenomenon with time and space. The chemical species is dissolved in a solvent or any gas mixture, for example as oxygen enrichment in the air. The evolution of species mass transfer depends on its concentration with respect to time and space. The gradient for diffusion occurs as a result of the motion of the molecules. Due to the kinetic energy of molecules, they can collide with each other in random directions. Then, flux Na for diffusion can be written as;

$$\mathbf{N\_{a}} = \ -\mathbf{D\_{AB}} \left[ \frac{\partial^2 \mathbf{C\_{a}}}{\partial \mathbf{x}} + \frac{\partial^2 \mathbf{C\_{a}}}{\partial \mathbf{y}} + \frac{\partial^2 \mathbf{C\_{a}}}{\partial \mathbf{z}} \right] \tag{2}$$

where mass transfer occurs as a result of diffusion DAB.

#### 2.2.3. Membrane Model

The unit cell is a small part of the membrane module that demonstrates the whole membrane module. The steady-state diffusion equation, which shows mass transport in a model with diffusion, can be written as;

$$\mathrm{D}\_{\mathrm{AB}} \left[ \frac{\partial^2 \mathrm{C}\_{\mathrm{a}}}{\partial \mathbf{x}} + \frac{\partial^2 \mathrm{C}\_{\mathrm{a}}}{\partial \mathbf{y}} + \frac{\partial^2 \mathrm{C}\_{\mathrm{a}}}{\partial \mathbf{z}} \right] = 0 \tag{3}$$

It can also be represented in terms of a nebula operator;

$$\mathbf{D} \cdot \mathbf{V}^2 \mathbf{C}\_4 = 0 \tag{4}$$

where Ca denotes concentration (mole/m3) and DAB is the diffusion coefficient of the diffusing species (m2/s). All boundaries are considered to be insulating

$$\mathrm{D}\_{\mathrm{AB}} \left[ \frac{\partial \mathcal{C}\_{\mathrm{a}}}{\partial \mathbf{x}} + \frac{\partial \mathcal{C}\_{\mathrm{a}}}{\partial \mathbf{y}} + \frac{\partial \mathcal{C} \mathbf{a}}{\partial \mathbf{z}} \right] = \mathbf{0} \tag{5}$$

$$\mathbf{D} \,\nabla \mathbf{C}\_{\mathbf{a}} \cdot \mathbf{n} = 0 \tag{6}$$

The two faces are applied as boundary conditions where the concentration of the two components is set from high to low. Boundary Condition 1 (B.C.1) is considered to be high concentrations on the feed side and Boundary Condition 2 (B.C.2) is represented on the reject side in low concentrations.

B.C.1

$$\mathbf{C} = \mathbf{C}\_0 \tag{7}$$

B.C.2

$$\mathbf{C} = \mathbf{C}\_{0,1} \tag{8}$$

#### 2.2.4. Membrane Flux

Diffusion through the membrane was represented in terms of the initial boundary condition C0 and the final boundary condition C0,1.

$$\mathbf{N\_{4}} = \frac{\mathbf{D}}{\delta} (\mathbf{C\_{0}} - \mathbf{C\_{0,1}}) \tag{9}$$

where <sup>D</sup> <sup>δ</sup> is a barrier with a corresponding thickness.

The permeability (P) of gas can be defined as a product of the diffusion coefficient (D) and solubility coefficient (S),

$$\mathbf{P} = \mathbf{D} \cdot \mathbf{S}.\tag{10}$$

The relationship between partial pressure and concentration is defined as

$$\mathbf{C} = \mathbf{S} \cdot \mathbf{p}.\tag{11}$$

The flux through the membrane can be presented in terms of permeability and partial pressure,

$$\mathbf{N\_{a}} = \frac{\frac{\mathbf{p}}{\xi}}{\delta} (\mathbf{S} \cdot \mathbf{p} \mathbf{f} - \mathbf{S} \cdot \mathbf{p} \mathbf{h}) \tag{12}$$

where feed pressure pf and permeate pressure ph are used for calculating gradient across the membrane. Solubility is obtained using a ratio of the concentration gradient and partial pressure difference for the binary mixture. The solubility S can be calculated as

$$\mathbf{S} = \frac{\Delta \mathbf{p}}{\Delta \mathbf{C}}.\tag{13}$$

This membrane model was used to study gas separation in flat sheet and spiral wound membrane modules.

#### *2.3. Geometry*

A sketch of the geometries of flat sheet and spiral wound membrane modules is shown in Figure 2.

**Figure 2.** Schematic diagram of the flat sheet and spiral wound membrane module.

#### *2.4. Meshing*

The subdomains are often called elements or cells, and the collection of all elements or cells is called a mesh (Figure 3). In this process, a physics-controlled mesh was applied. An extremely fine grid element was used, while a further increase did not affect the model results. The mesh consisted of 1,324,604 domain elements, and 46,750 boundary elements, and 900 edge elements were used for the flat sheet membrane module. The mesh consisted of 3,257,454 domain elements, 413,420 boundary elements, and 3439 edge elements used for the spiral wound membrane module.

**Figure 3.** Meshing of membrane module. (**a**) flat sheet; (**b**) spiral wound.

#### **3. Results**

CFD analysis was performed for the flat sheet and spiral wound membrane modules. The simulation was carried out using COMSOL Multiphysics software 5.2a. Fluxes of gases and concentration variations were found across the length of membrane modules. Various parameters like feed pressure, concentrations across the length of the module, permeability, and feed were considered in this study.

#### *3.1. CH4*/*C2H6 Separation*

The flat sheet membrane module is usually used for pilot plant testing or lab scale testing. The separation of methane and ethane was considered for the flat sheet membrane module in this study. The simulation of a flat sheet membrane module was performed to check the concentration change on the reject and permeate sides. The feed gas entered the membrane module and the permeate collected at the bottom of the module. A cross-flow model with specific boundary conditions was used. The boundary conditions included a high feed concentration on the right side and low feed concentration on the left side. Figure 4a shows the concentration variation of CH4 on the feed and permeate sides. The contour shows the concentration variation of gas on both sides of the membrane. The slice centration shows the gas variations in the center of the module. Figure 4b shows the differing concentration variation from the feed side to the permeate side. The gas moved through the flat sheet module and permeate collected at the bottom.

Figure 4c shows the concentration gradient of CH4 present in the flat sheet membrane module. Streamlines were used for concentration gradient representation in Figure 4c. The lines moving from high to low and passing through the membrane represent the gas diffusion through the membrane. The membrane is located in the center of the module and the permeate is shown on the bottom surface. The results verified that a gradient was present, and that gas diffused through the membrane. The flux was calculated for the given parameters using COMSOL Multiphysics software. The contour in Figure 4d indicates the flux variations in the flat sheet membrane module. The color bar shows the flux magnitude in the module. The flux calculated in the module can be explained by Fick's law.

**Figure 4.** (**a**) CH4 gas concentration variation in the flat sheet membrane module; (**b**) slice shows CH4 gas variation in the center of the module; (**c**) line shows the concentration gradient variation in the membrane module and (**d**) diffusive flux variation for CH4 in a flat sheet membrane module.

#### *3.2. CH4*/*CO2 Separation*

The simulation of a spiral wound membrane module was presented to show the concentration of feed CO2 on the reject and permeate sides. The feed gas entered the membrane module from the central tube and moved through to the end. The module shows that the concentration of CO2 changed throughout Figure 5a. The high concentration was applied to the first wrapped sheet as a boundary condition, permeate collected in the permeate channel, and a low concentration was applied to the second plate of the membrane module.

The contour shows that the gas concentration varied along the membrane on the permeate and reject sides. Figure 5b shows the slice concentration of CO2 in a spiral wound membrane module. The result shows the concentration variation in the center of the module. It verifies that the gas diffused through the membrane from the feed channel and collected in the permeate channel. Figure 5c shows the concentration gradient present in the spiral wound membrane module for CH4/CO2 separation. The arrows show that the gas passes through the membrane from the feed side to the permeate collector of the module. The results prove that a gradient was present, and that the gas diffused through the membrane. Figure 5d shows the diffusive flux magnitude in the spiral wound membrane module for CO2. The diffusive flux magnitude in the membrane module was obtained using concentration and permeability parameters. The flux through the membrane was calculated in the gas separation process and depended on gas diffusion in the perm-selective membrane, due to the pressure and concentration gradient. It illustated that mass transport occurred through the membrane in the spiral wound membrane module.

**Figure 5.** (**a**) CO2 gas concentration variation in the spiral wound membrane module; (**b**) slice shows CO4 gas variation in the center of the module; (**c**) line shows concentration gradient variation in the membrane module and (**d**) diffusive flux variation for CO4 in the spiral wound membrane module.

#### *3.3. Parametric Study of the Spiral Wound Membrane Module*

Feed Pressure and Permeate Pressure

The membrane length was considered for the spiral wound membrane module, and the results were verified via the literature. Figure 6 shows that permeate pressure was a function of membrane length for different values. Basically, the permeate pressure was less than the feed pressure. Therefore, it was necessary to calculate the changes in permeate pressure according to length. Increasing the membrane module's length decreased the permeate pressure because of the gradient changes throughout the module from the feed side to the reject side. Therefore, less mass transfer occurred as a result of less gradient in the module. It was observed that feed pressure had a reverse effect on permeate pressure. An increase in the feed pressure produced more gradient across the membrane, which resulted in higher mass transfer through the membrane. Figure 7 represents the feed pressure variation with methane recovery in the spiral wound module. The results show that higher amounts of methane permeate were obtained in the end because a high gradient was produced, and gas diffusion through the membrane was very high.

**Figure 6.** Permeate pressure variation with module length in a spiral wound membrane module compared with published values [31].

**Figure 7.** CH4 gas variation with feed pressure in the spiral wound membrane module compared with published values [31].

#### *3.4. Parametric Study for the Flat Sheet Membrane Module*

#### Feed Pressure and Length

It was observed that the permeate concentration changed along the length of the module. The more permeable components passed through the membrane while other components were left as rejects. The CH4 gas was more suitable for permeation through the membrane. Figure 8 demonstrates that the mole fraction of CH4 increased on the permeate side across the length of the module and decreased on the reject side, with respect to length. The ethane gas increased on the reject side because the permeation of ethane through the membrane was very low. Therefore, the maximum possible purity of the gas was achieved.

Figure 9 shows that increasing the feed pressure causes the gradient to increase. Due to high feed pressure, the permeation of components through the membrane increased when the driving force increased. Therefore, it was estimated that the higher the feed pressure the better the separation. This was also true for methane recovery on the permeate side because increasing feed pressure resulted in more pressure gradient and turbulence for mass transfer. Therefore, the gas passed through the membrane. The permeate values calculated by the numerical model for the spiral wound and flat sheet membrane modules are shown in Table 2 for comparison with published data. A maximum difference of 10.8% and 8.7% was found for the spiral wound and flat sheet membrane modules, respectively.

**Figure 8.** Mole fraction of methane with feed pressure in the flat sheet membrane module compared with published values [32].

**Figure 9.** Mole fraction of methane change with a length of the flat sheet membrane module compared with published values [32].


**Table 2.** Comparison of the model results with the published data for CH4 (permeate) at different values of selectivity.

#### **4. Discussion**

Computational fluid dynamics simulations were used to study mass transport in spiral wound and flat sheet membrane modules. Three-dimensional geometrics were considered to analyze the membrane-based gas separation for the binary gas mixture. Counter-current and cross-flow membrane models were used to verify the results. The simple Fick's law was used to explain mass flux transport of the binary gas mixture through the membrane. The membrane model was defined using COMSOL Multiphysics software for the separation of the binary mixture. The membrane model was applied as a thin diffusion barrier, which allowed certain species to pass through the membrane. The effect of molar flux for species in mass transport through a membrane was considered. In the present model, the different parameters were investigated to verify the models in the literature. Different gas mixtures, like CO2/CH4 and CH4/C2H6, were investigated in different membrane modules.

A spiral wound membrane module was discussed for CH4/CO2 membrane separation. The module geometry consisted of three flat sheets that wrapped around a central tube. The membrane collector was present at the center of the feed and reject sides. The cross-flow model was applied to verify the literature results. In the spiral wound membrane module, the increase in membrane length showed a considerable decrease in concentration polarization. When the length of the membrane increased, an increase in the residual mole fraction was observed. An increase in the permeate purity was also observed, which indicated that the concentration polarization was negligible.

A flat sheet membrane module was considered for the separation of CH4/C2H6. The model showed a variation of membrane results that have been discussed. This was novel work describing the flow profiles of gases in different membrane modules. The contour also showed the total flux, concentration gradient, and diffusive flux magnitude of different membrane modules. The investigated parameters were compared with published results. It was observed that in a flat sheet membrane module with increasing feed pressure, the pressure gradient also increased, which resulted in higher flux, higher permeation, and maximum purity of the permeate. The concentration polarization was observed to be negligible. Furthermore, with increasing module length in the flat sheet membrane module, a decrease in concentration polarization was observed because the increase in the module length resulted in more permeation of the desired component and an increase in permeate purity. In the spiral wound membrane module, increasing the membrane length resulted in a considerable decrease in concentration polarization. When the length of the membrane was increased, an increase in the residual mole fraction was observed. An increase in the permeate purity was also observed, which indicated that the concentration polarization was negligible.

#### **5. Conclusions**

The membrane modules were modeled using CFD to obtain the maximum possible value of the desired gas in the permeate. The effect of concentration polarization on gas separation performance was negligible. Different parameters were studied in the membrane modules, such as feed pressure, module length, permeate pressure, and feed concentration. Increasing feed pressure in the membrane modules caused the pressure gradient to increase. Therefore, maximum mass transfer occurred through

the membrane. Moreover, the length was considered to show gas variation in the permeate. In the end, the contours of gas flow profiles in the membrane modules were successfully reported. Modeling predictions were compared with the published data and validated, and it was found that there was good agreement between them for the different values of selectivity. The comparison indicated that the membrane modules were very efficient in terms of the separation of the desired gas at a higher pressure gradient.

**Author Contributions:** A.H. proposed the main idea and methodology; S.Q. performed the CFD analysis and wrote the manuscript; M.A. and A.H. provided key suggestions and improved the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **List of symbols**


#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
