NON-ISOTHERMAL KINETIC MODEL OF THE METHANE HYDRATE DISSOCIATION PROCESS AT TEMPERATURES BELOW ICE MELTING POINT Текст научной статьи по специальности «Физика»

Non-Isothermal Kinetic Model of the Methane Hydrate Dissociation Process at Temperatures Below Ice Melting Point

I.G. Donskoy1*, S.Ya. Misyura2

1Melentiev Energy Systems Institute of Siberian Branch of Russian Academy of Sciences, Irkutsk, Russia 2Kutateladze Institute of Thermophysics of Siberian Branch of Russian Academy of Sciences, Novosibirsk, Russia

Abstract — The paper proposes a new model of gas hydrate particles dissociation at high heat fluxes. In this model, the process of hydrate decomposition occurs under the conditions of competition between heterogeneous kinetics and gas filtration. An analysis of the experimental data gives new values of the kinetic coefficients for the hydrate dissociation at low temperatures. The calculation results make it possible to reproduce experimental data on the dynamics of methane hydrate powder dissociation (including dissociation under the conditions of gas burning above the surface) and to describe the phenomenon of self-conservation in terms of changes in the pore structure of the ice crust. The submodel of dissociation of a single particle is embedded in the mathematical model of transport processes in the powder layer, which allows analyzing the heterogeneity of heating and the collective effects of dissociation.

Index Terms — natural gas, gas hydrates, dissociation, kinetics, heat and mass transfer, numerical modeling.

I. Introduction Gas hydrates are crystalline compounds formed by the interaction of gases with water at low temperatures and high pressures [1]. The wide interest in gas hydrates is due to their formation in natural and technical systems. For example, gas hydrate formation can be a problem during gas production in the cold season and its transport over long distances [2-4]. Estimates of the reserves of combustible gases bound in the hydrated form are very large, but the specificity of gas hydrates (dispersion, inaccessibility) does not allow them to be extracted.

* Corresponding author. E-mail: donskoy.chem@mail.ru

http://dx.doi.org/10.38028/esr.2020.01.0003 Received February 28, 2020. Revised March 19, 2020. Accepted May 30, 2020. Available online July 31, 2020.

This is an open access article under a Creative Commons Attribution-NonCommercial 4.0 International License.

© 2020 ESI SB RAS and authors. All rights reserved.

According to the energy forecasts in [5], [6], cost-effective technology for the extraction of natural gas from hydrate deposits will not emerge until 2040. In a recent review of the prospects for the energy use of gas hydrates [7], the authors estimate the start of commercial production no earlier than 2030. This is due to some technical problems, for example, low gas flow rates. Despite the widespread occurrence of hydrates on the seabed and in permafrost (estimates of reserves vary, but even pessimistic methods give multiple excesses of the "traditional" natural gas reserves), their volume fraction in bearing rocks is usually too small to provide a stable high gas flow rate [8]. An important issue is also the safety of production since during gas extraction hydrate-bearing formation changes its mechanical characteristics, which means that it can lose stability. Nevertheless, gas hydrate programs are intensively developed abroad. Note that the USSR and Russia have successful long-term experience in the extraction of methane from gas hydrate deposits in northern conditions [9]. At present, the average cost of methane gas hydrate exceeds the cost of traditional natural gas by 2-3 times. However, research and industrial projects in this area are quite intensive [10-12]. In addition to underwater drilling, the possibility of extracting hydrate from permafrost zones is being investigated.

Another important aspect of the energy applications of gas hydrates is storage and transport. The bulk density of the gas in hydrated form is higher than that for compressed gas, and storage conditions are softer compared to liquefied gas [13]. Under certain conditions, gas hydrates can exhibit abnormal stability, which is called self-preservation: at atmospheric pressure and slight cooling, the hydrate particles get covered with a crust of ice, which prevents the outflow of gas from the inner layers [14]. Maintaining hydrates under self-preservation conditions requires much less energy, and fire and explosion hazard is significantly reduced. Gas hydrate transport can be cost-effective at short distances and small volumes of supply [15, 16], for example, for supply to small power systems. Other hydrate applications may include thermal energy storage, desalination, gas treatment [17], etc. The stability of gas hydrates is an important issue, in particular, in the studies

of the global climate. Methane is a greenhouse gas, and its natural emission includes, among other things, the decomposition of gas hydrates at seabed and permafrost [18]. To reliably predict the emissions, it is necessary to study the physicochemical features of hydrate dissociation under natural conditions. One way to use gas hydrates in the energy sector is to store CO2 in gas hydrate deposits. Under a suitable combination of conditions, methane can be displaced from its hydrate form by carbon dioxide [19, 20]. Such technologies will reduce CO2 emissions and control CH4 emissions which, according to some geological models, can be dangerous (i.e. clathrate gun hypothesis [21]).

The existing problems of using hydrates are somehow related to an insufficiently complete understanding of the physics and chemistry of the processes that occur during the formation and dissociation of the hydrate phase [22, 23]. This, in turn, is associated with a variety of phenomena taking place during hydrate processing. Usually, the theoretical consideration of such processes employs simplified approaches. They simulate the propagation of the dissociation front with various moving boundary conditions [24, 25], or use diffusion models [26]. More detailed models take into account the formal kinetics of physicochemical transformations [27], heat and mass transfer in particles [28] and rocks containing hydrates [29].

One of the features of the dissociation of gas hydrates is a wide range of scales: physical and chemical processes occur at the level of hydrate cells (10-20 angstroms), clusters (of the order of nm), pores and cracks (0.1-10 microns), micrograms (0.1-1 mm), granules (1-10 mm), agglomerates (1-10 cm), and, finally, hydrate deposits (1100 m).

Studies of the hydrate decomposition kinetics are possible only when sufficiently small (microns) particles are used [30]. There are also known studies, in which the molecular dynamics of the processes of decomposition and formation of gas hydrates is simulated [31-33]. In reality, however, hydrate particles are much larger, therefore, in the study of their dissociation, transport processes must be taken into account. In this case, the nature of the diffusion process remains uncertain, since dislocation diffusion in the ice layer allows one to explain the dissociation rate only in deep self-preservation. At the same time, it is known that the ice formed in the process of hydrate dissociation differs in properties from the ice obtained by freezing pure water. As a rule, it has a developed microporous structure, which results from the difference in the molar volumes of ice and hydrate [34]. Therefore, it is natural to assume that the key process responsible for the release of methane from the particle is filtration.

Previously, the filtration process was often considered as limiting on the scale of the gas flow in a porous rock containing gas hydrate. In this case, the decomposition of the hydrate occurs when the corresponding thermobaric

conditions are reached and is limited by processes on the scale substantially larger than the particle size [35-37], or kinetics of dissociation is given that does not take into account the structure of the particle [38].

Mathematical models of porous reactive particles are widely used in problems of chemical technology and combustion theory [39-43]. Chemical transformations in porous particles take place in catalytic processes and during pyrolysis and combustion of solid fuels. Despite the formal differences, there is a certain analogy between these processes and the decomposition of hydrates, which can be used in the theoretical analysis.

In the present research, a hydrate particle dissociation model is developed, taking into account the heat transfer of the particle with the environment and the filtration of the gas released during dissociation through the particle pores. The process of hydrate conversion into ice and gas is considered as a one-stage irreversible chemical reaction with the kinetic coefficients that are evaluated from an analysis of various experimental data.

II. Filtration-Kinetic Model of Dissociation of

Thermally Uniform Gas Hydrate Particle We consider a hydrate particle placed in a heating gas medium. Under suitable conditions (increasing temperature, lowering pressure), the gas hydrate reaches the stability boundary and releases gas (Fig. 1). Depending on the external conditions, the limiting factor may be heat transfer or mass transfer.

The mathematical model is built taking into account several simplifications and assumptions, namely:

• there are no impurities in the hydrate, and the dissociation products are pure gas and porous ice;

• the particle is considered spherically symmetric;

• gas hydrate dissociation is a one-stage irreversible reaction (association is excluded);

• dissociation front gradually moves from the surface towards the center of a particle;

• the temperature field in the particle is uniform;

• the gas pressure field in the particle pores is quasistationary;

• the impact of the Joule-Thomson effect is neglected. At the initial moment, the mass of the particle is

equal to the mass of the hydrate phase. The mass fraction of methane in the hydrate is denoted by B (for example, for methane, one can take B = 0.12. Then we can make the material balance equations for the components of the particle:

H

(1)

dmH r

dt - B

dm' 1 - B

dt — 1 B

dmG

= r - j J i

dt

(2) (3)

r

s

Here m is the mass, kg; indices H, I, G refer to hydrate, ice, and gas, respectively; r is the rate of gas formation during decomposition of the hydrate, kg/s; jg is the loss of gas mass due to mass transfer with the environment, kg/s. The thermal conductivity of a hydrate particle is of the order of 0.5 W/m/K [44]. Therefore, at moderate heating rates, it can be assumed with good accuracy that the Biot number for particles of the order of 1 mm is quite small. To be accurate, it is also necessary to take into account the contribution of the endothermic process, however, assessment of these disturbances can be made later. For small Biot numbers, the particle can be considered as thermally homogeneous, i.e. one can assume that the temperature gradient in it is not significant. In this case, the problem can be simplified without considering the process of thermal conductivity, and we can write the equation of heat and mass balance in the particle:

dm

= -g (4)

dt

dT

Cp = cHp m + cpm' + cGm

(6)

The density and heat capacity of ice and hydrate are taken from published data [45, 46]. The mass and linear dimensions of the particle are interconnected by obvious geometric relationships.

The fraction of unreacted hydrate can be found from the gas mass balance:

rdt

X = 1 --

Bmn

(7)

From equation (7), we can immediately obtain the expressions for the mass of the solid components of the particle: mH = m0X; m' = m0(1 - 5)(1 - X). In differential form, equation (7) can be written as follows:

dX ___

dt Bm0

Equations (4) and (5) are supplemented by the expression for the rate of the hydrate decomposition and gas flow from the surface:

(7a)

CP — = aS(Tout - T) - Qr - cG (T - t0 )jg (5)

In equation (5), Cp is the heat capacity of the particle, cpG is averaged specific heat capacity of gas, J/kg/K; m is the particle mass, kg; T is the temperature of the particle, K; Tout is the temperature of the surrounding medium, K; a is the particle heat transfer coefficient, W/m2/K; S is the particle outer surface area, m2; Q is the thermal effect of the hydrate dissociation, J/kg.

The last term on the right-hand side of (5) is associated with a change in enthalpy due to a decrease in the particle mass (dH = d(hm) = mdh + hdm). The change in particle mass is determined by the gas flow rate through the surface. We assume that the specific heat capacity Cp is additive in terms of components:

r = k0e~RTSr (Peq - P)

jK = pAp -

d

(8)

(9)

In equations (8) and (9), k0 and Ea are the pre-exponent factor and activation energy of the hydrate dissociation reaction, kg/m2/Pa/s and J/mol, respectively; Sr is the area of the dissociation front in particle, m2; P, Peq, Pout is the current (intraporous), equilibrium and external gas pressure, Pa; pg is gas density in the pores of the particle, kg/m3; kD is the permeability coefficient of porous ice; m is the dynamic gas viscosity, Pas; 6 is the effective thickness of the ice crust, m.

The equilibrium pressure of methane above the hydrate is usually expressed by an exponential function of temperature:

lg (p-dp)=a - a

(10)

Fig. 1.The scheme of dissociation of a spherical gas hydrate.

Coefficients A1 and A2 can be determined from experimental data or theoretical estimation [47]. In this paper, the following values are accepted: A1 = 4.77047; A2 = 920.701; DP = 0.5. These values give characteristic features for the dissociation of methane hydrate: equilibrium temperature of -80 °C at 1 atm and equilibrium pressure of 25 atm at 0 °C.

To determine the gas density and pressure, it is necessary to know the pore volume in the particle. The hydrate porosity is usually very low, and due to the difference in densities during dissociation, ice is formed with a porosity of the order of 0.1. The volume of gas in the particle is very small, and the contribution of gas to the density of the heat capacity is almost invisible. With external heating, the dissociation of hydrate particles usually occurs layer by layer, and the decomposition front moves from the surface to the center. Given this, the surface of the dissociation front and the outer surface of the particle are connected by a relationship arising from geometric considerations: if we denote by X the fraction of unreacted hydrate (with respect

t

e

to the initial amount), then the shrinking core model gives the expression for spherical particle:

S = So 3

(11)

Similarly, for the thickness of the ice crust, we obtain the expression (at linear approximation, i.e. for thin crust):

d (' - -13 )

(12)

In equation (12), d0 is the particle diameter. Thus, during dissociation, the area of decaying hydrate decreases, while at the same time the filtration resistance increases due to the growth of the ice crust.

Equation (9) models the Darcy filtration equation for the transfer of gas forced by pressure difference. In this case, a simplified form is adopted: the filtration rate is proportional to the average pressure difference in the particle, and the denominator in the expression for the pressure gradient is the thickness of the ice crust. Equation (5) is applicable to small X, when the change in the curvature of the dissociation front can be neglected. In a wide range of conversion degrees, we can consider the quasi-stationary problem [48, 49], then we can write:

d = d

1

(13)

2 X"13 -1

For a numerical study of the solutions to the system of equations (1-5), it is necessary to know the values of thermophysical and kinetic parameters. The most difficult parameter to evaluate is the permeability coefficient of porous ice. This parameter depends primarily on the size distribution of pores and their density. Following [50], it can be assumed that each pore can be represented as a cylindrical channel with an equivalent diameter equal to the average diameter of the pore. Then, the permeability coefficient for ice can be estimated as the total drag coefficient of pores per unit surface area of a particle. Assuming the laminar gas flow in small channels, we can write:

k = nl

D 32

(14)

In equation (14), dp is the average pore diameter, and n is the surface porosity (total pore cross-sectional area per unit area of ice surface). The average pore diameter can be estimated from the studies where the hydrate surface was scanned using microscopy [14, 51] (hereafter, it is taken equal to 0.5 ^m). In addition, an assumption is introduced that the average pore size does not change with the radius of the particle, which, of course, requires additional verification, but is applicable as a first approximation. Surface porosity is equal to the product of the surface pore density and the cross-sectional area of a single channel:

ji ^ n = a — dp 4 p

(15)

In equation (15), c is the surface pore density, m-2. This parameter varies both under the influence of thermobaric conditions and under the action of mechanical stresses in the particle; its current value most likely depends on the particle heating dynamics. Evaluation of pore density is the most difficult task. At the same time, this parameter is most important for modeling processes such as self-preservation. In terms of the proposed mathematical model, this process is accompanied by a sharp drop in pore density. Thus, for closure, the model requires an additional equation for the pore density, which relates this quantity to the conditions of the dissociation process (temperature and pressure). Obtaining such an equation (even in the form of an empirical formula) is of scientific and practical interest. Since we do not have them at the moment, we have to use a simplified approach, namely, setting a piecewise constant pore density for different parts of the experimental curves. This approach is essentially adjustable, however, it allows obtaining several interesting results.

At the initial point X = 1, equation (9) is singular, therefore, for correct calculations, we assume that at the initial time moment X0 = 1-e, where e is an arbitrarily small number (for example, of the order of the relative weighing error in experiments).

Preliminary studies have shown that to apply the model to explain the observed effects, it is also necessary to take into account in detail the conditions for the heat transfer of particles with external heat sources. Typically, in experiments, particles are heated by natural or forced convection of the airflow or combustion products (in the case of gas hydrates combustion), or by the heat flux passing through the wall of the vessel filled with hydrate. In both cases, there are characteristic effects of heat transfer in two-phase systems with phase transitions, as well as heat transfer in granular beds. The heat flux to a hydrate particle cannot be generally estimated simply by the temperature difference between the particle and the external environment: both the heat transfer coefficient and the local temperature difference (especially for the underlying layers of particles) in an unobvious way depend on the conditions specified in the experiment.

III. Low-temperature kinetics of methane hydrate dissociation

An important problem in hydrate dissociation modeling is the selection of reliable kinetic coefficients. At present, the published studies on the experimental measurement of the hydrate dissociation rate coefficient are very limited. This is due to some problems arising from such measurements: the mass transfer restrictions, the need to maintain isothermal conditions, self-preservation, etc. The results obtained by Bishnoi et al. [29, 52], where kinetic coefficients were measured for the dissociation of particles at temperatures above 273 K (which helped to avoid self-preservation) are most widely used. For negative temperatures, these coefficients, however, give a

significant error, which may be due to differences between the mechanisms of "low temperature" (t < 0°C) and "high temperature" (t > 0°C) dissociation. The following is an attempt to clarify the values of kinetic coefficients for the region of negative temperatures and to generalize the known experimental data on this topic.

The present mathematical model of gas hydrate dissociation is based on the model of a shrinking-core and irreversible kinetics of dissociation in the Bishnoi form. It is supposed that the particle is isothermal. In this case, one can write:

dX dt

k>

(Peq - P) S0X

(16)

Bmn

Here X is the fraction of unreacted hydrate, the ratio of the gas contained in the hydrated form to its initial amount; kd is the preexponential factor of the kinetic coefficient, kg/Pa/m2/ s; Ea is the activation energy of the dissociation process, J/mol; P is the gas pressure in the pores of the particle; Peq is the equilibrium gas pressure outside the hydrate particle; S0 is the surface of the particle; B is the mass fraction of gas in the hydrate; m0 is the initial mass of the particle.

For a spherical particle, we can write:

dX dt

1 o he

(Peq - P)pd20 X

Bpc

pd3

(17)

6k CeRT

(Peq - P)

BPo dc

X3

For isothermal conditions, the entire expression written before X2'3 can be considered constant, i.e. the specific dissociation rate increases with a decrease in the particle size, which is consistent with empirical observations. Complete decomposition of the particle is achieved at a finite time t = 3/K, where

K=

--L

6kadeRT (Peq - P)

BPo do

The half-life will be equal to:

= 3 (tf - ^05)

T(05 = K

The average half-life rate is: 0.5K

^0.5 =■

(18)

(19)

(20)

3 ( tf - ^05 )

To compare the data obtained under different conditions, it is necessary to take into account the temperature dependence of the reaction rate both for the kinetic coefficient (related to activation energy) and for the equilibrium pressure. The superposition of the temperature dependences of the kinetic coefficient and the equilibrium pressure can lead to a wide range of possible apparent

temperature coefficients of the dissociation rate. Therefore, it is necessary to separate these factors with sufficient accuracy.

It is worth noting that the values of the kinetic coefficients are very sensitive to the choice of the equation for the equilibrium curve.

The kinetic coefficient can be determined in two ways: by the differential method (by measuring the rate of the process and comparing its values with kinetic models) and by the integral method (by measuring the characteristic conversion time).

Differential method

From equation (17) we have:

2

X3

K = —

DX/

(21)

"At

where the numerator is the value of the kinetic function for the selected value of X, and the denominator is the average rate for a given moment in time, calculated as the difference approximation. For a constant heating rate, we have:

K = -

DTX3 ßDX

(22)

where b is the constant rate of temperature change.

Integral method

Knowing the characteristic time, for example, the time of complete conversion or half-life (which can be determined more precisely), we can write:

* = 3

k =

VT -

(23)

(24)

After the coefficient K is found, it is possible to isolate the component that corresponds directly to dissociation, according to (18):

kd = kde RT = '

BPodo

rK.

(25)

6 (Peq - P)

The pre-exponent and activation energy can be determined by constructing a regular Arrhenius plot in the coordinates ln(kd) ~ 1/T.

The two main sets of measurements that will be considered here are the data from [53] and [54]. Authors of [53] measured the relative intensity of the methane-water bond signal in samples of hydrate powder of different fractions at a constant heating rate of 1 K/min. The low rate of temperature change makes it possible to consider the conditions as isothermal between two consecutive measurements. Therefore, in this case, we can use the differential method, i.e. equation (21). In [54], the authors measured the average half-life rate of methane hydrate samples. The authors report possible temperature inhomogeneities, however, to a first approximation, we consider their experiments to be isothermal. The integral

e

E

2

t

t

E

1000/T

Fig. 2. Arrhenius plot constructed while processing the kinetic data [53, 54](the thick line is linear regression, the thin line is the kinetics of hydrate dissociation above 273 K [29]).

Fig. 3. The dependence of the residualfunction 0 on the values of the Arrhenius coefficients for the selected data set.

method, namely, equation (24), is well suited for processing these data.

Some experimental points are not suitable for analysis since, given constant P, equation (17) does not take into account self-preservation. According to the well-known literature data, this process begins at temperatures of the order of 240 K. Therefore, to isolate the kinetic coefficients without the influence of transfer processes, we have to exclude the temperature range from 240 K and higher. Nevertheless, there are still enough experimental points: 10 points from experiments [54] and about 20 from [53]. The processed results are shown in Fig. 2.

The average mass content of methane in the hydrate is assumed to be 0.12. The density of the hydrate is 910 kg/ m3. The average particle size for the data [53] is taken equal to the arithmetic mean for the boundaries of the indicated ranges. In [54], the samples were prepared by pressing a hydrated powder (an average size of 200-300 ^m) into cylinders with an average porosity of 0.29. In this case, the hydrate particles most likely agglomerate, although they do not turn into monolithic samples (as indicated by porosity measurements). Therefore, the size of 10 granules of the initial powder was adopted as the characteristic size d0.

Figure 2 shows that, despite significant differences in the experimental methods, the obtained kinetic coefficients quite well fall into a single dependence. The effective parameters of the Arrhenius form of the kinetic coefficient are equal to: k0 = 3.25x10-5 kg/Pa/m2/s; Ea = 34 kJ/mol. These values differ significantly from the kinetic coefficients from [29, 51], which may be due to the difference in the mechanism of dissociation. Interestingly, Clarke-Bishnoi (high-temperature) kinetic coefficient becomes larger than the obtained (low-temperature) values at a temperature close to 240 K (self-preservation starts here).

Fig. 4. Comparison of experimental and calculated kinetic curves of dissociation of methane hydrate fractions at a heating rate of 1 K/min [53]: solid lines — calculation, markers — measurements.

Fig. 5. Comparison of experimental and calculated data on the kinetics of methane hydrate dissociation (data [54]): (a) kinetic dissociation curves; (b) average integral half-life rate (solid lines — calculation, markers — measurements).

The estimation of kinetic coefficients from the measurement results is a typical inverse problem, the solution to which is not always unambiguous and stable. In this case, the solution was found by the least-squares method (search for the minimum of the non-negative residual function

min

in0 = ^raj (YJxp -Yjdc)2

(26)

Since the direct problem is nonlinear (especially in the non-isothermal case), and the set of experimental data includes experiments performed under various conditions, it is necessary to carefully select the weight coefficients for the residues obtained by approximating the kinetic coefficients for different options (heating rate, particle size, etc.).

The residual function includes discrepancies that arise when comparing the experimental and calculated data for four experimental sets: integral data [54]; isothermal kinetic curves [54]; non-isothermal kinetic curves with a constant heating rate [53]; non-isothermal kinetic curves with a varying heating rate [55]. In total, the residual function <P includes about fifty terms. The results of optimizing the residual function in the space of coefficients k0 and Ea are shown in Fig. 3.

One can see a region of parameters, where the pre-exponent and activation energy, varying over a fairly wide range, give approximately the same residues (the so-called compensation effect is observed). Therefore, for certainty in the calculations, the activation energy is taken to be 34 kJ/mol, as the generalized Arrhenius plot gives. The preexponential coefficient under these conditions is 0.019 kg/Pa/m2/s. For these values, calculations were performed according to the data of [53, 54], the results of which are shown in Figs. 4 and 5. The calculated curves are in a very good agreement with the experiment for small fractions (up to about 200 mm). With a further increase in size under the experimental conditions, transport processes most likely begin to affect the kinetics of decomposition (for example, the resistance of the ice crust increases, which is one of the mechanisms of self-preservation [56]). Therefore, the decomposition of particles after reaching a degree of conversion of 50% is slower than expected from the kinetic model.

IV. Effect of particle sizes of gas hydrates on the rate of their dissociation

The granule size of a gas hydrate determines its outer surface through which heat and mass transfer processes occur, and the total amount of stored gas. From the simplest geometric considerations, it is possible to evaluate the qualitative dependence of the hydrate dissociation rate on the particle size, namely, small particles dissociate faster, dX/dt ~ 1/d0. Such a simple dependence, however, will be observed only under isothermal conditions. If a gas hydrate dissociates at changing temperature, the dependence can change significantly due to nonlinear effects.

Equations (1-5) were numerically solved using a combination of explicit and implicit methods for ODEs. The equation for the kinetics of dissociation (17) is solved explicitly:

X* = Xk -1 - tK (Tk-1 ) ( Xk -1 )2/3 (27)

Equation (27) gives the gas production rate rk, the position of the dissociation front, and the ice porosity n. Knowing these quantities, we can solve the difference equation of the gas mass balance (3):

k k-1 k k n S mk -mk,1 =trk-t-Ç-—

j d

The pore volume VG is geometric relationships: in the

m,

Pg

œ RTk-1mG _ p

k out

V

MVk

(28)

determined from the "shrinking core" model

for the sphere VG = '' n(l - X). The gas density pG is p

specified iteratively (in most cases, however, an explicit approximation is sufficient). After all the components of the material balance are determined, the heat balance equation (5) is solved:

T - t" = ta (ToUI -T) - tQrk -1 J'C (T-T) (29)

p p p

Equation (27) is an explicit Euler method formula, and equations (28) and (29) are written as an implicit form of Euler method. The numerical scheme was obtained by splitting the initial system into steps related to a chemical reaction, convective mass transfer, and heat transfer. Typically, such schemes have the first order of accuracy related to t. In the calculations presented below, the time step was 10-4 s: such accuracy is sufficient for research purposes.

When considering the process of dissociation of a thin layer of gas hydrate, we can take the approximation of the reference particle [57, 58] and simulate the dynamics of the layer by the dynamics of dissociation of an individual particle. To test our model, experiments on the dissociation of a thin hydrate particles layer were carried out (the experimental technique is described in [55, 59]). Synthetic gas hydrate powder (mass fraction of methane 12%), poured into a metal tank with a layer with a height of 0.6 to 30 mm, was heated at atmospheric pressure and ambient temperature. In some experiments, the produced methane was ignited. Separate fractions (0.6 and 6 mm) and their mixtures were used.

The heat transfer conditions in the particle layer are rather difficult to establish with the necessary accuracy. In the experiment, the surface of the powder was controlled while the powder was heated by heat conduction through the walls of the tank, by the gas flow and the liquid water drainage, as well as flame radiation (in the experiments with hydrate combustion). Therefore, the heating of the layer was modeled by a constant heat flux. It was determined by the constant temperature difference between the particle surface and the gas at the heat transfer coefficient (q = aDT), which was estimated according to the recommendations for granular materials from [50].

Fig. 6. Gravimetric curves of methane hydrate powder. Solid lines: calculation at different values of dT; points: experiment (triangles- with gas burning, circles- without gas burning).

0.0006 0.0005 0.0004 i? 0.0003 0.0002 0.0001 0

1 '.

■0.5 1 -2 -3

0

2000 4000

Time, s

6000

8000

Fig. 7. The gas evolution rate during the dissociation of methane hydrate particles (numbers are average particle diameters, mm).

By varying the DT value in the calculations, we can find a suitable value for the experimental conditions. Under such conditions, the heat flux turns out to be a fitting parameter. The results of such a selection are shown in Fig. 6.

Hydrate dissociation without combustion of gas above particles is expected to occur more slowly due to a smaller heat flux. During dissociation with gas combustion above a layer, the temperature difference estimated using a mathematical model is about 100 K. It can be assumed that heat transfer occurs mainly between water vapor and melting ice on the surface of the particles. A high heat flux leads to intense heating of the powder; therefore, the hydrate passes rather quickly through the temperature region of self-preservation (240-270 K). The dissociation process without combustion is much slower. Therefore, in the experiment (circles in Fig. 6) a stepwise change in mass is observed: the hydrate particles remain in the self-preservation region for a time long enough to form stable ice crust on the surface. A model of the self-preservation phenomenon will be described in the next section.

The dynamics of methane hydrate dissociation with different particle sizes (reduced to the same mass) is shown in Fig. 7. Firstly, for the particles, the heating time to the dissociation temperature differs significantly: for small particles, heating takes time on the order of minutes, and for large particles, it takes tens of minutes. The peak dissociation rate also differs by an order of magnitude. Therefore, the dissociation and combustion of powders with different sizes of gas hydrate granules will occur in stages. Small particles are most reactive, however, their share in the powder depends on the method of hydrate production and processing (grinding, pressing, etc.). The production, storage, and transport of small particles are associated with additional energy costs and losses but the dissociation of large particles is much slower. Therefore, the question of the optimal fractional composition for different methods of gas hydrate processing is of interest. If the released methane ignites above the layer, then small particles could act as the initiator of dissociation.

v. Modeling the phenomenon of self-preservation of gas hydrate

The boundaries of the existence of gas hydrates are known with fairly high accuracy from direct measurements and models of chemical thermodynamics [1]. Self-preservation of gas hydrates is the abnormal stability of the hydrate in the region below the dissociation curve. The mechanisms of self-preservation are not completely understood [12]. It is known that the cause of abnormal stability is often an impermeable ice crust on the surface of particles. This crust can be formed both due to the liquid water contained in the system and due to the crystallization of supercooled water, which is formed during the hydrate decomposition [60, 61]. Surface microscopy shows that during self-preservation, the surface of particles is rearranged from a porous structure to a smooth continuous crust [62].

Using the molecular dynamics method, the authors of [63], [64] showed that the decomposition of hydration cells can be difficult due to the difference in the thermal expansion coefficients of ice and hydrate.

In [49], [65-68], the phenomenon of self-conservation is considered as a diffusion process with a damping diffusion coefficient (see also simplified diffusion-kinetic model with apparent mass transfer coefficient [69]). The author constructed an isothermal theory of self-preservation for particles of different geometries. Similar ideas were proposed in our studies [70, 71] of non-isothermal processes. A decrease in the transfer coefficients in a porous medium often means changes in its structure, namely, in the sizes and connectivity of cracks and pores of the ice crust. In this regard, it is interesting to study the relationship between the size and distribution of pores, and the dynamics of their change. The studies in [72, 73] proposed models for a developing pore system in carbonaceous fuels. It seems reasonable to use the available models for the inverse process associated with narrowing and blocking of pores [74, 75]. For example, exponent in time damping function used in [66], [68] is to have a

1

0.9 0.8 0.7

cO °r 0.6

CT

e

/.5

c0 CL

Ol 0.4 0.3 0.2 0.1 0

0.2

0.4

0.6

0.8

0

0.2

0.4

0.6

0.8

Fig. 8. The dependence of the solution to equation (30) on the Fig. 9.The dependence of the solution to equation (31) on the fraction of unreacted hydrate X and the ratio kF/kR (values are fraction of unreacted hydrate X and the ratio kF/kR (values are

shown on the corresponding curves).

(30) and 1/8 for (31). Therefore, when using the thin crust approximation, we obtain peak pressure in the pores after the hydrated particle dissociates by about 70%. Given the spherical curvature of the particle, we obtain peak pressure in the later stages of decomposition (about 90%).

The self-preservation range may be considered as regions of low permeability curves in Figs. 8 and 9, when the pressure in the pores is close to the equilibrium value for a given temperature most of the conversion time. It can be said that the dissociation with self-preservation (at least in one of the variants of this phenomenon) is limited by filtration. With a sharp change in permeability, the solution moves from one curve to another.

The permeability of the ice crust, however, is uncertain. Its direct measurement is challenging, therefore, it is necessary to look for suitable assessment methods. As suggested above in (8), permeability can be considered (as a first approximation) in the form of a product of pore density and hydraulic resistance of a single pore. Then we need to measure the average pore diameter and their density per unit of surface. Both quantities, however, can change in the process of dissociation. Based on the surface microscopy studies of hydrates, the following values were taken: an average pore diameter equal to 1 ^m and a pore density equal to 1011 m-2 (pores occupy about 8% of the surface).

In the region of self-preservation, the density and size of pores change quite sharply due to multiple local temperature inhomogeneities that cause phase transitions. Thus, a simulation model of self-preservation can be proposed: when a particle enters a chosen temperature range, its permeability characteristics drop to small (compared to the initial) values. In the calculations, this drop occurs stepwise, after which the permeability begins to increase (due to the development of cracks and thermomechanical stresses).

shown on the corresponding curves).

physical sense related to thermomechanical properties of the ice crust.

Model (1-9) can be used to construct a qualitative theory of self-preservation based on the idea of varying permeability of the ice crust. Under the constant thermobaric conditions, we can write the equations for the stationary pressure in the pores of a particle (using thin crust approximation (11)):

kR (Peq - P)X23 - kF (P - Pu) * = 0 (30)

1 - X/s

For different values of the fraction of unreacted hydrate X, this equation gives different values of quasi-stationary pressure level. The natural parameter in this equation is the ratio of the kinetic coefficients kR and kF, which determine the rates of dissociation and filtration processes, respectively.

Applying stationary flow approximation (12), we obtain the following equation:

2/ 1

kR (Peq -P)X23 -kF (P -Pout= 0 (31)

XV3 -1

The dependence of the dimensionless pressure in the pores of the particle on the ratio kF/kR and the fraction of unreacted hydrate X is shown in Fig. 8 for (30) and Fig. 9 for (31). The qualitative behavior of the solution does not depend on the approximation for the filtration resistance. With a high permeability of the crust, the pressure in the pores practically does not differ from the ambient pressure, however, even with comparable dissociation and filtration rate coefficients, the pressure increase becomes significant. The position of the maximum pressure can be determined by differentiating the expression for P with respect to X and equating the derivative to zero: the value of X at the peak of pressure does not depend on the coefficients kR and kF. Maximum pressure point Xmaxis equal to 8/27 for

0

Fig. 10. Kinetic curves of dissociation of gas hydrates: 1 and 3 - natural gas hydrate; 2 and 4 - synthetic gas hydrate (lines - calculation, markers - measurements).

1.00E-07

200 400 600 1.00E-09 -

300 700 1100 1500 1900 2300 2700 1.00E-09 -

1.00E-11

1.00E-13 -

1.00E-15 -

1.00E-17

Fig. 11. Gas emission rate during dissociation of methane hydrate, kg/s: a - synthetic gas hydrate, b - natural gas hydrate (solid lines - calculation, markers - measurements).

To validate the mathematical model, dissociation of methane hydrate of different origins (synthetic and natural) was studied under the conditions of gas combustion above a layer of powder [76]. The average particle size was 1.2 mm for synthetic hydrate and 2.5 mm for natural one. The calculations were carried out for the experimental conditions, given permeability as a function of time. The comparison of the calculated and experimental data is shown in Figs. 10 and 11.

500 1000 1500 2000 2500 3000

Time, s

0 500 1000 1500 2000 2500 3000

Time, s

Fig. 12. Calculated parameters of dissociating artificial hydrate particle: a - pressure, b - temperature.

The change in permeability is specified through a change in the pore density (we assume that the pore size varies slightly under the experimental conditions). The initial density is 1011 m-2, but when it goes into the region of self-preservation, it drops to 103 m-2. The synthetic gas hydrate remains in the self-preservation conditions for a relatively short time, therefore, pore density in the calculation decreases and grows stepwise, which allows a fairly close reproduction of the measured kinetic curve (Fig. 10). Particles of natural gas hydrate are 2 times larger, so the duration of its decay, including the time spent in the self-preservation conditions, lasts longer. A stepwise change in pore density, in this case, gives high deviations; therefore, attempts were made to find a better approximation for the dynamics of permeability reduction. As the calculations show, the best agreement with the experiment is provided by a linear increase in the logarithm of pore density, i.e. exponential growth in the number of pores in time (Fig. 11). This dependence can be explained in terms of the development of random cracks.

The change in temperature and gas pressure in the particle can be seen in Fig. 12: when the site of self-preservation occurs, the pressure in the pores increases sharply, reaching equilibrium values, and dissociation

T

Fig. 13. Hydrate powder layer dissociation model.

practically stops. After returning to the initial permeability, the mass-loss rate becomes high due to the high rate of transfer processes after depressurization in the pores of the particle. In the case of natural hydrate, the rate of dissociation, taking into account self-preservation, slows down more smoothly, due to a larger particle size and, hence, higher filtration resistance of the ice crust (Fig. 11). Thus, using a fairly simple model with dynamic permeability, we can reproduce the phenomenon of methane hydrate self-preservation and study some details of this process associated with the internal state of the particle (apparent kinetics of dissociation, gas pressure in the pores).

VI. Mathematical modeling of dissociation of a gas hydrate layer in a one-dimensional approximation

The above calculations were based on the approximation of the thermal uniformity of particles and the possibility of representing the entire powder layer as a single particle. This approximation is applicable to thin layers, and with an increasing thickness (which is especially important for technological applications) it is necessary to take into account the inhomogeneity of heating along the layer height. It is especially important for the combustion of gas hydrates when temperature gradients could be large [77, 78]. In this section, we consider a problem in the one-dimensional approximation, where the spatial coordinate is the layer height. The problem statement is presented in Fig. 13.

A layer of powder of height L on a thermally insulated surface is considered. At the initial moment, the temperature of the powder is Tin, and the gas pressure in the porous space is atmospheric (POut). Hydrate is heated by convective heat transfer from the ambient air with a temperature TOut. When the dissociation conditions are reached, the gas hydrate decomposes, and the released gas flows between the particles and leaves the powder layer. The perturbations introduced by the gas flow over the powder layer, the effect of self-preservation, and the influence of gravity are neglected. We use a single-temperature approximation (heat

transfer between gas and particles occurs faster than heat distribution due to heat conduction and gas movement). Then the system of heat and mass transfer equations inside the layer is written as follows:

C rn+c; r I1 -n)fT =

dt

=— W. ru dT+Q dr

dz Ç eff dz J pP dz * dt

U=-kD dP

ß dz

r=- K (P - P ) r

P -Pend

y

(32)

(33)

(34)

v p Pend ,

Here Peq is the equilibrium gas pressure above the hydrate at a given temperature, and Kd is the apparent dissociation kinetic coefficient, which is determined by the expression:

Kd = ^expl--^

RT0 Bdr

(35)

The boundary and initial conditions are as follows:

P (0, t) = Pu, (36)

-l

ÔP dz

eff

5T_

dZ

= a (T0u, - T)

= 0;

ST dz

= 0

P ( z,0) =Pout ; T ( z, 0) =Tin

(37)

(38)

(39)

The calculations were carried out with the following parameters: a layer height is 20 mm, an average particle size is 1 mm, porosity is 0.4, the initial temperature of the powder is 190 K, and the outer temperature is 473 K. The heat transfer coefficient of the layer surface is 10 W/m2/K. The effective coefficient of thermal conductivity of powder was calculated according to the recommendations [50]. Hydrate and ice phases have a thermal conductivity of about 0.2-0.5 W/m/K but heating of the whole layer occurs through porous space filled with gas and contact areas between particles, i.e. Biot number for a single particle is still small but Biot number for a layer of particles could be very large.

A numerical algorithm used for the calculations was previously developed for a mathematical model of pyrolysis of woody fuel particles [58]. Each moment is divided into several calculation stages. First, a kinetic step is made, gas sources and heat sinks are calculated. In the next step, the pressure and velocity fields are calculated. Finally, the problem of thermal conductivity is solved. The calculation results are shown in Figs. 14 and 15. It is possible to averagely divide the layer height into sections corresponding to particles and construct the dynamics of heating and dissociation of particles in these sections. Since the layer height is 20 particle diameters, the curves are plotted for these sections with reference to the layer (particle) number, and the calculation is carried out starting from the bottom of the tank.

z=0

Time, sec

Fig. 14. Dynamics of temperature changes in a layer of gas hydrate particles.

\ N ' ■ \ N ' - \ N \ S \ N ~i(T

\ \ \ \ ' • . 15

\ \

\ * 18 \20 ^ \ I \ 1 1 I I 1

I 1 1 1 1

0 50 100 150 200 250 300 350 400 450 500 Time, sec

Fig. 15. Dynamics of unreacted hydrate fraction changes in a layer of powder.

To simulate the heating of a layer from the air, it is necessary to take into account the melting of ice at temperatures above 273 K. Therefore, in addition to the equations of thermal conductivity and filtration, the Stefan condition for the phase transition is included in the system. Considering the scheme of layer dissociation similar to that used in [79], we can neglect the effect of water drainage.

Ice melting occurs as a monotonic process of heat absorption at a constant temperature: one can see in Fig. 14 the areas, where the temperature remains constant for long periods. In contrast to melting, dissociation is characterized by significant temperature fluctuations, including those associated with collective effects. Dissociation begins when thermobaric stability limits are reached but it is kinetically limited at low temperatures. Therefore, the heating of the lower layers is faster than the complete dissociation of the upper layers. When dissociation in the 18th particle begins, the surface still does not have time to decompose completely (see Fig. 15). As a result, heat sink intensifies, and sharp minima are observed on the temperature curves. Even the melting of ice on the surface (time point of about 250 s) is interrupted due to the intense absorption of heat by dissociating inner layers. The developed mathematical model will allow us to calculate the variance in the particle decomposition rate in different sections of the layer, find reasonable averaging methods, and compare them with the ones used above.

VII. Conclusions 1. The paper proposes a mathematical model for the dissociation of gas hydrates, which is based on the hypothesis of a filtration mechanism of self-preservation. The model takes into account the kinetics of heterogeneous decomposition and gas transport in the pores of a particle. The presented model is shown to distinguish several conditions of decomposition of

hydrate particles. The establishment of these conditions depends on the characteristics of the porous structure of the ice crust. With a high pore density, a kinetic decomposition occurs, and at a low pore density, filtration is involved. For intermediate values of pore density during the decomposition process, the conditions change.

2. Based on the experimental data analysis conducted by different scientific teams, new kinetic coefficients of dissociation of methane hydrate at low temperatures were obtained. The activation energy of the process of low-temperature dissociation (34 kJ/mol) is much lower than the known activation energy of dissociation at temperatures above ice melting, which may be due to a change in the dissociation mechanism.

3. The developed model and new kinetic coefficients were used to reproduce the experimental kinetic features of the methane hydrate dissociation in a thin powder layer: the dependence of the dissociation kinetics on time, particle size, and heat fluxes.

4. A one-dimensional unsteady mathematical model of dissociation of a layer of gas hydrate particles was developed taking into account heat transfer and filtration in a porous volume. The model allows us to investigate the features of the collective behavior of gas-emitting particles, the absorption of dissociation heat, and to evaluate inhomogeneities along the layer height.

Acknowledgments

The theoretical calculations were carried out within the framework of the Melentiev Energy Systems Institute SB RAS Basic Research Program (project No. 11.17.1.2, reg. No. AAAA-A17-117030310448-0). The experimental studies were carried out within the framework of the state assignment for the Kutateladze Institute of Thermophysics SB RAS (Project No. AAAA-A17-117080810004-8).

References

[1] E.D. Sloan and C.A. Koh, Clathrate hydrates of natural gases (3rd ed.). Boca Raton, USA: CRC Press, Taylor and Francis Group, 2008.

[2] E.A. Bondarev, VI. Vasilyev, A.F. Voevodin, N.N. Pavlov and A.P.Shadrina, Thermohydrodynamics of gas mining and transport systems. Novosibirsk: Nauka, 1988. (in Russ.)

[3] Yu.F. Makogon, Hydrates of hydrocarbons. Tulsa, Oklahoma: Pennwell Books, 1997.

[4] N.G. Musakaev, R.R. Urazov, and V.Sh. Shagapov, "Hydrate formation kinetics in piped natural-gas flows, "Thermophys. Aeromech., vol. 13, no. 2, pp. 275-281, 2006. DOI: 10.1134/S0869864306020090

[5] A.A. Makarov (ed.), Forecast of the world and Russian energy sector. Energy Research Institute of RAS, Analytical Center of RF Government. Moscow, 2016.

[6] The International Energy Agency World Energy Outlook 2016. IEA Publications, 2016.

[7] Z. Yin and P. Linga, "Methane hydrates: A future clean energy resource," Chinese J. Chem. Eng., vol. 27, no. 9, pp. 2026-2036, 2019. DOI: 10.1016/j. cjche.2019.01.005

[8] Gas hydrates: extraction technologies and development prospects. Analytical Center of RF Government, 2013.

[9] Y. Cui, C. Lu, M. Wu, Y. Peng, Y. Yao and W. Luo, "Review of exploration and production technology of natural gas hydrate," Adv. Geo-Energy Res., vol. 2, no. 1,pp. 53-62, 2018. DOI: 10.26804/ager.2018.01.05

[10] T. Yu, G. Guan, and A. Abudula, "Production performance and numerical investigation of the 2017 offshore methane hydrate production test in the Nankai Trough of Japan," App. Energy, vol. 251, no. 113338, 2019. DOI: 10.1016/j.apenergy.2019.113338

[11] Y. Sun, X. Ma, W. Guo, R. Jia, and B. Li, "Numerical simulation of the short- and long-term production behavior of the first offshore gas hydrate production test in the South China Sea," J. Petr. Sci. Eng., vol. 181, no. 106196, 2019. DOI: 10.1016/j.petrol.2019.106196

[12] Y.F. Makogon, "Natural gas hydrates - A promising source of energy," J. Nat. Gas Sci. Eng. , vol. 2, no. 1, pp. 49-59, 2010. DOI: 10.1016/j.jngse.2009.12.004

[13] Z.R. Chong, S.H.B. Yang, P. Babu, P. Linga, and X.-S. Li, "Review of natural gas hydrates as an energy resource: Prospects and challenges," App. Energy, vol. 162, pp. 1633-1652, 2016. DOI: 10.1016/j. apenergy.2014.12.061

[14] A. Falenty, W.F. Kuhs, M. Glockzin and G. Rehder, "Self-preservation of CH4 hydrates for gas transport technology: pressure-temperature dependence and ice microstructures," Energy Fuels, vol. 28, pp. 62756283, 2014. DOI: 10.1021/ef501409g

[15] H.P. Veluswamy, A. Kumar, Y. Seo, J.D. Lee, P. Linga, "A review of solidified natural gas (SNG) technology for gas storage via clathrate hydrates," App. Energy, vol. 216, pp. 262-285, 2018. DOI: 10.1016/j. apenergy.2018.02.059

[16] G. Rehder, R. Eckl, M. Elfgen, A. Falenty, R. Hamann, N. Kahler, W.F. Kuhs, H. Osterkamp and C. Windmeier, "Methane hydrate pellet transport using the self-preservation effect: a techno-economic analysis," Energies, vol. 5, pp. 2499-2523, 2012. DOI: 10.3390/ en5072499

[17] C. Xu, X. Li, K. Yan, X. Ruan, Z. Chen and Z. Xia, "Research progress in hydrate-based technologies and processes in China: A review," Chinese J. Chem. Eng., vol. 27, no. 9, pp. 1998-2013, 2019. DOI: 10.1016/j. cjche.2018.12.002

[18] S.Y. Misyura, "The influence of porosity and structural parameters on different kinds of gas hydrate dissociation," Sci. Rep., vol. 6, no. 30324, 2016. DOI: 10.1038/srep30324

[19] V.Sh. Shagapov, M.K. Khasanov, N.G. Musakaev, and N.H. Duong, "Theoretical research of the gas hydrate deposits development using the injection of carbon dioxide," Internat. J. Heat Mass Trans., vol. 107, pp. 347-357, 2017. DOI: 10.1016/j. ijheatmasstransfer.2016.11.034

[20] V.Sh. Shagapov, G.R. Rafikova, and M.K. Khasanov, "The Theory of the Replacement of Methane by Carbon Dioxide in Gas Hydrates," Theor. Found. Chem. Eng., vol. 53, pp. 64-74, 2019. DOI: 10.1134/ S0040579518060118

[21] E.D. Sloan, "Fundamental principles and applications of natural gas hydrates," Nature, vol. 426, pp. 353-363, 2003. DOI: 10.1038/nature02135

[22] A.Yu. Manakov, N.V. Penkov, T.V. Rodionova, A.N. Nesterov, and E.E. Fesenko, "Kinetics of formation and dissociation of gas hydrates," Russ. Chem. Rev., vol. 86, p. 845, 2017. DOI: 10.1070/RCR4720

[23] W. Ke, T.M. Svartaas, and D. Chen, "A review of gas hydrate nucleation theories and growth models," Journal of Natural Gas Science and Engineering, vol. 61, pp. 169-196, 2019. DOI: 10.1016/j.jngse.2018.10.021

[24] Z. Yin, Z.R. Chong, H.K. Tan, and P. Linga, "Review of gas hydrate dissociation kinetic models for energy recovery," J. Nat. Gas Sci. Eng., vol. 35B,pp. 13621387, 2016. DOI: 10.1016/j.jngse.2016.04.050

[25] S.D. Sleptsov and V.E. Dontsov, "Modelling of decomposition of a single gas hydrate particle in the water behind the shock wavefront," Thermophysics and Aeromechanics, vol. 17, no. 3, pp. 375-381, 2010. DOI: 10.1134/S086986431003008X

[26] S. Takeya, T. Ebunima, T. Uchida, J. Nagao, and H. Narita, "Self-preservation effect and dissociation rates of CH4 hydrate," J. Crystal Growth, vols. 237-239, pp. 379-382, 2002. DOI: 10.1016/S0022-0248(01)01946-7

[27] V.A. Vlasov, "Formation and dissociation of gas hydrate in terms of chemical kinetics," React. Kin. Mech. Cat., vol. 110, no. 1, pp. 5-13, 2013. DOI: 10.1007/s11144-013-0578-x

[28] T. Bar-Kohany and W.A. Sirignano, "Transient combustion of a methane-hydrate sphere," Comb. Flame, vol. 163, pp. 284-300, 2016. DOI: 10.1016/j. combustflame.2015.10.004

[29] H. Liang, Y. Song and Y. Chen, "Numerical simulation for laboratory-scale methane hydrate dissociation by depressurization," Energy Conv. Manag., vol. 51, pp. 1883-1890, 2010. DOI: 10.1016/j. enconman.2010.02.018

[30] M. Clarke and P.R. Bishnoi, "Determination of the activation energy and intrinsic rate constant of methane gas hydrate decomposition," Can. J. Chem. Eng., vol. 79, pp. 143-147, 2001. DOI: 10.1002/cjce.5450790122

[31] S.A. Bagherzadeh, S. Alavi, J.A. Ripmeester, and P. Englezos, "Evolution of methane during gas hydrate dissociation," Fluid Phase Equil., vol. 358, pp. 114120, 2013. DOI: 10.1016/j.fluid.2013.08.017

[32] R.V. Belosludov, Y.Y. Bozhko, R.K. Zhdanov, O.S. Subbotin, Y. Kawazoe and VR. Belosludov, "Hydrogen hydrates: Equation of state and self-preservation effect," Fluid Phase Equil., vol. 413, pp. 220-228, 2016. DOI: 10.1016/j.fluid.2015.11.031

[33] J. Kondori, S. Zendehboudi, and L. James, "New insights into methane hydrate dissociation: Utilization of molecular dynamics strategy," Fuel, vol. 249, pp. 264-276, 2019. DOI: 10.1016/j.fuel.2019.02.125

[34] C. Windmeier and L.R. Oellrich, "Visual observation of the methane hydrate formation and dissociation process," Chem. Eng. Sci., vol. 109, pp. 75-81, 2014. DOI: 10.1016/j.ces.2014.01.018

[35] R.I. Nigmatullin, VSh. Shagapov and VR. Syrtlanov, "Self-similar problem of decomposition of gas hydrates in a porous medium upon depression and heating," J. App. Mech. Tech. Phys., vol. 39, pp. 421-427, 1998. DOI: 10.1007/BF02468125

[36] M.J. Castaldi, Y. Zhou and T.M. Yegulalp, "Down-hole combustion method for gas production from methane hydrates," J. Petr. Sci. Eng., vol. 56, pp. 176-185, 2007. DOI: 10.1016/j.petrol.2006.03.031

[37] Y. Liu, Y. Chen, Y. Bai and S. Li, "Similarity theory for the physical simulation of natural gas hydrate reservoir development," Mining Sci. Tech., vol. 20, pp. 782-788, 2010. DOI: 10.1016/S1674-5264(09)60281-7

[38] Y. Liu and I.K. Gamwo, "Comparison between equilibrium and kinetic models for methane hydrate dissociation," Chem. Eng. Sci., vol. 69, pp. 193-200, 2012. DOI: 10.1016/j.ces.2011.10.020

[39] E.W. Thiele, "Relation between Catalytic Activity and Size of Particle," Ind. Eng. Chem., vol. 31, no. 7, pp. 916-920, 1939. DOI: 10.1021/ie50355a027

[40] Ya.B. Zeldovich, Selected works. Volume I: Chemical Physics and Hydrodynamics. Princeton, New Jersey: Princeton Univ. Press, 1992.

[41] M.G. Slin'ko, V.A. Kirillov, I.A. Mikhailova, and S.I. Fadeev, "Mathematical Model of a Catalytic Process on a Porous Grain in a Gas-Liquid-Solid Three-Phase System," Dokl. Chem., vol 376, pp. 42-45, 2001. DOI: 10.1023/A:1018816527288

[42] O. Senneca, P. Bareschino. M. Urciuolo and R. Chirone, "Prediction of structure evolution and fragmentation phenomena during combustion of coal: Effects of heating rate," Fuel Proc. Tech., vol. 166, pp. 228-236,

2017. DOI: 10.1016/j.fuproc.2017.06.010

[43] S.L. Singer and A.F. Ghoniem, "Comprehensive gasification modeling of char particles with multi-modal pore structures," Comb. Flame, vol. 160, pp. 120-137, 2013. DOI: 10.1016/j.combustflame.2012.09.007

[44] E. Chuvilin and B. Bukhanov, "Thermal Conductivity of Frozen Sediments Containing Self-Preserved Pore Gas Hydrates at Atmospheric Pressure: An Experimental Study," Geosciences, vol. 9, no. 65, 2019. DOI: 10.3390/geosciences9020065

[45] Y.P. Handa, "Compositions, enthalpies of dissociation, and heat capacities in the range 85 to 270 K for clathrate hydrates of methane, ethane, and propane, and enthalpy of dissociation of isobutane hydrate, as determined by a heat-flow calorimeter," J. Chem. Thermodynamics, vol. 18, pp. 915-921, 1996. DOI: 10.1016/0021-9614(86)90149-7

[46] S.A. Aromada, B. Kvamme, N. Wei, and N. Saeidi, "Enthalpies of Hydrate Formation and Dissociation from Residual Thermodynamics," Energies, vol. 12, no. 4726, 2019. DOI: 10.3390/en12244726

[47] A.H.S. Dehaghani and B. Karami, "A new predictive thermodynamic framework for phase behavior of gas hydrate," Fuel, vol. 216, pp. 796-809, 2018. DOI: 10.1016/j.fuel.2017.11.128

[48] D.B. Spalding, Combustion and mass transfer. Oxford: Pergamon Press, 1979.

[49] V.A. Vlasov, "Diffusion model of gas hydrate dissociation into ice and gas: Simulation of the self-preservation effect," Internat. J. Heat Mass Trans., vol. 102, pp. 631-636, 2016. DOI: 10.1016/j. ijheatmasstransfer.2016.06.057

[50] M.E.Aerov, O.M.Todes, and D.A.Narinskii, Stationary granular beds: hydraulics and heat transfer. Leningrad: Khimiya, 1979. (in Russ.)

[51] A. Fujimoto and T. Sugahara, "Scanning electron microscopic studies on the methane hydrate decomposition using the freeze-fracture replica method," Bull. Glaciolog. Res., vol. 35, pp. 39-45, 2017. DOI: 10.5331/bgr.17R02

[52] H.C. Kim, P.R. Bishnoi, R.A. Heidemann, and S.S.H. Rizvi, "Kinetics of methane hydrate decomposition," Chem. Eng. Sci., vol. 42, no. 7, pp. 1645-1653, 1987. DOI: 10.1016/0009-2509(87)80169-0

[53] S. Takeya, T. Uchida, J. Nagao, R. Ohmura, W. Shimada, Y. Kamata, T. Ebunima, and H. Narita, "Particle size effect of CH4 hydrate for self-preservation, "Chem. Eng. Sci., vol. 60, no. 5, pp. 1383-1387, 2005. DOI: 10.1016/j.ces.2004.10.011

[54] L.A. Stern, S. Circone, S.H. Kirby and W.B. Durham, "Anomalous preservation of pure methane hydrate at 1 atm," J. Phys. Chem. B, vol. 105, pp. 1756-1762, 2001. DOI: 10.1021/jp003061s

[55] S.Ya. Misyura, "Effect of heat transfer on the kinetics of methane hydrate dissociation, "Chem. Phys. Lett., vol. 583, pp. 34-37, 2013. DOI: 10.1016/j. cplett.2013.08.010

[56] E. Chuvilin, B. Buhanov, O. Guryaeva, V. Istomin, S. Takeya, and A. Hachikubo, "Experimental study of self-preservation mechanisms during gas hydrate decomposition in frozen sediments," Proc. 7th Internat. Conf. Gas Hydrate (ICGH 2011). Edinburgh, Scotland, UK, July 17-21 2011.

[57] E. Ranzi, M. Corbetta, F. Manenti and S. Pierucci, "Kinetic modeling of the thermal degradation and combustion of biomass," Chem. Eng. Sci., vol. 110, pp. 2-12, 2014. DOI: 10.1016/j.ces.2013.08.014

[58] I.G. Donskoy, "Mathematical modeling of woody particles pyrolysis in a fixed bed, " Computational Technologies, vol. 23, no 6, pp. 14-24, 2018. DOI: 10.25743/ICT.2018.23.6.003 (in Russ.)

[59] S.Y. Misyura and V.S. Morozov, "Influence of the granule size and composition uniformity on methane hydrate dissociation," J. Phys. Conf. Ser., vol. 1128, no. 012072, 2018. DOI: 10.1088/17426596/1128/1/012072

[60] VP. Mel'nikov, A.N. Nesterov and A.M. Reshetnikov, "Formation of supercooled water upon dissociation of propane hydrates at T < 270 K," Dokl. Phys. Chem., vol. 417, pp. 304-307, 2007. DOI: 10.1134/ S0012501607110048

[61] VP. Mel'nikov, A.N. Nesterov, L.S. Podenko, A.M. Reshetnikov, and V.V. Shalamov, "Metastable states of gas hydrates under pressures lower than equilibrium pressures of the ice-hydrate-gas system," KriosferaZemli, vol. 15, no. 5, pp. 80-83, 2011. (in Russ.)

[62] A. Falenty and W.F. Kuhs, "Self-preservation of CO2 gas hydrates - surface microstructure and ice perfection," J. Phys. Chem. B, vol. 113, pp. 1597515988, 2009. DOI: 10.1021/jp906859a

[63] O.S. Subbotin, V.R. Belosludov, T. Ikeshoji, E.N. Brodskaya, E.M. Piotrovskaya, V. Sizov, R.V Belosludov, and Y. Kawazoe, "Modeling the Self-Preservation Effect in Gas Hydrate/Ice Systems," Mat. Trans., vol. 48, no. 8, pp. 2114-2118, 2007. DOI: 10.2320/matertrans.N-MRA2007866

[64] VR.Belosludov, Y.Y.Bozhko, and R.K. Zhdanov, "Self-preservation effect modeling in hydrate systems using Lattice Dynamic methods," J. Phys.: Conf. Ser., vol. 1128, no. 012086, 2018. DOI: 10.1088/17426596/1128/1/012086

[65] V.A. Vlasov, "Phenomenological diffusion theory of the formation of gas hydrate from ice powder," Theor. Found. Chem. Eng., vol. 46, pp. 576-582, 2012. DOI: 10.1134/S0040579512060243

[66] A.A. Sizikov, V.A. Vlasov, A.S. Stoporev, A.Y. Manakov, "Decomposition Kinetics and Self-Preservation of Methane Hydrate Particles in Crude Oil Dispersions: Experiments and Theory," Energy Fuels, vol. 33, no. 12, pp. 12353-12365, 2019. DOI: 10.1021/ acs.energyfuels.9b03391

[67] V.A. Vlasov, "Mathematical Model of the Effect of Self-Preservation of Gas Hydrates," J. Eng. Phys. Thermophys., vol. 92, no. 6, pp. 1406-1414, 2019. DOI: 10.1007/s10891-019-02057-8

[68] V.A. Vlasov, "Diffusion model of gas hydrate dissociation into ice and gas that takes into account the ice microstructure," Chem. Eng. Sci., vol. 215, no. 115443, 2020. DOI: 10.1016/j.ces.2019.115443

[69] W. Liu, Q. Li, Y. Song, L. Zhang, M. Yang, L. Wang, and Y. Chen, "Diffusion Theory of Formation of Gas Hydrate from Ice Powder without Melting," Energy Procedia, vol. 61, pp. 513-522, 2014. DOI: 10.1016/j. egypro.2014.11.1161

[70] I.G. Donskoy and S.Ya. Misyura, "Dissociation model of gas hydrate particle with consideration of methane filtration in the pores of ice crust," Modern Science. Researches, Ideas, Results, Technologies, no. 1 (16), pp. 241-245, 2015.

[71] S.Y. Misyura and I.G. Donskoy, "Dissociation of natural and artificial gas hydrate," Chem. Eng. Sci., vol. 148, pp. 65-77, 2016. DOI: 10.1016/j.ces.2016.03.021

[72] G.R. Galavas, "A random capillary model with application to char gasification at chemically controlled rates," AIChE J., vol. 26, no. 4, pp. 577-585, 1980. DOI: 10.1002/aic.690260408

[73] S.K. Bhatia and D.D. Perlmutter, "A random pore model for fluid-solid reactions: II. Diffusional and transport effects," AIChE J., vol. 27, no. 2, pp. 247254, 1981. DOI: 10.1002/aic.690270211

[74] I.G. Donskoi, "Process simulation ofthe co-gasification of wood and polymeric materials in a fixed bed," Solid Fuel Chem., vol. 52, no. 2, pp. 121-127, 2018. DOI: 10.3103/S0361521918020027

[75] I.G. Donskoy, "Mathematical modeling of the agglomeration in a reactive porous medium with variable permeability," Computational Technologies, vol. 25, no. 2, 2020. DOI: 10.25743/ICT.2020.25.2.003 (in Russ.)

[76] S.Y. Misyura, "Dissociation of natural and artificial methane hydrate," EPJ Web of Conf., vol. 110, no. 01040, 2016. DOI: 10.1051/201611001040

[77] G. Cui, S. Wang, Z. Dong, X. Xing, T. Shan, and Z. Li, "Effects of the diameter and the initial center temperature on the combustion characteristics of methane hydrate spheres," App. Energy, vol. 257, no. 114058, 2020. DOI: 10.1016/j.apenergy.2019.114058

[78] S.Y. Misyura, "Non-stationary combustion of natural and artificial methane hydrate at heterogeneous dissociation," Energy, vol. 181, pp. 589-602, 2019. DOI: 10.1016/j.energy.2019.05.177

[79] F.H. Wu, R.E. Padilla, D. Dunn-Rankin, G.B. Chen, and Y.C. Chao, "Thermal structure of methane hydrate fueled flames," Proc. Comb. Inst., vol. 36, no. 3, pp. 4391-4398, 2017. DOI: 10.1016/j.proci.2016.06.012

Igor Donskoy received the Ph.D. degree in energy systems in 2014. He is currently a senior researcher in the Laboratory of Thermodynamics at Melentiev Energy Systems Institute SB RAS (Irkutsk). His main research interests include mathematical modeling, fuel processing, and thermal engineering.

Sergey Misyura received the Ph.D. degree in thermophysics in 2013. He is currently a researcher in the Laboratory of Transfer Processes at Kutateladze Institute of Thermophysics SB RAS (Novosibirsk). His main interests include thermophysical experiment, heat and mass transfer during phase transitions, and gas hydrates processing.