Available online at www.sciencedirect.

com

Applied Thermal Engineering 28 (2008) 532–546

www.elsevier.com/locate/apthermeng

Modelling and simulation of fluid flow and heat transfer

in the convective zone of a power-generation boiler
Antonio Gómez a, Norberto Fueyo a,*
, Luis Ignacio Dı́ez b

a
Fluid Mechanics Group, CPS, Universidad de Zaragoza, Marı́a de Luna 3, 50018 Zaragoza, Spain
b
Fundación CIRCE, Marı́a de Luna 3, 50018 Zaragoza, Spain

Received 26 September 2006; accepted 30 April 2007

Available online 18 May 2007

Abstract

The convective zone of a power-station boiler is a complex piece of engineering equipment, which comprises a multiplicity of inter-
connected heat-exchanging elements.
This paper presents a mathematical model of the convective zone, which allows for the calculation of the shell-side flow and the shell-
side, tube-side and tube-wall, thermal fields, and of the shell-tube heat-exchange. The model is solved using conventional CFD tech-
niques, in which the individual tubes are treated as sub-grid features. A geometrical model is used to describe the multiplicity of
heat-exchanging structures, and the interconnections among them.
The CFD model is first validated by comparison with simple heat exchanger geometries which can be solved with analytical methods,
and then applied to an actual 350 MW(e) power-station boiler.
Ó 2007 Elsevier Ltd. All rights reserved.

Keywords: Heat exchangers; Power-station boiler; Convective zone; Heat transfer; CFD; Modelling

1. Introduction ciency and durability. For the modelling to be accurate

enough, a desirable feature of the model should be the abil-
The shell-and-tube heat exchanger is a ubiquitous piece of ity to calculate simultaneously the (spatially-varying) tem-
equipment in the power-generation and processes industries. perature of both the shell and the tube fluids.
In this type of heat exchanger, heat is transferred to (or The modelling of shell-and-tube heat exchangers using
from) a fluid flowing through a tube bundle from (or to) Computational Fluid Dynamics techniques dates back to
another fluid which flows around the tubes through the shell. the early days of the discipline. It is difficult therefore to
Both fluids are everywhere separated by the tube wall (usu- trace it to a seminal paper; however, the technique
ally metal-made) and never become in direct contact. described in the one by Patankar and Spalding [1] was per-
The convective zone of a power-station boiler is a com- haps the first one capable of calculating simultaneously the
plex assembly of inter-connected heat exchangers. Around thermal fields in both fluids and in the tube metal, and can
50% of the cycle heat generated in the boiler may be trans- be regarded as the framework from which successive con-
ferred in this zone, and hence the interest of the modelling tributions were developed.
of fluid flow and heat transfer. An accurate computer A literature survey, however, shows that most of the
model can provide convenient and economical assistance published CFD-type simulations of large, shell-and-tube
in the design and in the operation of the convective zone, heat exchangers pay much attention to the prediction of
and can contribute significantly to the improvement of effi- the aerodynamics (such as flow distribution or pressure
losses), but fewer contributions can be found in which
*
Corresponding author. the thermal fields and heat transfer are also considered.
E-mail address: Norberto.Fueyo@unizar.es (N. Fueyo). This is the case in the early paper by Rhodes and Carlucci

1359-4311/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved.
doi:10.1016/j.applthermaleng.2007.04.019
 A. Gómez et al. / Applied Thermal Engineering 28 (2008) 532–546 533

Nomenclature

At cross-sectional surface-area of tubes in a cell ui ith component of velocity vector (m/s)

(m2) utE tube-side-fluid velocity across the east cell-face
bi inlet mass flow rate through element i (if a mani- (m/s)
fold) (kg/s) utW tube-side-fluid velocity across the west cell-face
Cpt specific heat of tube-side fluid (J/kg K) (m/s)
Cl constant in the k– turbulence model U absolute value of velocity vector (m/s)
Di inner tube diameter (m) v~s shell-side-fluid velocity-vector (m/s)
Do outer tube diameter (m) v~t tube-side-fluid velocity-vector (m/s)
Fi ith component of friction term (kg/m2 s2) VP volume of cell P (m3)
~
F friction term accounting for shell-side pressure-
losses in the tube bank (kg/m2 s2) Greek symbols
g tube-side mass flow rate per cell (kg/s) at!w volumetric heat-transfer coefficient from the
hs, ht, hw specific enthalpy of shell-side fluid, tube-side tube-side fluid to the tube wall (W/m3 K)
fluid and tube wall, respectively (J/kg) as!w volumetric heat-transfer coefficient from the
htE specific enthalpy of tube-side fluid in the east cell shell-side fluid to the tube wall (W/m3 K)
00
(J/kg) as!w heat-transfer coefficient from the tube-side fluid
htP specific enthalpy of tube-side fluid in the current to the tube wall (W/m2 K)
00
cell (J/kg) at!w heat-transfer coefficient from the tube-side fluid
htT specific enthalpy of tube-side fluid at the previ- to the tube wall (W/m2 K)
ous time-step (J/kg) dij connectivity parameter (=1 if element i feeds ele-
htW specific enthalpy of tube-side fluid in the west ment j and 0 otherwise)
cell (J/kg) Dt time interval (s)
hwt specific enthalpy of tube-side fluid at the tube-  dissipation rate of the turbulent kinetic energy
wall temperature (J/kg) (m2/s3)
k turbulent kinetic energy (m2/s2) s fraction of local space occupied by the shell-side
kt thermal conductivity of tube-side fluid (W/m K) fluid
l characteristic length scale in the tube bundle (m) t fraction of local space occupied by the tube-side
mi mass flow rate through element i (kg/s) fluid
Nu Nusselt number (–) w fraction of local space occupied by the tube wall
nj number of elements fed from element j (–) Cs effective diffusivity in the shell-side enthalpy
p pressure (Pa) equation (kg/ms)
P tube pitch in the cross-flow direction (m) ls molecular (dynamic) or effective (molecular plus
P average pitch in the tube bundle (m) turbulent) viscosity of the shell-side fluid (kg/
Prs, Prw Prandtl numbers of shell-side fluid calculated ms)
at the local fluid temperature and at the wall lT turbulent viscosity of the shell-side fluid (kg/ms)
temperature, respectively qs, qt, qw density of shell-side fluid, tube-side fluid and
Prt Prandtl number of tube-side fluid calculated at tube wall, respectively (kg/m3)
the local fluid temperature nui friction factor for the ith component of velocity
Rfs shell-side fouling resistance (m2 K/W) vector (m1)
Rft tube-side fouling resistance (m2 K/W)
Rs shell-side convection resistance (m2 K/W) Subscripts
Rt tube-side convection resistance (m2 K/W) E east cell
Res, Ret local Reynolds number of the shell-side and f fouling
tube-side fluids P current cell
S/V ratio of tube surface-area to volume in the bank s shell-side fluid
(m2/m3) t tube-side fluid
t time (s) w tube wall
Ts, Tt, Tw shell-side-fluid temperature, tube-side-fluid W west cell
temperature and tube-wall temperature (K)

[2], which reports the first aerodynamic calculation of a tation and discretisation, due to the limitations imposed by
model heat exchanger, and a comparison with experimen- the computers available at the time, the mathematical
tal measurements. Despite the crude geometrical represen- model had many of the features still used nowadays for
534 A. Gómez et al. / Applied Thermal Engineering 28 (2008) 532–546

calculating the tube-bundle aerodynamic performance, putationally impractical for complex arrays of intercon-
such as the use of porosities, friction factors and effective nected banks. The present model, on the other hand,
viscosities to represent the effect of the tube bundle on allows for the simultaneous calculation of the shell-side,
momentum transport. The problem was later revisited by tube-side and metal thermal fields even when several inter-
Theodossiou et al. [3]. connected heat exchangers have to be considered, as it is
The simulation of power-plant condensers, which for the the case in power-station boilers.
purpose of this paper can be regarded as a derivative of the
heat exchanger, has attracted much attention in recent 2. Conservation equations
years. One of the most comprehensive three-dimensional
models is that used by Zhang et al. [4], Zhang and Zhang The present mathematical model can be best viewed as
[5] and Zhang and Bokil [6] for power-plant condensers, representing three ‘‘superimposed’’ spaces or ‘‘worlds’’,
albeit the tube-side fluid is treated in a less general manner which coexist in physical space: the shell-side of the heat
than that employed in this paper. Ormiston et al. [7,8] pro- exchanger, the tube side and the metal. Conservation equa-
vide a literature survey of recent CFD models of power- tions are solved for the relevant magnitudes in these three
plant condensers, as well as some modifications to the spaces, and include the appropriate exchange terms
solution algorithm which can lead to better convergence between them (e.g., heat transfer).
behaviour. Prithiviraj and Andrews [9,10] reformulate the Thus, for the shell-side, the relevant equations are con-
finite-volume equations for the tube bundle using control- tinuity, momentum and enthalpy conservation. The conti-
volume integrals of the model equations. The treatment nuity equation is
of the tube-side enthalpy is again rather simple, since the
os qs
condensers considered in their paper are single-pass ones. þ r  ðs qs~
vs Þ ¼ 0 ð1Þ
No previous multi-dimensional simulations of the fluid ot
flow and heat transfer in the convective zone of a power-sta- where s is the void fraction, i.e., the fraction of local space
tion boiler, addressing both the shell-side and the tube-side available to the shell-side fluid; qs is the shell-fluid density;
fields have been traced in the open literature. The paper by and ~vs is the shell-fluid interstitial-velocity vector.
Tu et al. [11] addresses the CFD modelling of the fly-ash The shell-side-momentum equation is
flow (without heat transfer) in power-utility boilers (includ- vs Þ
oðs qs~
ing both the furnace and convective section) using Eulerian– þ r  ðs qs~
vs~ vs Þ ¼ s rp  ~
vs Þ  r  ðs ls r~ F ð2Þ
ot
Eulerian models. However, the emphasis is placed on the
multi-phase, gas-particle model, and the representation of where ls is either the molecular viscosity (if the flow is lam-
the convective section appears to be a simplified one. Coelho inar) or an effective (molecular plus turbulent) one if the
[12] has reported results from a CFD model of a boiler sim- flow is turbulent and an eddy-viscosity model is used; p is
ilar to the present one, but, unlike this work, the modelling is the pressure; and ~F is a friction term (specified below) that
restricted to the shell side, and either the tube-side tempera- accounts for the pressure losses in the tube bank.
ture is assumed or the overall heat transfer in a tube bank is The shell-side-enthalpy (hs) equation reads
known and distributed along the bank. Diez et al. [13] oðs qs hs Þ
address the modelling of the same boiler as featured in this þ r  ðs qs v~s hs Þ  r  ðs Cs rhs Þ ¼ as!w ðT w  T s Þ
ot
work, but using a zonal method for the gas flow and
ð3Þ
lumped-parameter heat exchanger network for the tube-
side, together with CFD calculations of the combustion where Cs is a diffusion coefficient; as!w is a volumetric heat
chamber to provide inlet conditions for the zonal model. transfer coefficient from the shell fluid to the tube wall,
Such a simplified model, it is argued, has the attractive merit with units W/m3 K; Ts is the local shell-fluid temperature;
of being usable on line. A similar method is used by Niu and and Tw is the local tube-wall temperature. The use of a vol-
Wong [14]. Huang et al. [15] model the combustion chamber umetric heat-transfer coefficient in this equation is required
and heat exchanger of a domestic boiler; although few for dimensional reasons; it is calculated from the usual,
details of the mathematical model and numerical procedure per-unit-surface one using the ratio of tube surface to shell
are given, the method employed is probably similar to the volume as detailed below. Other heat sources, such viscous
one used in the present work. heating or radiative heat-exchange, have been neglected in
This paper presents a mathematical model that can be this equation.
used for the CFD simulation of fluid flow and heat transfer If the flow in the shell side is turbulent, as it is usually the
in the convective section of power-station boilers. The case, the above equations can be interpreted as Favre-
model, which builds on the method by Patankar and averaged. Then, if suitable eddy-viscosity turbulence model
Spalding [1], represents the tubes as sub-grid features; the is used to close the unknown correlations, ls can be regarded
meshing of individual tubes or periodically-repeating as an effective shell-side viscosity, and Cs as an effective shell-
tube-arrays, while useful to improve the design of tube side Prandtl number for hs. Both can be calculated from the
banks or to calculate their heat-exchanging or aerodynamic model parameters. The flow in a utility boiler is clearly tur-
performance, such as pressure losses (see, e.g., [16]), is com- bulent; in the application shown in Section 6 below, the Rey-
 A. Gómez et al. / Applied Thermal Engineering 28 (2008) 532–546 535

nolds number based on an average velocity and a boiler F i ¼ nui qs ui U ð6Þ

transversal dimension is in excess of 500,000. Turbulence
in such a highly complex physical domain is correspond- where ui is the ith component of the velocity vector, U is the
ingly a complex issue. Given the impossibility to mesh the absolute value of the velocity vector, and the friction fac-
individual tubes in the banks, we have resorted to the use tors nui are calculated similarly to Rhodes and Carlucci
of a k– turbulence model, with some modifications to [2] as
  
impose a length scale within the tube bank. While more fk 1  s
sophisticated turbulence models are available, their use for nui ¼ 2 ð7Þ
P 1
the problem addressed in this paper does not justify the sub-
stantial additional computational expense, since the tube for the parallel-flow direction and
bank is nevertheless represented as a sub-grid feature.   2  
In the tube side, the flow is considered (locally) one- f? P s 1  s
nui ¼ 2 ð8Þ
directional (viz., in the direction of the tube orientation), P P  Do 1
which allows to obtain the local mass flow rate by simple for the cross-flow direction.
geometrical calculations (see below). Regarding the tube- In Eqs. (7) and (8), the factors fk and f? are functions of
side momentum equations, their only use would be to cal- the local Reynolds number [2]; P is the tube pitch in the
culate tube-side head losses. However, if head losses in the cross-flow direction; P is an average pitch; Do is the tube
tube-side circuit are quantity of interest, they can be calcu- outer diameter; s is the local void fraction, and  is the void
lated using simpler alternative methods, e.g., pipe-flow for- fraction in the local tube-bank. These both can be readily
mulae. Thus, in the current application, the tube-side calculated from geometrical data.
momentum equations are not required. The only relevant An alternative means of combining the contributions
conservation equation for the tube side is therefore the from the principal directions to the pressure loss is the ten-
one for enthalpy (ht), which reads sorial approach suggested by Butterworth [18].
oðt qt ht Þ
þ r  ðt qt~
vt ht Þ ¼ at!w ðT w  T t Þ ð4Þ
ot 3.2. Heat-transfer coefficients
where t is the fraction of local space occupied by the tube-
side fluid; ~vt is the (one-directional) tube-side-fluid velocity For the calculation of the volumetric heat transfer coef-
vector; at!w is a volumetric heat-transfer coefficient (from ficient at!w in Eq. (4), we first compute the (per-unit-sur-
the shell fluid to the tube wall), with units W/m3 K; and Tt face) heat-transfer coefficient a00t!w by considering the
is the (bulk) tube-side-fluid temperature. heat-transfer resistances due to heat convection and
Streamwise heat conduction (molecular or turbulent) fouling:
has been neglected in the above equation, since in all cases 1
it will be much less important than convection. (Cross- ¼ Rt þ Rft ð9Þ
a00t!w
stream conduction is accounted for in the volumetric
heat-transfer coefficient at!w.) The convection resistance Rt is obtained using a Nusselt-
Finally, for the metal (the tube wall), the only relevant number correlation:
equation is the energy-conservation one, which, neglecting
1 Nuk t
axial heat-conduction, reads ¼ ð10Þ
Rt Di
oðw qw hw Þ
¼ as!w ðT s  T w Þ þ at!w ðT t  T w Þ ð5Þ with kt being the tube-side-fluid thermal-conductivity and
ot
Nu ¼ 0:023Re0:8
t Pr t
0:4
[19]. Transport and thermodynamic
where w is the fraction of local space occupied by the tube properties for this correlation are calculated at the local
wall. tube temperature, which is an output from the model.
An alternative framework to that presented here for Finally, the dimensions of at!w in Eqs. (4) and (5), i.e.,
addressing the space-sharing feature of the flow under con- W/m3 K, are recovered from those of a00t!w (W/m2 K) by
sideration is the use of different computational spaces for multiplying the coefficient by the surface-to-volume ratio
each of the ‘‘co-existing worlds’’. Beale and Zhubrin [17] in the bank, i.e., by the local ratio of heat-exchanging sur-
apply this concept in their simulation of Solid Oxide Fuel face S to physical volume V:
Cells.
S
at!w ¼ a00t!w ð11Þ
3. Auxiliary equations V
The shell-to-wall heat-transfer coefficient as!w is calculated
3.1. Pressure losses in a similar manner:

Pressure losses (~F in Eq. (2)) are calculated by means of 1

¼ Rs þ Rfs ð12Þ
friction factors n, as follows: a00s!w
536 A. Gómez et al. / Applied Thermal Engineering 28 (2008) 532–546

with the resistance coefficient Rs now obtained from the heat cannot be regarded as constant, and for an accurate
correlation proposed by Zhukauskas [20] for the Nusselt calculation of shell-to-tube heat transfer they must be com-
number: puted as a function of the local temperature and pressure.
 1=4 In this work, the IAPWS-IF97 formulation for water and
Prs steam [28] has been attached to the simulation code for this
Nu ¼ CRems Pr0:36
s ð13Þ
Prw purpose.
Here, the Reynolds number is calculated with the tube
external diameter and the maximum gas velocity within 4. Geometrical description
the tube bank. The gas Prandtl numbers Prs and Prw are
calculated respectively at the local gas temperature and at One of the main difficulties in simulating heat exchang-
the wall temperature, both obtained from the model. ers and convective zones is often their geometrical com-
The constants C and m depend on tube-bank configura- plexity, since they are frequently composed of a large
tion and Reynolds number; see [20] for details and for the number of interconnected tube banks. This is certainly true
range of validity for the correlation. of the convective zone of a power-station boiler.
It is difficult to reliably calculate the shell-side fouling- The approach embodied in the equations introduced
resistance from other flow parameters, despite past and above implies the treatment of tubes as sub-grid features.
on-going research efforts on fouling modelling and predic- This somehow alleviates the problem, since it circumvents
tion [21–26]. In the present work, empirical values, calcu- the need for meshing individual tubes; but some significant
lated from plant data, will be used. difficulties still remain, and this section describes how they
have been addressed in the present model.
The geometrical description used in this model considers
3.3. Turbulence essentially two types of element: tube banks and manifolds.
The tube bank is the main heat-exchanging item. They can
The modification of the turbulence in the presence of be ‘‘two-dimensional’’ (i.e., a single plane of tubes) or
tube bundles has been addressed by several authors in the ‘‘three-dimensional’’. In either case, they can additionally
past. When eddy-viscosity models are used, a common be present in the form of ‘‘arrays’’ of identical structures
approach is to modify the turbulent viscosity in the tube regularly arranged in space. Tube banks are characterised
bundle, e.g., by using a length scale (typically proportional by geometrical, physical and topological properties.
to the clearance between tubes) and a velocity fluctuation Among the geometrical properties, the main relevant ones
to obtain an effective viscosity or thermal diffusivity [2]. are pitch, external and internal diameters, position and
Another approach (see, e.g. [9]) modifies the source term dimensions. The physical properties include all the tube-
in the k and  equations to account for the enhanced pro- side-fluid properties, the tube material (and hence its prop-
duction and dissipation of turbulence in the tube bundle. erties), and the mass flow rate per tube. The all-important
In this work, a version of the first approach is used. The topological properties are two: the flow direction within the
gas effective-viscosity for points located within the tube tubes, by reference to the simulation axes, and the bank
bundle is calculated considering that the tube clearance is connectivity, i.e., to which other bank(s) each bank is con-
a characteristic length scale l. For the fluctuating velocity, nected and how.
k1/2 is used. Thus, One type of such connection (but not the only one) is the
lT ¼ C l qs lk 1=2 ð14Þ manifold (or header). This is defined as a virtual element
(i.e., without a geometrical representation in the model)
where Cl is a standard constant in the k– model [27]; the where the tube-side flow from one or several banks is col-
turbulent kinetic energy k is calculated from the standard lected, and probably mixed and distributed to other banks.
equation [27]; and the characteristic length-scale l is set The ‘‘inlet’’ manifold is a special type of manifold, which is
within the tube bundle as l ¼ P  Do , with P being the used to distribute to the appropriate tube banks the tube-
average local pitch. side fluid arriving from other plant devices not included
Outside the tube bundle, the effective viscosity is calcu- in the simulation (e.g., the steam fed to the reheater from
lated as usual with the Prandtl–Kolmogorov formula and the turbine). Several inlet manifolds are allowed in the
the k– model. model, so that several independent (not interconnected)
The effective diffusivity Cs in Eq. (3) is obtained as steam and water circuits can be present in the calculation,
Cs = lT/Prs. The turbulent Prandtl number for the shell- as is the case in power-generation boilers.
side enthalpy is taken as Prs = 0.9.
5. Solution
3.4. Steam properties
5.1. Simultaneous treatment of the three phases
The tube-side of the convective zone often consists of
steam in greatly-differing thermodynamic conditions. Thus, Eqs. (1)–(5) can be discretised and solved using any
properties such as density, thermal conductivity or specific CFD technique. The solution procedure is however compli-
 A. Gómez et al. / Applied Thermal Engineering 28 (2008) 532–546 537

cated by the presence of three distinct phases: the shell-side

gas-phase; the tube-side vapor- (or liquid-) phase; and the
metal, solid phase. An economical and easy-to-implement
method can be used when quasi-steady-state is assumed
in the metal-side enthalpy equation (5). In these circum-
stances, the left-hand side of the equation vanishes, and
the equation can be used to calculate the metal
temperature:
as!w T s þ at!w T t
Tw ¼ ð15Þ
as!w þ at!w
This value for Tw can then be used in Eqs. (3) and (4).
A further simplification in the numerical solution can be
achieved if the tube-side flow in the heat-exchanging banks
is always considered to be one-directional, and aligned with
any of the mesh directions. Then Eq. (4) for ht can be dis-
cretised as
VP
ðq htP  qtT htT Þ  max½At utw qtW htW ; 0
Dt tP
Fig. 2. Schematic of the two-dimensional heat exchanger, showing the
þ max½At utw qtP htP ; 0 þ max½At ute qtE htE ; 0 main dimensions and parameters.
at!w
 max½At ute qtP htP ; 0 ¼ V P ðhwtP  htP Þ ð16Þ
Cpt
solved alongside the discretised versions of the shell-side
where it has been assumed (without loss of generality) that equations (1)–(3), which in general of course do have a con-
the tube is oriented in the x (West-to-East) direction (see vection term.
Fig. 1); the max function expresses the upwinding practice
for the enclosed quantity (i.e., the value convected into the 5.2. Tube-side mass-conservation
cell through the cell face is the value in the upwind cell cen-
ter); At is the total cross-sectional area of the tubes in the The model conservation equations do not include a
cell; g = AtutPqtW and similar terms are the tube-side mass three-dimensional, tube-side continuity-equation because
flow rate per cell, constant along the cells in the same tube the flow in each tube is one-directional and confined. The
bank on mass-conservation grounds; and VP is the cell vol- geometrical model thus ensures mass conservation within
ume. Compared to the differential equation (4), the right- the same tube bank, and also when a bank is directly linked
hand side uses enthalpies, rather than temperatures. Thus, to another one. However, some provisions must be made to
ht = CptTt is the tube-side-fluid enthalpy, and hwt ¼ Cpt T w calculate the flow distribution when a manifold is fed by,
is the tube-side-fluid enthalpy at the wall temperature. and in turn feeds, several tube banks.
Since the tube-side mass flow rate g is known for each The mass flow rate mi through element i (be it a tube
cell in the bank from the geometrical model, Eq. (16) can bank or a manifold) is given by the equation
be cast into simple, source-only, linear equation in the X dij
form: mi ¼ m j þ bi ð18Þ
j6¼i
nj
0 ¼ S 1  S 2 htP ð17Þ
where mj is the mass flow rate through element j, dij is 1 if
with S1 and S2 being appropriate coefficients, variable from
element i feeds element j and 0 otherwise, nj is the number
cell to cell, but not explicitly dependent on htP. After this
of (equal) elements fed from j, and bi is the mass flow rate if
manipulation, the convective term in Eq. (16) has been
i is an inlet manifold, and 0 otherwise.
transferred to the source term in (17); and, since the latter
Eq. (18) is cast into a matrix form [A][m] = [b], which is
does not (formally) have a convection term, it can be easily
solved by LU decomposition to calculate the vector of the
unknown mass flow rates [m].

5.3. Numerical details

The equations and algorithms described above have

been coded by the authors into the finite-volume CFD-
code PHOENICS ([29]); version 2.2 has been used. For
Fig. 1. Three-cell schematic showing the discretisation of the tube-side accuracy and efficiency, structured, Cartesian grids are
enthalpy equation (for tubes aligned with the x direction). employed, and porosities are selectively used were needed
538 A. Gómez et al. / Applied Thermal Engineering 28 (2008) 532–546

Table 1
Validation cases: all cases solved with LMTD method, except case 5 with effectiveness-NTU
Case L (m) NT (–) NL (–) PT (m) PL (m) D (m) V1 (m/s) T in
s (°C) T out
s (°C) T in
t (°C) T out
t (°C)
1 0.75 20 15 0.038 0.038 0.019 8 325 341.95 375 375
2 1.0 12 8 0.0285 0.057 0.019 7 293.15 303.8 328.15 328.15
3 1.0 14 14 0.015 0.015 0.01 5 298.15 344.21 373.15 373.15
4 1.0 5 20 0.02 0.02 0.01 10 1200 615 400 400
5 4.0 10 10 0.05 0.05 0.025 5 800 625 300 345.04

Table 2
Validation cases: comparison of theoretical and simulation results
Case T out
s (K) T out
t (K) W (kW) a00s!w (W/m2 K)
Theoretical Simulation Theoretical Simulation Theoretical Simulation Theoretical Simulation
1 341.95 341.21 Constant Constant 91.80 91.64 167.0 168.5
2 303.80 303.38 Constant Constant 33.41 33.98 199 203.6
3 344.21 341.88 Constant Constant 62.39 62.26 210 211.9
4 615 641 Constant Constant 201.1 241.3 148.5 150.8
5 625.0 634.4 345.0 343.5 885.0 910.8 81.0 83.3

Fig. 3. Schematic of the simulated convective zone. Dotted lines indicate the circulation of the tube-side fluid between heat-exchanging elements. Triangles
are inlet manifolds. Circles are intermediate manifolds. See text for nomenclature.
 A. Gómez et al. / Applied Thermal Engineering 28 (2008) 532–546 539

to account of the fraction of local space occupied by tubes

and thus unavailable to the flue gas. Second-order interpo-
lation schemes are used for shell-side velocities, turbulence
statistics and enthalpy.
A typical discretisation for a convective zone, such as
that shown in Results section below, involved around
700,000 cells, and takes around 10 h CPU time to converge
on a Pentium IV at 3 GHz. A parametric analysis was per-
formed with a finer grid of 1.1 million cells to confirm that
results are grid independent.

6. Results

6.1. Validation tests

Validation tests have been performed mainly in respect

of two different aspects of the model. The first aspect is
the correct working of the several features of the geometri-
cal model, such as the different connectivity options and
the conservation of overall tube-side-fluid mass. These tests
will not be reported here.
The second validation initiative has consisted in predict-
ing of the overall heat-transfer rates in geometrically-sim-
ple configurations. The results of this validation exercise
are reported in this section.
The simple configurations considered are single-pass,
cross-flow heat exchangers, which can be analysed alterna-
tively using either the log mean temperature difference
(LMTD) method or the effectiveness-number of transfer
units (NTU) method [30]. While these alternative theoreti-
cal methods are also approximations, the similarity of the
results with both the numerical and theoretical approaches
adds confidence in the accuracy of the former.
Fig. 4. Geometry of the convective zone, as represented in the model.

Table 3 A general schematic of the configuration is represented

Operating conditions for different loads (mfr = mass flow rate)
(in cross-section) in Fig. 2, where the main heat exchanger
Load 60% 75% 85% 95% 100% and flow parameters are shown. Table 1 presents the values
Flue-gas mfr (kg/s) 256 313 360 402 440 of these parameters for the different cases used in the vali-
Flue-gas inlet-temperature 1300 1300 1300 1300 1300 dation. The test cases include a range of tube diameters,
(K)
pitch values, gas velocities and gas and tube-side inlet tem-
Water mfr at economiser 201 240 275 292 329
inlet (kg/s) peratures. Cases with constant and variable tube-side tem-
Water temperature at 503 516 523 526 542 perature are both considered.
economiser inlet (K) Validation results are presented in Table 2, and show a
Water pressure at 166 170 173 179 183 good agreement between both calculation methods. For
economiser inlet (bar)
instance, discrepancies in outlet temperatures in both cal-
Steam mfr at steam-drum 148.2 193.6 230 268 284
outlet (kg/s) culation methods are normally around 1%, with only one
Steam temperature at steam- 623 626 626 626 630 case (case 4) yielding differences around 4%.
drum outlet (K)
Steam pressure at steam- 165 171 172 178 180
drum outlet (bar)
6.2. Convective-zone modelling
Steam mfr at reheater inlet 142 180.2 208 232 266
(kg/s)
Steam temperature at 595 612 620 627 640 The model has been applied to the simulation of the
reheater inlet (K) convective zone of a power-station boiler. This is one of
Steam pressure at reheater 24 30 34.5 40.5 42.2 the three 350 MW(e), wall-fired boilers of the Teruel power
inlet (bar)
station located at Andorra (Teruel, Spain). The unit is a
540 A. Gómez et al. / Applied Thermal Engineering 28 (2008) 532–546

sub-critical coal-fired utility boiler, with natural circulation bine through the reheater (RH) and sent back to the power
and balanced draft. The boiler possesses three stages of cycle. The other one is the main water–steam circuit. The
superheat, a single reheat a two economisers. The combus- water is fed into this circuit through the lower economiser
tion and heat transfer in the furnace zone has been already (LE) and then through the upper economiser (UE). From
investigated numerically in the past by one of the present there, the hot water is sent to the drum separator, where
authors [31]. liquid water and steam are at equilibrium conditions. The
A schematic of the convective zone is shown in Fig. 3. liquid phase from the drum circulates through down-com-
After leaving the furnace, the flue gas flows along the radi- ers to the furnace water-walls where the phase change from
ant superheater (wing walls) and final superheater (FS). liquid to steam takes place using the heat released in the
Downstream of the final superheater, the gas enters the combustion zone, and then returns as a two-phase mixture
convective zone separated in two parallel passes. Two inde- to the drum. From there, dry steam is distributed through
pendent water–steam circuits can be found in the boiler. In the tubes that line the ceiling and walls of the convective
the shorter one, steam is drawn from the high-pressure tur- zone (which therefore serve as a preliminary superheater),

1199
1203 1201
1231

1189
1199

84
11
10
09

1203
1258

870

1120
828
939

74

731
4

68
9
0
72

773 70
3
731
759

73
1

615
62
634

0
648

5
62

Fig. 5. Shell-side temperature and shell-side flow-field on a middle plane for 60% load.
 A. Gómez et al. / Applied Thermal Engineering 28 (2008) 532–546 541

and then is collected in a manifold and circulated succes- to regulate the gas temperature at each of air preheater
sively through the primary superheater (PS), the wing walls inlets.
(WW) and the final superheater (FS). From this last heat Some plant measurements in the convective zone are
exchanger, the steam is sent to the high-pressure turbine. available, and will be used for comparison with computed
The gas flow-rates through the parallel passes of the con- results in this section. This type of utility boilers are con-
vective section are adjusted to maintain the steam temper- ventionally instrumented, providing of course enough mea-
ature at reheater outlet by means of regulating dampers, surements for their safe operation and control. However,
located under the lower economiser. The reheater is located such information is neither complete nor reliable enough
in the rear pass of the boiler, whereas the primary super- to be used in the validation of thermo-fluid-dynamical
heater and upper economiser are located in the frontal models. In part to remedy this situation, a plan for the
one. Aside from the main gas outlet at the cold end of improvement and extension of instruments and measure-
the lower economiser, there is a secondary extraction of ment procedures was designed and implemented in this
hot gases (between the reheater and the lower economiser) case-study boiler. The rationale and details of this plan

1217

1201
1189

1106
120
9

1199

10
22

1258
911
1134
20
11

898

12
64
983

80

814
0

759

793
3
828

77
81

81
4

79
3
683

68
731

9
70

703
3

Fig. 6. Shell-side temperature and shell-side flow-field on a middle plane for 100% load.
542 A. Gómez et al. / Applied Thermal Engineering 28 (2008) 532–546

can be found in [32]. As for the variables used for valida- 640 K. Table 3 presents the operating conditions for differ-
tion in the present paper, discussed in this section, a grid ent loads which are used as input data in the simulations
of eight gas thermocouples was installed at the lower carried out to validate the model. Most of these data are
economiser outlet, and six water and steam thermocouples directly-available measurements, but some of them are
were installed at upper economiser, primary superheater indirectly estimated using other operating data. This is
and reheater outlets; the former were installed to reduce the case of gas mass-flow-rate (obtained by closing the
the uncertainty due to temperature stratification in the flow overall mass balance in the boiler), steam mass-flow-rate
and the latter to obtain more reliable data by taking addi- leaving the drum (mass balance in the main steam circuit)
tional measurements at ends of the external headers. Seri- and steam mass-flow-rate at reheater inlet (mass balance
ous technical problems arose when attempting to considering the high-pressure extractions to the water heat-
instrument headers located inside the boiler enclosure. ers in the cycle).
Thermowells equipped with armoured outfit can be For these convective-zone calculations, tube-side heat-
installed in this situation, but this possibility was discarded transfer resistances are neglected. For the shell side, values
owing to doubts about the integrity of the pressure circuit. calculated from the application of semi-empirical heat-
Therefore, the water temperature at lower economiser out- transfer models are used. The direct measurement of ash
let (i.e., upper economiser inlet) and the steam temperature fouling rates in this type of equipment is technically unfea-
at primary superheater inlet were not measured. sible, and therefore one has to resort to indirect thermal
The operating conditions change considerably depend- calculations. In this case, a semi-empirical approach was
ing on the load, or fraction of the nominal power at which proposed in [33], aiming at the on-line estimation of these
the boiler operates. For instance, the pressure at econo- thermal resistances. The method basically consists in the
miser inlet varies from 166 bar to 183 bar, and the range formulation of a lumped model for large heat exchangers,
of temperatures at the reheater inlet from 595 K to such as those found in utility boilers, coupled with energy

Fig. 7. Tube-side temperature on a middle plane for 60% load.

 A. Gómez et al. / Applied Thermal Engineering 28 (2008) 532–546 543

balances in a separated section-by-section analysis. Typical dimensional view of the mid-plane for 60% and 100% load,
mean values for the shell-side heat-transfer resistances, cal- respectively, showing how the gas temperature gradually
culated in this way to account for ash fouling, have been decreases along the flow path. The flow pattern is very sim-
used in the CFD simulations: for the re-heater Rfs = ilar in both cases; and the temperature map is also similar,
0.00459 K m2/W; for the primary superheater Rfs = albeit the temperatures are lower for the lower load as
0.004163 K m2/W; for the lower economiser Rfs = 0.00173 expected.
K m2/W; for the upper economiser Rfs = 0.002062 K m2/W. Figs. 7 and 8 are the corresponding tube-side tempera-
Fig. 4 is a three-dimensional view of the geometry as tures. The interconnection of the tube banks in both steam
generated directly from the geometrical model (tube banks circuits is apparent from the gradual change in the tube-
are represented as solid boxes for simplicity, but fluid flow side temperature even when the interconnected banks are
is of course allowed through them in the model). The sim- not contiguous in physical space.
ulation domain starts at the end of the combustion zone, The graph in Fig. 9 compares the total heat exchanged
just under the ‘‘nose’’, where uniform velocity and temper- in the convective zone (reheater, primary superheater,
ature profiles are presumed. (The convective-zone model lower and upper economiser) with the measured one for
can be coupled with models for the combustion zone, several loads. The agreement is generally good, but the case
which would allow this presumption to be relaxed by using for 100% load shows a more significant discrepancy, with
calculated results in the absence of detailed experimental the model tending to follow more closely the linear behav-
data.) iour found for lower loads. Fig. 10 is the corresponding
Figs. 5 and 6 represent the shell-side temperature and graph for the calculated and measured average outlet
(superimposed) the shell-side velocity vectors in a two- temperatures of water or steam in the upper economiser,

Fig. 8. Tube-side temperature on a middle plane for 100% load.

544 A. Gómez et al. / Applied Thermal Engineering 28 (2008) 532–546

Fig. 9. Comparison between calculated and measured exchange heat for different loads.

Fig. 10. Comparison between calculated and measured mean temperatures at the outlets of RH, PS, UE and the gas exit.

reheater and primary superheater outlets, and the com- isfactory. These results are indicative of the capability of
puted and measured mean temperature of the gas at the the model to deal with the different operating conditions
outlet of the convective zone. There is a good agreement of a real plant and its usefulness in the study of the perfor-
in both figures between calculated and experimental data. mance of convective zones.
The model described here can be extended with relative
7. Conclusions and future work ease to account for additional physical processes. Firstly,
heat transfer by radiation is bound to relevant in the con-
The present paper has introduced a means of simulta- vective surfaces which are exposed to the flame. Heat trans-
neously calculating the shell-side flow-field, the shell-side fer by radiation can be accounted for in the model either as
and tube-side thermal fields and the tube-wall temperature a direct flame-to-surface phenomenon (e.g., by use of view
in the convective zone of a power-station boiler. The model factors), or by considering the gas as a participating media,
allows for several arbitrarily-interconnected heat-exchang- by use of a full radiation model. In the absence of other
ing elements to be simulated in a flexible manner. The data, incoming radiation from the flame zone to the con-
model has been validated with the simulations of a real vective one can be estimated from furnace-combustion
power-station convective zone for different loads, and the models. These can be usefully coupled to the present con-
agreement between calculated and plant data has been sat- vective-zone model to obtain other relevant boundary con-
 A. Gómez et al. / Applied Thermal Engineering 28 (2008) 532–546 545

ditions which are difficult to obtain experimentally (e.g., [10] M. Prithiviraj, M. Andrews, Three-dimensional numerical simulation
velocity or temperature distributions on the convective- of shell-and-tube heat exchangers. Part II: Heat transfer, Numerical
Heat Transfer, Part A 33 (8) (1998) 817–828.
zone inlet-plane). A two-way coupling would allow, among [11] J.Y. Tu, C.A.J. Fletcher, M. Behnia, Numerical modelling of three-
other features, to account for the possible effect of pressure dimensional fly-ash flow in power utility boilers, International
losses in the convective zone on the flow pattern in the Journal for Numerical Methods in Fluids 24 (1997) 787–807.
furnace. [12] P.J. Coelho, Mathematical modeling of the convection chamber of a
Shell-side fouling has been modelled in the present work utility boiler – an application, Numerical Heat Transfer, Part A 36
(1999) 411–428.
using experimental values for the associated heat-transfer [13] L.I. Diez, C. Cortes, A. Campo, Modelling of pulverized coal boilers:
resistance. Research work, such as that cited in the intro- review and validation of on-line simulation techniques, Applied
duction, can be with some effort integrated in the model Thermal Engineering 25 (2005) 1516–1533.
to simulate the motion of the fly-ash in the convective sec- [14] Z. Niu, K.-F.V. Wong, Adaptive simulation of boiler unit perfor-
tion, the build-up of deposits, and their impact on heat mance, Energy Conversion Management 39 (13) (1998) 1383–1394.
[15] L.Y. Huang, J.X. Wen, T.G. Karayiannis, R.D. Matthews, Numer-
transfer. Work along these lines has been recently pub- ical prediction of high efficiency boiler heat exchanger performance,
lished by Kær et al. [26]. Applied Thermal Engineering 18 (1998) 1089–1099.
In respect of the tube-side sub-models, an extension can [16] A. Horvat, M. Leskovar, B. Mavko, Comparison of heat transfer
be devised to calculate head losses in the tube-side circuits. conditions in tube bundle cross-flow for different tube shapes,
The flow distribution from the manifolds to the individual International Journal of Heat and Mass Transfer 49 (2006) 1027–
1038.
tubes can then be calculated, if required, in a more accurate [17] S.B. Beale, S.V. Zhubrin, A distributed resistance analogy for solid
manner using this information. oxide fuel cells, Numerical Heat Transfer, Part B 47 (2005) 573–591.
[18] D. Butterworth, The development of a model for three-dimensional
Acknowledgements flow in tube bundles, International Journal of Heat and Mass
Transfer 21 (1977) 253–256.
[19] ASHRAE Handbook, Fundamentals, American Society of Heating,
This work was partly supported by the European Union Refrigerating and Air Conditioning Engineers, Inc., Atlanta, GA,
under Contract ECSC 7220/ED/096 to Grupo ENDESA. 1989.
The cooperation of Dr. Juan Francisco Gonzalez and Mess [20] A. Zhukauskas, Heat transfer from tubes in crossflow, in: J.P.
Enrique Perez and Mariano Lacarta of Grupo ENDESA is Harnett, T.F. Irvine (Eds.), Advances in Heat Transfer, vol. 8,
gratefully acknowledged. We also thank CHAM (London, Academic Press, New York, 1972, pp. 93–160.
[21] T. Erickson, S. Allan, D. McCollor, J. Hurley, S. Srinivasachar, S.
UK), for granting us the use of their PHOENICS CFD Kang, J. Baker, M. Morgan, S. Johnson, R. Borio, Modelling of
code, in which the models outlined in this paper have been fouling and slagging in coal-fired utility boilers, Fuel Processing
coded by the authors. Technology 44 (1995) 155–171.
[22] S. Zibas, S. Idem, Boiler heat transfer modeling using CEMS data
with application to fouling analysis, in: ASME Joint Power Gener-
References ation Conference, vol. 30, 1996, pp. 655–661.
[23] G. Bergeles, D. Bouris, M. Yianneskis, S. Balabani, A. Kravaritis, S.
[1] S.V. Patankar, D.B. Spalding, A calculation procedure for the Itskos, Effects of fouling on the efficiency of heat exchangers in lignite
transient and steady-state behaviour of shell-and-tube heat exchang- utility boilers, Applied Thermal Engineering 17 (8–10) (1997) 739–
ers, in: A.A. Afgan, E.U. Schluner (Eds.), Heat Exchangers: Design 749.
and Theory Sourcebook, McGraw Hill, New York, 1974, pp. 155–176. [24] H. Wang, J. Harb, Modeling of ash deposition in large-scale
[2] D.B. Rhodes, L.N. Carlucci, Predicted and measured velocity combustion facilities burning pulverized coal, Programmes in Energy
distributions in a model heat exchanger, in: International Conference Combustion Science 23 (1997) 267–282.
on Numerical Methods in Engineering, Canadian Nuclear Society/ [25] J. Isdale, P. Mercier, J. Grillot, A. Mulholland, J. Gomatam,
American Nuclear Society, 1983, pp. 935–948. Integrated modelling of process heat transfer with combustion and
[3] V.M. Theodossiou, A. Sousa, L. Carlucci, Flow field predictions in a fouling, Applied Thermal Engineering 17 (8–10) (1996) 751–762.
model heat exchanger, Computational Mechanics 3 (1988) 419–428. [26] S. Kær, L. Rosendahl, L. Baxter, Towards a CFD-based mechanistic
[4] C. Zhang, A.C.M. Sousa, J.E.S. Venart, The numerical and exper- deposit formation model for straw-fired boilers, Fuel 85 (2006) 833–
imental study of a power plant condenser, Journal of Heat Transfer 848.
115 (1993) 435–445. [27] B.E. Launder, D.B. Spalding, Mathematical Models of Turbulence,
[5] C. Zhang, Y. Zhang, A quasi-three-dimensional approach to predict Academic Press, New York, 1972.
the performance of steam surface condensers, Journal of Energy [28] W. Wagner, J. Cooper, A. Dittmann, J. Kijima, H.-J. Kretzschmar,
Resources Technology 115 (1993) 213–220. A. Kruse, R. Mareš, K. Oguchi, H. Sato, I. Stocker, O. Šifner, Y.
[6] A. Zhang, C. Bokil, A quasi-three-dimensional approach to simulate Takaishi, I. Tanishita, J. Trübenbach, T. Willkommen, The IAPWS
the two-phase fluid flow and heat transfer in condensers, Interna- Industrial Formulation 1997 for the thermodynamic properties of
tional Journal of Heat Mass Transfer 40 (15) (1997) 3537–3546. water and steam, Journal of Engineering for Gas Turbines and Power
[7] A. Zhang, C. Bokil, A quasi-three-dimensional approach to simulate 122 (2000) 150–182.
the two-phase fluid flow and heat transfer in condensers, Interna- [29] D.B. Spalding, A general purpose computer program for multi-
tional Journal of Heat Mass Transfer 40 (15) (1997) 3537–3546. dimensional one and two phase flow, Journal of Mathematics and
[8] S.J. Ormiston, G.D. Raithby, L.N. Carlucci, Numerical modeling of Computers in Simulation 23 (1981) 265–276.
power station steam condensers. Part 2: Improvement of solution [30] F.P. Incropera, D.P.D. Witt, Fundamentals of Heat and Mass
behavior, Numerical Heat Transfer, Part B 27 (1) (1995) 103–125. Transfer, second ed., John Wiley and Sons, New York, 1985.
[9] M. Prithiviraj, M. Andrews, Three-dimensional numerical simulation [31] N. Fueyo, J. Ballester, C. Dopazo, An Eulerian–Eulerian model of
of shell-and-tube heat exchangers. Part I: Foundation and fluid coal combustion, NOx formation and reburning, in: 12th Annual
mechanics, Numerical Heat Transfer, Part A 33 (8) (1998) 799–816. International Pittsburgh Coal Conference, 1995, pp. 1113–1116.
546 A. Gómez et al. / Applied Thermal Engineering 28 (2008) 532–546

[32] L.I. Dı́ez, C. Cortés, I. Arauzo, A. Valero, Combustion and heat [33] C. Cortés, L.I. Dı́ez, A. Campo, Modeling large-size boilers as a set of
transfer monitoring in large utility boilers, International Journal of heat exchangers: tips and tricks, ASME HTD Combustion and
Thermal Science 40 (2001) 489–496. Energy Systems 369 (4) (2001) 41–48.