International Journal of Thermal Sciences 48 (2009) 1184–1195

Contents lists available at ScienceDirect

International Journal of Thermal Sciences

www.elsevier.com/locate/ijts

Optimal design of a plate heat exchanger with undulated surfaces

A.G. Kanaris, A.A. Mouza, S.V. Paras ∗
Laboratory of Chemical Process and Plant Design, Department of Chemical Engineering, Aristotle University of Thessaloniki, Univ. Box 455, GR 54124, Thessaloniki, Greece

a r t i c l e i n f o a b s t r a c t

Article history: The purpose of this study is to suggest a general method for the optimal design of a plate heat exchanger
Received 17 January 2008 (PHE) with undulated surfaces that complies with the principles of sustainability. A previously validated
Received in revised form 12 September CFD code is employed to predict the heat transfer rate and pressure drop in this type of equipment.
2008
The computational model is a three-dimensional narrow channel with angled triangular undulations in a
Accepted 5 November 2008
Available online 28 November 2008
herringbone pattern, whose blockage ratio, channel aspect ratio, corrugation aspect ratio, angle of attack
and Reynolds number are used as design variables. To limit the number of simulations needed, the Box–
Keywords: Behnken technique is employed. An objective function that linearly combines heat transfer augmentation
Plate heat exchanger with friction losses, using a weighting factor that accounts for the cost of energy, is employed for the
Sustainable design optimization procedure using response surface methodology (RSM). New correlations are provided for
CFD predicting Nusselt number and friction factor in such PHEs. The results are in very good agreement with
Optimization published data. Finally, optimal design specifications are suggested for a range of Re for two values of the
Response surface methodology
weighting factor.
 2008 Elsevier Masson SAS. All rights reserved.

1. Introduction to the flow and increases the friction losses. It is therefore appar-
ent that, in order to come up with an optimal design, a ‘compro-
The need for designing process equipment that complies with mise’ between heat transfer and pressure drop must be set. As the
the principles of economic and ecological sustainability (sustain- plate patterns greatly influence both hydraulic and thermal perfor-
able design) acted as a driving force towards the evolution in the mance, the final design is certainly dependent on the initial choice
design of plate heat exchangers (PHE). The advantages of a PHE of the plate pattern.
over conventional heat transfer equipment also increased their ac- The majority of the performance data on this type of equipment
ceptance in two-phase flow applications. They are commonly used are considered proprietary, in a highly competitive PHE market.
in process and power industries for a wide range of temperatures Thus, the lack of data in open literature has held back the devel-
due to their compactness, close temperature approach and ease on opment of a more generic design model for the corrugated plates
inspection and cleaning [1,2]. The plates of these heat exchangers of a PHE, even though there are studies referring to the effect
comprise some form of near-sinusoidal corrugations in a herring- of several geometrical parameters on the heat transfer coefficient
bone (or chevron) pattern, a design commonly used for PHEs as and the friction factor. For instance, Focke et al. [5] and Heavner
it is considered the most successful type [3]. It is known [4] that et al. [6] provide empirical equations based on experimental data
there are two mechanisms for the augmentation of heat transfer, regarding the effect of the corrugation inclination angle on the per-
which are also accompanied by pressure loss increase; inducing formance of an industrial PHE. Martin [3], in a theoretical approach
flow separation and reattachment, and increasing the wall skin using the generalized Lévêque equation, provide correlations based
friction. When two of these plates are arranged and placed abut- on chevron angle for predicting friction factor and Nusselt num-
ting, a channel with complicated passages is formed. Due to the
ber in typical PHEs, with limited accuracy (±40%) [4]. Numerical
breakup and reattachment of boundary layers, secondary flows and
work has also been performed; either on representative elements
small hydraulic diameter of the flow passages, high heat transfer
or in complete conduits. For example, Mehrabian and Poulter [7]
coefficients are achieved. This allows a small requirement in sur-
study the flow inside a furrow of a PHE with sinusoidal corru-
face area, up to one third that of a shell-and-tube heat exchanger
gations, while Hossain and Sadrul Islam [8] study similar control
for a given duty [4], thus reducing the cost, overall volume and
volumes for three different corrugation shapes. Moreover, Asako et
space requirement for the exchanger. Nevertheless, the series of
al. [9] study the case of rounding the corrugations in wavy pas-
periodic changes in flow direction induces a significant resistance
sages with constant plate spacing using a numerical approach for
low Reynolds numbers.
* Corresponding author. There is an increasing interest on the study of various geomet-
E-mail address: [email protected] (S.V. Paras). rical parameters of the modulated pattern of the plates with the

1290-0729/$ – see front matter  2008 Elsevier Masson SAS. All rights reserved.
doi:10.1016/j.ijthermalsci.2008.11.001
 A.G. Kanaris et al. / International Journal of Thermal Sciences 48 (2009) 1184–1195 1185

Nomenclature

A PHE channel area . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m2 y+ y-plus

Af mean cross flow section . . . . . . . . . . . . . . . . . . . . . . . . . . . . m2 z width of corrugations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m
BR blockage ratio, d/ H ∆ T ln logarithmic temperature difference . . . . . . . . . . . . . . . . . . K
ChanAR channel aspect ratio, H / W ∆P pressure drop . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Pa
CorAR corrugation aspect ratio, d/ z
cp specific heat capacity . . . . . . . . . . . . . . . . . . . . . . . J kg−1 K−1 Greek
d height of corrugations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m β weighting factor of the objective function
Dh hydraulic diameter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m βc weighting factor correction accounting for corrugation
F objective function shape
f friction factor η function
fc friction factor correction accounting for corrugation
θ angle of attack . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . deg
shape
µ dynamic viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Pa s
H half of the channel height . . . . . . . . . . . . . . . . . . . . . . . . . . . m
ρ density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . kg m−3
h heat transfer coefficient . . . . . . . . . . . . . . . . . J m−2 K−1 s−1
Φ area enlargement factor
k thermal conductivity . . . . . . . . . . . . . . . . . . . . J m−1 K−1 s−1
L length of the plate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m Subscripts
N number of flow channels per stream
c cold-fluid stream
Nu overall Nusselt number, based on hydraulic diameter
f friction-related term
P wetted perimeter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m
Pr Prandtl number h hot-fluid stream
Re Reynolds number based on hydraulic diameter Nu thermal-related term
T temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K 0 smooth surface-based term
u mean entrance velocity . . . . . . . . . . . . . . . . . . . . . . . . . . m s−1 i inlet
W half of the channel width . . . . . . . . . . . . . . . . . . . . . . . . . . . m o outlet

intention to propose an optimum design of this type of equipment. mate goal is to suggest general correlations for the design of a PHE
Kim and Kim [10] seek the optimum design in the case of a rectan- with corrugations in a herringbone pattern.
gular rib-roughened channel placed vertically to the flow direction,
while Kim and Lee [11] also use rectangular ribs, but having them 2. ‘Experimental’ procedure
placed in a crossed pattern. Zhang et al. [12] check the effect of the
fin waviness and spacing for low Reynolds numbers and suggest an 2.1. Selection of design parameters
optimal design based on these parameters. Wang and Sundén [13]
A series of single three-dimensional channels with triangular
suggest an iterative procedure for determining the optimal design
ribs in a V-shaped (i.e., herringbone) pattern on both of channel
of PHEs with and without pressure drop limitations, taking into ac-
sides is considered for this study. The plates are identical and are
count the effect of the corrugation angle. Nevertheless, to the best
placed in a way that a crossed pattern is created, as presented in
of author’s knowledge, a correlation combining all the design pa-
Fig. 1.
rameters of a PHE is not available in the open literature.
The main geometrical parameters (Fig. 2) which, as expected,
The type of flow inside the passages of a PHE is a significant
influence the hydrodynamic behavior of the system are:
issue. Focke and Knibbe [14], who performed flow visualization
experiments in narrow passages simulating a corrugated PHE, sug-
• the height (d) and width (z) of the corrugation,
gest that local flow structure controls the heat transfer process,
• the height (2H ) and width (2W ) of the channel, and
due to the existence of secondary swirling flow patterns, known
• the angle of attack, theta (θ ).
as Goertler vortices. The vortices and secondary flows increase tur-
bulence production and help convective processes for heat trans-
As the corrugation pattern is continuous, (i.e., no gap between the
fer augmentation by increasing advection of fluid from the center triangular corrugations exists) the width of the corrugation is also
towards the side walls and by producing high shear [4]. Hessel- the pitch of the repeated modulation pattern.
graves [15] also mention that in plain ducts with polygonal ge- To perform the simulations, five dimensionless groups are se-
ometries the transition between laminar and turbulent flows is lected as design variables, namely:
rather sharp. Moreover, Shah and Wanniarachchi [1] acknowledge
that there is evidence of turbulent flow even for the low Reynolds • the blockage ratio (BR), defined as d/ H , which expresses the
number range of 100–1500. The latter statement is also supported percentage of the entrance of the channel that is ‘blocked’
by Vlasogiannis et al. [16], who verify that the flow inside an in- with corrugations and is also referred to in literature as rib-
dustrial PHE is turbulent for Re > 650. to-channel height ratio [11];
Previous work conducted in this Laboratory [17,18] have proven • the channel aspect ratio (ChanAR), defined as H / W , which is
that CFD is a reliable tool for simulating the operation of PHE. a measure of how narrow the channel is;
Thus, instead of expensive and time consuming laboratory experi- • the corrugation aspect ratio (CorAR), defined as d/ z, which
ments, CFD simulations can be used for predicting the performance expresses the obtuseness or acuteness of the triangular cor-
of this type of equipment. The purpose of the present work is to rugation and is also directly related to the area enlargement
employ numerical simulation as an ‘experimental’ tool to exam- factor [13];
ine the combined effect of most of the geometrical parameters on • the sine of twice the angle of attack (sin 2θ ) which is com-
the operation of PHE and to check in what extent each of these monly used by researchers, i.e., Martin [3], in correlations con-
parameters affects the momentum and heat transfer rate. The ulti- cerning heat transfer and friction losses, and
1186 A.G. Kanaris et al. / International Journal of Thermal Sciences 48 (2009) 1184–1195

Fig. 1. Computational domain used for the simulations.

Table 1
Constrains of the design variables.

Parameter Lower limit Upper limit

BR 0.1 0.9
ChanAR 0.05 0.25
CorAR 0.05 0.40
θ , deg 25 75
Re 500 6000

tions, the hydraulic diameter can be calculated using the geomet-

rical design parameters, by the following equation:
 
2d 2
Dh = −1 (3)
ChanAR + Φ BR
where Φ is the plate area enlargement factor, defined for triangu-
lar corrugations (as in the present study) as:
Fig. 2. Geometric parameters of the computational domain. 0.5
Φ = 1 + 4(CorAR)2

(4)

• the Reynolds number (Re) defined as:

2.2. Assumptions and simplifications
u Dh ρ
Re = (1) The design specifications, which are commonly encountered in
µ
industrial PHEs, are used to select the bounds of the design vari-
where u is the mean entrance velocity, ρ and µ the density and ables values in this study, as presented in Table 1. The upper limit
of 0.9 on blockage ratio (BR) is imposed by restrictions when cre-
dynamic viscosity of the fluid, respectively, and D h is the hydraulic
ating the grid, as the small gap between the plates stiffens the
diameter of the conduit, defined as:
creation of appropriate grid cells in this area. As the height of the
corrugations in a typical PHE is approximately 2 mm, the space
4A f between the plates of a channel with BR = 0.9 is about 0.2 mm.
Dh = (2)
P Therefore, it can be assumed that for BR = 0.9, the channel plates
are nearly in contact. In the case of corrugation aspect ratio, most
where A f is the mean cross sectional area and P is the wetted PHEs do not typically implement an area enlargement factor, Φ ,
perimeter of the cross-section. In the case of triangular corruga- greater than 1.25 [19], which, according to Eq. (4), equals to a value
 A.G. Kanaris et al. / International Journal of Thermal Sciences 48 (2009) 1184–1195 1187

of 0.375 for CorAR. An upper limit of 0.40 for CorAR is set for this ements) is imposed in the vicinity of the walls. This procedure is
study. necessary for confined geometries [20] in order to create elements
It is a known fact that pressure drop in a PHE is very high with y + smaller than 1.0 near the corrugated wall, as imposed
compared to a shell-and-tube heat exchanger [1]. In the relevant by the SST model. As the flow is symmetrical with respect to
work of Wang and Sundén [13], allowable pressure drop does not the yz-plane and to keep the memory requirements low, the final
exceed 120 kPa. Preliminary simulation results also prove that the computational model represents half the actual channel (Fig. 1).
pressure drop inside the PHE channel increases above this value Steady-state condition is set for all simulations. Fig. 3 presents typ-
for the majority of the geometrical design-variable combinations ical results of temperature and pressure drop distribution in the
as Re increases, e.g., up to 350 kPa for Re = 7000. To keep the channel. The results are finally expressed in terms of local and
pressure drop within acceptable values, as suggested by Shah and overall Nusselt number and friction factor.
Wanniarachchi [1], the Re of the flow inside a single PHE channel
is set between 500 and 6000. 3. Quest for the optimal design
For the chevron-type PHE, the shape of the corrugation can vary
between the sharp shape of the triangular-type and the smooth 3.1. Optimization method
shape of the sinusoidal-type modulation, depending on the shape
of the fabricating tool used to create the corrugations during the The optimal solution for a V-shaped PHE is obtained using a
manufacturing process. The inclusion of the shape of corrugations polynomial-based Response Surface Method (RSM), a global op-
as a design variable would further increase the complexity of the timization method that includes a collection of techniques, i.e.,
optimization problem (in terms of computational demand and op- design of experiments (DOE), regression analysis and analysis of
timization method). Kanaris et al. [18] have proven that it is only variance (ANOVA) techniques [22]. The DOE techniques allow the
the pressure drop that is affected by the shape of the corrugations designer to extract as much information as possible from a limited
and not the heat transfer rate. Thus, the shape of the corrugations number of test cases. For example, for the case of five design fac-
is selected to be triangular for simplicity in constructing the grid tors, a three-level full factorial design would require 243 design (or
for each computational case, while the effect of the corrugation training) points compared to the 41 design points required by the
shape on the friction factor will be discussed later in this study. Box–Behnken, a common DOE technique. This method is ideal for
The working fluid is water and enters the conduit on the z-di- CFD cases, where a single simulation can take hours or even days
rection (Fig. 1), having an entrance temperature of 40 ◦ C in all sim- to run.
ulations. The temperature of the corrugated surfaces is set to a Using a DOE technique, a series of experiments or numerical
constant value of 20 ◦ C, while all other walls are considered adi- simulations are performed for a prescribed set of design points,
abatic. Although the simplification of constant wall temperature in order to construct a response surface of the measured quan-
seems unrealistic, it allows the designer to focus on the PHE opti- tity over the design space, i.e., a function η(x1 , x2 , . . . , xn ), where
mization without simulating a second channel for the cold stream, x1 , x2 , . . . , xn are the design variables with values constrained in
which would lead to an enormous increase in computational de- specific ranges. When a second-order polynomial is used, the re-
mands [18]. A commercial CFD code, namely ANSYS CFX® 10.0, is sponse surface equation is expressed as:
employed to simulate the flow and to obtain results concerning n
 n
 n

heat transfer characteristics and pressure drop inside the channel. η = a0 + a jxj + a j j x2j + a i j xi x j (5)
In the present calculations, the CFD code uses a high resolution j =1 j =1 i = j
advection scheme for the discretization of the momentum equa-
tions, while for the pressure-velocity decoupling, a method similar where αii are the unknown coefficients of the polynomial equation
to that by Rhie and Chow is applied [20]. Due to the increased determined using the calculated values of the measured quantity
computational demands and to the number of simulations needed, for the prescribed set of design points. This equation, η , is also re-
a high performance cluster (HPC) for parallel computing, consisting ferred to as a quadratic model. A set of training points is selected
of six 64-bit AMD processors with a total of 16 GB RAM, is used. using a Box–Behnken design and is presented in Table 2. The Box–
The HPC uses Gbit Ethernet connections and runs Gentoo Linux, Behnken design method is valid for use in RSM as it is proven [23]
while the commercial CFD code uses the MPICH protocol for mes- that it is sufficient to fit a quadratic model. Each independent vari-
saging between the nodes. The wall clock time needed for each able is placed at one of three equally spaced values – low, center
simulation run varied between 8 and 24 hours, depending on the and high value (as in a three-level full factorial design) – and the
mesh density of each computational model. ratio of the number of training points to the number of coefficients
Each computational model with different geometrical details of the quadratic model is kept between 1.5 and 2.6. A certain num-
is constructed using the parametric design features of ANSYS ber of factors are put through all combinations for the factorial
Workbench® 10.0. Relevant work in this Lab concerning PHE [17,18] design, while the other factors are kept to their center values. This
has proven that the SST turbulence model [21] is the most appro- type of DOE technique ensures that factors will not be set to their
priate for simulating the flow inside this type of conduit. More high levels simultaneously and thus avoids extreme factor combi-
specifically, the SST model activates k–ω model near the wall, nations [24].
where it performs best, while it switches to k–ε for the rest of The RSM has been successfully used in several attempts to op-
the flow. As the SST model is basically a two-equation turbulence timize a compact heat exchanger [10,11,25,26]. A common method
model, it provides the eddy diffusivity variable in order to express for testing the significance of a response surface model is the ad-
the turbulent fluctuation terms in the Reynolds-averaged transport justed coefficient of determination, R 2adj . It is suggested that, when
equations [20]. The quality of the mesh is examined and modified 0.9  R 2adj  1, the value of the objective function is accurately pre-
to speed up convergence and increase robustness, while maintain- dicted by the response surface model [10].
ing accuracy in the simulations.
The various grids used for the simulations are unstructured 3.2. Definition of the objective function
meshes consisting of tetrahedral and prism elements. The tetrahe-
dral elements are created by the ANSYS ICEM CFD® Mesher using To optimize the PHE performance, a function that incorporates
the Delaunay method [20]. In order to facilitate the boundary layer both the enhancement of heat transfer and the increase of fric-
calculations, a layer of prism elements (also known as ‘inflated’ el- tion losses must be formulated [27]. All the other cost factors, i.e.,
1188 A.G. Kanaris et al. / International Journal of Thermal Sciences 48 (2009) 1184–1195

Fig. 3. Typical pressure and temperature distributions across a corrugated plate.

64
equipment, maintenance and labor cost are not practically affected f0 = , Re < 2300 (11)
by the geometry of the PHE. Consequently, a thermal-related func- Re
tion (ηNu ) that gives the heat transfer augmentation with respect
f 0 = (1.8 log Re − 1.5)−2 , Re  2300 (12)
to that for a smooth surface, is defined as:
 −1
Nu Each of the functions is expressed as a quadratic model of the
ηNu = (6)
Nu0 five design variables (Table 1), using the RSM. Finally, the afore-
mentioned two functions are linearly combined in an objective
The overall Nusselt number Nu is defined as:
function, F , which compromises between the heat transfer aug-
hD h mentation and the inevitable increase of friction losses [27]:
Nu = (7)
k
where D h is the hydraulic diameter of the conduit, k the average F = ηNu + β η f (13)
thermal conductivity and h the average heat transfer coefficient
predicted by the simulation. Nu0 is the Nusselt number for fully- where the parameter β is a weighting factor that accounts for the
developed flow in a smooth pipe and is calculated by the Dittus– pumping cost to thermal energy cost and can be estimated using
Boelter correlation [15]: data concerning the cost of a unit of heat produced by a common
Nu0 = 0.023Re0.8 Pr0.3 (8) fuel (e.g., natural gas) and the cost of the same amount of elec-
tric energy. The definition of these two functions is suggested [28]
where Pr is the Prandtl number. Accordingly, a friction-related func- and used [10,11] in the literature for this optimization approach.
tion (η f ), which expresses the increase of friction losses with re- The linear combination of two factors concerning the heat transfer
spect to that for a smooth surface is defined as: augmentation and the corresponding friction losses is also men-
  1/ 3 tioned in the fundamental work of Bejan et al. [29]. According to
f
ηf = (9) literature, β varies between 0 and 0.1 [10]. For the purpose of this
f0
study, results are presented for various Re and for two values of the
where the friction factor, f , is given by: weighting factor, i.e., 0.03 and 0.06, where 0.03 corresponds to the
present day conditions in Greece, and the effect of the weighting
∆P Dh
f = (10) factor on the values of the optimal design parameters is discussed.
L ρ u 2 /2 Therefore, to maximize the performance of the type of PHEs
where ∆ P / L is predicted by the simulation. The friction factor for in question, the objective function defined above should be mini-
smooth circular tube, f 0 , is calculated by [15]: mized.
 A.G. Kanaris et al. / International Journal of Thermal Sciences 48 (2009) 1184–1195 1189

Table 2
Prescribed set of design points using the Box–Behnken technique.

BR ChanAR CorAR θ, deg Re

0.5 0.15 0.225 50 3250
0.5 0.05 0.225 50 3250
0.5 0.05 0.225 50 3250
0.1 0.05 0.225 50 3250
0.5 0.15 0.400 50 3250
0.5 0.25 0.400 50 3250
0.1 0.15 0.225 50 3250
0.5 0.15 0.225 75 3250
0.5 0.25 0.225 25 500
0.5 0.15 0.225 25 500
0.5 0.25 0.225 50 6000
0.1 0.25 0.225 50 6000
0.5 0.15 0.050 75 3250
0.5 0.25 0.050 50 3250
0.5 0.05 0.400 50 3250
0.5 0.15 0.400 50 3250
0.9 0.15 0.050 50 500
0.5 0.15 0.225 75 500
0.9 0.15 0.400 50 6000
0.5 0.05 0.225 75 6000
0.9 0.15 0.225 75 3250
0.5 0.05 0.225 25 3250
0.1 0.15 0.225 25 3250
0.9 0.15 0.225 25 3250
0.1 0.15 0.400 50 3250 Fig. 4. Nu vs Re for different corrugation shapes.
0.9 0.05 0.225 50 3250
0.5 0.15 0.400 75 3250
0.5 0.25 0.225 50 3250 4.1. Nusselt number estimation
0.5 0.15 0.225 25 500
0.9 0.15 0.225 50 500
The validity of the proposed correlation for Nusselt number es-
0.5 0.15 0.050 50 6000
0.5 0.25 0.225 75 6000 timation is checked using limited experimental data [18]. It is also
0.1 0.15 0.225 75 500 compared with the most popular correlation available in the lit-
0.5 0.05 0.050 50 500 erature, proposed by Martin [3]. In Fig. 4, the proposed plotted
0.9 0.15 0.225 50 6000
correlation for the Nu refers to a PHE channel with the same ge-
0.5 0.15 0.050 50 6000
0.1 0.15 0.050 50 3250 ometrical configuration as the one used in the experimental setup
0.5 0.15 0.400 25 3250 by Kanaris et al. [18] and found to be in excellent agreement. Com-
0.5 0.15 0.050 25 3250 parison with the correlation by Martin, used for a given angle
0.9 0.25 0.225 50 3250 of attack (θ = 67◦ ), shows that the two correlations are in good
0.1 0.15 0.225 50 3250
agreement (±20%). It must be noted that Martin’s correlation [3]
takes into account only the angle of attack. An interesting observa-
Table 3 tion is that the difference on the corrugation shape does not seem
Values of the coefficients of the quadratic models for: (a) ηNu , (b) η f . to affect the Nu values, i.e., the proposed correlation refers to tri-
(a) angular, available experimental data refer to rounded trapezoidal
ηNu 1 BR ChanAR CorAR sin 2θ Re and Martin’s correlation refers to sinusoidal corrugations.
1 −0.00220
BR −0.05756 0.31230 4.2. Friction factor estimation
ChanAR −0.28845 −0.72577 6.10191
CorAR −1.41326 −1.40661 −1.06706 4.58153
sin 2θ 0.94572 0.10292 −0.59416 −0.27910 −0.57396
In a similar way, the validity of the proposed correlation for
Re 1.56E–05 −1.30E–05 5.51E–05 −7.67E–05 −1.64E–05 4.90E–09 friction factor is checked using available experimental data [18]
(b)
and is also compared with the correlation available in the liter-
ature, proposed by Martin [3]. In Fig. 5, four different groups of
ηf 1 BR ChanAR CorAR sin 2θ Re
data are presented. More specifically:
1 5.35675
BR 1.31185 3.34690
ChanAR 3.15423 −0.40463 −8.15968 • data from the proposed correlation, referring to triangular cor-
CorAR −8.86597 2.47967 −10.92780 2.14140 rugations,
sin 2θ −12.19970 −2.24776 −0.83843 12.66120 7.51322 • data from CFD simulations concerning trapezoidal corrugations
Re 7.92E–4 2.00E–4 −5.74E–5 1.22E–3 −1.05E–4 −8.28E–08
[18],
• experimental data corresponding to rounded trapezoidal cor-
4. Results and discussion rugations [18] and finally,
• data calculated by the correlation of Martin [3] which concerns
sinusoidal corrugations.
The coefficients of the quadratic models for the thermal-related,
ηNu , and the friction-related, η f , functions were calculated using It is obvious that the sharpness of the corrugation influences the
the RSM and are presented in Table 3. For both functions, the ad- friction factor values, i.e., the sharpest the corrugation, the higher
justed coefficient of determination, R 2adj , is found to be greater than the friction factor. However, all data groups are well fitted with a
0.9. The ηNu and η f functions are used for the calculation of Nus- power law function, as suggested by Martin [3], with the value of
selt number, Nu, and friction factor, f . −0.135 being the exponent of Re.
1190 A.G. Kanaris et al. / International Journal of Thermal Sciences 48 (2009) 1184–1195

Table 5
Optimal values of design geometrical parameters for β = 0.03 and 0.06.

β = 0.03 β = 0.06
Re BR ChanAR CorAR θ BR ChanAR CorAR θ
500 0.900 0.226 0.400 41.2 0.771 0.250 0.393 38.6
1000 0.900 0.222 0.400 41.2 0.755 0.242 0.386 39.1
2000 0.900 0.216 0.400 41.7 0.666 0.203 0.340 42.2
3000 0.900 0.212 0.400 41.7 0.641 0.189 0.326 42.7
4000 0.900 0.206 0.400 42.2 0.617 0.171 0.308 43.2
5000 0.900 0.200 0.400 42.7 0.593 0.163 0.298 43.7
6000 0.900 0.196 0.400 42.7 0.585 0.153 0.287 44.2

Fig. 5. Friction factor vs. Re for different corrugation shapes.

Table 4
Friction factor correction and weighting factor correction for each shape case.

K Friction correction Weighting factor

( f = K Re−0.135 ) factor, f c 1/3
correction, βc = f c
Triangular 11.26
Sinusoidal 7.69 0.68 0.88
Rounded trapezoidal 9.38 0.83 0.94
Sharp trapezoidal 10.19 0.90 0.97
(a)

Consequently, a correction factor, f c , calculated as a ratio of

friction factor for one of the aforementioned shapes (i.e., sinu-
soidal, rounded trapezoidal or sharp trapezoidal) over the friction
factor for the case of triangular shaped corrugations, can be used
to account for the corrugation shape. The proposed values for the
correction factor are presented in Table 4.

4.3. Optimization results

In general the values of the four geometrical design variables

and the Re number which correspond to the optimal design are
calculated by the aforementioned method. Keeping Re constant, the
objective function is repeatedly optimized with respect to the four
geometrical parameters used as design variables for various Re val-
ues. This allows the designer to use the correlations derived in this
study to check the effect of Re on the optimal parameter values
and to design the best PHE for each case. The aforementioned cor-
relations can also be used for the calculation of the optimal design
specification when the value of any of the geometrical parameters
is pre-set (e.g., BR ≃ 1 or θ = 45◦ ).
As previously mentioned, the optimization of Eq. (13) has been
performed for various values of Re and for two values of the (b)
weighting factor, β, 0.03 and 0.06, in order to find the geometrical
Fig. 6. Optimal values of design geometrical parameters using CFD, for: (a) β = 0.03
configuration that provides the most efficient design. The values of
and (b) β = 0.06.
the geometrical parameters that correspond to this optimal design
are presented in Fig. 6 and are also tabulated in Table 5.
It is obvious that the pumping cost, incorporated in the weight- ChanAR ≃ 0.20). However, for β = 0.06 (i.e., pumping cost is more
ing factor, plays a significant role. As shown in Fig. 6, when significant), the values of the design variables for the optimal de-
β = 0.03, it is apparent that the optimal values of the geometri- sign of a PHE, not only differ from the previous values, but they
cal parameters are practically constant for the range of Re studied are also a function of Re. Since, in this case, the friction losses
(i.e., the plates must be placed abutting, θ ≃ 400 , CorAR = 0.40 and have a greater effect on the objective function values, an optimal
 A.G. Kanaris et al. / International Journal of Thermal Sciences 48 (2009) 1184–1195 1191

(a)

(b)

Fig. 7. Boundary layer breakup and reattachment of the flow on a cross section of a channel (ChanAR = 0.15, CorAR = 0.225, θ = 75, Re = 3250) for: (a) BR = 0.9 and
(b) BR = 0.5.

design leads to greater distances between the plates (a behavior In the case of CorAR, the increase of the parameter value leads
also mentioned in the similar work by Kim and Kim [10]) and less to an increase on the reattachment length. This is apparent in
sharp corrugations (i.e., smaller values of BR and CorAR). Fig. 8a and c, where the type of flow is presented for two CorAR
The effect of the distance between the plates, which is ex- values, 0.40 and 0.15, for BR = 0.9, ChanAR = 0.25, θ = 75◦ and
pressed by BR, is presented in Fig. 7 for the case of two closely Re = 500. For low values of the weighting factor (e.g., β = 0.03)
placed plates (Fig. 7a, BR = 0.9) and for the case that the plates the effect of this parameter on the overall energy economy is in-
are located in some distance apart from each other (Fig. 7b, significant (Fig. 6a), whereas for β = 0.06 the optimal CorAR value
BR = 0.5), all other parameters being the same (ChanAR = 0.15, decreases with Re (Fig. 6b) and is always smaller than that for
CorAR = 0.225 and θ = 75◦ ). It is obvious that, when the plates β = 0.03 (Table 5). In the work of Kim and Kim [10], it is stated
are placed abutting, the type of flow inside the channel greatly af- that CorAR reaches its highest value as the weighting factor de-
fects the heat transfer augmentation and the friction losses, due to creases, in agreement with the results of the present study. In
Fig. 9, it is shown that friction losses increase linearly with CorAR.
the ‘restart’ of the boundary layer.
The upper limit of the design variable used for the calculations
As for the angle of attack, θ , while for low Re the optimal de-
is imposed by PHE construction constrains, as mentioned earlier.
sign leads to ‘soft’ plates (i.e., with small angles, θ = 30◦ –40◦ ),
Nevertheless, in cases of low pumping cost (i.e., low values of β)
for higher Re the optimal values are shifted to geometrical con-
the designer could explore the use of sharper corrugations.
figurations with greater angles (‘harder’ plates) on the corrugation
Optimal values of ChanAR, the variable that expresses how nar-
pattern. This is in agreement with Martin [3], who observes that,
row or wide the PHE channel is, decreases with Re (for the Re
for a given pressure drop, ‘harder’ plates result to higher temper-
range studied). This parameter represents the effect of the side
ature changes. In Fig. 8a and b, the type of flow inside a chan- walls on the performance of the PHE. It is reported that two
nel is presented for two different corrugation angles, i.e., θ = 75◦ kinds of flow occur in the PHE: the ‘crossing flow’ and the ‘zig-zag
and θ = 25◦ , respectively, while all the other geometrical param- flow’ [5]. During crossing flow (Fig. 10a), small sub-streams follow
eters retain the same values (i.e., BR = 0.9, ChanAR = 0.25 and the furrows of the corrugated plates, while, during ‘zig-zag flow’
CorAR = 0.40). It is shown that the flow reattachment length is (Fig. 10b), the fluid ‘walks’ over the corrugations, by rapidly chang-
greater for higher values of the angle of attack. It is noted that in ing direction. It is known [5] that the wall shear stress attains its
Fig. 8a and b, the different angles of attack and the fact that the peak values on the crests of the corrugations, resulting in increased
cross sections lie on a plane parallel to the flow create the illusion heat flux, and this study confirms this behavior (Fig. 11). There-
of different CorAR. It is obvious that on ‘harder’ plates the ‘restarts’ fore, when the flow is forced to follow the ‘zig-zag flow’, friction
of boundary layer are more frequent, contrary to the case of low losses increase. Existence of ‘crossing-flow’ in this type of equip-
angles of attack. ment is one of the results of the reflection of fluid flowing in the
1192 A.G. Kanaris et al. / International Journal of Thermal Sciences 48 (2009) 1184–1195

(a)

(b)

(c)

Fig. 8. Boundary layer breakup and reattachment of the flow on a cross section of a channel (ChanAR = 0.25, BR = 0.9, Re = 500) for: (a) CorAR = 0.40, θ = 75◦ , (b) CorAR =
0.40, θ = 25◦ and (c) CorAR = 0.15, θ = 75◦ .

furrows upon channel side walls. For a channel with low ChanAR shown that the geometrical parameter values that correspond to
(i.e., ‘wide’ channel) reflections on side walls become less frequent, the optimal design remain practically unaffected when β < 0.04.
because for long furrows, part of the flow inside them tends to If β exceeds 0.04, the friction-related term seems to become the
change direction and to follow the ‘zig-zag flow’ pattern. As the controlling step of the optimization, and leads to geometrical con-
weighting factor increases, the optimal value of ChanAR reaches figurations with lower pressure drop.
higher values, resulting in narrower channels, where pressure drop It must be pointed out that, the corrugation-shape correction
decreases.
factor ( f c ) is raised to the 1/3 power when it is included in
the objective function and thus its significance on the objective
4.4. Effect of weighting factor
function values is reduced (Eq. (9)). So finally a new weighting
The optimal values of the objective function have been also cal- factor can be defined, which is the product of β and βc , where
1/ 3
culated for various values of the weighting factor β and for two βc = f c (Table 4) and the objective function for optimization be-
typical Re (1000 and 5000) and are presented in Fig. 12. It is comes:
 A.G. Kanaris et al. / International Journal of Thermal Sciences 48 (2009) 1184–1195 1193

Fig. 9. Effect of CorAR on the values of the thermal and friction related functions, for BR = 0.5, ChanAR = 0.15, θ = 50◦ and Re = 6000.

(a) (b)

Fig. 10. Flow streamlines presented for different channel aspect ratio: (a) ChanAR = 0.25 (narrow channel), (b) ChanAR = 0.05 (wide channel).
1194 A.G. Kanaris et al. / International Journal of Thermal Sciences 48 (2009) 1184–1195

Fig. 11. Shear stress and heat flux distribution on the corrugated plate.

ology (RSM) is used as a global optimization method and results

are presented in terms of Nusselt number and friction factor. Box–
Behnken design was selected for the four geometrical design fac-
tors and the Reynolds number to construct the response surface.
Results from the CFD simulations show very good agreement with
correlations for this type of industrial PHE from the literature as
well as with available experimental data. An objective function is
appropriately defined as a linear trade-off between heat transfer
and pressure drop, combined with a weighting factor accounting
for the cost of energy. The correlations proposed for the calcu-
lation of heat transfer rate and friction losses can be considered
a valuable tool for the optimal design of a PHE. Additionally, the
designer can independently choose which of the five available de-
grees of freedom to use as design variables to optimize a PHE.
More specifically, it is shown that optimal performance is
achieved as the channel plates are coming closer to each other
and for less obtuse (i.e., sharper) corrugations, while the PHE per-
formance can be improved for lower values of channel aspect ratio
(i.e., wider channels) and for higher values of the angle of attack,
as Re increases. In all cases, the designer can be led to the optimal
geometrical configuration of a PHE with the use of values of de-
sign parameters that create intense flow inside the furrows, a case
where secondary flow due to the corrugations is dominant. This
Fig. 12. Optimal values of design geometrical parameters vs the weighting factor for
Re = 1000 and Re = 5000. secondary flow increases flow separation and reattachment, which
is a known mechanism that improves heat transfer augmentation.
F = ηNu + (β · βc )η f (14) As the weighting factor increases, meaning that the friction losses
play a very important role, the optimal design of a PHE dictates
It is interesting that even for the smoothest shape (i.e., si- greater distances between the plates and less sharp corrugations,
nusoidal) the effect of the correction factor would be about for high values of Re. More work can be done to check the effect
10%. This correction factor is also assumed to be independent of flow configuration in the geometrical parameter values that op-
of BR. timize the performance of a PHE complying with the principles of
ecological and economic sustainability.
5. Conclusions
Acknowledgements
A commercial CFD code validated previously by experimental
work in our Lab is used for the calculation of heat transfer rate Financial support by the General Secretariat for Research and
and friction losses in a channel of a PHE with undulated surfaces. Technology and the European Union (PENED 2003) is greatly ac-
This study extends the use of CFD for seeking the optimal design knowledged. The authors wish to thank Prof. V. Bontozoglou, Prof.
for this type of equipment. Results presented in this study consider S.G. Yiantsios and Prof. N. Sahinidis for their helpful comments and
the use of water as working fluid. The response surface method- suggestions.
 A.G. Kanaris et al. / International Journal of Thermal Sciences 48 (2009) 1184–1195 1195

References [14] W.W. Focke, P.G. Knibbe, Flow visualization in parallel-plate ducts with
corrugated walls, J. Fluid Mech. 165 (1986) 73–77.
[15] J.E. Hesselgraves, Compact Heat Exchangers: Selection, Design and Operation,
[1] R.K. Shah, A.S. Wanniarachchi, Plate heat exchanger design theory, in: J.M.
1st ed., Pergamon, New York, 2001.
Buchlin (Ed.), Industrial Heat Exchangers, von Karman Institute Lecture Series,
[16] P. Vlasogiannis, G. Karagiannis, P. Argyropoulos, V. Bontozoglou, Air–water
1991.
two-phase flow and heat transfer in a plate heat exchanger, Int. J. Multiphase
[2] W.M. Kays, A.L. London, Compact Heat Exchangers, 3rd ed., Krieger Publ. Co.,
Flow 28 (5) (2002) 757–772.
Florida, USA, 1998.
[17] A.G. Kanaris, A.A. Mouza, S.V. Paras, Flow and heat transfer in narrow channels
[3] H. Martin, A theoretical approach to predict the performance of chevron-type
with corrugated walls: a CFD code application, Chemical Engineering Research
plate heat exchangers, Chem. Eng. Proc. 35 (1996) 301–310.
and Design 83 (A5) (2005) 460–468.
[4] R.K. Shah, D.P. Sekulić, Fundamentals of Heat Exchanger Design, John Wiley and
[18] A.G. Kanaris, A.A. Mouza, S.V. Paras, Flow and heat transfer prediction in a
Sons, Inc., Hoboken, NJ, 2003.
corrugated plate heat exchanger using a CFD code, Chemical Engineering and
[5] W.W. Focke, J. Zachariades, I. Olivier, The effect of the corrugation inclination
Technology 29 (8) (2006) 923–930.
angle on the thermohydraulic performance of plate heat exchangers, Int. J. Heat
[19] A. Bejan, A.D. Kraus, Heat Transfer Handbook, John Wiley and Sons, Inc.,
Mass Transfer 28 (1985) 1469–1497.
Canada, 2003.
[6] H. Heavner, R.L Kumar, A.S. Wanniarachchci, Performance of an industrial plate
[20] CFX® Release 10.0 Manual, ANSYS Inc. International, Canonsburg, 2005.
heat exchanger: Effect of chevron angle, in: AIChE Symposium Ser. Heat Trans-
[21] Menter F., Esch T., Elements of industrial heat transfer predictions, in: 16th
fer, Am. Inst. Chem. Eng. Atlanta, GA, 1993.
Brazilian Congress of Mechanical Engineering (COBEM), Uberlandia, Brazil,
[7] M.A. Mehrabian, R. Poulter, Hydrodynamics and thermal characteristics of cor- 2001.
rugated channels: computational approach, Applied Mathematical Modeling 24 [22] R.H. Myers, D.C. Montgomery, Response Surface Methodology: Process and
(2000) 343–364. Product Optimization Using Designed Experiments, 2nd ed., John Wiley & Sons
[8] M.Z. Hossain, A.K.M. Sadrul Islam, Numerical investigation of fluid flow and Inc., Canada, 2002.
heat transfer characteristics in sine, triangular and arc-shaped channels, Ther- [23] G.E.P. Box, D.W. Behnken, Some new three level designs for the study of
mal Science 11 (1) (2007) 17–26. quantitative variables, Technometrics 2 (4) (1960) 455–475.
[9] Y. Asako, H. Nakamura, M. Faghri, Heat transfer and pressure drop characteris- [24] Matlab 6.5 Manual, The Mathworks, Inc., 2002.
tics in a corrugated duct with rounded corners, Int. J. Heat Mass Transfer 31 (6) [25] H.-M. Kim, K.-Y. Kim, Shape optimization of three-dimensional channel rough-
(1988) 1237–1245. ened by angled ribs with RANS analysis of turbulent heat transfer, Int. J. Heat
[10] H.-M. Kim, K.-Y. Kim, Design optimization of rib-roughened channel to enhance Mass Transfer 49 (2006) 4013–4022.
turbulent heat transfer, Int. J. Heat Mass Transfer 47 (2004) 5159–5168. [26] M.A. Ansari, K.-Y. Kim, Shape optimization of a micromixer with staggered
[11] K.-Y. Kim, Y.-M. Lee, Design optimization of internal cooling passage with V- herringbone groove, Chem. Eng. Sci. 62 (2007) 6687–6695.
shaped ribs, Numerical Heat Transfer A 51 (2007) 1103–1118. [27] M.E. Taslim, C.M. Wadsworth, An experimental investigation of the rib surface-
[12] J. Zhang, J. Kundu, R.M. Manglik, Effect of fin waviness and spacing on the averaged heat transfer coefficient in a rib-roughened square passage, Journal
lateral vortex structure and laminar heat transfer in wavy-plate-fin cores, Int. of Turbomachinery – Transactions of ASME 119 (2) (1997) 381–389.
J. Heat Mass Transfer 47 (2004) 1719–1730. [28] D.L. Gee, R.L. Webb, Forced convection heat transfer in helically rib-roughened
[13] L. Wang, B. Sundén, Optimal design of plate heat exchangers with and with- tubes, Int. J. Heat Mass Transfer 23 (1980) 1127–1136.
out pressure drop specifications, Applied Thermal Engineering 23 (2003) 295– [29] A. Bejan, Entropy Generation Through Heat and Fluid Flow, John Wiley and
311. Sons, Inc., New York, 1982.