A Distributed-Parameter Model of Black Liquor Falling Film

Ind. Eng. Chem. Res. 2003, 42, 1925-1937

1925

A Distributed-Parameter Model of Black Liquor Falling Film Evaporators. Part I. Modeling of a Single Plate Z. Stefanov and K. A. Hoo* Department of Chemical Engineering, Mail Stop 3121, Texas Tech University, Lubbock, Texas 79409-3121

Phenomenological modeling of black liquor evaporators, found in the pulp and paper industry, presents a very challenging task. The physicochemical phenomena that occur do not lend themselves readily to known mechanisms, and in many instances, data to support these hypotheses are difficult to obtain. In this work, a distributed-parameter model (a system of partial differential equations) based on first-principles knowledge about the fluid dynamics and heattransfer processes is developed for a falling film lamella type evaporator. Primarily, the model describes falling film evaporation on one lamella. The model is solved using orthogonal collocation on finite elements in the presence of scaling and disturbances in the mass feed rate, feed dry solids content, and wall temperature. Introduction The black liquor evaporation process is a very important part of every pulp mill. This process produces concentrated or strong black liquor. The black liquor is a residual liquid from the wood cooking process that consists of organic and inorganic compounds dissolved or suspended in water, usually with 10-14% dry solid content.1,2 The organic part comes from the wood and consists of a complex mixture of the products of the reaction of lignin with cooking chemicals and the destruction of the sugars. The actual composition of this mixture is not well-known. The inorganic part is mainly sodium hydroxide, sodium sulfide, sodium sulfate, and sodium carbonate that come from the cooking chemicals and other mineral salts natural to wood. The strong black liquor, which is the raffinate of the evaporator process, is used as fuel in the recovery boiler (see Figure 1). The recovery boiler, a combination of a power boiler and a chemical reactor, plays a vital role in mill operation because this is where the first step in the regeneration of the cooking chemicals takes place, along with the production of a significant amount of the energy that impacts the overall mill economy. Although many factors can lead to unsatisfactory performance of the recovery boiler, the elimination of inconsistencies in the feed composition (dry solids content) and feed rate can result in improved performance. Another reason for better understanding of the evaporation process is the impact the recovery boiler and the evaporator plant might have on the environment.3,4 (The evaporator plant usually consists of several evaporators or “effects” connected in series. The actual number of evaporators is a function of the optimal economic parameters.) The evaporators are one of the main sources of total reduced sulfur (TRS) compounds that are restricted by the U.S. EPA. The source of these compounds in the evaporator plant is the foul condensates. Unstable operation of the recovery boiler results in increased sulfur oxides and sodium carbonate and sulfate carryover. From an economic perspective, the carryover means significant chemical losses. * To whom correspondence should be addressed. E-mail: [email protected]. Tel.: (806) 742-4079. Fax: (806) 742-3552.

Figure 1. Simple illustration of the evaporator plant and the recovery boiler.

This motivates the development of a rigorous model of the evaporation process that can be used to study the impact of different operating conditions, optimize operating conditions, investigate design changes, and promote more sophisticated control strategies. The evaporator models that can be found in the open literature are mostly lumped-parameter, steady-state models. These models can be used cautiously for design purposes, but in the case of control applications, dynamic models are required. A first-principles-based distributed-parameter model will be better able to account for the various phenomena that occur in the evaporator. An example is the change in the extent of heat transfer along the length of the plate that results from changes in the black liquor parameters during evaporation. A lumped model cannot correctly represent this distributed phenomenon. Although the derivation of a distributed-parameter model might be straightforward for simple liquid systems with known properties, this is not the case for black liquor systems. As stated above, the exact composition of this liquor is not known. Indeed, the characteristics of the liquor (viscosity, density, surface tension, etc.) at the operating conditions in the evaporator (temperature, dry solids content) lead to Reynolds, Prandtl, and Kapiza numbers for which no experimental data are available. Moreover, no fundamental models that describe black liquor hydrodynamics accurately have been reported. This work focuses on the derivation of a distributedparameter model that can be used for control studies and design changes. Thus, the model is designed to

10.1021/ie020483a CCC: $25.00 © 2003 American Chemical Society Published on Web 03/29/2003

1926

Ind. Eng. Chem. Res., Vol. 42, No. 9, 2003

Figure 2. (Left) One effect (or evaporator) of a multiple-effect evaporator plant. (Right) One plate (or lamella) of the effect.

describe the most important parameters of the black liquor, namely, the mass flow rate and the dry solids content. It is outside the purview of this work to provide detailed modeling of black liquor hydrodynamics (wavy phenomena, breakup of the falling film, etc.) on the plate. However, there is wide interest in the determination of the operating conditions where the regime changes from turbulent to wavy laminar to laminar, as this is critical for proper calculation of the transport coefficients (e.g., the heat-transfer coefficient of evaporation). Many different types of evaporators are available such as long tube vertical (LTV), falling or rising film, horizontal, forced circulation, and falling film plate, to name a few of the more common ones.1,2 The falling film lamella type of evaporator is chosen for this work, as it is the most modern type used in the pulp and paper industry and its design from an energy point of view is superior to that of other types of evaporators.1,2 The paper is organized as follows. Section 2 briefly reviews the relevant transport and fluid dynamic phenomena of falling films. This review is followed by the development of the phenomenological model. Section 4 analyzes the simulated results that are obtained in the presence of some important disturbances, and section 5 discusses the limitations of the proposed model and future work. The final section summarizes the important findings. 2. Transport Phenomena in Falling Films This section first describes the evaporation process of a single plate (lamella) in the falling film evaporator. The remarks are restricted to the processes found in a pulp mill. This discussion is followed by a brief review of some important fluid dynamic phenomena of falling films as they relate to the development of the phenomenological model. 2.1. Black Liquor Falling Film Evaporation Process. The schematic on the left-hand side of Figure 2 represents an evaporator body, and that on the righthand side is a magnification of the important flow streams on a single plate. The black liquor enters at the top of the evaporator body and flows down on one side of the plate. The plate construction is dimpled to promote turbulent flow. When there are several plates, a uniform distribution of the liquor is achieved through the use of distribution devices. Steam is introduced at

the other side of the plate. As the water in the black liquor evaporates, it moves in a direction opposite to that of the black liquor. The concentrated black liquor exits the plate stack (see Figure 2) and is forced to leave at the bottom of the effect by a circulating pump. The vapor that exits the plate stack passes through the impingement separator and becomes the steam feed to the next effect. 2.2. Relevant Fluid Dynamic Phenomena. The transport phenomena that occur in falling films, especially when evaporation is present, are very complex. Numerous studies have been conducted on the fluid dynamics of falling films, starting with the classic Nusselt theory for laminar, flat falling films. However, real films are wavy laminar or turbulent in most industrial applications.5 Thus, it is not surprising that Nusselt’s theory is unable to predict important parameters such as the heat-transfer coefficient. The modeling presented in this paper follows the concepts and models presented by Alhusseini and Chen.5 The surface waves in the wavy laminar regime carry a significant part of the mass flow, and they change the nature of the flow, for instance, causing local recirculation. In the turbulent regime, the nature of the entire flow changes, and eddy momentum transport becomes significant. In both regimes, the heat and mass transport increase. To determine the transition from wavy laminar to fully turbulent flow, Chun and Seban proposed two correlations.6 The first is a function of the Prandtl (Pr) number, whereas the second is a function of the Kapitza number (Ka)6

Ret,t ) 5800Pr-1.06

(1)

Ret,t ) 0.215Ka-1/3

(2)

Ka ≡

µ4 g Fσ3

(3)

The Nomenclature section provides definitions of the symbols used. These two correlations appear to be used widely and are recommended, for instance, by Åsblad and Berntsson.7 However, calculations (see Table 1) at conditions for black liquor mixtures show large differences in their quantitative predictions, but qualitative agreement on the type of flow regime. This lack of quantitative

Ind. Eng. Chem. Res., Vol. 42, No. 9, 2003 1927 Table 1. Estimates of the Reynolds Number Using the Correlations of Chun and Seban6 no. 1 1 2 2 3 3 4 4 5 5

effect position entrance exit entrance exit entrance exit entrance exit entrance exit

Pr 17.42 37.41 24.34 53.20 38.69 91.56 48.13 328.48 127.72 18376

Ka 10-8

3.403 × 3.814 × 10-7 1.198 × 10-7 1.287 × 10-6 5.156 × 10-7 1.062 × 10-5 -

agreement is attributed to the differences between the simple liquids used by Chun and Seban6 to generate the correlations and the more complex black liquor mixture. Because no such correlations are available in the open literature for black liquor mixtures and because the two correlations in eqs 1 and 2 are consistent in their prediction of the type of regime, as a first approach to the phenomenological modeling of this type of falling film evaporator, the existence of a turbulent regime (recall the dimpling of the plates) will be assumed. Hence, the correlations of Chun and Seban6 will be used. One of the most important parameters in falling films, when evaporation is present, is the heat-transfer coefficient (h). Numerous studies attempting to determine this parameter have been published in the open literature. One of the earliest correlations was provided by Nusselt.8 Fundamental research by Wilke in 1962 provided a more precise correlation for this parameter.6,9 Limberg’s analysis of the energy balance equation for turbulent films resulted in a seventh-order eigenvectoreigenfunction solution.9 Limberg also provided a correlation that calculated the turbulent Nusselt number to within 1% of the exact solution. It is not clear, from the physics of the process, why this correlation predicts the Nusselt number. Indeed, Limberg provided no explanation. Recent experimental work has indicated that both regimes, wavy laminar and turbulent, can be present in simple to complex applications that involve a falling film. For a wavy laminar regime, numerical simulations by Miladinova et al.10 and Miyara11 showed the onset and development of wave formation. In a work by Chang et al.,12 the results showed the presence of small waves with aperiodic behavior. All of these studies solved the mass, momentum, and energy balance equations. A general conclusion is that falling film flow is very unstable and sensitive to flow disturbances. In the present study, turbulent behavior occurs for almost all of the available operating data used. Many models that predict the transport processes in these films are available. The early models assume that a turbulent falling film consists of a boundary layer at the wall and a turbulent layer. More recent acceptable theories hypothesize that these films consist of a boundary layer at the wall, a turbulent core, and a second boundary layer at the free surface.5 The presence of the additional boundary layer might explain the overprediction of the models by Chun and Seban6 and Limberg,9 because their models do not account for the damping of turbulence at the free surface. Alhusseini and Chen5 compared four models in terms of the presence of the boundary layer at the free surface and the eddy viscosity in the film.5 They provided a

Re

Ret eq 1

Ret eq 3

regime

3331 1039 1501 539 743 264 412 51 74 0.691

287 125 200 87 122 49 97 13 35 0.18

89 34 57 22 32 11 0

turbulent turbulent turbulent turbulent turbulent turbulent uncertain uncertain uncertain uncertain

{

model for the eddy diffusivity of momentum given by

1 1 + x1 + (2ky*DwF)2 0 e y* < y/i 2 2 (4) ) 1.199 × 10-16 Re2m δ* - y* / ν y e y* < δ* i 2/3 / 2 δ* A -

(

i

)

where /

Dw ) (1 - e-y*/Aw)2

/

/

F ) (e-cy /δ )2

The first expression in eq 4 covers the layer influenced by the wall and the turbulent core. The damping effect of the wall is represented by the van Driest13 damping function Dw, which models the damping action of a smooth wall. The parameter F represents the mixing length correction function in the core region; k is von Karman’s constant, which has a value of 0.41; y* is the dimensionless y coordinate (see Figure 4); and δ* is the dimensionless film thickness. The second expression in eq 4 covers the influence of the free surface. Here, m is an empirical function of the Kapitza number, and A/i is a different empirical function of the Kapitza number that corrects for representing the wavy interface as a smooth surface. In this study, the plates in the evaporators are dimpled to promote artificial turbulence. It is reasonable to assert that the wall is not smooth and remove the damping effect of the wall from the first expression in eq 4. This modification gives

1 1 ) - + x1 + (2ky*F)2 ν 2 2

(5)

Other models that do not account for the damping effect of the wall can be used. However, the model given by Alhusseini and Chen5 better describes the general system and appears to predict uniformly all of the transfer coefficients of interest: absorption, heating, and evaporation. The heat-transfer coefficient for evaporation is given by5

h/E )

δ*1/3Pr y* 1 dy* 0 1 1 + Pr Prt ν

∫

1.73 < Pr < 46.6 (6)

()

It is important to point out that most studies cover a very narrow range of Prandtl numbers. The exception is Limberg’s work,9 which, however, has been found to overpredict the real data. In the present study, the

1928


Figure 4. Differential volume element of the lamella turned on its side.

Figure 3. Calculated value of the eddy diffusivity of momentum at the entrances and exits of three consecutive effects. The changes in the value of /ν are a function of the dry solids concentration, which increases with each consecutive effect.

Prandtl numbers are extremely high at some operating conditions (see Table 1). However, the calculations show that, in the present study, turbulent flow exists for Prandtl numbers up to 90. Thus, the use of eq 6 is justified. The average velocity, u j , in the direction of flow in the film can be obtained from knowledge of the volumetric flow rate (V) and cross-sectional area (A). That is

∫0 ∫0 u*xδg dy dx a

u j)

V ) A

δ

δa

(7)

where u*, the velocity distribution, can be calculated from5

u* )

∫0

y*

(1 - δ*ζ*) dζ* (1 + ν)

(8)

As a check on the validity of this model, Figure 3 shows the eddy diffusivity of momentum (/ν) at operating conditions typical for black liquor falling film evaporation. The eddy diffusivity of momentum decreases with each successive effect (evaporator) because the kinematic viscosity increases as the dry solids concentration increases. Alhusseini and Chen5 presented similar results but for a water system at 373 K and a Reynolds number of 2500. The maximum value of /ν in that simple system is approximately 1. For black liquor operations, this is similar to the conditions at the exit of the third effect (labeled 3 out in Figure 3). It is reasonable to conclude that the present results provide some confidence that the proposed model describes the system properly because, at the operating conditions of interest, the eddy diffusivity of momentum is significant. This implies that the assumption of turbulent flow is valid. (Note that a value of unity means that the eddy viscosity and kinematic viscosity are equal. Thus, the viscous contribution to the total shear stress is the same as the contribution of the turbulent stresses.)

energy, and momentum balances. A differential volume approach, shown in Figure 4, is used. Assumption 1. The black liquor feed supplied to the plate (lamella) is at its boiling temperature. In the real case, there are two possibilities: (i) The liquor temperature is above the boiling point, and thus, flashing occurs when the liquor enters the effect. (ii) The liquor temperature is below the boiling point. In this case, the liquor must first be heated to its boiling point before evaporation occurs. In general, the evaporation is carried out at the liquor’s boiling point (Tb), which is a function of the pressure and a boiling point rise (∆Tb). The additional temperature increase is a function of the dry solids content. For black liquor mixtures, Tb is given by

Tb ) Tsw(P) + ∆Tb(Xs)

Assumption 2. The wall temperature is uniform along the length of the plate. In practice, this assumption might not hold. Assumption 3. No nucleate boiling occurs. The assumption is reasonable because high temperature differences are not favored in black liquor evaporation. When nucleate boiling is present at high temperatures, it is hypothesized that the accompanying high mixing rate increases collisions between the salt molecules, thus permitting the formation of crystals.14 This is the primary reason for evaporator scaling, which decreases heat transfer significantly (see section 4.2.4).2,14,15 Assumption 4. The parameters that describe the fluid flow and the liquor are constant in the horizontal direction (see Figure 4). 3.1. Energy Balance. Energy conservation applied to a differential volume element (see Figure 4) gives

E|t+∆t - E|t ) (Q - Qw)∆t

(10)

where E is the energy, t is the time, Q is the rate at which energy enters the differential volume through the plate, and Qw is the energy removed by the evaporation of water. Energy removal due to evaporation is given by

Qw ) W∆Hv

(11)

where W, defined as

3. Phenomenological Modeling This section describes the development of a model that is based on the fundamental principles of mass,

(9)

W≡

a∆zK(Tw - T) ∆Hv

(12)

Ind. Eng. Chem. Res., Vol. 42, No. 9, 2003 1929

is the mass of water evaporated from a surface of size a∆z and K, defined as

K≡

1 δw 1 + / λw h

(13)

E

is the heat-transfer coefficient, which is a function of the black liquor concentration, temperature, and mass flow rate. The amount of energy entering the differential volume is given by

Q ) a∆zK(Tw - T)

(14)

where Tw is the wall temperature on the steam side. Dividing eq 10 by ∆t and taking the limit as ∆t f 0 gives

dE ) Q - Qw dt

3.2.2. Balance on Dry Solids. At a given position z, eq 16 gives

∂(MXs) ∂Xs ∂Xs ∂G ∂M )M + Xs ) Gτ + Xsτ ∂t ∂t ∂t ∂t ∂t

Applying this equation from z to z + ∆z and using eq 18, the left-hand side of the eq 21 becomes

G ∂Xs Xs ∂G (GXs)|z - (GXs)|z+∆z + ) vz ∂t vz ∂t ∆z Taking the limit as ∆z f 0 results in the following partial differential equation

(

(15)

M|t+∆t - M|t ) G|z∆t - G|z+∆z∆t - W∆t

∂Xs XsvzW′ ∂Xs + vz ) ∂t ∂z G

(16)

3.2.1. Balance on Black Liquor. In the case of the black liquor component, eq 16 becomes

∂G ) G|z - G|z+∆z - W ∂t

(17)

where τ, the residence time in the differential volume, is assumed to be constant and can be expressed as a function of the velocity, vz, in the z direction. That is

τ)

∆z vz

(18)

(Note, vz ≡ u j .) Substituting this expression into eq 17 and defining W′ ≡ W/∆z gives

1 ∂G G|z - G|z+∆z ) - W′ vz ∂t ∆z

(19)

Taking the limit as ∆z f 0 then yields

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

(20)

(23)

To summarize, the balances are

∂G ∂G + vz ) -W′vz ∂t ∂z ∂Xs ∂Xs XsvzW′ + vz ) ∂t ∂z G This system is classified as quasilinear.16 The known initial conditions are

z)0

where M is the mass of the black liquor, G is the mass flow rate, and the other symbols are as before. Dividing through by ∆t and taking the limit as ∆t f 0 yields

τ

(22)

which, upon substitution of eq 20, gives

In practice, there is a small accumulation of energy, due to the boiling point rise. However, at the operating limits for a single evaporator, this change can be neglected. 3.2. Material Balance. Applying the principle of material conservation to a differential volume of size ∆z (see Figure 4) gives

∂M ) G|z - G|z+∆z - W ∂t

)

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

It is assumed that there is no accumulation of energy and that only a phase change occurs. Thus

Q ) Qw

(21)

{

G)G ∀ t g 0 X ) X0 s s,0

3.3. Numerical Scaling. Equations 20 and 23 can be made dimensionless by introducing the following variables

G G* ) G 0 Xs X* ) Xs,0 z z* ) L vz / vz ) vf,0 tvf,0 t* ) L W′ W′ ) W′ 0 LW′0 ) η G0

mass flow rate dry solids content space coordinate along the plate velocity along the plate time

(24)

Substitution into the mass balance equations above gives the following dimensionless system of PDEs

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

(25)

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

(26)

(

)

1930


Figure 5. Responses to a +5% change in feed mass flow rate: (top) response in time at z ) L, (bottom) response in space.

3.4. Black Liquor Properties. This section introduces the correlations used to calculate the model parameters. It should be noted that, although the correlations or modifications used here fit the process and the operating data, in every case, some corrections must be made as there are differences, some greater, some lesser, among different black liquors. 1. Density (G). The widely accepted correlation provided by Hough2 was used. It is given by

F ) 1007 - 0.495T + 6Xs

0 e Xs e 0.55

(27)

where T is the temperature and Xs is the dry solids content. 2. Kinematic Viscosity (ν). The correlation provided by Gullichsen and Paulapuro,1 using the coefficients for

softwood, was employed to calculate the kinematic viscosity

ln(ν) ) A(X) +

B(Xs)

T3 0 e Xs e 0.76 313 K e T e 388 K (28)

Here, A and B are third-order polynomial functions of the dry solids content. 3. Thermal Conductivity (λ). The experimental data presented in Ramamurthy et al.17 provide a reliable source for values of the thermal conductivity. These data are valid for 0 e Xs e 0.83 and 297 K e T e 367 K. 4. Specific Heat (Cp). The correlation for specific heat was obtained from Hough (originally Harvin and


Figure 6. Responses to a -5% change in feed mass flow rate: (top) response in time at z ) L, (bottom) response in space.

Brown)2

Cp ) 4106.2 + 3.822T - 3328.67Xs

(29)

It is valid in the range 0 e Xs e 0.6 and 297 K e T e 380 K. 5. Surface Tension (σ). The experimental data of Krishnagopalan et al.18 were employed. The data were obtained using a maximum bubble pressure technique, which ensures that certain specific requirements are satisfied for the proper testing of high-solids black liquor mixtures. The data are valid in the interval 0 e Xs e 0.65 and 298 K e T e 373 K. 6. Specific Heat of Evaporation (∆Hv). The specific heat of evaporation, as a function of pressure and temperature, is assumed to be the same as that of water.

The data for the specific heat were obtained from steam tables provided in Pavlov et al.19 7. Boiling Point Rise (∆Tb). The boiling point rise correlation was obtained from Hough (originally Hultin)2

∆Tb ) K

(

Xs 1 - Xs

)

6.5 < K < 7.5, 0 e Xs e 0.5 (30)

4. Results In this section, the simulated results of the phenomenological model when disturbances are introduced are presented and discussed. The method of orthogonal collocation on finite elements is used to solve numerically the system of PDEs given by eqs 25 and 26. As

1932


Figure 7. Responses to a +5% change in feed dry solids content: (top) response in time at z ) L, (bottom) response in space.

stated in ref 20, when the solution is expected to exhibit steep gradients, it is better to use trial functions, defined over only a part of the region. The orthogonal collocation method is also known to provide satisfactory solutions when the system to be solved is stiff (cf. finite difference methods) and to be less computationally expensive than the Galerkin method. The numerical solution to eqs 25 and 26 was obtained using two internal collocation points in each of 10 finite elements. This resulted in a system of 62 differentialalgebraic equations. There was no attempt to optimize the number of finite elements or the number of internal collocation points. 4.1. Simulated Dynamic Responses. The performance of the dynamic model is demonstrated here by applying step changes in the feed mass flow rate, the feed dry solids content, and the wall temperature. These

are some of the most important disturbances in this process. The impact of scaling (by steam) is also examined, as scaling is known to affect the heat-transfer process, thereby impacting the operation of the evaporator plant. The nominal conditions are a dry solids content of 0.3 kg/kg (kilograms of dry solids per kilogram of black liquor) and a mass flow rate of 2.2 kg/s. The temperature of the wall was chosen to be 370° K at these operating conditions. The black liquor mass flow rate and dry solids content are modified by (5% from the initial conditions, whereas the wall temperature is changed by (1% of its nominal value. To mimic scaling conditions, scale deposit thicknesses of 0.1 and 1 mm are used. The carbonate scale thermal conductivity is between 0.163 and 3.49 W/m‚K.19 An average value of 2 W/m‚K is used in this study. The


Figure 8. Responses to a -5% change in feed dry solids content: (top) response in time at z ) L, (bottom) response in space.

heat-transfer coefficient in this case is calculated by

K≡

1 δw δs 1 + + λw λs h/

(31)

E

In actual situations, the thermal conductivity in the presence of scaling can be determined experimentally, and thus, more precise calculations can be performed. 4.2. Analysis of the Simulated Results. 4.2.1. Effect of Change in Feed Mass Flow Rate. The temporal and spatial behaviors of the model when subjected to disturbances in the mass flow rate are shown in Figures 5 and 6. The time profiles are calculated at z ) L, the end of the plate. The results show that a (5% change in the mass flow rate produces

almost the same change in the output. The system is stable, and the mass flow rate attains a new steady state. In contrast, the steady-state value of the dry solids content remains unchanged. This is because the heattransfer coefficient changes in the same direction as the mass flow rate. Thus, the system is self-compensating for small changes in the mass flow rate. The temporal behavior shows a moving front that develops after the disturbance is introduced. The traveling front can also be observed in the spatial direction. The residence time of the black liquor on the heating surface is calculated to be around 7 s. Compared to the results of Adams,14 where the residence time for the evaporator studied was approximately 10 s, this result appears to be consistent, considering that, in his work, the dry solids content was higher.

1934


Figure 9. Responses to a +1% change in wall temperature: (top) response in time at z ) L, (bottom) response in space.

4.2.2. Effect of a Change in the Feed Dry Solids Content. Figures 7 and 8 show the simulated response for (5% change in the dry solids contents of the feed. The behavior is similar to that observed for mass flow rate disturbances. The process gain is on the order of unity, and there is no change in the output mass flow rate. This means that the system is not sensitive to small changes in the dry solid content. 4.2.3. Effect of a Change in the Wall Temperature. The change in the wall temperature is introduced uniformly over the entire plate. The response is shown in Figures 9 and 10. There is no transport delay in the response because the uniform change affects the end of the plate immediately. The time required to reach a new steady state is the same as in the previous two cases. The process gain with respect to the mass flow rate and the dry solids content is about unity.

4.2.4. Effect of Scaling. The results shown in Figure 11 indicate the effect of scaling on the operation of the evaporator. Scale with an assumed thickness of 0.1 mm is shown to reduce the output dry solids content by 10%, whereas a scale thickness of 1 mm reduces the output dry solids content by more that 50%. These results show that the condition of scaling has to be avoided. From a control perspective, the temperature of the wall has to be controlled to prevent scaling. 4.2.5. Meaning of η. A mathematical analysis of the dimensionless parameter η, defined in eq 24, shows that, when G0 f 0, η f ∞, and when G0 f ∞, η f 0. In general, for very high feed flow rates, there will be no evaporation, whereas for very low feed flow rates, the evaporation will be almost instantaneous. This is consistent with the physics of the process. Thus, W′0 is a weak function of G0.


Figure 10. Responses to a -1% change in wall temperature: (top) response in time at z ) L, (bottom) response in space.

Numerical results for a feed flow rate of approximately 100 kg/s (assuming that feeding this amount of liquor to a single plate is possible) showed negligible evaporation. For the opposite case, that of very low flow rates, the flow regime is not expected to be turbulent. Thus, the model should not be used to infer any behavior unless the flow regime is turbulent. In general, reasonable values of η, at reasonable operating conditions, can provide a measure of evaporation intensity. 4.2.6. Significance of the Results from a Control Perspective. Observe that the proposed model shows the system to be open-loop stable. Both the mass flow rate and the dry solids content reach their respective steady-state values after step changes in the positive and negative directions. Another important observation is that, in the specific range of operating conditions

investigated, the system is almost decoupled with respect to the mass flow rate and the dry solids content. That is, changes in the mass flow rate do not lead to significant changes in the dry solids content and vice versa. Finally, the process gain is almost unity in the case of the mass flow rate and the dry solids content. 5. Discussion It is prudent to point out some limitations of the proposed model. (i) The model represents only one type of evaporator, the falling film evaporator. It is possible that, with some minor modifications, a similar model could be used for the tubular type of falling film evaporators, because the tube curvature is negligible with respect to the film

1936


Figure 11. Spatial profiles of the mass flow rate and the dry solids content at steady state as a function of scale thickness.

thickness. Thus, the tube surface can be assumed to be locally flat, and the same geometry can be used. (ii) Because the model does not account for the direction of vapor flow, it can be used for evaporators in which the vapor traffic proceeds down the lamella. However, for other types of evaporators, especially the rising film type, it is suspected that a completely different approach will be needed because the processes are of a very different nature. (iii) The energy supplied was assumed to be uniform along the lamella. In practice, this state of uniform heating might not be realistic. (iv) The hydrodynamics of the flow reflects turbulent conditions. This assumption can be argued to be valid in the early evaporators, but in later ones and in the superconcentrators, the flow characteristics (Re, Pr, and Ka numbers) are not turbulent. Consequently, it is possible for the flow to be laminar or wavy laminar in these regions.6,9,11,12,21,22 Thus, the heat-transfer coefficient and other relevant flow-dependent parameters have to be determined using the proper flow models. (v) The correlations for black liquor physical properties also introduce uncertainties that affect the predictive nature of the model. With a black liquor dry solids content outside the range of 0.3-0.4, the empirical relationships must be reexamined. As mentioned previously, the black liquor will have very irregular properties from mill to mill, as it is a function of many things, including wood type, cooking conditions, etc. Thus, for every separate pulp mill, it is advised to carry out experiments to determine the exact correlation coefficients. The next steps in this study are to (1) account for other processes that impact the operation, including flashing or heating of the black liquors, mixing that occurs due to recirculation, and accumulation at the bottom of the effect; (2) assess the impact of nonuniform heating; (3) address the transitions among the different flow regime types; and (4) model the entire evaporator plant, addressing the numerical issues that come with higher dimensions and greater complexities. A few remaining remarks about the control of the entire evaporator plant are appropriate. In general, the control of a multiple-effect evaporator plant is not a

trivial task; for that matter, even the control of a single effect is not without serious challenges. There are measurement problems, especially with the dry solids content and a dearth of quantity and type of manipulated variables, and if both the mass flow rate and dry solids content are the controlled variables, the obvious choice of manipulated variables creates interactions. Finally, these systems are constrained in terms of both inputs and outputs. Thus, the control structure, whether simple or complex, must address these issues. In practice, the control of truly nonlinear distributedparameter systems (DPSs) is usually accomplished by disregarding the spatial nature of the system and applying lumped system methods. However, by so doing, important interactions due to the underlying distributed characteristics can be ignored. Zheng and co-workers proposed a novel method for developing a reduced-order model for implementable control solutions.23 In another approach to the synthesis of control structures, Cardoso24 applied robust control concepts that ignored any spatial interactions. We will investigate both approaches in the future. 6. Summary This study involved the development of a phenomenological model that describes the transport processes for a single plate in a falling film black liquor evaporator. To be precise, the plate is at the front end of the evaporator plant. The fluid flow regime was assumed to be turbulent, and material and energy balances were used to describe the physicochemical processes. The parameters for the model were obtained from wellknown correlations and validated by real operating data, when available. In cases where the correlations were inaccurate or not appropriate, correlations were developed using the data. The resulting nonlinear model (system of PDEs) was solved using the well-known method of orthogonal collocation on finite elements. The simulated responses at different loads, feed mass flow rate and feed solid dry content changes, wall temperature changes, and heat-transfer limitations due to the presence of scaling, were shown to yield responses that can be explained by process knowledge. Acknowledgment The first author was supported by special research assistantship funds provided by the College of Engineering at Texas Tech University, Lubbock, TX, and Petroleum Research Foundation Grant ACS-PRF 35789AC9. Both authors are grateful to Tembec (Quebec, Canada), Kraft Pulp Division Skookumchuck Operations personnel, and the Crestbrook (Cranbrook, BC, Canada) plant for operating data. Nomenclature A ) film cross-sectional area (m2) A/i ) empirical function for correction of inadequacies in the modeling of the wavy interface a ) plate length (m) Cp ) specific heat (J/kg) Dw ) van Driest’s damping function E ) energy (J) F ) mixing length correction function G ) black liquor mass flow rate (kg/s) G0 ) initial black liquor mass flow rate (kg/s)

Ind. Eng. Chem. Res., Vol. 42, No. 9, 2003 1937 G* ) dimensionless black liquor mass flow rate, G/G0 g ) gravitational acceleration, 9.81 (m/s2) h ) film heat-transfer coefficient (W/m2‚K) hE ) film heat-transfer coefficient for evaporation (W/m2‚K) / hE ) dimensionless film heat-transfer coefficient for evaporation, hE(ν2/g)1/3/λ K ) overall heat-transfer coefficient (W/m2‚K) K ) constant in ∆Tb equation k ) von Karman’s constant ) 0.41 Ka ) Kapitza number, µ4g/Fσ3 L ) plate length (m) M ) mass (kg) m ) exponent of Re in eddy diffusivity of the momentum equation for the outer boundary layer P ) pressure (Pa) Pr ) Prandtl number, Cpµ/λ Pri ) constant in Prt function, 0.86 Prt ) turbulent Prandtl number Q ) heat flux (J/s) Re ) film Reynolds number, 4Γ/µ Ret,t ) transition Reynolds number for turbulent regime T ) black liquor temperature (K) Tb ) black liquor boiling temperature (K) Tsw ) water saturation temperature (K) Tw ) wall temperature at steam side (K) t ) time (s) t* ) dimensionless time, tvf,0/L u* ) velocity distribution (m/s) u j ) average velocity (m/s) vf ) friction velocity, xδg vf,0 ) friction velocity at initial conditions, xδ0g V ) black liquor volumetric flow rate (m3/s) vz ) flow velocity along the plate (m/s) v/z ) dimensionless velocity along the plate, vz/vf,0 W ) evaporated water mass flow rate (kg/s) Xs ) black liquor dry solids content (kg/kg) Xs,0 ) initial black liquor dry solids content (kg/kg) X* ) dimensionless dry solids content y ) space coordinate perpendicular to the free surface (m) y* ) dimensionless space coordinate perpendicular to the free surface, yvf/ν y/i ) y* at the interface between the turbulent core and outer boundary layer z ) space coordinate along the plate (m) z* ) dimensionless space coordinate along the plate, z/L ∆Hv ) specific heat of evaporation (J/kg) ∆Tb ) boiling point rise (K) Greek Letters R0 ) smallest eigenvalue in Limberg’s model Rt ) thermal diffusivity (m2/s) δ ) film thickness (m) δ0 ) film thickness at initial conditions (m) δw ) plate wall thickness (m) δs ) scale thickness (m) δ* ) dimensionless film thickness,δuf/ν η ) dimensionless variable ) eddy viscosity (m2/s) λ ) film thermal conductivity (W/m‚K) λw ) wall thermal conductivity (W/m‚K) λs ) scale thermal conductivity (W/m‚K) µ ) dynamic viscosity (Pa‚s) ν ) black liquor kinematic viscosity (m2/s) F ) black liquor density (kg/m3) σ ) surface tension (N/m)

τ ) residence time (s) Γ ) black liquor mass flow rate per unit length (kg/s ‚ m)

Literature Cited (1) Gullichsen, J., Paulapuro, H., Eds. Chemical Pulping, 1st ed.; TAPPI Press: Atlanta, GA, 1999. (2) Hough, G. Chemical Recovery in the Alkaline Pulping Processes, 1st ed.; TAPPI Press: Atlanta, GA, 1985. (3) Technical Support Document for Best Management Practices for Spent Liquor Management, Spill Prevention and Control; Technical Report EPA-821-R-97-011; U.S. Environmental Protection Agency, U.S. Government Printing Office: Washington, DC, 1997. (4) Kraft Pulp Mill Compliance Assessment Guide (CAA, CWA, RCRA and EPCRA); Technical Report EPA/310-B-99-001; U.S. Environmental Protection Agency, U.S. Government Printing Office: Washington, DC, 1999. (5) Alhusseini, A. A.; Chen, J. C. Transport phenomena in turbulent falling films. Ind. Eng. Chem. Res. 2000, 39, 2091. (6) Chun, K. R.; Seban, R. A. Heat transfer to evaporating liquid films. J. Heat Transfer 1971, 93, 391. (7) Åsblad, A.; Berntsson, T. Surface evaporation of turbulent falling films. Int. J. Heat Mass Transfer 1991, 34, 835. (8) Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. Transport Phenomena, 2nd ed.; John Wiley & Sons: New York, 2000. (9) Limberg, H. Wa¨rmeu¨bergang an turbulente und laminare rieselfilme. Int. J. Heat Mass Transfer 1973, 16, 1691. (10) Miladinova, S.; Slavtchev, S.; Lebon, G.; Legros, J. Long wave instabilities of nonuniformly heated falling films J. Fluid Mech. 2002, 453, 153. (11) Miyara, A. Numerical analysis on flow dynamics and heat transfer of falling liquid films with interfacial waves. Heat Mass Transfer 1999, 35, 298. (12) Chang, H.; Demekhi, E. A.; Kaladin, E. Simulation of noise driven wave dynamics on a falling film. AIChE J. 1996, 42, 1553. (13) van Driest, E. R. On turbulent flow near a wall J. Aeronaut. Sci. 1956, 23, 1007. (14) Adams, T. N. Sodium salt scaling in black liquor evaporators and concentrators TAPPI J. 2001, 84, 1. (15) Severtson, S. J.; Duggirala, P. Y.; Carter, P. W.; Reed, P. E. Mechanism and chemical control of calcium carbonate scaling in the kraft process. TAPPI J. 1999, 82, 167. (16) Rhee, H. K.; Rutherford, A.; Amundson, N. R. First-Order Partial Differential Equations, 1st ed.; Prentice Hall: Englewood Cliffs, NJ, 1989; Vol. II. (17) Ramamurthy, P.; van Heiningen, A. R. P.; Kubes, G. J. Viscosity and thermal conductivity of black liquor. TAPPI J. 1993, 76, 175. (18) Krishnagopalan, J.; Stockel, I. H.; Pendse, H. P.; Kiran, E. Surface tension of kraft black liquor Chem. Eng. Technol. For. Prod. Process. 1988, 2, 97. (19) Pavlov, K. F.; Romankov, P. G.; Noskov, A. A. Examples and Problems on Unit Operations in the Chemical Technology, 10th ed.; Himia: Leningrad, USSR, 1987. (20) Finlayson, B. A. Nonlinear Analysis in Chemical Engineering; McGraw-Hill: New York, 1980. (21) Levich, V. G. Physicochemical Hydrodynamics; Prentice Hall: Englewood Cliffs, NJ, 1962. (22) Chun, K. R.; Seban, R. A. Performance prediction of fallingfilm evaporators. J. Heat Transfer 1972, 94, 432. (23) Hoo, K. A.; Zheng, D. Low-order control-relevant models for a class of distributed parameter systems. Chem. Eng. Sci. 2001, 56, 6683. (24) Zheng, D.; Hoo, K. A.; Piovoso, M. J. Low-order model identification of distributed parameter systems by a combination of singular value decomposition and the Karhunen-Loe`ve expansion. Ind. Eng. Chem. Res. 2002, 41, 801. (25) Cardoso, A. J. L. Modelizaç aö e Controlo Robusto de uma Cadeia de Evaporaç aö. M.S. Thesis, Universidade de Coimbra, Coimbra, Portugal, 1995.

Received for review June 28, 2002 Revised manuscript received October 18, 2002 Accepted February 19, 2003 IE020483A

A Distributed-Parameter Model of Black Liquor Falling Film

Recommend Documents