CFD modeling of coal pyrolysis in externally heated fixed-bed reactor

Fuel 233 (2018) 685–694

Contents lists available at ScienceDirect

Fuel journal homepage: www.elsevier.com/locate/fuel

Full Length Article

CFD modeling of coal pyrolysis in externally heated fixed-bed reactor Yanan Qian

a,b

c

, Yin Yu , Guangwen Xu

a,d

, Xiaoxing Liu

a,b,⁎

T

a

State Key Laboratory of Multiphase Complex Systems, Institute of Process Engineering, Chinese Academy of Sciences, Beijing 100190, China University of Chinese Academy of Sciences, Beijing 100049, China National Key Laboratory of Shock Wave and Detonation Physics, Institute of Fluid Physics, CAEP, Mianyang 621900, China d Institute of Industrial Chemistry and Energy Technology, Shenyang University of Chemical Technology, Shenyang 110142, China b c

A R T I C LE I N FO

A B S T R A C T

Keywords: Coal pyrolysis CFD Water evaporation and condensation Volatiles Heat transfer

A three dimensional transient model was developed to simulate the pyrolysis process of low-rank sub-bituminous coal in a laboratory-scale externally heated fixed-bed reactor. The water evaporation and the release of volatiles were taken into account in the model. The simulation was first validated by the central temperature evolution of reactor. The evolution of bed temperature, generation of gaseous products, variation of voidage and its influence on the flow behavior of gaseous products were further analyzed. The obtained results show that the phase change of moisture determined the heating process, delaying the temperature of inner coal zone increase above the boiling temperature and consequently the start of pyrolysis. The non-uniform radial distribution of voidage led to the phenomenon that gaseous products generated in the inner low temperature zone tended to flow radially towards heating wall before flowing out the reactor. Analyses indicate that in the temperature range of intense volatiles release the dominant heat transfer mechanism was radiation.

1. Introduction Pyrolysis for producing valuable chemicals and liquid fuels provides an efficient way to utilize coal resources with added product value. This is important especially for China which has relative abundant coal resources but lacks seriously crude oil and natural gas resources [1]. In recent years the increasing cost of petroleum and high dependence of China crude oil supply on petroleum import greatly reawakened the necessity of research and development for high efficient coal pyrolysis technologies [2]. According to heating strategy, coal pyrolysis technologies can be roughly classified into two types: direct-heating where hot solid or gas act as heat carrier [3,4], and indirect-heating where energy is supplied through the externally heated walls and the reactor is generally fixed/moving bed [5,6]. Compared with direct-heating, indirect-heating coal pyrolysis technology can produce liquid fuels with low dust intake and fuel gas free of dilution [1]. Its disadvantage is the low efficiency in heat transfer. Although in the past years extensive experimental studies on indirect-heating coal pyrolysis technology have been performed with the aim of improving the heating efficiency and regulating the flow of pyrolysis gas [1,6–8], the indirect-heating coal pyrolysis is still generally based on empirical development. The main reason is that, due to the harsh thermochemical and physicochemical conditions within the reactor, precise measurement of the temperature and products evolutions during continuous pyrolysis operation is

⁎

generally difficult, if not impossible [1,6]. With the rapid development of high-performance computer technology, Computational Fluid Dynamics (CFD) has been increasingly used as a supplementary of chemical experiments to explore the complex heat and mass transfer characteristics of the reactors. Multiphase CFD models incorporating physical and chemical reactions and interphase momentum, heat and mass transfer have been used to model different systems like biomass fast pyrolysis [9,10], coal combustion and gasification [11,12], and fluid catalytic cracking [13]. Researches focusing on coal pyrolysis in different types of reactors have also been realized in CFD simulations. Hart et al. [14] modeled the pyrolysis of a Loy Yang low-rank coal in a pressurized drop tube furnace. Based on the CFD simulation results, they discussed the influence of using single reaction and competing reaction models on the predicted mass loss. Iavarone et al. [15] evaluated the accuracies of different one-step reaction models in predicting the volatile yield of bituminous coal in an entrained flow reactor based on CFD modeling results. Shu et al. [16] modeled coal pyrolysis in a downer reactor and discussed the contributions of different heat transfer mechanisms to the temperature increase of particles. Yan et al. [17] simulated particle heating and devolatilization during rapid coal pyrolysis in a thermal plasma downer reactor using an Euler-Lagrange model. Through similar approach recently Ma et al. [18] modeled the hydropyrolysis of pulverized coal in magnetically rotating hydrogen plasma downer reactor. CFD studies

Corresponding author at: State Key Laboratory of Multiphase Complex Systems, Institute of Process Engineering, Chinese Academy of Sciences, Beijing 100190, China. E-mail address: [email protected] (X. Liu).

https://doi.org/10.1016/j.fuel.2018.06.100 Received 20 December 2017; Received in revised form 5 June 2018; Accepted 25 June 2018 0016-2361/ © 2018 Elsevier Ltd. All rights reserved.

Fuel 233 (2018) 685–694

Y. Qian et al.

W/m3 Sradi, g ,Sradi, s heat source due to gas and solid phase radiation, W/m3 Sreaction heat transfer due to coal pyrolysis, W/m3 Tg , Ts temperature of gas and solid phases, K temperature of liquid water and water vapor, K Tl , Tv Tref reference temperature, 298 K the saturation temperature of water, 373 K Tsat ⎯→ ⎯u , ⎯→ ⎯ gas and solid velocity vectors, m/s g us v yield of volatile products in a control volume v∞ potential ultimate yield of the volatile products Yg, i mass fraction of species i in gas phase, kg/kg mass fraction of species j in solid phase, kg/kg Ys, j

Nomenclature pre-exponential factor, s−1 absorption coefficient, 1/m empirical deformation parameter coefficients of condensation coefficients of evaporation specific heat, kJ/(kg·K) diffusion coefficient for gas i through a gas mixture, m2/s particle diameter, m apparent activation energy, J/mol radiation exchange factor gravity acceleration, m/s2 specific enthalpy of gas and solid phases, J/kg heat transfer coefficient between gas and solid phases, W/ (m3·K) hH2O the latent heat of vaporization, J/kg I radiation intensity, W/(m2·sr) ⎯⎯⎯→ Jg, i diffusion flux of species i in gas phase, kg/(m2·s) gas and solid thermal conductivities, W/(m·K) k g , ks ksr, e effective radiative conductivity of solid phase, W/(m·K) k g, e , ks, e effective thermal conductivity of gas and solid phases, W/ (m·K) devolatilization rate constant kv ṁ coal → product devolatilization rate of coal, kg/(m3·s) ṁ coal → species(i) generation rate of species i, kg/(m3·s) ṁ coal → volatile generation rate of total volatiles, kg/(m3·s) ṁ H2O mass transfer of moisture between phases, kg/(m3·s) ṁ reaction mass transfer flux due to coal pyrolysis, kg/(m3·s) ṁ sg mass transfer from solid phase to gas phase, kg/(m3·s) ṁ s jg i mass transfer from species j of solid phase to species i of gas phase, kg/(m3·s) ṁ s, jj′ mass flux between different species in solid phase, kg/ (m3·s) nusselt number, Nus Pr prandtl number p pressure, N/m2 Qreaction reaction heat, J/kg qradi, s radiative heat flux through solid phase, W/m2 R universal gas constant, 8.314 J/(mol·K) reynolds number Res SH2O heat transfer due to the phase transformation of moisture,

A0 a B Ccon Cev Cp Di, m dp E0 FE → g Hg , Hs hgs

Greek symbols

βgs Γ εg , εs εl , ε v εr κ μg Af ρg , ρs ρl , ρv σ τg Ω ω

gas-solid drag coefficient, kg/(m3·s) combination parameter voidage and volume fraction of solid phase volume fractions of liquid water and water vapor particle emissivity dimensionless parameter, κs / κ g dynamic viscosity, Pa·s dimensionless solid conductivity density of gas and solid phases, kg/m3 density of liquid water and water vapor, kg/m3 Stefan-Boltzmann constant, 5.67 × 10−8 W/(m2·K4) shear stress of gas phase, N/m2 solid angle, sr surface fraction parameter

Subscripts

e g l p radi ref s v

about coal pyrolysis in fixed-bed reactor are mainly come from the research region of coking. And research attentions are generally paid to the influence of formed plastic layer on the flow of gaseous products. Guo and Tang [19] investigated the temperature evolution, gas flow and voidage variation of solid phase during the coking process using PHOENICS CFD package. Polesek-Karczewska et al. [20] built a onedimensional three phase transient model for wet coal pyrolysis and discussed the effect of low permeability of coal in its plastic stage on the internal pressure peaks. Lin et al. [21] developed a two dimensional transient model of coal carbonization and modeled the heat transfer and fluid flow structure in a 6-m-height coke oven. Slupik et al. [22] modeled the carbonization process in a 1.2-m-height coke oven. Their simulation results indicate that the coke shrinkage nearly has no influence on heat transfer. Based on our previous studies [23] on the evaluations of different conduction and radiation heat transfer models, in this work the pyrolysis process of low-rank sub-bituminous coal in a laboratory-scale externally heated cylindrical reactor was modeled by taking the water evaporation and volatiles release into account. Our previous experimental results [1,7] show that the final yield and quality of the tar product can be notably improved by mounting specially designed

effective gas phase liquid water particle radiation reference temperature solid phase vapor

internals in the reactor. For the future optimization and scale-up of the new developed pyrolyzer, it is thus necessary to gain a thorough understanding and prediction of the heat transfer and reaction characteristics inside the reactor through CFD simulation. As a first step, in this work we will mainly test the predictive ability of our three-dimensional transient CFD model. For this purpose, the externally heated reactor without internals (Type A reactor in [7]) was modeled. Rather than treating the coal charge as porous medium [19–22], in our simulations the coal charge was approximated to continuous fluid. The temperature evolution of coal charge was simulated by solving the energy conservation equation of solids phase. The porosity variation was automatically considered by modeling the water evaporation and volatiles release. And the contributions of heat conduction and radiation to the effective thermal conductivity of coal bed were separately considered through the two parameters of effective thermal conductivity and radiative conductivity. The simulation was validated based on the temperature evolution experimentally measured in the center of reactor.

686

Fuel 233 (2018) 685–694

Y. Qian et al.

∂ ⎯ + ∇ ·(ε k ∇T )−h (T −T ) ⎯ H ) = ε ∂p + τ : ∇ ⎯→ (εs ρs Hs ) + ∇ ·(εs ρs ⎯→ u u s s s g s s s, e s gs g s ∂t ∂t

2. Mathematical model and simulation layout The applied mathematical model was based on the following assumptions:

+ Sradi, s−SH2O + Sreaction

(1) The gaseous mixture was treated as an ideal gas. (2) Moisture diffusion inside coal particles was neglected. (3) The following contents of gaseous pyrolysis products was assumed: H2, CO, CO2, CH4, C2H6, H2O and tar (4) The temperature of reactor wall was assumed to be the same of that of heating furnace. (5) Size variation of coal particle during pyrolysis was ignored. 2.1. Model development A two-fluid model was adopted to describe the mass and heat transfer process inside a laboratory-scale externally heated fixed-bed composed of coal particles and pyrolysis gas. The conservation equations of continuity, momentum, energy, and species were integrated with the conduction-radiation model, pyrolysis kinetics model, and moisture evaporation model to simulate the heat and mass transfer behaviors using Fluent 6.3.26.

Gas phase

Solid phase

∂ ⎯ ) = −ṁ u (εs ρs ) + ∇·(εs ρs ⎯→ s sg ∂t

(1)

(3)

⎯ − ⎯→ ⎯u ) ṁ sg ( ⎯→ u s g represents the momentum transfer caused by mass transfer between phases. βgs denotes the gas-solid drag coefficient due to momentum exchange between gas phase and solid phase. In this work, the Gidaspow drag correlation [24] which is one of the default drag models available in Fluent 6.3.26 was used to calculate the inter-phase drag interaction. For the investigated fixed bed system, solid phase was stationary. Thus the momentum conservation equation for the solid phase was ⎯ =→ u 0. ignored and the velocity of solid was set to be zero, ⎯→ s The energy conservation equation for gas phase is

∂ ⎯ Y )= u (εs ρs Ys, j ) + ∇·(εs ρs ⎯→ s s, j ∂t

M

∑

ṁ s, jj′− ∑ ṁ sjgi

j′ ≠ j

2.1.2.1. Zehner-Bauer-Schlünder (ZBS) model. According to ZBS model, the effective thermal conductivities of gas phase k g, e and solid phase ks, e can be formulated as

kg, e =

∂ ⎯u + ∇ ·(ε k ∇T ) ⎯u H ) = ε ∂p + τ : ∇ ⎯→ (εg ρg Hg ) + ∇ ·(εg ρg ⎯→ g g g g g g g, e g ∂t ∂t −hgs (Ts−Tg ) + Sradi, g + SH2O

ṁ s jg i

i=1

2.1.2. Constitutive relations for heat transfer According to Yaji and Kunni [25], the following mechanisms contribute the heat transfer in packed bed with the presence of flowing fluid: conduction through solid phase, radiation through solid phase, conduction through gas phase, radiation through gas phase, and interphase convection heat transfer. On the basis of our previous studies on the evaluation of various conduction and radiation heat transfer models for packed bed [23], in this work the thermal conductions through solid and gas phases were accounted for using the ZehnerBauer-Schlünder (ZBS) model [26,27], the effective radiative conductivity due to solid radiation was calculated through Breitbach and Barthels (B-B) model [28], the heat transfer due to gas radiation was described through discrete ordinates (DO) model [29], and the conventional Gunn's heat transfer model [30] was adopted to calculate the convection heat transfer coefficient. Details of these models are given below.

(2)

ṁ sg is the mass source term. It denotes the mass transfer between solid and gas phases due to chemical reaction or phase transformation. In our simulations, it was expressed as the sum of mass transfer from solid phase to species i in gas phase, ṁ sg = ∑i ṁ sg i . εg and εs are volume fractions of gas and solid phases, εg + εs = 1. The momentum conservation equation for the gas phase is ∂ ⎯u ) + ∇·(ε ρ ⎯→ ⎯ ⎯→ ⎯ → (εg ρg ⎯→ g g g ug ug ) = −εg ∇p + ∇ ·(εg τg ) + εg ρg g ∂t ⎯ − ⎯→ ⎯u ) + ṁ ( ⎯→ ⎯ ⎯→ ⎯ u + βgs ( ⎯→ s g sg us − ug )

M

∑

(7) j=1 ⎯⎯⎯→ where the diffusion flux of species i in gas phase is Jg, i = ρg Di, m ∇Yg, i , Di, m is the effective diffusion coefficient for gas i through a mixture of gases. ṁ s jg i represents mass transfer from species j of solid phase to species i of gas phase. ṁ s, jj′ is the mass flux between different species due to a homogeneous reaction in solid phase.

2.1.1. Governing equations The continuity equations for gas (g) and solid (s) phases are expressed as follows:

∂ ⎯u ) = ṁ (εg ρg ) + ∇ ·(εg ρg ⎯→ g sg ∂t

∂ ⎯u Y ) = −∇ ·(ε ⎯⎯⎯J→) + (εg ρg Yg, i ) + ∇ ·(εg ρg ⎯→ g g, i g g, i ∂t

(6)

Solid phase

Gas phase

(5)

where ks, e is the effective thermal conductivity of solid phase, Sradi, s is the heat source due to particle radiation. The heat source term Sreaction represents the heat source due to coal pyrolysis. Its value can be obtained from Sreaction = Qreaction·ṁ reaction . Parameter ṁ reaction is the mass transfer flux due to coal pyrolysis, and parameter Qreaction is the reaction heat. Note that in our simulations the reaction heat of pyrolysis was added into the solid phase. Following the strategy adopted by Shu et al. [16], in our simulations both the gas and solid phases were treated as mixture composed of certain species. The gas phase was assumed to be a mixture of volatile products and composed of H2, CO, CO2, CH4, C2H6, H2O (g) and tar. And the solid phase consisted of coal, char and H2O (l). Transport equations for species in gas phase (species i) and solid phase (species j) can be expressed as:

Γ= (4)

(1− 1−εg ) k g εg

, ks, e =

1−εg [ωκ + (1−ω)] k g 1−εg

⎤ 2 ⎡ κ−1 B κ B−1 1 ⎢ ln ⎛ ⎞− B − ( B+ 1)⎥ 2 B B κ ⎝ B ⎠ 1− 2 ⎥ 1− κ ⎢ 1− κ κ ⎦ ⎣

( )

(8)

(9)

10/9

where k g, e is the effective thermal conductivity of gas phase, hgs is the gas-solid heat transfer coefficient, Sradi, g is the heat source due to gas radiation. The heat source term SH2O represents the heat transfer due to the phase transformation of moisture, which can be obtained from SH2O = ṁ H2O·h H2O . ṁ H2O is the mass transfer of moisture between phases, and h H2O represents the latent heat of vaporization. The energy conservation equation for solid phase is

1−εg ⎞ B= 1.25 ⎜⎛ ⎟ ⎝ εg ⎠

(10)

κ = ks / k g . Parameters k g and ks are the thermal conductivities of gas mixture and solid particles, respectively. Parameter ω is a surface fraction parameter reflecting the heat transfer through particle contact and its value depends on many undetermined quantities such as surface 687

Fuel 233 (2018) 685–694

Y. Qian et al.

coal → aTar + bCH 4 + cC2 H6 + dCO + eCO2 + f H2 + g H2 O+ hChar

roughness and mechanical properties of particle[26]. Literature survey shows that the value of ω is usually set to be 7.26 × 10−3 [31,32]. This value was also used in our simulations.

(17) The pyrolysis products were divided into two groups: volatiles and solid char. Volatiles included tar and light gas, and the light gas was composed of H2, CO, CO2, CH4, C2H6 and H2O. It should be noted that both tar and light gas were in gaseous state at the operating temperature, and the coal and generated char were in the solid phase. According to Bradley et al. [42], the devolatilization rate of coal can be expressed as

2.1.2.2. Breitbach and Barthels (B-B) model. The effective radiative conductivity approach has been commonly used to describe the heat transfer due to solid radiation [33–35]. In this approach the radiative heat flux is formulated as a set of simple algebraic equations similar to the Fourier conduction law [33],

qradi, s = −ksr, e ∇Ts

(11)

ṁ coal → product =

ksr, e

is the effective radiative conductivity of solid phase. The general form of ksr, e is [33]

ksr, e = 4FE σdp Ts3

where σ is the Stefan-Boltzmann constant, dp is the particle diameter which was assumed to be constant during coal pyrolysis, FE is radiation exchange factor. The key of effective radiative conductivity approach is to get the appropriate formula of FE . Various empirical correlations of FE have been reported in literature, as have been nicely reviewed and evaluated by van Antwerpen et al. [33]. In this work, the model proposed by Breitbach and Barthels [27] was adopted, as shown below,

Af =

ks 4σdp Ts3

E kv = A 0exp ⎛− 0 ⎞ ⎝ RTs ⎠ ⎜

. Parameter εr is the surface emissivity of particle. The

empirical parameter B can be calculated through Eq. (10). Then the solid radiation source Sradi, s = −∇ ·qradi, s .

ṁ coal → species(i) = 2.1.2.3. Discrete ordinates (DO) model. The DO model available in Fluent has been widely used to describe the heat transfer due to gas radiation [36–39]. DO model considers the radiative transfer equation → at position r in the direction → s as a field equation

where I is the radiation intensity, a is the absorption coefficient of gas phase. The value of a could be determined through the gray gas model in which the absorption coefficient of the gas is assumed constant in the entire spectrum. The gray gas model has been commonly used in generic CFD programs for modeling radiative transfer [40]. In this work, following the works of Crnomarkovic et al. [38,39], the value of a was set to be 0.07 m−1. The gas radiation source Sradi, g = −∇ · ∫4π I (→ r,→ s )→ s dΩ. Ω is the solid angle.

6k g εs εg Nus

+

+ (1.33−2.4εg +

1.2εg2) Res0.7 Pr 1/3 (16)

Res is the Reynolds number and Pr is Prandtl number: Res =

i (20)

(21)

When the temperature of liquid water Tl is higher than the saturation temperature Tsat , evaporation of liquid water takes place. And condensation of water vapor occurs when its temperature Tv is lower than the saturation temperature. To model the phase transition of water, the evaporation-condensation model available in Fluent was

Nus is the Nusselt number Nus = (7−10εg +

∑

H2 O(l) ↔ H2 O(g )

(15)

0.7Res0.2 Pr 1/3)

i h

2.1.4. Moisture evaporation and condensation model The raw coal particles contain a certain amount of physical water. Ideally, during the heating process the coal particles are firstly heated to the water boiling point. The temperature will keep constant until the moisture content inside them is completely evaporated. The water vapor might condense again when it flows through the low temperature raw coal bed. Thus the moisture content of raw coal particles significantly influences temperature evolution of coal bed, and subsequently the pyrolysis behavior [1,43]. The phase transformation of physical water content in coal particles can be expressed as

2.1.2.4. Gunn's heat transfer model. According to this model, the gassolid heat transfer coefficient hgs can be formulated as [30]

5εg2)(1

dvi = ṁ coal → product dt

ṁ coal → species(i) is the generation rate of species i due to the pyrolysis reaction. The stoichiometric coefficients a ∼ h for devolatilization are given in Table 2. They are determined through the experimental data of the product distribution of Yilan coal pyrolyzing at 1173 K [7]. Following Shu et al. [16], in our simulations the pyrolysis reaction heat Qreaction was set to be −0.425 MJ/kg (coal).

(14)

dp2

(19)

i=a

4

hgs =

⎟

To get the values of the pre-exponential factor A 0 and apparent activation energy E0 , non-isothermal kinetics analyses of coal pyrolysis based on thermogravimetric analysis (TGA) experiments were conducted. The tested coal sample was Yilan coal, which was also the coal used in the experiments of Zhang et al. [7]. Table 1 lists the major characterization data for the used coal. According to the thermogravimetric (TG) analysis of the Yilan coal, the values of A 0 and E0 were 29.07 s−1 and 63067 J/mol, respectively. The entire devolatilization rate can be allocated to each component species respectively in proportion to its stoichiometric coefficient as [16]

(13)

σTg ∇ ·(I (→ r,→ s )→ s ) + aI (→ r,→ s)=a π

(18)

v∞ is the potential ultimate amount of the volatile products, v is the amount of volatile products at specific time t. According to the Arrhenius law the rate constant Kv can be expressed as

(12)

⎧ ⎫ 1−εg B+1 1 FE = (1− 1−εg ) εg + × · 1 ⎨ 2/ εr −1 B 1 + (2 / ε − 1) ⎬ r f ⎭ ⎩

dv = k v (v∞−v ) dt

dp | ug − us | ρg , μg

Cpg μg

Table 1 Properties of Yilan coal sample.

Pr = k . Parameter μg is the viscosity of gas phase and Cpg is the g specific heat of gas phase.

proximate analysis (wt%)

2.1.3. Pyrolysis kinetic models In this work, the coal pyrolysis mechanism was described using the following global devolatilization model [16,41]

ultimate analysis (wt%)

688

M 4.61 Cd 70.96

A 33.36 Hd 6.23

V 32.15 Nd 1.57

FC 29.88 Od 20.61

Sd 0.63

Fuel 233 (2018) 685–694

Y. Qian et al.

Table 2 Stoichiometric coefficients for devolatilization.

Table 3 Properties of solid and gas phases.

Species

tar

CH4

C2H6

CO

CO2

H2

H2 O

char

Properties of solid phase

value

label value

a 0.0564

b 0.0432

c 0.0195

d 0.0295

e 0.0457

f 0.0060

g 0.0859

h 0.7138

Coal mixture Average diametera, mm Density, kg/m3 Specific heat, J/(kg·K) Thermal conductivity, W/ (m·K) Voidagea Properties of coal Densitya, kg/m3 Specific heat, J/(kg·K)

coal char H2O (l) 2 volume-weightedb Ts mass-weightedc mass-weightedd

Thermal conductivity, W/ (m·K) Emissivity Properties of char Densitya, kg/m3 Specific heat, J/(kg·K) Thermal conductivity, W/ (m·K) Emissivity Properties of physical liquid water Density, kg/m3 Specific heat, J/(kg·K) Thermal conductivity, W/ (m·K) Properties of gas phasee Mixture of pyrolysis volatiles Density, kg/m3 Specific heat, J/(kg·K) Thermal conductivity, W/ (m·K)

Fig. 1. Schematic plot for experimental reactor of coal pyrolysis.

a b

adopted to calculate the evaporation and condensation rates [44].

For evaporation,

For condensation,

T −T Tl > Tsat, ṁ H2O = Cev ·⎡ρl εl l sat ⎤ ⎥ ⎢ Tsat ⎦ ⎣ T −T Tv < Tsat ṁ H2O = Ccon ·⎡ρv ε v v sat ⎤ ⎢ Tsat ⎥ ⎦ ⎣

(22)

1400 1150 Ts 573K ; 1150 + 2.03×(Ts -573)0.00155×(Ts −573)2 Ts 573K 0.19 Ts ⩽ 573K ; 0.19 + 0.00025×(Ts -573) Ts > 573K 0.9

[48]

[49] [50]

999 1381.6 0.37

[51] [52]

0.7

[50]

998.2 4182 0.6

tar CH4 C2H6 CO CO2 H2 H2O volume-weightedb mass-weightedc mass-weightedd

Experimental measurement. The volume-weighted density of the mixture is calculated by ρ =

1 Y ∑i i ρi

.

d The mass-weighted thermal conductivity of the mixture is calculated by k = ∑i Yi ki . e

(23)

The properties of each species of gas phase are in supplementary material.

inner diameter of 100 mm and effective volume of 1500 mL for coal samples. In the experiments the reactor loaded with coal was quickly placed into an electrically heating furnace with preset temperature being 1173 K. The temperature evolution of the coal bed was measured through a thermocouple with its free end placed at the center of the bed. And the exit of reactor was connected with a set of product collection devices to collect pyrolysis gas and tar. The cylindrical reactor was numerically discretized using hexahedral structured grids. The maximal grid size was 2 mm. And the grids were found to be fine enough to achieve grid independent simulation results. The initial temperatures of gas and solid phases were set to be 300 K. Constant temperature boundary (1173 K) was set for the heating wall of the reactor, whose temperature was assumed to be the same with that of furnace. And all other walls were set to be adiabatic. The exit at the top of reactor was set as atmosphere pressure outlet. To reduce numerical diffusion, the second-order upwind scheme was used to discretize the equations of momentum, energy and discrete ordinates. And mass conservation equations were solved by the QUICK scheme. Table 3 summarizes the main properties of gas and solid phases and also the parameter settings used in simulations. The initial species mass fraction of solid phase was determined by the proximate analysis of Yilan coal as shown in Table 3, where the coal mixture consisted of 4.61 wt% H2O(l), 95.39 wt% coal content (including ash), and zero percent of char content. All the required functions describing the temperature dependent material properties and the source terms

Tsat

∫ Cp, l dT + hH2O, ref /MH2O + ∫ Cp, v dT , where Cp, l and Cp, v are specific Tref

0.4

c The mass-weighted specific heat of the mixture is calculated by Cp = ∑i Yi Cp, i .

where ρl and ρv are densities of liquid water and water vapor, εl and ε v are volume fractions of liquid water and water vapor, and Cev and Ccon are evaporation frequency and condensation frequency, respectively. Cev (C(Tcon ) is a function of evaporation (condensation) surface, evaporation (condensation) coefficient, latent heat of water, etc (see supplementary material). As stated in Fluent theory guide, the evaporation (condensation) surface and evaporation (condensation) coefficient are usually not very well known. Thus Cev and Ccon generally serves as tuning parameters. According to the theory guide of Fluent, the latent heat h H2O at the saturation temperature (373 K at the atmospheric oph H2 O = − eration pressure) is calculated from Tsat

Refs.

Tref

heat of liquid water and vapor. is the reference temperature, 298 K. MH2O is the molecular weight of water, 18 kg/kmol. And hH2O, ref is the latent heat at the reference temperature. Following Shu et al. [16], its value was set as the difference between standard state enthalpies of vapor and liquid phase for 4.4 × 107 J/kmol under 0.1 MPa and 298 K conditions. 2.2. Simulation layout The simulated system has been experimentally investigated by Zhang et al. [7]. Fig. 1 gives a schematic diagram of the experimental apparatus. The cylindrical reactor was made of stainless steel with the 689

Fuel 233 (2018) 685–694

Y. Qian et al.

Fig. 2. (a) Comparison of the predicted and experimentally measured temperature evolutions at the center of coal charge, (b) evolutions of solids temperature at different radial positions and (c) contours of solids temperature at different moments.

radial positions. All the monitor points located at the axial position of 75 mm above the base wall. It can be seen that the heating rates of the coal particles varied greatly along the radial direction. The closer is the radial region to the heating wall, the higher is the heating rate. In the temperature range of 373 K and 773 K, the average heating rate was around 130 K/min for the coal zone 5 mm away from the heating wall, whereas it was only about 10 K/min for coal zone in the central region. In addition, in the region close to the heating wall, the temperature plateau of 373 K did not appear. But in the central region, there was a long temperature plateau which nearly took half of the whole heating time. The reason is that in the near wall region, only the moisture evaporation takes place, and the moisture content decreases rapidly, whereas in the central region re-condensation of water vapor dominated at the early stage of heating, as can be seen from Fig. 3. Fig. 2(c) displays contours of temperature at different moments. It can be seen that the temperature gradient along radial direction was very large at the early stage. And the temperature gradient decreased gradually with the ongoing of heating. The evolutions of the moisture content of the coal charge are shown in Fig. 3. The initial moisture content was 4.61 wt% (air dried). Fig. 3(a) displays the moisture content fields along a cross section of the investigated reactor at different moments. Due to the high temperature of reactor wall, a large amount of heat was quickly transferred into the coal charge adjacent to the heating wall. This led to the rapid moisture evaporation in that zone. The dry coal zone slowly broadened with the ongoing of heating. Meanwhile, the moisture content in the central zone firstly increased (t ≤ 20 min) and then decreased (t ≥ 30 min). The moisture evaporation and condensation can be more clearly seen from Fig. 3(b), which shows the evolutions of moisture content at different radial positions. It can be found that the water vapor released from outer zones of coal charge partially accumulated in the inner zones. This is a result of the transport of water vapor. The hot water

appeared in the transport equations that are not available in Fluent 6.3.26 were implemented into Fluent with user-define function (UDF).

3. Results and discussion As mentioned in Section 2.1.4, the temperature evolution of coal bed is closely related to the moisture evaporation and condensation. Unfortunately, precise determination of the values of the evaporization (Cev ) and condensation (Ccon ) frequencies of moisture is generally very difficult. The measurements conducted by Rubel and Gentry [45] suggested that accommodation coefficient of condensation was approximately 1.2 times evaporization. Thus, in this work we also kept Ccon = 1.2Cev as suggested by Zhang [46] and fitted the values of Cev and Ccon by comparing the simulation results with experimentally measured data, as shown in the figure embedded in Fig. 2(a). It can be seen that the case with Cev = 5 and Cev = 6 gives the best agreement with the experimental data. Thus, in all the following discussed simulation cases these two values were used. Fig. 2(a) gives the comparison of simulated and experimental measured temperature evolutions at the center of the investigated reactor. The experimentally observed constant temperature stage in the initial stage (t < 10 min) has not been correctly represented by our simulation. Nevertheless, the intermediate constant temperature stage and the later temperature increasing stage were successfully captured by the simulation and the predicted and measured temperature changes were in quantitatively agreement. The temperature plateau in between 20 min and 50 min represents the evaporation stage. Once the moisture had all evaporated (t > 50 min) the temperature of coal started to increase. The predicted time for central coal reaching 773 K was around 90 min, which is consistent with the actual pyrolysis time measured by Zhang et al. [7]. Fig. 2(b) gives the temperature evolutions of coal charge at different 690

Fuel 233 (2018) 685–694

Y. Qian et al.

Fig. 3. (a) Moisture contents in the axial-section of the pyrolyzer at different moments, (b) evolutions of moisture content at different radial positions (solid lines: temperature of coal charge, dash lines: moisture content).

vapor re-condensed when flowing through the coal charge zone which still remained at the temperature below 373 K, leading to the increase of moisture content. Similar phenomena have also been observed in many studies, e.g. [20,22]. Fig. 4(a) gives the distributions of volatiles generation rate at different moments. It can be seen that the radial position of most intense volatiles release slowly shifted to reactor center with the progressing pyrolysis. At the same time, the zone of volatiles release broadened. This can be attributed to the fact that the heating rate in the inner coal zone was smaller than that in the region close to heating wall, as manifested in Fig. 2(b). Fig. 4(b) presents the time evolution of volatiles generation rate at different radial positions, together with the temperature evolutions. The simulation results show that the zones of most intense volatiles release were in the temperature range between approximately 600 K and 800 K. It can also be found that, from heating wall to reactor center, the peak of volatiles generation rate first decreased and then increased. Such results have also been reported by Polesek-Karczewska et al. [20] who modeled the carbonization of coal in a coke oven using three phase transient model. The intensity of volatiles release is directly related to the heating rate of the coal bed. Based on the predicted temperature evolutions at different radial positions, it is found that from the heating wall to reactor center the heating rate decreased from 130 K/min (5 mm) to 8 K/min (30 mm) and then increased to 10 K/min (50 mm). The re-increase of heating rate may be attributed to that the moisture evaporation was complete all along the bed, leaving behind coal/char partly heated up. Such phenomenon has also been reported by [20].

Due to the evaporation of liquid water and the release of volatiles, the voidage of coal bed continually changed during the heating process. Fig. 5(a) presents the voidage evolution. For the coal zone adjacent to the heating wall, its voidage sharply increased at the early heating stage due to the high heating rate and thus the intense release of gaseous products (vapor and volatiles). Its voidage continually increased at the later heating stage but the increasing rate became notably smaller. The slow increase of voidage is due to the slow release of volatiles. Note that according to Eq. (18) the volatiles generation rate can never achieve zero. And the slow release of volatiles at the later heating stage can also be seen from Fig. 4(a). The variation of voidage in the inner coal zone shows several stages. The voidage first slightly decreased due to the condensation of water vapor, as illustrated in Fig. 3. As the moisture evaporation front moved forwards from the heating wall to reactor center, the voidage in the inner coal zone slowly increased. Once the coal charge drying process finished, the temperature of the coal charge quickly increased, leading to the intense release of volatiles, which resulted into the rapid increase of voidage. It should be addressed again that the Yilan coal considered in this work is a type of sub-bituminous coal with high volatile content. The softening and thermoplastic deformation which commonly appear during the heating of coking coal in coke oven [20,21] do not take place in our modeled system. The non-uniform distribution of voidage along the radial direction of the considered cylindrical reactor directly influenced the flow behavior of gaseous products. Fig. 5(b) shows the altitudinal direction (zdirection, see Fig. 1) velocity distributions and streamline patterns of gaseous products at different heating time. Due to the higher voidage of 691

Fuel 233 (2018) 685–694

Y. Qian et al.

Fig. 4. (a) Contours of volatiles generation rate at different moments, (b) evolutions of temperature (solid lines) and volatiles generation rates (dash lines) at different radial positions.

phenomenon happens, its influence on heat transfer is minor since in this region the dominant heat transfer mechanism is radiation. And most of the effective radiative conductivity models show that the effective radiative conductivity increases with the increase of voidage [33,34]. Qualitatively, the simulation results are consistent with the findings of Slupik et al. [22]. They found that the shrinkage of coke and thus the formation of gap between heating wall and coke charge had ignorable influence on the heat transfer. In the central region, the dominant heat transfer mechanism at the early heating stage (t ≤ 50 min) was conduction. As heating progressed, the temperature of coal charge increased, leading to the increase of effective radiative conductivity (see Fig. 6(b)). The effective thermal conductivity due to conduction decreased slightly due to the increase of voidage. Thus, at the later heating stage (t ≥ 70 min), the dominant heating mechanism changed to be radiation.

the zone near heating wall, the gaseous products generated in the relatively lower temperature inner coal zone tended to flow towards the higher temperature heating wall before flowing out of the reactor. In real coal pyrolysis system, due to the secondary reactions of primary tar, primary pyrolysis gas product passing through high temperature zone will lead to a significant decrease of tar yield [47]. Thus, for the purpose of increasing the tar yield, it is better to regulate the flow of primary pyrolysis gas from high temperature zone to low temperature zone. This is exactly the reason that Zhang et al. [7] set a central gas collection pipe connected to the reactor exit. In our next work, we’ll investigate the influence of different internals on the heating rate of coal charge, the flow behavior of gaseous products, and the yield of tar by including a secondary reaction model of tar in simulations. For the investigated fixed bed reactor, the coal charge gained heat mainly via two modes, i.e., conduction and radiation. It should be noted that for the case considered here, the temperature difference between gas and solid phases is negligible (see Fig. A1 in Supplementary material). Thus the interphase convective heat transfer is ignorable. Fig. 6(a) compares their relative contributions. It can be observed that in the region close to heating wall the main heat transfer mechanism was radiation. Around 80% heat flux was contributed by radiation. This is due to the high temperature of heating wall and also the relatively larger voidage in that region, leading to far larger effective radiative conductivity, as shown in Fig. 6(b). In real fixed bed reactor, the voidage of coal zone adjacent to reactor wall is generally larger than that of coal zone in inner zone. And the release of gaseous products and the possible shrinkage of coal particle will further increase the voidage. In some cases, for example, the coke oven system, this might lead to the formation of gap between reactor wall and char charge [20,21]. The simulation results shown in Fig. 6 suggests that, even if this

4. Conclusions In this work, a two-fluid CFD numerical framework was developed to simulate the coal pyrolysis process in a cylindrical fixed-bed reactor. The complex physical and chemical conversions during the pyrolysis of Yilan coal (a kind of sub-bituminous coal) was characterized by coal heating, moisture evaporation, and coal devolatilization. The CFD framework was successfully validated by comparing the predicted temperature evolution in the reactor center with experimentally measured data. The main findings of this study are summarized as below: 1. The simulation results show that the moisture content determines the coal pyrolysis process. The water vapor generated in the charge zone close to heating wall recondensed in the cooler coal zone when 692

Fuel 233 (2018) 685–694

Y. Qian et al.

Fig. 5. (a) Voidage contours at different moment, (b) altitudinal direction gas velocity distributions and streamlines at different moments.

Fig. 6. (a) The contributions of conduction and radiation; (b) evolutions of the effective conductivities.

flowing towards to the reactor center, leading to the delay of the temperature increase above the boiling point and consequently the start of pyrolysis. 2. From the heating wall to reactor center, the peak of volatiles generation rate shows a first decrease and later increase trend. Analyses indicate that the volatiles generation rate is directly related to the heating rate of coal bed. The re-increase of heating rate leads to the

later increase of volatiles generation, which might be attributed to the moisture have all evaporated. 3. Due to the evaporation of moisture and the release of volatiles, the voidage of coal zone adjacent to heating wall is larger than that of inner zone. Such non-uniform radial distribution of voidage leads to the phenomenon that gaseous products generated in the low temperature inner zone tend to flow towards heating wall before 693

Fuel 233 (2018) 685–694

Y. Qian et al.

flowing out of the reactor. 4. The simulation results indicate that in the temperature range of intense volatiles generation radiation is the dominant heat transfer mechanism.

heat transfer in coke ovens. Appl Therm Eng 2015;81:353–8. [22] Slupik L, Fic A, Bulinski Z, Nowak AJ, Kosyrczyk L, Labojko G. CFD model of the coal carbonization process. Fuel 2015;150:415–24. [23] Qian YN, Geng SJ, Zhan J-H, Xu GW. Liu XX. CFD modeling of heat transfer in packed beds: evaluation of models. Submitted to Int. J Heat Mass Transf 2017. [24] Gidaspow D. Multiphase flow and fluidization: continuum and kinetic theory descriptions. Academic Press; 1994. [25] Yagi S, Kunii D. Studies on effective thermal conductivities in packed beds. AIChE J 1957;3(3):373–81. [26] Zehner P, Schlünder EU. Wärmeleitfähigkeit von Schüttungen bei mäβigen Temperature. Chem Ing Tech 1970;42:933–41. [27] Bauer R, Schlünder EU. Part I: Effective radial thermal conductivity of packings in gas flow. Part II: Thermal conductivity of the packing fraction without gas flow. Int Chem Eng 1978;18:189–204. [28] Breitbach G, Barthels H. The radiant heat transfer in the high temperature reactor core after failure of the afterheat removal systems. Nucl Technol 1980;49:392–9. [29] Modest MF. Radiative heat transfer. Academic Press; 2003. [30] Gunn D. Transfer of heat or mass to particles in fixed and fluidised beds. Int J Heat Mass Transf 1978;21:467–76. [31] Yusuf R, Halvorsen B, Melaaen MC. An experimental and computational study of wall to bed heat transfer in a bubbling gas–solid fluidized bed. Int J Multiphas Flow 2012;42:9–23. [32] Patil DJ, Smit J, van Sint Annaland M, Kuipers JAM. Wall-to-bed heat transfer in gas–solid bubbling fluidized beds. AIChE J 2006;52:58–74. [33] Van Antwerpen W, du Toit CG, Rousseau PG. A review of correlations to model the packing structure and effective thermal conductivity in packed beds of mono-sized spherical particles. Nucl Eng Des 2010;240:1803–18. [34] Zhou CG, Yang WH. Effect of heat transfer model on the prediction of refuse-derived fuel pyrolysis process. Fuel 2015;142:46–57. [35] Yang YB, Ryu C, Khor A, Sharifi VN, Swithenbank J. Fuel size effect on pinewood combustion in a packed bed. Fuel 2005;84(16):2026–38. [36] Zavala-Guillen I, Xaman J, Salinas C, Ismail KAR, Hernandez-Perez I, HernandezLopez I. Optical thickness effect on natural convection in a vertical channel containing a gray gas. Int J Heat Mass Transf 2017;107:510–9. [37] Zhang J, Prationo W, Zhang L, Zhang ZX. Computational fluid dynamics modeling on the air-firing and oxy-fuel combustion of dried Victorian brown coal. Energy Fuels 2013;27(8):4258–69. [38] Crnomarkovic N, Sijercic M, Belosevic S, Stankovic B, Tucakovic D, Zivanovic T. Influence of forward scattering on prediction of temperature and radiation fields inside the pulverized coal furnace. Energy 2012;45:160–8. [39] Crnomarkovic N, Sijercic M, Belosevic S, Tucakovic D, Zivanovic T. Numerical investigation of processes in the lignite-fired furnace when simple gray gas and weighted sum of gray gases models are used. Int J Heat Mass Transf 2013;56:197–205. [40] Liu F, Becker HA, Bindar Y. A comparative study of radiative heat transfer modelling in gas-fired furnaces using the simple grey gas and the weighted-sum-of-greygases models. Int J Heat Mass Transf 1998;41:3357–71. [41] Mahalatkar K, Kuhlman J, Huckaby ED, O’Brien T. CFD simulation of a chemicallooping fuel reactor utilizing solid fuel. Chem Eng Sci 2011;66:3617–27. [42] Bradley D, Lawes M, Park H-Y, Usta N. Modeling of laminar pulverized coal flames with speciated devolatilization and comparisons with experiments. Combust Flame 2006;144:190–204. [43] Hu EF, Zeng X, Ma DC, Wang F, Yi XJ, Li Y, et al. Effect of the moisture content in coal on the pyrolysis behavior in an indirectly heated fixed-bed reactor with internals. Energy Fuels 2017;31(2):1347–54. [44] Lee WH. A pressure iteration scheme for two-phase flow modeling. Veziroglu TN, editor. Multiphase Transport Fundamentals, Reactor Safety, Applications, vol. 1. Washington, DC: Hemisphere Publishing; 1980. [45] Rubel GO, Gentry JW. Measurement of the kinetics of solution droplets in the presence of adsorbed mono-layers: determination of water accommodation coefficients. J Phys Chem 1984;88(14):3142–8. [46] Zhang XB, Li JF, Zhu JK, Qiu LM. Computational fluid dynamics study on liquefied natural gas dispersion with phase change of water. Int J Heat Mass Transf 2015;91:347–54. [47] Franklin HD, Peters WA, Howard JB. Mineral matter effects on the rapid pyrolysis and hydropyrolysis of a bituminous coal. 1. Effects on yields of char, tar and light gaseous volatiles. Fuel 1982;61(2):155–60. [48] Agroskin AA, Gonczarow EI, Makeev LA, Jakunin WP. Thermal capacity and heat of pyrolysis of Donbass coal. Koki Chimija 1970;5:8–13. [49] Agroskin AA. The change of heat and temperature transfer coefficient of coal during heating. Bergakad Freiberg 1957;9:177–86. [50] Bhattacharya SP, Wall TF. Development of emittance of coal particles during devolatilisation and burnoff. Fuel 1999;78:511–9. [51] Kobayashi H. Devolatilization of pulverized coal at high temperatures. Massachusetts Institute of Technology; 1976. [52] Gupta M, Yang J, Roy C. Specific heat and thermal conductivity of softwood bark and softwood char particles. Fuel 2003;82:919–27.

Acknowledgements The work presented in this paper is financially supported by National Natural Science Foundation of China, Grant No. 21576265, Nippon Steel Arts Foundation, Grant No. U1630105, and National Basic Research Program of China, Grant No. 2014CB744303. Xiaoxing Liu acknowledges the financial support from the “Hundred Talents Program” of Chinese Academy of Sciences. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.fuel.2018.06.100. References [1] Zhang C, Wu RC, Hu EF, Liu SY, Xu GW. Coal pyrolysis for high-quality tar and gas in 100 kg fixed bed enhanced with internals. Energy Fuels 2014;28:7294–302. [2] Zhang Y. Analysis on metamorphism of coal pyrolysis technology. Chem Ind 2012;30(3):23–5. [3] Guo S, Luo C, Zhang D, Han Z, Liu H, Kang X, et al. Experiment in pilot plant of new technology for lignite retorting using solid heat carrier. J Dalian Univ Technol 1995;1:46–50. [4] Zhang JW, Wu RC, Zhang GY, Yu J, Yao CB, Wang Y, et al. Technical review on thermochemical conversion based on decoupling for solid carbonaceous fuels. Energy Fuels 2013;27(4):1951–66. [5] Khan MR. Production of high quality liquid fuels from coal by mild pyrolysis of coallime mixtures. Fuel Sci Technol Int 1987;5(2):185–231. [6] Hu EF, Zeng X, Ma DC, Wang F, Li Y, Guo EW, et al. Characterization of coal pyrolysis in indirectly heated fixed bed based on field effects. Fuel 2017;200:186–92. [7] Zhang C, Wu RC, Xu GW. Coal pyrolysis for high-quality tar in a fixed-bed pyrolyzer enhanced with internals. Energy Fuels 2014;28:236–44. [8] Hu EF, Zhu CQ, Rogers K, Han X, Wang J, Zhao J, et al. Coal pyrolysis and its mechanism in indirectly heated fixed-bed with metallic heating plate enhancement. Fuel 2016;185:656–62. [9] Xiong Q, Aramideh S, Kong SC. Modeling effects of operating conditions on biomass fast pyrolysis in bubbling fluidized bed reactors. Energy Fuels 2013;27:5948–56. [10] Xiong Q, Aramideh S, Kong SC. Assessment of devolatilization schemes in predicting product yields of biomass fast pyrolysis. Environ Prog Sustain Energy 2014;33:756–61. [11] Hashimoto N, Kurose R, Shirai H. Numerical simulation of pulverized coal jet flame employing the TDP model. Fuel 2012;97:277–87. [12] Jeong HJ, Hwang IS, Park SS, Hwang J. Investigation on co-gasification of coal and biomass in Shell gasifier by using a validated gasification model. Fuel 2017;196:371–7. [13] Nikolopoulos A, Malgarinos I, Nikolopoulos N, Grammelis P, Karrelas S, Kakaras E. A decoupled approach for NOx−N2O 3-D CFD modeling in CFB plants. Fuel 2014;115:401–15. [14] Hart J, Al-Abbas AH, Naser J. Numerical investigation of pyrolysis of a Loy Yang coal in a lab-scale furnace at elevated pressures. Heat Mass Transfer 2013;49(12):1725–32. [15] Iavarone S, Galletti C, Contino F, Tognotti L, Smith PJ, Parente A. CFD-aided benchmark assessment of coal devolatilization one-step models in oxy-coal combustion conditions. Fuel Process Technol 2016;154:27–36. [16] Shu Z, Fan CG, Li SG, Wang JW. Multifluid modeling of coal pyrolysis in a downer reactor. Ind Eng Chem Res 2016;55:2634–45. [17] Yan BH, Cheng Y, Jin Y, Guo CY. Analysis of particle heating and devolatilization during rapid coal pyrolysis in a thermal plasma reactor. Fuel Process Technol 2012;100:1–10. [18] Ma J, Zhang M, Su BG, Wen GD, Yang YW, Yang QW, et al. Numerical simulation of the entrained flow hydropyrolysis of coal in magnetically rotating plasma reactor. Energy Covers Manage 2017;148:431–9. [19] Guo ZC, Tang HQ, Liu JL. Desulfurization of coke by recycling COG in coking process. Fuel 2005;84:893–901. [20] Polesek-Karczewska S, Kardas D, Cizminski P, Mertas B. Three phase transient model of wet coal pyrolysis. J Anal Appl Pyrol 2015;113:259–65. [21] Lin W, Feng YH, Zhang XX. Numerical study of volatiles production, fluid flow and

694