Modeling of Combustion Process in 600 MW Utility Boiler Using

Boiler Using Comprehensive Models and Its. Experimental Validation. Jianren Fan,* Ping Sun, Xudong Zha, and Kefa Cen. Department of Energy Engineering...
0 downloads 0 Views 148KB Size
Energy & Fuels 1999, 13, 1051-1057

1051

Modeling of Combustion Process in 600 MW Utility Boiler Using Comprehensive Models and Its Experimental Validation Jianren Fan,* Ping Sun, Xudong Zha, and Kefa Cen Department of Energy Engineering, Zhejiang University, Hangzhou, 310027, P. R. China Received February 15, 1999

In this work, an Eulerian/Lagrangian approach has been employed to numerically investigate flow characteristics, heat transfer, and combustion processes in tangentially fired furnaces. Comprehensive models have been described in detail. An RNG k- model and a new cell face velocity interpolation method based on a nonstaggered grid system are employed to gain accuracy. To avoid pseudo-diffusion that is significant in modeling tangentially fired furnaces, some attempts have been made at improving the differential volume scheme. Some new developments on turbulent diffusion of particles are taken into account for improving computational accuracy. Numerical results achieve reasonable agreement with experimental data.

1. Introduction

2. Global Numerical Models

There are large amounts of anthracites and lean coals in China, which provide parts of the generated electricity. At present, tangentially fired boilers are widely used to burn anthracites and low-volatile-matter coals. The predictions of boiler characteristics are important for the boiler design and operation. It is therefore very desirable to develop an analytical model of the boiler, which can serve as a tool to optimize operation and to predict the effects of changing conditions. Generally, simulations of power plant boilers include several indispensable parts: turbulent flow and turbulent transfer processes, coal-particle motion and turbulent diffusion, turbulent combustion, evaporation devolatilization and carbon combustion, and heat transfer between gases particle and the walls of boilers. Comprehensive models of combustion processes have been the subjects of many investigations. These comprehensive models are based on the numerical solution of multidimensional differential equations for conservation of mass energy and momentum, combined with rateprocess laws or correlation, physical and chemical properties, and coefficients from experimental data. Our model predicts three-dimensional flow, combustion and trajectories of coal particles, temperatures, major species concentrations, and radiate fluxes in the boiler. In our work, we focus on improving numerical accuracy by using newly developed nonstaggered grid technology, introducing a highly accurate, highly stable differential volume scheme, and incorporating an updated k- model into the flow field simulation. The data resulting from the present study may be used to enhance the understanding of combustion processes and also provide a useful basis for further researching comprehensive models of combustion processes and designing and operating boiler furnaces with high efficiency.

In this work, a Lagrangian/Eulerian approach has been employed for gas-solid two phases flow simulation. The gas phase is described by Navier-Stokes equations, coupled with appropriate equations for density and viscosity. For closure of the turbulence equations, we use the RNG k- model. The general Eulerian equation for the gas phase takes the following form:

* Author to whom correspondence should be addressed. Fax: +86571-7991863. E-mail: [email protected].

(

)

∂(Fu bjΦf) ∂Φf ∂ ) Γ + SΦ ∂xj ∂xj Φ ∂xj

(1)

where f is a flux coefficient, which is used to describe the characteristics of the platen. To make a precise simulation for platen and convective superheater, we use the method of multi-porosity combined with distributed resistance, that is to consider the platen as a sort of porous media, inside which the pressure gradient along the flow direction is redistributed according to thermal calculation. The correcting surface flux coefficient on its each surface as well as the volume flux coefficient are defined by its structure. In the upper equation, Φ represents the variables u, v, w, k, and , (coordinate velocities, turbulent fluctuation energy, and turbulent energy dissipation) as done in ref 1.1 A modified SIMPLER method2 has been employed to determine velocities and pressures using a nonstaggered grid system in Cartesian coordinates as demonstrated in the following part. To eliminate the oscillations with the present nonstaggered grid arrangement, a new treatment of the locally linear convection term2 was introduced. (1) Fan, J. R.; Liang, X. H.; Xu, Q. S.; Zhang, X. Y.; Cen, K. F. Energy 1997, 22, 847-857. (2) Rahman, M. M.; Miettinen, A.; Siikonen, T. Numer. Heat Transfer, Part B 1996, 30, 291-314.

10.1021/ef9900233 CCC: $18.00 © 1999 American Chemical Society Published on Web 08/20/1999

1052 Energy & Fuels, Vol. 13, No. 5, 1999

u*w )

u*W + u*P

{

2

λd Kw

|

∂fu ∂x

w

Fan et al.

-

[21(∂p* ∂x |

- Cw

W

+

| ]}

|)

∂p* ∂p* ∂x P ∂x w

(2)

where λd is a scaling factor with typical values of 0.5 to 1.0, Kw ) ∆xw/4(VP/AP - VW/AW), f u is the source term from which the pressure gradient term is excluded, and Cw ) Vw/Aw. Variables with subscript w represent those at cell face, while those with majuscule subscript W represent variables at nodes westward of controlled volume P. 2.1 Mathematical Models for Gas and Solid Flow. 2.1.1 Models for Gas Flow. How to reduce pseudodiffusion in numerical simulation has been a long-time concern in furnace simulations, especially in tangentially fired boiler furnaces where the angle between inlet jet direction and x or y coordinate is nearly 45° which is supposed to lead to the highest degree of pseudodiffusion. One way to reduce pseudo-diffusion is to adopt an adjustable grid system to keep inlet jet direction as parallel to the x or y coordinate as possible.3 However, that requires not only a lot of re-programming but also some reconsideration on solving method and wall function, etc. In this work, we focus on an improved differential volume method. To solve eq 1 numerically, it should be written in the numerical form

APΦp ) AEΦE + AWΦW + ANΦN + ASΦS + ATΦT + ABΦB + S (3) where A is a coefficient, and subscripts E, W, N, S, T, B represent the six neighboring points of control volume P. To apply an artificial viscosity method, the x-direction momentum equation is written into

(

∂(Φ) ∂(Φ) ∂(Φ) ∂2Φ ∂2Φ 1∂p + 2 + u +v +w )+ν ∂x ∂y ∂z F ∂x ∂x2 ∂y ∂2Φ ∂ ∂Φ ∂Φ 2 + νj cos A + νj cos A cos B + ∂x ∂x ∂y ∂z2 ∂ ∂Φ ∂Φ + νj cos A cos B νj cos A cos C + ∂z ∂y ∂x

) (

(

νj cos2 B

) (

)

∂Φ ∂Φ + + νj cos B cos C ∂y ∂z

∂ ∂Φ ∂Φ + νj cos C cos B + νj cos A cos C ∂z ∂x ∂y

)

∂Φ νj cos C + S (4) ∂z 2

in which, A, B, and C are the angles between Cartesian coordinate and flowing direction. A standard k- model performs well in predicting flows without swirling flow or without sharp change within calculated region. But for tangentially fired boiler furnaces, swirling flow is very marked and we must resort to other more valid, more efficient turbulent models to gain accuracy. In this paper, we try to use RNG k- model4 as an alternative substitution for a standard k- model. Generally, k and  equations are follows:

uj uj

∂ ∂k ) PK -  + ∂xj ∂xj

[( ) ] [( ) ] )

∂  2 ∂ ) C1 PK - C2 + ∂xj k k ∂xj

PK ) 2νTSijSij, Sij )

(

ν+

νT ∂k σT ∂xj

ν+

νT ∂ σT ∂xj

(6) (7)

k2 1 ∂ui ∂uj + , νT ) Cµ , 2 ∂xj ∂xi  k ) uiui/2,  ) ν

∂ui ∂ui (8) ∂xj ∂xj

For an RNG k- model Cµ ) 0.085, C1 ) 1.42 - η(-η/ η0)/(1 + βη3), C2 ) 1.68, σT ) 0.7179, C ) 0.7179, η ) h ij)1/2, β ) 0.015, and η0 ) 4.38. Sk/, S ) (2S h ijS 2.1.2 Models for Particle Flow. The equation of motion for a particle is

bp/dt) ) mpb g+ mp(dv 1 b -b v p)|v bg - b v p|Ap + F Breaction (9) FC (v 2 D g where F Breaction ) V B J(-dmνp/dt) is the force exerted when devolatilization occurs, the direction of the reaction force is modeled stochastically. When the particle has advanced in the furnace, reaction processes such as vaporization, devolatilization, and char combustion should be taken into consideration. In an attempt to model turbulent dispersion of particles, we take the stochastic model as recommended by Fan et al.5 2.2 Models for Heat Transfer. The heat transfer equation is

(FucpT)∆y∆z + (FvcpT)∆x∆z +(FwcpT)∆x∆y ) λxT∆y∆z + λyT∆x∆z + λzT∆x∆y + Qradiation 4KσT4∆V + Qreaction∆V (10)

where νj is artificial viscosity,

|V|∆ νj ) [cth(Re∆/2) - 2/Re∆] 2 Re∆ )

|V|∆ , ∆ ) (∆x + ∆y + ∆z)/3, ν |V| ) x(u2 + v2 + w2) (5)

(3) Patankar, S. V.; Spalding, D. B. Int. J. Heat Mass Transfer 1972, 15, 1787-1806.

where Qradiation is the radiation energy transfer, which is modeled by using a Monte Carlo method. Qreaction is the heat energy released during coal combustion, which is modeled as a source term after completing particle tracing. (4) Yakhot, V.; Orszag, S. A. J. Scientific Computing 1986, 1, 3951. (5) Fan, J. R.; Zhang, X. Y.; Chen, L. H.; Cen, K. F. Chem. Eng. J. 1997, 66, 207-215.

Modeling of Combustion Process in 600 MW Utility Boiler

Energy & Fuels, Vol. 13, No. 5, 1999 1053

Table 1. Analysis of Coal Used proximate analysis (%)

ultimate analysis (ar%)

MS

Ash

VM

FC

C

H

O

N

S

heating value (MJ/kg)

8.66

24.24

25.36

41.74

53.05

3.86

8.3

1.00

0.89

20.93

2.3 Models for Combustion. Diffusion-limited vaporization of moisture from the coal particle is described by

rw ) MwNumCgDwmAp(Xwp - Xwg)/dp(1 - Xwprp/rw) (11)

The conservation of particle energy is

mp(dhp/dt) ) Qc + Qr + Qw + Qv + Qh

(23)

dV/dt ) (V′daf - V)‚K‚exp(-E/RT)

(12)

where Qc, Qr, Qw, Qv, and Qh are referring to heat related to gas-particle convection, radiation, and heat release during evaporation, devolatilization, and charoxidation along the trajectory of the particle. The EBU-Arrhenius model models the gas-phase turbulent combustion process

V′daf ) QVdaf

(13)

Wfu ) min(Wfu,A, Wfu,T)

Fu et al.6 use the following equations to describe the formation of volatile matter from coal dust:

Char is produced in competition with volatile production as expressed by

rhm ) rv(1 - Ym)/Ym

(24)

(14)

 Wfu,T ) F min(Yfu, YO2)CE R

(15)

Equations with the form of dx/dt ) f(x,t) are integrated using a five-stage Runge-Kutta method as described in reference 8.1

The combustion of char is calculated as follows:

rhl ) rphrch/(rph + rch)

Wfu,A ) AFYfuYO2 exp(-E/RT)

3. Results and Discussions

where

rch ) A exp(-E/RT)

(16)

in which, A ) 0.08596 kg Pa-1 s-1 m-2, E ) 1.5 × 105 kj mol-1 kg-1

rph ) 4DMΦ/R′dp(Tp + Tg)

(17)

( )( )

D ) D0

P0 Tg P T0

1.75

(18)

in which D0 ) 3.13 × 10-4 m2/s. The ratio of CO/CO2 takes the following expression recommended by Mitchell:7

(

)

Eco (moles of CO) ) Aco exp R*Tp (moles of CO2)

(19)

where, Aco ) 2.630 × 1011, Eco ) 73200.0 cal/mol. The total reaction rate for the coal particle is

rp ) rv + rhl + rw

(20)

The carbon reaction rate is

rh ) rhm - rhl

(21)

The mass change of a particle follows

dmp/dt ) mw + mv + mh

(22)

where mw, mv, and mh are the rate of mass change when the coal particle is evaporated, devolatilized, and charoxidized, respectively. (6) Fu, W. B.; Zhang, Y. P.; Han, H. G.; Duan, Y. N. Combust. Flame 1987, 70, 253. (7) Mitchell, R. E. 22nd Symp. (Int) Combust., [Proc.] 1988, 69-78.

The boiler we studied is made by CE corporation. It serves 600 MW(e) output. The furnace dimensions are 19.558 m × 16.4325 m with 15 levels of burners. The analysis of the coal used is presented in Table 1. Using 16 different size classes represents a continuous distribution of coal-particle sizes. The mass fractions for various size classes are shown in Table 2. The numerical procedure for gas-phase flow is based on a well-known volume finite discretization of nonstaggered grids, which has been developed by Patankar and Spalding.3 A five-stage iteration method is used in integrating the conventional difference equations for the coal particles.1 The simulations are carried out considering half of the furnace with symmetry conditions at the middle plane using a grid comprising 35 × 37 × 143 control volumes. We start to make the calculation by solving the gas flow-field equations assuming that the particles are absent. Using the flow-field particles’ trajectories, we can determine particle temperatures and burnout histories. In the case of gas-solid flow, especially when the particle loading is high, the back-influence can be so large that the solution diverges. This can be explained by the coupling of the Lagrangian with the Eulerian part of the calculation, which is done by means of the particle source terms. The particles travel along a finite number of trajectories, leaving behind a field of source terms for the fluid equations. This source term field must be continuous in a flow domain in order to obtain a converged solution of the fluid equations. The mass, momentum, and energy source terms for each cell are calculated. The source terms are included in the gasphase equations and the flow field is then recalculated. The process is repeated until further repetition fails to change the solution. Thus, the mutual interaction of the gas and particles is accounted for.

1054 Energy & Fuels, Vol. 13, No. 5, 1999

Fan et al.

Table 2. Percentages of Various Size Classes of Pulverized Coal particle diameter (µm) percentage particle diameter (µm) percentage

5 20.22 160 1.42

20 28.04 180 0.94

40 17.92 200 0.62

60 11.61 220 0.41

80 7.58 240 0.27

100 4.97 260 0.18

120 3.27 280 0.12

140 2.16 300 0.25

Figure 2. Trajectories of particles.

Figure 1. Flow field.

Figure 1 shows the predicted results of flow field of a cross-section. In this work, we apply the pressure gradient which is calculated by thermal calculation into SIMPLER method to coupled calculate the flow field near the platen. Flow field shown in Figure 1b is taken from a cross-section of one of the secondary burners. As

seen from the figure, our method of porosity combined with distribution resistance is reasonable. Figure 2 shows the predicted tracks of four different particles with different sizes of 20, 50, 80, and 110 µm diameter. The resident time of a 20 µm diameter particle is 20 s, whereas the resident time of a 50 µm diameter particle is 14 s. From the figure, we can conclude that the smaller particles follow the gas flow better. Figure 3 shows that the RNG k- model is closer to experimental data than the standard k- model. As shown in Figure 3b, the gradient of the k curve obtained from the RNG model is bigger than that obtained from the standard k- model. This is because the coefficients of the RNG model are closer to the practice of the strongly swirling flow in the furnace. It has a better

Modeling of Combustion Process in 600 MW Utility Boiler

Energy & Fuels, Vol. 13, No. 5, 1999 1055

Figure 3. Comparison between standard k- model and RNG k- model. (A) Velocity profile on cross-section (----- RNG k- model, s standard k- model, and 2 experiment). (B) Turbulent kinetic energy profile on cross-section (s RNG k- model, s standard k- model).

Figure 5. Comparison between numerical and measured temperature.

Figure 4. Predicted temperature distribution for base case (T(K)) (a) vertical central cross section, (b) secondary burner plane.

effect than the standard k- model. That the RNG k- model is superior to the standard k- model in nonlinear flow has been validated in ref 4. In the RNG k- model, the introduction of scale η ) Sk/ enables the model to be efficient and effective for nonlinear flows featuring sharp change, swirling. The RNG k- model can be applied to both high Reynolds number flows and low Reynolds number flows. So there is no need to introduce wall-function near wall when the RNG k- model is applied. Figure 4 shows the distribution of the temperature field. It can be deduced by eq 10 that the temperature field is determined by Qradiation, Qconvection, and Qreaction, in which Qradiation is the radiation energy transfer, Qconvection is the convective energy transfer, and Qreaction is the heat energy released during coal combustion. Figure 5 shows the comparison between the predicted temperature and the experimental data on the central line of the cross-section of the furnace chamber. As we can see, the prediction yields reasonable results. Figure 6 shows the comparison between our calculated value and the measured value of the convective superheater. The superheater consists of 25 rows of tubes located near the furnace exit. It is obvious that the deviation along the direction of depth is caused by the deviation of the velocity field near outlet of the furnace. As seen from the figure, the results of calcula-

1056 Energy & Fuels, Vol. 13, No. 5, 1999

Fan et al.

Figure 6. Comparison between numerical and measured temperature of the connective superheater (°C).

Figure 8. Contour of CO2 volume fraction.

Figure 7. Contour of oxygen volume fraction (A) on central cross section, (B) on the cross-section of a burner plane.

tion and measurement show the same trend, while the introduction of the Monte Carlo method results in nonuniform temperature. The nonuniformity of tem-

perature resulted from the characteristic of Monte Carlo methods. Both sizes of grids and the numbers of energy beam are influential factors, especially the later. Further division of the energy beam can achieve accurate prediction, but it will make the prediction time-consuming. Figures 5 and 6 show that predicted temperature deviates from the experimental results somewhere. The deviation sometime can be as high as 15%-25%. That is partly due to the sort of method we use to predict the temperature field. The causes are not only the convective term but also the distribution of reactive heat. In this work, the reactive heat distribution is taken by finishing pulverized coal combustion and then summing the source term of heat. Obviously, detailed grouping of particles and multi-iteration are necessary to reach an accurate prediction. However, it will make the calculation process very time-consuming. Besides, there are also simplifications in the experiment; for example, we assume the temperature of the water-steam mixture in all tubes along the furnace width is same. In fact, a deviation in the water-side does exist. Figures 7, 8, and 9 show the contours of oxygen volume fraction and CO2, NOx concentration in the furnace. It is observed that the oxygen concentration is relatively high near the burners, whereas it is low at the outlet of the furnace. It is because oxygen is highly supplied through the burners and consumed quickly due to the process of combustion. Figure 8 shows the CO2 concentration profile, as we can see that it is relatively high in the middle of the furnace due to rapid reaction and became rather higher due to further oxidization of CO. As the reasons of the high oxygen concentration near the burners, the rate of reaction is high. Figure 9 shows the profile of the distribution of NOx that is closely dependent on the profile of O2. The rapid

Modeling of Combustion Process in 600 MW Utility Boiler

Energy & Fuels, Vol. 13, No. 5, 1999 1057

Foundation of China and Zhejiang Provincial Natural Science Foundation of China. Nomenclature

Figure 9. Contour of NOx volume fraction.

devolatilation process leads to the high volatile-NOx concentration near the burners. The high temperature and high oxygen density in the middle of the furnace lead to the high thermal-NOx concentration. 4. Conclusions We have made several attempts at improving accuracy in simulating the combustion process inside tangentially fired furnaces. To gain a reasonably good prediction of aerodynamics field, we made comparisons between the standard k- model and the RNG k- model. We found that the RNG k- model can make a relatively accurate prediction. We also introduced a stochastic turbulent diffusion model into particle tracing. As a result, our numerical results achieve reasonable agreement with experimental data. Acknowledgment. This present work has been financially supported by the National Natural Science

A ) coefficient in numerical equation Ap ) coal-particle area c ) specific heat CE ) coefficient in EBU model Cg ) gas molar concentration dp ) char particle diameter d ) particle diameter D0 ) binary diffusivity Dw ) binary diffusivity E ) activation energy f ) source term excluding pressure gradient F ) force h ) enthalpy g ) gravitational acceleration k) turbulent kinetic energy m ) mass M ) molecular weight Nu ) Nusselt number p ) pressure Q ) heat energy r ) particle reaction rate R ) thermodynamics constant rhl ) coefficient of reaction rate Re ) Reynolds number S ) source term t ) time T ) temperature u ) velocity along x-direction v ) velocity along y-direction ∆V ) cell volume V ) averaged velocity w ) velocity along z-direction W ) reactive rate x ) coordinate X ) mole fraction y ) coordinate Y ) concentration z ) coordinate Greek letters Φ ) general variables  ) turbulent dissipation µ ) viscosity νj ) artificial viscosity F ) density ∆ ) control volume width λ ) conducting coefficient Subscripts r ) radiation c ) convection daf ) Daf coal w ) evaporation v ) devolatilization h ) char-oxidation EF9900233