- Email: [email protected]
Direct numerical simulation (DNS) of alkali metals released during char combustion
Fuel 255 (2019) 115763
Contents lists available at ScienceDirect
Fuel journal homepage: www.elsevier.com/locate/fuel
Full Length Article
Direct numerical simulation (DNS) of alkali metals released during char combustion Sibo Qu, Changfu You
T
⁎
Key Laboratory for Thermal Science and Power Engineering of the Ministry of Education, Tsinghua University – University of Waterloo Joint Research Centre for Micro/ Nano Energy & Environment Technologies, Department of Energy and Power Engineering, Tsinghua University, Beijing 100084, China
A R T I C LE I N FO
A B S T R A C T
Keywords: Alkali metals release Char combustion Direct numerical simulation (DNS) Sodium chloride (NaCl)
Alkali metals in pulverized coal easily released in gas phase, cause significant fouling, atmospheric pollution and etc. during combustion and gasification. A TURE DNS method, based on fictitious domain strategy with finite volume (FV) scheme, was developed to describe the behavior of alkali metals release and transport in reacting fluid–solid flow at sub-particle level. Taking water-soluble sodium chloride (NaCl) of high harmfulness as an example, alkali metals release and transport in physical field of char combustion, as the main stage of coal combustion, was simulated. A 50 × 10−6 m single particle, and 150 1 × 10−3 m closely spaced particles combustion at 1273 K air environment were respectively calculated. Computational results revealed that continuous release of alkali metals depends heavily on the local high temperature from the stable combustion stage (approximately 85% of the whole combustion stage) which produces stable heat. And complicated movement of fluid medium has a great effect on the transport behavior of alkali metals in gas phase, especially the entrainment of vortexes, which results in the local enrichment of alkali metals in gas phase.
1. Introduction Alkali-rich coals, e.g., Zhundong Coal from Xinjiang, China, because of its huge deposit, high quality characterized by low ash, aluminum, iron, sulfur, and trace element contents [1], has wide prospects in industrial application and scientific research. However, the high content of alkali metals is the main limitation. Existed researches showed that Na content in ash is generally more than 5%, which is significantly more than the level in typical Chinese coal, ranging from 1 to 2 % [2]. Alkali metals have the strongest sublimation tendency, for the main mineral content [3] during combustion and gasification. Especially Na [4], is generally accepted as the primary cause of serious ash-related problems, such as fouling, slagging and corrosion, for many industrial systems [5]. Further, alkali metals contribute greatly to the formation of atmospheric pollutants, such as fine particulates (size less than 1 µm) due to homogeneous nucleation or heterogenous condensation [5–7]. A chemical fractionation method of successive leaching with H2O, CH3COONH4 and HCl solutions has been widely used to quantify alkali metals using terms such as, “H2O soluble”, “CH3COONH4 soluble”, “HCl soluble”, and “insoluble form” [8]. Among them, H2O solubles such as alkali metal chlorides [33], should be given a priority consideration, due to their proven harmfulness, high content and easy conversion into gas phase [9]. Taking Na as an example, experimental investigations on ⁎
its release and reacting dynamics evolved from offline to online measuring techniques using advanced laser diagnostics [10–12]. On the numerical side, a one-step Arrhenius model firstly proposed by van Eyk et al. [13] was used to study the Na release during char combustion. A two-step kinetics model for all stages of coal combustion was developed by Liu et al. [12] then, based on the Na release, surface temperature and diameter of particles. These models can quantitatively describe the evolution of Na in solid fuels, but can not fully deal with the multiphysical fields in real fluid-solid combustion. So far, few researches have been done on kinetic behavior of alkali metals in multi-physical fields in detail after they release from coal in gas phase. A method using one-dimensional (1D) premixed flamelet generated manifolds (FGM), was used to study Na emissions in an early-stage pulverized coal flame [14]. But, the solid particle was treated as Lagrangian point source, internal structure and boundary layers of which were not resolved. So that it can only be called “quasi- DNS” or “point-particle DNS”, since it lost physical integrity. Therefore, based on existed results, a TURE DNS method, was developed in this work, to describe the behavior of alkali metals release and transport in reacting fluid–solid flow at sub-particle level in detail, further, to study the characteristics of alkali metals during combustion of solid fuels, and to provide a guide for industrial design of high-efficiency clean combustion technology. Based on fluid dynamics [15], mass and heat transfer theory [29],
Corresponding author. Tel.:/fax: +86 10 62785669. E-mail address: [email protected] (C. You).
https://doi.org/10.1016/j.fuel.2019.115763 Received 28 March 2019; Received in revised form 23 June 2019; Accepted 3 July 2019 Available online 31 July 2019 0016-2361/ © 2019 Elsevier Ltd. All rights reserved.
Fuel 255 (2019) 115763
S. Qu and C. You
2.1. Reacting fluid-solid flow model
numerical simulation is superior to the experiment in obtaining the details of multi-physical fields. There are two main types of the fully resolved simulation (FRS) method for fluid-solid flow system at subparticle level: body conformal mesh method [16] represented by Arbitrary Lagrange-Euler (ALE) [17] method, and fixed mesh method. Body conformal mesh method which constructs body conformal meshes outside solid particle, can describe the boundary condition of solid particle well, so that the computational result is more accurate. However, a moving particle leads to meshes need to be reconstructed constantly. Especially in the case of multiple particles, meshes reconstruction is frequent, which makes computational calculation tedious. And in particles collision, too small spacing between two particles can easily cause meshes distortion. Thus, it is not suitable to simulate the movement of a large number of particles with body conformal mesh method. Fixed mesh method performs well, because it does not require frequent meshes reconstruction, which mainly includes level set method [18], physalis method [19], immersed boundary method (IBM) [20], and fictitious domain method created by Glowinski et al. [21]. Level set method is mainly used to calculate the movement and deformation of the interface, thus, it is more suitable for the calculation of gas-liquid flow. Physalis method requires complicated local analytic solutions as a medium to connect the surface of solid particle to background fluid meshes. Thus, it lacks applicability for some cases without local analytical solutions. IBM is a commonly used method for FRS with fixed meshes. This method was proposed by Peskin [22] to study the flow around flexible membranes (mainly for heart). Subsequently, the method was applied to calculate the flow around suspended particles and fixed cylinders [23,24]. IBM uses a series of Lagrange points to describe the solid particle boundary, and a specific strategy to connect them with the Euler points on the background fluid meshes to realize the coupling of fluid-solid phases. Because of its better completeness in describing physical process, higher computational efficiency and stability, IBM has been widely used. The solution strategy of fictitious domain method developed by Apte [25], Ghasemi [26], Haeri [27] and Zhang [28] et al. is similar to that of IBM, which has also been widely used and developed [30]. The main difference between these two is that fictitious domain method not only deals with the boundary, but also needs to adopt a certain numerical method in interior of particle to ensure the rigid body motion law. Therefore, with the basic advantages of IBM, fictitious domain method is more dominant in describing reacting fluid–solid flow with a high density ratio. A TURE DNS method, based on fictitious domain strategy with finite volume scheme was developed. Unified structured meshes were used to represent both fluid and solid particle regions in computational domain. Here, the solid particle regions adapted to the unified equation forms by adding different source terms, named fictitious source terms. Taking water-soluble sodium chloride (NaCl) of high harmfulness as an example, alkali metals release and transport in physical field of char combustion, as the main stage of coal combustion, was simulated. Section 2 describes the algorithm development process, including governing physical equations and mathematical methods for coupled equation solutions. Section 3 validates the present model with existed experimental results at first, and the calculation results of temperature and component distribution obtained under the present method follow then, which prove feasibility, clarify the effects of local parameters, and preliminarily study the characteristics of alkali metals release and transfer in physical field of char combustion.
Physical descriptions for different phases can be divided into two types: fluid solver, and solid particle solver. Fluid solver mainly comprises applicable physical equations, dynamic equations, presented in Section 2.1.1, species conservation equation for mass transfer, presented in Section 2.1.4, and energy equation for heat transfer, presented in Section 2.1.5. Particle solver mainly comprises dynamics equations of Newton's laws and the rigid constraint, presented in Section 2.1.2, energy equation, presented in Section 2.1.5 either. Descriptions of homogeneous and heterogeneous reactions occurring in physical field of reacting fluid-solid flow are presented in Section 2.1.3, which are used as source terms. 2.1.1. Fluid dynamics equations The Navier-Stokes equations are used to describe fluid movement. Here, according to the characteristics of the fluid medium and closure of flow field, whether to adopt Navier-Stokes equations in compressible form or not can be considered, further. According to the ideal gas hypothesis, the state equation is:
ρf =
P RT
(1)
where ρf is fluid density, the unit is: kg/m ; P is pressure, the unit is: Pa; T is temperature, the unit is: K; and R is the specific gas constant of local fluid, usually the value is: 8.314 J/(mol·K). Eq. (1) is used to describe the basic relation among fluid density, ρf , pressure, P and temperature, T. According to the basic idea of taking fluid element as control body and the law of mass conservation, the continuity (mass conservation) equation of fluid is: 3
∂ρf ∂t
+ ∇∙ (ρf Uf ) = 0
(2)
where Uf is the local fluid velocity, the unit is: m/s; t is physical time, the unit is: s. According to the law of conservation of momentum, the momentum equation of fluid is:
∂ (ρf Uf ) ∂t
+ ∇∙ (ρf Uf ∙Uf )
= −∇p + ∇∙μ [∇Uf + (∇Uf )T − 2/3∇∙Uf ] + ρf g + S
(3)
where μ is fluid viscosity, the unit is: Pa·s; g is gravity acceleration, usually the value is: 9.8 m/s2; and S , the fictitious domain source term, adopting different forms for different conditions, is discussed below. 2.1.2. Solid particle dynamics equations Particle movement can be described by Newton's motion laws with specific consideration given to buoyancy effects, due to heat and frequent particle collision [28]. The correlative equations are:
∂Ut 1 = [Ff + ∂t ρs Vs
∑ Fcoll + (ρs − ρf ) Vsg]
∂ (ωs × r ) = ρs Tf / Is − |ωs |2 × r ∂t
(4)
(5)
where Ut is translational particle velocity, the unit is: m/s; ρs is particle density, the unit is: kg/m3; Vs is particle volume, the unit is: m3; Ff is force from fluid, the unit is: N; Fcoll is force from particle and wall collision, the unit is: N; ωs is particle angular velocity, the unit is: rad/s; r is vector position from particle center to focal point; Is is the moment of particle inertia, the unit is: kg·m2; and Tf is torque from fluid, the unit is: N·s. The specific counting equations of fluid force, Ff and fluid torque Tf are:
2. Numerical method To deal with complicated multi-physical field problem of reacting fluid-solid flow, a set of governing equations is established, based on the basic conservation relation and transfer theory of relevant physical quantities, in Section 2.1. And fictitious domain strategy to link particle variables with background fluid variables for realizing interaction is discussed in detail in Section 2.2. 2
Fuel 255 (2019) 115763
S. Qu and C. You
Ff =
∫ n·τds = ∫ n·{−pI + μ [∇Uf + (∇Uf )T ]} ds
(6)
Tf =
∫ r × τds = ∫ r × {−pI + μ [∇Uf + (∇Uf )T ]} ds
(7)
Table 1 Chemical reaction rate and heat parameters of carbon. Reaction n
where n is a normal unit vector; τ is the stress tensor on the particle surface, the unit is: N/m2; I is a unit vector on the pressure action direction; and ds is the differential area on the particle surface, the unit is: m2. In order to consider the effect of particle deformation, the soft sphere hypothesis is adopted. According to Glowinski et al. [21], the specific counting equations of Fcoll (on particle i imposed by particle j) is:
Fcoll = d ij kR R c 2
qṅ
An
En
(I) 2C + O2 → 2CO
1.97 × 107
1.98 × 105
2.21 × 105
(II) C + CO2 → 2CO
1.29 × 105
1.91 × 105
− 1.73 × 105
(III) 2CO + O2 → 2CO2
1.30 × 108
1.26 × 105
5.56 × 105
Table 2 Chemical reaction rate and heat parameters of NaCl.
(8)
Reaction n
An
En
qṅ
(IV) NaCl (s ) → NaCl (g )
6.32 × 109
2.81 × 105
− 181.42 × 105
where d ij is the sphere center vector of two particles, the unit is: m, and kR is a repulsive factor, defined as:
kR = kC hd
ρs ρf
(9)
where kC is 0.2–2; h is a characteristic scale related to the mesh generation method; and d is the problem dimension, given as 2 or 3. R c is defined as:
R c = Ri + Rj + lc − |d ij|
(10)
where Ri is the radius of particle i; Rj is the radius of particle j; lc is a repulsive range; and lc = h is satisfactory where particles do not overlap. Their units are all: m. 2.1.3. Chemical reaction model Only irreversible reactions described by chemical kinetics are considered, since their accuracy has been proved repeatedly [28,34]. Assuming N reactions and M species exist in the system, the relation between reactants and products is given by:
∑m ∈M vmn(r ) = ∑m ∈M vmn(p)
(11)
where v is the stoichiometric number. According to the basic theory of chemical reaction kinetics, the production rate of species m , ωṁ (the unit is: m3/s) in fluid can be described by:
ωṁ =
1 Wm ̇ C 1 Wm ωm = ∙ ∙ 1000 ρf 1000 ρf
N
∑n=1 vmn(p) kn ∏m ∈M
[m]vmn
Fig. 1. Computational domains [28].
(r )
(12) ̇ C where Wm is the molecular weight, the unit is: g/mol; ωm is the molar production rate, the unit is: mol/s; [m] is the molar concentration of specie m , the unit is: mol/m3; and kn is the reaction rate constant of reaction n , which can be evaluated by the Arrhenius equation:
E kn = An exp ⎛− n ⎞ ⎝ RT ⎠
(13)
where An is pre-exponential factor of reaction n , whose unit differs according to specific reaction, and En is activation energy of reaction n , the unit is: J/mol. The heat production rate for total reactions, q̇ (the unit is: J/s) is given by: M ̇ C= q ̇ = − ∑ hm ωm m
N
∑n
qṅ
(14)
Fig. 2. Validation of carbon burning rates for different temperatures.
where hm is the enthalpy of species m , and qṅ is the heat production rate of reaction n . Their units are all: J/mol.
According to the law of conservation of mass, as the mass of carbon in the particle decreases, its density variation can be described by:
2.1.3.1. Carbon transformation. The carbon transformation can be summarized by the heterogeneous reactions of C and O2 to produce CO and C and CO2 to produce CO on particle surface, and by the homogeneous reaction of CO and O2 to produce CO2 in the fluid. Table 1 lists the correlative kinetic parameters for chemical reactions.
∂ρs Wc vcn (r ) WC = −∑ ρ ωṁ = − ∑ vcn (r ) ωṅ C n = I , II W v (p) f n = I , II 1000 ∂t m mn
(15)
where Wc is the molecular weight of carbon; and ωṅ C is the molar 3
Fuel 255 (2019) 115763
S. Qu and C. You
where cp, f is the heat capacity of the local field; λT , f is thermal conductivity of the local field; ST , the fictitious domain source term, adopting different forms for different conditions is discussed below. For solid particle, the specific equation for solid particle temperature, Ts , is:
ρs Vs cp, s
∂Ts =− ∂t
∫ ρf cp,f Tf Uf ·nds + ∫ λT,s ∇Tf ·nds + q̇
(18)
where cp, s is heat capacity and λT , s is the thermal conductivity of the solid particle. 2.2. Coupling method for fluid and solid particle Governing differential equations for fluid and solid particle are very different in form. And complicated coupling relations among physical quantities in space and time make the problem lack analytical solutions. Therefore, a special strategy needs to be designed to satisfy the discrete solving logic of computer and keep the calculation results of two phases synchronized. According to fictitious domain method, a numerical solving strategy for multi-physical field of reacting fluid–solid flow is discussed in Section 2.2.1. The method to discrete differential equations into formatted numerical matrices is discussed in Section 2.2.2. Mesh description of computational domain is discussed in Section 2.2.3. And systematic calculation steps are summarized in Section 2.2.4.
Fig. 3. Validation of NaCl release characteristics for temperature increase.
production rate of reaction n . 2.1.3.2. Alkali metal transformation. Current assumption is that Na is primarily released in the form of NaCl gas, which mainly derives from NaCl solid [31]. The gasification of NaCl, was primarily determined by chemical reaction kinetics and Yang [34] confirmed the accuracy of this method. Table 2 lists the correlative kinetic parameters.
2.2.1. Unified equations form The main idea of fictitious domain method is that, suppose the region occupied by solid particle is also filled with fluid in flow field, as fictitious domain, and then a fictitious force is applied to the fluid in it, so that the fluid motion behavior satisfies the constraints of the rigid body motion. The specific strategy is to add different forms of the source term, named fictitious source terms, into the fluid governing equations in the region occupied by solid particle, so that the numerical results are the same as those obtained by the solid particle dynamics equations. For the momentum equation of fluid, i.e., Eq. (3), to realize the description of particle movement, the specific form of S can be derived by:
2.1.4. Species conservation equation The variation of gas species distribution in flow field can be described by the Fick's diffusion law. The specific volume fraction equation of specie m (Ym ) is:
∂ (Ym) 1 + Uf ∙∇Ym = Dm ∇2 Ym + ∇ρf Dm ∙∇Ym + ωṁ + Sm ∂t ρf
(16) 2
where Dm is the diffusion coefficient of specie Ym , the unit is: m/s , and Sm , the fictitious domain source term, adopting different forms for different conditions is discussed below. 2.1.5. Energy conservation equation The variation of temperature distribution in flow field can be described by the basic heat transfer theory. Coupling the fluid dynamics equations, the specific equation for fluid temperature, Tf , is:
∂Tf + Uf ·∇Tf ⎞ = λT , f ∇2 Tf + ∇λT , f ·∇Tf + q ̇ + ST ρf cp, f ⎛ ⎝ ∂t ⎠ ⎜
S = ρs / ρf ∙1/ Vs [Ff + |ωs
∑ Fcoll + (ρs − ρf ) Vs g] + ρs Tf /Is − ρs
× r + ∇p − ∇∙μ [∇Uf + (∇Uf )T − 2/3∇∙Uf ] − ρf g
(19)
where, ρs / ρf is density correction of solid particle and fluid, 1/ Vs [Ff + ∑ Fcoll + (ρs − ρf ) Vs g ] + ρs Tf / Is − ρs |ωs |2 × r is used to add kinetic information of solid particle, while
⎟
(a) Temperature (K) field
|2
(17)
(b) NaCl mass fraction (%) field
(c) CO mass fraction (%) field
Fig. 4. Temperature, NaCl mass fraction and CO mass fraction at t = 0.02 s. 4
Fuel 255 (2019) 115763
S. Qu and C. You
(a) Temperature (K) field
(b) NaCl mass fraction (%) field
(c) CO mass fraction (%) field
Fig. 5. Temperature, NaCl mass fraction and CO mass fraction at t = 0.15 s.
(a) Temperature (K) field
(b) NaCl mass fraction (%)
(c) CO mass fraction (%)
Fig. 6. Temperature, NaCl mass fraction and CO mass fraction at t = 0.45 s.
(a) Temperature (K) field
(b) NaCl mass fraction (%)
(c) CO mass fraction (%)
Fig. 7. Temperature, NaCl mass fraction and CO mass fraction at t = 0.49 s.
5
Fuel 255 (2019) 115763
S. Qu and C. You
Sm =
Wm C ωm, sur ρf
(20-b)
C
where, ωm, sur can be obtained by the surface reaction of the solid particle, by Eqs. (11)–(14). Here, fictitious source term, Sm , is used to add the products of heterogeneous reaction to flow field. For the energy conservation equation, i.e., Eq. (18), to realize a description of particle temperature variation, the specific form of ST can be derived by:
ST
( ∫ρ c
= −
f p, f
+ Uf ·∇Tf
Tf Uf ·nds + ρf cp, f
)ρVc
s s p, s
Fig. 8. Variation of temperature and mass residual rate of NaCl.
2.2.2. Spatial and temporal discrete solving method Due to velocity, Uf and pressure, p of flow field are coupled by the momentum equation i.e., Eq. (3), it is difficult to solve them independently. Therefore, it is necessary to design a numerical solution method. According to the Semi-Implicit Method for Pressure Linked Equations (SIMPLE), firstly proposed by SV Patankar, 1972, solving method used in this work iterates the velocity, Uf and pressure, p between the continuity equation i.e., Eq. (2) and the momentum equation i.e., Eq. (3) in discrete form repeatedly, so that the numerical values of velocity, Uf and pressure, p are more closer to the analytical solution after each iteration. Mole fraction of specie m , Ym and temperature, T of flow field are independently described by species conservation equation i.e., Eq. (16) and energy conservation equation i.e., Eq. (17), so that they can be solved respectively. Here, finite volume method (FVM) is used to discretize the differential equations, based on its advantage for satisfying integral conservation of each physical quantity in whole computational domain.
∇p − ∇∙μ [∇Uf + (∇Uf )T − 2/3∇∙Uf ] − ρf g is used to delete kinetic information of fluid to the equation. Thus, the original fluid equation can be used to calculate the solid particle movement by fictitious source term, S . For the species conservation equation of fluid, i.e., Eq. (16), to realize the description of species behavior in the region occupied by solid particle, assuming no reaction within the particles, the specific form of Sm can be derived by: 1 ∇ρ Dm ·∇Ym + ω̇m ) ρf f
(21)
where (ρf cp, f )/(ρs Vs cp, s ) is heat capacity correction of solid particle and fluid, − ∫ ρf cp, f Tf Uf ·nds + ∫ λT , s ∇Tf ·nds is used to add heat information of solid particle, while − λT , f ∇2 Tf − ∇λT , f ·∇Tf + Uf ·∇Tf is used to delete heat information of fluid to the equation. Thus, the original fluid equation can be used to calculate the particle temperature by fictitious source term, ST . Although the specific form of the governing equations is different from the original one after adding fictitious source terms, the numerical results are theoretically consistent, that is, no information of solid particles and fluids is lost. Due to it is not necessary to deduce governing equation frequently in the boundary between fluid and solid particle, effected by the strong interaction, fictitious domain strategy has obvious advantages in solving fluid-solid flow problems. As a result, it improves efficiency and stability in solving key parameters, such as velocity in whole computational domain.
Fig. 9. Particles arrangement at t = 0 s.
Sm = Uf ·∇Ym − (Dm ∇2 Ym +
∫ λT ,s ∇Tf ·nds − λT ,f ∇2Tf − ∇λT ,f ·∇Tf
2.2.3. Description of solid particle and fluid After establishing unified equations in discrete form, structural equilibrium meshes can be used to describe both solid particle and fluid. Fig. 1 shows three types of computational domains: region only occupied by fluid (only red triangles); region only occupied by solid particle (only blue rectangles); and the region partially occupied by solid particles and fluid (both red triangles and blue rectangles). The background flow field comprises Euler meshes, which are applicable to the governing fluid equations (or unified equations without fictitious source terms). The field (Euler meshes) totally covered by solid particles (solid) is applicable to the governing rigid body equations (or unified equations with fictitious source terms). During numerical calculation, Euler staggered meshes are used to store information such as velocity, Uf and pressure, p of flow field. Representing particle surface, the boundary between fluid and solid particle comprises Lagrangian points (as shown by black nodes in
(20-a)
By adding fictitious source term, Sm , D (Ym)/ Dt = ∂ (Ym)/ ∂t + Uf ∙∇Ym is 0 in Eq. (16), since there is no component distribution in solid particle under consideration. While in the region near the solid particle boundary (on the surface of the solid particle), the specific form of Sm can be derived by: 6
Fuel 255 (2019) 115763
S. Qu and C. You
(a) Temperature (K) field
(b) NaCl mass fraction (%) field
(c) CO mass fraction (%) field
Fig. 10. Temperature, NaCl mass fraction and CO mass fraction at t = 0.05 s.
(a) Temperature (K) field
(b) NaCl mass fraction (%) field
(c) CO mass fraction (%) field
Fig. 11. Temperature, NaCl mass fraction and CO mass fraction at t = 0.1 s.
physical field. The classical SIMPLE can be used to solve velocity, Uf and pressure, p of pure flow field well. On this basis, strategy adopts to fictitious source terms to modify the governing equation forms of fluid, so that the behavior of solid particle and fluid in whole computational domain can be solved in a unified way. In addition, FVM, specific meshes construction, introducing more equations, and etc. make FRS more authentic in physics. The solving process is summarized as follows:
Fig. 1) that distribute Lagrangian source terms into Eulerian meshes outside the solid particle. The background meshes, including Lagrangian points, is the region influenced by both particles and fluid. This region can also be calculated by the unified equations with fictitious source terms adding modified parameter, φ , such that:
φ = Vs, E / Vf , E
(22)
where Vs, E is the particle volume in the Eulerian mesh and Vf , E is the fluid volume in the same Eulerian mesh. The value of modified parameter, φ ranges from 0 to 1.
A. Initialize physical parameters and set structured meshes for the computational domain at the beginning of the n-th time loop; B. Calculate the force, Fs and torque, Ts on solid particle by Eqs. (6)–(10), then solve Eqs. (4) and (5) to obtain the velocity, Us , position vector, Ss and acceleration velocity, as of solid particle; C. Divide the solid particle regions and fluid regions in whole
2.2.4. Solving process For FRS of reacting fluid–solid flow, the difficulty lies in the description and solution of the strong interaction between two phases and complicated coupling relations among physical quantities in the multi7
Fuel 255 (2019) 115763
S. Qu and C. You
(a) Temperature (K) field
(b) NaCl mass fraction (%) field
(c) CO mass fraction (%) field
Fig. 12. Temperature, NaCl mass fraction and CO mass fraction at t = 0.15 s.
(a) Temperature (K) field
(b) NaCl mass fraction (%) field
(c) CO mass fraction (%) field
Fig. 13. Temperature, NaCl mass fraction and CO mass fraction at t = 0.16 s.
D. E. F.
G.
consume) mass and heat, by Eqs. (11)–(14) and correct the physical parameters of solid particle and fluid, by thermodynamic equations and the mass conservation equation, e.g., Eq. (15); H. Analyze the specific fictitious source terms, Sm and ST by Eqs. (20-a), Eq. (20-b) and Eq. (21) and substitute them into Eqs. (16) and (17); I. Discretize Eqs. (16) and (17) by FVM and solve them to obtain the temperature, Tf and species, Ym distribution of flow field; J. Judge that, if the iteration results of key parameters such as velocity, Uf and pressure, p of flow field converge, terminate the computational iteration of SIMPLE and enter step. K. Otherwise, return to
computational domain. For the meshes including the boundary between solid particle and fluid, calculate the modified parameter, φ by Eq. (22); Calculate the temperature, Ts of solid particle by Eq. (18); Analyze the specific fictitious source term, S by Eq. (19) and substitute it into Eq. (3); Start the computational iteration of SIMPLE, discretize and solve Eqs. (2) and (3) to obtain the coupled solution of velocity, Uf and pressure, p distribution of flow field; Calculate chemical reaction source terms which produce (or 8
Fuel 255 (2019) 115763
S. Qu and C. You
(a) Temperature (K) field
(b) NaCl mass fraction (%) field
(c) CO mass fraction (%) field
Fig. 14. Temperature, NaCl mass fraction and CO mass fraction at t = 0.17 s.
(a) Temperature (K) field
(b) NaCl mass fraction (%) field
(c) CO mass fraction (%) field
Fig. 15. Temperature, NaCl mass fraction and CO mass fraction at t = 0.2 s.
coupled with fictitious domain strategy, into SIMPLE, the comprehensive algorithm system is established.
step. F; K. Update the physical parameters of solid particle and fluid by the physical equations, e.g., Eq. (1) and return to Step. A to commence the (n + 1)-th time loop, until the end of a given physical time.
3. Calculation results and analyses
Step. F to J represent the whole calculating process of solving flow field by SIMPLE. Step. C, E and H embody fictitious domain strategy of unified treatment of solid particle and fluid. Step. I represents the calculating of the temperature, Tf and species, Ym equations with integral discrete scheme by FVM, based on the conservation law of mass and energy. By embedding the discrete equations solving with FVM,
Gasification behavior of alkali metals coupled with carbon combustion is the key to this problem. Therefore, carbon combustion model and Na gasification model are validated at first in Section 3.1, to prove their accuracy. Then, discrete particle combustion and multiple particles combustion in air are simulated respectively in Sections 3.2 and 3.3. 9
Fuel 255 (2019) 115763
S. Qu and C. You
boundary condition of the computational domain is open and the other three sides are walls. It can be observed from the calculation results that particle starts to burn at t = 0 from which point there is a rapid temperature increase. Figs. 4 and 5 show that heating accelerates chemical reactions from t = 0.02 to 0.15 s. Fig. 5 shows that combustion stabilizes at t = 0.15 s. The combustion rate of carbon and the production rate of gas remain unchanged until the combustion reaches its final stage. Fig. 6 shows that the combustion rate of carbon decreases after t = 0.45 s. Fig. 6(a) shows that the particle temperature rapidly decreases due to air flow and insufficient heat supply from the combustion reactions, which leads to a rapid decrease in Na release, as shown in Fig. 6(c). However, there is still a considerable amount of Na in the particle. Fig. 7 shows that the carbon combustion process completes at about t = 0.49 s. Fig. 7(a) shows that the remaining CO is gradually removed by the air flow. Fig. 7(c) shows that there is a low Na release rate in the higher temperature region. In addition, Fig. 8 shows more obviously that, in the initial stage (approximately 10% of the whole combustion), temperature rises rapidly and Na release accelerates. During the stable stage (approximately 85% of the whole combustion), temperature is constant in a range slightly higher than that of ambient medium. At this time, combustion rate of particle is unchanged as a stable heat source, while Na release slows down slightly, but is not obvious. This is mainly due to the continuous decrease of the absolute mass of Na. At the end (approximately 5% of the whole combustion), particle stops the oxidation and gasification, and no heat is generated continuously, which leads to the rapid reduction of temperature, and Na release decreases rapidly, although more than 50% of it remains in particle. Therefore, it can be concluded that Na release depends heavily on the local high temperature, which is consistent with the phenomenon obtained by experimental investigations that the chemical transformation reactions of Na are closely related to temperature [35–38]. Na release mainly occurs during the stable combustion stage due to the need for a stable heat source to provide a sustained high temperature condition.
Fig. 16. Comparison of NaCl residue variation for a single particle against the average value of 150 particles.
3.3. Multiple particles combustion Fig. 17. Frequency distribution of NaCl residue for 150 particles at 0.1 s and 0.2 s.
Corresponding to dense particles combustion, such as in the area near the boiler burner, or in the dense flow zone of fluidized bed, multiple particles combustion moving on the background meshes is simulated. Especially when the particles volume fraction is large, there are frequent collisions between particles, and between particles and walls. These collisions significantly affect the molecular transport of gas products and thus, affect the specific combustion state of particles. Figs. 10–15 shows the calculation results for a multiple particulate system where particles are moving, colliding, burning and interacting with fluid. The particle has a diameter of 1 mm, density of 1000 kg/m3 and the mass fraction of Na is 0.5%. The ambient pressure is 1.01 × 105 Pa; the temperature is 1273 K; and the velocity of air flow from the bottom of the computational domain (0.35 × 0.1 m) is 2 m/s. The boundary conditions are identical to the previous cases. At the beginning of the calculation (i.e., t = 0, as shown in Fig. 12), the group of 150 (15 × 10) particles evenly arranged at the top of the computational domain are released. Fig. 10 shows that at t = 0.05 s, the incoming flow begins to disrupt the structure of the multiple particulate system. There is a gradual heat increase on the windward side of the particles and reaction I, II, IV commence on the surface of the particles; reaction III primarily occurs at the tail region of the particles. Fig. 11 shows that at t = 0.1 s, the particles settle in the middle of the computational domain. Particles collide on the windward side and small vortexs separate from the wall in the tail region of the particles. Fig. 12 shows that at t = 0.15 s, some particles reach the bottom and begin to bounce. Here, the motion patterns of the particles and fluid in the computational domain maintain a degree of symmetry. Fig. 13 shows that at t = 0.16 s, all particles
3.1. Model validation Fig. 2 shows a validation of char combustion model. Simulation results are obtained by taking the average value of the carbon mass release rate in the stable burning stage (i.e. linear variation is satisfied). The results here are compared with those from the Makino’s experiment [32], which comprised a graphite rod of 5 mm diameter to establish a two-dimensional stagnation flow combustion under an ambient temperature of 320 K with an air flow velocity of 4.125 m/s. Errors from the current simulation are within the acceptable range. Fig. 3 shows a validation of NaCl release model, established according to Yang [34]. Simulation results are obtained by calculating the release ratio of the NaCl mass, accounting for 68% of Na content in Zhundong coal, as the temperature increases from 673 to 1473 K. The model-based results are in agreement with the experimental data. 3.2. Single particle combustion Corresponding to sparse particles combustion, such as pulverized coal boiler combustion, a single particle combustion fixed in the background meshes is simulated. Figs. 4–7 shows the calculation results. The particle has a diameter of 50 × 10−6 m, density of 1000 kg/ m3, and the mass fraction of Na is 0.5%. Ambient pressure is 1.01 × 105 Pa; temperature is 1273 K; and the velocity of air flow from the bottom of the computational domain (0.002 × 0.001 m) is 1 m/s. Top 10
Fuel 255 (2019) 115763
S. Qu and C. You
produces a large number of vortexes, which seriously affects the transport behavior of Na in gas phase. Especially for the large vortex, the entrainment makes Na in gas phase locally enriched in the flow field.
reach the bottom and the volume fraction of the particles achieves maximum compression. Frequent collisions occur in the particle group and, thus, the flow field and particles no longer have symmetry. Fig. 14 shows that at t = 0.17 s, the particles diffuse back and collide frequently, which results in chaotic movement both for the fluid and solid particles. Also, large vortexes are more likely to occur. It is necessary to highlight here that where the large vortex exists, Na in gas phase is enriched. Fig. 15 shows that at t = 0.2 s, the rebounding particles are distributed across the entire flow field. The whole simulation shows that the fluid movement has a significant impact on the molecular transport of gas products, especially the emergence of vortexes, which promote the formation of local high gas concentration due to entrainment. As shown in Fig. 16, comparing Na residue of a single particle with the average Na residues of dense particles under the same combustion conditions and with the same particle characteristics, it can be found that the average Na release of dense particles is slightly larger than that of a single particle in the initial period (< 0.05 s) of combustion. As combustion proceeds, the release rate of dense particles slows down gradually, while Na release rate of a single particle does not change significantly. The preliminary guess is that, at the beginning of combustion, the oxygen is sufficient, and particles on the windward side begin to burn, which makes the heat accumulated in particles more than that carried away by the fluid, which results in local high temperature to make Na release rate at a very high level, and even makes the average release rate of all 150 particles is higher than that of a single particle. As combustion proceeds, oxygen is rapidly consumed and even insufficient to supply, which limits the local combustion rate, prevents the local high temperature of particles from continuing, and leads to the decrease of the average release rate of dense particles. The gap of Na release rate between dense particles and a single particle is larger. As shown in Fig. 17, from t = 0.1 to 0.2 s, Na residue of particles generally decreases. Further, the centrality of the Na release rate decreases and the peak value of particle Na residue surpasses 80% at 0.1 s, then decreases to less than 75% at 0.2 s. Very few particles, guessed to be located on the windward side, lose most of their Na in a short time. It can be predicted that, as combustion proceeds, the dense particulate system becomes chaotic and the peak value of the Na residue of particles shifts leftwards. The peak distribution is more uniform, either.
In this work, the feasibility of this TURE DNS method was confirmed. Its potential can be further exploited for numerically investigating the characteristics of alkali metals release and transfer in multi-physical fields in more details. Acknowledgement This work was financially supported by the National Natural Science Foundation of China (51761125012 and 51390493). References [1] Zhou JB, Zhuang XG, Alastuey A, Querol X, Li JH. Geochemistry and mineralogy of coal in the recently explored Zhundong large coal field in the Junggar basin, Xinjiang province, China. Int J Coal Geol 2010;82:51–67. [2] Yang ZC, Liu JL, He HG. Study on properties of Zhundong coal in Xinjiang region and type-selection for boilers burning this coal sort. Thermal Power Gener Chin 2010;39:38–44. [3] Sugawara K, Enda Y, Inoue H, Sugawara T, Shirai M. Dynamic behavior of trace elements during pyrolysis of coals. Fuel 2002;81:1439–43. [4] Li GD, Li SQ, Qian H, Qiang Y. Fine particulate formation and ash deposition during pulverized coal combustion of high-sodium lignite in a down-fired furnace. Fuel 2015;143:430–7. [5] Takuwa T, Naruse I. Emission control of sodium compounds and their formation mechanisms during coal combustion. Proc Combust Inst 2007;31:2863–70. [6] Takuwa T, Mkilaha ISN, Naruse I. Mechanisms of fine particulates formation with alkali metal compounds during coal combustion. Fuel 2006;85:671–8. [7] Naruse I, Kamihashira D, Khairil, Miyauchi Y, Kato Y, Yamashita T, Tominaga H. Fundamental ash deposition characteristics in pulverized coal reaction under high temperature conditions. Fuel 2005;84:405–10. [8] Benson SA, Holm PL. Comparison of inorganics in three low-rank coals. Ind Eng Chem Prod Res Dev 1985;24:145–9. [39] Wang CA, Jin X, Wang YK, Yan Y, Cui J, Liu YH, Che DF. Release and transformation of sodium during pyrolysis of zhundong coals. Energy Fuels 2015;29:78–85. [10] He Y, Zhu J, Li B, Wang ZH, Li ZS, Alden M, et al. In-situ measurement of sodium and potassium release during oxy-fuel combustion of lignite using laser-induced breakdown spectroscopy: effects of O2 and CO2 concentration. Energy Fuels 2013;27:1123–30. [11] Wang ZH, Liu YZ, Ronald W, Wan KD, He Y, Xia J, et al. Measurement of atomic sodium release during pyrolysis and combustion of sodium-enriched Zhundong coal pellet. Combust Flame 2017;176:429–38. [12] Liu YZ, He Y, Wang ZH, Wan KD, Xia J, Liu JZ, et al. Multi-point LIBS measurement and kinetics modeling of sodium release from a burning Zhundong coal particle. Combust Flame 2018;189:77–86. [13] Van Eyk PJ, Ashman PJ, Nathan GJ. Mechanism and kinetics of sodium release from brown coal char particles during combustion. Combust Flame 2011;158:2512–23. [14] Wan KD, Vervisch L, Xia J, Domingo P, Wang ZH, Liu YZ, et al. Alkali metal emissions in an early-stage pulverized-coal flame: DNS analysis of reacting layers and chemistry tabulation. Proc Combust Inst 2019;37:2791–9. [15] Fox RW, McDonald AT, Pritchard PJ. Introduction to fluid dynamics. New York: John Wiley and Sons; 1985. p. 354. [16] Hu HH, Patankar NA, Zhu MY. Direct numerical simulations of fluid–solid systems using the arbitrary Lagrangian-Eulerian technique. J Comput Phys 2001;169:427–62. [17] Lomtev I, Kirby RM, Karniadakis GE. A discontinuous Galerkin ALE method for compressible viscous flows in moving domains. J Comput Phys 1999;155:128–59. [18] Osher S, Fedkiw R. Level set methods and dynamic implicit surfaces. Surfaces 2002;44:77. [19] Sierakowski AJ, Andrea P. Resolved-particle simulation by the Physalis method: enhancements and new capabilities. J Comput Phys 2016;309:164–84. [20] Uhlmann M. An immersed boundary method with direct forcing for the simulation of particulate flows. J Comput Phys 2005;209:448–76. [21] Glowinski R, Pan TW, Hesla TI, Joseph DD. A distributed Lagrange multiplier/fictitious domain method for particulate flows. Int J Multiph Flow 1999;25:755–94. [22] Peskin CS. The immersed boundary method. Acta Numerica 2002;11. [23] Fogelson AL, Peskin CS. A fast numerical method for solving the three-dimensional stokes' equations in the presence of suspended particles. J Comput Phys 1988;79:50–69. [24] Lai MC, Peskin CS. An immersed boundary method with formal second-order accuracy and reduced numerical viscosity. J Comput Phys 2000;160:705–19. [25] Apte SV, Martin M, Patankar NA. A numerical method for fully resolved simulation (FRS) of rigid particle–flow interactions in complex flows. J Comput Phys 2009;228:2712–38. [26] Ghasemi A, Pathak A, Raessi M. Computational simulation of the interactions
4. Conclusion A TURE DNS method, based on fictitious domain strategy with finite volume (FV) scheme, was developed to describe the behavior of alkali metals release and transfer during char combustion. Coupling fluid dynamics equations, particles dynamics equations, chemical reaction models, species conservation equations and energy equations, the introduction of SIMPLE enabled the algorithm to perform better in terms of accuracy and computational efficiency. After verifying the method against existed experimental results, a single fixed particle combustion simulation corresponding to sparse particles combustion, and multiple particles combustion, corresponding to dense particles combustion were carried out respectively. The conclusions are summarized as follows. A. Na release depends heavily on the local high temperature in flow field, and mainly occurs during the stable combustion stage (approximately 85% of the whole combustion) due to the need for a stable heat source to provide a sustained high temperature condition. B. Due to rapid heat accumulation with sufficient oxygen, the average Na release rate of dense particles is higher than that of a single particle at the beginning of combustion. With oxygen consumes, average Na release rate of dense particles decreases gradually. C. Frequent collision of particles makes the flow field disordered and 11
Fuel 255 (2019) 115763
S. Qu and C. You
[27] [28] [29] [30] [31] [32] [33]
between moving rigid bodies and incompressible two-fluid flows. Comput Fluids 2014;94:1–13. Haeri S, Shrimpton JS. A new implicit fictitious domain method for the simulation of flow in complex geometries with heat transfer. J Comput Phys 2013;237:21–45. Zhang LH, Liu K, You CF. Fictitious domain method for fully resolved reacting gassolid flow simulation. J Comput Phys 2015;299:215–28. Patankar S. Numerical heat transfer and fluid flow. CRC Press; 1980. Yu ZS, Shao XM. Direct numerical simulation of particulate flows with a fictitious domain method. Int J Multiph Flow 2010;36:127–34. Lee SHD, Teats FG, Swiift WM, Banerjee DD. SHORT COMMUNICATION alkalivapor emission from PFBC of illinois coals. Combust Sci Technol 1992;86:327–36. Makino A, Namikiri T, Kimura K. Combustion rates of graphite rods in the forward stagnation field with high-temperature airflow. Combust Flame 2003;132:743–53. Yang YM, Zhang H, Wu YX, Lv JF. Release of alkali/alkaline earth metal species in zhundong coal at different ashing temperatures. J Combust Sci Technol Chin
2015;21:297–300. [34] Yang YM. Occurrence modes and release characteristics of alkali and alkaline earth metal species in zhundong coals. Beijing: Tsinghua University; 2017. [35] Kosminski A, Ross DP, Agnew JB. Transformations of sodium during gasification of low-rank coal. Fuel Process Technol 2006;87:943–52. [36] Li WQ, Wang LY, Qiao Y, Lin JY, Wang MJ, Chang LP. Effect of atmosphere on the release behavior of alkali and alkaline earth metals during coal oxy-fuel combustion. Fuel 2015;139:164–70. [37] Song GL, Song WJ, Qi XB, Lu QG. Transformation characteristics of sodium of zhundong coal combustion/gasification in circulating fluidized bed. Energy Fuels 2016;30:3473–8. [38] Wang ZS, Wang LY, Lin JY, Wang MJ, Chang LP. The influence of the addition of sodium on the transformation of alkali and alkaline-earth metals during oxy-fuel combustion. J Energy Inst 2018;91:502–12.
12
Contents lists available at ScienceDirect
Fuel journal homepage: www.elsevier.com/locate/fuel
Full Length Article
Direct numerical simulation (DNS) of alkali metals released during char combustion Sibo Qu, Changfu You
T
⁎
Key Laboratory for Thermal Science and Power Engineering of the Ministry of Education, Tsinghua University – University of Waterloo Joint Research Centre for Micro/ Nano Energy & Environment Technologies, Department of Energy and Power Engineering, Tsinghua University, Beijing 100084, China
A R T I C LE I N FO
A B S T R A C T
Keywords: Alkali metals release Char combustion Direct numerical simulation (DNS) Sodium chloride (NaCl)
Alkali metals in pulverized coal easily released in gas phase, cause significant fouling, atmospheric pollution and etc. during combustion and gasification. A TURE DNS method, based on fictitious domain strategy with finite volume (FV) scheme, was developed to describe the behavior of alkali metals release and transport in reacting fluid–solid flow at sub-particle level. Taking water-soluble sodium chloride (NaCl) of high harmfulness as an example, alkali metals release and transport in physical field of char combustion, as the main stage of coal combustion, was simulated. A 50 × 10−6 m single particle, and 150 1 × 10−3 m closely spaced particles combustion at 1273 K air environment were respectively calculated. Computational results revealed that continuous release of alkali metals depends heavily on the local high temperature from the stable combustion stage (approximately 85% of the whole combustion stage) which produces stable heat. And complicated movement of fluid medium has a great effect on the transport behavior of alkali metals in gas phase, especially the entrainment of vortexes, which results in the local enrichment of alkali metals in gas phase.
1. Introduction Alkali-rich coals, e.g., Zhundong Coal from Xinjiang, China, because of its huge deposit, high quality characterized by low ash, aluminum, iron, sulfur, and trace element contents [1], has wide prospects in industrial application and scientific research. However, the high content of alkali metals is the main limitation. Existed researches showed that Na content in ash is generally more than 5%, which is significantly more than the level in typical Chinese coal, ranging from 1 to 2 % [2]. Alkali metals have the strongest sublimation tendency, for the main mineral content [3] during combustion and gasification. Especially Na [4], is generally accepted as the primary cause of serious ash-related problems, such as fouling, slagging and corrosion, for many industrial systems [5]. Further, alkali metals contribute greatly to the formation of atmospheric pollutants, such as fine particulates (size less than 1 µm) due to homogeneous nucleation or heterogenous condensation [5–7]. A chemical fractionation method of successive leaching with H2O, CH3COONH4 and HCl solutions has been widely used to quantify alkali metals using terms such as, “H2O soluble”, “CH3COONH4 soluble”, “HCl soluble”, and “insoluble form” [8]. Among them, H2O solubles such as alkali metal chlorides [33], should be given a priority consideration, due to their proven harmfulness, high content and easy conversion into gas phase [9]. Taking Na as an example, experimental investigations on ⁎
its release and reacting dynamics evolved from offline to online measuring techniques using advanced laser diagnostics [10–12]. On the numerical side, a one-step Arrhenius model firstly proposed by van Eyk et al. [13] was used to study the Na release during char combustion. A two-step kinetics model for all stages of coal combustion was developed by Liu et al. [12] then, based on the Na release, surface temperature and diameter of particles. These models can quantitatively describe the evolution of Na in solid fuels, but can not fully deal with the multiphysical fields in real fluid-solid combustion. So far, few researches have been done on kinetic behavior of alkali metals in multi-physical fields in detail after they release from coal in gas phase. A method using one-dimensional (1D) premixed flamelet generated manifolds (FGM), was used to study Na emissions in an early-stage pulverized coal flame [14]. But, the solid particle was treated as Lagrangian point source, internal structure and boundary layers of which were not resolved. So that it can only be called “quasi- DNS” or “point-particle DNS”, since it lost physical integrity. Therefore, based on existed results, a TURE DNS method, was developed in this work, to describe the behavior of alkali metals release and transport in reacting fluid–solid flow at sub-particle level in detail, further, to study the characteristics of alkali metals during combustion of solid fuels, and to provide a guide for industrial design of high-efficiency clean combustion technology. Based on fluid dynamics [15], mass and heat transfer theory [29],
Corresponding author. Tel.:/fax: +86 10 62785669. E-mail address: [email protected] (C. You).
https://doi.org/10.1016/j.fuel.2019.115763 Received 28 March 2019; Received in revised form 23 June 2019; Accepted 3 July 2019 Available online 31 July 2019 0016-2361/ © 2019 Elsevier Ltd. All rights reserved.
Fuel 255 (2019) 115763
S. Qu and C. You
2.1. Reacting fluid-solid flow model
numerical simulation is superior to the experiment in obtaining the details of multi-physical fields. There are two main types of the fully resolved simulation (FRS) method for fluid-solid flow system at subparticle level: body conformal mesh method [16] represented by Arbitrary Lagrange-Euler (ALE) [17] method, and fixed mesh method. Body conformal mesh method which constructs body conformal meshes outside solid particle, can describe the boundary condition of solid particle well, so that the computational result is more accurate. However, a moving particle leads to meshes need to be reconstructed constantly. Especially in the case of multiple particles, meshes reconstruction is frequent, which makes computational calculation tedious. And in particles collision, too small spacing between two particles can easily cause meshes distortion. Thus, it is not suitable to simulate the movement of a large number of particles with body conformal mesh method. Fixed mesh method performs well, because it does not require frequent meshes reconstruction, which mainly includes level set method [18], physalis method [19], immersed boundary method (IBM) [20], and fictitious domain method created by Glowinski et al. [21]. Level set method is mainly used to calculate the movement and deformation of the interface, thus, it is more suitable for the calculation of gas-liquid flow. Physalis method requires complicated local analytic solutions as a medium to connect the surface of solid particle to background fluid meshes. Thus, it lacks applicability for some cases without local analytical solutions. IBM is a commonly used method for FRS with fixed meshes. This method was proposed by Peskin [22] to study the flow around flexible membranes (mainly for heart). Subsequently, the method was applied to calculate the flow around suspended particles and fixed cylinders [23,24]. IBM uses a series of Lagrange points to describe the solid particle boundary, and a specific strategy to connect them with the Euler points on the background fluid meshes to realize the coupling of fluid-solid phases. Because of its better completeness in describing physical process, higher computational efficiency and stability, IBM has been widely used. The solution strategy of fictitious domain method developed by Apte [25], Ghasemi [26], Haeri [27] and Zhang [28] et al. is similar to that of IBM, which has also been widely used and developed [30]. The main difference between these two is that fictitious domain method not only deals with the boundary, but also needs to adopt a certain numerical method in interior of particle to ensure the rigid body motion law. Therefore, with the basic advantages of IBM, fictitious domain method is more dominant in describing reacting fluid–solid flow with a high density ratio. A TURE DNS method, based on fictitious domain strategy with finite volume scheme was developed. Unified structured meshes were used to represent both fluid and solid particle regions in computational domain. Here, the solid particle regions adapted to the unified equation forms by adding different source terms, named fictitious source terms. Taking water-soluble sodium chloride (NaCl) of high harmfulness as an example, alkali metals release and transport in physical field of char combustion, as the main stage of coal combustion, was simulated. Section 2 describes the algorithm development process, including governing physical equations and mathematical methods for coupled equation solutions. Section 3 validates the present model with existed experimental results at first, and the calculation results of temperature and component distribution obtained under the present method follow then, which prove feasibility, clarify the effects of local parameters, and preliminarily study the characteristics of alkali metals release and transfer in physical field of char combustion.
Physical descriptions for different phases can be divided into two types: fluid solver, and solid particle solver. Fluid solver mainly comprises applicable physical equations, dynamic equations, presented in Section 2.1.1, species conservation equation for mass transfer, presented in Section 2.1.4, and energy equation for heat transfer, presented in Section 2.1.5. Particle solver mainly comprises dynamics equations of Newton's laws and the rigid constraint, presented in Section 2.1.2, energy equation, presented in Section 2.1.5 either. Descriptions of homogeneous and heterogeneous reactions occurring in physical field of reacting fluid-solid flow are presented in Section 2.1.3, which are used as source terms. 2.1.1. Fluid dynamics equations The Navier-Stokes equations are used to describe fluid movement. Here, according to the characteristics of the fluid medium and closure of flow field, whether to adopt Navier-Stokes equations in compressible form or not can be considered, further. According to the ideal gas hypothesis, the state equation is:
ρf =
P RT
(1)
where ρf is fluid density, the unit is: kg/m ; P is pressure, the unit is: Pa; T is temperature, the unit is: K; and R is the specific gas constant of local fluid, usually the value is: 8.314 J/(mol·K). Eq. (1) is used to describe the basic relation among fluid density, ρf , pressure, P and temperature, T. According to the basic idea of taking fluid element as control body and the law of mass conservation, the continuity (mass conservation) equation of fluid is: 3
∂ρf ∂t
+ ∇∙ (ρf Uf ) = 0
(2)
where Uf is the local fluid velocity, the unit is: m/s; t is physical time, the unit is: s. According to the law of conservation of momentum, the momentum equation of fluid is:
∂ (ρf Uf ) ∂t
+ ∇∙ (ρf Uf ∙Uf )
= −∇p + ∇∙μ [∇Uf + (∇Uf )T − 2/3∇∙Uf ] + ρf g + S
(3)
where μ is fluid viscosity, the unit is: Pa·s; g is gravity acceleration, usually the value is: 9.8 m/s2; and S , the fictitious domain source term, adopting different forms for different conditions, is discussed below. 2.1.2. Solid particle dynamics equations Particle movement can be described by Newton's motion laws with specific consideration given to buoyancy effects, due to heat and frequent particle collision [28]. The correlative equations are:
∂Ut 1 = [Ff + ∂t ρs Vs
∑ Fcoll + (ρs − ρf ) Vsg]
∂ (ωs × r ) = ρs Tf / Is − |ωs |2 × r ∂t
(4)
(5)
where Ut is translational particle velocity, the unit is: m/s; ρs is particle density, the unit is: kg/m3; Vs is particle volume, the unit is: m3; Ff is force from fluid, the unit is: N; Fcoll is force from particle and wall collision, the unit is: N; ωs is particle angular velocity, the unit is: rad/s; r is vector position from particle center to focal point; Is is the moment of particle inertia, the unit is: kg·m2; and Tf is torque from fluid, the unit is: N·s. The specific counting equations of fluid force, Ff and fluid torque Tf are:
2. Numerical method To deal with complicated multi-physical field problem of reacting fluid-solid flow, a set of governing equations is established, based on the basic conservation relation and transfer theory of relevant physical quantities, in Section 2.1. And fictitious domain strategy to link particle variables with background fluid variables for realizing interaction is discussed in detail in Section 2.2. 2
Fuel 255 (2019) 115763
S. Qu and C. You
Ff =
∫ n·τds = ∫ n·{−pI + μ [∇Uf + (∇Uf )T ]} ds
(6)
Tf =
∫ r × τds = ∫ r × {−pI + μ [∇Uf + (∇Uf )T ]} ds
(7)
Table 1 Chemical reaction rate and heat parameters of carbon. Reaction n
where n is a normal unit vector; τ is the stress tensor on the particle surface, the unit is: N/m2; I is a unit vector on the pressure action direction; and ds is the differential area on the particle surface, the unit is: m2. In order to consider the effect of particle deformation, the soft sphere hypothesis is adopted. According to Glowinski et al. [21], the specific counting equations of Fcoll (on particle i imposed by particle j) is:
Fcoll = d ij kR R c 2
qṅ
An
En
(I) 2C + O2 → 2CO
1.97 × 107
1.98 × 105
2.21 × 105
(II) C + CO2 → 2CO
1.29 × 105
1.91 × 105
− 1.73 × 105
(III) 2CO + O2 → 2CO2
1.30 × 108
1.26 × 105
5.56 × 105
Table 2 Chemical reaction rate and heat parameters of NaCl.
(8)
Reaction n
An
En
qṅ
(IV) NaCl (s ) → NaCl (g )
6.32 × 109
2.81 × 105
− 181.42 × 105
where d ij is the sphere center vector of two particles, the unit is: m, and kR is a repulsive factor, defined as:
kR = kC hd
ρs ρf
(9)
where kC is 0.2–2; h is a characteristic scale related to the mesh generation method; and d is the problem dimension, given as 2 or 3. R c is defined as:
R c = Ri + Rj + lc − |d ij|
(10)
where Ri is the radius of particle i; Rj is the radius of particle j; lc is a repulsive range; and lc = h is satisfactory where particles do not overlap. Their units are all: m. 2.1.3. Chemical reaction model Only irreversible reactions described by chemical kinetics are considered, since their accuracy has been proved repeatedly [28,34]. Assuming N reactions and M species exist in the system, the relation between reactants and products is given by:
∑m ∈M vmn(r ) = ∑m ∈M vmn(p)
(11)
where v is the stoichiometric number. According to the basic theory of chemical reaction kinetics, the production rate of species m , ωṁ (the unit is: m3/s) in fluid can be described by:
ωṁ =
1 Wm ̇ C 1 Wm ωm = ∙ ∙ 1000 ρf 1000 ρf
N
∑n=1 vmn(p) kn ∏m ∈M
[m]vmn
Fig. 1. Computational domains [28].
(r )
(12) ̇ C where Wm is the molecular weight, the unit is: g/mol; ωm is the molar production rate, the unit is: mol/s; [m] is the molar concentration of specie m , the unit is: mol/m3; and kn is the reaction rate constant of reaction n , which can be evaluated by the Arrhenius equation:
E kn = An exp ⎛− n ⎞ ⎝ RT ⎠
(13)
where An is pre-exponential factor of reaction n , whose unit differs according to specific reaction, and En is activation energy of reaction n , the unit is: J/mol. The heat production rate for total reactions, q̇ (the unit is: J/s) is given by: M ̇ C= q ̇ = − ∑ hm ωm m
N
∑n
qṅ
(14)
Fig. 2. Validation of carbon burning rates for different temperatures.
where hm is the enthalpy of species m , and qṅ is the heat production rate of reaction n . Their units are all: J/mol.
According to the law of conservation of mass, as the mass of carbon in the particle decreases, its density variation can be described by:
2.1.3.1. Carbon transformation. The carbon transformation can be summarized by the heterogeneous reactions of C and O2 to produce CO and C and CO2 to produce CO on particle surface, and by the homogeneous reaction of CO and O2 to produce CO2 in the fluid. Table 1 lists the correlative kinetic parameters for chemical reactions.
∂ρs Wc vcn (r ) WC = −∑ ρ ωṁ = − ∑ vcn (r ) ωṅ C n = I , II W v (p) f n = I , II 1000 ∂t m mn
(15)
where Wc is the molecular weight of carbon; and ωṅ C is the molar 3
Fuel 255 (2019) 115763
S. Qu and C. You
where cp, f is the heat capacity of the local field; λT , f is thermal conductivity of the local field; ST , the fictitious domain source term, adopting different forms for different conditions is discussed below. For solid particle, the specific equation for solid particle temperature, Ts , is:
ρs Vs cp, s
∂Ts =− ∂t
∫ ρf cp,f Tf Uf ·nds + ∫ λT,s ∇Tf ·nds + q̇
(18)
where cp, s is heat capacity and λT , s is the thermal conductivity of the solid particle. 2.2. Coupling method for fluid and solid particle Governing differential equations for fluid and solid particle are very different in form. And complicated coupling relations among physical quantities in space and time make the problem lack analytical solutions. Therefore, a special strategy needs to be designed to satisfy the discrete solving logic of computer and keep the calculation results of two phases synchronized. According to fictitious domain method, a numerical solving strategy for multi-physical field of reacting fluid–solid flow is discussed in Section 2.2.1. The method to discrete differential equations into formatted numerical matrices is discussed in Section 2.2.2. Mesh description of computational domain is discussed in Section 2.2.3. And systematic calculation steps are summarized in Section 2.2.4.
Fig. 3. Validation of NaCl release characteristics for temperature increase.
production rate of reaction n . 2.1.3.2. Alkali metal transformation. Current assumption is that Na is primarily released in the form of NaCl gas, which mainly derives from NaCl solid [31]. The gasification of NaCl, was primarily determined by chemical reaction kinetics and Yang [34] confirmed the accuracy of this method. Table 2 lists the correlative kinetic parameters.
2.2.1. Unified equations form The main idea of fictitious domain method is that, suppose the region occupied by solid particle is also filled with fluid in flow field, as fictitious domain, and then a fictitious force is applied to the fluid in it, so that the fluid motion behavior satisfies the constraints of the rigid body motion. The specific strategy is to add different forms of the source term, named fictitious source terms, into the fluid governing equations in the region occupied by solid particle, so that the numerical results are the same as those obtained by the solid particle dynamics equations. For the momentum equation of fluid, i.e., Eq. (3), to realize the description of particle movement, the specific form of S can be derived by:
2.1.4. Species conservation equation The variation of gas species distribution in flow field can be described by the Fick's diffusion law. The specific volume fraction equation of specie m (Ym ) is:
∂ (Ym) 1 + Uf ∙∇Ym = Dm ∇2 Ym + ∇ρf Dm ∙∇Ym + ωṁ + Sm ∂t ρf
(16) 2
where Dm is the diffusion coefficient of specie Ym , the unit is: m/s , and Sm , the fictitious domain source term, adopting different forms for different conditions is discussed below. 2.1.5. Energy conservation equation The variation of temperature distribution in flow field can be described by the basic heat transfer theory. Coupling the fluid dynamics equations, the specific equation for fluid temperature, Tf , is:
∂Tf + Uf ·∇Tf ⎞ = λT , f ∇2 Tf + ∇λT , f ·∇Tf + q ̇ + ST ρf cp, f ⎛ ⎝ ∂t ⎠ ⎜
S = ρs / ρf ∙1/ Vs [Ff + |ωs
∑ Fcoll + (ρs − ρf ) Vs g] + ρs Tf /Is − ρs
× r + ∇p − ∇∙μ [∇Uf + (∇Uf )T − 2/3∇∙Uf ] − ρf g
(19)
where, ρs / ρf is density correction of solid particle and fluid, 1/ Vs [Ff + ∑ Fcoll + (ρs − ρf ) Vs g ] + ρs Tf / Is − ρs |ωs |2 × r is used to add kinetic information of solid particle, while
⎟
(a) Temperature (K) field
|2
(17)
(b) NaCl mass fraction (%) field
(c) CO mass fraction (%) field
Fig. 4. Temperature, NaCl mass fraction and CO mass fraction at t = 0.02 s. 4
Fuel 255 (2019) 115763
S. Qu and C. You
(a) Temperature (K) field
(b) NaCl mass fraction (%) field
(c) CO mass fraction (%) field
Fig. 5. Temperature, NaCl mass fraction and CO mass fraction at t = 0.15 s.
(a) Temperature (K) field
(b) NaCl mass fraction (%)
(c) CO mass fraction (%)
Fig. 6. Temperature, NaCl mass fraction and CO mass fraction at t = 0.45 s.
(a) Temperature (K) field
(b) NaCl mass fraction (%)
(c) CO mass fraction (%)
Fig. 7. Temperature, NaCl mass fraction and CO mass fraction at t = 0.49 s.
5
Fuel 255 (2019) 115763
S. Qu and C. You
Sm =
Wm C ωm, sur ρf
(20-b)
C
where, ωm, sur can be obtained by the surface reaction of the solid particle, by Eqs. (11)–(14). Here, fictitious source term, Sm , is used to add the products of heterogeneous reaction to flow field. For the energy conservation equation, i.e., Eq. (18), to realize a description of particle temperature variation, the specific form of ST can be derived by:
ST
( ∫ρ c
= −
f p, f
+ Uf ·∇Tf
Tf Uf ·nds + ρf cp, f
)ρVc
s s p, s
Fig. 8. Variation of temperature and mass residual rate of NaCl.
2.2.2. Spatial and temporal discrete solving method Due to velocity, Uf and pressure, p of flow field are coupled by the momentum equation i.e., Eq. (3), it is difficult to solve them independently. Therefore, it is necessary to design a numerical solution method. According to the Semi-Implicit Method for Pressure Linked Equations (SIMPLE), firstly proposed by SV Patankar, 1972, solving method used in this work iterates the velocity, Uf and pressure, p between the continuity equation i.e., Eq. (2) and the momentum equation i.e., Eq. (3) in discrete form repeatedly, so that the numerical values of velocity, Uf and pressure, p are more closer to the analytical solution after each iteration. Mole fraction of specie m , Ym and temperature, T of flow field are independently described by species conservation equation i.e., Eq. (16) and energy conservation equation i.e., Eq. (17), so that they can be solved respectively. Here, finite volume method (FVM) is used to discretize the differential equations, based on its advantage for satisfying integral conservation of each physical quantity in whole computational domain.
∇p − ∇∙μ [∇Uf + (∇Uf )T − 2/3∇∙Uf ] − ρf g is used to delete kinetic information of fluid to the equation. Thus, the original fluid equation can be used to calculate the solid particle movement by fictitious source term, S . For the species conservation equation of fluid, i.e., Eq. (16), to realize the description of species behavior in the region occupied by solid particle, assuming no reaction within the particles, the specific form of Sm can be derived by: 1 ∇ρ Dm ·∇Ym + ω̇m ) ρf f
(21)
where (ρf cp, f )/(ρs Vs cp, s ) is heat capacity correction of solid particle and fluid, − ∫ ρf cp, f Tf Uf ·nds + ∫ λT , s ∇Tf ·nds is used to add heat information of solid particle, while − λT , f ∇2 Tf − ∇λT , f ·∇Tf + Uf ·∇Tf is used to delete heat information of fluid to the equation. Thus, the original fluid equation can be used to calculate the particle temperature by fictitious source term, ST . Although the specific form of the governing equations is different from the original one after adding fictitious source terms, the numerical results are theoretically consistent, that is, no information of solid particles and fluids is lost. Due to it is not necessary to deduce governing equation frequently in the boundary between fluid and solid particle, effected by the strong interaction, fictitious domain strategy has obvious advantages in solving fluid-solid flow problems. As a result, it improves efficiency and stability in solving key parameters, such as velocity in whole computational domain.
Fig. 9. Particles arrangement at t = 0 s.
Sm = Uf ·∇Ym − (Dm ∇2 Ym +
∫ λT ,s ∇Tf ·nds − λT ,f ∇2Tf − ∇λT ,f ·∇Tf
2.2.3. Description of solid particle and fluid After establishing unified equations in discrete form, structural equilibrium meshes can be used to describe both solid particle and fluid. Fig. 1 shows three types of computational domains: region only occupied by fluid (only red triangles); region only occupied by solid particle (only blue rectangles); and the region partially occupied by solid particles and fluid (both red triangles and blue rectangles). The background flow field comprises Euler meshes, which are applicable to the governing fluid equations (or unified equations without fictitious source terms). The field (Euler meshes) totally covered by solid particles (solid) is applicable to the governing rigid body equations (or unified equations with fictitious source terms). During numerical calculation, Euler staggered meshes are used to store information such as velocity, Uf and pressure, p of flow field. Representing particle surface, the boundary between fluid and solid particle comprises Lagrangian points (as shown by black nodes in
(20-a)
By adding fictitious source term, Sm , D (Ym)/ Dt = ∂ (Ym)/ ∂t + Uf ∙∇Ym is 0 in Eq. (16), since there is no component distribution in solid particle under consideration. While in the region near the solid particle boundary (on the surface of the solid particle), the specific form of Sm can be derived by: 6
Fuel 255 (2019) 115763
S. Qu and C. You
(a) Temperature (K) field
(b) NaCl mass fraction (%) field
(c) CO mass fraction (%) field
Fig. 10. Temperature, NaCl mass fraction and CO mass fraction at t = 0.05 s.
(a) Temperature (K) field
(b) NaCl mass fraction (%) field
(c) CO mass fraction (%) field
Fig. 11. Temperature, NaCl mass fraction and CO mass fraction at t = 0.1 s.
physical field. The classical SIMPLE can be used to solve velocity, Uf and pressure, p of pure flow field well. On this basis, strategy adopts to fictitious source terms to modify the governing equation forms of fluid, so that the behavior of solid particle and fluid in whole computational domain can be solved in a unified way. In addition, FVM, specific meshes construction, introducing more equations, and etc. make FRS more authentic in physics. The solving process is summarized as follows:
Fig. 1) that distribute Lagrangian source terms into Eulerian meshes outside the solid particle. The background meshes, including Lagrangian points, is the region influenced by both particles and fluid. This region can also be calculated by the unified equations with fictitious source terms adding modified parameter, φ , such that:
φ = Vs, E / Vf , E
(22)
where Vs, E is the particle volume in the Eulerian mesh and Vf , E is the fluid volume in the same Eulerian mesh. The value of modified parameter, φ ranges from 0 to 1.
A. Initialize physical parameters and set structured meshes for the computational domain at the beginning of the n-th time loop; B. Calculate the force, Fs and torque, Ts on solid particle by Eqs. (6)–(10), then solve Eqs. (4) and (5) to obtain the velocity, Us , position vector, Ss and acceleration velocity, as of solid particle; C. Divide the solid particle regions and fluid regions in whole
2.2.4. Solving process For FRS of reacting fluid–solid flow, the difficulty lies in the description and solution of the strong interaction between two phases and complicated coupling relations among physical quantities in the multi7
Fuel 255 (2019) 115763
S. Qu and C. You
(a) Temperature (K) field
(b) NaCl mass fraction (%) field
(c) CO mass fraction (%) field
Fig. 12. Temperature, NaCl mass fraction and CO mass fraction at t = 0.15 s.
(a) Temperature (K) field
(b) NaCl mass fraction (%) field
(c) CO mass fraction (%) field
Fig. 13. Temperature, NaCl mass fraction and CO mass fraction at t = 0.16 s.
D. E. F.
G.
consume) mass and heat, by Eqs. (11)–(14) and correct the physical parameters of solid particle and fluid, by thermodynamic equations and the mass conservation equation, e.g., Eq. (15); H. Analyze the specific fictitious source terms, Sm and ST by Eqs. (20-a), Eq. (20-b) and Eq. (21) and substitute them into Eqs. (16) and (17); I. Discretize Eqs. (16) and (17) by FVM and solve them to obtain the temperature, Tf and species, Ym distribution of flow field; J. Judge that, if the iteration results of key parameters such as velocity, Uf and pressure, p of flow field converge, terminate the computational iteration of SIMPLE and enter step. K. Otherwise, return to
computational domain. For the meshes including the boundary between solid particle and fluid, calculate the modified parameter, φ by Eq. (22); Calculate the temperature, Ts of solid particle by Eq. (18); Analyze the specific fictitious source term, S by Eq. (19) and substitute it into Eq. (3); Start the computational iteration of SIMPLE, discretize and solve Eqs. (2) and (3) to obtain the coupled solution of velocity, Uf and pressure, p distribution of flow field; Calculate chemical reaction source terms which produce (or 8
Fuel 255 (2019) 115763
S. Qu and C. You
(a) Temperature (K) field
(b) NaCl mass fraction (%) field
(c) CO mass fraction (%) field
Fig. 14. Temperature, NaCl mass fraction and CO mass fraction at t = 0.17 s.
(a) Temperature (K) field
(b) NaCl mass fraction (%) field
(c) CO mass fraction (%) field
Fig. 15. Temperature, NaCl mass fraction and CO mass fraction at t = 0.2 s.
coupled with fictitious domain strategy, into SIMPLE, the comprehensive algorithm system is established.
step. F; K. Update the physical parameters of solid particle and fluid by the physical equations, e.g., Eq. (1) and return to Step. A to commence the (n + 1)-th time loop, until the end of a given physical time.
3. Calculation results and analyses
Step. F to J represent the whole calculating process of solving flow field by SIMPLE. Step. C, E and H embody fictitious domain strategy of unified treatment of solid particle and fluid. Step. I represents the calculating of the temperature, Tf and species, Ym equations with integral discrete scheme by FVM, based on the conservation law of mass and energy. By embedding the discrete equations solving with FVM,
Gasification behavior of alkali metals coupled with carbon combustion is the key to this problem. Therefore, carbon combustion model and Na gasification model are validated at first in Section 3.1, to prove their accuracy. Then, discrete particle combustion and multiple particles combustion in air are simulated respectively in Sections 3.2 and 3.3. 9
Fuel 255 (2019) 115763
S. Qu and C. You
boundary condition of the computational domain is open and the other three sides are walls. It can be observed from the calculation results that particle starts to burn at t = 0 from which point there is a rapid temperature increase. Figs. 4 and 5 show that heating accelerates chemical reactions from t = 0.02 to 0.15 s. Fig. 5 shows that combustion stabilizes at t = 0.15 s. The combustion rate of carbon and the production rate of gas remain unchanged until the combustion reaches its final stage. Fig. 6 shows that the combustion rate of carbon decreases after t = 0.45 s. Fig. 6(a) shows that the particle temperature rapidly decreases due to air flow and insufficient heat supply from the combustion reactions, which leads to a rapid decrease in Na release, as shown in Fig. 6(c). However, there is still a considerable amount of Na in the particle. Fig. 7 shows that the carbon combustion process completes at about t = 0.49 s. Fig. 7(a) shows that the remaining CO is gradually removed by the air flow. Fig. 7(c) shows that there is a low Na release rate in the higher temperature region. In addition, Fig. 8 shows more obviously that, in the initial stage (approximately 10% of the whole combustion), temperature rises rapidly and Na release accelerates. During the stable stage (approximately 85% of the whole combustion), temperature is constant in a range slightly higher than that of ambient medium. At this time, combustion rate of particle is unchanged as a stable heat source, while Na release slows down slightly, but is not obvious. This is mainly due to the continuous decrease of the absolute mass of Na. At the end (approximately 5% of the whole combustion), particle stops the oxidation and gasification, and no heat is generated continuously, which leads to the rapid reduction of temperature, and Na release decreases rapidly, although more than 50% of it remains in particle. Therefore, it can be concluded that Na release depends heavily on the local high temperature, which is consistent with the phenomenon obtained by experimental investigations that the chemical transformation reactions of Na are closely related to temperature [35–38]. Na release mainly occurs during the stable combustion stage due to the need for a stable heat source to provide a sustained high temperature condition.
Fig. 16. Comparison of NaCl residue variation for a single particle against the average value of 150 particles.
3.3. Multiple particles combustion Fig. 17. Frequency distribution of NaCl residue for 150 particles at 0.1 s and 0.2 s.
Corresponding to dense particles combustion, such as in the area near the boiler burner, or in the dense flow zone of fluidized bed, multiple particles combustion moving on the background meshes is simulated. Especially when the particles volume fraction is large, there are frequent collisions between particles, and between particles and walls. These collisions significantly affect the molecular transport of gas products and thus, affect the specific combustion state of particles. Figs. 10–15 shows the calculation results for a multiple particulate system where particles are moving, colliding, burning and interacting with fluid. The particle has a diameter of 1 mm, density of 1000 kg/m3 and the mass fraction of Na is 0.5%. The ambient pressure is 1.01 × 105 Pa; the temperature is 1273 K; and the velocity of air flow from the bottom of the computational domain (0.35 × 0.1 m) is 2 m/s. The boundary conditions are identical to the previous cases. At the beginning of the calculation (i.e., t = 0, as shown in Fig. 12), the group of 150 (15 × 10) particles evenly arranged at the top of the computational domain are released. Fig. 10 shows that at t = 0.05 s, the incoming flow begins to disrupt the structure of the multiple particulate system. There is a gradual heat increase on the windward side of the particles and reaction I, II, IV commence on the surface of the particles; reaction III primarily occurs at the tail region of the particles. Fig. 11 shows that at t = 0.1 s, the particles settle in the middle of the computational domain. Particles collide on the windward side and small vortexs separate from the wall in the tail region of the particles. Fig. 12 shows that at t = 0.15 s, some particles reach the bottom and begin to bounce. Here, the motion patterns of the particles and fluid in the computational domain maintain a degree of symmetry. Fig. 13 shows that at t = 0.16 s, all particles
3.1. Model validation Fig. 2 shows a validation of char combustion model. Simulation results are obtained by taking the average value of the carbon mass release rate in the stable burning stage (i.e. linear variation is satisfied). The results here are compared with those from the Makino’s experiment [32], which comprised a graphite rod of 5 mm diameter to establish a two-dimensional stagnation flow combustion under an ambient temperature of 320 K with an air flow velocity of 4.125 m/s. Errors from the current simulation are within the acceptable range. Fig. 3 shows a validation of NaCl release model, established according to Yang [34]. Simulation results are obtained by calculating the release ratio of the NaCl mass, accounting for 68% of Na content in Zhundong coal, as the temperature increases from 673 to 1473 K. The model-based results are in agreement with the experimental data. 3.2. Single particle combustion Corresponding to sparse particles combustion, such as pulverized coal boiler combustion, a single particle combustion fixed in the background meshes is simulated. Figs. 4–7 shows the calculation results. The particle has a diameter of 50 × 10−6 m, density of 1000 kg/ m3, and the mass fraction of Na is 0.5%. Ambient pressure is 1.01 × 105 Pa; temperature is 1273 K; and the velocity of air flow from the bottom of the computational domain (0.002 × 0.001 m) is 1 m/s. Top 10
Fuel 255 (2019) 115763
S. Qu and C. You
produces a large number of vortexes, which seriously affects the transport behavior of Na in gas phase. Especially for the large vortex, the entrainment makes Na in gas phase locally enriched in the flow field.
reach the bottom and the volume fraction of the particles achieves maximum compression. Frequent collisions occur in the particle group and, thus, the flow field and particles no longer have symmetry. Fig. 14 shows that at t = 0.17 s, the particles diffuse back and collide frequently, which results in chaotic movement both for the fluid and solid particles. Also, large vortexes are more likely to occur. It is necessary to highlight here that where the large vortex exists, Na in gas phase is enriched. Fig. 15 shows that at t = 0.2 s, the rebounding particles are distributed across the entire flow field. The whole simulation shows that the fluid movement has a significant impact on the molecular transport of gas products, especially the emergence of vortexes, which promote the formation of local high gas concentration due to entrainment. As shown in Fig. 16, comparing Na residue of a single particle with the average Na residues of dense particles under the same combustion conditions and with the same particle characteristics, it can be found that the average Na release of dense particles is slightly larger than that of a single particle in the initial period (< 0.05 s) of combustion. As combustion proceeds, the release rate of dense particles slows down gradually, while Na release rate of a single particle does not change significantly. The preliminary guess is that, at the beginning of combustion, the oxygen is sufficient, and particles on the windward side begin to burn, which makes the heat accumulated in particles more than that carried away by the fluid, which results in local high temperature to make Na release rate at a very high level, and even makes the average release rate of all 150 particles is higher than that of a single particle. As combustion proceeds, oxygen is rapidly consumed and even insufficient to supply, which limits the local combustion rate, prevents the local high temperature of particles from continuing, and leads to the decrease of the average release rate of dense particles. The gap of Na release rate between dense particles and a single particle is larger. As shown in Fig. 17, from t = 0.1 to 0.2 s, Na residue of particles generally decreases. Further, the centrality of the Na release rate decreases and the peak value of particle Na residue surpasses 80% at 0.1 s, then decreases to less than 75% at 0.2 s. Very few particles, guessed to be located on the windward side, lose most of their Na in a short time. It can be predicted that, as combustion proceeds, the dense particulate system becomes chaotic and the peak value of the Na residue of particles shifts leftwards. The peak distribution is more uniform, either.
In this work, the feasibility of this TURE DNS method was confirmed. Its potential can be further exploited for numerically investigating the characteristics of alkali metals release and transfer in multi-physical fields in more details. Acknowledgement This work was financially supported by the National Natural Science Foundation of China (51761125012 and 51390493). References [1] Zhou JB, Zhuang XG, Alastuey A, Querol X, Li JH. Geochemistry and mineralogy of coal in the recently explored Zhundong large coal field in the Junggar basin, Xinjiang province, China. Int J Coal Geol 2010;82:51–67. [2] Yang ZC, Liu JL, He HG. Study on properties of Zhundong coal in Xinjiang region and type-selection for boilers burning this coal sort. Thermal Power Gener Chin 2010;39:38–44. [3] Sugawara K, Enda Y, Inoue H, Sugawara T, Shirai M. Dynamic behavior of trace elements during pyrolysis of coals. Fuel 2002;81:1439–43. [4] Li GD, Li SQ, Qian H, Qiang Y. Fine particulate formation and ash deposition during pulverized coal combustion of high-sodium lignite in a down-fired furnace. Fuel 2015;143:430–7. [5] Takuwa T, Naruse I. Emission control of sodium compounds and their formation mechanisms during coal combustion. Proc Combust Inst 2007;31:2863–70. [6] Takuwa T, Mkilaha ISN, Naruse I. Mechanisms of fine particulates formation with alkali metal compounds during coal combustion. Fuel 2006;85:671–8. [7] Naruse I, Kamihashira D, Khairil, Miyauchi Y, Kato Y, Yamashita T, Tominaga H. Fundamental ash deposition characteristics in pulverized coal reaction under high temperature conditions. Fuel 2005;84:405–10. [8] Benson SA, Holm PL. Comparison of inorganics in three low-rank coals. Ind Eng Chem Prod Res Dev 1985;24:145–9. [39] Wang CA, Jin X, Wang YK, Yan Y, Cui J, Liu YH, Che DF. Release and transformation of sodium during pyrolysis of zhundong coals. Energy Fuels 2015;29:78–85. [10] He Y, Zhu J, Li B, Wang ZH, Li ZS, Alden M, et al. In-situ measurement of sodium and potassium release during oxy-fuel combustion of lignite using laser-induced breakdown spectroscopy: effects of O2 and CO2 concentration. Energy Fuels 2013;27:1123–30. [11] Wang ZH, Liu YZ, Ronald W, Wan KD, He Y, Xia J, et al. Measurement of atomic sodium release during pyrolysis and combustion of sodium-enriched Zhundong coal pellet. Combust Flame 2017;176:429–38. [12] Liu YZ, He Y, Wang ZH, Wan KD, Xia J, Liu JZ, et al. Multi-point LIBS measurement and kinetics modeling of sodium release from a burning Zhundong coal particle. Combust Flame 2018;189:77–86. [13] Van Eyk PJ, Ashman PJ, Nathan GJ. Mechanism and kinetics of sodium release from brown coal char particles during combustion. Combust Flame 2011;158:2512–23. [14] Wan KD, Vervisch L, Xia J, Domingo P, Wang ZH, Liu YZ, et al. Alkali metal emissions in an early-stage pulverized-coal flame: DNS analysis of reacting layers and chemistry tabulation. Proc Combust Inst 2019;37:2791–9. [15] Fox RW, McDonald AT, Pritchard PJ. Introduction to fluid dynamics. New York: John Wiley and Sons; 1985. p. 354. [16] Hu HH, Patankar NA, Zhu MY. Direct numerical simulations of fluid–solid systems using the arbitrary Lagrangian-Eulerian technique. J Comput Phys 2001;169:427–62. [17] Lomtev I, Kirby RM, Karniadakis GE. A discontinuous Galerkin ALE method for compressible viscous flows in moving domains. J Comput Phys 1999;155:128–59. [18] Osher S, Fedkiw R. Level set methods and dynamic implicit surfaces. Surfaces 2002;44:77. [19] Sierakowski AJ, Andrea P. Resolved-particle simulation by the Physalis method: enhancements and new capabilities. J Comput Phys 2016;309:164–84. [20] Uhlmann M. An immersed boundary method with direct forcing for the simulation of particulate flows. J Comput Phys 2005;209:448–76. [21] Glowinski R, Pan TW, Hesla TI, Joseph DD. A distributed Lagrange multiplier/fictitious domain method for particulate flows. Int J Multiph Flow 1999;25:755–94. [22] Peskin CS. The immersed boundary method. Acta Numerica 2002;11. [23] Fogelson AL, Peskin CS. A fast numerical method for solving the three-dimensional stokes' equations in the presence of suspended particles. J Comput Phys 1988;79:50–69. [24] Lai MC, Peskin CS. An immersed boundary method with formal second-order accuracy and reduced numerical viscosity. J Comput Phys 2000;160:705–19. [25] Apte SV, Martin M, Patankar NA. A numerical method for fully resolved simulation (FRS) of rigid particle–flow interactions in complex flows. J Comput Phys 2009;228:2712–38. [26] Ghasemi A, Pathak A, Raessi M. Computational simulation of the interactions
4. Conclusion A TURE DNS method, based on fictitious domain strategy with finite volume (FV) scheme, was developed to describe the behavior of alkali metals release and transfer during char combustion. Coupling fluid dynamics equations, particles dynamics equations, chemical reaction models, species conservation equations and energy equations, the introduction of SIMPLE enabled the algorithm to perform better in terms of accuracy and computational efficiency. After verifying the method against existed experimental results, a single fixed particle combustion simulation corresponding to sparse particles combustion, and multiple particles combustion, corresponding to dense particles combustion were carried out respectively. The conclusions are summarized as follows. A. Na release depends heavily on the local high temperature in flow field, and mainly occurs during the stable combustion stage (approximately 85% of the whole combustion) due to the need for a stable heat source to provide a sustained high temperature condition. B. Due to rapid heat accumulation with sufficient oxygen, the average Na release rate of dense particles is higher than that of a single particle at the beginning of combustion. With oxygen consumes, average Na release rate of dense particles decreases gradually. C. Frequent collision of particles makes the flow field disordered and 11
Fuel 255 (2019) 115763
S. Qu and C. You
[27] [28] [29] [30] [31] [32] [33]
between moving rigid bodies and incompressible two-fluid flows. Comput Fluids 2014;94:1–13. Haeri S, Shrimpton JS. A new implicit fictitious domain method for the simulation of flow in complex geometries with heat transfer. J Comput Phys 2013;237:21–45. Zhang LH, Liu K, You CF. Fictitious domain method for fully resolved reacting gassolid flow simulation. J Comput Phys 2015;299:215–28. Patankar S. Numerical heat transfer and fluid flow. CRC Press; 1980. Yu ZS, Shao XM. Direct numerical simulation of particulate flows with a fictitious domain method. Int J Multiph Flow 2010;36:127–34. Lee SHD, Teats FG, Swiift WM, Banerjee DD. SHORT COMMUNICATION alkalivapor emission from PFBC of illinois coals. Combust Sci Technol 1992;86:327–36. Makino A, Namikiri T, Kimura K. Combustion rates of graphite rods in the forward stagnation field with high-temperature airflow. Combust Flame 2003;132:743–53. Yang YM, Zhang H, Wu YX, Lv JF. Release of alkali/alkaline earth metal species in zhundong coal at different ashing temperatures. J Combust Sci Technol Chin
2015;21:297–300. [34] Yang YM. Occurrence modes and release characteristics of alkali and alkaline earth metal species in zhundong coals. Beijing: Tsinghua University; 2017. [35] Kosminski A, Ross DP, Agnew JB. Transformations of sodium during gasification of low-rank coal. Fuel Process Technol 2006;87:943–52. [36] Li WQ, Wang LY, Qiao Y, Lin JY, Wang MJ, Chang LP. Effect of atmosphere on the release behavior of alkali and alkaline earth metals during coal oxy-fuel combustion. Fuel 2015;139:164–70. [37] Song GL, Song WJ, Qi XB, Lu QG. Transformation characteristics of sodium of zhundong coal combustion/gasification in circulating fluidized bed. Energy Fuels 2016;30:3473–8. [38] Wang ZS, Wang LY, Lin JY, Wang MJ, Chang LP. The influence of the addition of sodium on the transformation of alkali and alkaline-earth metals during oxy-fuel combustion. J Energy Inst 2018;91:502–12.
12