Direct Numerical Simulation of Hydrate Dissociation in Homogeneous Porous Media by Applying CFD Method: One Example of CO2 Hydrate

Abstract

Computational fluid method (CFD) is popular in either large-scale or meso-scale simulations. One example is to establish a new pore-scale (m~μm) model of laboratory-scale sediment samples for estimating the dissociation rate of synthesized CO2 hydrate (CO2H) reported by Jeong. It is assumed that CO2H formed homogeneously in spherical pellets. In the bulk flow, concentration and temperature of liquid CO2 in water flow was analyzed by CFD method under high-pressure state. Finite volume method (FVM) were applied in a face-centered packing in unstructured mesh. At the surface of hydrate, a dissociation model has been employed. Surface mass and heat transfer between hydrate and water are both visualized. The initial temperature 253.15K of CO2H pellets dissociated due to ambient warm water flow of 276.15 and 282.15K and fugacity variation, ex. 2.01 and 1.23 MPa. Three tentative cases with porosity 74, 66, and 49% are individually simulated in this study. In the calculation, periodic conditions are imposed at each surface of packing. Numerical results of this work show good agreement with Nihous’ model. Kinetic modeling by using 3D unstructured mesh and CFD scheme seems a simple tool, and could be easily extended to determine complex phenomena coupled with momentum, mass and heat transfer in the sediment samples.

1. Introduction

Computational fluid method (CFD) is popular in either large-scale or meso-scale simulations. One example is to establish a new pore-scale (m~μm) model of laboratory-scale sediment samples for estimating the dissociation rate of synthesized CO₂ hydrate (CO₂H) reported by [1]. To decrease the CO₂ concentration in the air, carbon dioxide capture and storage (CCS) is regarded to be an effective way. One concept of CCS is to store CO₂ in gas hydrate in sub-seabed geological formation, as was illustrated by [6]. Besides, many studies about the formation and dissociation of CO₂ hydrate (CO₂H) while stored in the deep ocean or geologic sediment have been introduced. In particular, flow and transport in sediment is multidisciplinary science including the recovery of oil, groundwater hydrology and CO₂ sequestration. It reported the measurements of the dissociation rate of well-characterized, laboratory-synthesized carbon dioxide hydrates in an open-ocean seafloor [5]. The pore effect in the phase equilibrium mainly due to the water activity change was discussed in [7]. The reactive transport at the pore-scale to estimate realistic reaction rates in natural sediments was discussed in [3]. This result can be used to inform continuum scale models and analyze the processes that lead to rate discrepancies in field applications. Pore-scale model is applied to examine engineered fluids [4]. Unstructured mesh is well suited to pore-scale modeling because of adaptive sizing of target unit with high mesh resolution and the ability to handle complicated geometries [17, 18]. Particularly, it can easily be coupled with computational fluid dynamics (CFD) methods, such as finite volume method (FVM) or finite element method (FEM). Unstructured tetrahedral mesh used to define the pore structure is discussed in [19]. Another case includes a numerical simulation of laminar flow based on FVM with unstructured meshes was used to solve the incompressible, steady Navier-Stokes equations through a cluster of metal idealized pores by [20].

The objective of this work is to develop a new pore-scale model for estimating the dissociation rate of CO₂H in homogeneous porous media. To cooperate with molecular simulation and field-scale simulators, we aim at establishing pore-scale modeling to analyze the simultaneous kinetic process of CO₂H dissociation due to non-equilibrium states. Major assumptions in this study are listed as below:

1. Only dissociation occurred at the surface, no any formation occurred immediately with dissociation.

2. CO₂ dissociated at the surface is assumed to be dissolved into liquid water totally without considering the gas nucleated.

3. The surface structure does not collapse with the dissociation of CO₂H at the surface of pellets.

4. Homogeneous face-centered packing of multi-CO₂H pellets.

5. Single phase flow coupled with mass, heat, and momentum transfers.

2. Dissociation modeling at the surface

In this study, the dissociation flux ($F_1$) is assumed to be proportional to the driving force, the free energy difference ($∆\mu$) introduced by [6], presented as

$$F_1=k_{\mathrm{bl}}\text{RTln}\left(\frac{C_{\mathit{\text{Hsol}}}}{C_I}\right)$$E1

where $k_{\mathit{bl}}$ is the rate constant $\left[{\mathrm{mol}}^2{\mathrm{J}}^{-1}{\mathrm{s}}^{-1}{\mathrm{m}}^{-2}\right]$ of dissociation. According to [21], $k_{\mathit{bl}}$ is listed as below:

$$k_{\mathit{bl}}=\text{exp}(-\frac{\mathrm{11,729}}{T}+26.398)$$E2

where $C_{\mathit{\text{Hsol}}}$ is the mole fraction of CO₂ in the aqueous solution at equilibrium state with hydrate, and $C_I$ means surface concentration in the ambient aqueous solution at the surface of the hydrate $C_I$.

3. Basic transport equations

Flow in the porous media around CO₂H is governed by the continuity and the Navier-Stoke’s equations. The advection-diffusion equations of non-conservative type for mass concentration C and temperature T are also considered.

$$\nabla \cdot \mathbf{u}=0$$E3

$$\frac{\partial \mathbf{u}}{\partial t}+\nabla \cdot \left(\mathbf{uu}\right)=-\nabla P+\frac{1}{\text{Re}}\nabla \cdot \left[\nabla \mathbf{u}+{\left(\nabla \mathbf{u}\right)}^T\right]+\frac{{\rho }_w}{{\mathrm{Fn}}^2}g$$E4

$$\frac{\partial C}{\partial t}+\mathbf{u}\cdot \nabla C=\frac{1}{\text{ReSc}}{\nabla }^{2}C$$E5

$$\frac{\partial T}{\partial t}+\mathbf{u}\cdot \nabla T=\frac{1}{\text{RePr}}{\nabla }^{2}T$$E6

where the viscosity, diffusivity, and thermal conductivity of pure water are included in dimensionless parameters such as the Reynolds number, the Schmit number, and the Prandtl number, which are interpolated as functions of temperature and are updated at every computational time step as summarized in Table 1. U and d (=0.001 m) are the velocity of inflow and diameter of hydrate pellet.

No.	Function	Definition
1	$D=7.4\times {10}^{-12}\frac{{\left(\phi M_B\right)}^{1/2}T}{{\eta }_L{V}_A^{0.6}}$	D: diffusion coefficient of CO2 in water; $\phi$(=2.6): association parameter for the solvent water; $M_B\left(=18\mathrm{g}{\mathrm{mol}}^{-1}\right)$: molecular weight of water; $V_A\left(=3.4\times {10}^{-5}{\mathrm{m}}^3{\mathrm{mol}}^{-1}\right)$: molar volume of CO2; ${\eta }_L$[$\mathrm{mPa}\bullet \mathrm{s}$]: viscosity of water
2	${\nu }_L=\left(8.8286\times {10}^{-10}\right)T^2-\left(5.3886\times {10}^{-7}\right)T+8.314\times {10}^{-5}$	${\nu }_L\left[\mathrm{m}{\mathrm{s}}^{-2}\right]$: kinematic viscosity
3	${\lambda }_L=487.85\text{ln}\left(T\right)-2173.8$	${\lambda }_L\left[\mathrm{WK}{\mathrm{m}}^{-1}\right]$: heat conductivity of water
4	${\alpha }_L=\frac{{\lambda }_L}{{\rho }_W{C}_p}$	${\alpha }_L\left[\mathrm{m}{\mathrm{s}}^{-2}\right]$: the thermal diffusivity of aqueous phase; ${\rho }_W\left(=997.1\frac{\mathrm{kg}}{{\mathrm{m}}^3}\right)$: density of water; $C_p\left(=\mathrm{4,180}\frac{\mathrm{J}}{\mathrm{kg}}\bullet \mathrm{K}\right)$: isobaric specific heat, quoted from “Chemical Engineering Handbook”, Japan (1985)
5	$P_{\mathit{eq}}=\text{exp}\left(\alpha +\frac{\beta }{T}\right)\times {10}^3$	$\alpha$=44.580 and $\beta$ = − 10246.28
6	$C_{\mathit{\text{Hsol}}}=a\bullet \text{exp}\left(b\bullet P\times {10}^{-6}+1.321\times {10}^{-4}T-2.292\times {10}^{-2}\right)\bullet 1.8\times {10}^{-5}$	$C_{\mathit{\text{Hsol}}}\left[\text{mole}\bullet {\mathrm{m}}^{-3}\right]$: solubility of hydrate (275.15 K < T < 281.15 K); $a=0.0016{\left(T-273.15\right)}^{0.9211}$; $b=-0.0199\text{log}\left(T-273.15\right)+0.0942$; by Aya et al. [10], Yang et al. [11], and Servio and Englezos [12]

Table 1.

Parameters used in this study.

4. Mass transfer

To rewrite Eq. (1), the flux at the surface of the hydrate can be discretized as

$$k_{\mathit{bl}}\mathit{RT}\text{ln}\left(\frac{C_{\mathit{\text{Hsol}}}}{C_I}\right)=D\nabla C=D\frac{C_I-C^{\prime }}{h_I}$$E7

where $C_I$ is the varying surface concentration calculated locally at each surface cell, $C^{\prime }$ is the centroid concentration, and $h_I$ is the thickness of centroid surface cell, as shown in Figure 1.

Figure 1.

Schematic image of discretized surface concentration.

5. Heat transfer

The equation of energy balance at the surface of CO₂H is given by

$${\dot{Q}}_H+{\lambda }_H\nabla T_H={\lambda }_L\nabla T_L$$E8

where ${\dot{Q}}_H$ ($=H_L{F}_1$, where $H_L$ is the latent heat of hydrate dissociation) is the dissociation heat transferred to the CO₂H, ${\lambda }_H$ is the thermal conductivity of hydrate. Dissociation heat per mole of hydrate, $H_L,$ is interpolated from [2] as

$$H_L=\mathrm{207,917}-530.97\times T_I$$E9

where $T_I$ is the surface temperature. Then, we have

$$T_I=\frac{{\lambda }_L{h}_H{T}_L+{\lambda }_H{h}_L{T}_H-\mathrm{207,917}\bullet F_{\mathit{cal}}{h}_L{h}_H}{{\lambda }_L{h}_H+{\lambda }_H{h}_L-530.97\bullet F_{\mathit{cal}}{h}_L{h}_H}$$E10

where $T_L$ and $T_H$ are the temperatures defined at the centroids cell in the aqueous phase and solid hydrate, respectively; $h_L$ and $h_H$ are half widths of centroid in the aqueous phase and solid hydrate, respectively. Besides, the temperature in the pellet is calculated by using the heat conductivity equation.

$$\frac{\partial T}{\partial t}={\alpha }_H\frac{{\partial }^{2}T}{\partial x^2}$$E11

where ${\alpha }_H=\frac{{\lambda }_H}{{\rho }_H{C}_p^{\prime }}$ is the thermal diffusivity of CO₂H. These relative properties of CO₂H are quoted from [9], the thermal diffusivity of aqueous phase (${\alpha }_H$) of $1.38\times {10}^{-7}\mathrm{m}{\mathrm{s}}^{-2}$, the heat capacity of hydrate ($C_p^{\prime }$) of $2080.0\mathrm{J}{\mathrm{K}}^{-1}{\mathrm{kg}}^{-1}$, and the thermal conductivity (${\lambda }_H$) of $0.324\mathrm{W}{\mathrm{K}}^{-1}{\mathrm{m}}^{-1}$. The density of CO₂H (${\rho }_H$) is given as $1116.8\mathrm{kg}{\mathrm{m}}^{-3}$.

6. Computational conditions

Two types of cells, tetrahedrons and triangular prisms, are applied in the present unstructured grid system, as introduced in Figure 2. In detailed, the surface of hydrate uses prism. Both the flow field and inside the pellet are tetrahedral meshes. Upward is the inflow where initially the uniform velocity profile is adopted. Prism mesh and no-slip condition are imposed at the surface of the pellet. To analyze more detailed mass and heat transfer simulatneously, one cell-layer of the prisms that attached to the CO₂H surface is divided into at least five very thin layers as referred in [8] for high Prandtl or Schmidt number. The basic parameters of computation are denoted in Table 1. The initial values of dimensionless parameters are listed in Table 2 at the temperatures from 276.15 to 283.15 K. Reynolds number, Schmidt number, and Prandtl number function of the temperature or pressure are listed in Table 2. The minimum grid size of this computational model is listed in Table 3. Lm, Lc, and LT are the applied mesh thicknesses. δm, δc, and δT are the thickness of momentum, concentration, and thermal boundary layers, respectively. The relationship between δm, δc, and δT quoted from the theory of flat plate boundary layer is listed below:

$${\delta }_m=\frac{5.48}{\sqrt{\mathit{Re}}}\frac{d}{2}$$E12

$${\delta }_c=\frac{{\delta }_m}{1.026·{\mathit{Sc}}^{1/3}}$$E13

$${\delta }_T=\frac{{\delta }_c}{{\mathit{Pr}}^{1/3}}$$E14

Figure 2.

Description of mesh in unstructured grid system of 67,104 cells. (a) Overview, (b) surface of CO₂H, (c) water, and (d) pellets of CO₂H.

Case	Reynold number	Prandtl number	Froude number	Schmidt number	Porous ratio	Temperature of water (K)	Fugacity of equilibrium (MPa)	Fugacity (MPa)
1	50	10	0.023	755	74%	282.15	3.89	3
2		12		880		276.15	1.77	1
3		10		755	66%	282.15	3.89	3
4		12		880		276.15	1.77	1
5		10		755	49%	282.15	3.89	3
6		12		880		276.15	1.77	1
7	100	10	0.046	755		282.15	3.89	3
8	135		0.062

Table 2.

Calculation conditions of this work.

Cell number (porosity)	Sc*	δm*	δc&T*	Lm	VTL number	Lc&T
53,440 (49%)	755	2.358E-04	2.252E-05	1.0E-05	5	2.0E-06
67,104 (66%)	755	2.345E-03	2.252E-04	6.0E-05	5	1.2E-05
77,432 (74%)	755	2.594E-03	2.276E-04	6.0E-05	5	1.2E-05

Table 3.

The thicknesses of boundary layers, δ_m and δ_c&T, and grid sizes, L_m and L_c&T (unit: meter).

*

Values are quoted from “Chemical Engineering Handbook”, Japan (1985)

To follow [22] of Eq. (14), the boundary layer’s thickness for temperature, ${\mathrm{\delta }}_{\mathrm{T}}$ is assumed as the same size as that for mass concentration, ${\mathrm{\delta }}_{\mathrm{c}}$. For the initial temperature of the CO₂H pellet, T_ini is assumed as a constant value of 253.15 K.

7. Verification

The in-house code originally developed by [7] has been applied to determine the intrinsic dissociation rate of methane hydrate. The numerical results verified by experimental results are successfully used in calculating one pellet of hydrate in a slow flow rate of high pressure without considering the collapse of hydrate and the nucleation of bubbles referring to [6, 20] as well.

8. Results of case study

In this study, cases with porosity of 74, 66, and 49 are individually discretized as face-centered unstructured packing of hydrate in sediment. CO₂H pellets with initial temperature of 253.15 K dissociate due to the variation of driving force, ex. 0.89 and 0.77 MPa, under thermal stimulation of ambient warm water, ex. 282.15 and 276.15 K. Comparative small driving forces selected here is due to the assumption of no surface’s collapse. Computational conditions are listed in Table 2. Result of flux at the surface is the converge value as shown in Figure 3. In Figure 4(a)–(c) at time 0.16(s) show velocity vector of case 1 in three specific sections. In Figure 5, the distributions of concentration at 0.16 s are presented. Slight CO₂ discharges at the surface. Relative temperature distributions are indicated in Figure 6. As time increases, the dissociation heat of CO₂ hydrate results in water temperature drop significantly as shown in Figure 7(a)–(c). Relative concentration distribution in center section is shown in Figure 8. The heat of water conducts to the solid-side pellet rapidly in few seconds, and slow mass transfer at the surface dominates the dissociation rate rather than fast heat transfer at the surface.

Figure 3.

Converge value of case 1.

Figure 4.

Velocity vector in the vicinity of pellets (case 1: T = 282.15 K, T_ini = 253.15 K, Re = 50, Sc = 755, Pr = 10, VTL = 5, and time = 0.16 s). (a) Front section, (b) Center section and (c) Back section.

Figure 5.

Concentration profile in three specific sections of cubic unit (case 1: time = 0.16 s, unit: mole/m² s). (a) Front section, (b) Center section and (c) Back section.

Figure 6.

Temperature profile in three specific sections of cubic unit (case 1, time = 0.16 s, unit: K). (a) Front section, (b) Center section and (c) Back section.

Figure 7.

Temperature versus time in center section of cubic unit (case 1, time = 0.16 s, 0.27 s, and 0.54 s, unit: K). (a) Front section, (b) Center section and (c) Back section.

Figure 8.

Concentration versus time in center section of cubic unit (case 1, time = 0.16 s, 0.27 s, and 0.54 s, unit: mole/m² s). (a) Front section, (b) Center section and (c) Back section.

9. Discussion

To follow the modeling illustrated in [14]:

$$F_3=k_D\left(f_{\mathit{eq}}-f_g\right)=k_{D0}\text{exp}\left(-\frac{\mathrm{\Delta }E}{\mathit{RT}}\right)\left(f_{\mathit{eq}}-f_g\right)$$E15

where $k_{D0}\left[\mathrm{mol}{\mathrm{Pa}}^{-1}{\mathrm{s}}^{-1}{\mathrm{m}}^{-2}\right]$ is the rate constant, $f_g\left(\mathrm{Pa}\right)$ is the fugacity of gaseous CO₂, and $f_{\mathit{eq}}$ is the fugacity of the quadruple equilibrium. They obtained $k_{D0}$ and $∆E$ for CO₂H as $1.83\times {10}^8\mathrm{mol}{\mathrm{Pa}}^{-1}{\mathrm{s}}^{-1}{\mathrm{m}}^{-2}$ and $102.88\mathrm{kJ}{\mathrm{mol}}^{-1}$, respectively, at temperature and pressure ranging from 274.15 to 281.15 K and from 1.4 to 3.3 MPa. However, new modified value of $∆E$, if considered the real case in the ocean quoted from [16], is $96.49\mathrm{kJ}{\mathrm{mol}}^{-1}$. The order of Reynolds number based on the size of a particle, about 16 μm, is calculated as 50. Clarke et al. [28] determined the dissociation rate of CO2H by measuring the nucleated methane gas in V-L state [14]. The comparison of three models is listed in Table 4. The results of dissociation flux are summarized in Figure 9. Higher water temperature induces higher dissociation flux at the surface of hydrate. Data correlated by [14] show much lower level than both Nihous’ model and new model. The numerical results in this work marked as new model in Figure 9 show consistent result compared with Nihous’ model. The dissociation flux in various flow rates in cases 5, 7, and 8 are listed in Figure 10. Here, it is noted that porosity is not considered in both Clarke’s and Nihous’ models, and these two models are only function of temperature and fugacity as presented in Eq. (15). The trend of flux becomes saturated in the figure due to the surface dissociation flux that becomes slow due to the fast mass transfer in bulk flow at Reynolds number over 100.

Item	Modeling	Intrinsic rate of dissociation	Ref.
Clarke’s Model	Eq. (15)	KD0 = 1.83 * 10^8mole/Pa s m; ΔE = 102.89 kJ/mole	[14]
Nihous’ Model	Eq. (15)	KD0 = 1.83 * 10^8mole/Pa s m; ΔE = 96.49 kJ/mole	[16]
New Model	Eq. (1)	Kbl = exp(−11,729/T + 26,398)	[21]

Table 4.

The comparison of three models.

Figure 9.

Results of simulation compared to existing two models.

Figure 10.

Results of flux in cases 5, 7, and 8.

10. Conclusions

The objective of this work is to establish a new pore-scale model for estimating the dissociation rate of CO₂H in laboratory-scale sediment samples. It is assumed that CO₂H formed homogeneously in spherical pellets. In the bulk flow, concentration and temperature of liquid CO₂ in water flow was analyzed by computational fluid dynamics (CFD) method without considering gas nucleation under high-pressure state. In this work, finite volume method (FVM) was applied in a face-centered regular packing in unstructured mesh. At the surface of hydrate, a dissociation model has been employed. Surface mass and heat transfer between hydrate and water are both visualized. The initial temperature 253.15 K of CO₂H pellets dissociated due to ambient warm water flow of 276.15 and 282.15 K and fugacity variation, ex. 2.01 and 1.23 MPa. Three tentative cases with porosity 74, 66, and 49% are individually simulated in this study. In the calculation, periodic conditions are imposed at each surface of packing. Additionally, the flux at CO₂H’s surfaces is compared to Clarke and Bishnoi [13] and Nihous and Masutani [15] ‘s correlations at Reynolds number of 50. Numerical results of this work show good agreement with Nihous’ model. Kinetic modeling by using 3D unstructured mesh of regular cubic unit and CFD scheme seems to be a simple tool to estimate the dissociation rate of CO₂H in laboratory-scale experiments, and could be easily extended to determine complex phenomena coupled with momentum, mass, and heat transfer in the sediment samples.

Acknowledgments

This work was supported by DOIT, Ministry of Science and Technology under contract No. MOST 106-3113-M-002-006. The authors also wish to acknowledge Professor Toru SATO for the valuable advices and guidance.

Nomenclature

C	volumetric molar concentration of CO2 in the ambient water molm−3
CH	volumetric molar concentration of CO2 in the aqueous solution equilibrated with the stable hydrate phase molm−3
C′	volumetric molar concentration of CO2 in water at the centroid of a cell attaching to the hydrate surface molm−3
CI	volumetric molar concentration of CO2 in the ambient aqueous solution at the surface of the hydrate ball molm−3
CX	average molar volumetric concentration of CO2 in the ambient water flow for a given cross section of water flow molm−3
d	diameter of the CO2 hydrate ball [m]
D	diffusion coefficient of CO2 in water ms−2
E	activation energy Jmol−1
F	dissociation rate flux mols−1m−2
feq	fugacity of the quadruple equilibrium Pa.
fg	fugacity of gaseous CO2 Pa
G	molar Gibbs free energy Jmol−1
HL	latent heat of hydrate dissociation Jmol−1
hL	length of the water layer attached to the hydrate surface [m]
kD0	intrinsic dissociate rate constant based on Clarke-Bishnoi model molPa−1s−1m−2
kbl	dissociation rate constant based on new model mol2J−1s−1m−2
L	thickness of computational cell [m]
MB	molecular weight of water gmol−1
P	thermodynamic pressure [Pa]
Peq	quadruple equilibrium pressure for CO2 hydrate as a function of T [Pa]
Q	volumetric flow rate of the ambient water m3s−1
Q̇H	the rate at which the latent heat is transferred to the CO2 hydrate by dissociation Jm−2s−1
R	gas constant, 8.314 JK−1mol−1
T	absolute temperature [K]
TL	temperature at the centroids of a cell in the solid hydrate [K]
TH	temperature at the centroids of a cell in the aqueous phase [K]
u	velocity vectors [m/s]
x	mole fraction of CO2 [-]
xeq	solubility of CO2 in the aqueous solution in equilibrium with the stable hydrate phase [-]
xI	mole fraction of CO2 in the aqueous phase at the surface of the hydrate ball [-]
αL	thermal diffusivity in the aqueous phase ms−2
αH	thermal diffusivity in the hydrate ball ms−2
Δr	thickness of the computational cell [m]
δ	thickness of the boundary layer [m]
Δμ	chemical potential difference Jmol−1
ρw	density of the ambient water kgm−3
φ	the association parameter for the solvent water
ηL	the viscosity of water [mPas]
VA	the molar volume of CO2 m3mol−1
νL	kinematic viscosity of water ms−2
λL	heat conductivity of water WK−1m−1
λH and λL	the heat conductivities in the hydrate and water WK−1m−1