Parameters Variation of Natural Gas Hydrate with Thermal Fluid Dissolution Based on Multi-field Coupling under Pore-Scale Modeling

1. Introduction

Using geothermal energy for mining natural gas hydrates under the seabed is an efficient and low-carbon method[1]. As the fluid injected into the deeper geothermal reservoir, the heated water conveys the heat and is reinjected into the fractured hydrate reservoir[2]. During the hydrate dissolving in the reservoir, the thermal fluid enhances the porosity of structure, thereby augmenting the flow capacity and convective heat transfer capacity of the reservoir[3]. Accordingly, when heat flows through the medium of hydrate, a condition of fully multi-field coupling must be considered with every physical field of thermal-flowing-mechanic-chemical (TFMC)[4]. The simulation must therefore analyze these interactions between the fields. It is necessary to analyze the mutual influence of key parameters between multi-field coupling during simulation.

According to the reviews of literatures, there must be functions or laws to express the coupling ways of each field in TFMC. The coupling ways between flowing and mechanic fields are the Biot equation and Darcy law[5]. The consolidation coefficient was used to describe the harmonic relationship between the solid stress field and the fluid pressure field. The pressure field changes the fluid velocity field by Darcy law. However, the porosity changes greatly due to hydrate dissolution, and the Biot equation is no longer applicable. The coupling ways of flowing and thermal field are the convection heat transfer equation and viscosity-temperature relationship[6]. The change of temperature has little effect on the viscosity, and the flow in the microporous structure is laminar, and the average size of the pores is generally smaller than the thickness of the boundary layer. Therefore, the convection heat transfer equation is directly used to simulate the internal energy transfer in the pores. Its thermal diffusion coefficient should include heat conduction and heat convection, while heat conduction is retained in the solid. The coupling ways between flowing field and chemical field are the Kim-Bishnoi (K-B) equation in pressure form and the kinetics of methane hydrate decomposition[7]. The pressure distribution in the flowing field will affect the chemical reaction rate based on the K-B equation. The energy consumption caused by chemical reactions at the microscale will lead to kinetic energy transfer in the flowing field, but it will not be too obvious. The coupling way of stress field and chemical field is porosity change. Hydrate dissolution will increase porosity and reduce the elastic modulus of the solid skeleton. The elastic modulus decreasing will lead to changes in the pressure field, which in turn affects the reaction rate. The coupling ways of mechanic and thermal field are thermal expansion and effective heat conductivity. When hydrates dissolving, the thermal expansion deformation of the solid structure is very tiny, but the porosity changes greatly, which in turn affects the effective heat conductivity, which has a significant impact on convective heat transfer. The coupling way of thermal and chemical field is K-B equation in thermal form. The reaction rate coefficient of the K-B equation is related to temperature, and the equilibrium pressure is also related to temperature, but the effect is not significant. The chemical reaction has a great influence on the internal energy transfer. The dissolution enthalpy of the chemical reaction will absorb a lot of heat and affect the internal energy diffusion rate[8].In summary, when comprehensively considering the multi-field coupling of TFMC, the unsteady Darcy's law, the unsteady convective heat transfer equation, and the unsteady K-B equation are important governing equations. Therefore, the average permeability, effective thermal conductivity and reaction rate coefficient are critical physical parameters. At microscale simulation, Darcy's law becomes the N-S equation, and the K-B equation is a macroscale equation and is not applicable. Therefore, it is necessary to solve the above key physical parameters through pore-scale simulation of the N-S equation, convective heat transfer equation, and fluid-solid boundary changes, and verify Darcy's law and K-B equation of macro-scale condition.

The lattice Boltzmann method (LBM) is a discretely REV scale simulation method based on the Boltzmann distribution equation. It can not only characterize the effects of microscale simulation, but also has perfect convergence and stability. Using the dual distribution function (DDF-LBM) method, the distribution of velocity and temperature field at the microscale modelling can be more accurately calculated[9]. At the same time, the calculation efficiency is relatively faster and it has better self-stability for the fluid-solid boundary behaviors. Huo et al. used the DDF-LBM to simulate the solid-liquid boundary change during the dissolution of phase change of silicon materials in a square cavity and verified the feasibility of this method under natural thermal convection conditions[10]. He et al. reviewed and introduced the detailed method of LBM simulation of single fluid solid-liquid phase change[11]. The literature listed the expression of phase change latent heat in the convection heat transfer equation in the DDF method, and introduced MRT to correct the nonlinear coefficient of the equation. However, this paper did not describe the analysis from REV to macroscale, and achieve no physical meaning overlap from micro to macro. Zhang et al. used the DDF-LBM to simulate the hydrate dissolution process at the pore scale and calculated the variation of macroscale permeability with hydrate saturation[12]. However, the literature did not conduct sufficient research on the effective thermal conductivity, hydrate activation reaction and phase change latent heat, resulting in inaccurate calculation results. Ji et al. used the MRT-LBM to study the two-phase curves of gas-water flow during hydrate dissolution, but did not consider the influence of multi-field coupling on the dissolution law[13]. Ross-Jones et al. used the DDF-LBM to study the simulation results of the effective thermal conductivity of porous structures at the nanoscale[14]. The calculation results verified that the DDF-LBM can be directly used to simulate convective heat transfer at the nanoscale, but the model did not consider the effects of forced convection and phase change on flow simulation. In summary, LBM can be applied to the microscale simulation of porous media, but the key parameters of the methods need to be modified to make the calculation results more satisfied with really physical phenomena such as convective heat transfer and phase change entropy absorption.

In this paper, the influence of hydrate activation reaction on flow-heat transfer will be considered, based on DDF-LBM. And a simulation calculation method of heat flow and solidification will be designed, based on fully multi-field coupling, by TFMC. According to the calculated results, the permeability, thermal conductivity, activation reaction rate and other parameters of the models under multi-field coupling conditions would be analyzed, which can provide more accurate reference for mining natural gas hydrates using geothermal fluids.

2. Methodology

Many literatures have shown that natural gas hydrate reservoirs currently exist in different modes at microscale[12]. In order to better characterize the hydrate dissolution mechanism, we designed two different microstructures, the Pore-filling type and the Grain-coating type, based on two common hydrate reservoir modes. The model size Length × Width = 250×110, where the size of hydrate and rock particles is 11×11. The square particles of rock and hydrate make the fluid flowing in the pores more similar to the flowing in slender pipes. This design ensures that the models have a porosity closer to the real reservoir (less than 15%) while allowing thermal fluid to flow smoothly[2].

Figure 1. Two different microstructures of typical hydrate reservoirs.

Figure 1. Two different microstructures of typical hydrate reservoirs.

According to the literatures[5], when simulating the thermal fluid flowing through porous media with hydrates, we need to set the Mass Equation (1), Momentum Equation (2), and Energy Equation (3) under the microscale as the controlling functions of the model. The controlling functions used in LBM are as follows.

$${\displaystyle \frac{\partial \rho }{\partial t}}+\nabla \cdot \left(\rho \mathit{u}\right)=0$$

(1)

$${\displaystyle \frac{\partial \rho \mathit{u}}{\partial t}}+\nabla \cdot \left(\rho \mathit{u}\mathit{u}\right)=-\nabla \cdot p+\nabla \cdot \left(\upsilon \nabla \mathit{u}\right)$$

(2)

$${\displaystyle \frac{\partial \rho c_{p}T}{\partial t}}+\nabla \cdot \left(\rho c_{p}T\mathit{u}\right)=\nabla \cdot \left(\lambda \nabla T\right)+S_T,S_T=\Delta H\cdot R_h$$

(3)

where $\rho$, $\mathit{u}$, $p$, $T$ stand for the density, velocity, pressure and temperature. $\upsilon$ is the viscosity of fluid. $c_p$ is the heat capability of fluid or solid. $\lambda$ is the heat conductivity. $S_T$ is the heat source term. $\Delta H$ is the heat absorption entropy of hydrate dissolution. $R_h$ is the hydrate dissolution rate.

Some literatures have shown the experiments results under microscale, with the actual pressure and temperature conditions of natural gas hydrate reservoirs[15,16]. The simulations have shown that, even if hydrates release gases, such as methane after dissolution, these types of gas will be difficult to gather as continuous bubbles in water due to the high pressure in storage and will basically be mixed in the flowing water as isolated particles. Therefore, only two structural grids, fluid and solid, were established in the LBM of this paper. The influence of continuous bubbles on the flow were not considered, and only a little correction to the viscosity and density of the liquid, based on mixed bubbles.

2.1. Lattice Boltzmann Method for Convection Simulation

The simulation fields in DDF-LBM include density field and temperature field. The Distribution Equation (4) and Equilibrium Distribution Function (5) of the density field and the Distribution Equation (6) and Equilibrium Distribution Function (7) of the temperature field in the D2Q9 model are as follows.

$$f_a\left(x+\Delta x,t+\Delta t\right)-f_a\left(x,t\right)=-{\Omega }_f\left[f_a\left(x,t\right)-{f_a}^{eq}\left(x,t\right)\right]$$

(4)

$${f_a}^{eq}=\rho {\omega }_a\left[1+3e_a\cdot u+{\displaystyle \frac{9}{2}}{\left(e_a\cdot u\right)}^2-{\displaystyle \frac{3}{2}}{u}^2\right]$$

(5)

$$g_a\left(x+\Delta x,t+\Delta t\right)-g_a\left(x,t\right)=-{\Omega }_T\left[g_a\left(x,t\right)-{g_a}^{eq}\left(x,t\right)\right]-{\omega }_a{S}_T\Delta t$$

(6)

$${g_a}^{eq}=\rho {\omega }_a\left[1+3e_a\cdot u\right]$$

(7)

$$\rho ={\displaystyle \sum f_a},\rho u={\displaystyle \sum f_a{e}_a},P=\rho RT={\displaystyle \frac{1}{3}}\rho ,T={\displaystyle \sum g_a}$$

(8)

where ${\Omega }_f$ and ${\Omega }_T$ are the relaxion frequency of density and temperature. In solid heat conduction problems, ${\Omega }_T$ is generally a function of the thermal diffusivity $\alpha$. However, in conjugate convection, energy transfer is not only heat conduction, but also the flow affects the movement of internal energy, which is forced heat convection. Therefore, in the DDF-LBM simulation for conjugate convection, the thermal diffusion rate of fluid should be the combined effect of heat conduction and energy migration, which is, ${\Omega }_T={\Omega }_{\alpha }+{\Omega }_f$, where ${\Omega }_{\alpha }={\displaystyle \frac{1}{3\alpha +0.5}}$ and ${\Omega }_f={\displaystyle \frac{1}{3\upsilon +0.5}}$. $\alpha$ stands for the thermal diffusivity of fluid at stagnation[17]. Using the modified relaxation frequency ${\Omega }_T$ in the fluid temperature grid can more accurately calculate the coupling of fluid convection and solid heat transfer in porous media under microscale.

2.2. Lattice Boltzmann Method for Phase Change

Since the heat source term $S_T$ in the temperature field distribution equation only exists at the equilibrium temperature and pressure, it is zero under other conditions. Therefore, in the LBM simulation codes, the processing of source term $S_T$ can be carried out according to the following equations.

$$\left\{\begin{array}{cccccc}when& T<T_{eq}& \Rightarrow & S_T=0& \Rightarrow & f_H\to f_H,g_H\to g_H\\ when& T>T_{eq}& \Rightarrow & \Delta T_{\Delta t}=T\left(t+\Delta t\right)-T\left(t\right)& \Rightarrow & {\displaystyle \sum \Delta T_{\Delta t}}\\ when& {\displaystyle \sum \Delta T_{\Delta t}}<{\displaystyle \frac{\Delta H}{c_H{M}_H}}& \Rightarrow & T\left(t+\Delta t\right)=T_{eq}& \Rightarrow & f_H\to f_H,g_H\to g_H\\ when& {\displaystyle \sum \Delta T_{\Delta t}}>{\displaystyle \frac{\Delta H}{c_H{M}_H}}& \Rightarrow & f_H\to f_F,g_H\to g_F& \Rightarrow & S_T=0\end{array}\right.$$

(9)

where $T_{eq}$ is the equilibrium temperature. $f_H$, $g_H$ stand for the lattice of hydrate as solid. $f_F$, $g_F$ stand for the lattice of hydrate as fluid. $\Delta T_{\Delta t}$ is the temperature corresponding to the cumulative internal energy of the solid hydrate lattice reaching equilibrium. $M_H$ is the molecular weight of the hydrate. Common Type I and Type II natural gas hydrate $C{H}_4\cdot 5.5H_{2}O$ have a molecular weight of 115g/mol[18]. Then source term in Equation (9) can be calculated as follow.

$${\displaystyle \frac{\Delta H}{c_H{M}_H}}={\displaystyle \frac{52\mathrm{kJ}/\mathrm{mol}}{2.1\mathrm{kJ}/\mathrm{kg}\cdot \mathrm{K}}}\cdot {\displaystyle \frac{1}{0.115\mathrm{kg}/\mathrm{mol}}}\approx 215\mathrm{K}$$

(10)

It means that, hydrate in equilibrium needs to absorb the internal energy corresponding to 215K before it can dissolve from the solid phase into fluid phase.

2.3. Lattice Boltzmann Method for Hydrate Reaction

The macroscopic model of hydrate dissolution and activation can be characterized by the K-B model, as shown in the following[7]:

$$R_h={\displaystyle \frac{\Delta n_h}{\Delta t}}=k_d\cdot A_{hs}=A_{hs}{k}_{d0}\exp \left({\displaystyle \frac{-\Delta E}{RT}}\right)\left(f_{eq}-f\right)$$

(11)

where $n_h$ is the number of moles of dissolved hydrate, $A_{hs}$ is the reaction contact area, $\Delta E$ is the activation energy, $k_{d0}$ is the intrinsic activation rate, $f_{eq}$ is the equilibrium pressure. The K-B model takes into account the influence of temperature-pressure equilibrium curves during hydrate dissolution and the influence of reaction contact area.

According to the literatures, the K-B model is based on a statistical experiment of each phase of hydrate at the equilibrium state in a sealed container under approximately steady-state conditions, by reducing pressure or increasing temperature. And calculated the dissolution rate of hydrates within steady-state conditions. The methods did not fully consider the effect of the thermal fluid injected into the hydrate. Therefore, the reaction rate $k_{d0}$ in the literature is a fixed value after fitting, but in fact the reaction rate changes with the temperature gradient, which cannot be calculated from the experiments. Therefore, DDF-LBM can be used to intuitively count the contact area during dissolving, and according to the average pressure of the lattices during dissolving, the change of reaction rate $k_d$ as the hydrate saturation decreases can be given, which means that, this change has relationship with the hydrate saturation in this model. The dissolved and dissolving lattice can be illustrated in the Figure 2, including the reacting boundaries.

The undissolved grids shown in Figure 2 are the equations shown in the first-row of Eq. (9), that is, the current grids temperature have not reached the equilibrium condition. The dissolving grids shown in Figure 2 are those in the second and third rows of Eq. (9), that is, the current grids satisfied the equilibrium condition, but the accumulated energy does not satisfy the phase transition condition. These grids maintain solid properties, and their contact boundaries with the flow grids in the pores are the activation contact areas $A_{hs}$, as shown by the red lines in the Figure 2. The dissolved grids shown in Figure 2 are those in the fourth-row of Eq. (9), that is, the current accumulated energy of these grids was over the limit of heat absorption entropy. Then those grids will alter from solid into fluid, and being parts of collision calculation of fluids. Then the fluid-solid boundary evolutions with the hydrate activation reaction are realized.

2.4. Lattice Boltzmann Method for Unstable Simulation

The simulation of hydrate dissolution in porous media is a process with approximately steady-state flowing and non-steady-state heat transferring, so the statistical method of the model calculation results needs to be revised[19]. The average permeability $K_a$ of the porous structure can be calculated using the differential form of Darcy's law, as shown in Eq. (12). The efficient thermal conductivity ${\lambda }_a$ can be calculated using the differential form of the energy equation, as shown in Eq. (13). The model in this paper can be verified using experimental results in the literatures, and the Eq. (14) (K-B model) can be used to solve the reaction rate $k_d$.

$$K_a={\displaystyle \frac{{\mu }_w{L}_x\cdot \left(\overline{u_{in}}+\overline{u_{out}}\right)}{2\left(\overline{P_{in}}-\overline{P_{out}}\right)}}$$

(12)

$${\lambda }_a=\left[{\displaystyle \frac{{\displaystyle \sum {\rho }_i{c}_{i}T}}{\Delta t}}+{\displaystyle \frac{{\rho }_w{c}_w\left(\overline{u_{in}}+\overline{u_{out}}\right)}{2}}\cdot {\displaystyle \frac{\left(\overline{T_{in}}-\overline{T_{out}}\right)}{L_x}}\right]/{\displaystyle \frac{\left(\overline{T_{in}}+\overline{T_{out}}-\overline{T_{mid}}\right)}{L_x^2}}$$

(13)

$$k_d={\displaystyle \frac{R_h}{A_{hs}}}=k_{d0}\exp \left({\displaystyle \frac{-\Delta E}{RT}}\right)\left(P_{eq}-\overline{P_r}\right)$$

(14)

In the simulation of dissolution process, the pressure, velocity and temperature fields at the entire physical microscale are non-uniform. Therefore, it is impossible to directly using import and export statistics results to calculate. It is necessary to divide the inhomogeneous physical fields into regions so that the calculation results can be more approximately stable. In this study, the simulation area will be divided into eight statistical areas of equal size. And use the above equations to calculate the permeability and efficient thermal conductivity in each area, then average them. This can improve the impact of unsteady state in the physical field within the models, making the calculation results more accurate.

3. Results and Discussions

Set the left side as the hot fluid inlet and the right side as the hot fluid outlet. Other parameters and boundaries can be set as Tab. 1, including initial conditions and different lattices[20,21]. Set the hydrate totally dissolved as the ending point. Simulate the pressure, velocity and temperature at different times in the porous model. Fully consider the impact of fluid-solid boundary behaviors after hydrate dissolution. The calculation results of the structural field, velocity field, and temperature field at some time nodes are as follows.

Table 1. Parameters used in the DDF-LBM model with TFMC.

Parameters	Symbols	Value with Units
Kinetic viscosity	$\upsilon$	1×10-6 m2/s
Heat diffusivity of liquid	$\alpha$	1.31×10-7 m2/s
Thermal conductivity of rock	${\lambda }_s$	2.0 W/m/K
Thermal conductivity of hydrate	${\lambda }_H$	9.0 W/m/K
Thermal conductivity of liquid	${\lambda }_w$	0.55 W/m/K
Density of hydrate	${\rho }_h$	920 kg/m3
Density of rock	${\rho }_s$	2650 kg/m3
Specific heat of rock	$c_s$	1.0 kJ/kg/K
Specific heat of hydrate	$c_h$	2.08 kJ/kg/K
Inlet temperature	$T_{in}$	293 K
Phase change temperature	$T_{eq}$	274 K
Original temperature	$T_w$	271 K
Inlet velocity	$u_{in}$	0.1 m/s
Original pressure	$P_w$	2.3 MPa
Length scale	$L_r$	2×10-7 m
Time scale	$t_r$	2×10-8 s
Height	$H$	10×10-7 m
Activation energy	$\Delta E$	81 kJ/mol
Reaction heat absorption	$\Delta H$	52 kJ/mol
Phase equilibrium curve		$P_{eq}=\exp \left(38.98-8533.8/T_{eq}\right)$

Table 1. Parameters used in the DDF-LBM model with TFMC.

Parameters	Symbols	Value with Units
Kinetic viscosity	$\upsilon$	1×10-6 m2/s
Heat diffusivity of liquid	$\alpha$	1.31×10-7 m2/s
Thermal conductivity of rock	${\lambda }_s$	2.0 W/m/K
Thermal conductivity of hydrate	${\lambda }_H$	9.0 W/m/K
Thermal conductivity of liquid	${\lambda }_w$	0.55 W/m/K
Density of hydrate	${\rho }_h$	920 kg/m3
Density of rock	${\rho }_s$	2650 kg/m3
Specific heat of rock	$c_s$	1.0 kJ/kg/K
Specific heat of hydrate	$c_h$	2.08 kJ/kg/K
Inlet temperature	$T_{in}$	293 K
Phase change temperature	$T_{eq}$	274 K
Original temperature	$T_w$	271 K
Inlet velocity	$u_{in}$	0.1 m/s
Original pressure	$P_w$	2.3 MPa
Length scale	$L_r$	2×10-7 m
Time scale	$t_r$	2×10-8 s
Height	$H$	10×10-7 m
Activation energy	$\Delta E$	81 kJ/mol
Reaction heat absorption	$\Delta H$	52 kJ/mol
Phase equilibrium curve		$P_{eq}=\exp \left(38.98-8533.8/T_{eq}\right)$

Figure 3. Cloud maps of structural field, velocity field, and temperature field at 0.24mS and 1.80mS.

Figure 3. Cloud maps of structural field, velocity field, and temperature field at 0.24mS and 1.80mS.

3.1. Simulation Verifying by Reaction Rate

Changing the hydrate storage state may not change the dissolving rate much. Figure 4 shows the variation of residual hydrate saturation over time for two different pore structures. $S_h$ is the absolute saturation of hydrate in the pores and $S_h=V_h/V_P$. $V_h$ and $V_P$ are the solid hydrate volume and pore volume respectively.

Although the initial hydrate saturations of the grain-coating and the pore-filling types are different, their dissolution rates change in the same way because of the same boundary conditions of inlet velocity and temperature. By linear fitting the curves in Figure 4, the hydrate dissolution rates of the two models can be approximately described by the following equations.

$$\left\{\begin{array}{l}{n}_h=0.98\times {10}^{-3}\cdot t+5.9\times {10}^{-7}\left(\mathrm{Grain}-\mathrm{coating}\right),R^2=0.97\\ n_h=0.99\times {10}^{-3}\cdot t+5.8\times {10}^{-7}\left(\mathrm{Pore}-\mathrm{filling}\right),R^2=0.98\end{array}\right.$$

(15)

According to Eq. (15), the fitting results of the hydrate dissolution rate $R_h$ of the two models can be directly obtained as below.

$$R_h={\displaystyle \frac{\Delta n_h}{\Delta t}}=\left\{\begin{array}{l}0.98\times {10}^{-3}\mathsf{\mu }\mathrm{mol}/\mathrm{s}=3.52\mathsf{\mu }\mathrm{mol}/\mathrm{h}\left(\mathrm{Grain}-\mathrm{coating}\right)\\ 0.99\times {10}^{-3}\mathsf{\mu }\mathrm{mol}/\mathrm{s}=3.55\mathsf{\mu }\mathrm{mol}/\mathrm{h}\left(\mathrm{Pore}-\mathrm{filling}\right)\end{array}\right.$$

(16)

Substituting the hydrate dissolution rate $R_h$ into Eq. (14), we can obtain the hydrate activation reaction rate per contact area simulated by LBM, as known as $k_d$. The activation reaction area $A_{hs}$ can be obtained by counting the solid-liquid phase field boundary in Figure 2. The temperature difference $\Delta T$ in Eq. (17) comes from the temperature difference between the inlet $T_{in}$ and the boundary $T_w$ in Tab. 1.

$$\begin{array}{l}{k}_d={\displaystyle \frac{R_h}{A_{hs}}}=\left\{\begin{array}{l}4.06\times {10}^4{\mathrm{mol}/\mathrm{h}/\mathrm{m}}^2\left(\mathrm{Grain}-\mathrm{coating}\right)\\ 4.87\times {10}^4{\mathrm{mol}/\mathrm{h}/\mathrm{m}}^2\left(\mathrm{Pore}-\mathrm{filling}\right)\end{array}\right.\\ k_d=6.464\times {\left(\Delta T\right)}^{2.05}=3.65\times {10}^3{\mathrm{mol}/\mathrm{h}/\mathrm{m}}^2,Kamath et al. \left(1987\right)\end{array}$$

(17)

The simulation results in Eq. (17) are familiar with the experimental results in the literatures[22]. However, the experimental condition in the literature is to directly heat the upper surface of the hydrate in a sealed container, the thermal fluid is not injected into the porous structure with hydrate. Therefore, the hydrate dissolution reaction rate caused by the injection of hot fluids simulated by LBM must be greater than the experimental results, but not much more. It is further proved that the DDF-LBM simulation method that comprehensively considers TFMC in this paper can approximately and accurately describe the hydrate dissolution process when thermal fluid is injected into porous media.

3.2. Relationship between Saturation and Permeability

Firstly, the solid boundary state in the model under different saturation conditions will be retained, as shown in the structure map in Figure 3. Then, the velocity and pressure field distributions under different solid boundary conditions were calculated using the BLGK method with 8 partitions method. The approximate stable average permeability $K_a$ is calculated using Eq. (12). The average permeability of the grain-coating model increased from 14mD at undissolved state to 116mD at totally dissolved state. The average permeability of the pore-filling model increased from 14mD at undissolved state to 107mD at totally dissolved state. Then the average permeability was processed as dimensionless relative permeability $K_r$, and it was found that the hydrate saturation $S_h$ and relative permeability $K_r$ basically conformed to the power exponential relationship. According to relevant literatures, there are two typical models: Masuda model (18) and KC model (19), which can describe the power exponential relationship[23,24]. The fitting results based on the calculation results are shown in the following figures and equations.

$$K_r={\displaystyle \frac{K}{K_a}}={\left(1-S_h\right)}^n,Masuda et al.\left(1997\right)\left\{\begin{array}{l}n=2.92,R^2=0.99\left(\mathrm{Grain}-\mathrm{coating}\right)\\ n=3.09,R^2=0.99\left(\mathrm{Pore}-\mathrm{filling}\right)\end{array}\right.$$

(18)

$$K_r={\displaystyle \frac{K}{K_a}}={\displaystyle \frac{{\left(1-S_h\right)}^3}{{\left(1-B\times S_h\right)}^2}},Kozeny-Carman\left(2014\right)\left\{\begin{array}{l}B=0.052,R^2=0.99\left(\mathrm{Grain}-\mathrm{coating}\right)\\ B=-0.049,R^2=0.99\left(\mathrm{Pore}-\mathrm{filling}\right)\end{array}\right.$$

(19)

Figure 5. Fitting results between hydrate saturation and relative permeability by Masuda model and KC model.

Figure 5. Fitting results between hydrate saturation and relative permeability by Masuda model and KC model.

From the fitting results above, the average permeability calculation results of different saturations can be obtained by using the steady-state LBM. The fitting results of the Masuda model and KC model are basically the same, which is consistent with the conclusions in the literatures. This method can be used to accurately given the expression of permeability by saturation and provide accurate results for other coupled parameters.

3.3. Relationship between Saturation and Heat Conductivity

The effective thermal conductivity can only be calculated using a non-steady-state method, because Eq. (3) is a partial differential equation. Eq. (13) can be used to calculate the effective heat conductivity at different hydrate saturation. According to the literatures[25,26], if the medium of the porous is single or has a simple distribution, the average thermal conductivity can be directly calculated using the multi-component percentage weight of hydrate, rock, and fluid. However, the porous structure of hydrate reservoirs is relatively complex. So, the non-steady-state calculation results and component average calculation results of different structures are compared in the figure below.

Figure 6. Calculation results between hydrate saturation and effective heat conductivity.

Figure 6. Calculation results between hydrate saturation and effective heat conductivity.

The simulation results were shown in Figure 5, although the effective thermal conductivity fluctuates as the saturation decreases, the overall result gradually increases as the saturation decreases. The effective thermal conductivity of the grain-coating hydrate is smaller than that of the interval hydrate. Fitting the above calculation results, it is found that they can be approximately fitted by a power exponential equation, and the results are shown as follow.

$${\lambda }_a=a\cdot {S_h}^b\left\{\begin{array}{l}a=4.49,b=-0.062\left(\mathrm{Grain}-\mathrm{coating}\right),R^2=0.42\\ a=5.40,b=-0.082\left(\mathrm{Pore}-\mathrm{filling}\right),R^2=0.78\end{array}\right.$$

(20)

$${\lambda }_a={\varphi }_h\cdot {\lambda }_h+{\varphi }_s\cdot {\lambda }_s+{\varphi }_f\cdot {\lambda }_f\left(\mathrm{Heat} \mathrm{Flux} \mathrm{Parallel}\right)$$

$${\lambda }_a=\left({\varphi }_h\cdot {\lambda }_h+{\varphi }_s\cdot {\lambda }_s\right){\displaystyle \frac{2\left(1-{\varphi }_f\right)\left({\varphi }_h\cdot {\lambda }_h+{\varphi }_s\cdot {\lambda }_s\right)+\left(1+2{\varphi }_f\right){\lambda }_f}{\left(2+{\varphi }_f\right)\left({\varphi }_h\cdot {\lambda }_h+{\varphi }_s\cdot {\lambda }_s\right)+\left(1-{\varphi }_f\right){\lambda }_f}} \left(\mathrm{Sphere} \mathrm{Accumulation}\right)$$

Compared with the results of heat flux parallel and sphere accumulation, the rules are just the opposite. The reason is that, after the hydrates are reduced, the permeability increases, resulting in an increase in flow rate. And then the heat convection rate increases, resulting in an increase in the effective thermal conductivity.

3.4. Relationship between Saturation and Reaction Rate

Calculate the average pressure $\overline{P_r}$ at the activation reaction area, what is shown in Figure 2, then we can get the intrinsic activation reaction rate at different dissolution time or hydrate saturation by Eq. (14). The simulation results are shown as follow.

It can be seen from Figure 7 that the intrinsic reaction rate gradually decreases with the decrease of saturation. Comparing the simulation results of two different structures, the intrinsic reaction rate of the grain-coating type hydrate is smaller than that of the pore-filling type. By fitting the results, the expressions of the intrinsic reaction rate can be obtained.

$$k_{d0}=c\cdot e^{d\cdot S_h}\left\{\begin{array}{l}c=3.31\times {10}^4,d=2.65,R^2=0.83\left(\mathrm{Grain}-\mathrm{coating}\right)\\ c=3.99\times {10}^4,d=3.08,R^2=0.94\left(\mathrm{Pore}-\mathrm{filling}\right)\end{array}\right.$$

(21)

$$\begin{array}{l}{k}_{d0}=3.6\times {10}^4{\mathrm{mol}/\mathrm{m}}^2/\mathrm{Pa}/\mathrm{s},Clarke et al. \left(2001\right)\\ k_{d0}=1.24\times {10}^5{\mathrm{mol}/\mathrm{m}}^2/\mathrm{Pa}/\mathrm{s},Kim et al. \left(1987\right)\end{array}$$

Kim et al. used experimental methods to obtain the intrinsic activation reaction rate $k_{d0}$. However, Clarke et al. used a combination of experimental and numerical simulation methods to prove that Kim's results may be too large[27,28]. The fitting results in Eq. (21) show that if thermal fluid injection dissolution was considered, $k_{d0}$ will be larger and related to the saturation $S_h$.

Compared with the fitting results obtained above, the K-B model results are too large because the actual reaction contact area under micro-scale conditions is not fully considered. Not all contact boundaries between solid hydrates and fluids, only the contact boundaries between hydrates and fluids in the activation reaction are the actual reaction contact areas $A_{hs}$.Clarke's simulation results obtained a fixed value because it did not fully consider the influence of temperature gradient under micro-scale conditions. Although the inlet temperature $T_{in}$ remains unchanged during the dissolution process, the front of the activation reaction gradually advances as the hydrate dissolves. By comparing the structural field and temperature field results at different times in Figure 3, it is shown that the distance between the dissolving lattices (green) and the entrance gradually increases as the saturation decreases, so the temperature gradient of the model decreases as the saturation decreases. Therefore, the intrinsic reaction rate of hydrates is actually correlated with the hydrate saturation, and the K-B model can be further modified.

4. Conclusions

In this paper, the DDF-LBM was used to simulate the dissolution of natural gas hydrate caused by the thermal fluid injection under micro-scale porous model. The simulation results were solved by unsteady difference method and it was found that the dissolution rate of thermal fluid injection was higher than that of surface heat absorption only. According to the simulated velocity field, pressure field, temperature field and structure field results, the key physical parameters under different hydrate saturation conditions were calculated, and their exponential fitting equations were given. The fitting results were compared with those from other literatures, which further proved that the DDF-LBM was feasible and realized the TFMC fully multi-field coupled dynamic simulation of thermal fluid injection to dissolve natural gas hydrates. Based on the above results, the following conclusions can be drawn.

(1) When DDF-LBM simulates conjugate convection, the relaxation frequency of the thermal lattice in the pores needs to be modified to consider the simultaneous existence of thermal diffusion and thermal convection. The modification is necessary when simulating forced convection heat transfer, otherwise the calculated results will be much smaller. When describing the changes in fluid-solid boundaries, we must first determine whether the current solid hydrate lattices have reached the P-T phase equilibrium condition, and then determine whether the cumulative internal energy absorbed by the current lattice exceeds the heat required for its reaction. The lattices that have reached phase equilibrium conditions but do not satisfy the cumulative endothermic conditions are dissolving solid hydrate lattices, and their contact boundaries with the flowing lattices are the reacting contact areas in the K-B model. The reaction contact area can be directly obtained from the structural field simulation results, so that the effect of dissolution endothermicity on internal energy transfer in the convective heat transfer process can be described at micro-scale.

(2) By comparing the results of heating dissolution experiments in closed containers in other literature, the hydrate dissolution rate under the thermal fluid injection conditions simulated in this paper is slightly higher than the result of only heat absorption on surface. This is because the hydrate reaction rate calculated from the experiments in the literature was based on surface heat diffusion only, and did not take into account the actual situation of thermal fluid flowing under the porous structure of the reservoir. Therefore, the activation reaction rate calculated by DDF-LBM is more consistent with the real reservoir situation. Compared the dissolution rates of the grain-coating type and the pore-filling type, although the hydrate structures and initial saturations of the two models are different, their dissolution rates are basically the same due to the same boundary conditions of injection flow velocity and temperature. It is shown that the DDF-LBM simulation with TFMC multi-field coupling in this paper can approximately and accurately describe the hydrate dissolution process when thermal fluid injected into porous media.

(3) The modified expressions of the key physical parameters of TFMC multi-field coupling can be obtained, by fitting the hydrate saturation with the average permeability, effective thermal conductivity and intrinsic reaction rate. The fitting results show that the parameters’ trending of two different storage structures are the same. The average permeability increases as the saturation decreases, and the relative permeability has an exponential relationship with saturation. The effective thermal conductivity increases slightly with decreasing saturation and increases substantially when most of the hydrates are dissolved. Therefore, there is an approximate linear or exponential relationship between the effective thermal conductivity and the saturation, and a power exponential equation was used for fitting in this paper. The intrinsic reaction rate decreases with decreasing saturation, which is different from the fixed value given by the K-B model in the literature. This is because the model does not fully consider the impact of temperature gradient and reaction contact area on the results at micro-scale, so the K-B model can be further modified. The above modified expressions of key parameters can provide necessary references for simulating the real-time full coupling of thermal fluid injection into hydrates under macros-scale conditions.