Distributed Parameter Model of Black Liquor Falling-Film Evaporators

In the first of two parts of this work, a distributed parameter model, based on first-principles knowledge about the fluid dynamics and heat-transfer ...
1 downloads 0 Views 627KB Size
Ind. Eng. Chem. Res. 2004, 43, 8117-8132

8117

Distributed Parameter Model of Black Liquor Falling-Film Evaporators. 2. Modeling of a Multiple-Effect Evaporator Plant Z. I. Stefanov† and K. A. Hoo* Department of Chemical Engineering, Texas Tech University, Lubbock, Texas 79409-3121

In the first of two parts of this work, a distributed parameter model, based on first-principles knowledge about the fluid dynamics and heat-transfer processes, was developed for the evaporation side of a single plate (lamella) of a falling-film evaporator. In the second part, the work expands on the single-plate model to develop a fundamental model of a falling-film evaporator that accounts for the condensation side of the plate, the heating/flashing section at the evaporator entrance, the evaporator inventory, mixing, and recirculation flows. On the basis of the single evaporator model, a multiple evaporator plant that consists of one superconcentrator (three falling-film evaporators as one unit) and four falling-film evaporators is modeled. The model, its numerical solution, and its behavior at different operating conditions are presented and discussed. 1. Introduction The first part of this work presented a first-principlesbased model of the evaporation side of one plate (lamella) in a falling-film type of evaporator. Two simplifying assumptions were made in the development: uniform heating along the plate wall and the absence of a heating section at the beginning of the plate. In reality, there is generally nonuniform heating along the plate wall, and there may be a nonnegligible heating section. In this work, these assumptions are removed. The plate wall is heated by saturated steam whose condensation is not uniform because the condensate film parameters are functions of the plate length. Additionally, the corresponding heat-transfer coefficient changes along the length of the plate. A heating section is introduced at the beginning of the plate if the temperature of the fluid feed to the plate is below the boiling point that corresponds to the current saturation pressure. If the temperature of the fluid feed to the plate is above the boiling point, flashing rather than heating occurs. Figure 1 is a schematic of a falling-film evaporator. On the basis of experimental investigations by Patterson et al.,1 the wall dynamics are assumed to be very fast (the wall thickness is ∼1.5 mm). When the liquor leaves the plate stack, it enters the evaporator inventory (see Figure 1), which is controlled. Thus, the evaporator inventory can be considered a buffer and a sink. The evaporator inventory is modeled as a well-mixed thermal mixing tank. The recirculation pump mixes the feed with the contents of the evaporator inventory. This pump is modeled as a steady-state thermal mixer with the assumptions that the mixing dynamics are extremely fast and that the pump is an ideal mixer. The pump discharge is split into two flow streams; one flow stream is returned to the plate stack, and the other is the product stream. Two other assumptions made are as follows: heat loss to the surroundings is negligible, and the saturation * To whom correspondence should be addressed. E-mail: [email protected]. † Currently affiliated with Tembec Inc., St. Francisville, LA.

Figure 1. Schematic representation of a falling-film evaporator.

pressures at the condensation and evaporation sides of the plate are time-invariant. When multiple evaporators (MEVs) are considered, the pressure dynamics will be explored more closely. The paper is organized as follows. Section 2 describes the development of the model for one falling-film evaporator detailing the transport phenomena that occur (heating, flashing, and condensation) as well as accounting for the evaporator inventory and an energy balance at the wall. Section 3 develops a fundamental model for an evaporator plant that contains four fallingfilm plate-type evaporators and one superconcentrator. The model development pays particular attention to the pressure dynamics and unit-unit interdependencies.

10.1021/ie049611g CCC: $27.50 © 2004 American Chemical Society Published on Web 11/11/2004

8118

Ind. Eng. Chem. Res., Vol. 43, No. 25, 2004

Section 4 provides a discussion of the results of the MEV plant and the validation of the steady-state results with real plant data. The final section, section 5, summarizes the work and provides some concluding remarks. 2. Model Development of One Falling-Film Plate-Type Evaporator 2.1. Hydrodynamics. With any change in the fluid properties and operating conditions, the falling film will transition from a turbulent to a wavy-laminar and eventually to a laminar regime. Determination of the transition point has been studied by several researchers.2-4 However, the data are for a limited number of liquids with properties that are very different from those of the one of interest in this work. A more quantitative approach such as that provided by Chun and Seban5 may be considered. However, it was shown in ref 6 that the predictions obtained from these correlations were not correct for the fluid of interest. In light of this, the turbulent model of Alhusseini and Chen,7 together with an empirical correlation for the wavy-laminar regime,8 is used. The proposed expression, eq 1, uses the calculated heat-transfer coefficients for both turbulent and wavy-laminar regimes, thus providing a smooth transition

hE,total ) (hE,l5 + hE,t5)1/5

(1)

where hE,l is the laminar heat-transfer coefficient of evaporation and hE,t is the turbulent heat-transfer coefficient of evaporation. This equation is used in this work to calculate the heat-transfer coefficients for evaporation, condensation, and sensible heating over the entire range of operating conditions. 2.2. Evaporation. The heat transfer of evaporating or condensing films, in the laminar region, is a wellstudied subject.4,5,9 A number of accurate numerical solutions to the fundamental momentum, mass, and energy balances for wavy films were developed in refs 10 and 11. A general characteristic of these studies is the occurrence of large irregular (e.g., varying amplitude, asymmetric) waves on the film surface. Experimental research in refs 12 and 13 confirmed this observation. The reason for the increase in the heat transfer is attributed to the transport of a large part of the fluid in the observed solitary waves. Thus, the average film thickness decreases, hence the increase in the value of the heat-transfer coefficient. The following correlation was obtained experimentally by Alhusseini et al.8,14

hE,wl ) 2.65Re-0.158Ka0.0563

(2)

where Re is the film Reynolds number and Ka is the Kapitza number and is valid in the Prandtl number range of 1.7-47. According to Miyara,15 the influence of the Prandtl number is significant in the range of 0.110. For higher values, the influence is much smaller. This equation will be used in this study to obtain the heat-transfer coefficient of evaporation and condensation in the laminar region. 2.3. Sensible Heating. The information on the sensible heating of wavy-laminar falling films is considerable less than the information on condensation and

Figure 2. Experimental design.

evaporation. Heating without evaporation of wavy films is not well studied either experimentally or theoretically. One of the few reliable sources is the work of Kapitza.2 Kapitza determined experimentally and theoretically that the heat-transfer coefficient of sensible heating is 1.21 times greater than that for pure laminar conditions. If this information is to be used, it is necessary to calculate the heat-transfer coefficient of sensible heating for the pure laminar regime. This problem already has been studied extensively by Pohlhausen.16 The analytical solution is based on solving the differential energy balance of incompressible flow without heat generation.4,17,18 The general expression to calculate the heattransfer coefficient is given by

hH,wl ) 0.332Re1/2Pr1/3

(3)

To account for a wavy interface influence, eq 3 is multiplied by the factor determined by Kapitza (1.21).

hH,wl ) 0.402Re1/2Pr1/3

(4)

2.4. Model Dimension. The problem of model dimensionality is important especially when the system is distributed. Mathematical models describing falling films are usually two-dimensional, nonlinear, and very complex.10,11,19 In this work, it is enough to represent the spatial (distributed) behavior of the film, but the dimension of its distributed nature (axial or plate direction or both axial and transverse directions) must be established. Thus, an experiment is built and tested to determine the number of spatial dimensions. The experimental device, shown in Figure 2, consists of a tin plate with a dimple. The tin plate is mounted onto a glass plate and is heated from the opposite side by a 1200 W heating element. The temperature of the plate is controlled to within (1 K deviation from the set point. The dimensions of the tin plate and dimple are shown in Figure 3. A flow distributor is placed at the top of the plate to create a uniform falling film. The liquid is recirculated using a centrifugal pump. A tracer is used to determine if there are any significant velocity components in the transverse direction. The construc-

Ind. Eng. Chem. Res., Vol. 43, No. 25, 2004 8119

Figure 3. Experimental setup to determine the flow distribution. The dimensions of the dimpled plate are in millimeters. The size of the dimple is comparable with the size of the dimples found on the plate. Table 1. Experimental Conditions To Determine Model Dimensions temp, °C

fluid

Γ, kg/s‚m

µ, mPa‚s

Re

24 60 24 60

water water SCMC SCMC

0.819 0.819 0.916 0.916

0.921 0.470 5237 2821

3557 6970 0.699 1.298

tion permits the introduction of the tracer in either a heated or a nonheated falling film. The experiment is performed at two different temperatures, 24 and 60 °C. Two different fluids are used, water and a water solution of sodium carboxymethylcellulose (SCMC). The solutions of SCMC are very viscous even at low concentrations (about 10 000 mPa‚ s for one or two solutions). The experimental conditions for each of the four tests are listed in Table 1. The viscosity of the water is calculated from tabulated data, while that of the SCMC is measured using a CannonFenske A919-500 viscometer at room temperature. The viscosity of SCMC at 60 °C is determined using the data provided in ref 20. The results of the experiment are shown in Figure 4. The results indicate that the velocity component in the transverse (horizontal) direction is negligible in the case of SCMC; evidence of the tracer can be found only in the vertical direction. In the case of water, at a Reynolds number of 3557, there is a small dispersion of the tracer in the transverse direction. At a Reynolds number of 6970, there is significant dispersion of the tracer in the transverse direction especially when the tracer is introduced at the edge of the dimple. Fortunately, such high Reynolds numbers are not encountered in the system studied. On the basis of the experimental results, it is concluded that a one-dimensional fallingfilm model is satisfactory. 2.5. Evaporator Design. Falling-film plate evaporators consist of a shell, a heating element, and a vapor space that can be either above or below the heating element.21-23 In this study, the vapor space is positioned over the heating element (the plate stack).

Fresh feed enters the evaporator at the recirculation pump, where it is mixed with the evaporator inventory. A portion of the liquid product stream is recycled to the evaporator by way of a distributor, which ensures full usage of the available heat-transfer surface. As the fluid flows down the plates, the liquid in the fluid evaporates as a result of heating supplied by steam that condenses on the opposite side of the plate. The released vapors travel upward and exit the plate stack. To prevent entrainment in the vapor stream, a separator is installed. The concentrated fluid exits the plate stack and collects in the evaporator inventory. There is a significant dependence of the working fluid parameters (viscosity, heat capacity, enthalpy, specific heat of evaporation, etc.) on the solids concentration and the temperature. These parameters are not constant along the plate length. This change is accounted for by using experimental correlations found in refs 23-26. More details on the actual correlations can be found in ref 27. From the evaporator design, the following distinct sections can be identified: (i) liquor distributor; (ii) plate stack; (iii) evaporator inventory; (iv) inline mixing before recirculation; (v) stream splitting after the recirculation pump. 2.6. Liquor Distributor. The main purpose of the liquor distributor is mechanical; the real process of interest is the possibility of flashing by the liquor. If flashing occurs, it is assumed to be instantaneous. The amount of evaporated water and the changes in the fluid parameters (solids concentration, temperature, and mass flow rate) must be known to determine the input conditions to the plate stack. A steady-state mass balance at the distributor is given by

Gout ) Gin - Wf

(5)

where Gout is the mass flow rate of the liquor leaving the distributor, Gin is the recirculation mass flow rate, and Wf is the mass flow rate of the water vapors. Assuming that the vapor is saturated, the steadystate energy balance is given by

Gincp,inTin ) Goutcp,outTboil + Wf(∆Hevap + cp,wTboil) where cp,in, cp,out, and cp,w are the heat capacities of the recirculating and flashed fluid streams and water, respectively, Tin is the temperature of the recirculating fluid stream, Tboil is the boiling temperature of the fluid stream, and ∆Hevap is the heat of evaporation of the fluid stream at Tboil. The heat capacities of the fluid are functions of the solids concentration and the stream temperature. The boiling temperature is a function of the solids concentration and the saturation temperature. Determination of the fluid stream’s parameters after flashing requires solving a system of nonlinear algebraic equations. Because the fluid of interest is black liquor, the boiling point rise is small and can be neglected.23 The saturation temperature can be used in place of the boiling temperature. In the countercurrent process studied here, flashing is likely to occur only in the evaporators in which the feed stream is supplied, and in these evaporators, the boiling point rise is negligible because the solids concentration in the fluid is low. The change

8120

Ind. Eng. Chem. Res., Vol. 43, No. 25, 2004

Figure 4. Tracer patterns. Water at 24 °C: Re ) 3557. Water at 60 °C: Re ) 6970. SCMC at 24 °C: Re ) 0.699. SCMC at 60 °C: Re ) 1.298.

in the heat capacity is negligible; thus, the steady-state energy balance becomes

GincpTin ) GoutcpTsat + Wf(∆Hevap + cp,wTboil) (6) 2.7. Plate Stack. 2.7.1. Evaporation. A fundamental model of the evaporation process that occurs on a single plate was presented in part 1 of this work.6 For completeness, they are repeated here.

∂G ∂G + vz ) -W′vz ∂t ∂z ∂Xs ∂Xs Xs + vz ) W′vz ∂t ∂z G

(7)

where G is the fluid mass flow rate, Xs is the solids concentration in the fluid, vz is the fluid average velocity, and W′ is defined as W′ ) W/∆z, in which W is the vapor mass flow rate and ∆z is the length of the differential volume element. 2.7.2. Heating. The differential volume element with the important flows and variables is shown in Figure 5. The following assumptions have been made: 1. The average film velocity, vz, in the vertical direction is independent of the horizontal position parallel to the plate. Thus,

∂vz/∂x ) 0

(8)

Figure 5. Differential volume of the heating film.

2. There is no nucleate boiling.6 3. The parameters of the fluid of interest, black liquor, do not change significantly in the differential volume. 4. Ideal mixing is assumed. Thus, the heat of dilution of the fluid can be neglected. It has been shown that if the heat of dilution is neglected, the total error in the evaporation energy requirements in the concentration of softwood black liquor (from 20% to 65-90% solids concentration) is about 4-6%.25,26 Because the change

Ind. Eng. Chem. Res., Vol. 43, No. 25, 2004 8121

in the solids concentration for a single evaporator is about 10%, the heat of dilution can be neglected. 5. The pressures at the vapor and steam sides of the plate do not change significantly with respect to friction losses. According to Perry et al.’s28 section 11-110, the pressure drop in a falling-film evaporator is negligible and the pressure over the falling film is essentially the same as the vapor head pressures. Energy Balance. The energy balance of the differential volume element is given by

(| [|

dEdv 1 ) Q - Ws + Gl|z hl z + vz2 + gz dt 2 1 Gl|z+∆z hl z+∆z + vz2 + g(z + ∆z) (9) 2

)

]

where Gl and hl are the mass flow rate and specific enthalpy of the fluid, respectively, g is the gravitational constant, and Ws is the shaft work. Using the specific enthalpy relation Figure 6. Differential volume of the condensing film.

hl ) cp,l(T - Tr) where Tr ) 0 is the reference temperature and cp,l is the heat capacity of the process fluid and assuming that density changes in the differential volume element are small, shaft work is negligible, and the kinetic and potential energy terms can be neglected. Therefore, eq 9 becomes

Energy Balance. The overall energy balance of the differential volume is given by

( |

dEdv 1 ) -Q - Ws + Gc|z hc z + vz2 + gz dt 2 1 1 Gc|z+∆z hc z+∆z + vz2 + g(z + ∆z) + Wc hst + vv,z2 + 2 2

[ |

) [

]

]

(10)

g(z + ∆z) (15)

The mass flow rate of the liquor is Gl ) aδFvz, where δ is the film thickness, a is the plate width, and F is the density of the fluid. The total energy of the differential volume can be assumed to be equal to its internal energy. The internal energy is given by

A material balance of the differential volume is given by

dEdv/dt ) Q + Glcp,l(T|z - T|z+∆z)

dMdv/dt ) Gc|z - Gc|z+∆z + Wc

(16)

(11)

Substitution of the above equation into eq 15 and neglecting the kinetic and potential terms give

where Mdv is the mass of the differential volume element. The heat, Q, entering the differential volume can be expressed as

dMdv dEdv ) -Q + hc|z + Gc|z+∆z(hc|z - hc|z+∆z) + dt dt Wc(hst - hc|z) (17)

U ) Mdvcv,l(T - Tr) ) aFδ∆zcv,lT

Q ) hHa∆z(Tw - T)

(12)

where hH is the heat-transfer coefficient of sensible heating and Tw is the wall temperature. Using eqs 11 and 12, eq 10 becomes

aFδ∆zcv,l

dT ) hHa∆z(Tw - T) + aδFcp,lvz(T|z dt T|z+∆z) (13)

Taking the limit as ∆z f 0 gives

∂T hH(Tw - T) ∂T + vz ) ∂t ∂z Fδcp,l

(14)

where it is assumed that cv,l ≈ cp,l. (Recall that cp - cv ) R, where R is the universal gas constant.) 2.7.3. Steam Condensation. The model describing the condensation process is similar to that developed for the evaporation process. Figure 6 shows the differential volume with the important flows.

The specific enthalpy is given by

h ) cp(T - Tr) Substitution of the above equation into eq 17 and assuming that the material accumulation term is negligible give

dEdv/dt ) -Q + Gc|z+∆z(cp,cT|z - cp,cT|z+∆z) + Wc(hst - hc|z) (18) where cp,c is the specific heat capacity of the condensate. Only a phase change occurs along the plate; thus, the temperature of the film does not change. The heat of condensation (evaporation), ∆Hc, is given by the difference hst - hw|z. It then follows that the energy balance is given by

dEdv/dt ) -Q + Wc∆Hc

(19)

As in the case of the evaporation model, the energy balance is solved at steady state. Thus

8122

Ind. Eng. Chem. Res., Vol. 43, No. 25, 2004

Q ) Wc∆Hc

(20)

Q ) hca∆z(T - Tw,c)

(21)

where

Material Balance. Applying the principle of material conservation to a differential volume of size, ∆z (see Figure 6), gives

M|t+∆t - M|t ) Gc|z∆t - Gc|z+∆z∆t + Wc∆t where M is the mass of the condensate, Gc is the mass flow rate, W is the mass flow rate of the condensing steam, and ∆t is the time difference. Dividing through by ∆t and taking the limit as ∆t f 0 yield

∂M/∂t ) Gc|z - Gc|z+∆z + Wc

(22)

It is assumed that the residence time τ of the condensate in the differential volume is constant. Thus, the residence time can be expressed as

τ ) ∆z/vz

(23)

Figure 7. Wall differential volume element.

is achieved instantaneously. It follows that eq 26 becomes

Qcw ) Qwe

Substitution of the above equation into eq 22 gives

τ

∂Gc ∆z ∂Gc ) ) Gc|z - Gc|z+∆z + Wc ∂t vz ∂t

(24)

which represents the mass balance with respect to the condensate mass flow rate. Define

(27)

The following equations describe the heat transferred at the wall due to condensation/heating and condensation/evaporation, respectively.

hC(Tf,c - Tw,c) ) hH(Tw,h - Tf,h) )

δw (T - Tw,h) λw w,c (28)

hC(Tf,c - Tw,c) ) hE(Tw,e - Tf,e) )

δw (T - Tw,e) λw w,c (29)

W′c ) Wc/∆z Dividing eq 24 by ∆z and taking the limit as ∆z f 0 yield

∂Gc ∂Gc + vz ) W′cvz ∂t ∂z

(25)

2.7.4. Heat Transfer at the Wall. Because the wall does not exchange any mass with its surroundings, only an energy balance is necessary (see Figure 7). The boundaries through which the wall differential volume can exchange heat are as follows: (i) heat flux from the condensate film to the element, Qcw, (ii) heat flux from the element to the evaporating film, Qwe, (iii) heat flux to the top adjacent element, and (iv) heat flux to the bottom adjacent element. The fluxes to the adjacent elements can be neglected (the temperature gradient in the z direction is negligible because of the phase change on the both sides of the wall). The significant heat fluxes are those from and to the films at both sides of the wall. Thus, the overall energy balance of a differential volume at the wall, neglecting potential and kinetic energy terms and shaft work, is given by

dEw/dt ) Qcw - Qwe

These equations are used to calculate the surface temperatures at the wall. To summarize, the equations that describe heating of the fluid are given by eqs 14, 25, and 29, and those that describe evaporation are given by eqs 7, 25, and 28. The above mathematical model can be transformed into a dimensionless form by introducing the following dimensionless variables:

T T0

T* )

z* )

z L

T/we )

v/z )

Twe T0

vz vf,0

G* )

/ vc,z )

vc,z vf,0

G G0

G/c )

t* )

Gc Gc,0 X* )

Xs Xs,0

W′* )

W′ W′0

tvf,0 L

(26)

In the current study, the plate is made of stainless steel with a thickness of 1.5 mm. According to ref 1, the thermal response of thin (1.15 mm) metal (copper) walls is on the order of 0.01 s. Thus, it is possible to solve the energy balance at the wall assuming that equilibrium

W′/c )

W′c W′c,0

η)

LW′0 G0

ηc )

LW′c,0 Gc,0

ηT )

h HL Fcpδv/f,0

where T* is the dimensionless fluid temperature, T/we is the dimensionless wall temperature on the evapora-

Ind. Eng. Chem. Res., Vol. 43, No. 25, 2004 8123

tion side of the plate, G* and G/c are the dimensionless mass flow rate of the fluid and condensate, respectively, X* is the dimensionless solids concentration, x* is the / are the dimensionless spatial variable, v/z and vc,z dimensionless fluid and condensate velocities, respectively, and t* is dimensionless time. Substitution of these dimensionless variables in the model equations leads to the following dimensionless system:

dL NGp - Gt,out L dF ) dt FAt FAt dt

(38)

A mass balance of the solids in the fluid is given by N

dMs/dt )

∑i Gs,p - Gs,out

(39)

∂T* ∂T* + v/z ) ηT(T/w - T*) ∂t* ∂z*

(30)

∂G/c ∂G/c + v/z ) W′c/v/z ηc ∂t* ∂z*

where Ms is the mass of the solids in the tank, Gs,p is the mass flow rate of the solids exiting each plate, and Gs,out is the mass flow rate of the solids exiting the tank. By analogy with the overall mass balance, the above can be rewritten as

(31)

dMs/dt ) NGs,p - Gs,out

∂G* ∂G* + v/z ) -W′*v/z η ∂t* ∂z*

(32)

∂X* X* ∂X* + v/z ) W′*v/z η ∂t* ∂z* G*

(33)

An analysis of the dimensionless parameter η shows that when G0 f 0, η f ∞ and when G0 f ∞, η f 0. Thus, at high feed flow rates, evaporation will be negligible, while at very low feed flow rates, evaporation will be almost instantaneous. This is consistent with the physics of the process. In the case of ηc, the mass flow rate at the top of the plate is zero because there is no condensate at the top. In this case, Gc,0 is a sufficiently small number (0.001 kg/s). 2.8. Evaporator Inventory. The evaporator inventory is located at the bottom of the evaporator and can be considered a tank with a constant area and volume. The level is well controlled to maintain smooth operation of the centrifugal pump (see Figure 1). 2.8.1. Mass Balance. The overall mass balance is given by N

dM/dt )

∑i Gp,i - Gout

(34)

where M is the mass of the liquor in the tank, Gp,i is the mass flow rate at the output of each plate, Gout is the mass flow rate at the exit of the tank, and N is the number of plates. The outflow from each plate is assumed to be identical; thus, the above balance can be rewritten as

dM/dt ) NGp,i - Gt,out

(35)

where N is the number of plates in the stack. A material balance of the tank is given by

(40)

The mass of the liquor can be defined in terms of its density, M ) FAtL. Define X to be the mass fraction of the solids. Using eq 35, it then follows that eq 40 can be rewritten as

dX NGp ) (X - X) dt FAtL p

(41)

where Xp is the solids concentration of the fluid exiting the plates. 2.8.2. Energy Balance The energy balance of the liquor in the evaporator inventory is given by

dE d(FAtcpT) ) ) Qin - Qout ) NGpcpTp - Gt,outcpT dt dt (42) Expanding the derivative at the right-hand side of this equation yields

dT NGpTp - Gt,outT T dF T dL T dcp ) (43) dt FLAt F dt L dt cp dt where the rate of change of the heat capacity is given by

dcp dT dX dX ) 3.822X + 3.822T - 3322.67 dt dt dt dT

(44)

To summarize, the mathematical model describing the evaporator inventory is given by eqs 38, 41, and 43. 2.9. Inline Mixer. Mixing of the fresh feed with the recycled fluid that occurs before the recirculation pump is assumed to be instantaneous and ideal (no heat effects due to mixing). Also, fluid holdup in the inline mixer is assumed to be negligible. The overall mass and energy balances of the mixer are given by

Gm,out ) Gfeed + Gt,out

(45)

(36)

Gm,outTm,outcp,m,out ) GfeedTfeedcp,feed + Gt,outTt,outcp,t,out (46)

where At is the tank cross-sectional area and F is the density of the fluid in the tank

respectively. A material balance of the solids is given by

d(LAtF)/dt ) NGp - Gt,out

dF dT dX ) -0.495 + 600 dt dt dt

(37)

Expanding the time derivative at the left-hand side of eq 36 and rearranging the equation give

Gm,outXm,out ) GfeedXfeed + Gt,outXt,out

(47)

2.10. Splitter. Because the holdup is negligible, only the steady-state mass balance is necessary to represent the division of the flow at the splitter. Thus

8124

Ind. Eng. Chem. Res., Vol. 43, No. 25, 2004

Gspl ) Grec + Gprod

(48)

3. Model Development of a Multiple-Effect Evaporator Plant 3.1. Description of the MEV Plant. Usually more than one evaporator is necessary to concentrate the solids to the desirable value. Thus, a MEV plant consists of evaporators connected in series with respect to their vapor lines.21,23 The vapor (steam) produced during the evaporation process by one evaporator is used as the heating source for the next evaporator in the series. Fresh steam is usually supplied to the first evaporator in the MEV. The driving force for the operation of the MEV plant is determined by the overall temperature difference between the steam temperature and that of the heat sink (cooling water).21 The overall temperature difference is reduced by the boiling point rise. The latter occurs because of an increase in the solids concentration. An important measure of the MEV plant efficiency is a parameter called steam economy,21,23 which is defined as the mass of water evaporated per unit mass of steam or, in the general case, the mass of liquid evaporated per unit of heat supplied. In this work, steam economy (kilograms of evaporated water per kilograms of steam) will be used to evaluate the MEV plant performance. Steam usage is affected by many factors including the evaporator type, process fluid and steam characteristics, condensate flashing, and liquid-vapor flow sequence. There are two basic liquid-vapor sequences: forward (cocurrent) and backward (countercurrent) flow.21,23 In the forward flow scheme, the cold feed is the input to the evaporator with the highest temperature. The advantage of this scheme is that pumps are not needed to move the process fluid from evaporator to evaporator because the fluid is transported by the natural pressure gradient in the MEV plant. However, significant sensible heating should be applied to the process fluid to reach the boiling temperature. In the case of backward flow, the feed stream does not require significant preheating; however, the process fluid has to be pumped from evaporator to evaporator. Pumping the process fluid is not a disadvantage in the case of a falling-film-type evaporator because the recirculation pumps are inherently part of the design. Therefore, the reverse flow scheme is preferred when falling-film evaporators are used. When the working fluid is concentrated, it is common practice to include superconcentrators as part of the evaporator plant. In these unit operations, the highest solids concentration is achieved. The superconcentrators are usually two or three falling-film evaporators mounted in one unit and connected in series with respect to the flow path of the working fluid. Each evaporator in the superconcentrator has its own fresh steam supply, but they all share the same steam header and secondary vapor space. A typical modern black liquor multiple-effect fallingfilm evaporator plant has a sequence of evaporators and at least one superconcentrator. In this study, the MEV consists of four falling-film plate-type evaporators in series followed by one superconcentrator that has three falling-film plate-type evaporators mounted in a single unit. A diagram of the MEV plant is shown in Figure 8. The feed stream is split equally between evaporators E-4 and E-5. This is normal practice in a black liquor

evaporation plant23 in an attempt to maintain an even distribution of the throughput. A total of 55 lb. (55 psig) of steam is used to heat the superconcentrator. Vacuum conditions are achieved by a surface condenser followed by a steam ejector (EJ-1) for removal of the noncondensable gases. The auxiliary flash units, F-1-F-4, are common in a MEV plant. Their purpose is to improve steam economy by flashing the condensates of the secondary vapors produced by the evaporators. The condensates that come from the fresh steam sources are returned to the boiler because flashing these condensates results in increased boiler feedwater usage and hence higher utilities. In this work, all condensates (primary and secondary) are flashed in order to achieve the highest possible steam usage. Modeling the MEV plant requires models of each evaporator unit, the pressure dynamics, and condensate flashing processes. Models of a single evaporator and condensate flashing were already developed and presented. A model that describes the pressure dynamics follows. As pointed out in section 2.5, the working fluid parameters are not constant. To account for the varying parameter values among the evaporators in the MEV plant, the same correlations as those used in the case of a single evaporator also are employed. 3.2. Pressure Dynamics Model. The model should describe the pressure variations in each vapor/condensate space shared between two evaporators. This space comprises the vapor space of a given evaporator and the condensate space of the next evaporator in the series. Therefore, the pressure in this closed space depends on the rates of evaporation and condensation in the corresponding spaces. The following assumptions are made: (i) Water vapor behaves as an ideal gas. (ii) Pressure losses are negligible during vapor transport between the bodies. (iii) Condensation due to heat losses to the surroundings is negligible. With these assumptions, the mass balance of the system is given by

dM/dt ) W - Wc

(49)

where W is the mass flow rate of the evaporated water in the source evaporator and Wc is the mass flow rate of the condensing vapors in the next evaporator. Using the ideal gas law, the pressure dynamics is given by

(

M dP d ) RT dt dt MmV

)

(50)

where Mm is the molecular weight of water (18 g/mol), V is the volume of the system (a design parameter) and is equal to the sum of the volume of the vapor space of the source evaporator, the volume of the piping between the evaporators, and the volume of the condensate space of the next evaporator. The temperature is assumed to be the saturation temperature (Ts) of water vapor at the given pressure. The volume of the vapor space is assumed constant; thus, expansion of the right-hand side of the above equation gives

(

dTs dP R dM dP ) Ts + M dt MmV dt dP dt

)

(51)

The term dTs/dP is very small (∼10-5) and can be omitted from the above equation. Substituting eq 49 into

Ind. Eng. Chem. Res., Vol. 43, No. 25, 2004 8125

Figure 8. Multiple-effect falling-film evaporator plant. There are four falling-film evaporators (E-1-E-4) and one superconcentrator with three evaporators (SC-1-SC-3). Table 2. Evaporator Operating Conditions

Table 3. Evaporator Sensitivity to Disturbances

variable

value

change, %

feed mass flow rate, Gf feed solids, Xf feed temperature, Tf vapor pressure, Pvap steam pressure, Ps

58 kg/s 0.145 kg/kg 303.15 K 24 000 Pa 47 000 Pa

(5 (5 (5 (9.355 (5

the above equation gives the following expression for the pressure dynamics:

RTs dP ) (W - Wc) dt MmV

(52)

4. Results and Discussion 4.1. Single Evaporator. The partial differential equations (PDEs) are discretized using orthogonal collocation on finite elements. The scheme used has 10 finite elements with 2 internal collocation points, which gives a total of 31 spatial points. In 21 points, the solution is calculated by integration of 21 ordinary differential equations for each PDE. In the case of heating, there are 93 differential-algebraic equations (DAEs), and for evaporation, there are 124 DAEs. The solver used to solve the discretized systems of PDEs is the LSODI (Livermore solver for ordinary differential equations, implicit systems) solver, included in the package ODEPACK.29-31 The execution time is approximately 14 h on a 1.2 GHz AMD Athlon node, which is equivalent to about 1 h of real time. The operating conditions of the evaporator are listed in Table 2. The solids concentration of the product stream is the most important variable in the evaporation process. The sensitivity of this parameter to unmeasured disturbances and changing operating conditions is shown in Figures 9-11. Define the sensitivity of the evaporator to the disturbances as

sensitivity ) ∆SS/∆d

(53)

where ∆d is the percent change in the disturbance and ∆SS is the percent change in the solids concentration with respect to the steady-state value. Table 3 lists the sensitivity values to each disturbance. 4.1.1. Feed Flow Rate. The evaporator is less sensitive to changes in the feed flow rate when com-

sensitivity disturbance

∆d > 0

∆d < 0

feed mass flow rate, Gf feed solids, Xf feed temperature, Tf vapor pressure, Pvap steam pressure, Ps

-0.444 1.105 1.431 -0.945 0.911

0.523 -1.086 -1.188 1.157 -0.750

pared to changes in the other variables listed in Table 2. The direction of the response is as expected (see Figure 9): an increase in the production rate leads to a decrease in the product stream solids concentration because more feed is being processed but the amount of energy supplied is unchanged. The magnitude of the response is similar but not the same for the increase and decrease in the mass flow rate. 4.1.2. Feed Concentration. The evaporator performance is affected by changes in the feed solids concentration (see Figure 9). The solids concentration in the product stream follows the same direction as increases and decreases in the feed solids concentration. 4.1.3. Feed Temperature. Changes in temperature have the largest effect on the performance of the evaporator, especially increases in the temperature (see Figure 9). The nonlinearity of the system is well expressed because the changes in the product solids concentration differ by about 30% for increases and decreases in the feed temperature. 4.1.4. Steam Pressure and Vapor Pressure. Figures 10 and 11 show the responses to changes in steam and vapor pressure, which can be explained from the physics of the process. An increase in the steam pressure leads to larger heat transfer (greater temperature difference across the plate) and therefore greater evaporation, which leads to an increase in the solids concentration of the product. Changes in the vapor pressure follow a similar logic. For all operating changes, the simulated responses of the evaporator are smooth. There are no significant process lags from the time the operating condition changes are applied to the time the changes affect the product stream. This is because the mixer, pump, and splitter sections contain no significant dynamics (see Figure 1). The long-term behavior (low frequency) reflects the heat- and mass-transfer phenomena that are affected by the changes in the operating conditions.

8126

Ind. Eng. Chem. Res., Vol. 43, No. 25, 2004

Figure 9. Product solids concentration: (a) 5% increase and (b) 5% decrease in the feed mass flow rate; (c) 5% increase and (d) 5% decrease the feed solids concentration.

The difference between the effect of the throughput change and the feed concentration change follows from how these changes affect the heat-transfer coefficient of sensible heating. For the operating conditions provided in Table 2, the change in the heat-transfer coefficient of sensible heating is approximately doubled (4% as compared to 2%) for the same changes in the feed solids concentration and the throughput. 4.2. MEV Plant. The model equations are integrated using the same solver as that in the case of a single evaporator. Because of the computational burden, parallel processing (message passing interface or MPI) is used. Each evaporator is assigned to a node (1.2 GHz AMD Athlon of a Rocks cluster.2 [Rocks is an opensource high-performance Linux cluster solution.] A total of 1 h of real operation takes 24 h of computation on the cluster.

The operating conditions of the MEV plant are listed in Table 4. The diameter of each body varies between 4 and 7 m, and their height is about 12 m. Each body contains between 130 and 200 plates. [The actual parameters were supplied by Tembec, Inc. (Montreal, Quebec, Canada). For confidential reasons, the true values cannot be published.] The levels of evaporators E-2-E-5 and the superconcentrator SC-1 are controlled using single-loop feedback proportional-integral (PI) controllers. The responses in the face of disturbances are shown in Figures 12-15.All of the responses converge in a smooth fashion. The levels of all evaporators (not shown) are regulated to within (0.2% of their nominal values. In each case, the MEV plant reaches a new steady state in approximately 1 h after the disturbance is introduced. In the case of E-5, the new steady state is achieved over a longer

Ind. Eng. Chem. Res., Vol. 43, No. 25, 2004 8127

Figure 10. Product solids concentration: (a) 5 increase and (b) 5% decrease in the feed temperature; (c) 5% increase and (d) 5% decrease in the steam pressure. Table 4. Operating Conditions of the MEV Plant parameter

value

Gfeed Xfeed Tfeed Psteam (SC-1-SC-3) Pcondenser

116 kg/s, 50/50 split 0.145 kg/kg 303.15 K 410 kPa 21.4 kPa

length of time as compared to the time observed in the responses of a single evaporator (about 20 min; see Figures 9-11). This is attributed to process feedback by the secondary vapors that are absent in the single evaporator case. The effect of the disturbances on the solids concentration at SC-1 is about 5 times greater than that observed at E-5. This behavior will be explained in the following section (also see Figure 16).

4.2.1. Sensitivity to Disturbances. The sensitivity of the MEV to the set of disturbances is assessed with respect to the product solids concentration and steam economy. Mass Balance. The sensitivity of the product solids concentration is provided in Table 5. The values are calculated using eq 53. In the case of a single evaporator, it was found that the feed temperature disturbance had the strongest effect on the product solids concentration. In the case of the MEV, the largest sensitivity is associated with changes in the feed mass flow rate (throughput). One explanation for this behavior is the nonlinear dependence of the product solids concentration on the amount of water in the process fluid. Consider the following example where the feed stream has a 0.2 kg/kg solids concentration. Figure 16 shows the amount of water that must be evaporated to achieve

8128

Ind. Eng. Chem. Res., Vol. 43, No. 25, 2004

Figure 11. Response of the product solids concentration to changes in the secondary vapor saturation temperature. Top: 9.355% increase. Bottom: 9.335% decrease.

Figure 12. Response of the MEV plant to changes in the feed flow rate. Top: 5% decrease. Bottom: 5% increase.

a different product solids concentration for two different feed mass flow rates. Observe that this relationship is nonlinear. When the product solids concentration is in the range of 0.27-0.3 kg/kg, the amount of water to be evaporated for each feed rate is not significantly differ-

ent. However, if the range of solids concentration is between 0.6 and 0.7 kg/kg, then the amount of water to be evaporated is significantly different between the two feed rates. For a 8.9% difference between the feed rates, there is a 13.25% increase in the amount of water to be

Ind. Eng. Chem. Res., Vol. 43, No. 25, 2004 8129

Figure 13. Response of the MEV plant to changes in the feed solids concentration. Top: 5% increase. Bottom: 5% increase.

evaporated for a 0.6 kg/kg solids concentration. It can be concluded that a larger final solids concentration means a greater sensitivity to throughput changes. Energy Balance. The sensitivities are calculated as the differences between the steam economy at nominal (ideal) conditions and the steam economy after the disturbances (nonideal) are introduced. Steam economy (SE) sensitivity is calculated as follows: SC-1

SE )

W ∑ E-5

SC-1

(54)

∑ Wc

SC-3

where W is the mass flow rate of the evaporated water and Wc is the mass flow rate of the condensate. The sensitivity values of the steam economy to each disturbance are provided in Table 6. The energy balance shows the largest sensitivity to changes in the feed temperature. It is observed that a decrease in the steam economy occurs when the feed temperature decreases because more steam is consumed to heat the liquor to its boiling point. Also the steam economy decreases with an increase in the feed mass flow rate because there is more throughput to be processed. The disturbance in the heat transfer does not appear to affect the steam economy. However, changes in the

heat-transfer capability result in a lower product solids concentration. Thus, if the same solids concentration in the product stream is to be achieved, then the steam economy must decrease because more steam is consumed. Although changes in the steam economy may appear to be small, the analysis must consider the effect over the entire duration of the operation. For example, if the MEV plant is processing 116 kg/s of process fluid (solids feed concentration of 0.145 to a product solids concentration of 0.73 kg/kg), the amount of steam is correspondingly 24.92 kg/s when the feed temperature is decreased by 5% with respect to the base case and 22.45 kg/s when the feed temperature is increased. If the cost of steam is assumed to be $5 USD/ton, then the per year difference is ∼$400 000 USD, which is not an insignificant number. 4.2.2. Model Validation. Industrial steady-state data for a seven-effect falling-film evaporator plant with similar physical characteristics, dimensions, and operating conditions were provided by a mill that is owned and operated by Tembec, Inc. (Montreal, Quebec, Canada). On the basis of these parameters, simulation of the mathematical model described above produced a product solids concentration that differed from the data by ∼0.006 kg/kg. Although this is not definitive, it demonstrates the potential of the model to predict the product concentration for different operating conditions. Also, the model assumptions and empirical and semiempirical correlations appear to be reasonable. Because no

8130

Ind. Eng. Chem. Res., Vol. 43, No. 25, 2004

Figure 14. Response of the MEV model to changes in the feed temperature. Top: 5% decrease. Bottom: 5% increase.

Figure 15. Response of the MEV plant to a 10% decrease in the heat transfer of the superconcentrator.

dynamic data are available for further validation, further comparisons cannot be made at this time. However, the steady-state validation together with the fact that the observed responses follow the physics of the process indicates the potential of the fundamental model to be used for other studies such as design modifications, optimization, process monitoring, and real-time control. 5. Summary This work developed a fundamental distributed parameter model of a multiple-effect falling-film evapora-

tor plant. The model describes the important phenomena of evaporation, heating, and condensation for different hydrodynamic regimes as well as the pressure dynamics of the MEV plant. Experimental investigations were performed to determine the model dimensions. Because the recovery boiler requires the concentrated fuel (product of the MEV) to have a consistent solids concentration within the specified range,27 the sensitivity to different but expected disturbances and process conditions was investigated using the fundamental model. It was found that the single evaporator and the

Ind. Eng. Chem. Res., Vol. 43, No. 25, 2004 8131

ful for the data and discussions with Tembec, Inc., at St Francisville (U.S.) and Crestbrook (Canada). Nomenclature

Figure 16. Response of the MEV model to a 10% decrease in the heat transfer of the superconcentrator. Table 5. Sensitivity of the Product Solids Concentration sensitivity disturbance

d>0

d