Numerical Modeling of a Discontinuous ... - ACS Publications

Ind. Eng. Chem. Res. 2005, 44, 3159-3177

3159

Numerical Modeling of a Discontinuous Incineration Process with On-line Validation Davide Manca* and Maurizio Rovaglio Dipartimento di Chimica, Materiali e Ingegneria Chimica “G. Natta”, Politecnico di Milano, Piazza Leonardo da Vinci n.32, I-20133 Milano, Italy

This paper describes the dynamic model of the hot section of an incineration plant with steam production. The model is based on a first-principles approach that starting from material, energy, and momentum balances leads to the integration of a differential-algebraic system of 132 equations. The hot section of the plant comprises a furnace, a postcombustion chamber, and a waste heat boiler. The reciprocating grate kiln represents one of the most important features of the process, where the heterogeneous combustion of waste takes place. The waste, coming from diverse collections, is transferred from the pit into a vertical hopper and enters into the furnace pushed by a discontinuous feeding grate. Also the burning waste is discontinuously moved through the primary combustion chamber by reciprocating grates. These discontinuous features are mainly responsible for the high oscillations of process variables. The model was implemented for control purposes with the aim of reducing oscillations and offsets through a model-based control system. The paper describes thoroughly the key points that were addressed in structuring and solving the numerical detailed model. Finally, a validation procedure sketches a quantitative comparison between experimental measured data and simulated values. The measured process variables were acquired on-line from the real-time database running on the distributed control system. 1. Introduction

Table 1. Final Destination of Municipal Wastes in Europe (Cagiano de Azevedo, 2002)

Waste disposal represents a major problem for modern society. The exponential increase of production, technology, and commodities in the past 30 years has produced, as a consequence, a parallel increase of exhaust, old, unused scraps. Packaging dimensions have also increased with a direct effect on the quantity of wastes. In the past, wastes were mainly disposed of in landfills and municipal dumps. In contrast, incineration is nowadays growing in importance due to environmental and economic reasons. Table 1 reports an overview of the dichotomy between existing incinerators and municipal dumps in European countries. It is possible to observe that incineration is prevailing over dumping in some countries, although dumping is still dominant in most nations. This is mainly due to political choices made in past decades. Nonetheless, waste incineration is becoming the dominant option when new plants are involved. The reduction of environmental impact and the production of electric power represent the main reasons for this wisdom reform. Italy is slowly moving toward the incineration option. Several incinerators are being designed and built up throughout the Italian territory. Latest data from FederAmbiente, the Italian environmental federation, show that 35 plants, having incineration capacities from 30000 to 300000 t/y, burn about 2400000 t/y of solid wastes.1 More stringent specifications for correct operation, adopted throughout the European community countries, call for a proper optimal control of those units. On the economic side, the plant efficiency must be stressed to * To whom correspondence should be addressed. Mailing address: CMIC Department - Politecnico di Milano, 20133 Milano, Italy. Phone: +39 02 23993271. Fax: +39 02 70638173. E-mail: [email protected].

European countries

municipal waste production 103 t/y

Austria Belgium Denmark Finland France Germany Greece Ireland Italy Luxembourg Holland Portugal Sweden United Kingdom

2509 1646 2610 2100 20000 47392 3200 1550 27000 218 8430 3500 3200 20000

incineration 103 t/y % 410 499 1466 50 9759 12468 1 0 1400 126 2192 0 1300 2500

16.3 30.3 56.2 2.4 48.8 26.3 0.0 0.0 5.2 57.8 26.0 0.0 40.6 12.5

municipal dump 103 t/y % 1381 903 586 1500 10025 33197 2970 1432 24000 83 2870 3080 1200 14000

55.0 54.9 22.5 71.4 50.1 70.0 92.8 92.4 88.9 38.1 34.0 88.0 37.5 70.0

increase thermal efficiency and power production. Revenues come from both the quantity of burnt wastes and the steam production. Steam generates electric power in the turbine engine. Exhaust steam can also be used in cogeneration plants for civil purposes (winter heating and hot water). These articulate profiles about environment and economics must meet together in the plant operation, which is supervised by a control unit. Online control of heterogeneous combustion processes can be a challenge due to intrinsic problems such as inverse response, time delay, hot spots, and unknown feeding conditions as extensively described in refs 2-5. Conventional control loops based on multiple input and multiple output (MIMO) techniques are ineffective in solving such issues. A specifically developed system, based on advanced control algorithms, is the winning answer. The family of model-based control algorithms represents an effective solution to this problem. Model

10.1021/ie0495473 CCC: $30.25 © 2005 American Chemical Society Published on Web 03/24/2005

3160

Ind. Eng. Chem. Res., Vol. 44, No. 9, 2005

predictive control (MPC) is one of the most diffused techniques. It is based on the availability of a numerical dynamic model of the process that is able to predict the plant behavior, subject to some manipulated variables. This paper deals with the modeling of a heterogeneous combustion process for control purposes. The hot section of an incineration plant is described, analyzed, and dynamically modeled. The MPC has been applied to a municipal incineration plant in the northern part of Italy comprising two equivalent process lines capable of burning a total of 55000 t/y of solid waste. The attention will be focused on the dynamic modeling of the process only. A future work will deal with the MPC structure, on-line implementation, and testing of such an advanced control technique. The most peculiar aspect of this paper is the numerical modeling of a primary combustion chamber. Since it is a reciprocating grate kiln, the discontinuous feeding and movement of the waste, throughout the chamber, are the main features that are discussed and solved. The approach to the problem is based on first-principles equations. The model validation is performed through the acquisition of on-line process data from a distributed control system (DCS) device. Once a detailed dynamic model of the process is developed, it can be reduced in dimension and complexity or even used for identification purposes. Actually, model-based control consists of an optimization procedure that queries the dynamic model for prediction purposes. If the number of predictions is relevant, the detailed numerical model will probably not be able to keep pace with control timings. In this case, a more efficient model is required. Such a model can be obtained from the detailed one by reduction in complexity. Equivalently, it is possible to identify either a linear or nonlinear simplified model (ARX, ARMAX, NARX;6 ANN7) using the detailed one as a substitute for the real process. One might argue that the real combustion process could be directly used to perform a step-test procedure for on-line identification purposes. Unfortunately, when waste incineration processes are involved, it is not practical to perform step-tests on input variables such as waste flow rate, heat of combustion, and composition since the measurement and characterization of waste are not feasible on-line. At the same time, some step-tests performed on process variables, such as number of strokes, primary and secondary air flow rates, and air splitting under each grate, would be confused with other side effects produced by undesired and unforeseen disturbances (waste density, waste composition, waste heat of combustion; blower fluctuations; flow rate flickering; air leakages; oscillating pressure drops). One of the most important attributes of an incineration plant is the geometry of the primary combustion chamber. There are three main chamber layouts: drum kiln, grate kiln, and rotary kiln. When municipal and hospital wastes are involved, the choice usually goes to either a drum or a grate kiln. When either industrial toxic waste or sludge is available, a rotary kiln is often recommended.8-10 With reference to drum and grate kilns, the most significant difference is waste movement along the combustion bed. In the former, movement is almost continuous while in the latter waste advances discontinuously. The characterization and numerical modeling of the discontinuous motion of waste are the key points of this paper. An on-line validation of the dynamic model is also proposed and discussed. Actually,

Figure 1. Hot section of the incineration plant.

the operation of grate kilns is an unsteady process. No steady state can be reached, even theoretically, due to discontinuous pushes produced by reciprocating grates. 2. The Numerical Model With reference to Figure 1 the hot section of the municipal incinerator, which has been studied, modeled, and validated, comprises three distinct units: a primary kiln, a postcombustion chamber, and a waste heat boiler. The plant burns only solid waste of different typologies. The nominal capacity of the plant is 55000 t/y split into two parallel lines. An orange-peel bucket collects the waste from the pit and throws it into a loading hopper. Waste is then pushed into the primary kiln by a feeder that operates discontinuously. Four descending and reciprocating grates move the burning waste forward. The primary combustion air is fed through the holes of the reciprocating grates, while the secondary combustion air enters over the waste bed. The hot gases pass to the postcombustion chamber where oxidation is completed and finally arrive at the boiler. Feedwater is preheated in the economizer (see Figure 1) and then vaporizes in the cross-flow tube banks. Finally, from the cylindrical drum, steam goes to the superheater before reaching the turbine to produce electric power. 2.1. Primary Kiln. The heterogeneous combustion of solid waste over the grates of the primary kiln produces hot gases that continue oxidizing in the upper homogeneous volume before passing to the postcombustion chamber. The multifaceted profile of heterogeneous combustion, where exsiccation, ignition, burning, and extinction are spatially sequential, suggests adopting a cascade of perfectly mixed reactors: one for each of the reciprocating grates. In contrast, the hot gases, produced by waste combustion, are collected to a homogeneous volume where the

Ind. Eng. Chem. Res., Vol. 44, No. 9, 2005 3161

three distinct temperature values spatially distributed in the upper homogeneous section of the furnace. Consequently, the cascade of perfectly mixed heterogeneous reactors is characterized by a system of separate material balance equations while an overall heat balance equation is enough for the whole furnace. 2.1.1. Heterogeneous Combustion. The material balance equations for the solid phase on each reciprocating grate of the primary kiln can be expressed as

dMW,i in out ) FW,i - FW,i - RW,i i ) 1,...,NG dt

(1)

dMI,i in out ) FI,i - FI,i i ) 1,...,NG dt

(2)

in out in out FW,i ) FW,i-1 , FI,i ) FI,i-1 i ) 2,...,NG

(3)

in FW,1 ) FW(1 - ωI,0 - ωM,0)

(4)

in ) FWωI,0 FI,1

(5)

with

Figure 2. Schematization of the primary combustion chamber.

complex fluid dynamics produced by secondary air inlets recalls a perfectly mixed reactor. Figure 2 shows the modeling approach and spatial discretization of the primary combustion chamber. Besides primary and secondary air flow rates, some air leakages must also be accounted for. Main air leakages are present across the moving grates and are due to a permanent depression imposed on the furnace to avoid puffs of smokes. A fraction of these air leakages participates in the heterogeneous combustion of waste over each grate, while the complement bypasses the grates and directly arrives at the above homogeneous region as if it entered with a secondary air flow rate. Actually, the nonideal mixing conditions within the homogeneous combustion volume can be accounted for by introducing an additional bypass fraction of the gas flow rate produced by heterogeneous combustion. In Figure 2, flow rate #9 shows how a portion of hot gases bypasses the combustion volume without further oxidation. This assumption allows numerically modeling of the experimental presence of carbon monoxide (CO) entering the postcombustion chamber. We assume that CO in stream #8 is oxidized to CO2, while CO in stream #9 remains unoxidized. Both kinetics and reaction schemes will be discussed in the following. Furnace temperature is assumed constant throughout the chamber. This assumption, although quite strong, is reasonable and finds good agreement with the experimental thermocouple data that measure on-line

Equations 4 and 5 suggest that an elementary characterization of waste is measured experimentally or is numerically evaluated on-line by means of a Reconciliation procedure10,11 in terms of inert and moisture fractions as well as C,H,S,O,Cl,N mass fractions. Moreover, eq 4 assumes that the waste exsiccation takes place rapidly over the first grate. This hypothesis can be justified experimentally by observing the flames produced over that grate. Actually, waste first exsiccates and then ignites and starts burning. Given the elementary characterization of waste, an overall oxidation mechanism can be envisaged to describe the heterogeneous combustion: µO ,i 2

CnHmSpOqClkNy 98nCO + kHCl + pSO2 + xiNO + m-k y - x1 N2 + H2O (6) 2 2 The assumption is that C is oxidized to CO2, Cl to HCl, residual H goes to H2O, S to SO2, while atomic N in part produces N2 and in part NO. The splitting factor of the nitrogen fraction is xi ) ΨNO,iy and is evaluated through Bowman’s theory12 about heterogeneous combustion. More details are given in ref 10. The waste combustion is rate-determined by the diffusion of oxygen toward the flame front. The chemical reaction is almost instantaneous if compared to diffusion. Consequently, the combustion rate is

RW,i )

kx,ixO2,i

(7)

µO2,i

and the stoichiometric coefficient µO2,i can be obtained from eq 6:

µO2,i )

1 m-k -q n + 2p + xi + 2 2

(

)

(8)

The solid phase over each grate i is modeled by assuming the hypothesis of a granular bed and it is characterized by an equivalent particle diameter dp,i

3162


that allows calculating the mass-transfer coefficient: mix kx,i ) ShiDO ac 2,i i tot,i

where the mean waste thickness, sB,i, is deduced from chamber geometry:

(9)

Since the solid particles can be considered spherical, it follows that

6(1 - i) ai ) dp,i

(10)

sB,i )

MW,i + MI,i FWlGwG

(16)

By going back to the cascade of perfectly mixed reactors, which discretize the waste bed combustion, it is possible to write the following balance equations for the gas phase:

Sherwood number, Shi, is evaluated from the Colburn’s analogy

Shi ) jD,iReiSci1/3

{

with

Rei ) Sci )

(11)

FAvA,i µA a i µA

(12)

mix FADO 2,i

mix where jD,i and DO can be found in the literature.13 2,i The hypothesis of a granular bed, for the moving waste over the grates, looks adequate since the operator to the bucket, which discharges the waste into the loading hopper, has the task of homogenizing the rubbish in the pit. Equation 7 requires the evaluation of the exchange area (encapsulated in term kx,i) that can be brought back to the correlations for packed systems, provided that it is corrected by a proper coefficient. Such a coefficient takes into account the waste movement over the grates. Actually, waste movement is not integral with the grates but relative movements make it advance along the chamber. Therefore, surface renewal, achieved by discontinuous strokes, increases the exchange area and, as a consequence, the combustion rate. The correction coefficient Γi allows one to explicitly bind the exchange area, AE,i, to the number of strokes given in the time unit, NGS:

ai MW,i + MI,i AE,i ) Γi fW,i FW (1 - i)

(13)

Γi ) 1 + δ(1 - e-NGS,i/λ)

(14)

with fW,i being the waste mass fraction, over the ith grate: fW,i ) MW,i/(MW,i + MI,i). It is worth observing that Γi assumes an asymptotic value when the number of grate strokes increases. This is sound since, over a given NGS value, no further increase in the combustion rate is achieved. The two adaptive parameters, δ and λ, represent respectively the maximum correction reached in asymptotic conditions and the design value of NGS under regular operation. The waste void fraction, i, in eqs 10 and 13, can be evaluated from pressure drops averaged out from plant data measured over each grate. Consequently, the Ergun’s equation, for granular beds, can be adopted:

∆Pi 150µAvA,i(1 - i)2 1.75FAvA,i2(1 - i) ) + sB,i d 2 3 d 3 p,i

i

p,i i

(15)

where ΨNO,i is given elsewhere in the literature.10,12 2.1.2. Waste Movement. As aforementioned, one of the noteworthy aspects of grate kiln modeling is represented by the characterization of the waste movement. The mechanical system adopted in the primary combustion chamber comprises four reciprocating grates. Each grate can be operated separately by means of discontinuous strokes. A sudden and rather fast stroke moves the waste forward. Sometimes two or three consecutive strokes are given to make the waste advance. In our opinion, it is nonetheless better to operate single strokes in order to smoothen the combustion process while avoiding either stockpiles or voids. The kinematics of the grates has been schematized in a simplified way. In practice, the grate is represented as a moving shelf over which a uniform bed of waste is laid. Every stroke makes the waste advance a known segment, ∆xG. When the grate returns to the original position, the volume of waste belonging to that moving segment falls upon the following grate due to mass inertia (see Figure 3). The mass of waste that falls with a stroke is

MW,F,i ) FW∆xGwGsB,i

(18)

Given the number of strokes in the time unit, NGS,i, it is possible to evaluate the outlet mass flow rate from each section of the primary kiln: out ) MW,F,ifW,iNGS,i FW,i

(19)

out ) MW,F,i(1 - fW,i)NGS,i FI,i

(20)


Figure 4. Functional schematization of a grate stroke. Figure 3. Schematization of the waste movement over a reciprocating grate.

Doing so, we obtain a mean value of the waste flow rate over each grate. As a matter of fact, the discontinuity feature cannot be disregarded when modeling the combustion dynamics of the primary kiln. When the main process variables are observed, a peculiar oscillating trend can be pointed out. Such an aspect can only be explained by considering the discontinuous operation of the plant. When the grate is idle, the primary combustion air enters the corresponding perfectly mixed reactor while the hot gases are discharged from it. The inlet and outlet waste flow rates should be considered null as far as the grate is idle. In contrast, when a grate stroke is issued, the solid waste could be thought of falling instantaneously and numerically modeled by a Dirac delta, whose integral is given by eq 18. Actually, such an assumption would be both incorrect and impractical. The waste falls during the whole length of the stroke that can last a few seconds. A more correct formulation suggests the adoption of a bell-shaped function, such as the normal distribution, able to suitably describe the phenomenon and furthermore capable of numerically smoothing the discontinuity. Starting from the normal distribution,

1

f(z) )

x2πσ

(

exp -

)

(z - z0)2 2σ2

(21)

σ2 being the variance, z0 the central value, and +∞ f(z) dz ) 1, it is possible to adapt the distribution ∫-∞ to our modeling purposes. The independent variable, z, becomes the process time, t; the central value is t0,i + θi and MW,F,i is a multiplicative quantity. Eventually, we have out ) MW,F,ifW,i FW,i

1

x2πσ

out FI,i

(

exp -

2

) MW,F,i(1 - fW,i)

1

x2πσ2

)

correlated to the time length of the waste fall. Given a positive h value, the higher the variance, the lower the +h f(z) dz < 1. For small values of σ2 integral value of ∫zz00-h and high values of h, the following expression +h f(z) dz = 1 is quite correct. ∫zz00-h In our case, h is determined from the time interval between two consecutive grate strokes:

hi =

0.5 NGS,i

(24)

As far as σ2 is concerned, it must be determined so to satisfy the following equation:

∫0

3600

(

)

[t - (t0,k + θ)]2 1 MW,F,0 dt exp 2σ2 x2πσ = 3600 F h W t0,k < t < t0,k +

1 (25) NGS,0

where F h W is the desired mean mass flow rate of waste that should be fed to the primary kiln and t0,k is the initial time for the kth stroke of the feeding grate. The second adaptive parameter, θ, can be physically regarded as the time interval between the mechanical start of grate and the effective fall of waste, which is over it, on the following grate. Numerically, θ is chosen to avoid a sudden discontinuity between two consecutive values of waste flow rates due to previous and following strokes. Since the FW,0(t0) value is practically null before a new sudden stroke, it follows that θ is evaluated by the equation

θ ) kxσ2

(26)

2

[t - (t0,i + θi)] 2σ2

(

exp -

(22)

)

[t - (t0,i + θi)]2 2σ2

(23)

where t0,i is the instant when grate i starts moving (see also Figure 4). Specific attention should be devoted to parameters σ2 and θ. σ2, the variance of the curve, is responsible for the shape of the Gaussian distribution and is physically

where k is obtained through a variable substitution by setting

Z)-

t0 - (t0 + θ)

xσ2

Doing so, it is possible to determine the normal distribution that allows evaluating

θ)Z h xσ2 w k ) Z h )-

t0 - (t0 + θ)

xσ2

(27)

3164


2.1.3. Homogeneous Combustion. As discussed before, the gas phase, within the primary kiln, can be considered as a perfectly mixed system, given the high turbulence produced by inlet convective fluxes and the significant material and heat fluxes for the combustion of solid waste. Given the flow rates and compositions of gaseous streams produced by every heterogeneous section (physically bound to each reciprocating grate), it is possible to evaluate the inlet flow rate to homogeneous section that receives contributions of both secondary air and secondary leakage streams. Consequently, the inlet gas flow rates are determined by the following equations:

)

∑ i)1

dnHom,CO2 dt

out WG,j,i

dnHom,j out out in ) -(Wout Hom + WL )xHom,j + WHom,j dt j ) SO2,H2O,HCl dt

j ) CO,SO2,HCl,NO,H2O

out out ) -(Wout Hom + WL )xHom,CO2 + in RHom,CO2VHom + WHom,CO 2

dnHom,N2

NG

in WHom,j

The following reported here are the material balance equations that describe the homogeneous combustion within the furnace:

1 out out ) -(Wout Hom + WL )xHom,N2 - RHom,NOVHom + 2 in WHom,N 2

NG

in ) WHom,O 2

out WG,O ∑ ,i + (WA,S + WA,LS)xO ,A i)1 2

NG

WHom,N2 ) in xHom,j

)

(28)

out WG,N ∑ i + (WA,S + WA,LS)xN ,A i)1 2

in WHom,j in WHom,j

The reactions that take place in the homogeneous combustion section are

1 CO + O2 f CO2 2

(29)

N2 + O2 f 2NO

(30)

The kinetics of eq 29 is reported in ref 14 where a simplified expression is proposed:

dcCO 25000 ) 1.8 × 1010cCOxcO2 exp dt RT

(

)

(31)

Given the operating temperature (higher than 1000 °C) and the residence time (higher than a second) of the primary kiln, eq 31 predicts an almost complete conversion of CO to CO2. However, when we experimentally measure the CO content of the exhaust gas that leaves the furnace and enters the postcombustion chamber, we see non-zero values. Assuming that kinetics of eq 31 is correct, we have to introduce a fluid dynamics term to explain the relevant presence of CO. A simplistic but effective hypothesis suggests introducing a bypass stream that does not participate in homogeneous combustion but goes directly to postcombustion section. Such a bypass flow rate accounts for the presence of CO in the exhaust gas and allows describing quantitatively the imperfect mixing within the primary combustion chamber. As far as eq 30 is concerned, it is described by the Zeldovich kinetic mechanism15 that consists of the following reactions:

{

dt

)

-(Wout Hom

+

out Wout L )xHom,O2

O + N2 f N + NO N + O2 f O + NO N + OH f H + NO

(32)

(33) 1 - RHom,NOVHom 2

1 in R V + WHom,O 2 2 Hom,CO2 Hom

2

j ) O2,N2,CO,SO2,HCl,NO,H2O

NC

∑ j)1

dnHom,O2

2

dnHom,CO* out out ) -(Wout Hom + WL )xHom,CO* dt in (1 - fbypass) RHom,CO2VHom + WHom,CO dnHom,CO** out out ) -(Wout Hom + WL )xHom,CO** + dt in fbypass WHom,CO out xHom,CO )

nHom,CO* + nHom,CO** NC+1

nHom,j ∑ j)1

where CO* is the carbon monoxide fraction that reacts in the homogeneous phase, while CO** is the carbon monoxide portion that does not react since it belongs to the bypass stream whose fraction is equal to fbypass. The evaluation of mass flow rates is performed by the Bernoulli equation that, applied to the specific geometry of the furnace,16 is

v12 P1 - P2 ) γF1 2

(34)

where γ is a geometry coefficient that depends on the inlet section. Two separate cases should be considered when applying eq 34. In the first case an ambient flow rate enters the furnace while in the second, a combustion gas flow rate exits the primary kiln. If PPK < Pamb, we have

h x∆P Win L ) k

(35)

If PPK > Pamb, we have

ctot,g h Wout L ) k ctot,amb

x

Famb x∆P Fg

(36)

where k h can be determined if the air leakage flow rate


is known under steady-state conditions. Actually, when the incineration plant works properly, under design conditions, a negative relative pressure is kept within the process units to avoid the combustion gases being ejected from the seals. Consequently, during regular operation, eq 35 applies with respect to eq 36. Once Win L has been evaluated, we assume that it can be evenly divided into two contributions: the former refers to the air leakages under the reciprocating grates while the latter quantifies the secondary air leakages. Analogously, the gas flow rate that passes from the primary kiln to the postcombustion chamber is equal to

Wout Hom )

x

2PPK |PPK - PPC| out ‚ A γ MWgRT

(37)

Since the chamber pressure is almost equal to the ambient one and the furnace temperature is high, it is correct assuming the hot gas ideal and therefore applying the ideal gas equation for evaluating the volumetric quantities of eq 37. 2.1.4. Energy Balance. We chose to write one energy balance for the whole furnace. Such a decision was based on the assumption that the system is thermally homogeneous and solid and gas phases are in equilibrium. For sure, this assumption is quite strong since thermographic studies17,18 have shown qualitatively and quantitatively that the temperature profile of walls, gases, and waste bed are nothing but constant or in equilibrium. Consequently, a suitable description of temperature profiles would require the adoption of a detailed approach, like the zonal method,19 which takes into account the radiative heat transfer between the discretization sections of the furnace. Nonetheless, such a detail of description was regarded as too cumbersome for the numerical model that must be implemented on-line for monitoring, control, and supervision purposes. The energy balance over the whole furnace is then

dU ) H˙ in - H˙ out - Q˙ disp - Q˙ V,PK dt

it can be demonstrated that eq 38 becomes NC

(MScˆ p,S + M ˙ in S

∑ i)1

dT

nic˜ v,i)

NC

)-

dt

ˆ p,s dT in c

RT

d

NC

(39)

By applying the following formulas that are valid for ideal gases

(40)

and by assuming for the solid phase

cˆ v,S ≈ cˆ p,S w US ≈ HS

(42)

NC

in (WHom,i ∫T ∑ i)1

M ˙ in S

T

∫TT

in

in

c˜ p,i dT) ) (WA,PK + WA,L)

∫TT

amb

c˜ p,A dT (43)

cˆ p,S dT ) Fin ˆ p,S(T - Tin) + W(1 - ωH2O)c Fin WωH2O MWH2O

∫TT

eb

v c˜ p,H dT (44) 2O

From eq 44 it can be outlined that the solid heat capacity is assumed to be constant and the water contribution appears only as vapor sensible heat since the standard net heat of combustion accounts for liquid heating and phase change of the water fraction. • Reaction heat rate: For the homogeneous reactions we have NR

R ∆HRj (T)Rj(T) ) ∆HCOfCO RCO ∑ j)1 2

2

R + ∆HN R 2fNO NO

(45)

The solid combustion adopts a modified standard net heat of combustion, WLHC*, to account for partial combustion of the burning waste fraction that produces CO and partial conversion of the nitrogen waste fraction to NO:

(

RW,i

∑ i)11 - ω

- ω H 2O

WLHC +

)

ωC ωN R R ∆HCOfCO + ∆H (1 - φNO,i) (46) N2fNO 2 MWC MWN

niUg,i + MSUS ∑ i)1

Hg,i ) Ug,i + PV ˜ ) Ug,i + RT

NC

(

Each term in eq 42 can be applied and extended to the furnace context as follows: • Power for heating the reactants:

NG

NC

c˜ p ) c˜ v + R

c˜ p,i dT) -

∑ni) - Q˙ disp,PK - Q˙ V,PK dt i)1

I

WiHg,i + M ˙ S HS ∑ i)1

in

∆HRj (T)Rj(T) + RWWLHC + ∑ j)1

(38)

where

H˙ )

T

NR

∫TT

RWWLHC* )

U)

in (WHom,i ∫T ∑ i)1

(41)

• Power for increasing the total number of moles:

RT

d NC

∑ ni

dti)1

(47)

This term is present only for the gas phase. As shown in eq 41, the solid enthalpy and internal energy can be considered almost equal. Consequently, the contribution of eq 47 becomes null for the solid phase. • Dissipated power: The furnace dissipates energy by means of radiative and convective heat transfer. The model considers both contributions through the Hottel and Sarofim ap-

3166


proach.19 The net radiative flux exchanged between the hot gas and the furnace walls is

I Ti,W : insulator temperature of the ith discretization layer of the vertical wall

Q˙ rad ) σ(AGfWTg4 - AWfGTW4)

R : refractory temperature of the ith Ti,C discretization layer of the ceiling

(48)

Biardi and co-workers,20 starting from the pioneer work of Senkara,21 have shown that, for industrial furnaces, by adopting the one-gas-zone hypothesis, the following approximation holds:

AGfWTg = AWfGTW

(49)

Consequently, eq 48 becomes

TW 1 - k3 Q˙ rad ) σAGfW(Tg4 - TW4) k) 4 Tg 1-k

AGW 1 1 + -1 g R

(50)

(51)

According to the Hottel and Sarofim approach, the convective contribution can be expressed as a pseudoradiative flux: 4

4

The difference between vertical walls and ceiling is due to the different fluid dynamic regime. Actually, a forced convective regime prevails on the vertical walls while a natural convection characterizes the ceiling. Consequently, the only difference in the heat balance can be attributed to the laminar heat-transfer coefficient. Starting with the partial differential equation that describes the Fourier’s conduction problem through a one-dimensional plain wall,

∂2T ∂T ) kW 2 ∂t ∂z

hint

Q˙ conv ) σAGfW(Tg - TW )AGW 4σTGfW3

Tg + TW TGfW ) (52) 2

The convective heat coefficient, hint, is different for either vertical or horizontal walls. A forced convective regime applies to vertical walls while a natural regime characterizes the ceiling of the furnace (for further details see ref 22). • Absorbed power of the production of steam: The furnace walls, besides being isolated by a refractory and insulator thickness, have some boiling tubes that pass inside them to heat the process water flow rate and reduce thermal dissipations. Through design specifications, the power absorbed to preheat process water can be evaluated as a fraction, βPK, of energy dissipated by the kiln:

Q˙ V,PK ) βPKQ˙ disp,PK

(53)

2.1.5. Evaluation of the Refractory/Insulator Temperatures. A dedicated section is devoted to the evaluation of the refractory and insulator temperatures since thermal dissipation must be analyzed more deeply. Actually, the unknowns of the problem are as follows: int TS,W : skin temperature of the refractory wall in contact with the combustion gas int : TS,C

ext : skin temperature of the external casing wall TS,W in contact with the ambient air ext TS,C : skin temperature of the external casing ceiling in contact with the ambient air

where

AGfW )

I Ti,C : insulator temperature of the ith discretization layer of the ceiling

skin temperature of the refractory ceiling in contact with the combustion gas

R : refractory temperature of the ith Ti,W discretization layer of the vertical wall

(54)

it is possible to spatially discretize it by introducing a proper number of layers, NL: W Ficˆ p,i Vi

dTi int out ) q˘ in ˘ out i ) 1,...,NL (55) i Ai - q i Ai dt

int ) Ai+1 , as Since walls are flat, we have that Aout i shown in Figure 5. A further simplification can be done by evaluating the inlet and outlet heat fluxes from the steady-state Fourier equation, with boundary conditions:

T(zi) ) Ti

(56)

T(zi+1) ) Ti+1

Doing so, we obtain the inlet flux that leaves ith layer:

q˘ out ) -kW i

|

dT dz

i+1/2

) kW

Ti - Ti+1 zi+1 - zi

(57)

As a matter of fact, walls are not perfectly homogeneous. There is at least one discontinuity passing from refractory to insulator since thermal conductivity changes. If we assume that the interface between such materials is in steady-state conditions, it is possible to equal inlet and outlet fluxes, therefore obtaining

TRNRL - TI1 q˘ RI ) sI sR + 2kR 2kI

(58)

Another discontinuity source is located between combustion gases and internal refractory walls of the


Figure 5. Discretization of the refractory/insulator walls.

furnace. Once again, we assume a steady-state condition and by equating heat fluxes we obtain

hint(Tg,PK - Tint W) )

2kR int (T - TR1 ) sR W

(59)

As aforementioned, hint assumes different values if referred either to vertical walls or to the ceiling. Finally, the external skin temperature is again determined by applying the same approach to heat fluxes:

hext(Text S - Tamb) )

2kI I (T - Text S ) sI NIL

(60)

As described before, the heat-transfer coefficient can be evaluated adopting Hottel and Sarofim approach:

hext )

NukA + 4σmetTmet3 L

Tamb + Text S Tmet ) 2

(61)

2.2. Postcombustion Chamber. The postcombustion chamber of an incineration plant, by receiving the hot gas stream, produced in the furnace, has to carry out the following tasks: (1) keep the combustion temperature over a minimum value to assist the oxidation kinetics; (2) reach a suitable residence time to complete the oxidation of most resistant pollutant molecules; (3) produce a high gas turbulence to ensure a deep mixing between combustible and comburent fractions; (4) keep the oxygen content over a minimum value to increase combustion rate. At the inlet of the postcombustion chamber an auxiliary burner, fed by a methane stream, keeps the temperature of the chamber over a minimum value. Usually, the burner is switched off since the waste combustion is autothermal. Conversely, it is switched on during startup and shutdown operations in order to

avoid steep temperature changes. An inlet air flow rate keeps the outlet oxygen content over the limit prescribed by law. The chamber geometry comprises a rectangular section with two adjacent vertical ducts (upward and downward) connected by two consecutive elbows (see section 6 of Figure 1). Walls are made of two layers of refractory and insulating materials. Inside the refractory layer a series of boiling tubes absorbs a fraction of dissipated heat and produces a steam flow rate that feeds the cylindrical drum. Since turbulence and residence time of the postcombustion chamber are high and the process stream is gaseous, we chose to model both material and energy balances in terms of a perfectly mixed approach. Such a hypothesis allows one to highly simplify the numerical model by reducing the number of equations that describe the combustion process. 2.2.1. Material Balance. There is an analogy between equations wrote in section 2.1.3 to describe the material balances of the homogeneous gas phase within the primary kiln, and those describing the postcombustion chamber:

dnCO2 dt

out in in ) -Wout PC xCO2,PC + RCO2VPC + WPCxCO2,PC

dnSO2

(62)

out in in ) -Wout PC xSO2,PC + WPCxSO2,PC

(63)

out in in ) -Wout PC xH2O,PC + WPCxH2O,PC

(64)

dnHCl out in in ) -Wout PC xHCl,PC + WPCxHCl,PC dt

(65)

dt dnH2O dt

3168


dnN2

1 out in in ) -Wout PC xN2,PC - RNOVPC + WPCxN2,PC dt 2

dnO2 dt

(66)

1 1 out ) -Wout R V + PC xO2,PC - RNOVPC 2 2 CO2 PC in Win PCxO2,PC (67)

dnCO out in in ) -Wout PC xCO,PC - RCO2VPC + WPCxCO,PC dt

(68)

dnNO out in in ) -Wout PC xNO,PC + RNOVPC + WPCxNO,PC dt

(69)

As far as the pressure drop across the postcombustion chamber is concerned, the following equation evaluates it as a function of outlet flow rate:

PPC - PPK ) ∆Pdis + ∆Ploc

(70)

The distributed pressure drop is

∆Pdis )

2fDFgvg2LPC Dh,PC

(71)

where the drag coefficient, fD, is evaluated through Chu and Churchill formula:23

fD )

1 7 0.9 4 log 0.27 + Dh Re

[ (

( ) )]

2

(72)

Local pressure drops across elbows and in the inlet and outlet sections are

Fgvg2 ∆Ploc ) γ 2

(73)

where parameter γ depends on the geometry of the section. Once again, the equation of state for the hot gas is the ideal one since the operating temperature is high and the pressure is low. 2.2.2. Energy Balance. The energy balance over the postcombustion chamber is analogous to that described in paragraph 2.1.4 for the furnace. The overall energy balance

dU ) H˙ in - H˙ out - Q˙ disp,PC - Q˙ V,PC dt

(74)

can be made explicit in accordance with eq 42. The only exception is that there is one homogeneous phase; therefore, the solid contributions are missing: NC

(

dT

NC

in nic˜ v,i) ) -∑(Wi,PC ∫T ∑ dt i)1 i)1

T

NR

∑ j)1

∆HRj (T)Rj(T)

+ RT

in

c˜ p,i dT) NC

d

(

dt

ni) - Q˙ disp,PC - Q˙ V,PC ∑ i)1 (75)

The same approach proposed in section 2.1.4 applies also for the evaluation of terms Q˙ disp,PC and Q˙ V,PC. 2.3. Waste Heat Boiler. The postcombustion chamber is directly connected to the heat recovery section.

This section comprises a water tube boiler, a superheater, and an economizer (see sections 7-10 in Figure 1). Steam is superheated at a pressure of 10 bar and then expanded in a turbine that produces electrical energy. The exhaust steam is condensed and then recycled to the heat recovery units. Recovered enthalpy is used for cogeneration purposes such as hot water for civil use. 2.3.1. Boiling Section. As shown in Figure 1, the hot gases coming from the postcombustion chamber directly pass to the boiler that belongs to the same building as the primary and secondary kilns. The boiler can be schematized as a parallelepiped whose walls are filled by a raw of embedded tubes. Water absorbs a fraction of heat dissipated by the hot gas stream and starts vaporizing. A natural convection mechanism moves water upside, toward the cylindrical drum where actual vaporization takes place. The saturated steam, produced in the cylindrical drum, is then fed to superheater and eventually to turbine. The tube banks within the boiler section are licked by hot gases that from a mean temperature of about 1100 °C become quenched to 500 °C, which is the inlet temperature to superheater. Once again, the modeling of a waste heat boiler is performed by means of a simplified approach since a detailed model would be too cumbersome in terms of simulation times. In fact, a rigorous approach would require the description of (a) the biphasic motion of the liquid/vapor stream within boiling tubes and (b) the radiative heat transfer between gas and tubes. Subject (b) could be illustrated by the zonal method that would require the evaluation of view factors and direct exchange areas. Moreover, a detailed approach would entail the evaluation of energy stored in either liquid or vapor phases, considering the hydraulic head and both the material and thermal gradients across the liquid-vapor interface. In contrast, we adopted a simplified approach based on semiempirical relations24 that, nonetheless, allows obtaining good results in accordance with experimental data. 2.3.1.1. Cylindrical Drum. As mentioned before, we assumed that the boiling tubes are filled with water while the cylindrical drum is in part filled by a boiling water holdup and the remaining volume by saturated steam. Further hypotheses are that the inlet water flowrate is instantly mixed with water holdup in the vessel and that liquid and vapor phases are in equilibrium. Consequently, boiling water and steam holdups are at the same temperature and pressure. The material balance equations become

dML in - G˙ LfV ) G˙ H 2O dt

(76)

dMV ) G˙ LfV - G˙ V dt

(77)

Outlet steam flow rate, G˙ V, is adjusted by a valve installed on the vapor line that connects the cylindrical drum to superheater and finally to turbine. According to plant conditions, it looks correct assuming the operating pressure of turbine as constant. Actually, this assumption is partially correct since condensation pressure is the only measured value and it is downstream the turbine. That pressure value depends on the external operating conditions of condenser.


The distributed pressure drops along superheater tubes are

LSH ∆PSH ) 2fDFVvV2 Dh,SH

The hot gases material balance is trivial:

dng,B in out - Wg,B ) Wg,B dt

(78)

where the drag coefficient, fD, can be evaluated through conventional literature formula. Since

G˙ V ) FVAVvV

(79)

The fluid dynamic equation that evaluates pressure drops of the gas stream across the boiler has an analogous structure to eq 78 and therefore will be omitted. Conversely, the outlet gas flow rate from the boiler is

by combining eqs 78 and 79, we obtain

G˙ V ) AV

x

(80)

As far as the cylindrical drum pressure is concerned, it can be computed through the ideal gas equation by adopting VV as the effective volume occupied by steam inside the vessel. The pressure-temperature relationship is regulated by a thermodynamic equilibrium between liquid and vapor phases:25

(

P ) exp 73.649 -

7258.2 - 7.3037 log (T) + T

)

4.1653 × 10-6‚T2 (81) The liquid flow rate that is vaporized depends on the heat flux absorbed by tubes. Due to system inertia, which is bound to significant liquid and vapor holdups within the drum, it would not be correct assuming that heat absorbed by hot gases is instantaneously transformed into a steam flow rate that leaves the vessel. Actually, from a detailed point of view, vapor bubbles produced inside the boiling tubes move up toward the cylindrical drum. They go past the water head and finally leave the vessel. Although all these steps require a non-negligible amount of time, inertia can be disregarded with respect to the overall hot section dynamics. This assumption allows writing of the following equation:

G˙ LfV )

in in cˆ (Teq - TH ) G˙ Babs - G˙ H 2O p,L 2O

∆Hev(Teq)

(

1 - TR 1 - TR,ref

)

TR )

T (83) TC

in , depends on a Finally, the inlet water flow rate, G˙ H 2O dynamic equation that controls the liquid level inside the vessel. For the sake of simplicity, by adopting a conventional proportional-integral controller, we obtain

(

in in,0 ) G˙ H + kp c + G˙ H 2O 2O

1 τI

∫0t c dt

)

(84)

2.3.1.2. Boiler. With the term “boiler” we mean the vertical tube banks that are crossed by hot gases coming from the postcombustion chamber. A boiling water flow rate passes inside tubes while the hot gases are shellside.

(86)

mix MWg,B

Some words should be spent for the energy balance since it represents the main equation and quantifies heat exchange and therefore steam production. Due to high turbulence generated by hot gases that cross the tube banks, it is sound considering the whole boiler as a perfectly mixed volume. This means that gas temperature in the boiler is uniform and equal to the outlet value. Moreover, the high temperature of gas allows considering the heat transfer radiative. In addition, the ideal gas heat capacity is quite low. All these points lead to a pseudo-steady-state assumption for enthalpy balance:

˙ out ˙ V,B ) 0 H˙ in B - H B - Q

(87)

The Orrok and Hudson hypothesis24,26 assumes that only a fraction, m, of the inlet enthalpy to the boiler is absorbed by tube banks:

1

m)

155.3xq hg

1+ q)

H˙ in B Aexc B

hg )

∫T

Tin g ref

(88)

mix cˆ p,g dT

Consequently, the exchanged thermal power is

Q˙ V,B ) mH˙ in ) mW ˙ g,B

(82)

with

∆Hev(T) ) ∆Hev(Tref)

FgABvg

out ) Wg,B

FVDh,SH∆PSH 2fDLSH

(85)

∫TT

ing,B

ref

mix c˜ p,g dT

(89)

An alternative approach to the Orrok and Hudson formula is given by Konakov’s method27 that, although being based on an analogous energy balance, it uses different empirical parameters as shown in the following equation:

(

)

rad,out Tg,B 100

4

(

+ 239‚q‚

)

rad,out Tg,B -1 )0 hg Tref + mix cˆ p,g

(90)

rad,out Once the outlet gas temperature, Tg,B , has been determined by solving eq 90, it is trivial to evaluate the heat exchanged in the radiative section of the waste heat boiler. Unfortunately, in most cases the two approaches from Orrok-Hudson and Konakov do not agree and produce differing results. Annaratone24 proposed a synthesis of aforementioned methods that rad,out . Startalways produces intermediate values for Tg,B

3170


ing from Orrok and Hudson formulation, the exchanged fraction, m, of eq 88 is slightly modified:

Q˙ V,B ) mH˙ in )

1 in Wg,B 461.9 xq h0.15 g 1+ hg

∫TT

in g,B

ref

mix c˜ p,g dT

(91)

hrad ) σAGfT

It is then possible to determine the outlet gas temperature from the following equation: rad,out ) Tref + Tg,B

1000 mix 2.52 0.15 cˆ p,g hg + hg xq

in out ) Wg,B + WL,SH Wg,SH

∫TT

in g,SH

ref

mix out c˜ p,g dT ) Wg,B

∫TT

out g,B

ref

(93)

mix c˜ p,g dT +

WL,SH

∫TT

L,SH

ref

c˜ p,A dT (94)

Both heat exchangers are of the shell and tube typology. The hot gas is shell-side while either a steam or a water flow rate is tube-side. Since these exchangers differ for the physical properties of cold stream, only one model will be described in the following. The reader will have to adapt physical properties to such process streams. The heat-transfer coefficient for a hot stream is a function of gas velocity that depends on the geometry of tube banks. The tube layout can be either square/ rectangular or triangular. The shell-side mean gas velocity is

uinst if 2(sd - dt) g (st - dt) um ) st - dt um )

uinst 2(sd - dt)

(

)

T hH + T hC 2

3

(96)

Finally, the overall heat-transfer coefficient, U, for a shell and tube exchanger is

(92)

Actually, to account for thermal dispersion across the boiler walls, the inlet gas enthalpy, H˙ in B , is corrected by an adaptive coefficient, Rdisp,B, and the effective value used in eq 87. 2.3.2. Superheater and Economizer. An additional heat recovery section is installed downstream the boiler. A superheater raises the sensible enthalpy of saturated steam produced in the cylindrical drum. An economizer preheats the process water flow rate that subsequently flows in tubes inside the refractory walls of both the furnace and the postcombustion chamber. The hot gases coming from the boiler section are slightly quenched by a variable air leakage flow rate, WL,SH, due to the imperfect sealing of the connecting duct. The material and enthalpy balances that model the air leakage quenching are

in Wg,SH

of the exchanged heat is via a convective regime. Consequently, the radiative heat-transfer coefficient assumes a formulation analogous to the convective one (see also eq 51):

U)

1 1 1 + + ffou hint hconv + hrad

(97)

where ffou is the fouling factor. This quantity is usually not negligible: it depends mainly on the quality of waste fed to the furnace. For example, when animal flours, coming from the treatment of possible mad cows, are mixed and burnt with other types of wastes, to reduce their decidedly high heat of combustion, it takes only a few days for the tube banks to become encrusted. Whenever this happens, an operator has to manually remove the fouling layer over the tube banks by means of a steam jet directly withdrawn from plant production. This operation is quite expensive since it takes about 2 h and must be performed at low-temperature values, i.e., reduced steam production. Finally, the outlet exchanger temperatures are evaluated through the efficiency approach:22

{

1 - exp(-NTU(1 - r)) if r e 0.98 1 - r‚exp(-NTU(1 - r)) NTU η) if 0.98 < r e 1 1 + NTU

η)

r)

{

cMIN cMAX

cMIN ) Min{FHcˆ p,H,FCcˆ p,C} cMAX ) Max{FHcˆ p,H,FCcˆ p,C} NTU )

(98)

UAGT cMIN

If the steam flow rate is the stream having the minimum heat capacity and is tube-side, then in in in Tout C ) TC + η(TH,eff - TC )

(99)

in in in Tout H ) TH,eff - ηr(TH,eff - TC )

(95)

if 2(sd - dt) g (st - dt)

where sd ) xsl2+(st/2)2 for a square/rectangular tube pitch. We suggest referring to scientific literature22 for the correlations about heat-transfer coefficients of the shellside gases. The same approach should be devoted to the evaluation of heat-transfer coefficients for tube-side streams: either steam or water. In this case, the radiative heat term cannot be neglected, although most

otherwise in in in Tout C ) TC + ηr(TH,eff - TC )

(100)

in in in Tout H ) TH,eff - η(TH,eff - TC )

Once again, the effective inlet gas temperature, in TH,eff , is a function of heat dissipated by the heat exchanger. Therefore, a multiplicative coefficient (Rdisp,SH < 1 or Rdisp,EC < 1) should correct the inlet gas in enthalpy and as a consequence TH,eff .

Ind. Eng. Chem. Res., Vol. 44, No. 9, 2005 3171 Table 2. Enumeration of First Principle Equations Modeling the Hot Section Dynamics of the Incineration Plant material energy momentum balance balance balance total

section primary kiln postcombustion chamber boiler section overall hot section

55 10 6 71

29 15 9 53

2 1 5 8

86 26 20 132

Pressure drops across the heat exchanger depend on the number of longitudinal tubes, nT, present in the bank: 2

∆Pexc ) fDFgvg,MAX nT

(101)

Green and Perry28 report several correlations for the drag coefficient, fD, as a function of the tube pitch geometry. 3. Model Validation The whole hot section of the incineration plant, which has been modeled in section 2, is numerically summarized in Table 2 that lists the number of equations of the primary kiln, postcombustion chamber, and boiler. One hundred and thirty-two differential and algebraic equations (DAE system) describe the hotsection dynamics of the plant. The numerical system is intrinsically stiff due to coexistence of both differential and algebraic equations. The numerical solver adopted is DASPK29 that with version 3.0 has gained some advantages over DASSL30 mainly in terms of robustness and features. However, it should be underlined that the preconditioned Krylov feature (PK) was not exploited due to the rather small dimensions of the system (only 132 equations). For the sake of clarity, Krylov method and preconditioning are methodologies adopted when solving linear systems, Ax ) b, of large dimensions through iterative methods. If large linear systems with thousands of equations are involved, the classical solution based on factoring the coefficient matrix, A, becomes impractical in terms of both CPU time and memory allocation. In this case, an iterative procedure based on the Generalized Minimum Residual method (Krylov method) is viable and effective.31 The overall memory allocation required by matrix A can be further decreased by pre-multiplying it (and also vector b) by a suitable preconditioning matrix, P. Then the Krylov method is applied to this altered, but equivalent, linear system. The extra iterations resulting from this approximation are offset by the lower storage and linear system solution costs for a reduced PA matrix. Incidentally, a linear system solution is required whenever the order of the DAE solver, the step length, or the Jacobian matrix change. The solver was able to integrate the DAE system in almost all process transients. The intrinsically oscillating behavior of the furnace avoids the system reaching any steady-state condition. Consequently, the CPU time taken for a conventional plant simulation is significant. About 2 min of CPU time is required to simulate a real plant dynamics of 30-40 min, adopting respectively as absolute and relative tolerances the following values: 1 × 10-10 and 1 × 10-8. DASPK, as any other solver, is unable to increase significantly the integration step since process variables have an unsteady and bumpy behavior (see also Figures 6 and 7).

Table 3. Pseudo-stationary Conditions Describing the Mean Operation of Incineration Plant Input Variables number of strokes of the -1 feeding grate [h ] inlet waste flow rate [kg/h] number of strokes of the first reciprocating grate [h-1] number of strokes of the second reciprocating grate [h-1] number of strokes of the third reciprocating grate [h-1] number of strokes of the fourth reciprocating grate [h-1] primary air flow rate under the grates [Nm3/h] secondary air flow rate to the furnace [Nm3/h] Output Variables outlet gas temperature from the furnace [°C] outlet gas temperature from the postcombustion chamber [°C] O2 molar fraction in the outlet stream from the postcombustion chamber [%] CO content in the outlet stream from the postcombustion chamber [mg/Nm3] steam flow rate [t/h]

40 4000 23 20 13 23 12500 6500

1035 1115 8.5 11 12

3.1. Identification of the Adaptive Parameters. As discussed in section 2 of this paper, some adaptive parameters were introduced in the first principle model to account for the uncertainties of the process. The main adaptive parameters are: δ and λ in eq 14; fbypass in eq 33; βPK in eq 53; βPC for term Q˙ disp,PC in eq 75; ffou in eq 97. Their evaluation is performed once for all by considering the mean operating condition of the incineration plant. Table 3 reports the values of main input and output process variables referring to a pseudostationary condition. Once determined by means of a nonlinear regression procedure, the adaptive parameters are assumed constant throughout the dynamic simulation of the plant.11 The reader could argue that these adaptive parameters may not be constant if the plant operates at different regimes. However, it should be observed that the dynamic simulation lasts a few days. During this short time interval, the variability of adaptive parameters is well below that of unmeasurable process variables whose evaluation will be described in the following. To be more precise, ffou in eq 97 has the highest degree of variability within the set of adaptive parameters. However, as described in paragraph 2.3.2, the periodic removal of fouling by a steam jet procedure keeps its variability moderate. The other adaptive parameters depend mainly on the geometry, fluid dynamics, and mechanical movement of the reciprocating grates. Consequently, they can be regarded as constant. Another quite important aspect is represented by waste characterization. The aforementioned parameter identification is performed assuming a mean ultimate analysis of waste fed to the furnace. Table 4 reports such a characterization. As will be discussed in the following, the extreme variability of waste, deposited in the pit, is one of the main aspects responsible for the unsteady behavior of the incineration plant. The waste is heterogeneous since it comes from municipal, hospital, agricultural, and zootechnical collections. Moreover, the same waste typology is a function of wet mass fraction and apparent density that change with atmospheric

3172


Figure 6. Pendulum-like trend of the steam flow rate. Comparison between experimental data (a) and simulated values (b).

Figure 7. The adoption of increasing values for the moving average of process variables allows smoothing of the local noisy oscillations. At the same time it introduces a delay time into the process dynamics. Table 4. Waste Characterization (Ultimate Analysis) Adopted for the Identification of Adaptive Parameters H mass fraction C mass fraction O mass fraction N mass fraction S mass fraction Cl mass fraction inert mass fraction wet mass fraction standard net heat of combustion [kJ/kg]

0.0166 0.2622 0.0754 0.0010 0.0008 0.0065 0.2675 0.3700 12500

conditions (such as rain and sun) and permanence within the pit. The model validation is performed by comparing predicted numerical values with experimental ones. The plant control structure comprises a Distributed Control

System (DCS) that collects on-line the measured values and stores them in a proprietary database. An external computer may be connected to the DCS to acquire in real time the measured data. The input/output exchange protocol is based on the OPC (OLE for Process Control) layer (www.opcfoundation.org) that has the significant feature of abstracting from both hardware and software implemented in the DCS structure. Also working with a high level programming language, it is quite simple to exchange input/output data by connecting the acquisition client to the OPC server running on the DCS. The acquisition computer and DCS are simply connected through an ethernet cable and exchange data by means of TCP/IP protocol. The validation procedure required acquiring about 40 process variables. We chose a sampling time of 2 s to have a high degree of precision for most sensitive variables that often exhibit fast dynamics. Once some days of on-line data were acquired, we were able to study off-line the plant behavior and to mimic its evolution through the numerical simulator. The first benchmark was the comparison of simulated data with the oscillating behavior of main process variables, produced by discontinuous feeding conditions. Figure 6 shows a comparison between steam flow rate, measured on-line along a 1 h interval, and corresponding simulated curve. The pendulum-like motion of the plant is well-reproduced in terms of both oscillation amplitudes and almost periodic intervals. Actually, Figure 6 reports a pseudo-stationary evolution of the incineration process. Steam flow rate oscillates around an almost constant mean value and the curve period recalls the grate strokes that push waste into the furnace. It should be underlined that the 1 h interval refers to a manual


Figure 8. Diagrams describing the variability of waste heat of combustion, over an 8 h interval, inferred through the HCl measurement after the boiler section and the reconciled total inlet air flow rate.

Figure 9. Measured and inferred input data for the validation procedure.

operation condition. The operator manually set the feeding and reciprocating grate frequencies and kept them constant throughout the 1 h interval. Before the introduction of an on-line Model Predictive Control

(MPC) structure (which was made possible by the dynamic numerical model described in this paper) it was common use to operate the feeding grate by means of three consecutive and sudden strokes. Doing so, the

3174


Figure 10. Validation procedure: comparison between experimental data and simulated output values.

operator was sure that a significant amount of waste fell into the furnace and, with the total number of strokes in the time unit (1 h) almost constant, he could reduce his monitoring activity to one-third of that necessary if only one single stroke had been issued at a time. Three consecutive strokes were responsible for both high amplitudes and high pendulum-like periods. By keeping the overall frequency and therefore waste flow rate constant, and, at the same time, by equally spacing each single stroke, it was possible to eventually reduce both amplitude and period of the steam flow rate curve. But this is another story because it is connected to the implementation of an advanced control strategy (MPC) that will be described in a future paper. 3.2. The Moving Average Problem. A secondary problem that we had to face was the high-frequency oscillations produced by discontinuity sources of the process (mainly the reciprocating grates). Since each grate has an independent pulsing frequency, process variables have a noisy local behavior that is superimposed to the distributed dynamics produced by external disturbances. The local high-frequency oscilla-

tions, reported in Figure 7, can be visually reduced in amplitude by evaluating the moving average of the corresponding process variable. Doing so, the overall dynamic transient can be easily detected but a time delay is introduced in the system. As shown in Figure 7, high values of the characteristic moving average time, tma (defined as a multiple of the sampling time), reduce local oscillations but increase the induced sluggishness of the system. Conversely, low tma values increase the simulated process noise that gets summed to on-line noise produced by external flickering of disturbances and of measuring devices. Both asymptotic conditions should be avoided for control purposes. Consequently, values of tma equal to either 1 or 10 min should be avoided, while tma ) 5 min looks to be a good compromise as shown in Figure 7. 3.3. Measurable and Unmeasurable Input Process Variables. As mentioned before, the validation procedure is based on comparing on-line data with corresponding simulated values. The numerical model requires a number of input data to simulate the dynamics of the incineration process. Some input data, such


as primary and secondary air flow rates to the furnace and the postcombustion chamber, are easily and continuously measured on-line and stored in the DCS database. Other input variables cannot be measured because either the measuring device is missing (i.e., CO fraction in the outlet gas from furnace) or the quantity is unmeasurable on-line (i.e., waste heat of combustion). As mentioned before, waste is only partially known. The operator usually sets the feeding grate frequency but it is not possible to measure the inlet waste flow rate since the orange-peel bucket is not equipped with a dynamometer. To correctly identify waste characteristics, we based our approach on two distinct assumptions. On one side, we adopted a steady-state reconciliation procedure to infer the waste inlet flow rate to the furnace. On the other side, the waste heat of combustion was evaluated dynamically by stating that the municipal fraction has an almost constant quality while fluctuations are mainly due to the hospital fraction that is mixed in the pit.10,11 An indirect measure of hospital waste content is given by the HCl fraction measured before the separation units. Actually, the plastics content of hospital waste is higher than the municipal one (also due to the differentiated collection performed by municipality). Figure 8a shows that both trends of HCl fraction (after the hot section) and waste heat of combustion are comparable when put on the same time scale. The coaptation procedure, produced by a nonclassical approach to steady-state reconciliation,3,11 if dynamically iterated, produces good results as far as the total amount of inlet air flow rate is taken into account (see Figure 8b). Actually, the total inlet air flow rate, by summing up leakage contributions, is responsible for diluting or concentrating HCl in the hot gas flow rate, regardless of waste heat of combustion. This point explains why some portions of the two diagrams of Figure 8a are not synchronously overlapping. A sequence of equally spaced steady-state reconciliation procedures, run every 5 min over an 8 h interval, identified the unmeasurable input data that, coupled through over-sampling to those acquired on-line every 2 s from the DCS, represented the overall validation input data set. Figure 9 reports some measured or inferred input data adopted for validating the numerical model of the incineration hot section. The plant dynamics was integrated by the DAE solver with 2 s output times, tout (the same as the sampling time). Once the tout value was reached, a new set of input data was assigned as the initial condition and integration was iterated for two more seconds. The overall 8 h interval allows taking into account the typical evolution and disturbances of the incineration plant. Figure 10 shows good agreement between measured and simulated output values. Not only the main process variable, i.e., steam flow rate, but also furnace and postcombustion outlet temperatures and also oxygen molar fraction after the postcombustion chamber have comparable trends. The accordance can be observed in terms of both absolute values and response times. As far as steam flow rate is concerned, it is possible to observe a maximum difference of half a ton between experimental and simulated data. Actually, such a difference is quite low and process operators consider it negligible. In our opinion, this is a good result that represents also a necessary condition for implementing a model-based control of the incineration plant.

4. Conclusions The good agreement with experimental data, shown by the first principle numerical model, described in this paper, allows consideration of it as a sound candidate for model-based control. The whole hot section of an incineration plant, comprising a reciprocating grate kiln, a postcombustion chamber, and a boiler, was modeled in detail. The resulting differential-algebraic system, once solved, describes the dynamics of those units. Such a model becomes the engine for predicting the plant response necessary for a model predictive control algorithm to optimize steam production. The optimizing routine also has to satisfy both law and process constraints that are correctly described by the numerical model. The validation procedure, described in previous paragraphs, was proven reliable for this task. The unmeasurable input variables, such as waste characteristics, which are required to run the dynamic simulation of the hot section, can be inferred through an online steady-state reconciliation procedure that is run at equally spaced time intervals. Through that procedure, unknown quantities such as waste heat of combustion, ultimate composition, and air leakages can be determined and used as input values together with on-line measured input variables. Acknowledgment The authors want to acknowledge the precious work done by Luigi Borsini and Pietro Manzoni in their master degree thesis. A special thank you goes to the technical staff of AGAC incineration plant in Reggio Emilia, Italy. Notation A ) area, m2 a ) specific surface per unit volume, m-1 AE ) waste exchange area, m2 AGfT ) gas-tube direct exchange area, m2 AGfW ) gas-wall direct exchange area, m2 AWfG ) wall-gas direct exchange area, m2 c ) gas molar concentration, kmol/m3 cp ) heat capacity at constant pressure, kJ/(kmol K) or kJ/ (kg K) cv ) heat capacity at constant volume, kJ/(kmol K) or kJ/ (kg K) d ) equivalent diameter of waste particles, m Dh ) hydraulic diameter, m DO2 ) oxygen diffusivity, m2/s dt ) tube diameter, m F ) mass flow rate, kg/s f ) mass fraction fbypass ) bypass fraction fD ) drag coefficient ffou ) fouling factor, (m2 K)/W G˙ ) water or steam flow rate, kg/s H ) enthalpy, kJ/kmol or kJ/kg h ) heat-transfer coefficient, W/(m2 K) h ) time interval between two consecutive grate strokes, s hg ) specific enthalpy, kJ/kg i ) grate or reactor or zone or layer counter j ) chemical species or reaction counter jD ) Colburn factor k ) thermal conductivity, W/(m K) kp ) proportional gain of the controller kx ) mass-transfer coefficient, kmol/(m2 s) L ) tube length, m

3176


l ) length, m M ) mass holdup, kg m ) enthalpy fraction MW ) molecular weight, kg/kmol n ) number of moles, kmol NC ) number of components NG ) number of reciprocating grates NGS ) number of grate strokes in the time unit, h-1 NIL ) number of insulator layers NL ) number of layers NR ) number of reactions NRL ) number of refractory layers nT ) number of tubes NTU ) number of transfer units Nu ) Nusselt number P ) pressure, Pa Q ) heat, kJ q ) heat flux, kW/m2 R ) gas-law constant, J/(mol K) Re ) Reynolds number RHom,X ) specific reaction rate for gas species X, kmol/(m3 s) Rj ) specific reaction rate for reaction jth, kmol/(m3 s) RW ) waste reaction rate, kg/s s ) thickness, m Sc ) Schmidt number sd ) square/rectangular tube pitch, m Sh ) Sherwood number sl ) longitudinal pitch, m st ) transversal pitch, m T ) temperature, K t ) time, s t0 ) initial time of the grate stroke, s U ) internal energy, kJ/kmol (∼) or kJ/kg (/\) U ) overall heat-transfer coefficient, W/(m2 K) uin ) inlet gas velocity before the tube banks, m/s um ) mean gas velocity through the tube banks, m/s V ) volume, m3 v ) gas velocity, m/s W ) molar flow rate, kmol/s w ) width, m WLHC ) waste heat of combustion, kJ/kg x ) molar fraction z ) axial coordinate Greek Letters R ) adaptive coefficient β ) adaptive parameter for the absorbed heat by process water Γ ) correction coefficient γ ) geometry coefficient for the evaluation of the pressure drops ∆Hev ) heat of vaporization, kJ/kg ∆HR ) heat of reaction, kJ/kmol ∆P ) pressure drop, Pa ∆xG ) axial movement of the reciprocating grate, m δ ) adaptive parameter for the waste exchange area ) emissivity ) waste void fraction ) wall roughness c ) distance between the controlled variable and its setpoint η ) heat-exchanger efficiency θ ) central time value of the waste normal distribution, s λ ) adaptive parameter for the waste exchange area µ ) gas viscosity, kg/(m s) µO2 ) stoichiometric molar ratio between oxygen and waste φNO ) partial conversion coefficient for the waste N f NO reaction F ) mass density, kg/m3 σ ) Boltzmann constant, W/(m2 K4)

σ2 ) variance τi ) integral gain of the controller ψNO ) split factor for nitrogen of Bowman’s theory ω ) mass fraction Subscripts, Superscripts, and Accents A ) air abs ) absorbed amb ) ambient B ) boiler B ) waste bed C ) ceiling C ) cold gas C ) critical conv ) convective dis ) distributed disp ) dissipated EC ) economizer eq ) equilibrium exc ) heat exchanger ext ) external F ) falling G ) grate g ) gas GT ) gas-tube GW ) gas-wall G f W ) gas to wall I ) inert I ) insulator H ) hot gas Hom ) homogeneous H2O ) liquid water (for the waste heat boiler section) in ) inlet or inner int ) internal L ) leakage L ) liquid loc ) local LS ) secondary leakage L f V ) vaporizing flow rate (liquid to vapor) M ) moisture ma ) moving average MAX ) maximum met ) metal MIN ) minimum mix ) mixture out ) outlet or outer PC ) postcombustion chamber PK ) primary kiln, furnace R ) reaction R ) reduced R ) refractory rad ) radiative ref ) reference condition S ) secondary S ) skin SH ) superheater tot ) total V ) steam v ) vapor W ) wall W ) waste 0 ) inlet to the first reciprocating grate, initial A h ) mean value of quantity A A˙ ) flux value for quantity A

Ind. Eng. Chem. Res., Vol. 44, No. 9, 2005 3177 A ˜ ) molar value for quantity A A ˆ ) mass value for quantity A

Literature Cited (1) Ricci, M.; Tornavacca, A.; Francia, C. Gestione Integrata dei Rifiuti Urbani: Analisi Comparata dei Sistemi di Raccolta; FederAmbiente: Monza, 2003. (2) Manca, D.; Rovaglio, M.; Pazzaglia, G.; Serafini, G. Inverse Response Compensation and Control Optimization of Incineration Plants with Energy Production. Comput. Chem. Eng. 1998, 22 (12), 1879-1896. (3) Rovaglio, M.; Manca, D.; Rusconi, F. Supervisory Control of a Selective Catalytic Reactor for NOx Removal in Incineration Plants. J. Waste Manage. 1998, 18, 525-538. (4) Ibanez, R.; Andres, A.; Viguri, J. R. Characterisation and Management of Incinerator Wastes. J. Hazard Mater. 2000, 79 (3), 215-227. (5) Zabaniotou, A.; Giannoulidis, N. Incineration of Municipal Solid Waste with Electricity Production and Environmental Safety: The Case of a Small Capacity Unit in Greece. Energy Sources 2002, 24 (2), 115-126. (6) Ljung, L. System Identification: Theory for the User; Prentice Hall: Upper Saddle River, NJ, 1999. (7) Manca, D.; Bonanomi, E.; Ciapparelli, F.; Rovaglio, M. Model Predictive Control by Neural Networks Identification. AIDIC Conference Series, ERIS, ISBN 0391-7401, Pierucci, S., Ed.; Milano, Italy, 1999; pp 201-210. (8) Park, S. I.; Kyong, N. H.; Park, Y. J.; Lee, S. K., NumericalSimulation to Control Rotary-Kiln Incineration of Municipal SolidWaste. Energy 1994, 19 (2), 179-186. (9) Cook, C. A.; Cundy, V. A.; Larsen, F. L.; Lighty, J. S. Comprehensive heat transfer model for rotary desorbers. Can. J Chem. Eng. 1996, 74 (1), 63-76. (10) Rovaglio, M.; Manca, D.; Biardi, G. Dynamic Modeling of Waste Incineration Plants with Rotary Kilns: Comparisons between Experimental and Simulation Data. Chem. Eng. Sci. 1998, 53 (15), 2727-2742. (11) Rovaglio, M.; Manca, D. Reconciliation and Model Identification as a First Step for On-line Optimization of Incineration Processes. Proceedings of International Conference on Incineration & Thermal Treatment Technologies, IT3, Cohen, L., Ed.; Irvine, CA, 1997; pp 587-596. (12) Bowman, C. T. Kinetics of Pollutant Formation and Destruction in Combustion. Prog. Energy Combust. Sci. 1975, 1, 33. (13) Bird, R. B.; Steward, W. E.; Lightfoot, E. N. Transport Phenomena; John Wiley: New York, 2001. (14) Niessen, W. R. Combustion and Incineration Processes: Application in Environmental Engineering; Marcel Dekker Inc.: New York, 2002. (15) Zeldovich, Y. B. Acta Physicochim. URSS 1946, 21, 557.

(16) Douglas, J. F.; Gasiorek, J. M.; Swaffield, J. A. Fluid Mechanics; Pitman Press: Bath, 1980. (17) Manca, D.; Rusconi, F.; Arrigoni, M.; Cattaneo, F.; Rovaglio, M. Temperature Mapping of a Combustion Chamber by Infrared Thermographic Image Recognition. Proceedings of the Combustion Meeting, XXII Event, Salatino, P., Ed.; Napoli, Italy, 2-5 May, 1999; pp 175-178. (18) Manca, D.; Rovaglio, M. Infrared Thermographic Image Processing for the Operation and Control of Heterogeneous Combustion Chambers. Combust. Flame 2002, 130 (4), 277-297. (19) Hottel, H. C.; Sarofim, A. F. Radiative Transfer; McGrawHill: New York, 1967. (20) Biardi, G.; Pellegrini, L.; Grottoli, M. G.; Rovaglio, M. Chemical Furnaces under Transient Conditions; Proceedings of CEF’87, Ferraris, G. B., Ed.; Taormina: Milano, Italy, 1987; pp 559-564. (21) Senkara, T. Wa¨ rmetechnische Rechnungen fu¨ r gas-und o¨ lbehizte Wa¨ rmeo¨ fen; Vulkan-Verlag: Essen, 1977. (22) Incropera, F. P.; De Witt, D. P. Fundamentals of Heat and Mass Transfer; John Wiley: New York, 2001. (23) Chu, H. N. S.; Churchill, S. W. The Development and Testing of a Numerical Method for Computation of Laminar Natural Convection in Enclosures. Comput. Chem. Eng. 1977, 1, 103-108. (24) Annaratone, D. A new formula for the calculation of the exhaust gases temperature at the exit of the combustion chamber. Termotecnica 1984, 12, 57-61. (25) Daubert, T. E.; Danner, R. P. D.I.P.P.R., a Data Bank Compilation Table of Properties of Pure Compounds; DIPPRPress: New York, 1985. (26) Gumz, W. Die Vorga¨nge in den Feuerungen bei hohen Temperaturen. Mitt. Wgb H. 1955, 38. (27) Konakov, P. K. Waermeuebertragung in Kesselfeuerungen; Acad. Wiss. USSR, 1952. (28) Green, D. W.; Perry, R. H.; Maloney, J. O. Perry’s Chemical Engineers’ Handbook; McGraw-Hill: New York, 1997. (29) Brown, P. N.; Hindmarsh, A. C.; Petzold, L. R. A Description of DASPK: A Solver for Large-Scale Differential-Algebraic Systems; Lawrence Livermore National Report UCRL; Lawrence Livermore National Laboratory: Livermore, CA, 1992. (30) Petzold, L. R. A Description of DASSL: a Differential/ Algebraic System Solver, Scientific Computing; Stepleman, R. S., et al., Eds.; North-Holland: Amsterdam, 1983; pp 65-68. (31) Brown, P. N.; Hindmarsh, A. C.; Petzold, L. R. Using Krylov Methods in the Solution of Large-Scale DifferentialAlgebraic Systems. SIAM J. Sci. Comput. 1994, 15, 1467-1488.

Received for review May 26, 2004 Revised manuscript received February 17, 2005 Accepted February 18, 2005 IE0495473

Numerical Modeling of a Discontinuous ... - ACS Publications

Recommend Documents