- Email: [email protected]
On the application of different turbulence models for the computation of fluid flow and combustion processes in small scale wood heaters
Experimental Thermal and Fluid Science 21 (2000) 99±108
www.elsevier.nl/locate/etfs
On the application of dierent turbulence models for the computation of ¯uid ¯ow and combustion processes in small scale wood heaters H. Knaus, S. Richter, S. Unterberger *, U. Schnell, H. Maier, K.R.G. Hein Institute of Process Engineering and Power Plant Technology (IVD), University of Stuttgart, Pfaenwaldring 23, D-70550 Stuttgart, Germany Received 20 June 1999; received in revised form 15 November 1999; accepted 1 December 1999
Abstract In order to optimise the combustion process in small scale wood heaters with a thermal capacity of up to 15 kW using natural wood logs an important task is to investigate the mixing and combustion processes within the reaction zones of such ®ring systems. For a basic understanding of the formation and decomposition of dierent gas components and the dependency on mixing intensity and therefore on turbulent behaviour of the reacting ¯ow numerical studies were carried out using the simulation program AIOLOS. A comparison of results obtained during the numerical modelling studies applying the k,e-, low Reynolds-number k,e- and the Reynolds-stress turbulence model on dierent grids with detailed experimental data shows a good correspondence both for isothermal conditions as well as during combustion. The main aim is to validate the numerical model enabling its use as an engineering tool to investigate the in¯uence of parameters such as mixing conditions, combustion air distribution and furnace geometry on the combustion process within the wood heater on a complete combustion for reduced emission of unburnt components, such as CO, hydrocarbons and soot. Ó 2000 Elsevier Science Inc. All rights reserved. Keywords: Wood combustion; Modelling studies; Turbulence models; Homogeneous gas phase reactions
1. Introduction One way of reducing the global CO2 emission problem is to substitute wood for fossil fuels in heat and power generation to a certain extent, an option also appropriate for small scale combustion systems for domestic heating purposes. However, the combustion of natural wood logs in small scale stoves with a thermal capacity in the range of up to 15 kW for domestic heating purposes often involves higher emissions of CO, unburnt hydrocarbons and particles. It is estimated that in Germany, for example, the share of small scale wood combustion systems contributing to the emissions of CO, hydrocarbons and soot is between 16% and 40%, although the amount of the total energy consumption is only about 1%. In order to improve the emission behaviour of small scale wood combustion systems a basic understanding of the combustion process is essential. Facilitating the in*
Corresponding author. Tel.: +49-711-685-3572; fax: +49-711-6853491. E-mail address: [email protected] (S. Unterberger).
sight into the complex inter-linked phenomena of chemical reactions and turbulent ¯ow ®eld behaviours as well as for the investigation of parametric eects, the simulation program AIOLOS is used as an eective tool for analysis and further optimisation of the combustion process. The numerical modelling studies were based on boundary conditions determined by detailed pro®le measurements of gas concentrations, temperatures, velocities and turbulence intensities within the burnout zone of the selected test stove [1]. 2. Design of the test stove A tile stove heating insert with a thermal capacity of 10 kW representing the state-of-the-art was selected as a test stove. The design of the stove and the area investigated by the mentioned pro®le measurements is shown in Fig. 1. The commercially available stove consists of two spatially separated reaction zones. In the ®rst reaction zone, both the wood fuel is gasi®ed and the gases released by the wood logs start to combust. In the secondary reaction zone, preheated air is added to the gases for almost complete burnout. The injection of the
0894-1777/00/$ - see front matter Ó 2000 Elsevier Science Inc. All rights reserved. PII: S 0 8 9 4 - 1 7 7 7 ( 9 9 ) 0 0 0 5 9 - X
100
H. Knaus et al. / Experimental Thermal and Fluid Science 21 (2000) 99±108
Fig. 1. Investigated tile stove heating insert, thermal capacity 10 kW.
secondary air is realised in two ways: by means of three nozzles before the throat and by a slot in the middle of the throat.
3. Experimental investigations 3.1. Experimental set-up The experimental test ®eld consists of measurement devices for gas concentrations, temperatures and velocities. The gas analysis system is designed to measure the gas concentrations of the main components CO2 , CO and O2 in the ¯ue gas as well as in the burnout zone at the same time. Additionally H2 O and hydrocarbons can be measured either in the ¯ue gas or in the burnout zone. The suction probes for the pro®le measurements are arranged at dierent points in the secondary reaction zone. The arrangement of the measuring points is shown in Fig. 2. Pro®le measurements were carried out by measuring successively at the positions 1±11 for a duration of approximately 1 min per single probe. Gas temperature measurements were done by using Ni±CrNi
thermocouples which were installed inside the gas sampling suction probes. The measured gas temperatures were corrected according to a formula obtained by using a suction pyrometer. Wall temperatures were measured by thermocouples at two measuring points in the burnout combustion chamber as well as the temperature of the preheated secondary air before entering the nozzles and the slot. The Laser-Doppler Velocimeter (LDV) system used for the ¯ow ®eld characterisation consisted of a 2 W argon-ion laser used with a DANTEC/invent 2D FiberFlow system facilitating the simultaneous measurement of two perpendicular velocity components in x- and z-directions, respectively. 3.2. Determination of boundary conditions and pro®le measurements In a ®rst step the boundary conditions at the dierent inlets of the burnout zone were determined at isothermal and combustion conditions, respectively. Velocity measurements were carried out on several pro®les covering the total depth (270 mm) of the channel between the gasi®cation/combustion zone and the burnout zone as
Fig. 2. Measuring positions and locations for comparison of predictions and experimental data.
H. Knaus et al. / Experimental Thermal and Fluid Science 21 (2000) 99±108
well as in the outlet area of the slot and the three nozzles for obtaining the combustion air volume ¯ows added to the reacting main ¯ow. In contrast to the isothermal case, the determination of velocity boundary conditions during combustion was done by a reduced number of measuring points due to the non-stationary combustion process. In order to get reliable data during combustion, results of several dierent burn cycles were used to calculate mean values of velocities and turbulence intensities. In addition to the above-mentioned investigations by LDV, gas concentration measurements of CO2 , CO, O2 and hydrocarbons were done at the inlet of the burnout zone. Temperature measurements were carried out for the entering combustible gases, the added burnout air and at two wall sections in order to get complete data sets characterising the combustion conditions in the investigated test stove during the dierent burn phases mentioned above. Furthermore, velocity pro®le measurements were carried out at isothermal conditions and for combustion tests, respectively. Gas concentration and temperature pro®le measurements at dierent positions in the burnout zone according to Fig. 2 were carried out for the determination of combustion conditions during a complete burn cycle. 4. Numerical modelling studies The present calculations were performed with the ®nite-volume code AIOLOS [2]. This simulation program has been developed for the numerical calculation of stationary, turbulent reacting ¯ows. The SIMPLECmethod [3] is applied to compute the pressure. The pressure-weighted interpolation method (PWIM) [4] is used to prevent the decoupling of velocities and pressure on the non-staggered grid. The convection ¯uxes are approximated by the upwind or higher-order dierencing schemes [5]. In the following the physical models used for the simulation of the wood heater are shortly described. 4.1. Kinetic model of chemical reactions The combustible volatiles are converted by the following reactions {C1}±{C3}: hm ni a O2 fC1g : Cm Hn 2 4 hn ni n ! mCO 1 a H2 a H2 O 2 4 2 1 fC2g : CO O2 ! CO2 2
2
1 3 fC3g : H2 O2 ! H2 O: 2 The distribution of water and hydrogen as products of the hydrocarbon combustion described by reaction {C1} is calculated from the water-gas-shift equilibrium.
101
All reactions are treated as irreversible reactions and the kinetic rates of the reactions {C1}±{C3} are taken from [6,7]. 4.2. Eddy dissipation concept (EDC) The EDC is a general concept for treating the interaction between turbulence and chemistry in ¯ames [8]. The method is based on a detailed description of the dissipation of turbulent eddies. In the EDC the total space is subdivided into a reaction space, called the `®ne structures' and the surrounding ¯uid. In the presented reaction scheme the reactions {C1}±{C3} are treated as taking place only in these ®ne structures which means that the reactions only take place in the smallest turbulent length scales. 4.3. Thermal radiation The discrete ordinates method in an S4-approximation [9] is used to solve the radiation transport equation. Since the intensity of radiation depends on absorption, emission and scattering characteristics of the medium passed through, a detailed representation of the radiative properties of a gas mixture would be very complex and currently beyond the scope of a 3D-code for the simulation of industrial combustion systems. Thus, contributing to the numerical eciency, simpli®cations are introduced, even at the loss of some accuracy. The absorption coecient of the gas phase is assumed to have a constant value of 0.2/m. 4.4. Transport equations and turbulence modelling In the following the conservation equations of mass, momentum, scalar quantities and turbulence quantities and their discretisation are described. All formulas and modelling assumptions presuppose high Reynoldsnumbers and steady-state ¯ow conditions. It is assumed that the ¯ow ®eld is weakly compressible which means that the density depends only on temperature and ¯uid composition but not on pressure. 4.5. Finite volume discretisation A body-®tted grid is used to represent the slope of the furnace and especially of the burner quarls exactly. Therefore the conservation equations have to be formulated in general curvilinear co-ordinates. Regarding general curvilinear co-ordinates, a huge number of different formulas exist, in particular for the momentum equation [10]. In this paper the following equations proposed by Peric [11] are used, which provide a strong conservative form of the conservation equations: Continuity: o Uj q 0: onj
4
102
H. Knaus et al. / Experimental Thermal and Fluid Science 21 (2000) 99±108
Momentum: o l oui j ouk m j Uj quj ÿ Bm bi bk pbji onm onj J onm bjk qu0i u0k J qgi : Scalar: o C/ o/ j Uj q/ ÿ B Js/ : J onm m onj
o leff oe j e Uj qe ÿ B Ce1 Pkk ÿ J qCe2 e; re J onm m onj k 5
6
The contravariant velocity components Uj and the tensor Bij are de®ned by Uj J bjm um ;
Bij bjk bik ;
7
where J denotes the Jacobian determinant and the elements of the inverse Jacobian matrix bij represent the coordinate transformation of Cartesian into general curvilinear co-ordinates bij
oxi : oyj
8
The components uj of the velocity vectors are expressed in terms of the spatially invariant Cartesian basis. 4.6. k,e-turbulence model For the k,e-model a ®rst-order turbulence model closure according to Boussinesq is used. In analogy to laminar ¯ows the Reynolds-stresses are assumed to be proportional to the gradients of the mean velocities l oui n oum 1 2 ouk k 2 bm bi ÿ dmi bk ÿ dmi qk 9 smi t J onn on1 onk 3 3 introducing lt as a proportional factor of the turbulent viscosity. The model is based on the Prandtl±Kolmogorov approach mt Cm k 1=2 L
10
and the following equation for e, resulting from a dimensional analysis: k 3=2 11 e Ce L where L and k 1=2 are the characteristic length scales and velocities of turbulence. Eqs. (10) and (11) result in an equation for the turbulent viscosity dependent on the turbulent kinetic energy k and the turbulent dissipation e k2 mt q; leff lt l1 ; 12 e with the dynamic viscosity ll . For the turbulent kinetic energy and the turbulent dissipation, transport equations can be developed from the Navier±Stokes equations assuming an isotropic turbulence: o l ok j Uj qk ÿ eff Bm Pkk ÿ J qe; 13 rk J onm onj lt Cl q
P ÿlt
oui n ouj m b b onn j onm i
oui n b onn j
ÿ
14
lt oq i b q2 oni j
op i b; oni j
15
where P is the term for the production of turbulence by gradients of the mean velocities. The parameters used in the turbulence model equations are chosen according to Launder and Spalding [12]. 4.7. Low Reynolds-number k,e-turbulence model The standard k,e- and the Reynolds-stress model are restricted to high Reynolds-number applications. In the test stove investigated, however, Reynolds-numbers might become low, especially in the upper part of the burnout section. Moreover, the in¯uence of the walls on free ¯ow characteristics may not be adequately modelled by standard wall functions as used for the k,e-model and the Reynolds-stress model in such narrow geometries. The k,e-model has thus been extended to enable its use at low Reynolds-numbers and to describe more accurately the ¯ow close to a solid wall. Based on the ®ndings of Patel et al. [13] the approach of Lam and Bremhorst [14] was selected from the various low Reynolds-number models available. The Lam and Bremhorst model combines good performance with simple numerical formulation and it diers from the standard k,e-model by the inclusion of the following functions f to modify the constants C. Turbulent viscosity equation (12): ClLRM Cl fl :
16
Turbulent dissipation rate equation (14): LRM Ce1
LRM Ce2 Ce2 f2 :
Ce1 f1 ;
Reynolds-number functions: 2 20:5 1 ; fl 1 ÿ exp ÿ 0:0165 Ry RT f1 1 k2 q ; RT le
0:05 fl
3 ;
f2 1 ÿ exp ÿR2T ;
p kyq ; Ry l
17
18
19
20
where y is the distance normal to the wall. At high Reynolds-numbers RT or Ry , the functions ``f'' tend to unity and the turbulent quantities are thus determined analogously to the standard k,e-model. For the low Reynolds-number model the transport equations of k and e are integrated up to the wall attempting to model relaminarisation eects by the functions f described above. In Fig. 3, the functions fl , f1 and f2 are
H. Knaus et al. / Experimental Thermal and Fluid Science 21 (2000) 99±108
oui m ouj m 0 0 0 0 Pij ÿ quj uk b qui uk b onm k onm k lt oq k op k oq k op k b b ÿ b b : ÿ qJ onk j onk i onk j onk j
Fig. 3. Low Reynolds-number functions f1 , f2 and fl over the wall distance y .
depicted over the wall distance y de®ned by the following equation: y
Cl qk 1=2 y : l1
21
It can be seen that the functions f1 and f2 only have a strong in¯uence for wall distances y less than 5, whereas the function fl has to be considered in the whole boundary layer. Function f2 has a steep gradient in a wall distance y between 0 and 5 and for this reason the grid must be extremely ®ne especially in this part of the boundary layer to get a sucient resolution. 4.8. Dierential Reynolds-stress model (RSM) Flows with non-isotropic turbulence cannot adequately be described by a two-equation model like the k,e-model. Under such ¯ow conditions turbulence should be described by dierential transport equations for all Reynolds-stress components o 2 Uk qu0i u0j Dij Pij Wij ÿ J dij qe; onk 3
24
The expression in the ®rst brackets describes the interaction between the stresses and the rate of the mean strain. These terms are exact and need not be modelled. The expressions in the second brackets arise only under non-uniform density conditions and originate from Favre-averaging. Interactions between pressure and strain leading to a redistribution process within the Reynoldsstresses towards an isotropic state [16] are represented by the term Wij and are modelled as follows: e 2 1 qu0i u0j ÿ dij qk ÿ C2 Pij ÿ dij Pkk ; Wij ÿ JC1 k 3 3 25 where Wij consists of two distinct parts which describe the relations between purely turbulent quantities and between turbulence and mean strain. The dissipation rate e of the Reynolds-stresses is modelled using transport equation (14) as described for the standard k,e-model. It implies the isotropy of the dissipation process at smallest length scales. A summary of the constants appearing in the above equations is given in [5].
5. Comparison of experimental results and simulation studies Two body-®tted grids with dierent ®neness, one with 60 000 and the other with 200 000 nodes, were used to model the complex geometry of the secondary reaction zone. In Fig. 4 the coarse grid is shown. The ®ne
22
where Dij , Pij and Wij represent the diusion, production and pressure±strain processes, respectively. The diusion of the Reynolds-stresses Dij is approximated by ! o leff ou0i u0j k B : Dij 23 onk rk J onm m This expression is strongly simpli®ed compared to the generalised gradient-diusion hypothesis of Daly and Harlow [15]. The replacement of the detailed non-isotropic viscosity by the isotropic viscosity drastically reduces the complexity of the diusion term and promotes the numerical stability of the model. The generation term of the Reynolds-stress transport equations is given by
103
Fig. 4. Coarse computational grid.
104
H. Knaus et al. / Experimental Thermal and Fluid Science 21 (2000) 99±108
grid oers a sucient resolution of the boundary layer needed for the low Reynolds-number k,e-model. For most of the cell centres adjacent to the wall y is less than 5 and the largest distance is 40 for isothermal as well as for combustion conditions. It was not possible, however, to put a constant number of cells in the laminar boundary sublayer because of the complex geometry and the complex ¯ow situation within the reaction zone. 5.1. Isothermal conditions Because of the non-stationary processes of wood log combustion and the resulting measuring problems the ®rst simulations were done without combustion. Results of the k,e-model and the RSM are compared on the two grids to judge the grid dependency of the solutions. Comparing the measured and computed ¯ow ®eld a qualitatively good agreement is observed. For the simulation using the coarse grid and the k,e-model the location of the centre of the recirculation zone is well predicted (Fig. 5), but it is striking that this zone in the upper part of the furnace is more noticeably shaped in the measurements. On the ®ne grid using the k,e-model the recirculation zone is predicted as too strong and its centre is located more in the middle of the furnace. The best results are given by the RSM. The location of the recirculation zone as well as its strength is qualitatively in good accordance with the measurements. To validate the model quantitatively additional comparisons of velocity measurements and simulation results were done along dierent lines which are depicted in Fig. 2. Looking at the height of z ÿ20 mm the prediction of the w-velocity component is good for both grids using the k,e-model and the RSM (Fig. 6). However, on the ®ne grid the peak at x 10 mm is better predicted by the RSM than by the k,e-model. At a height of z 50 mm the k,e-model can predict neither on the coarse nor on
the ®ne grid as good as the RSM. Especially on the ®ne grid the shape and the location of the w-velocity maximum agrees well with the measurements for the RSM. Because of these results it was decided to continue with the simulations only on the ®ne grid and additionally to implement the low Reynolds-number k,emodel. This model turned out to give the best predictions. At z ÿ20 mm the simulation results of all turbulence models are very similar (Fig. 7). However the predictions of the magnitude and the gradient of the w-velocity at z 50, 100 and 180 mm are superior for the low Reynolds-number k,e-model in comparison to the other models. The conditions are slightly better predicted by the RSM than by the k,e-model. The ¯ow conditions seem to be more in¯uenced by the walls than by anisotropic turbulence eects and therefore the low Reynolds-number k,e-model gives better results than the RSM and the k,e-model. The grid is ®ne enough in the boundary layer to get signi®cant dierences between the low Reynolds-number k,e- and the standard k,emodel. 5.2. Combustion conditions Further numerical and experimental investigations were carried out during combustion in order to get more insight into the formation and decomposition processes of dierent gas components. The experimental velocity data shown in Fig. 8 were obtained during two reference burn cycles BC1 and BC2. Pro®le measurements were carried out during the main burning phase in a period between 20 and 30 min after reloading a new batch of wood logs. In Fig. 8 a comparison of calculated and measured velocities is shown for the isothermal investigations on dierent heights at z ÿ20, 40, 100 and 180 mm. The characteristics of the measured velocity pro®les are qualitatively well predicted. The maximum of the mean w-velocities on the ®rst three measuring levels is
Fig. 5. Comparison of measured and calculated ¯ow ®elds at isothermal conditions.
H. Knaus et al. / Experimental Thermal and Fluid Science 21 (2000) 99±108
105
Fig. 6. Comparison of dierent grids at isothermal conditions: z ÿ20, z 50 mm.
Fig. 7. Comparison of dierent turbulence models at isothermal conditions: z ÿ20, z 50, z 100 and z 180 mm.
underpredicted by both the k,e-model and the RSM. As for the isothermal case the maximum level of the w-velocity component is much better predicted using the low Reynolds-number k,e-model. At a height of z 180 mm
a good agreement of measurements and simulation results is achieved for all turbulence models. In a further step the modelling studies for combustion conditions were used for obtaining information
106
H. Knaus et al. / Experimental Thermal and Fluid Science 21 (2000) 99±108
Fig. 8. Dierent turbulence models at combustion: z ÿ20, z 40, z 100 and z 180 mm.
Fig. 9. Calculated mixture fraction, temperature and CO concentration in the burnout zone.
about the mixing intensity of the burnout air distribution by the slot and the three nozzles. Fig. 9 shows results of the calculated mixing ratio between combustible gases and burnout air. The penetration depth of the burnout air injected by the nozzles is approximately half of the combustion chamber width revealing locally high mixing intensity. The lower momentum of the burnout air by the slot which is due to the relatively large outlet area results in a lower penetration depth and a lower mixing intensity. The injected air is de¯ected to the left-
hand side of the burnout zone without complete mixing. Therefore, a streak is formed on the right-hand side of the burnout zone, retaining approximately the initial concentration. The characterisation of the combustion process within the burnout zone of the test stove is done by studying the measured and calculated gas concentration and temperature ®elds. In Fig. 9 contour plots are shown for the computed CO-concentrations and temperatures using the low Reynolds-number k,e-turbulence
H. Knaus et al. / Experimental Thermal and Fluid Science 21 (2000) 99±108
107
Fig. 10. Calculated and measured CO concentrations and temperatures in the burnout zone.
model. Additionally, the positions for gas and temperature measurements are depicted. For these positions, measured values are opposed to the computed values in Fig. 10. The CO-concentration ®eld, computed in the middle of the burnout zone (Fig. 9) shows a rapid reduction of the initially high values due to the injection of secondary air through the slot and the nozzles and the high temperatures within this part of the burnout zone. Related to the asymmetric arrangement of the air distribution and therefore to the insucient mixing intensity the high CO-content of the streak at the right-hand side of the burnout chamber cannot be reduced while passing the throat of the burnout zone. In the upper part of the burnout zone, lower temperatures and a reduced mixing intensity prevent a further oxidation of remaining CO. The comparison of the numerical simulations with measurement data shows that the characteristics of the computed gas concentration and temperature ®elds are in good accordance with the measurement data. The ¯uctuations of measured temperatures at the dierent measuring positions are due to changes in the ¯ow ®eld and combustion conditions in the burnout zone depending on the burn course and on dierent boundary conditions. Dierences between measured and calculated data, especially for the temperature measurements, can be put down to the fact that the predicted size of the recirculation zone formed in the upper left part of the burnout zone is slightly too small. The aected measuring points are the positions 1, 2, and 3, which are also the positions where the biggest dierences between measured and computed temperatures are found. The reaction of CO is predicted as too slow in the lower part of the burnout zone.
6. Summary The turbulent ¯ow and combustion processes within the burnout zone of a commercially available tile stove heating insert have been investigated. Numerical modelling studies based on detailed experimental boundary conditions were carried out in order to characterise the
¯ow ®eld at isothermal conditions as well as the gas concentration, temperature and ¯ow ®eld in the reaction zone for combustion of natural wood logs. Three different turbulence models, namely the k,e-, the low Reynolds-number k,e- and the Reynolds-stress model have been used. A good agreement between the computed results and measurement data was found especially for the low Reynolds-number k,e-model. However, the grid nodes in the boundary layer are not well enough distributed for a good resolution in this area. The results of the low Reynolds-number k,e-model might be further improved using a grid adaption technique. The use of higher-order discretisation schemes could also contribute to a further improvement. To get better predictions of temperatures and concentration of the species a more detailed reaction scheme instead of the simple global reaction scheme used would be advantageous. The processes of ¯uid ¯ow and combustion are well captured by the numerical model. Therefore it can be used as an engineering tool for optimising the air distribution arrangement and to improve burnout and eciency in order to reduce the environmental impact of the investigated test stove. Acknowledgements This work is partly funded by the European Commission DG XII within the non-nuclear JOULE-III programme under the contract JOR3-CT95-0056. References [1] S. Unterberger, H. Knaus, J. Sutinen, H.G. Heller, H. Maier, U. Schnell, K.R.G. Hein, Optimisation of the combustion and emission behaviour of wood burning systems ± An experimental and numerical approach, in: 10th European Conference and Technology Exhibition on Biomass for Energy and Industry, W urzburg, Germany, June 1998, pp. 1456±1459. [2] U. Schnell, R. Schneider, H.-C. Magel, B. Risio, J. Lepper, K.R.G. Hein, Numerical simulation of advanced coal-®red combustion systems with in-furnace NOx control technologies, in: Third International Conference on Combustion Technologies for a Clean Environment, Lisbon, Portugal, 1995, pp. 336±343.
108
H. Knaus et al. / Experimental Thermal and Fluid Science 21 (2000) 99±108
[3] J.P. Van Doormaal, G.D. Raithby, Enhancements of the simple method for predicting incompressible ¯ows, Numerical Heat Transfer 7 (1984) 147±163. [4] C.M. Rhie, W.L. Chow, Numerical study of the turbulent ¯ow past an airfoil with trailing edge separation, AIAA Journal 21 (11) (1983) 1525±1532. [5] R. Schneider, B. Risio, J. G orres, U. Schnell, K.R.G. Hein, Application of a dierential Reynolds-stress turbulence model to the numerical simulation of coal-®red utility boilers, in: Third International Symposium on Coal Combustion, Beijing, China, 1996, pp. 1±11. [6] V.L. Zimont, Y.M. Trushin, Total combustion kinetics of hydrocarbon fuels, Combustion, Explosion, and Shock Waves 5 (4) (1969) 391±394. [7] J.B. Howard, G.C. Williams, D.M. Fine, Kinetics of carbon monoxide oxidation in Institute, Pittsburg, PA, 1972, pp. 975±986. [8] B.F. Magnussen, The Eddy dissipation concept, in: XI Task Leaders Meeting ± Energy Conversion in Combustion, IEA, 1989. [9] W.A. Fiveland, Three dimensional radiative heat transfer solutions by the discrete ordinates method, AIAA Journal of Thermophysics 2 (4) (1988) 309±316.
[10] I.A. Demirdzic, A Finite Volume Method for Computation of Fluid Flow in Complex Geometries, Ph.D. Thesis, Mech. Eng. Dept., Imperial College of Science and Technology, London, 1982. [11] M.A. Peric, A ®nite volume method for the prediction of threedimensional ¯uid ¯ow in complex ducts, Ph.D. Thesis, Mech. Eng. Dept., Imperial College of Science and Technology, London, 1985. [12] B.E. Launder, D.B. Spalding, The numerical computation of turbulent ¯ows, Computer Methods in Applied Mechanics and Engineering 3 (1974) 269±289. [13] C. Patel, W. Rodi, G. Scheurer, Turbulence models for near-wall and low Reynolds number ¯ows: A review, AIAA Journal 23 (9) (1984) 1308±1319. [14] C.K.G. Lam, K.A. Bremhorst, Modi®ed form of the k-e-model for predicting wall turbulence, Journal of Fluids Engineering 103 (1981) 456±460. [15] B.J. Daly, F.H. Harlow, Transport equations in turbulence, Physics of Fluids 13 (1970) 2634±2649. [16] S. Hogg, M.A. Leschziner, Second-moment-closure calculation of strongly swirling ¯ow with large density gradients, International Journal of Heat and Fluid Flow 10 (1) (1989) 16±27.
www.elsevier.nl/locate/etfs
On the application of dierent turbulence models for the computation of ¯uid ¯ow and combustion processes in small scale wood heaters H. Knaus, S. Richter, S. Unterberger *, U. Schnell, H. Maier, K.R.G. Hein Institute of Process Engineering and Power Plant Technology (IVD), University of Stuttgart, Pfaenwaldring 23, D-70550 Stuttgart, Germany Received 20 June 1999; received in revised form 15 November 1999; accepted 1 December 1999
Abstract In order to optimise the combustion process in small scale wood heaters with a thermal capacity of up to 15 kW using natural wood logs an important task is to investigate the mixing and combustion processes within the reaction zones of such ®ring systems. For a basic understanding of the formation and decomposition of dierent gas components and the dependency on mixing intensity and therefore on turbulent behaviour of the reacting ¯ow numerical studies were carried out using the simulation program AIOLOS. A comparison of results obtained during the numerical modelling studies applying the k,e-, low Reynolds-number k,e- and the Reynolds-stress turbulence model on dierent grids with detailed experimental data shows a good correspondence both for isothermal conditions as well as during combustion. The main aim is to validate the numerical model enabling its use as an engineering tool to investigate the in¯uence of parameters such as mixing conditions, combustion air distribution and furnace geometry on the combustion process within the wood heater on a complete combustion for reduced emission of unburnt components, such as CO, hydrocarbons and soot. Ó 2000 Elsevier Science Inc. All rights reserved. Keywords: Wood combustion; Modelling studies; Turbulence models; Homogeneous gas phase reactions
1. Introduction One way of reducing the global CO2 emission problem is to substitute wood for fossil fuels in heat and power generation to a certain extent, an option also appropriate for small scale combustion systems for domestic heating purposes. However, the combustion of natural wood logs in small scale stoves with a thermal capacity in the range of up to 15 kW for domestic heating purposes often involves higher emissions of CO, unburnt hydrocarbons and particles. It is estimated that in Germany, for example, the share of small scale wood combustion systems contributing to the emissions of CO, hydrocarbons and soot is between 16% and 40%, although the amount of the total energy consumption is only about 1%. In order to improve the emission behaviour of small scale wood combustion systems a basic understanding of the combustion process is essential. Facilitating the in*
Corresponding author. Tel.: +49-711-685-3572; fax: +49-711-6853491. E-mail address: [email protected] (S. Unterberger).
sight into the complex inter-linked phenomena of chemical reactions and turbulent ¯ow ®eld behaviours as well as for the investigation of parametric eects, the simulation program AIOLOS is used as an eective tool for analysis and further optimisation of the combustion process. The numerical modelling studies were based on boundary conditions determined by detailed pro®le measurements of gas concentrations, temperatures, velocities and turbulence intensities within the burnout zone of the selected test stove [1]. 2. Design of the test stove A tile stove heating insert with a thermal capacity of 10 kW representing the state-of-the-art was selected as a test stove. The design of the stove and the area investigated by the mentioned pro®le measurements is shown in Fig. 1. The commercially available stove consists of two spatially separated reaction zones. In the ®rst reaction zone, both the wood fuel is gasi®ed and the gases released by the wood logs start to combust. In the secondary reaction zone, preheated air is added to the gases for almost complete burnout. The injection of the
0894-1777/00/$ - see front matter Ó 2000 Elsevier Science Inc. All rights reserved. PII: S 0 8 9 4 - 1 7 7 7 ( 9 9 ) 0 0 0 5 9 - X
100
H. Knaus et al. / Experimental Thermal and Fluid Science 21 (2000) 99±108
Fig. 1. Investigated tile stove heating insert, thermal capacity 10 kW.
secondary air is realised in two ways: by means of three nozzles before the throat and by a slot in the middle of the throat.
3. Experimental investigations 3.1. Experimental set-up The experimental test ®eld consists of measurement devices for gas concentrations, temperatures and velocities. The gas analysis system is designed to measure the gas concentrations of the main components CO2 , CO and O2 in the ¯ue gas as well as in the burnout zone at the same time. Additionally H2 O and hydrocarbons can be measured either in the ¯ue gas or in the burnout zone. The suction probes for the pro®le measurements are arranged at dierent points in the secondary reaction zone. The arrangement of the measuring points is shown in Fig. 2. Pro®le measurements were carried out by measuring successively at the positions 1±11 for a duration of approximately 1 min per single probe. Gas temperature measurements were done by using Ni±CrNi
thermocouples which were installed inside the gas sampling suction probes. The measured gas temperatures were corrected according to a formula obtained by using a suction pyrometer. Wall temperatures were measured by thermocouples at two measuring points in the burnout combustion chamber as well as the temperature of the preheated secondary air before entering the nozzles and the slot. The Laser-Doppler Velocimeter (LDV) system used for the ¯ow ®eld characterisation consisted of a 2 W argon-ion laser used with a DANTEC/invent 2D FiberFlow system facilitating the simultaneous measurement of two perpendicular velocity components in x- and z-directions, respectively. 3.2. Determination of boundary conditions and pro®le measurements In a ®rst step the boundary conditions at the dierent inlets of the burnout zone were determined at isothermal and combustion conditions, respectively. Velocity measurements were carried out on several pro®les covering the total depth (270 mm) of the channel between the gasi®cation/combustion zone and the burnout zone as
Fig. 2. Measuring positions and locations for comparison of predictions and experimental data.
H. Knaus et al. / Experimental Thermal and Fluid Science 21 (2000) 99±108
well as in the outlet area of the slot and the three nozzles for obtaining the combustion air volume ¯ows added to the reacting main ¯ow. In contrast to the isothermal case, the determination of velocity boundary conditions during combustion was done by a reduced number of measuring points due to the non-stationary combustion process. In order to get reliable data during combustion, results of several dierent burn cycles were used to calculate mean values of velocities and turbulence intensities. In addition to the above-mentioned investigations by LDV, gas concentration measurements of CO2 , CO, O2 and hydrocarbons were done at the inlet of the burnout zone. Temperature measurements were carried out for the entering combustible gases, the added burnout air and at two wall sections in order to get complete data sets characterising the combustion conditions in the investigated test stove during the dierent burn phases mentioned above. Furthermore, velocity pro®le measurements were carried out at isothermal conditions and for combustion tests, respectively. Gas concentration and temperature pro®le measurements at dierent positions in the burnout zone according to Fig. 2 were carried out for the determination of combustion conditions during a complete burn cycle. 4. Numerical modelling studies The present calculations were performed with the ®nite-volume code AIOLOS [2]. This simulation program has been developed for the numerical calculation of stationary, turbulent reacting ¯ows. The SIMPLECmethod [3] is applied to compute the pressure. The pressure-weighted interpolation method (PWIM) [4] is used to prevent the decoupling of velocities and pressure on the non-staggered grid. The convection ¯uxes are approximated by the upwind or higher-order dierencing schemes [5]. In the following the physical models used for the simulation of the wood heater are shortly described. 4.1. Kinetic model of chemical reactions The combustible volatiles are converted by the following reactions {C1}±{C3}: hm ni a O2 fC1g : Cm Hn 2 4 hn ni n ! mCO 1 a H2 a H2 O 2 4 2 1 fC2g : CO O2 ! CO2 2
2
1 3 fC3g : H2 O2 ! H2 O: 2 The distribution of water and hydrogen as products of the hydrocarbon combustion described by reaction {C1} is calculated from the water-gas-shift equilibrium.
101
All reactions are treated as irreversible reactions and the kinetic rates of the reactions {C1}±{C3} are taken from [6,7]. 4.2. Eddy dissipation concept (EDC) The EDC is a general concept for treating the interaction between turbulence and chemistry in ¯ames [8]. The method is based on a detailed description of the dissipation of turbulent eddies. In the EDC the total space is subdivided into a reaction space, called the `®ne structures' and the surrounding ¯uid. In the presented reaction scheme the reactions {C1}±{C3} are treated as taking place only in these ®ne structures which means that the reactions only take place in the smallest turbulent length scales. 4.3. Thermal radiation The discrete ordinates method in an S4-approximation [9] is used to solve the radiation transport equation. Since the intensity of radiation depends on absorption, emission and scattering characteristics of the medium passed through, a detailed representation of the radiative properties of a gas mixture would be very complex and currently beyond the scope of a 3D-code for the simulation of industrial combustion systems. Thus, contributing to the numerical eciency, simpli®cations are introduced, even at the loss of some accuracy. The absorption coecient of the gas phase is assumed to have a constant value of 0.2/m. 4.4. Transport equations and turbulence modelling In the following the conservation equations of mass, momentum, scalar quantities and turbulence quantities and their discretisation are described. All formulas and modelling assumptions presuppose high Reynoldsnumbers and steady-state ¯ow conditions. It is assumed that the ¯ow ®eld is weakly compressible which means that the density depends only on temperature and ¯uid composition but not on pressure. 4.5. Finite volume discretisation A body-®tted grid is used to represent the slope of the furnace and especially of the burner quarls exactly. Therefore the conservation equations have to be formulated in general curvilinear co-ordinates. Regarding general curvilinear co-ordinates, a huge number of different formulas exist, in particular for the momentum equation [10]. In this paper the following equations proposed by Peric [11] are used, which provide a strong conservative form of the conservation equations: Continuity: o Uj q 0: onj
4
102
H. Knaus et al. / Experimental Thermal and Fluid Science 21 (2000) 99±108
Momentum: o l oui j ouk m j Uj quj ÿ Bm bi bk pbji onm onj J onm bjk qu0i u0k J qgi : Scalar: o C/ o/ j Uj q/ ÿ B Js/ : J onm m onj
o leff oe j e Uj qe ÿ B Ce1 Pkk ÿ J qCe2 e; re J onm m onj k 5
6
The contravariant velocity components Uj and the tensor Bij are de®ned by Uj J bjm um ;
Bij bjk bik ;
7
where J denotes the Jacobian determinant and the elements of the inverse Jacobian matrix bij represent the coordinate transformation of Cartesian into general curvilinear co-ordinates bij
oxi : oyj
8
The components uj of the velocity vectors are expressed in terms of the spatially invariant Cartesian basis. 4.6. k,e-turbulence model For the k,e-model a ®rst-order turbulence model closure according to Boussinesq is used. In analogy to laminar ¯ows the Reynolds-stresses are assumed to be proportional to the gradients of the mean velocities l oui n oum 1 2 ouk k 2 bm bi ÿ dmi bk ÿ dmi qk 9 smi t J onn on1 onk 3 3 introducing lt as a proportional factor of the turbulent viscosity. The model is based on the Prandtl±Kolmogorov approach mt Cm k 1=2 L
10
and the following equation for e, resulting from a dimensional analysis: k 3=2 11 e Ce L where L and k 1=2 are the characteristic length scales and velocities of turbulence. Eqs. (10) and (11) result in an equation for the turbulent viscosity dependent on the turbulent kinetic energy k and the turbulent dissipation e k2 mt q; leff lt l1 ; 12 e with the dynamic viscosity ll . For the turbulent kinetic energy and the turbulent dissipation, transport equations can be developed from the Navier±Stokes equations assuming an isotropic turbulence: o l ok j Uj qk ÿ eff Bm Pkk ÿ J qe; 13 rk J onm onj lt Cl q
P ÿlt
oui n ouj m b b onn j onm i
oui n b onn j
ÿ
14
lt oq i b q2 oni j
op i b; oni j
15
where P is the term for the production of turbulence by gradients of the mean velocities. The parameters used in the turbulence model equations are chosen according to Launder and Spalding [12]. 4.7. Low Reynolds-number k,e-turbulence model The standard k,e- and the Reynolds-stress model are restricted to high Reynolds-number applications. In the test stove investigated, however, Reynolds-numbers might become low, especially in the upper part of the burnout section. Moreover, the in¯uence of the walls on free ¯ow characteristics may not be adequately modelled by standard wall functions as used for the k,e-model and the Reynolds-stress model in such narrow geometries. The k,e-model has thus been extended to enable its use at low Reynolds-numbers and to describe more accurately the ¯ow close to a solid wall. Based on the ®ndings of Patel et al. [13] the approach of Lam and Bremhorst [14] was selected from the various low Reynolds-number models available. The Lam and Bremhorst model combines good performance with simple numerical formulation and it diers from the standard k,e-model by the inclusion of the following functions f to modify the constants C. Turbulent viscosity equation (12): ClLRM Cl fl :
16
Turbulent dissipation rate equation (14): LRM Ce1
LRM Ce2 Ce2 f2 :
Ce1 f1 ;
Reynolds-number functions: 2 20:5 1 ; fl 1 ÿ exp ÿ 0:0165 Ry RT f1 1 k2 q ; RT le
0:05 fl
3 ;
f2 1 ÿ exp ÿR2T ;
p kyq ; Ry l
17
18
19
20
where y is the distance normal to the wall. At high Reynolds-numbers RT or Ry , the functions ``f'' tend to unity and the turbulent quantities are thus determined analogously to the standard k,e-model. For the low Reynolds-number model the transport equations of k and e are integrated up to the wall attempting to model relaminarisation eects by the functions f described above. In Fig. 3, the functions fl , f1 and f2 are
H. Knaus et al. / Experimental Thermal and Fluid Science 21 (2000) 99±108
oui m ouj m 0 0 0 0 Pij ÿ quj uk b qui uk b onm k onm k lt oq k op k oq k op k b b ÿ b b : ÿ qJ onk j onk i onk j onk j
Fig. 3. Low Reynolds-number functions f1 , f2 and fl over the wall distance y .
depicted over the wall distance y de®ned by the following equation: y
Cl qk 1=2 y : l1
21
It can be seen that the functions f1 and f2 only have a strong in¯uence for wall distances y less than 5, whereas the function fl has to be considered in the whole boundary layer. Function f2 has a steep gradient in a wall distance y between 0 and 5 and for this reason the grid must be extremely ®ne especially in this part of the boundary layer to get a sucient resolution. 4.8. Dierential Reynolds-stress model (RSM) Flows with non-isotropic turbulence cannot adequately be described by a two-equation model like the k,e-model. Under such ¯ow conditions turbulence should be described by dierential transport equations for all Reynolds-stress components o 2 Uk qu0i u0j Dij Pij Wij ÿ J dij qe; onk 3
24
The expression in the ®rst brackets describes the interaction between the stresses and the rate of the mean strain. These terms are exact and need not be modelled. The expressions in the second brackets arise only under non-uniform density conditions and originate from Favre-averaging. Interactions between pressure and strain leading to a redistribution process within the Reynoldsstresses towards an isotropic state [16] are represented by the term Wij and are modelled as follows: e 2 1 qu0i u0j ÿ dij qk ÿ C2 Pij ÿ dij Pkk ; Wij ÿ JC1 k 3 3 25 where Wij consists of two distinct parts which describe the relations between purely turbulent quantities and between turbulence and mean strain. The dissipation rate e of the Reynolds-stresses is modelled using transport equation (14) as described for the standard k,e-model. It implies the isotropy of the dissipation process at smallest length scales. A summary of the constants appearing in the above equations is given in [5].
5. Comparison of experimental results and simulation studies Two body-®tted grids with dierent ®neness, one with 60 000 and the other with 200 000 nodes, were used to model the complex geometry of the secondary reaction zone. In Fig. 4 the coarse grid is shown. The ®ne
22
where Dij , Pij and Wij represent the diusion, production and pressure±strain processes, respectively. The diusion of the Reynolds-stresses Dij is approximated by ! o leff ou0i u0j k B : Dij 23 onk rk J onm m This expression is strongly simpli®ed compared to the generalised gradient-diusion hypothesis of Daly and Harlow [15]. The replacement of the detailed non-isotropic viscosity by the isotropic viscosity drastically reduces the complexity of the diusion term and promotes the numerical stability of the model. The generation term of the Reynolds-stress transport equations is given by
103
Fig. 4. Coarse computational grid.
104
H. Knaus et al. / Experimental Thermal and Fluid Science 21 (2000) 99±108
grid oers a sucient resolution of the boundary layer needed for the low Reynolds-number k,e-model. For most of the cell centres adjacent to the wall y is less than 5 and the largest distance is 40 for isothermal as well as for combustion conditions. It was not possible, however, to put a constant number of cells in the laminar boundary sublayer because of the complex geometry and the complex ¯ow situation within the reaction zone. 5.1. Isothermal conditions Because of the non-stationary processes of wood log combustion and the resulting measuring problems the ®rst simulations were done without combustion. Results of the k,e-model and the RSM are compared on the two grids to judge the grid dependency of the solutions. Comparing the measured and computed ¯ow ®eld a qualitatively good agreement is observed. For the simulation using the coarse grid and the k,e-model the location of the centre of the recirculation zone is well predicted (Fig. 5), but it is striking that this zone in the upper part of the furnace is more noticeably shaped in the measurements. On the ®ne grid using the k,e-model the recirculation zone is predicted as too strong and its centre is located more in the middle of the furnace. The best results are given by the RSM. The location of the recirculation zone as well as its strength is qualitatively in good accordance with the measurements. To validate the model quantitatively additional comparisons of velocity measurements and simulation results were done along dierent lines which are depicted in Fig. 2. Looking at the height of z ÿ20 mm the prediction of the w-velocity component is good for both grids using the k,e-model and the RSM (Fig. 6). However, on the ®ne grid the peak at x 10 mm is better predicted by the RSM than by the k,e-model. At a height of z 50 mm the k,e-model can predict neither on the coarse nor on
the ®ne grid as good as the RSM. Especially on the ®ne grid the shape and the location of the w-velocity maximum agrees well with the measurements for the RSM. Because of these results it was decided to continue with the simulations only on the ®ne grid and additionally to implement the low Reynolds-number k,emodel. This model turned out to give the best predictions. At z ÿ20 mm the simulation results of all turbulence models are very similar (Fig. 7). However the predictions of the magnitude and the gradient of the w-velocity at z 50, 100 and 180 mm are superior for the low Reynolds-number k,e-model in comparison to the other models. The conditions are slightly better predicted by the RSM than by the k,e-model. The ¯ow conditions seem to be more in¯uenced by the walls than by anisotropic turbulence eects and therefore the low Reynolds-number k,e-model gives better results than the RSM and the k,e-model. The grid is ®ne enough in the boundary layer to get signi®cant dierences between the low Reynolds-number k,e- and the standard k,emodel. 5.2. Combustion conditions Further numerical and experimental investigations were carried out during combustion in order to get more insight into the formation and decomposition processes of dierent gas components. The experimental velocity data shown in Fig. 8 were obtained during two reference burn cycles BC1 and BC2. Pro®le measurements were carried out during the main burning phase in a period between 20 and 30 min after reloading a new batch of wood logs. In Fig. 8 a comparison of calculated and measured velocities is shown for the isothermal investigations on dierent heights at z ÿ20, 40, 100 and 180 mm. The characteristics of the measured velocity pro®les are qualitatively well predicted. The maximum of the mean w-velocities on the ®rst three measuring levels is
Fig. 5. Comparison of measured and calculated ¯ow ®elds at isothermal conditions.
H. Knaus et al. / Experimental Thermal and Fluid Science 21 (2000) 99±108
105
Fig. 6. Comparison of dierent grids at isothermal conditions: z ÿ20, z 50 mm.
Fig. 7. Comparison of dierent turbulence models at isothermal conditions: z ÿ20, z 50, z 100 and z 180 mm.
underpredicted by both the k,e-model and the RSM. As for the isothermal case the maximum level of the w-velocity component is much better predicted using the low Reynolds-number k,e-model. At a height of z 180 mm
a good agreement of measurements and simulation results is achieved for all turbulence models. In a further step the modelling studies for combustion conditions were used for obtaining information
106
H. Knaus et al. / Experimental Thermal and Fluid Science 21 (2000) 99±108
Fig. 8. Dierent turbulence models at combustion: z ÿ20, z 40, z 100 and z 180 mm.
Fig. 9. Calculated mixture fraction, temperature and CO concentration in the burnout zone.
about the mixing intensity of the burnout air distribution by the slot and the three nozzles. Fig. 9 shows results of the calculated mixing ratio between combustible gases and burnout air. The penetration depth of the burnout air injected by the nozzles is approximately half of the combustion chamber width revealing locally high mixing intensity. The lower momentum of the burnout air by the slot which is due to the relatively large outlet area results in a lower penetration depth and a lower mixing intensity. The injected air is de¯ected to the left-
hand side of the burnout zone without complete mixing. Therefore, a streak is formed on the right-hand side of the burnout zone, retaining approximately the initial concentration. The characterisation of the combustion process within the burnout zone of the test stove is done by studying the measured and calculated gas concentration and temperature ®elds. In Fig. 9 contour plots are shown for the computed CO-concentrations and temperatures using the low Reynolds-number k,e-turbulence
H. Knaus et al. / Experimental Thermal and Fluid Science 21 (2000) 99±108
107
Fig. 10. Calculated and measured CO concentrations and temperatures in the burnout zone.
model. Additionally, the positions for gas and temperature measurements are depicted. For these positions, measured values are opposed to the computed values in Fig. 10. The CO-concentration ®eld, computed in the middle of the burnout zone (Fig. 9) shows a rapid reduction of the initially high values due to the injection of secondary air through the slot and the nozzles and the high temperatures within this part of the burnout zone. Related to the asymmetric arrangement of the air distribution and therefore to the insucient mixing intensity the high CO-content of the streak at the right-hand side of the burnout chamber cannot be reduced while passing the throat of the burnout zone. In the upper part of the burnout zone, lower temperatures and a reduced mixing intensity prevent a further oxidation of remaining CO. The comparison of the numerical simulations with measurement data shows that the characteristics of the computed gas concentration and temperature ®elds are in good accordance with the measurement data. The ¯uctuations of measured temperatures at the dierent measuring positions are due to changes in the ¯ow ®eld and combustion conditions in the burnout zone depending on the burn course and on dierent boundary conditions. Dierences between measured and calculated data, especially for the temperature measurements, can be put down to the fact that the predicted size of the recirculation zone formed in the upper left part of the burnout zone is slightly too small. The aected measuring points are the positions 1, 2, and 3, which are also the positions where the biggest dierences between measured and computed temperatures are found. The reaction of CO is predicted as too slow in the lower part of the burnout zone.
6. Summary The turbulent ¯ow and combustion processes within the burnout zone of a commercially available tile stove heating insert have been investigated. Numerical modelling studies based on detailed experimental boundary conditions were carried out in order to characterise the
¯ow ®eld at isothermal conditions as well as the gas concentration, temperature and ¯ow ®eld in the reaction zone for combustion of natural wood logs. Three different turbulence models, namely the k,e-, the low Reynolds-number k,e- and the Reynolds-stress model have been used. A good agreement between the computed results and measurement data was found especially for the low Reynolds-number k,e-model. However, the grid nodes in the boundary layer are not well enough distributed for a good resolution in this area. The results of the low Reynolds-number k,e-model might be further improved using a grid adaption technique. The use of higher-order discretisation schemes could also contribute to a further improvement. To get better predictions of temperatures and concentration of the species a more detailed reaction scheme instead of the simple global reaction scheme used would be advantageous. The processes of ¯uid ¯ow and combustion are well captured by the numerical model. Therefore it can be used as an engineering tool for optimising the air distribution arrangement and to improve burnout and eciency in order to reduce the environmental impact of the investigated test stove. Acknowledgements This work is partly funded by the European Commission DG XII within the non-nuclear JOULE-III programme under the contract JOR3-CT95-0056. References [1] S. Unterberger, H. Knaus, J. Sutinen, H.G. Heller, H. Maier, U. Schnell, K.R.G. Hein, Optimisation of the combustion and emission behaviour of wood burning systems ± An experimental and numerical approach, in: 10th European Conference and Technology Exhibition on Biomass for Energy and Industry, W urzburg, Germany, June 1998, pp. 1456±1459. [2] U. Schnell, R. Schneider, H.-C. Magel, B. Risio, J. Lepper, K.R.G. Hein, Numerical simulation of advanced coal-®red combustion systems with in-furnace NOx control technologies, in: Third International Conference on Combustion Technologies for a Clean Environment, Lisbon, Portugal, 1995, pp. 336±343.
108
H. Knaus et al. / Experimental Thermal and Fluid Science 21 (2000) 99±108
[3] J.P. Van Doormaal, G.D. Raithby, Enhancements of the simple method for predicting incompressible ¯ows, Numerical Heat Transfer 7 (1984) 147±163. [4] C.M. Rhie, W.L. Chow, Numerical study of the turbulent ¯ow past an airfoil with trailing edge separation, AIAA Journal 21 (11) (1983) 1525±1532. [5] R. Schneider, B. Risio, J. G orres, U. Schnell, K.R.G. Hein, Application of a dierential Reynolds-stress turbulence model to the numerical simulation of coal-®red utility boilers, in: Third International Symposium on Coal Combustion, Beijing, China, 1996, pp. 1±11. [6] V.L. Zimont, Y.M. Trushin, Total combustion kinetics of hydrocarbon fuels, Combustion, Explosion, and Shock Waves 5 (4) (1969) 391±394. [7] J.B. Howard, G.C. Williams, D.M. Fine, Kinetics of carbon monoxide oxidation in Institute, Pittsburg, PA, 1972, pp. 975±986. [8] B.F. Magnussen, The Eddy dissipation concept, in: XI Task Leaders Meeting ± Energy Conversion in Combustion, IEA, 1989. [9] W.A. Fiveland, Three dimensional radiative heat transfer solutions by the discrete ordinates method, AIAA Journal of Thermophysics 2 (4) (1988) 309±316.
[10] I.A. Demirdzic, A Finite Volume Method for Computation of Fluid Flow in Complex Geometries, Ph.D. Thesis, Mech. Eng. Dept., Imperial College of Science and Technology, London, 1982. [11] M.A. Peric, A ®nite volume method for the prediction of threedimensional ¯uid ¯ow in complex ducts, Ph.D. Thesis, Mech. Eng. Dept., Imperial College of Science and Technology, London, 1985. [12] B.E. Launder, D.B. Spalding, The numerical computation of turbulent ¯ows, Computer Methods in Applied Mechanics and Engineering 3 (1974) 269±289. [13] C. Patel, W. Rodi, G. Scheurer, Turbulence models for near-wall and low Reynolds number ¯ows: A review, AIAA Journal 23 (9) (1984) 1308±1319. [14] C.K.G. Lam, K.A. Bremhorst, Modi®ed form of the k-e-model for predicting wall turbulence, Journal of Fluids Engineering 103 (1981) 456±460. [15] B.J. Daly, F.H. Harlow, Transport equations in turbulence, Physics of Fluids 13 (1970) 2634±2649. [16] S. Hogg, M.A. Leschziner, Second-moment-closure calculation of strongly swirling ¯ow with large density gradients, International Journal of Heat and Fluid Flow 10 (1) (1989) 16±27.