A heat transfer model for the steady-state simulation of climbing-falling-film plate evaporators

Journal of Food Engineering 54 (2002) 309–320 www.elsevier.com/locate/jfoodeng

A heat transfer model for the steady-state simulation of climbing-falling-film plate evaporators ~o Andrade C.P. Ribeiro Jr., M.H. Can

*

Departamento de Engenharia Quımica, Universidade Federal de Minas Gerais, Rua Espırito Santo, 35, Centro, 30160-030 Belo Horizonte, MG, Brazil Received 2 March 2001; accepted 3 November 2001

Abstract As an economical alternative for the concentration of heat sensitive materials, the climbing-falling-film plate evaporator was used for the first time in 1958. Nevertheless, there seems to be little information in the literature about such equipment. In this paper, an algorithm for the steady-state simulation of this evaporator is presented, which takes into account a general unit with n sets of four plates, whose feed can be at a temperature equal, inferior or superior to its saturation value for the operating pressure. The heating fluid, steam, can be fed either in a saturated or superheated state, with the possibility of condensate subcooling within the equipment channels. Initially, the algorithm was tested for simulating a sample unit at different operating conditions, and the physical consistency of the obtained results was analysed. The efficiency of the developed model was then demonstrated by the successful simulation of six industrial units used for concentrating milk in a Brazilian milk powder plant. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Plate evaporator; Algorithm; Simulation; Milk concentration

1. Introduction Evaporation is a unit operation in which a solution is concentrated by removing part of the solvent in the form of vapour. In 99% of the industrial cases, the solvent is water and the latent heat of evaporation is supplied by condensing steam, whose energy is transmitted to the solution by indirect heat transfer through metallic surfaces (Standiford, 1963). Although evaporators with tubular heating surfaces are by far the most common choice in the chemical industry, the problems of chemical degradation which arise in the concentration of heat-sensitive materials, such as milk, fruit juices, coffee, tea and whey, have led to the development of a new type of equipment known as plate evaporator, which has actually been used to concentrate liquid foods for at least 40 years (Carter & Kraybill, 1966; Kumar, 1993). In a plate evaporator, instead of tubes, metal plates are employed to separate the process and heating fluids. Such units require low headroom, present low installation costs (Usher, 1965) and are flexible, since plates can *

Corresponding author. Fax: +55-31-32381789. E-mail address: [email protected] (M.H. Ca~ no Andrade).

be added or removed for changes in the processing capacity. Besides, one can easily inspect and clean all heattransfer surfaces. The plates arrangement for a climbing-falling-film plate evaporator is shown in Fig. 1. Plates are grouped in sets of four with a climbing and a falling-film product section, between which there is a steam passage. The specific number of plate sets is mounted in a frame to give the necessary heat-transfer surface. All sets in the frame are fed in parallel with the material to be evaporated, and the concentrated product, together with the vapour produced, leaves from the head of the frame for separation in a centrifugal separator. A more detailed description of the equipment is given by Usher (1965). Despite their introduction in 1958, there appears to be little information in the literature about the concentration of solutions in plate evaporators. Only a few papers on this issue have been published, which are rather general and not much technical in nature (Brotherton, 1994; Gray, 1971; Hellstr€ om, 1994; Kumar, 1993; Licha, Valentin, Wersel, & Witte, 1989; Morgenroth, Jayatilaka, & Punter, 1997; Usher, 1965). To the authors’ knowledge, neither the mathematical modelling nor the steady-state simulation of climbing-falling-film plate evaporators has ever been addressed before.

0260-8774/02/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 2 6 0 - 8 7 7 4 ( 0 1 ) 0 0 2 1 7 - 5

310

C.P. Ribeiro Jr., M.H. Ca~no Andrade / Journal of Food Engineering 54 (2002) 309–320

Nomenclature heat-transfer area, m2 plates heat-transfer area per length unit, m distance between two consecutive plates, m
1 boiling number, qG
1 DHcon , dimensionless boiling point rise, K constants in Eq. (41), dimensionless 0:8 convection number, ð1
qualÞ
0:5 0:5 qv ðqual ql Þ , dimensionless Cp heat capacity, J kg
1 K
1 Deq channel equivalent diameter, m F total mass flow rate fed to the evaporator, kg s
1 Ffl fluid-dependent parameter, dimensionless Frl Froude number with all flow as fluid,
1 G2 ðq2l gDeq Þ , dimensionless G mass flux per unit area in the channel, kg m
2 s
1 g acceleration due to gravity, m s
2 H stream specific enthalpy, J kg
1 DHcon enthalpy of condensation, J kg
1 DHsuph enthalpy change associated with the cooling of superheated steam to its saturation temperature, J kg
1 h heat-transfer coefficient, W m
2 K
1 k plate thermal conductivity, W m
1 K
1 L plate length, m l plate thickness, m LMTD logarithmic-mean temperature difference, K m stream mass flow rate in the channels, kg s
1 n number of four-plate sets in the evaporator, dimensionless Nu Nusselt number, hDeq j
1 , dimensionless P pressure, Pa Pr Prandtl number, Cp lj
1 , dimensionless q heat flux, W m
2 Q heat exchanged between fluids, W qual steam quality, dimensionless Re Reynolds number, GDeq l
1 , dimensionless Rel liquid phase Reynolds number, Gð1
qualout Þ Deq l
1 , dimensionless T temperature, K solvent saturation temperature, K Ts A a b Bo BPR c1 –c5 Co

In this paper, an algorithm for the steady-state simulation of such evaporators is presented. The algorithm is based on the simultaneous resolution of the global heat and mass balance equations, together with the heattransfer rate equations, for each channel of the fourplate sets.

U w X

overall heat-transfer coefficient, W m
2 K
1 channel width, m total solids mass fraction, dimensionless

Subscripts 1–4 climb conc cond ev evc evf evfl fall feed flash h i in l oper out s sol SP stm TP v w

regions 1–4 of the evaporating solution path climbing-film region concentrated solution condensate evaporated solvent evaporated solvent in the climb-film region evaporated solvent in the falling-film region evaporated solvent during the flash process falling-film region feed stream flash process of the superheated feed stream heating of the subcooled solution to its saturation temperature each region of the evaporating solution path inlet condition liquid phase operating value outlet condition saturated state solution single-phase condensing steam two-phase vapour phase steam side of the plates

Superscripts ev sol stm w

evaporated solvent solution condensing steam steam side of the plates

Greek symbols a; b q l j c

constants in Eq. (40), dimensionless density, kg m
3 viscosity, Pa s fluid thermal conductivity, W m
1 K
1 parameter defined by (37), dimensionless

2. Mathematical modelling Fig. 2 shows the important variables in the four-plate set of a climbing-falling-film plate evaporator. In order to elaborate the mathematical model, the following assumptions are made:

C.P. Ribeiro Jr., M.H. Ca~no Andrade / Journal of Food Engineering 54 (2002) 309–320

311

evaporator. The overall mass balances for the climbing and falling-film sections are written as follows: mfeed ¼ mclimb þ mevc ;

ð1Þ

mclimb ¼ mfall þ mevf ;

ð2Þ

where mfeed is the feed flow rate in the climbing-film channel; mclimb and mfall are the solution mass flow rates at the outlet of each section; mevc and mevf are the evaporated solvent mass flow rates for the climbing and falling-film sections, respectively. Considering assumption (g), one can write the component mass balances below: Fig. 1. A climbing-falling-film plate evaporator comprising three fourplate sets.

mfeed Xfeed ¼ mclimb Xclimb ;

ð3Þ

mclimb Xclimb ¼ mfall Xfall ;

ð4Þ

where Xfeed is the feed total solids concentration; Xclimb and Xfall are, respectively, the total solids contents of the solutions at the outlet of the climbing and falling-film sections. As for the heat balances, assumptions (c), (d), (e), (h) and (j) lead to:

Fig. 2. Schematic representation of a four-plate set.

(a) steady-state operation; (b) uniform flow distribution in the channels; (c) each fluid is ideally mixed in the direction perpendicular to the flow; (d) heat conductivity in the direction of flow is negligible; (e) channels are of equal length and heat-transfer area; (f) constant cross-section of each channel on its whole length; (g) no solution entrainment in the vapour; (h) no chemical reactions taking place in the channels; (i) no pressure drop in the channels; (j) adiabatic operation. On account of assumption (b), the modelling of a four-plate set corresponds to the modelling of the whole

mfeed Hfeed
ðmclimb Hclimb þ mevc Hevc Þ þ2AUclimb LMTDclimb ¼ 0;

ð5Þ

mclimb Hclimb
ðmfall Hfall þ mevf Hevf Þ þ2AUfall LMTDfall ¼ 0;

ð6Þ

where Hfeed is the specific enthalpy of the feed; Hclimb and Hfall are the specific enthalpies of the solutions at the outlet of the climbing and falling-film sections, respectively; Hevc and Hevf are the specific enthalpies of the evaporated solvent in the climbing and falling-film sections, respectively; Uclimb and Ufall are the overall heat-transfer coefficients associated with each of these sections; A is the plates heat-transfer area and LMTD is the logarithmic-mean temperature difference, defined as follows, for either saturated or superheated steam LMTDclimb ¼ LMTDfall ¼

ln




ðTsstm

Tclimb
Tfeed ;
Tfeed Þ=ðTsstm
Tclimb Þ

Tfall
Tclimb
 ln ðTsstm
Tclimb Þ=ðTsstm
Tfall Þ

ð7Þ ð8Þ

being Tfeed the feed inlet temperature; Tclimb and Tfall the solution temperatures at the outlet of the climbing and falling-film sections, respectively, and Tsstm the steam saturation temperature. The analysis of Eqs. (7) and (8) shows that the thermal potential derived from the use of steam as heating fluid is always calculated at the saturation temperature, even when superheated steam is employed. This is due to the fact that, once superheated steam condenses, a condensate film at Tsstm is formed on the plates surface and acts as an additional resistance to heat transfer (Kreith, 1973; Westphalen, 1999).

312

C.P. Ribeiro Jr., M.H. Ca~no Andrade / Journal of Food Engineering 54 (2002) 309–320

For the steam channels, supposing there is no condensate subcooling and considering assumptions (c), (i) and (j), the heat balance is equal to   mstm DHsuph þ ðqualin
qualout ÞDHcon þ AðUclimb LMTDclimb þ Ufall LMTDfall Þ ¼ 0:

ð9Þ

In the previous equation, mstm is the steam mass flow rate in the channels; qual is the steam quality, that is, the mass fraction of mstm which is actually vapour (Smith, Van Ness, & Abbott, 1996); DHcon is the heat of condensation and DHsuph is the enthalpy change of cooling superheated steam, at pressure Pstm , from Tstm to its saturation temperature, Tsstm
 DHsuph ¼ Hstm Tsstm ; Pstm
Hstm ðTstm ; Pstm Þ: ð10Þ Finally, from assumption (i), the solution outlet temperature for each section can be written as a function of the operating pressure, Poper , and total solids content, using the concept of boiling point rise (Billet, 1988):
 Tclimb ¼ Ts Poper þ BPRðXclimb Þ; ð11Þ
 ð12Þ Tfall ¼ Ts Poper þ BPRðXfall Þ; where Ts is the solvent saturation temperature at Poper and BPR is the boiling point rise, which is a function of the total solids content. It should be emphasised that, if BPR is so small that its effect may be neglected, the solution temperature will only depend upon the pressure in the channels. From assumption (i), it follows that, in this case, the solution temperature will be constant throughout the channels, so that LMTD can be replaced by the simple difference Tsstm
Ts ðPoper Þ in Eqs. (5), (6) and (9). An analysis of Fig. 2 shows that there are 18 independent variables for only nine equations and, therefore, nine design parameters must be specified for a unique solution to be obtained. The same set of equations can be used either for designing a new evaporator or for predicting the performance of an existing one, depending upon the specified parameters. In the former case, the desired outlet solids concentration ðXfall Þ is given and the necessary heattransfer area is calculated, whereas the opposite procedure is adopted for performance prediction problems. Whichever the case, both the feed properties and the operating pressure are always known.

been recommended for solving the equations of multiple effect evaporation systems (Holland, 1975), was chosen to solve all equations simultaneously. The developed algorithm, implemented in FORTRAN language, enables one to predict the performance of a generic climbing-falling-film plate evaporator comprising n four-plate sets, to which solution may be pumped subcooled, saturated or superheated in relation to the equipment’s operating pressure. Superheated or saturated steam is accepted as heating fluid, with the possibility of condensate subcooling. The flowchart for the developed algorithm is presented in Fig. 3. As input data, the feed streams properties (flow rate, temperature, pressure and total solids mass fraction), the operating pressure and the plates number and dimensions are required. The output data include total solids content, evaporated flow rate, solution temperature and overall heat-transfer coefficients for both the climbing and falling-film sections, together with the condensate outlet temperature and, when this is not equal to the saturation value, the fraction of the equipment’s heat-transfer area which is used for cooling the condensate. Once the input data are read, the streams flow rates in the channels which compose the evaporation units are quantified. From Fig. 1 and assumption (b) mstm ¼

Fstm ; 2n þ 1

ð13Þ

mfeed ¼

Ffeed ; n

ð14Þ

where Fstm and Ffeed are, respectively, the total flow rates of steam and solution pumped to the evaporator. Before commencing the calculation itself, the steam saturation condition is checked, being the algorithm aborted whenever the values provided for Tstm and Pstm do not have thermodynamic consistency. As shown in Fig. 3, the calculation sequence, as well as the equations actually solved, depends upon the feed temperature. The main features of the calculation procedure are highlighted in the next sections. Regardless of the feed condition, as a convergence criterion for all non-linear equation systems, the residues for all functions must be less than 10
6 . Furthermore, the computed variables are always kept within the range of physical significance values.

3. Algorithm description

3.1. Saturated feed stream

The general set of equations associated with an evaporator can be solved simultaneously by some generalised nonlinear numerical method or, alternatively, by linearising the equations and using integrated loops of iteration (Lambert, Joye, & Koko, 1987). In this work, the Newton–Raphson method, which has already

For a saturated solution, initial values of Xclimb and Xfall are generated considering total condensation of steam in the channels. Let L be the plates length and a their heat-transfer area per length unit; substituting Eqs. (1) and (3) into (5), for the climbing-film region, and Eqs. (2)–(4) into (6), for the falling-film one, the heat

C.P. Ribeiro Jr., M.H. Ca~no Andrade / Journal of Food Engineering 54 (2002) 309–320

313

Fig. 3. Flowchart of the developed algorithm.

balance equations are rewritten as follows, using a functional form for the enthalpies:  Xfeed mfeed Hsol ðTfeed ; Xfeed Þ
Hsol ðTclimb ; Xclimb Þ Xclimb   
 Xfeed
1
Hevc Tclimb ; Poper Xclimb ð15Þ þ 2LaUclimb LMTDclimb ¼ 0;  1 1 Hsol ðTclimb ; Xclimb Þ
Hsol ðTfall ; Xfall Þ mfeed Xfeed Xclimb Xfall   
 1 1

Hevf Tfall ; Poper Xclimb Xfall þ 2LaUfall LMTDfall ¼ 0:

ð16Þ

This new arrangement reduces the six original equations to only two, in which Xclimb and Xfall are the independent variables. Assuming that no liquid is entrained in the steam fed to the evaporator, Eq. (9) becomes

  mstm DHsuph þ ð1
qualout ÞDHcon þLaðUclimb LMTDclimb þ Ufall LMTDfall Þ ¼ 0:

ð17Þ

Furthermore, since all energy received by the solution in the channels is necessarily transferred through the plates by conduction, the following relations can be written (Choudhary, Ali, & Das, 1996):  
1 1 l Uclimb LMTDclimb
þ hclimb kclimb Tclimb
Tfeed  ¼ 0; 
w ð18Þ w ln ðTclimb
Tfeed Þ=ðTclimb
Tclimb Þ 


1 1 l þ hfall kfall Tfall
Tclimb  ¼ 0; 
w w ln ðTfall
Tclimb Þ=ðTfall
Tfall Þ

Ufall LMTDfall


ð19Þ

where l is the plate thickness; hclimb and hfall are the heattransfer coefficients on the solution side for the climbing and falling-film regions, respectively; kclimb and kfall are

314

C.P. Ribeiro Jr., M.H. Ca~no Andrade / Journal of Food Engineering 54 (2002) 309–320

the thermal conductivity of the plate in these regions, whose values are not necessarily the same, due to the influence of the total solids concentration on the soluw w tion saturation temperature; Tclimb and Tfall are the plate surface temperatures on the steam side for the climbing and falling-film sections, which, however close to Tsstm , can never be equal to it, or no heat transfer would take place. By solving Eqs. (15)–(19) simultaneously, the outlet steam quality, the total solids contents and the surface temperatures for each evaporation section are all computed. If the steam quality is superior or equal to zero, outlet stream flow rates are quantified and the obtained results are, then, printed (Fig. 3). Nonetheless, provided that the overall heat-transfer coefficients are sufficiently high, the amount of energy possible to be transferred to the solution will be superior to the available one, resulting in a negative value for the outlet steam quality (Eq. (17)), which means that a steam flow rate greater than the one specified could be completely condensed. In practice, as the steam flow rate is fixed, part of the equipment heat-transfer area is utilised for cooling the condensate. 3.2. Condensate subcooling If there is condensate subcooling, the energy balances presented previously in Section 2 do not apply, since such event was not considered during their formulation. Therefore, a new set of equations must be written. These new equations derive from a more complex model which takes condensate subcooling into account. As presented in Table 1, this new model assumes that the evaporating solution path can be divided into four regions, each of them having its own overall heattransfer coefficient ðfUi g; 1 6 i 6 4Þ. In this case, the variables to be calculated are the total solids content at the outlet of each region ðfXi g 1 6 i 6 4Þ, the plate length associated with each region ðfLi g; 1 6 i 6 4Þ and the condensate outlet temperature, Tcond , summing nine variables, which, hence, would require nine independent equations. However, since the same steam channel provides energy for both the climbing and falling-film regions, the following relations are valid: L1 ¼ L4 ;

ð20Þ

Table 1 Features of the four regions of the evaporating solution path for the model taking condensate subcooling into account Region

Description

1 2 3 4

Climbing film, condensate subcooling Climbing film, condensing steam Falling film, condensing steam Falling film, condensate subcooling

L2 ¼ L3 ;

ð21Þ

L1 þ L2 ¼ L3 þ L4 ¼ L:

ð22Þ

Consequently, only the length of one region is an independent variable from which the values associated with the other regions can be calculated, reducing to six the number of necessary equations. For each evaporation region, the following material and heat balances can be written: sol ev msol i
1 ¼ mi þ mi ;

ð23Þ

sol msol i
1 Xi
1 ¼ mi Xi ;

ð24Þ

msol i
1 Hsol ðTi
1 ; Xi
1 Þ

2Ui Li aLMTDi þ


ev ev
msol Ti ; Poper ¼ 0; i Hsol ðTi ; Xi Þ þ mi H

ð25Þ

where msol i
1 is the solution mass flow rate at the entrance of region i; msol i is the mass flow rate of the concentrated solution which leaves region i; mev i is the mass flow rate of the evaporated solvent for region i; Xi
1 and Xi are the total solids contents of the solutions at the entrance and outlet of region i, respectively; Ti
1 and Ti are the temperatures of the solutions which enter and leave region i, respectively. The index zero refers to the feed stream. Upon substituting (23) and (24) into (25), one obtains four independent equations similar to Eqs. (15) and (16). The two equations left are the heat balances for the heating fluids, that is, steam and condensate:
 mstm DHsuph þ DHcon þ ð L
L1 ÞaðU2 LMTD2 þ U3 LMTD3 Þ ¼ 0;


mstm Hcond ðTcond Þ
Hcond Tsstm þ L1 aðU1 LMTD1 þ U4 LMTD4 Þ ¼ 0

ð26Þ

ð27Þ

in which Hcond is the specific enthalpy of the condensate. 3.3. Superheated feed stream When the feed stream is superheated, it is assumed that part of it vaporises instantaneously or ‘flashes’, thereby reducing the feed temperature to its saturation value in relation to the evaporator operating pressure. The solution total solids concentration after the flash, Xflash , can be quantified by the equation below, obtained from the mass and heat balances  Xfeed Hsol ðTfeed ; Xfeed Þ
Hsol ðTflash ; Xflash Þ Xflash   
 Xfeed þ 1
Hevfl Tflash ; Poper ¼ 0 ð28Þ Xflash in which Tflash is the solution temperature after the instantaneous vaporisation, and Hevfl is the specific enthalpy of the solvent vapour formed on account of the flash.

C.P. Ribeiro Jr., M.H. Ca~no Andrade / Journal of Food Engineering 54 (2002) 309–320

The solution mass flow rate for this new concentration, mflash , is given by the component mass balance mflash

mfeed Xfeed ¼ : Xflash

ð29Þ

The result of this flash process is a saturated solution whose flow rate and total solids content are equal to mflash and Xflash , respectively, and, hence, once these parameters are computed, the calculation follows the steps previously described for the saturated solution, as shown in Fig. 3. 3.4. Subcooled feed stream If the evaporator feed is subcooled, part of the equipment heat-transfer area is used for heating the solution to its saturation temperature. In this case, apart from the variables connected to the simpler model for the saturated feed, the plate length required for the solution heating, Lh , must be determined, which leads to a total of six independent variables. Considering that the equipment purpose is to concentrate the feed stream, the use of more than half of its heat-transfer area for heating the solution would be of little practical interest. Therefore, it was assumed that all heating without phase transition takes place in the climbing-film region. As a result, if it is further assumed that the plate temperature on the steam side is constant, Eqs. (16), (18) and (19) apply directly. The heat balance equations for the saturated solution in the climbing-film region and condensing steam are rewritten as follows: 
 Xfeed mfeed Hsol Tsfeed ; Xfeed
Hsol ðTclimb ; Xclimb Þ Xclimb   
 Xfeed
1
Hevc Tclimb ; Poper Xclimb þ 2ð L
Lh ÞaUclimb LMTDclimb ¼ 0; abLh Uh LMTDh þ ð L
Lh ÞUclimb LMTDclimb þ LUfall LMTDfall c þ mstm DHsuph

þ ð1
qualout ÞDHcon ¼ 0:

ð30Þ

ð31Þ

The new equation in this case is the heat balance for the solution region where no phase change occurs


mfeed Hsol ðTfeed ; Xfeed Þ
Hsol Tsfeed ; Xfeed þ 2Lh aUh LMTDh ¼ 0:

ð32Þ

Since the overall heat-transfer coefficient for the solution single-phase region is considerably smaller than the one associated with the evaporating region, the possibility of condensate cooling for the subcooled feed was not considered in this work.

315

3.5. Computation of outlet streams Once the appropriate set of equations is solved, the evaporator outlet streams can be computed. For the concentrate flow rate, Fconc , depending on the condensate temperature, we have, from Fig. 1  nmfall if Tcond ¼ Tsstm ; Fconc ¼ ð33Þ nmsol if Tcond < Tsstm : 4 Moreover, for the evaporated flow rate, Fev , the situations analysed in the previously described calculation strategy give rise to the following relations: 8 nðmevc þ mevf Þ if Tfeed 6 Tsfeed and Tcond ¼ Tsstm ; > > > P4 > if Tfeed ¼ Tsfeed and Tcond < Tsstm ; n i¼1 mev > i > > < nðmevfl þ mevc þ mevf Þ stm Fev ¼ Tsfeed and ð34Þ >

if Tfeed >  Tcond ¼ Ts ; > P > 4 ev > > n mevfl þ i¼1 mi > > : if Tfeed > Tsfeed and Tcond < Tsstm : 3.6. Calculation of the overall heat-transfer coefficient Whichever the feed temperature, the overall heattransfer coefficient is always estimated by adding up the individual resistances to the heat-transfer process 
1  1 1 l U¼ þ ; ð35Þ þ hsol hstm k where hsol and hstm are the heat-transfer coefficients on the solution and steam sides, respectively, which are computed using literature correlations. For the steam side, the equation developed by Wang and Zhao (1993) is utilised  0:248  0:983 Rel ql 0:33 Nu ¼ 0:00115 Prl ; ð36Þ c qv
 Cpl Tsstm
Tw


c¼ ð37Þ
1 DHcon 1 þ 0:68Cpl Tsstm
Tw DHcon in which ql and qv are the liquid and vapour densities, respectively; Cpl is the liquid heat capacity; Nu is the Nusselt number; Prl is the liquid phase Prandtl number and Rel is the liquid phase Reynolds number, defined as Rel ¼

Gð1
qualout ÞDeq ; ll

ð38Þ

where ll is the liquid viscosity, G is the mass flow rate per unit area in the channel and Deq is the channel equivalent diameter Deq ¼

2wb wþb

ð39Þ

being w the channel width and b the distance between two consecutive plates.

316

C.P. Ribeiro Jr., M.H. Ca~no Andrade / Journal of Food Engineering 54 (2002) 309–320

For the solution side, when no phase transition takes place, the heat-transfer coefficient is estimated by the equation proposed by Buonopane and Troupe (1969) jl hSP ¼ a Reb Pr0:40 ð40Þ Deq in which jl is the liquid thermal conductivity and a; b are constants whose values depend upon the geometry of the ribs on the plate surface. No specific equation for the two-phase heat transfer to water or aqueous solutions in ribbed channels was found in the literature. Although true that some correlations for the evaporation of refrigerants in plate-type heat exchangers have already been proposed (Azarskov, Danilova, & Zemskov, 1981; Nakaoka & Uehara, 1988; Uehara, Stuhltr€ ager, Miyara, Murakami, & Miyazaki, 1997; Yan & Lin, 1999), such correlations were not regarded as suitable for the problem considered here on account of two reasons. First, the physical properties of the liquids employed for obtaining these correlations (ammonia, refrigerant R-134a, freon 22) differ considerably from the ones related to the fluids focused in this work; second, plates used in refrigerant evaporators are quite different from the ones associated with the vacuum concentration of heat-sensitive materials, insofar as the corrugations on their surfaces, which strongly influence heat-transfer characteristics, are concerned (Kumar, 1993). However, probably as a result of the industrial hegemony of tubular units, a vast number of equations for predicting the boiling heat-transfer coefficient in tubes can be found in the literature, from which the ones presented by Bjorge, Hall, and Rohsenow (1982), Chen (1966), Chun and Seban (1971), G€ usewell, Behrendt, Hoffmann, and Klostermann (1987), Kandlikar (1990), Klimenko (1988) and Klimenko (1990) are just a few examples. Kandlikar (1990), in particular, after correlating over 5000 experimental data on boiling heattransfer of 10 different fluids, obtained the following formula for the ratio between the two- and single-phase heat-transfer coefficients in tubes: hTP ¼ c1 Coc2 ð25Frl Þc5 þ c3 Boc4 Ffl ; hSP

Eq. (41) for computing the heat-transfer coefficient on the side of the evaporating solution. For all correlations, the channel equivalent diameter is taken as the characteristic dimension, being the solution dimensionless groups calculated with the average values of physical properties at inlet and outlet conditions.

4. Results and discussion Since no previous work on the steady-state simulation of climbing-falling-film plate evaporators was found in the literature, the validation of the developed algorithm by comparing its results with other authors’ ones could not be performed. Moreover, although a set of typical operating data for a plate evaporation unit was provided by Morgenroth et al. (1997), this could not be used for validating the algorithm either, as it lacked some of the calculation input data, such as feed temperature, operating pressure and total number of plates. Consequently, in order to provide a means for analysing the veracity of the algorithm results, an example calculation was performed, varying some operating conditions, namely the feed temperature and the steam flow rate, so as to cover all the possibilities the algorithm takes into account. A milk concentrating evaporator, operating at 30.4 kPa with stainless steel plates, was chosen as the sample unit. Its main features are listed in Table 2. For performing the calculations, water thermophysical properties were estimated using the models given by Daubert and Danner (1985). As regards milk, its density was calculated by means of the equation presented by Cosme, Guerrero, and Velez (1997); for its heat capacity, the model proposed by Fernandez-Martın (1972) was employed. The viscosity of milk was computed utilising the equation of Ruiz and Canovas (1997). In addition, the model developed by Fernandez-Martın and Montes (1972) was chosen for estimating the thermal conductivity of milk, whereas, for the boiling point

ð41Þ

where the constants c1 –c5 depend on the boiling regime (convective or nucleate boiling), while the parameter Ffl is fluid-dependent. If one now assumes that, despite the differences for the individual heat-transfer coefficients in cylindrical and rectangular channels, the ratio between the two- and single-phase coefficients is approximately the same for both geometries, Eq. (41) may be employed to calculate the film coefficient on the solution side of a plate evaporator, provided that hSP is estimated by means of a correlation associated with ribbed rectangular channels. Therefore, in this work, Eq. (40) was substituted into

Table 2 Characteristics of the sample unit simulated Feed flow rate, Ffeed (kg/h) Feed solids content, Xfeed (w/w) Steam temperature, Tstm ð°CÞ Steam pressure, Pstm (kPa) Operating pressure, Poper (kPa) Plate heat-transfer area, a ðm2 =mÞ Plate length, L (m) Plate width, w (m) Plate thickness, l (mm) Distance between plates, b (mm) Plate material Number of four-plate sets, n

4000.0 0.1300 75.0 38.56 30.40 0.30 0.80 0.40 1.50 5.00 AISI 316 steel 20

C.P. Ribeiro Jr., M.H. Ca~no Andrade / Journal of Food Engineering 54 (2002) 309–320

317

Table 3 Simulation results for the different operating conditions of the sample unit Case

Tfeed (°C)

Fstm (kg/h)

Xout

Tcond (°C)

L1 =L

Lh =L

1 2 3 4 5 6 7 8 9 10

69.4 69.4 69.4 69.4 70.4 72.4 74.4 68.4 64.4 60.4

789.0 770.0 700.0 600.0 789.0 789.0 789.0 789.0 789.0 789.0

0.16075 0.15990 0.15675 0.15235 0.16111 0.16182 0.16252 0.15890 0.15397 0.15135

75.0 74.2 72.0 70.4 74.7 74.2 73.8 75.0 75.0 75.0

0.0000 0.0167 0.0790 0.1680 0.0052 0.0154 0.0249 0.0000 0.0000 0.0000

0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0673 0.2630 0.3995

rise, a linear regression of the experimental data provided by Fergusson (1989) was carried out. As shown in Table 3, 10 different operating conditions regarding feed temperature and steam flow rate were investigated. For each condition, the calculated values of the following variables are presented: total solids content of the concentrate, condensate temperature and fractions of the plates length used for cooling the condensate ðL1 =LÞ and heating the feed ðLh =LÞ. Case 1 in Table 3, with saturated feed ðLh ¼ 0 mÞ, total condensation and no condensate subcooling ðL1 ¼ 0 mÞ, corresponds to the simplest situation and will, therefore, be the base of comparison for the other conditions. Going from cases 1–4, the steam flow rate is progressively reduced for a constant feed temperature, which means a lower amount of available latent heat. Hence, a smaller degree of concentration is expected. The results presented in Table 3 for the solids content of the outlet stream do confirm such expectation. Furthermore, simulation results indicate condensate subcooling for cases 2–4, whose extent augments as the steam flow rate decreases. This observation is also in agreement with theoretical expectations and can be understood as follows: convection on the steam side is the lowest resistance to heat-transfer in the equipment, so that no appreciable diminution in the overall heattransfer coefficient is verified when the steam flow rate is reduced. Being U approximately the same, the lower the amount of steam available, the smaller the heat-transfer area required for its total condensation, or, in other words, the greater the portion of the plates area used to cool the condensate. Such fact is clearly shown in Table 3 by the augmentation of the L1 =L ratio, which indicates that, for case 4, almost 17% of the plates area is utilised for condensate subcooling. It should be noted that this increase in the sensible heat transferred to the evaporating solution is by no means enough to compensate for the diminution in the latent heat, in view of the fact that the former is considerably inferior to the later, and, consequently, a lower solids concentration at the outlet is obtained.

If one now focuses his attention on cases 5–7, comparing them with case 1, an analysis of the calculation mode related to superheated feed streams can be performed. In this case, for a fixed steam flow rate, the feed temperature, initially equal to the saturation value, is continuously increased. As the equipment operating pressure is fixed, an increase in the feed temperature beyond its saturation value will lead to flash vaporisation, being vapour formed at the expense of the liquid, which results in a higher solids concentration at the outlet. The higher the feed temperature, the greater the amount of solvent evaporated in this process, justifying, thus, the results presented in Table 3 for the outlet solids content of the solution in cases 5–7. Another interesting outcome of the use of a superheated feed is a decrease in the condensate outlet temperature, as shown in Table 3. On account of the flash vaporisation, the amount of evaporated solvent in the channels rises, which, in turn, leads to an increase in the heat-transfer coefficient on the solution side (Kandlikar, 1990; Klimenko, 1988; Yan & Lin, 1999). Convection on the solution side is an important resistance to heat transfer in the equipment and, therefore, the overall heat-transfer coefficient for both the climbing and falling-film sections also increases. With a higher U value, the same steam flow rate can be totally condensated with a smaller heat-transfer area, leaving, thereby, part of the plates area for condensate subcooling to take place. The three last situations in Table 3 refer to subcooled feed streams. For cases 8–10 the feed temperature was systematically reduced in comparison to the reference value, that is, case 1. As no evaporation can occur until saturation temperature is reached, a fraction of the plates’ area has to be utilised to heat the solution, accordingly reducing the amount of solvent evaporated in the unit, which is precisely the pattern exhibited by the results presented in Table 3. The colder the feed stream, the bigger the plates’ length associated with the solution heating. In case 10, for instance, due to a 13% decrease in the feed temperature, the single-phase region accounts for almost 40% of the climbing-film section, which is

318

C.P. Ribeiro Jr., M.H. Ca~no Andrade / Journal of Food Engineering 54 (2002) 309–320

Table 4 Operating conditions of the industrial plate evaporators simulated Variables

EV-101

EV-102

EV-103

EV-104

EV-301

EV-302

Feed flow rate, Ffeed (kg/h) Feed solids content, Xfeed (w/w) Feed temperature, Tfeed (°C) Steam flow rate, Fstm (kg/h) Steam temperature, Tstm (°C) Number of four-plate sets, n Pressure, Poper (kPa)

12,090 0.125 70.6 4579 80.0 54 26.14

7249 0.209 63.6 1201 65.5 37 20.87

5871 0.258 58.2 1238 61.0 45 16.11

4466 0.339 54.4 1260 56.0 48 9.32

6054 0.125 77.1 2114 93.0 20 26.14

3561 0.213 62.3 1881 66.0 21 8.00

highly prejudicial to the equipment performance, reducing the evaporate flow rate in almost 27%. The results presented so far indicate that all the calculation modules of the developed algorithm provide answers in accordance with theoretical reasoning. However, in order to prove the practical usefulness of the algorithm, it should be employed in the simulation of at least one real unit. Therefore, with the aid of the developed algorithm, the steady-state operation of six different industrial climbing-falling-film plate evaporators was simulated. The chosen units are part of the milk powder plant from Embare Ind ustrias Alimentıcias S.A., a dairy industry located in Minas Gerais, Brazil. There are two concentrating systems at Embare’s plant, which together produce 1900 kg/h of milk powder when operating at full capacity. The first concentrating system (CS1) comprises four evaporation effects, EV101 to EV-104, and is responsible for processing two thirds of the total amount of raw milk fed to the plant, whereas there are only two evaporation effects in the

second concentrating system (CS2), EV-301 and EV302. These evaporators are APV units with stainless steel (AISI 316) plates showing staggered hemispherical ribs on their surface. Their operating conditions are outlined in Table 4, whilst Table 5 presents some geometrical characteristics of the plates. In these simulations, water and milk thermophysical properties were estimated by means of the same models specified previously for the sample unit simulation. However, for the thermal conductivity of milk, due to the validation range of the equation presented by Fernandez-Martın and Montes (1972), when the total solids content was superior to 0.37, this property was Table 5 Geometrical characteristics of the industrial evaporator plates Length, L (cm) Width, w (cm) Thickness, l (mm) Distance between plates, b (mm) Heat-transfer area, a ðcm2 =cmÞ

92.5 53.0 1.6 6.0 53.1

Table 6 Comparison between calculated values and operating data for the CS1 evaporators Variable

Real

Simulation

Error (%)

EV-101 Concentrate flow rate, Fconc (kg/h) Evaporated flow rate, Fev (kg/h) Outlet solids content, Xout (w/w) Condensate temperature, Tcond (°C)

6648 4704 0.205 70.0

6668 4680 0.204 66.2

0.3 )0.5 )0.5 )5.4

EV-102 Concentrate flow rate, Fconc (kg/h) Evaporated flow rate, Fev (kg/h) Outlet solids content, Xout (w/w) Condensate temperature, Tcond (°C)

5410 1238 0.252 62.0

5431 1214 0.251 63.2

0.4 )1.9 )0.4 1.9

EV-103 Concentrate flow rate, Fconc (kg/h) Evaporated flow rate, Fev (kg/h) Outlet solids content, Xout (w/w) Condensate temperature, Tcond (°C)

4150 1260 0.328 58.5

4164 1245 0.327 58.0

0.3 )1.2 )0.3 )0.8

EV-104 Concentrate flow rate, Fconc (kg/h) Evaporated flow rate, Fev (kg/h) Outlet solids content, Xout (w/w) Condensate temperature, Tcond (°C)

2840 1310 0.480 51.0

2847 1302 0.479 48.1

0.2 )0.6 )0.2 )5.7

C.P. Ribeiro Jr., M.H. Ca~no Andrade / Journal of Food Engineering 54 (2002) 309–320

319

Table 7 Comparison between calculated values and operating data for the CS2 evaporators Variable

Real

Simulation

EV-301 Concentrate flow rate, Fconc (kg/h) Evaporated flow rate, Fev (kg/h) Outlet solids content, Xout (w/w) Condensate temperature, Tcond (°C)

3403 2272 0.200 67.3

3457 2213 0.197 67.7

1.6 )2.6 )1.5 0.6

EV-302 Concentrate flow rate, Fconc (kg/h) Evaporated flow rate, Fev (kg/h) Outlet solids content, Xout (w/w) Condensate temperature, Tcond (°C)

1420 1983 0.480 54.5

1432 1970 0.476 45.4

0.8 )0.7 )0.8 )16.7

estimated using the equation proposed by More and Prasad (1988). A comparison between simulation results and the field operating data for the units is given in Tables 6 and 7. It can be seen that the calculated and real values agree quite well. For the CS2, milk solids contents at the outlet of all units are predicted with a maximum deviation of 1.5%, whereas, for the CS1, this value is as low as 0.5%. Similarly, with regard to the water evaporation rates, the highest difference between simulated and operating data was equal to only 2.6%. As for the condensate temperature, which is associated with the higher errors, the mean deviation for all units is 5.2%. It should be pointed out that the condensate temperature deviation for the EV-302, 16.7%, is close to the mean deviation value reported by Kandlikar (1990) for his correlation, 15.9%.

5. Conclusion An algorithm for the steady-state simulation of climbing-falling-film plate evaporators, based on the simultaneous resolution of the heat and mass balances associated with the heat-transfer rate equation, was developed and implemented in FORTRAN language. As no previous work on the simulation of such evaporator was found in the literature, the adequacy of the calculation procedure was demonstrated by the simulation of a sample unit operating at different conditions. The algorithm was then successfully utilised for the performance prediction of six different industrial evaporators employed in the concentrating step of a Brazilian milk powder plant.

Acknowledgements The authors would like to thank FAPEMIG for the financial support that enabled the development of this

Error (%)

work, and also the company Embare Ind ustrias Alimentıcias S.A., for the industrial data provided.

References Azarskov, V. M., Danilova, T. N., & Zemskov, B. B. (1981). Heat transfer in plate evaporators of different geometries. Kholodil’naia Tekhnika, 4, 25–31 (in Russian). Billet, R. (1988). Evaporation. In W. Gerhartz (Ed.), Ullmann’s encyclopedia of industrial chemistry, Vol. B3 (5th ed., pp. 3.1–3.35). New York: VCH. Bjorge, R. W., Hall, G. R., & Rohsenow, W. M. (1982). Correlation of forced convection boiling heat transfer data. International Journal of Heat and Mass Transfer, 25(6), 753–757. Brotherton, F. (1994). Evaporation in plate type heat exchangers. Heat Recovery Systems & CHP, 14(5), 555–561. Buonopane, R. A., & Troupe, R. A. (1969). A study of the effects of internal rib and channel geometry in rectangular channels – part II: heat transfer. AIChE Journal, 15(4), 592–596. Carter, A. L., & Kraybill, R. R. (1966). Low pressure evaporation. Chemical Engineering Progress, 60(2), 99–110. Chen, C. J. (1966). Correlation for boiling heat transfer to saturated fluids in convective flow. Industrial and Engineering Chemistry Process Design and Development, 5(3), 322–329. Choudhary, R., Ali, S., & Das, H. (1996). A heat transfer model for performance evaluation of single effect multistage mechanical vapour recompression evaporator. Journal of Food Science and Technology, 33(4), 308–312. Chun, K. R., & Seban, R. A. (1971). Heat transfer to evaporating liquid films. Journal of Heat Transfer, 93(4), 391–396. Cosme, A. L., Guerrero, J. A., & Velez, J. F. (1997). Evaluacion de propriedades fisicoquimicas de leche concentrada. Informacion Tecnologica, 8(1), 35–40. Daubert, T. E., & Danner, R. P. (1985). Data compilation tables of properties of pure compounds. American Institute of Chemical Engineers. Fergusson, P. H. (1989). Developments in the evaporation and drying of dairy products. Journal of the Society of Dairy Technology, 42(4), 94–101. Fernandez-Martın, F. (1972). Influence of temperature and composition on some physical properties of milk and milk concentrates. I. heat capacity. Journal of Dairy Research, 39(65), 65–73. Fernandez-Martın, F., & Montes, F. (1972). Influence of temperature and composition on some physical properties of milk and milk concentrates. III. thermal conductivity. Milchwissenschaft, 27(12), 772–776. G€ usewell, M., Behrendt, D., Hoffmann, M., & Klostermann, W. (1987). Evaporators. In S. Weiß (Ed.), Verfahrenstechnische

320

C.P. Ribeiro Jr., M.H. Ca~no Andrade / Journal of Food Engineering 54 (2002) 309–320

Berechnungsmethoden Teil 1 – W€arme€ubertrager (pp. 125–198). Weinheim: VCH (in German). Gray, R. M. (1971). The plate evaporator. Journal of Applied Chemistry and Biotechnology, 21(12), 359–362. Hellstr€ om, J. (1994). Industrial plate evaporators. Heat Recovery Systems & CHP, 14(5), 549–554. Holland, C. D. (1975). Fundamentals and modelling of separation processes. New Jersey: Prentice-Hall. Kandlikar, S. G. (1990). A general correlation for saturated two-phase flow boiling heat transfer inside horizontal and vertical tubes. Journal of Heat Transfer, 112, 219–228. Klimenko, V. V. (1988). A generalised correlation for two-phase forced flow heat transfer. International Journal of Heat and Mass Transfer, 31(3), 541–552. Klimenko, V. V. (1990). A generalised correlation for two-phase forced flow heat transfer – second assessment. International Journal of Heat and Mass Transfer, 33(10), 2073–2088. Kreith, F. (1973). Principles of heat transfer. New York: Intext Educational. Kumar, H. (1993). Evaporation in plate heat exchangers. AIChE Symposium Series, 295(89), 211–222. Lambert, R. N., Joye, D. D., & Koko, F. W. (1987). Design calculations for multiple-effect evaporators. 1. Linear method. Industrial and Engineering Chemistry Research, 26(1), 100–104. Licha, H., Valentin, P., Wersel, M., & Witte, G. (1989). The plate evaporator – a new method in evaporation technology. Zuckerindustrie, 114(10), 785–798. More, G. R., & Prasad, S. (1988). Thermal conductivity of concentrated whole milk. Journal of Food Process Engineering, 10, 105–112.

Morgenroth, B., Jayatilaka, D., & Punter, G. (1997). Development of plate evaporator technology, the market place and the choice for the sugar technologist. Zuckerindustrie, 122(9), 691–699. Nakaoka, T., & Uehara, H. (1988). Performance test of a shell-andplate type evaporator for OTEC. Experimental Thermal and Fluid Science, 1, 283–291. Ruiz, J. F. V., & Canovas, G. V. B. (1997). Effects of concentration and temperature on the rheology of concentrated milk. Transactions of the ASAE, 40(4), 1113–1118. Smith, J. M., Van Ness, H. C., & Abbott, M. M. (1996). Introduction to chemical engineering thermodynamics (5th ed., p. 199). Singapore: McGraw-Hill. Standiford, F. C. (1963). Evaporation is a unit operation. Chemical Engineering, 70, 158–176. Uehara, H., Stuhltr€ager, E., Miyara, A., Murakami, H., & Miyazaki, K. (1997). Heat transfer and flow resistance of a shell and plate type evaporator. Journal of Heat Transfer, 119, 160– 164. Usher, J. D. (1965). Plate evaporators – design and application. Chemical and Process Engineering, 46(5), 229–234. Wang, Z., & Zhao, Z. (1993). Analysis of performance of steam condensation heat transfer and pressure drop in plate condensers. Heat Transfer Engineering, 14(4), 32–41. Westphalen, D.L. (1999). Modelling, simulation and optimisation of evaporation systems. Ph.D. Thesis, UNICAMP, Campinas, 228 pp (in Portuguese). Yan, Y. Y., & Lin, T. F. (1999). Evaporation heat transfer and pressure drop of refrigerant R-134a in a plate heat exchanger. Journal of Heat Transfer, 121(2), 118–127.