Application of Population Balance Theory for Dynamic Modeling of

Subscriber access provided by UNIVERSITY OF TOLEDO LIBRARIES

Fossil Fuels

Application Of Population Balance Theory For Dynamic Modeling Of Methane And Ethane Hydrate Formation Processes Mohammad Saeed Khatami, and Akbar Shahsavand Energy Fuels, Just Accepted Manuscript • DOI: 10.1021/acs.energyfuels.8b01351 • Publication Date (Web): 12 Jun 2018 Downloaded from http://pubs.acs.org on June 12, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Application Of Population Balance Theory For Dynamic Modeling Of Methane And Ethane Hydrate Formation Processes Mohammad S. Khatami and Akbar Shahsavand*a a

Chemical Engineering Department, Faculty of Engineering, Ferdowsi University of Mashhad

Abstract: Natural gas hydrates are considered as a promising fuel source in a relatively near future. Kinetic modelling of various steps of natural gas hydrate formation process, such as dissolution, nucleation and growth processes have been received numerous attentions. A novel mechanism is introduced for the entire nucleation and growth steps and a proper mathematical model is presented to estimate the gas consumption rate in a constant temperature and pressure process. The proposed model covers the entire dissolution, nucleation and growth stages. The combined Lax-Wendroff/Crank-Nicholson method is employed to solve the population balance equation for estimation of total surface area of evolved hydrate particles and corresponding particle size distributions. A special class of artificial neural network (known as Regularization network) is used to predict the solid-liquid equilibria. The proposed model is successfully validated using experimental data borrowed from literature for both methane and ethane hydrate formation processes. The simulation results indicate that for both methane and ethane species, the mole fractions in the bulk of liquid are often close to the corresponding concentrations at solid-liquid interfaces and decrease over the time during hydrate growth processes. It is clearly demonstrated that the overall resistance shifts from the nucleation reaction to mass transfer as the hydrate formation progresses. Keywords: Kinetic modelling, methane and ethane hydrates, population balance equation, Lax-Wendroff/Crank-Nicholson, artificial neural network

1 ACS Paragon Plus Environment

Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1. Introduction Nowadays, natural gas plays a major role in energy supply. Global demand for natural gas experiences a monotonically increasing trend due to growing request for clean energy. Worldwide proven natural gas reserves are around 200 trillion standard cubic meters while natural gas hydrates (NGH) reserves are estimated over 15 quadrillion standard cubic meters.1 The present estimates of methane (as the main constituent of NG) stored in the form of hydrate suggest that these clathrates have the potential to produce NG for at least 300 years based on the current global consumption rates.2 Each volume of NGH contains around 180 volumes of natural gas, at the standard conditions.3 Kinetic modellings of hydrate formation and dissociation processes have great significance in the investigation and analysis of many industrial applications which involves natural gas storage and transmission 4 or producing gas from methane hydrate reservoirs.5 Other applications of NGH such as CO2 sequestration 6, use of ionic hydrate in refrigeration systems7, well control and flow assurance in subsea transportation pipelines 8,9 are also reported in the literature. Most of models previously presented for the gas hydrate formation process are thermodynamically oriented and kinetic modelling studies have been received much less attention in the literature. Glew and Hagget pioneered in kinetic modelling of gas hydrate formation process.10 The hydrate production rate was modelled by considering both convective and conductive heat-transfer resistances. They reported that the exothermic hydrate formation reaction rate is directly proportional to the temperature difference between the reactor and its cooling bath. Vysniauskas and Bishnoi 11 studied the methane hydrate formation in a 500 cm3 volume semibatch stirred reactor operating at constant pressure. The gas consumption rate was experimentally measured for different operating temperatures and pressures. A three-step reaction mechanism was proposed and then used to compute the methane consumption rate during hydrate growth via a semi-empirical approach. Englezos et al. 12 also conducted several kinetically oriented experiments for the formation of methane and ethane hydrates in the same experimental set-up as of Vysniauskas and Bishnoi.11 The experiments were conducted at various temperatures and pressures in the range of 274282 K and 0.636-8.903 MPa, respectively. A theoretical model with one adjustable parameter was developed by combining the crystallization in the bulk of liquid phase and mass transfer at the gas-liquid interface. The following three steps mechanism was considered:


Page 2 of 41

Page 3 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

a) Diffusion of guest molecules from the gas-liquid interface to the bulk of liquid phase, b) Travel of these molecules from bulk of liquid to the hydrate–solution interface, and c) Reaction of water and guest molecule at the hydrate–solution interface. The collected experimental data were used to validate the model predictions. They also extended the above model later for hydrate formation of methane-ethane mixture.13 Their proposed model was practically used in the next decade as the most complete kinetic approach available in the open literature for hydrate formation process. Skovborg and Rasmussen

14

denoted that the Englezos et al. model over-predicts the gas

consumption rate for prolonged times. In particular, they pointed out that for both operating conditions considered by Englezos et al., systematically increasing deviations exist between calculated and experimental gas consumption rates. Skovborg and Rasmussen assumed that the entire mass transfer resistance during the hydrate formation lay in the diffusion of dissolved gas from the gas-liquid interface to the bulk of liquid. Hence, the difference between the mole fractions of solute at the gas-liquid interface and bulk of liquid was considered as the overall mass transfer driving force. They showed that such simple film theory can successfully reproduce the available laboratory data.14 Herri et al. agreed with the important role of the gas-liquid interface mass transfer resistance in the kinetics of hydrate formation process, as suggested by Skovborg and Rasmussen. However, they remarked that a more pronounce model should be developed for accurate description of the entire hydrate formation process, especially the particle growth section. They also conducted several kinetic experiments at constant pressures in a 1000 cm3 volume stirred reactor and used light scattering technique to measure the corresponding particle size distribution. The collected experimental data indicated that the stirring speed strongly affect both particle mean diameter and the corresponding concentration in the reactor.15 In their next article, a new kinetic model was presented which considered both nucleation and growth stages without validating with their previously obtained experimental data.16 Kashchiev and Firoozabadi have investigated the nucleation induction time and the corresponding driving force in terms of chemical potential of a hydrate building block. A general relation was

developed for super-saturation of a single component as the driving force in gas hydrate formation process. They reported that the thermodynamic conditions of gas dissolution into the bulk of liquid phase have significant impact on the driving force (∆𝜇).17

In their next article, two expressions were mathematically derived for the nucleus size of hydrate particles and nucleation rate as functions of super saturation (∆𝜇) in both homogeneous


Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

and heterogeneous nucleation's. The proposed correlations were successfully employed to predict the nucleation of methane hydrate in the presence of various additives.18 A year later in 2003, Kashchiev and Firoozabadi determined the induction time during the gas hydrate formation process and predicted the consumed moles of gas component in this period, in the presence of various kinetic inhibitors. Both above parameters were described in terms of super saturation and the entire proposed relations were validated for methane, ethane and cyclopropane hydrate formation processes.19 From 2008 to 2010, Bergeron and Servio conducted a series of experimental and theoretical studies to predict the consumption rates of various gases such as methane, propane and carbon dioxide during the hydrate formation process.20-23 In the first study, a combination of mass transfer resistances in the bulk of liquid and in the layer around the hydrate particles is used along with the hydrate growth resistance to estimate the propane consumption rate during constant pressure hydrate formation. They concluded that the mass transfer resistance around the hydrate particles can be ignored against the other two resistances.20 In the second study, the reaction rate constant (kr) was determined for the carbon dioxide hydrate formation process by neglecting the mass transfer resistance in the bulk of liquid.21 In their next study, Bergeron and Servio measured the gas concentration in the bulk of liquid at the onset of the hydrate growth stage and afterwards, for carbon dioxide and methane species. The concentrations of carbon dioxide and methane in the bulk of liquid phase was essentially remained constant during the entire growth stage. Furthermore, the corresponding mole fractions increased monotonically at the onset of the growth stage by increasing and decreasing the corresponding pressure and temperature, respectively. 22 In their last study, they reported that the methane hydrate reaction rate constant (kr) was essentially independent of the operating pressure. Three approaches, including experimental particle size distribution (PSD) measurements, semi experimental technique and theoretical method were compared to predict the total surface area of hydrate particles. Evidently, the experimental method performed more adequately than the other two methods on predicting the gas consumption rate.23 In a relatively comprehensive study, ZareNezhad et al. have been experimentally and theoretically investigated the kinetics of gas hydrate formation processes for pure methane, ethane and carbon dioxide at isochoric and isobaric conditions in the presence and absence of kinetic additives.24 A combination of growth rate and mass transfer around the hydrate particles were used. Similar to Herri and Kwaterski Langmuir type model, many system parameters such as driving force and free energy change were expressed in terms of liquid phase composition.25 They surprisingly concluded that "kinetic rate constant which is the only adjustable parameter 4 ACS Paragon Plus Environment

Page 4 of 41

Page 5 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

of the proposed model is a weak function of the temperature and a strong function of the surfactant concentration".24 Khosharay, Varaminian and Roosta conducted several constant volume experiments regarding methane and carbon dioxide hydrate formation. They essentially used the Skovborg and Rasmussen model with some modifications and successfully validated their proposed model with the collected experimental data.26 The interface concentrations were calculated by resorting to interfacial properties and parachor model instead of traditionally used vapor-liquid equilibrium (VLE) approach. The current study directly solves the population balance equation by employing the combined Lax-Wendroff/Crank-Nicholson method to determine the particle size distribution of clathrates. Mass transfer resistance in the bulk of liquid phase and the combination of nucleation and growth processes are simultaneously considered to predict the gas consumption rates during methane and ethane hydrate formation. A special class of artificial neural networks known as Regularization Network is used to compute the solid-liquid equilibrium methane concentration in the bulk of liquid phase. Various experimental data sets from literature are used to successfully validate the proposed model.12,27,28 2. Mathematical background Previous section reviewed various articles regarding kinetics and modelling of gas hydrate formation processes and the corresponding influential parameters. In this section, mathematical details of kinetic modeling for three individual parts of dissolution, nucleation and growth will receive proper attention. The actual hydrate formation process is considered to be the combination of both latter phenomena. A brief review of the solution method (known as LaxWendrof/Crank Nickelson) will be also included.

2.1. Kinetic Model The kinetic model of the hydrate formation process starts with the dissolution of gas into the liquid (water) phase. The following rate equation is used to predict the amount of mass transfer due to the dissolution process:16 𝑑𝑛 𝑑𝑡

= 𝑘𝐿 𝑎𝐶𝑤0 𝑉𝐿 (𝑥𝑖𝑛𝑡 − 𝑥𝑏𝑢𝑙𝑘 )

(1)


Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 41

𝑚2

Where, 𝑘𝐿 (𝑚⁄𝑠) is the mass transfer coefficient, 𝑎(𝑚3 ) is the specific mass transfer area, 𝐶𝑤0 (

𝑔𝑟𝑚𝑜𝑙⁄ 3 𝑚3 ) is initial water concentration in the liquid phase with a volume of 𝑉𝐿 (𝑚 ).

Furthermore, 𝑥𝑖𝑛𝑡 and 𝑥𝑏𝑢𝑙𝑘 denote mole fractions of the dissolved gas at the gas-liquid interface and the corresponding bulk of liquid phase. The value of 𝑘𝐿 𝑎 parameter could be found by resorting to solubility experimental data collected at the transient condition. Evidently, the interface concentration remains constant at fixed temperature and gas pressure, while the bulk concentration increases with time and will approach to the interface concentration, in the absence of hydrate formation. The final value of 𝑥𝑏𝑢𝑙𝑘 is restricted by the corresponding fugacity when hydrate formation is plausible. In other words, the fugacity of dissolved gas in the bulk of liquid and the prevailing hydrate phase should be equal at the start of nucleation process. At this step, a remarkable jump occurs in the bulk concentration of the liquid phase, as the hydrate formation initiates. The new gas concentration in the liquid phase can be easily computed by equating the consumption rates for three different resistances. As Figure 1 depicts, the first resistance in the liquid phase is due to the interface-bulk driving force and the second one belongs to the mass transfer resistance around the solid particles. The last resistance is due to the reaction rate at the solid-liquid interface.29 Evidently, for pure gases, the gas phase mass transfer resistance is essentially zero. By recognizing the three distinct resistances of Figure 1, the following equation is valid by assuming pseudo steady state condition during the hydrate formation: 𝑑𝑛 = 𝑘𝐿 𝑎 𝐶𝑤0 𝑉𝐿 (𝑥𝑖𝑛𝑡 − 𝑥𝑏𝑢𝑙𝑘 ) 𝑑𝑡 = 𝑘𝑑 𝐴𝑝+ 𝐶𝑤0 (𝑥𝑏𝑢𝑙𝑘 − 𝑥𝑏𝑝𝑎𝑟𝑡𝑖𝑐𝑙𝑒 ) = 𝑘𝑟 𝐴𝑝 (

𝑝𝑎𝑟𝑡𝑖𝑐𝑙𝑒 𝑒𝑞 𝜌𝐻𝑦𝑑 𝑥 −𝑥 ⁄𝑀𝑤 )( 𝑏 𝑥 𝑒𝑞 𝑏 ) 𝐻𝑦𝑑 𝑏

(2)

Where, 𝐴𝑝+ & 𝐴𝑝 denote the outer and inner areas of the liquid film around the hydrate particles. By neglecting the small thickness of the liquid film around the particle, then both above areas can be considered equal together (𝐴𝑝+ = 𝐴𝑝 ). Furthermore, 𝑘𝑑 (𝑚⁄𝑠) and 𝑘𝑟 (𝑚⁄𝑠( are mass transfer coefficient for liquid film around the particle and hydrate formation rate constant, respectively, 𝜌𝐻𝑦𝑑 (kg/m3) and 𝑀𝑤𝐻𝑦𝑑 are hydrate density and the corresponding molecular weight. Also, 𝑥𝑖𝑛𝑡 , 𝑥𝑏𝑢𝑙𝑘 , 𝑥𝑏𝑝𝑎𝑟𝑡𝑖𝑐𝑙𝑒 and 𝑥𝑏𝑒𝑞 denote various mole fractions of 6 ACS Paragon Plus Environment

Page 7 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

dissolved gas in the liquid phase at gas-liquid interface, bulk of liquid, outer and inner surfaces of the liquid film covering the solid particles. Evidently, it is assumed that 𝑥𝑏𝑒𝑞 is in equilibrium with hydrate at the prevailing gas pressure for pure gases. The following equation represents the relation between three individual resistances and the overall resistance (R):

𝑅=

1

1 𝐴𝑃

1

(𝑘 + 𝑘 ) + 𝑘 𝑑

𝑟

1

(3)

𝐿 𝑎𝑉𝐿

The mass transfer coefficient for liquid film around the particle could be computed from the following relation: 𝑘 𝐿

𝑆ℎ𝑝 = 𝐷𝑑

(4)

𝑒𝑓

Where, 𝐷𝑒𝑓 represents the effective diffusivity of dissolved methane in water at prevailing temperature and L is the diameter of the growing hydrate particles. In the absence of severe fluid motion and neglecting the Grashof number, the Sherwood number for mass transfer 2

coefficient over a spherical particle will reduce to 2.30 Since 𝐷𝑒𝑓 value is around 10−10 𝑚 ⁄𝑠, therefore the order of magnitude of mass transfer coefficient for liquid film around the particle (𝑘𝑑 ) will be around 10−1 − 10−2 𝑚⁄𝑠. On the other hand, the order of magnitude of reaction rate constant (𝑘𝑟 ) lie in the range of 10−8-10−11 𝑚/𝑠.23,25 Evidently,

1 𝑘𝑑

is practically negligible

when compared to other term in the parenthesis of the right hand side of equation (3). In this case, the mass transfer driving force in the bulk of liquid phase tends to zero, which practically 𝑝𝑎𝑟𝑡𝑖𝑐𝑙𝑒

means (𝑥𝑏

= 𝑥𝑏𝑢𝑙𝑘 ).

Figure 2 depicts the schematic representation of the overall mass resistance which is the combination of mass transfer resistance in the bulk of liquid phase and the corresponding hydrate formation reaction resistance. Assuming pseudo steady state condition, then the mass accumulation can be ignored and mass transfer rate would be equal to the reaction consumption rate. With this assumption, the equation (s) 2 can be rearranged to the following form for explicit computation of gas concentration in the bulk of water:

𝑥𝑏𝑢𝑙𝑘 =

𝜌 𝑘𝐿 𝑎 𝐶𝑤0 𝑉𝐿 𝑥𝑖𝑛𝑡 +𝑘𝑟 𝐴𝑝 ( 𝐻⁄𝑀𝑤 ) 𝐻

𝑘𝐿 𝑎 𝐶𝑤0 𝑉𝐿 +(

𝑘𝑟 𝐴 𝑝

𝜌 ⁄ 𝑒𝑞 )∗( 𝐻⁄𝑀𝑤𝐻 ) 𝑥𝑏

(5)

Measurement of real 𝑘𝐿 𝑎 value is difficult in the presence of hydrate formation process. The corresponding value measured in the solubility experiments is usually used instead. 7 ACS Paragon Plus Environment

Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 41

2.2. Population Balance Equation A proper method should be used in kinetic modeling of the hydrate process to predict the area of the hydrate particles at each instance. The population balance approach which combines a set of equations describing the conservation laws, rate equations and fluid mechanics of the system is very popular to model the nucleation and growth processes in particulate systems such as crystallization process.31 The first and only kinetic model which uses population balance method to simulate the hydrate formation process is presented by Englezos et al.12 They preferred to introduce and use various momentum concepts instead of direct solution of population balance equation for pure methane and ethane hydrate formations. In this work, direct solution of PBE will be used to predict the incremental changes of hydrate particle sizes with respect to time during hydrate formation of methane at constant pressure and temperature. Garside et al. presented the following general equation to compute the particle size distribution (𝑓: 𝑁𝑢𝑚 𝑜𝑓 𝑝𝑎𝑟𝑡𝑖𝑐𝑙𝑒/𝑚4 ) of the crystal particles for a unit volume of particulate solution.31 𝜕𝑓(𝑟.𝑡) 𝜕𝑡

+

𝜕[𝐺(𝑟.𝑡)𝑓(𝑟.𝑡)] 𝜕𝑟

= 𝐵𝑖𝑟𝑡ℎ − 𝐷𝑒𝑎𝑡ℎ

(6)

Where, 𝐺(𝑟. 𝑡) is the continuous growth function or the rate of growth during crystallization or hydrate formation process with units of velocity (i.e. m/s). Birth represents new particles formed through nucleation, agglomeration or breakage and Death represents destruction of particles through agglomeration or breakage. The birth and death terms are both discrete functions, but are often treated as continuous functions.32 These terms are neglected in this work during the growth process. The rate of particle growth during the hydrate formation process can be simplified as the following relation. 𝐺 = 𝑘𝑟 (𝑆 − 1) 𝑔

(7)

In the absence of secondary nucleation, the primary nucleation is assumed to be responsible for the entire generation (birth) process which is computable from the following equation: 𝐵0 = 𝑘𝑏 (𝑆 − 1)𝑏

(8)

In the above equations, S represents the super-saturation ratio which is defined as (𝑥𝑏𝑢𝑙𝑘 /𝑥𝑏𝑒𝑞 ). Growth and birth rate constants (kr & kb) and the corresponding exponents should be found by fitting the proposed model to the available experimental measurements. This issue will receive more attention in section 4. Evidently, both B0 and kb have the units of [number/(m3s)]. For 8 ACS Paragon Plus Environment

Page 9 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

the sake of simplicity, the value of the g exponent is assumed to be unity. The following equations describe both the initial and boundary conditions required for the equation (6). 𝑓(0. 𝑟) = 0

(𝑓𝑜𝑟 𝑎𝑛𝑦 𝑟 > 0)

𝑓(𝑡. 𝑟𝑚𝑖𝑛 ) = 𝐵0 /𝐺

(9)

(𝑓𝑜𝑟 𝑎𝑛𝑦 𝑡 > 0)

(10)

Where, rmin is the size of initial hydrate particle generated at the primary nucleation process. The above set of equations 6 to 10 should be numerically solved to find the discrete versions of both continuous functions f and G. To enhance the solution process, a maximum anticipated size (rmax) can be introduced and an additional constraint of 𝑓(𝑡. 𝑟𝑚𝑎𝑥 ) = 0 may be applied. After evaluation of PSD, the area of the entire hydrate particles required for the estimation of 𝑥𝑏𝑢𝑙𝑘 from equation (5), can be easily computed from the subsequent relation.33 ∞

𝐴𝑝 = 𝑉𝑙𝑖𝑞 ∫0 4𝜋𝑟 2 𝑓(𝑟. 𝑡) 𝑑𝑟

(11)

2.3. Numerical solution algorithm As previously mentioned, the combination of equations 2 to 11 should be solved after the dissolution process to predict the gas consumption rate during the nucleation and growth processes of the hydrate formation phenomenon at the prevailing temperature and pressure. The following steps represent the computation algorithm for performing our in-house method which involves a relatively complex task. 1. Compute minimum concentration of gas in the bulk of liquid phase (𝑥𝑏𝑒𝑞 ) at the system pressure and temperature which is necessary for hydrate formation. Evidently at this point the fugacity of the gas molecules in the bulk of liquid phase should be equal to its fugacity in the 𝑙 𝑠 hydrate structure (i.e. 𝑓𝐶𝐻 = 𝑓𝐶𝐻 ). The combination of experimental data provided by Yanga 4 4

et al., Wanjun et al. 34,35 and our fully optimized Kernel type artificial neural network known as Regularization network is used in this work to compute best value of 𝑥𝑏𝑒𝑞 for dissolution of methane gas in pure water at the prevailing constant temperature and pressure. 2. Equation (1) should be solved iteratively for the entire dissolution process to compute 𝑘𝐿 𝑎 value which fits most properly the experimental dissolution rate data for pre-specified values of initial water concentration (𝐶𝑤0 ), the liquid phase volume (𝑉𝐿 ), initial bulk concentration of the solute (i.e 𝑥𝑏𝑢𝑙𝑘 = 0) and corresponding interface concentration at the prevailing temperature and pressure. 3. Using the 𝑘𝐿 𝑎 value obtained from previous step, the gas consumption rate and the corresponding concentration in the bulk of liquid phase should be updated until the latter 9 ACS Paragon Plus Environment

Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

concentration (𝑥𝑏𝑢𝑙𝑘 ) reaches the 𝑥𝑏𝑒𝑞 computed from first step. At this point, the dissolution process terminates and the nucleation process initiates. 4. At the start of the nucleation initial hydrate cages will appear in the liquid phase. To simulate the combination of nucleation and growth processes, the PBE (eqn. 6) should be solved simultaneously with the corresponding boundary conditions and its constraint of 𝑓(𝑡. 𝑟𝑚𝑎𝑥 ) = 0 to determine particle size distribution ( f ) and the matching total surface area of hydrate particles (Ap) using equation (11). As previously mentioned, both growth and birth rate constants (kr & kb) and the corresponding exponents should be found by fitting the proposed model to the available experimental measurements. The detailed description of the recruited solution technique for solving PBE via Lax-Wendroff / Crank-Nicholson method will be thoroughly discussed in section 2.3.1. 5. Assuming pseudo-steady state condition for the entire hydrate formation process, then no accumulation exists and the mass transfer rate from interface to the bulk of liquid should be equal to the consumption rate of dissolved gas in the hydrate formation (nucleation and growth) process. In such case, the new value of 𝑥𝑏𝑢𝑙𝑘 can be readily updated via equation (5). Then equation 2 would be used to compute the rate of gas consumption during both nucleation and growth processes. 6. Since for each mole of methane or ethane molecules, 𝑛𝑤 moles of water molecules are required to construct the corresponding hydrate cages in the form of (CH4 & nwH2O or C2H6 & nwH2O), therefore by multiplying the gas consumption rate computed in the previous step by 𝑛𝑤 , the total moles of consumed water will be computed and finally the remaining moles or volume of the water in the liquid phase can be calculated. The value of 𝑛𝑤 is reported to vary between 5.75 and 6.3 for methane hydrate and 7.67 for ethane. 12,36 The entire computation algorithm described in the above steps is condensed and schematically is shown in Figure 3. To compute the particle size distribution (f) from PBE, the combined Lax-Wendroff/Crank-Nicholson method is used in this work, as described in the next section. 2.3.1 The combined Lax-Wendroff/Crank-Nicholson method This is a highly efficient numerical algorithm of a finite-element type with time discretization that leads to a stationary partial differential equation. These integral terms are implemented for every time increment in the entire nucleation and growth processes.37 Assuming extremely small value for the time increment (t), the most recent value of the growth rate (G) from the last time step can be used and equation (6) reduces as follows in the absence of death.


Page 10 of 41

Page 11 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

𝜕𝑓

𝑟

𝜕𝑓

+ 𝐺 𝜕𝑟 = ∫0 𝑚𝑎𝑥 𝐵0 𝛿(𝑟 − 𝑟𝑚𝑖𝑛 )𝑑𝑟 𝜕𝑡

(12)

where, 𝛿 is the Dirac delta function and the entire right hand side reduces to B0|at r=rmin. The Lax-Wendroff assumes a quadratic approximation of f with respect to time using the second order Taylor expansion: 𝜕2 𝑓 ∆𝑡 2

𝜕𝑓

𝑓𝑙𝑖+1 = 𝑓𝑙𝑖 + 𝜕𝑡 ∆𝑡 + 𝜕𝑡 2

(13)

2

Replacing the above approximation for f into simplified PBE, making proper arrangement and replacing from equation 12, the following equation evolves: 𝜕2 𝑓 𝜕𝑡 2

𝜕

𝜕𝑓

+ 𝐺 𝜕𝑡 (𝜕𝑟 ) = 0

𝜕2 𝑓

→

𝜕𝑡 2

𝜕2 𝑓

= 𝐺 2 𝜕𝑟 2

(14)

Now, the semi-implicit Crank-Nicholson discretization technique can be used for first and second order differential with respect to the spatial dimension (r), as the following relations: 𝜕𝑓 𝜕𝑟

𝜕2 𝑓 𝜕𝑟 2

1

𝑖+1 𝑖+1 − 𝑓𝑙−1 𝑓𝑙+1

2

2∆𝑟

= (

+

𝑖 𝑖 − 𝑓𝑙−1 𝑓𝑙+1

2∆𝑟

1

𝑖+1 𝑖+1 𝑖 𝑖 ) = 4∆𝑟 ( 𝑓𝑙+1 − 𝑓𝑙−1 + 𝑓𝑙+1 − 𝑓𝑙−1 )

1

𝑖+1 𝑖+1 𝑖 𝑖 = 2(∆𝑟)2 ( 𝑓𝑙+1 − 2 𝑓𝑙𝑖+1 + 𝑓𝑙−1 + 𝑓𝑙+1 − 2 𝑓𝑙𝑖 + 𝑓𝑙−1 )

(15) (16)

Finally, substituting equations 12,14,15 and 16 into equation (13) and performing some mathematical manipulations, leads to: 𝑖+1 𝑖+1 (−𝐴 − 𝐵) 𝑓𝑙+1 + (1 + 2𝐴) 𝑓𝑙𝑖+1 + (𝐵 − 𝐴) 𝑓𝑙−1 =

i=1,2,…,N ,

𝑖 𝑖 𝐵0 𝑡 + (𝐴 + 𝐵) 𝑓𝑙+1 + (1 − 2𝐴) 𝑓𝑙𝑖 + (𝐴 − 𝐵) 𝑓𝑙−1

Where 𝐴 =

𝐺 2 ∆𝑡 2 4∆𝑟 2

l=1,2,…,M (17)

𝐺∆𝑡

and = − 4∆𝑟 .

Equation (17) can be represented in the matrix for as φ= ϕ, where  is the unknown vector 𝑇 for any desired time interval (i.e. 𝑓 = [𝑓2 . 𝑓3 … . 𝑓𝑀−1 ]1×𝑀−2 ), ϕ is the right hand side column

vector of known values in the form of: (𝑨 − 𝑩)𝒇𝟏 𝟐 + (𝑨 + 𝑩)𝒇𝟑 𝟏 + (𝟏 − 𝟐𝑨)𝒇𝟐 𝟏 + (𝑨 − 𝑩)𝒇𝟏 𝟏 (𝑨 + 𝑩)𝒇𝟒 𝟏 + (𝟏 − 𝟐𝑨)𝒇𝟑 𝟏 + (𝑨 − 𝑩)𝒇𝟐 𝟏 𝛟= (𝑨 + 𝑩)𝒇𝟓 𝟏 + (𝟏 − 𝟐𝑨)𝒇𝟒 𝟏 + (𝑨 − 𝑩)𝒇𝟑 𝟏 ⋮ [(𝑨 − 𝑩)𝒇𝑴 𝟐 + (𝑨 + 𝑩)𝒇𝑴 𝟏 + (𝟏 − 𝟐𝑨)𝒇𝑴−𝟏 𝟏 + (𝑨 − 𝑩)𝒇𝑴−𝟐 𝟏 ]

and φ is an M-2M-2 tri-diagonal square matrix with the following form:


Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

φ=

To solve the above set of linear equations, the square matrix φ should be inverted. For small or moderate values of M, the inv or pinv command of MATLAB software can be used to compute the corresponding inverse. The computational demand for performing this task is O(M3) order of magnitude.38 To avoid excessive mathematical operations and appreciably reduce the execution time from O(M 3) to O(M ) for excessively large matrices, Thomas method may be used to take into account the special structure of the three band matrix .39 3. Case Study description The experimental data of Englezos et al. for methane hydrate formation at 276 K and constant pressures of 4.86 MPa & 7.09 MPa are used here to validate our proposed method.12 In these isothermal and isobaric tests, the hydrate formation kinetics were observed by contacting methane gas with water in the reactor chamber at constant temperature and pressure. Sufficient gas is constantly injected into hydrate formation chamber (which initially was filled with 300 cc de-ionized water) to maintain constant pressure. The consumption rate (shown in Figure 4) can be readily computed from the variations of the injected gas rate fed to the reactor over the time. The methane solid-liquid equilibria (𝑥𝑏𝑒𝑞 ) is amongst the most important parameters required in the kinetic modeling of gas hydrate formation process. As mentioned in the previous sections, at the start of the nucleation process, the equality of methane fugacity in both liquid 𝑙 𝑠 and solid phases should be satisfied (𝑓𝑢𝑔𝐶𝐻 = 𝑓𝑢𝑔𝐶𝐻 ). Instead of the thermodynamic 4 4

approach, a specially trained radial basis artificial neural network (known as Regularization Network (RN)) is used in the present study to predict the equilibrium methane concentration in the bulk of liquid phase (𝑥𝑏𝑒𝑞 ) at any prevailing temperature and pressure.40 The tabulated experimental data borrowed from the available literature is used to train the optimal RN (Table 1).


Page 12 of 41

Page 13 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Figure 5 illustrates the 3D generalization performance of the optimally trained Regularization network for the variations of 𝑥𝑏𝑒𝑞 versus the operating temperature and pressure. The cubic plus association equation of state (CPA-EOS) is also used to compute the methane equilibrium concentration at the interface of gas-liquid (xint). The CPA-EOS which is initially presented by Kontogeorgis et al. is actually a hybrid method that combines Soave–Redlich– Kwong EOS with the Wertheim association theory.41 The pressure dependency via CPA for multicomponent mixtures is usually represented as following form:42 𝑃=𝑉

𝑅𝑇

𝑚 −𝑏

−𝑉

𝑎(𝑇)

1 𝑅𝑇

𝑚 (𝑉𝑚 +𝑏)

− 2 𝑉 (1 + 𝜌

𝜕ln (𝑔)

𝑚

𝜕𝜌

) ∑𝑖 𝑥𝑖 ∑𝐴𝑖(1 − 𝑋𝐴𝑖 )

(18)

Where 𝜌 (= 1/𝑉𝑚 ) is the total molar density of the multicomponent mixture, 𝑋𝐴𝑖 represents the fraction of sites A on molecule i which do not form bonds with other active sites and 𝑥𝑖 is the mole fraction of component i in the interface of gas-liquid phases. 𝑋𝐴𝑖 is related to the so called association strength (∆𝐴𝑖𝐵𝑗 ) between two sites belonging to two different molecules (e.g. site A on molecule i and site B on molecule j) and should be computed from the following relation: 𝑋𝐴𝑖 =

1 𝐴 𝐵 ∑ ∑ 1+𝜌 𝑗 𝑥𝑗 𝐵𝑗 𝑋𝐵𝑗 ∆ 𝑖 𝑗

(19)

The association strength ∆ 𝐴𝑖𝐵𝑗 is usually expressed as: 𝐴 𝐵 𝜀 𝑖 𝑗

∆ 𝐴𝑖 𝐵𝑗 = 𝑔(𝜌) [exp (

𝑅𝑇

) − 1] 𝑏𝑖𝑗 𝛽 𝐴𝑖 𝐵𝑗

(20)

Where, parameters ε, β and bij's are tabulated in Table 2 and the radial distribution function (𝑔) is described as: 1

𝑔(𝜌) = 1−1.9𝜂

1

, 𝜂 = 4 𝑏𝜌

, 𝑏𝑖𝑗 =

𝑏𝑖 +𝑏𝑗 2

(21)

Finally, the energy parameter required in equation (18) is often given by a Soave-type temperature dependency: 𝑎(𝑇) = 𝑎0 (1 + 𝑐1 (1 − √𝑇𝑟 ))2

(22)

The combination equation (18) and tradition equilibrium criteria (𝑃𝑦𝑖 𝜑𝑖 𝑣 = 𝑃𝑥𝑖 𝜑𝑖 𝑙 ) will provide the gas equilibrium concentration at the interface of gas-liquid (xint). 4. Simulation results & discussion The population base model described via equations 6 to 11 and schematically depicted via the flowchart represented in Figure 3 is solved by resorting to the combinatory method of LaxWendroff/Crank-Nicholson technique to provide the hydrate formation rate for the 13 ACS Paragon Plus Environment

Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 41

experimental case studies borrowed from Englezos et al. 12 The kinetic rate constants shown in Table 3 are used in the entire simulations to predict the nucleation and growth rates via equations 7 & 8. As described in more details in previous works24, the experimental data extracted from Figure 4 are used along with a trial error procedure to minimize the sum of squares of the prediction errors and estimate the appropriate values of kr and kb, for the given data set. Furthermore, the reaction rate constant (𝑘𝑟 ) usually depends on various operating conditions, such as temperature, pressure and agitator speed. In practice and in the absence of reliable experimental data, the effects of the latter parameters are almost always neglected and an Arrhenius type dependency is considered for the temperature.23 As it was clearly mentioned in the beginning of section 3, the current work investigates the hydrate formation kinetics at constant temperature for two different pressures. For this reason, the temperature dependency of 𝑘𝑟 is not addressed throughout this article To ensure proper accuracy and robustness of the simulated results, the proposed model predictions should be initially verified and finally validated with experimental data. 4.1 Verification For verification purpose, both radial step size (dr) and the corresponding time step (dt) should be selected in a manner that the final solution of the population balance equation (PBE) would be independent of those step lengths. The PBE should be integrated between minimum radius of the hydrate formed (r min) at the nucleation process and a pre-specified maximum radius. To be on a safe side, the latter radius is assumed to be 1200 nm during the hydrate formation process. Assuming that initially generated hydrate particles are solid spheres containing two cages of 512 and six cages of 51262, then the value of rmin can be simply computed from Table 4 and the following volumetric average relation: 2

6

3

2

6

𝑉𝑚𝑒𝑎𝑛 = (8) 𝑉𝑠𝑚𝑎𝑙𝑙 𝑐𝑎𝑔𝑒 + (8) 𝑉𝑙𝑎𝑟𝑔𝑒 𝑐𝑎𝑔𝑒 → 𝑟𝑚𝑖𝑛 = √(8) (0.395)3 + (8) (0.433)3 = 0.4241 𝑛𝑚

Figure 6 compares four different particle size distributions (f's) computed from the solution of PBE for the overall hydrate formation time of 10 minutes after the dissolution process. In all cases, the radial step size is kept fixed at dr=0.8nm during the solution procedure but different time steps (dt's) of 10s, 1s, 0.1s and 0.05s are used. As can be seen, the predicted distribution dramatically changes when the time step varies from 10 second to 0.1s. On the


Page 15 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

other hand, the final solution (distribution) is practically independent of time step when the corresponding step size is smaller than 0.1s for dr=0.8nm. Obviously, the value of dt=0.1s is used in all next simulations to minimize the required execution time. Evidently, both excessive oscillations and negative values encountered in particle size distributions computed with the time steps of 1 and 10 seconds are due to lack of convergence and instability of the computed final solutions. To verify the proper value of radial step length, a similar procedure is repeated by keeping the time step fixed at 0.1s and using different radial step sizes of 5, 3, 0.8 and 0.5 nanometers. Figure 7 clearly shows that the predicted distribution exceedingly changes when dr varies from 5nm to 0.8nm. Once again, the ultimate distribution is almost the same when radial step size is less than 0.8nm. All the next simulation results will be carried out by assuming dr=0.8nm and dt=0.1s. 4.2 Simulation predictions and final validations Figures 8a & 8b provide various plots of particle size distributions computed at four different times after the start of the nucleation process for two distinct pressures of 4.86 and 7.9 MPa, respectively. As can be seen, the maximum size of hydrate particles produced at the lower pressure of 4.86MPa, is less than 800nm after 7400s, while the corresponding value rapidly exceeds 900nm after only 2200s. This is due to higher dissolution rates at larger gas pressures which will lead to higher driving forces inside the liquid phase and provide higher growth rates. To the best of our knowledge, no previous research has provided PSD plots for the hydrate formation process, therefore the simulation results of Figure 8 can't be validated directly. Fortunately, the experimental results reported in the literature indicate that maximum size of hydrate particles obtained in the laboratories are around 1000 nm after reasonable time passed from the start of nucleation process.23 As Figure 9a depicts, the hydrate particles surface area increases monotonically with time for both pressures due to continuous growth of the corresponding clathrates. For higher pressures, both nucleation and growth processes occur more intensely and rapidly due to higher methane concentrations in the liquid phase. Figure 9b confirms the above statement and shows that the primary nucleation rate starts with an intense peak and then reduces rapidly due to its positive Gibb's free energy compared to the growth process which has negative free energy as depicted in Figure 10. The initial primary nucleation rates (𝐵0𝑖 ) for both pressures are computed as 2.8 × 1014 and 6 × 1016 (𝑁𝑢𝑚⁄ 3 ), respectively which is completely in accordance with our 𝑚 𝑠

previous discussion.


Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The final primary nucleation rates (𝐵0𝑓 ) are estimated as 2∙ 9 × 108 and 0∙ 3(𝑁𝑢𝑚⁄𝑚3 𝑠) at the end of simulation time for the pressures of 4.86 and 7.09 MPa, respectively. As can be seen in Figure 10, the primary hydrate nucleation rate drastically drops for any pressure, due to exponential increase in Gibb's free energy of the primary nucleation process. This effect is more pronounce for larger pressures because at sufficiently high pressures, the number of previously formed hydrate particles are extremely higher (around 100 times) and less nucleation sites are required. In other words, the growth process (with negative Gibb's free energy, as depicted in Figure 10) would be more dominant for larger pressures as the hydrate formation time increases. Figure 10 clearly illustrates that the overall Gibbs free energy (weighted sum of both nucleation and growth processes) experiences a distinct maximum at critical particle radius (rcr) and afterwards drops monotonically and continuously. Evidently, the non-spontaneous nucleation phenomenon (ΔG> 0) becomes less probable as time goes by, especially for larger pressures. Figures 11 show our model predictions for the variations of methane concentration in the bulk liquid phase (xbulk) during the entire stages of gas hydrate formation process at two pressures of P=4.86 and 7.9 MPa. In both cases, the value of xbulk increases monotonically during the dissolution stages due to the positive accumulation of methane molecules in the water, as previously discussed. The first nuclei spontaneously form when the fugacity of the methane at the bulk of liquid phase reaches the minimum fugacity required for the initial hydrate generation. The negative accumulation of methane molecules inside water is due to the methane consumption via the nucleation process and tends to decrease the value of xbulk. However, each hydrate particle requires around 6 water molecules per each methane molecule and therefore leads to a sudden increase in xbulk value. This increase is larger for higher pressures, because the nucleation rate increases more rapidly with pressure. Obviously, the concentration of methane in the bulk of liquid phase can't exceed the corresponding interface concentration at the prevailing temperature and pressure. This observation can be also explained from another point of view. At early stages of hydrate formation, the surface area of solid particles (𝐴𝑝 ) are extremely small and therefore, the rate of mass transfer in the bulk of liquid phase exceeds the methane consumption due to both nucleation and growth reactions. This process tends to increase the value of xbulk. As time goes by, the value of 𝐴𝑝 grows asymptotically and the corresponding bulk concentration tends to decrease to its equilibrium value (xbeq).


Page 16 of 41

Page 17 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Figures 12 compare our simulation results with the several experimental data of Englezos et al. at 276K regarding the variations of consumed methane moles versus time at two working pressures of 4.86 MPa and 7.09 MPa. In both cases, the value of growth exponent (g) in equation (7) is set to unity. The simulation results show close agreement with experimental data collected at lower pressure (P=4.86MPa), but perform impressively good for the higher pressure of 7.09MPa. To improve the model performance for the former case, the value of g is increased and Figure 13 illustrates excellent agreement between the simulation results and experimental data for g=2. Any remaining discrepancies between the simulation results and experimental measurements may be due to ignoring secondary nucleation, breakage and death terms. The following relation should be used instead of equation 5 to compute the methane mole fraction in the bulk of liquid phase (xbulk), when the value of g is equal to 2. 𝑓(𝑥𝑏𝑢𝑙𝑘 ) = 𝑘𝑟 𝐴𝑝 (

𝑒𝑞 2 𝜌𝐻𝑦𝑑 𝑥𝑏𝑢𝑙𝑘 −𝑥 ⁄𝑀𝑤 ) ( 𝑥 𝑒𝑞 𝑏 ) − 𝑘𝐿 𝑎 𝐶𝑤0 𝑉𝐿 (𝑥𝑖𝑛𝑡 − 𝑥𝑏𝑢𝑙𝑘 ) = 0 𝐻𝑦𝑑 𝑏

(23)

To investigate that our proposed model can be actually extended to other cases, two experimental data are borrowed from Clarke and Bishnoi for the ethane and carbon dioxide hydrate formation processes.27,28 Figure 14 clearly illustrates the successful performance of the proposed model on prediction the entire three steps ethane hydrate formation process. This was anticipated, because the dissolution process for both methane and ethane constituents are exactly the same. The computed values for birth and growth rate constants (kb & kr) and corresponding exponents (b & g) are calculates as 31016 number/(m3.s), 510-8 m/s, 4.0 and 1.0, respectively. Our further investigations revealed that the proposed model actually failed to successfully reproduce the carbon dioxide hydrate formation data as presented by Clarke and Bishnoi (results not included). This was not also very surprising, because in previous cases, the dissolution of both methane and ethane species were molecular type and no ionization was happened. In contrast, the dissolution of carbon dioxide in water leads to partial ionization and some of the CO2 is practically converts to carbonate ion. Evidently, this fourth step is not included in this model. In other words, our new model performs very satisfactorily for light hydrocarbons hydrate formation processes. Figure 15 illustrates large differences between the long term simulation results obtained by our model and those predicted using the Englezos et al. model for hydrate formation at constant pressure of 7.09MPa. As can be seen, the model presented by Englezos et al. always over17 ACS Paragon Plus Environment

Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

predicts the methane consumption rates, especially at very long times. As Figure 15 depicts and Skovborg and Rasmussen have also previously pointed out, in all kinetically oriented models (such as Englezos et al.) the overall hydrate formation rate is directly proportional to the surface area of the entire hydrate particles and the methane consumption rates predicted via such models show a monotonically increasing trend with growing slope. In our more complex model which combines the hydrate formation reaction kinetics with mass transfer resistance inside the liquid phase, the trend is also monotonically increasing but as shown in Figure 15, the growth rate is continuously decreasing due to increase in mass transfer resistance inside the liquid phase during excessively long term hydrate formation processes. To be more precise, Figure 15 clearly shows that according to Englezos et al. model predictions, the amount of methane consumed should exceed 3 moles after 16 hours. Regarding the previously mentioned structure for the hydrate (CH4,6H2O), the approximate quantity of water required to construct these cages should be around 18 moles, while, the entire water initially charged to the vessel was around 17 grmol (300 cm3). Evidently, the predicted rate via Englezos et al. model contradicts with the reality and could not be true. 5. Conclusion A new model was presented which combines the mass transfer in liquid phase with hydrate nucleation and growth processes. The population balance theory was applied to the entire hydrate formation process to estimate the final rate of hydrate production and corresponding total surface area at constant pressure operations. The combination of Lax-Wendroff/CrankNicholson method was used to solve the partial differential equation derived from this approach. The proposed model was initially verified and then successfully validated by resorting to the experimental data of Englezos et al. The predictions of our model for the final size of the hydrate particles at the prevailing temperature and pressure (around 800 nanometers) were in close agreement with the previously reported values measured in several empirical studies. A steep jump in methane concentration in the bulk of liquid phase (xbulk) was encountered at the end of dissolution process to force the nonspontaneous nucleation process. Afterwards, the spontaneous growth process with negative Gibb's free energy accelerated the methane consumption and lead to gradual decrease in xbulk value. According to our simulation results, the reaction resistance due to the nucleation process dominates the mass transfer resistance at early stage of hydrate formation process. As the growth process evolves, the mass transfer


Page 18 of 41

Page 19 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

controls the overall rate of clathrate production. The proposed model was also successfully validated by experimental data of Clarke and Bishnoi for ethane hydrate formation process.

AUTHOR INFORMATION Corresponding Author: *

E-mail: [email protected]

Notes The authors have no conflicts of interest to declare.


Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

References [1] Linga, P., Haligva, C., Nam, S.C., Ripmeester, J.A., Englezos, P. Gas hydrate formation in a variable volume bed of silica sand particles. Energy Fuels, 2009, 23, 5496-5507. [2] Hester, K.C., Brewer, P.G. Clathrate hydrates in nature. Annu. Rev. Mar. Sci., 2009, 1, 30327. [3] Sun, C., Li, W., Yang, X., Li, F., Yuan, Q., Mu, L. Progress in research of gas hydrate. Chin J. Chem Eng., 2011, 19, 151-162. [4] Kim, N. J., Hwan Lee, J., Cho, Y. S., Chun, W. Formation enhancement of methane hydrate for natural gas transport and storage. Energy, 2010, 35, 2717-2722. [5] Heeschen, K.U., Abendroth, S., Priegnitz, M., Spangenberg, E., Thaler, J., Schicks, J. M. Gas Production from Methane Hydrate: A Laboratory Simulation of the Multistage Depressurization Test in Mallik, Northwest Territories, Canada. Energy Fuels, 2016, 30, 62106219. [6] Boswell, R., Schoderbek, D., Collett, T.S., Ohtsuki, S., White, M., Anderson, B.J. The Iġnik Sikumi Field Experiment, Alaska North Slope: Design, Operations, and Implications for CO2-CH4 Exchange in Gas Hydrate Reservoirs. Energy Fuels, 2017, 31, 140-153. [7] Ogawa, T., Ito, T., Watanabe, K., Tahara, K., Hiraoka, R., Ochiai, J., Ohmura, R., Mori, Y. H. Development of a novel hydrate-based refrigeration system: A preliminary overview. Applied Thermal Engineering, 2006, 26, 2157-2167. [8] Lorenzo, M. D., Aman, Z. M., Kozielski, K., Norris, B. W. E., Johns, M. L., May, E. F. Underinhibited hydrate formation and transport investigated using a single-pass gas-dominant flowloop. Energy Fuels, 2014, 28, 7274-7284. [9] Wang, Z. Y., Zhao, Y., Sun, B. J., Chen, L. T., Zhang, J. B., Wang, X. R. Modeling of hydrate blockage in gas-dominated systems. Energy Fuels, 2016, 30, 4653-4666. [10] Ribeiro Jr., C.P, Lage, P.L. Modelling of hydrate formation kinetics: State-of-the-art and future directions. Chemical Engineering Science, 2008, 63, 2007-2034. [11] Vysniauskas, A., Bishnoi, P. Kinetic study of methane hydrate formation. Chemical Engineering Science, 1983, 38, 1061-1072. [12] Englezos, P., Kalogerakis, N., Dholabhai, P.D, Bishnoi, P.R. Kinetics of formation of methane and ethane gas hydrates. Chemical Engineering Science, 1987a, 42, 2647-2658. [13] Englezos, P., Kalogerakis, N.E., Dholabhai, P.D., Bishnoi, P.R. Kinetics of gas hydrate formation from mixtures of methane and ethane. Chemical Engineering Science, 1987b, 42, 2659-2666. 20 ACS Paragon Plus Environment

Page 20 of 41

Page 21 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

[14] Skovborg, P., Rasmussen, P. A mass transport limited model for the growth of methane and ethane gas hydrates. Chemical Engineering Science, 1994, 49, 1131-1143. [15] Herri, J.M., Pic, J.S., Gruy, F., Cournil, M. Interest of in situ turbidimetry for the characterization of methane hydrate crystallization: application to the study of kinetic inhibitors. Chemical Engineering Science, 1999a, 54, 1849-1858. [16] Herri, J.M., Pic, J.S., Gruy, F., Cournil, M. Methane hydrate crystallization mechanism from in situ particle sizing. American Institute of Chemical Engineers, 1999b, 45, 590-602. [17] Kashchiev, D., Firoozabadi, A. Driving force for crystallization of gas hydrates. Journal of Crystal Growth, 2002a, 241, 220-230. [18] Kashchiev, D., Firoozabadi, A. Nucleation of gas hydrates. Journal of Crystal Growth, 2002b, 243, 476-489. [19] Kashchiev, D., Firoozabadi, A. Induction time in crystallization of gas hydrates. Journal of Crystal Growth, 2003, 250, 499-515. [20] Bergeron, S., Servio, P. Reaction rate constant of propane hydrate formation. Fluid Phase Equilibria, 2008, 265, 30-36. [21] Bergeron, S., Servio, P. Reaction Rate Constant of CO2 Hydrate Formation and Verification of old Premises Pertaining to Hydrate Growth Kinetics. American Institute of Chemical Engineers, 2008, 54, 2964-2970. [22] Bergeron, S., Servio, P. CO2 and CH4 mole fraction measurements during hydrate growth in a semi-batch stirred tank reactor and its significance to kinetic modeling. Fluid Phase Equilibria, 2009, 276, 150-155. [23] Bergeron, S., Servio, P., Beltrán, J. Reaction rate constant of methane clathrate formation. Fuel, 2010, 89, 294-301. [24] ZareNezhad, B., Mottahedin, M., Varaminian, F. A new approach for determination of single component gas Hydrate formation kinetics in the absence or presence of kinetic promoters. Chemical Engineering Science, 2015, 81, 28-37. [25] Herri, J.M., Kwaterski, M. Derivation of a Langmuir type of model to describe the intrinsic growth rate of gas hydrates during crystallization from gas mixtures. Chemical Engineering Science, 2012, 137, 447-457. [26] Khosharay, S., Roosta, H., Varaminian, F. Investigation on the kinetics of methane and carbon dioxide hydrates by using a modified kinetic model. Journal of Natural Gas Science and Engineering, 2015, 26, 587-594. [27] Clarke, M., Bishnoi, P. R. Determination of the intrinsic rate of ethane gas hydrate decomposition. Chemical Engineering Science, 2000, 55, 4869-4883. 21 ACS Paragon Plus Environment

Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

[28] Clarke, M., Bishnoi, P. R. Determination of the intrinsic kinetics of CO2 gas hydrate formation using in situ particle size analysis. Chemical Engineering Science, 2005, 60, 695709. [29] Mohebbi, V., Naderifar, A., Behbahani, R.M., Moshfeghian, M. Investigation of kinetics of methane hydrate formation during isobaric and isochoric processes in an agitated reactor. Chemical Engineering Science, 2012, 76, 58-65. [30]. Cussler, E. L. Diffusion: Mass Transfer in Fluid Systems. Cambridge University Press, 1997, USA. [31] Garside, J. Industrial crystallization from solution. Chemical Engineering Science, 1985, 40, 3-26. [32] Meisler, K. T. Multi-dimensional population balance models of crystallization processes. Ph.D. Thesis, 2014, DTU, Denmark. [33] Jafari, D. Experimental Investigation and Modeling of the Semi-Batch supercritical antisolvent crystallization process for the production of pharmacuetical particles. Ph.D. Thesis, 2016, Ferdowsi University, Iran. [34] Yanga, S.O., Choa, S.H., Lee, C.S., Lee, H. Measurement and prediction of phase equilibria for water + methane in hydrate forming conditions. Fluid Phase Equilibria, 2001, 185, 53-63. [35] Wanjun, L, Burruss, R., MingChou, I. Determination of methane concentrations in water in equilibrium with SI methane hydrate in the absence of a vapor phase by in situ Raman spectroscopy. Geochimica et Cosmochimica Acta, 2008, 72, 412-422. [36] Circone, S, Kirby, S.H., Stern, A.L. Direct measurement of methane hydrate composition along the hydrate equilibrium boundary. The Journal of Physical Chemistry B, 2005, 109, 9468-9475. [37] Bennett, M.K., Rohani, S. Solution of population balance equations with a new combined Lax-Wendroff / Crank-Nicholson method. Chemical Engineering Science, 2001, 56, 66236633. [38]. Hastie, T, Tibshirani, R, Buja, A. Generalized Additive Models. CRC Press, 1990, USA. [39] Higham, N.J. Accuracy and stability of numerical algorithms. Society for Industrial and Applied Mathematics, 2002, England. [40] Shahsavand, A, PourafshariChenar, M. Neural networks modeling of hollow fiber membrane processes. Journal of Membrane Science, 2007, 297, 59-73. [41] Kontogeorgis, G.M., Voutsas, E., Yakoumis, I., Tassios, D.P. An equation of state for associating fluids. Ind. Eng. Chem. Res. 1996, 35, 4310-4318. 22 ACS Paragon Plus Environment

Page 22 of 41

Page 23 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

[42] Folas, G.K., Kontogeorgis, G.M. Thermodynamic models for industrial applications: from classical and advanced mixing rules to association theories. John Wiley & Sons Ltd, 2010, USA. [43] Abolala, M., Varaminian, F. Modeling the solubility of light reservoir components, HCFCs and HFCs in water using the CPA and sPC-SAFT equations of state. Journal of Molecular Liquids, 2013, 187, 359-367.


Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

List of Tables: Table 1. Variations of methane mole fraction in the bulk of liquid phase at equilibrium with solid for different temperature and pressures.34,35

Table 2. CPA parameters for methane and water species 43

Table 3. Kinetic constants of the nucleation and growth processes.

Table 4. Information of SI cages.


Page 24 of 41

Page 25 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Table 1. Variations of methane mole fraction in the bulk of liquid phase at equilibrium with solid for different temperature and pressures. 34,35 Temperature(K)

Pressure(bar)

xbeq

273.1 273.1 273.1 273.1 273.1 273.1 273.1 273.1 278.1 278.1 278.2 278.2 278.2 278.1 278.2 278.1 278.1 278.1 278.1 276.55 279.55 284.55 276.55 279.55 284.55 289.55 276.45 279.55 284.45 289.55 294.05 276.55 279.55 284.55 289.55

49.8 52.0 78.5 84.2 116.3 122.8 135.0 148.1 57.9 81.2 88.9 104.4 111.8 115.3 137.6 160.2 165.4 192.9 193.5 100 100 100 200 200 200 200 300 300 300 300 300 400 400 400 400

0.000775 0.000751 0.000770 0.000776 0.000763 0.000752 0.000807 0.000765 0.00114 0.00103 0.00104 0.00113 0.00101 0.00106 0.00106 0.00100 0.000953 0.000954 0.000960 0.001254 0.001503 0.002043 0.001205 0.001462 0.001980 0.002689 0.001204 0.001363 0.001849 0.002550 0.003454 0.001060 0.001236 0.001728 0.002358

Table 2. CPA parameters for methane and water species 43 25 ACS Paragon Plus Environment

Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 41

Component

𝒂𝟎 (𝒎𝟑 𝒎𝒐𝒍−𝟐 )

𝒃 ∗ 𝟏𝟎𝟓 (𝒎𝟑 𝒎𝒐𝒍−𝟏 )

𝒄𝟏

𝜺(𝑱𝒎𝒐𝒍−𝟏 )

𝜷

Water Methane

0.121324 0.230621

1.45000 2.90502

0.670000 0.456099

16599.96 -

0.07093 -

Table 3. Kinetic constants of the nucleation and growth processes

Pressure (MPa)

b

kb (𝑵𝒖𝒎⁄ 𝟑 ) 𝒎 𝒔

kr (𝒎⁄𝒔)

7.09

6

4 ∗ 1016

1 ∗ 10−7

4.86

6

3 ∗ 1015

7 ∗ 10−10

Table 4. Information of SI cages Cavity

Number of cavities

Average cavity radius (nm)

Number of water molecules

512 51262

2 6

0.395 0.433

46


Page 27 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Figure 1. Various mass transfer and reaction resistances between gas-liquid & hydrate phases.29 132x93mm (96 x 96 DPI)

ACS Paragon Plus Environment

Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 2. Schematic representation of various resistances involved in the hydrate formation process. 107x19mm (96 x 96 DPI)


Page 28 of 41

Page 29 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Figure 3. Flowchart of our in-house computation algorithm. 173x235mm (96 x 96 DPI)


Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 4. Experimental data for consumption of methane gas during the entire hydrate formation process versus time at T=276 K .12 143x69mm (96 x 96 DPI)


Page 30 of 41

Page 31 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Figure 5. The 3D plot of xbeq, predicted via the optimally trained RN versus T and P. 154x68mm (96 x 96 DPI)


Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 6. Particle size distribution curves obtained from the solution of PBE by using various time steps at fixed radial step length of dr=0.8 nm for P=7.09 MPa. 178x88mm (96 x 96 DPI)


Page 32 of 41

Page 33 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Figure 7. Particle size distribution curves obtained from the solution of PBE by using various radial step lengths at fixed time step of dt=0.1s for P=7.09 MPa.

182x88mm (96 x 96 DPI)


Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 8. Particle size distributions at different times after onset of nucleation for a) P=4.86 MPa and b) P=7.09 MPa

146x130mm (96 x 96 DPI)


Page 34 of 41

Page 35 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Figure 9. Variations of a) hydrate particle surface areas and b) primary nucleation rates vs. time after the onset of nucleation for two gas pressures. 178x119mm (96 x 96 DPI)


Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 10. Gibb's free energy profile during various steps of hydrate formation process 73x79mm (96 x 96 DPI)


Page 36 of 41

Page 37 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Figure 11. Variations of methane mole fraction in the bulk of liquid phase at a) P=4.86 MPa, b) P=7.09 MPa

133x135mm (96 x 96 DPI)


Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 12. Validation of our simulation results for two pressures of a) P=4.86 MPa and b) P=7.09 MPa using experimental data at 276K and unit growth exponent (g=1). 141x141mm (96 x 96 DPI)


Page 38 of 41

Page 39 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Figure 13. Improvement of the simulation results for the lower pressure case (P=4.86MPa) by considering a quadratic form for the growth rate (g=2). 139x67mm (96 x 96 DPI)


Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 14. Successful validation of our simulation results with unit growth exponent (g=1), for the entire three steps of ethane hydrate formation process at P=0.837 MPa and T=2.6°C. 139x69mm (96 x 96 DPI)


Page 40 of 41

Page 41 of 41 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Energy & Fuels

Figure 15. Significant difference between our simulation results and those predicted by Englezos et al. model for long term hydrate formation process at T=276 K and P=7.09 MPa. 141x70mm (96 x 96 DPI)


Application of Population Balance Theory for Dynamic Modeling of

Recommend Documents