Modelling and Numerical Simulation of
Fluid Flow and Heat Transfer in
Thermoplates
zur Erlangung des akademischen Grades eines
DOKTORS DER INGENIEURWISSENSCHAFTEN (Dr.-Ing.)
der Fakult¨at f¨ur Maschinenbau
der Univesit¨at Paderborn
genehmigte
DISSERTATION
von
Dipl.-Ing. Boban Maleti´c
aus Bijeljina, Bosnien und Herzegowina
Tag des Kolloquiums: 20. Februar 2009
Referent: Prof. Dr.-Ing. Jovan Mitrovi´c
Korreferent: Prof. Dr.-Ing. Eugeny Kenig
Contents
List of Publications iii
List of Symbols ix
1 Introduction 1
2 State of the Art 3
2.1 Thermoplates ................................. 3
2.2 Parallelplatechannel............................. 5
2.3 Two-dimensional wavy channel . . . . . . . . . . . . . . . . . . . . . . . 7
2.4 Three-dimensional wavy channel . . . . . . . . . . . . . . . . . . . . . . . 15
3 Heat Transfer in a Parallel Plate Channel 19
3.1 Physical model and governing equations . . . . . . . . . . . . . . . . . . 19
3.2 Normalised temperature function and analytical solution . . . . . . . . . 20
3.2.1 Characteristic temperature distributions . . . . . . . . . . . . . . 23
3.3 TheNusseltNumber ............................. 27
3.3.1 Positions of the vertical asymptote and the adiabatic point . . . . 35
4 Fluid Flow and Heat Transfer in Thermoplates 39
4.1 Problemformulation ............................. 39
4.2 Thephysicalmodel.............................. 43
4.3 Processquantities............................... 47
4.3.1 Symmetric heated channels . . . . . . . . . . . . . . . . . . . . . 47
4.3.2 Asymmetric heated channels . . . . . . . . . . . . . . . . . . . . . 50
5 Results 53
5.1 Symmetric heated channels, steady state flows . . . . . . . . . . . . . . . 53
5.1.1 Thevelocityfield........................... 53
5.1.2 The temperature field and the heat flux at the thermoplate surface 57
5.1.3 Influence of the geometrical parameters . . . . . . . . . . . . . . . 62
5.1.3.1 Welding spot spacing in the flow direction . . . . . . . . 62
5.1.3.2 Welding spot spacing in the cross direction . . . . . . . . 65
5.1.3.3 Maximal plate distance . . . . . . . . . . . . . . . . . . 68
5.1.3.4 Nusselt number correlation . . . . . . . . . . . . . . . . 71
i
Contents
5.1.3.5 The thermo-fluid characteristic . . . . . . . . . . . . . . 72
5.1.4 Welding spots pattern and thermoplate surface shape . . . . . . . 75
5.2 Symmetric heated channels, transient flows . . . . . . . . . . . . . . . . . 78
5.3 Asymmetric heated channels . . . . . . . . . . . . . . . . . . . . . . . . . 84
5.4 Validation through Experiments . . . . . . . . . . . . . . . . . . . . . . . 88
6 Conclusions 95
Bibliography 97
A Numerical Method 103
A.1 The finite volume method . . . . . . . . . . . . . . . . . . . . . . . . . . 103
A.1.1 Steady state calculations . . . . . . . . . . . . . . . . . . . . . . . 104
A.1.1.1 Differencing schemes . . . . . . . . . . . . . . . . . . . . 105
A.1.1.2 The SIMPLE algorithm . . . . . . . . . . . . . . . . . . 107
A.1.1.3 Convergence criteria . . . . . . . . . . . . . . . . . . . . 109
A.1.2 The finite volume method for time dependent flows . . . . . . . . 112
A.2 Discretised integrals for the process quantities . . . . . . . . . . . . . . . 113
B Definition of the Model for StarCD 115
C POSDAT.f 125
D Mesh Test and Other Numerical Details 135
E Constants 137
F Numerical Results 139
Acknowledgements 147
ii
List of Publications
1. Mitrovi´c, J. and Maleti´c, B., 2005, Effect of thermal assymetry on heat transfer in
a laminar annular flow, Chemical Engineering and Technology, 28, 1144-1150.
2. Mitrovi´c, J. and Maleti´c, B., 2006, Effect of thermal asymmetry on heat laminar
forced convection heat transfer in a porous annular channel, Chemical Engineering
and Technology, 29, 750-760.
3. Mitrovi´c, J., Maleti´c, B. and Baˇcli´c, B.S., 2006, Some peculiarities of the asymetric
Graetz problem, International Journal of Engineering Science, 44, 436-455.
4. Mitrovi´c, J. and Maleti´c, B., 2006, Numerical simulation of fluid flow and heat
transfer in thermoplates, Proceedings of the 13th International Heat Transfer Con-
ference IHTC-13, HEX-10, Sydney, Australia.
5. Mitrovi´c, J. and Maleti´c, B., 2007, Heat transfer with laminar forced convection in
a porous channel exposed to a thermal asymmetry, International Journal of Heat
and Mass Transfer, 50, 1106-1121.
6. Mitrovi´c, J. and Maleti´c, B., 2007, Numerical simulation of fluid flow, heat transfer
and pressure drop in thermoplates, Proceedings of the 5th International Conference
on Heat Transfer, Fluid Mechanics and Thermodynamics HEFAT2007, MJ3, Sun
City, South Africa.
7. Mitrovi´c, J. and Maleti´c, B., 2007, Numerical simulation of fluid flow and heat
transfer in thermoplates, International Journal of Heat Exchangers, in Press.
8. Maletic Maleti´c, B. and Mitrovi´c, J., 2008, Influence of the thermoplate geometry
on the heat transfer, Proceedings of the 5th European Thermal-Sciences Confer-
ence EUROTHERM 2008, HEX27, Eindhoven, The Netherlands.
iii
List of Publications
iv
List of Figures
2.1 Fluid flow arrangement in thermoplates (DEG-Engineering [64]) . . . . . 3
2.2 Different types of thermoplates, Behrend [3] . . . . . . . . . . . . . . . . 4
2.3 Thermoplate assembly (OMEGA heat transfer technology [66]) . . . . . . 4
2.4 The use of thermoplate heat exchangers as a head condenser, M¨uhlthaler
[39]....................................... 5
2.5 Plate-and-frame heat exchanger, (GEA FlatPlate Inc. [65]) . . . . . . . . 6
2.6 Fluid flow through a parallel plate channel . . . . . . . . . . . . . . . . . 6
2.7 Distribution of Nusselt number allong a parallel plate channel, Hatton
andTurton[19]................................ 7
2.8 Different types of wavy channels . . . . . . . . . . . . . . . . . . . . . . . 9
2.9 Axial velocity profiles at two different Reynolds numbers for a/Hmax =
0.175 and L/Hmax = 1.5, Mahmud et al. [28] . . . . . . . . . . . . . . . . 10
2.10 Mean streamlines at Re = 500 visualised by using PIV method, Kim [24] 11
2.11 Distribution of the local Nusselt number for Re = 500 and Pr = 6.93, [69] 12
2.12 Distribution of the average Nusselt number for Re = 500 and Pr = 6.93,
[69]....................................... 13
2.13 Distribution of the local Nusselt number for a/L = 0.2 and Pr = 6.93, [69] 13
2.14 Streamlines predicted numerically by Guzm´an and Amon [18] for different
Reynoldsnumbers............................... 14
2.15 Dimple plate geometry numerically analysed by Witry et al. [73] . . . . . 16
2.16 Internal flow velocity vectors in the plate radiator calculated by Witry et
al.[73] .................................... 17
2.17 Internal flow pressure distribution in the plate radiator calculated by
Witryetal.[73]................................ 17
2.18 Internal flow temperature distribution in the plate radiator calculated by
Witryetal.[73]................................ 17
2.19 Channel with two-dimensional corrugations, Sawyers et al. [56] . . . . . . 18
2.20 Channel with three-dimensional corrugations, Sawyers et al. [56] . . . . . 18
3.1 Illustration of the physical model . . . . . . . . . . . . . . . . . . . . . . 19
3.2 Temperature field at specified boundary conditions, ΘW1= 0 and ΘW2=
0.4, with hydrodynamically developed flow, Equation (3.3) . . . . . . . . 24
3.3 Temperature profiles for some fixed values of Z∗. At Z∗= 0.6816 the
upperplateisadiabatic............................ 24
v
List of Figures
3.4 The temperature slope at the plates for ΘW1= 0 and ΘW2= 0.4 . . . . . 24
3.5 Temperature field at specified boundary conditions, ΘW1= 0 and ΘW2=
1, with hydrodynamically developed flow, Equation (3.3) . . . . . . . . . 25
3.6 Temperature profiles for some fixed values of Z∗for ΘW1= 0 and ΘW2= 1 25
3.7 The temperature slope at the plates for ΘW1= 0 and ΘW2= 1 . . . . . . 25
3.8 Temperature field at specified boundary conditions, ΘW1= 0 and ΘW2=
2, with hydrodynamically developed flow, Equation (3.3) . . . . . . . . . 26
3.9 Temperature profiles for some fixed values of Z∗for ΘW1= 0 and ΘW2= 2 26
3.10 The temperature slope at the plates for ΘW1= 0 and ΘW2= 1 . . . . . . 26
3.11 Distribution of the local Nusselt numbers Nu(Z∗) along the flow direction
for ΘW1= 0 and ΘW2= 0.4 ......................... 29
3.12 Effect of boundary conditions on Nusselt number for hydrodynamically
developed flow for ΘW1= 0 at different ΘW2................ 30
3.13 Effect of boundary conditions on Nusselt number for hydrodynamically
developed flow for ΘW1= 0 at different ΘW2................ 31
3.14 Transitions and asymptotical behaviour of the Nusselt numbers for chosen
ΘW2...................................... 33
3.15 Effect of velocity distribution on the position of the vertical asymptote
for ΘW2= 2/5. (a) w= ¯w=const,(b)w=3
2¯w(1 −X2), (c) Pr = 32/3,
w= ¯wf(X, Z), Equation (3.1), (d) Pr = 0.7, w= ¯wf(X, Z), Equation
(3.1) ...................................... 34
3.16 Transitions and asymptotical behaviour of the Nusselt numbers for the
region ΘW2<1 ................................ 37
4.1 Geometry and dimensions of the three-dimensional simulated domain; a)
welding spots arrangend in a staggered manner, Equation 4.1, b) parallel
welding spots, Equation (4.2) with C=sL/2, c) parallel welding spots,
Equation (4.2) with C=sL/4, d) parallel welding spots, Equation (4.2)
with C= 0, and e) thermoplate surface defined using the exponential
function(4.4) ................................. 40
5.1 Velocity field in the symmetry plane y= 0 at different Reynolds numbers 54
5.2 Velocity profiles in the welding spots cross-sections at Re = 4000 . . . . . 55
5.3 Pressure field in the symmetry plane y= 0 at Re = 4000 and its linkage
withthevelocityfield............................. 56
5.4 Temperature field in the symmetry plane y= 0 at different Reynolds
numbers.................................... 58
5.5 Wall heat flux at different Reynolds numbers . . . . . . . . . . . . . . . . 59
5.6 Local, spanwise averaged, heat flux, ˙qW= ˙qW(x), and the mean heat
flux, ˙qmW = ˙qmW (x), according to Equations (4.31) and (4.32) at different
Reynoldsnumbers............................... 60
5.7 Local, spanwise averaged, heat transfer coefficient, hW=hW(x), and the
mean heat transfer coefficient, hmW =hmW (x), according to the Equation
(4.33) at different Reynolds numbers . . . . . . . . . . . . . . . . . . . . 61
vi
List of Figures
5.8 Heat transfer coefficients, hW=hW(x) and hmW =hmW (x), according to
the Equation (4.33), with fully developed region in a 500mm long channel
at Re =50................................... 61
5.9 Average heat transfer coefficient, hmW of 500 mm long channels at dif-
ferent Re as a function of longitudinal welding spots pitch, sL;sL→ ∞:
parallelplates................................. 63
5.10 Average Nusselt number, NumW of 500 mm long channels at different Re
as a function of longitudinal welding spots pitch, sL;sL/L → ∞: parallel
plates ..................................... 63
5.11 Average Nusselt number, NumW of 500 mm long channels as a function
of Re at different sL/L;sL/L → ∞: parallel plates . . . . . . . . . . . . 64
5.12 Pumping power, P, Equation (4.45), needed for forced convection through
a 500 mm long channel as a function of the longitudinal welding spots
pitch, sL, at different Reynolds numbers, Re;sL→ ∞: parallel plates . . 64
5.13 Average heat transfer coefficient, hmW of 500 mm long channels at differ-
ent Re as a function of transversal welding spots pitch, sT........ 66
5.14 Average Nusselt number, NumW of 500 mm long channels at different Re
as a function of transversal welding spots pitch, sT............ 66
5.15 Average Nusselt number, NumW of 500 mm long channels as a function
of Re at different sT/L ............................ 67
5.16 Pumping power, P, Equation (4.45), needed for forced convection through
a 500 mm long channel as a function of the transversal welding spots pitch,
sT, at different Reynolds numbers, Re .................... 67
5.17 Average heat transfer coefficient, hmW of 500 mm long channels at differ-
ent Re as a function of maximal plate distance, δ............. 69
5.18 Average Nusselt number, NumW of 500 mm long channels at different Re
as a function of maximal plate distance, δ................. 69
5.19 Average Nusselt number, NumW of 500 mm long channels as a function
of Re at different δ/L ............................. 70
5.20 Pumping power, P, Equation (4.45), needed for forced convection through
a 500 mm long channel as a function of maximal plate distance, δ, at
different Reynolds numbers, Re ....................... 70
5.21 Thermo-fluid characteristic of a 500 mm long thermoplate channel as
function of geometry at selected Reynolds numbers; sL/L → ∞: par-
allelplates .................................. 73
5.22 Thermo-fluid characteristic of the thermoplate compared with the char-
acteristic of the flat channel at selected Reynolds numbers . . . . . . . . 74
5.23 Distribution of the local heat transfer coefficient, Equation (4.33), for
different models, Figure 4.1a, c and e, at Re =50 ............. 75
5.24 Distribution of the local heat transfer coefficient, Equation (4.33), for
different models, Figure 4.1a, c and e, at Re =2000............ 76
5.25 Distribution of the local heat transfer coefficient, Equation (4.33), for
different models, Figure 4.1b, c and d, at Re =2000............ 76
vii
List of Figures
5.26 Average Nusselt number at different Reynolds numbers for different mod-
els,Figure4.1................................. 77
5.27 Pumping power number at different Reynolds numbers for different mod-
els,Figure4.1................................. 77
5.28 Thermo-fluid characteristic at different Reynolds numbers for different
models,Figure4.1 .............................. 78
5.29 Development of the velocity field in the symmetry plane y= 0 with time 80
5.30 Development of the temperature field in the symmetry plane y= 0 with
time ...................................... 81
5.31 Development of the heat flux at the thermoplate wall with time . . . . . 82
5.32 Heat transfer coefficient averaged in space for a 120 mm long model,
steadyandunsteadyflow........................... 83
5.33 Local heat transfer coefficient for a 120 mm long model, at different time
steps...................................... 84
5.34 Different inlet fluid velocities . . . . . . . . . . . . . . . . . . . . . . . . . 85
5.35 Average heat transfer coefficient (averaged in space) of a 300 mm long
model for different inlet fluid velocities . . . . . . . . . . . . . . . . . . . 85
5.36 Nusselt number distributions at different ΘW2; ΘW1= 0, Re = 100. The
thin lines represent the corresponding distributions for the parallel plate
channel,Figure3.13 ............................. 86
5.37 Nusselt number distributions at different ΘW2; ΘW1= 0, Re = 100. The
thin lines represent the corresponding distributions for the parallel plate
channel .................................... 87
5.38 Geometry of the thermoplate used in experiments by Mitrovi´c and Peter-
son,[38] .................................... 88
5.39 Flow field and temperature field distribution across the thermoplate; Mar-
lotherm oil at Re = 1860, ˙mIN = 510kg/h, TIN = 328.2K, TW= 333.2K . 91
5.40 Flow field and temperature field distribution across the thermoplate; Mar-
lotherm oil at Re = 5416, ˙mIN = 1500kg/h, TIN = 328.1K, TW= 337.2K 92
5.41 Nusselt number comparison for a thermoplate strip . . . . . . . . . . . . 93
5.42 Nusselt number comparison for whole thermoplate . . . . . . . . . . . . . 93
5.43 Pressure drop comparison for whole thermoplate . . . . . . . . . . . . . . 94
A.1 A three-dimensional cell with the node Pand neighbouring nodes W,E,
S,N,Band T,[68]..............................104
A.2 Contours of constant φat different P´eclet numbers [68] . . . . . . . . . . 106
A.3 The upwind differencing scheme for one-dimensional flows [68] . . . . . . 107
A.4 The two-dimensional, backward staggered grid [68] . . . . . . . . . . . . 110
A.5 The SIMPLE algorithm [68] . . . . . . . . . . . . . . . . . . . . . . . . . 111
A.6 The transient SIMPLE algorithm, [68] . . . . . . . . . . . . . . . . . . . 112
A.7 Definitions of the elementary volume and surface . . . . . . . . . . . . . 114
D.1 Illustration of the grid used in the numerical simulations . . . . . . . . . 136
viii
List of Figures
x
List of Tables
3.1 Some eigenvalues λnof the asymmetrical Graetz problem . . . . . . . . . 23
4.1 Parameter combinations for the numerical experiments. . . . . . . . . . . 42
4.2 Extended parameter combinations for the numerical experiments. . . . . 42
4.3 Different model geometries. . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.4 Inlet velocity for time dependent fluid flows. . . . . . . . . . . . . . . . . 43
4.5 Thermally asymmetric boundary conditions. . . . . . . . . . . . . . . . . 43
5.1 Constants used in Equation (5.2). . . . . . . . . . . . . . . . . . . . . . . 71
5.2 Constants used in Equation (5.2). . . . . . . . . . . . . . . . . . . . . . . 71
5.3 Volume flow rate comparison for different models, Figure 4.1. . . . . . . . 75
5.4 Boundary conditions for the comparison with the experiments of Mitrovi´c
andPeterson. ................................. 89
E.1 Constants used in Equation (5.3) with I=J=K=L= 2 and M=
m=0. ....................................137
F.1 Simulation results for model parameters given in Tables 4.1 and 4.2. . . 139
xi
List of Tables
xii
List of Symbols
Latin symbols
Symbol Unit Meaning
Am2wall surface
An−constants
am amplitude, plate width
ai,j, aP, anb different coefficients in discretised equations
bm plate spacing
C−constant
Cr−Courant number
cpJ·kg−1·K−1specific heat capacity at constant pressure
cvJ·kg−1·K−1specific heat capacity at constant volume
D, Dn−constants
dm average plate distance
dhm hydraulic diameter
En−constants
1F1(x, y, z)−hypergeometric function
fm·s−2gravitational force vector
Hm height
H(x) - Heaviside step function
Hx(y)−Hermite polynomial
hW·m−2·K−1heat transfer coefficient
hmW·m−2·K−1average heat transfer coefficient
i, j, k, l, m - summation indices
iJ·kg−1specific internal energy
kW·m−1·K−1thermal conductivity, iteration number
Lm length, half of the wave-length
lm cell length
Mkg·kmol−1molar mass
M−number of boundary cell faces
mkg mass
˙mkg·s−1, kg·h−1, mass flow rate
N−total number of welding spots, number of numerical cells
Nu −Nusselt number
xiii
List of Symbols
Symbol Unit Meaning
nm coordinate normal to the surface
n−normal vector
PW pumping power
Pe −P´eclet number
Pr −Prandtl number
pPa, bar pressure
˙
QW heat flow rate
˙qW·m−2heat flux vector
Rm welding spot radius
Re −Reynolds number
Sm2cross-sectional flow area
sLm distance between welding points in x-direction
sTm distance between welding points in x-direction
TK thermodynamic temperature
ts time
Um circumference
¯um·s−1average velocity component in x-direction
u, v, w m·s−1components of the velocity vector w
˜u∗, ˜v∗m·s−1initial guess velocity components
Vm3Volume
˙
Vm3·s−1volume flow rate
wm·s−1velocity vector
¯wm·s−1average velocity component in z-direction
x, y, z m coordinates
X, Y, Z, Z∗−dimensionless coordinates
x0, z0m coordinates of the welding spots centres
ZA∗−dimensionless coordinate of the vertical asymptote
Z0∗−dimensionless coordinate of the adiabatic point
Greek symbols
α−weighting parameter, 0 ≤α≤1
δm maximal plate distance
εm width of the numerical cell in x-direction
Θ−dimensionless temperature
¯
Θ−average dimensionless temperature
ϑ◦C temperature
¯
ϑ◦C average temperature
κm2·s−1thermal diffusivity
xiv
Symbol Unit Meaning
λ−tolerance
λn−eigenvalues
µPa·s dynamic viskosity
νm2/s kinematic viskosity
ξm coordinate
π= 3.141... −Ludolph’s constant
ρkg·m−3density
τPa stress tensor
Φ s−2dissipation funktion
φ−general scalar variable
ϕn(X)−eigenfunctions
ϕ◦phase shift angle
ωn−nth rooth of the equation tan(ω) = ω
Subscripts
Avertical asymptote
i, I summation index
IN value at the inlet
j, J summation index
Llongitudinal, spacing in x-direction
maverage value
nsummation index
nb neighbouring cells
OUT value at the outlet
Sa narrow strip of the thermoplate surface y=f(x, z)
SW1 a narrow strip of the lower thermoplate surface
SW2 a narrow strip of the upper thermoplate surface
Ttransversal, spacing in z-direction
Wvalue at the wall
W1 value at the lower wall, X=−1 or y=−f(x, z)
W2 value at the upper upper wall, X= +1 or y= +f(x, z)
0 adiabatic point
∞value in the fully developed region
xv
List of Symbols
Symbol Unit Meaning
Superscripts
0 value at starting time level t
∗initial guessed value
xvi
1 Introduction
Thermoplates are efficient heat transfer devices, which are used in several branches of
engineering practice, e.g. as condensers or evaporators in thermal process technology
and cooling technique. In comparison to shell-and-tube heat exchanges, their installation
is relatively simple, and the periphery is drastically reduced.
A thermoplate consists of two metallic sheets, which are spot-welded over the whole
surface according to an appropriate pattern, whereas the edges - except for connecting
tubes - are continuously seam-welded. By applying a hydro-form technique, a chan-
nel having a complex shape is established between the sheets. One fluid is conducted
through this channel, the other one through the channel created by two neighbouring
thermoplates that are assembled in parallel at certain spacing thus making a thermoplate
heat exchanger. The complex geometry of the inside channel of such a plate provides
the flowing fluid a pronounced three-dimensional character.
In the simulations, the geometrical parameters, such as the welding pattern, the
surface shape, the streamwise welding spots pitch, sL, the transversal pitch, sT, and the
maximal distance between the metallic sheets, δ, have been varied, whereas the radius,
R, of the welding spots and the channel length, L, were kept constant. The Reynolds
number, Re, that is, the fluid inlet velocity, uIN , and the distance, δ, has also been
varied. The fluid inlet temperature, TIN , and the wall temperature, TW, were taken to
be constant.
A parallel plate channel, and from some points of view, its equivalents, duct flow and
annular gap, can be considered as the simplest thermoplate having even walls. Although
a fabrication of such a channel without welding spots seems scarcely possible in sizes
required in common practice, the parallel plate channel will serve in this work as a model
for comparison with thermoplates.
The fluid flow was considered to be laminar, incompressible, at steady-state and
three-dimensional. Water of constant physical properties is adopted for the numerical
experiments. The velocity and the temperature fields are governed by the equations of
continuity, momentum and energy.
The objective of the numerical experiments is to numerically obtain the sets of the
geometrical parameters that, in interaction with process parameters, should pave the way
for optimal heat transfer of the inside fluid which is assumed to pass the thermoplate as
a single phase.
1
1 Introduction
For the case of fluid flow and heat transfer in a parallel plate channel it was possible
to find an analytical solution to the energy equation. In this connection the Mathemat-
ica software was used. In the numerical simulations of fluid flow and heat transfer in
thermoplates, the commercial software StarCD was employed. The calculations have
largely been performed on the PC2Paderborn center for parallel computing.
2
2 State of the Art
2.1 Thermoplates
A thermoplate considered in this work consists of two metallic sheets of the same thick-
ness, which are spot-welded over the whole surface area according to an appropriate
pattern, whereas the edges - except for connecting tubes - are continuously seam-welded.
By applying a hydro-form technique, a channel having a complex shape is established
between the sheets, Figure 2.1. One fluid (working fluid) is conducted through this
channel, the other one (process fluid) through the channel confined by two neighbouring
thermoplates.
Figure 2.1: Fluid flow arrangement in thermoplates (DEG-Engineering [64]).
Different types of thermoplates were defined by Behrend [3], Figure 2.2. If the
metallic sheets are of the same thickness, the so-called double embossed thermoplate
is generated. The single embossed thermoplate is developed using the metallic sheets
of different thicknesses, where only the thinner sheet is deformed due to hydroforming,
while the other, thicker one, remains even. More than two sheets can also be used
thereby producing thermoplates of very complex geometry, types triple and quadro. In
the case that the firm sheets are to be used, or have to be thicker for some reasons,
the hydro-form technique is not suitable anymore, and the plates are to be pre-formed,
i.e. dimpled, before they are welded. Depending on its shape and use, thermoplates are
named by some producers: immersion plates, clamp-on plates, pillow plates, plate coils
and so on.
3
2 State of the Art
Figure 2.2: Different types of thermoplates, Behrend [3].
A number of thermoplates can be assembled in parallel at certain spacing thus
making a thermoplate bank assembly or a thermoplate heat exchanger, Figure 2.3.
Such heat exchangers are encountered in several branches of engineering practice,
e.g. as condensers or evaporators in thermal process technology and cooling technique
and in recirculating systems where automatically controlled temperatures are needed.
They are also used in any industries, where cooling or heating in a restricted space is
required, such as: dairy industry, chemical industry, wine or chocolate industry, etc.
Because of very low pressure drop of the fluid, usually vapour, that flows around the
thermoplates, their use as head condensers of distillation plants is especially suitable,
Figure 2.3: Thermoplate assembly (OMEGA heat transfer technology [66]).
4
2.2 Parallel plate channel
see the contribution of M¨uhlthaler [39], Figure 2.4.
Figure 2.4: The use of thermoplate heat exchangers as a head condenser, M¨uhlthaler
[39].
Both fluid spaces in the heat exchanger of this type are completely separated. The
thermoplate heat exchangers, by some manufacturers also called plate-and-shell heat
exchangers, are similar to the plate-and-frame heat exchangers, Figure 2.5, by name,
but from the constructional point of view they are more similar to the shell-and-tube
heat exchangers. Unlike other plate heat exchangers (Figure 2.5), the plates of the
thermoplate heat exchangers do not touch each other, they are self-supporting and do
not cause any forces, that have to be carried by the next plate or housing, Figure 2.3.
Because of fully welded construction there are usually no gasket problems, which makes
this type of heat exchangers suitable for high pressure applications. In comparison to
shell-and-tube heat exchangers, their installation is relatively simple, and the periphery
is drastically reduced for the same heat exchanger duty.
One big disadvantage of thermoplates lies in the fact that their thermo-fluid char-
acteristic is scarcely explored, which prevents a reliable design of the corresponding
heat exchangers. All literature contributions concerning thermoplates can be classi-
fied into three main groups: parallel plate channel, two-dimensional models and three-
dimensional models.
2.2 Parallel plate channel
A parallel plate channel, and from some points of view, its equivalents, duct flow and
annular gap, can be considered as the simplest thermoplate having even walls. Although
a fabrication of such a channel without welding spots seems scarcely possible in sizes
5
2 State of the Art
Figure 2.5: Plate-and-frame heat exchanger, (GEA FlatPlate Inc. [65]).
required in common practice, the parallel plate channel will serve in this work as a
model for comparison with thermoplates. The fluid flow in this case is assumed to be a
single-phase laminar forced convection.
Figure 2.6: Fluid flow through a parallel plate channel.
The interest in this classical problem originates from an article by Graetz [15], who
assumed a Poiseuille flow through a duct and a constant wall temperature and found
an analytical solution of this heat transfer situation that was later called the symmetric
Graetz problem. There are several monographs dealing with various formulations of the
classical Graetz problem, e.g., by Petukhov [50], Shah and London [59], and Kaka¸c et
al [23]. As can be taken from these monographs, the plane asymmetric Graetz problem
when one plate is isothermal and the other one adiabatic was the topic of many papers.
Both thermally developed and developing flows were analysed and positive and negative
Nusselt numbers were observed. Such conditions, however, are rarely encountered in
common practice.
There are only few papers dealing with a developing laminar flow bordered by
walls of constant, but different temperatures, what is called the asymmetric Graetz
6
2.3 Two-dimensional wavy channel
Figure 2.7: Distribution of Nusselt number allong a parallel plate channel, Hatton and
Turton [19].
problem in this work, Figure 2.6. Hatton and Turton [19] solved the energy equation by
separation of variables for eigenvalues in the thermally developing region of the laminar
flow through a parallel plate channel, where the walls are kept at different temperatures,
with the Poiseuille’s velocity field. In the paper, it is shown that the Nusselt number
of one plate can jump from infinite positive to infinite negative becoming zero further
downstream and again positive, Figure 2.7. However, the authors did not show how
the positions where the Nusselt number is infinite and zero depend on the boundary
conditions. Furthermore, neither the conditions for this behaviour to establish nor the
case of a small thermal asymmetry have been investigated. Such data are unavoidable for
calculating the heat transfer under these conditions. A detailed study of the asymmetric
Graetz problem which includes the effects, that have not been analysed in the literature
yet, will be shown in Chapter 3. Some important details of that study were published
by Mitrovi´c, Maleti´c and Baˇcli´c [32].
The same heat transfer situation in an annular gap was treated numerically by
Mitrovi´c and Maleti´c [33], where not only the boundary conditions in form of the fluid
and wall temperatures, but also the radii ratio of the annulus has a big influence on the
heat transfer.
2.3 Two-dimensional wavy channel
A considerable amount of contributions on methods that improve the heat and mass
transfer in compact exchange devices can be found in the literature. The objective of
these methods is mostly to interrupt the boundary layer formed on the wall surface by
bringing the fresh fluid from the core into the near wall region, thereby increasing the
temperature and concentration gradients at the wall. Some examples of such techniques
are off-set fins, louvers, vortex generators, communicating channels, acoustic excitation
7
2 State of the Art
of the flow, oscillatory inflow, use of wavy wall channels etc. The optimal solution
is the one with the lowest pressure drop but the largest heat and mass transfer rate.
However, also other criteria such as ease of production and scaling have to be considered
in particular cases [70].
Channels with wavy walls are easy to fabricate and provide significant heat and
mass transfer enhancement if operated in an appropriate range of the Reynolds number.
Therefore, the wavy wall channels and pipes, Figure 2.8, have been considered in several
studies as means of the heat and mass transfer enhancement, both numerically and
experimentally. It is observed, that the heat transfer enhancement can not be achieved
at low Reynolds numbers. On the other hand, if the flow is made unsteady, either
through external forces or through natural transition to an unsteady state, a significant
increase of heat transfer is possible [70].
The flow of a viscous fluid through axially symmetric pipes and symmetrical wavy
channels under Stokes flow conditions was first treated analytically by Burns and Parks
[6]. The solution, used to calculate the flux through the channel wall for a given wave
velocity or pressure gradient, was obtained by expressing the stream function as a Fourier
series. Using a simple perturbation approach and collocation method, Selvarajan et al.
[57] developed a numerical procedure for the analysis of the incompressible flow through
a channel whose walls describe a travelling wave motion. The flow asymmetry was
observed even for very low Reynolds numbers. Deiber and Schowalter [10] analysed flow
of an incompressible Newtonian fluid through a wavy walled tube. Good agreement
between their numerical and experimental predictions of friction factor as a function of
Reynolds number was achieved. Similar results are given by Ralph [51].
Flow characteristics in a symmetrical channel with wavy walls were investigated
both numerically and experimentally by Nisimura et al. [46]. Up to Re = 350 the
symmetrical flow was observed. Stable twin vortices are formed at the maximum cross
section, where the convective mixing between the mainstream and the vortex hardly
occurs and the separation and reattachment points are fixed. The situation is com-
pletely different at larger Reynolds numbers, where the mainstream is disturbed by the
unsteady vortex motion. For the analysed channel geometry, the authors noted that
laminar flow exists at Reynolds numbers less than 350, and an increase of the Reynolds
number causes turbulent flow to develop. Separation characteristics in a channel of the
same geometry studied Mahmud et al. [28] numerically. A critical Reynolds number,
from which the streamlines start to separate from the wavy wall surface for the given
geometrical parameters, was proposed. According to their results, the critical Reynolds
number falls rapidly with the increase of the surface waviness. Some profiles of the
axial velocity component for two different Reynolds number are qualitative compared
in Figure 2.9. At Re = 200 velocity profiles with negative gradient are observed, wich
indicates a flow reversal. Kim [24] visualised the flow through the wavy channel using
PIV (Particle Image Velocimetry) method, Figure 2.10.
Nishimura et al. [42] investigated experimentally the flow patterns and mass transfer
characteristics in symmetrical two-dimensional channels with sinusoidal and arc-shaped
8
2.3 Two-dimensional wavy channel
Figure 2.8: Different types of wavy channels.
walls, Figure 2.8a and b. The flow through arc-shaped channels is characterised by flow
separation in each furrow, which leads to an earlier transition from laminar to turbu-
lent flow in comparison to wavy channels. The flow separation implies a renewal of the
boundary layers thus yielding the mass transfer enhancement. The arc-shaped wall has
a larger mass transfer rate than the sinusoidal one because of an earlier transition to tur-
bulence. Although the both channel geometries, the sinusoidal and the arc-shaped one,
are taken to be two-dimensional, a regular three-dimensional flow is observed. Niˇceno
and Nobile [40] investigated numerically fluid flow and heat transfer in channels of the
same geometries as used by Nishimura et al. [42] and confirmed their experimental
results. Using a time-accurate, unstructured covolume method, they found that the
heat transfer augmentation was only possible in unsteady regimes. Compared to the
sinusoidally shaped channel, the heat transfer increase for the arc-shaped passage was
9
2 State of the Art
Figure 2.9: Axial velocity profiles at two different Reynolds numbers for a/Hmax =
0.175 and L/Hmax = 1.5, Mahmud et al. [28].
higher. However, the friction factor was also higher than that of the sinusoidal channel.
Reports on similar experiments in a sinusoidal wavy walled tube were made in a further
paper by Nishimura et al. [43]. As they found, the laminar-like flow and turbulent-like
flow alternatively appear in transitional flow regime, indicating an intermittent flow be-
haviour. The increment in mass transfer in this flow regime was larger for the wavy tube
than for the wavy channel.
Stone and Vanka [63] analysed a wavy channel consisting of 14 waves. They solved
the Navier-Stokes equations governing a time-dependent flow and found out that the
flow was steady for low Reynolds numbers, but increasing it beyond a modest value,
Figure 2.10: Mean streamlines at Re = 500 visualised by using PIV method, Kim [24].
10
2.3 Two-dimensional wavy channel
the flow became unsteady. In this case the increased mixing between the core and
near-wall fluid resulted in enhanced heat transfer and pressure drop. Hossain and Islam
[21] studied sine-shaped wavy channels numerically by solving two-dimensional Navier-
Stokes and energy equations implementing periodic boundary conditions. They found
that the flow was steady up to a critical Reynolds number depending on the channel
geometry, which confirms the results of Stone and Vanka [63]. Their results show further
that decreasing channel height and increasing amplitude cause the flow to become more
unstable resulting in higher heat transfer and pressure drop, but variation of wavelength
has weak effect. The best heat transfer enhancement is achieved in the transitional
regime. Furthermore Hossain and Islam [22] investigated numericaly in the same way
the channels with sinusoidal, triangular and arc-shaped walls, Figure 2.8a, c and b,
and compared them. Their results show that in the arc-shaped channel the transition
between laminar and turbulent flow starts at a relatively low Reynolds number, for the
sine-shaped channel it is relatively high, and for triangular channel it is in between
these two. The comparison of all three channels was not possible in the same region
of the Reynolds number, since the solutions for arc and triangular channels become
divergent after certain Reynolds number, with the code they used. Similar study of the
converging-diverging tube conducted Sparrow and Prata [61]. They adopted periodic
boundary conditions in numerical experiments and achieved a very good agreement
with experiments. In comparison to the straight tube, the heat transfer enhancement
in corrugated tube was possible only for Prandtl numbers greater than 2.5. Wirtz et al.
[72] investigated air flow though a grooved channel experimentally. The Nusselt number
of the windward face was approximately twice that of the leeward face.
Russ and Beer [53] investigated a pipe with sinusoidal wavy surface, both numer-
ically and experimentally in laminar, transitional and turbulent flow regime. The nu-
merical and the experimental results are in good agreement, except for the transitional
flow regime. The positions of the minima and maxima of the Sherwood and Nusselt
numbers nearly coincide with the flow separation and the reattachment point, respec-
tively. Saniei and Dini [55] conducted an experimental study of convective heat transfer
in the turbulent regime in a wavy channel. The highest magnitude of the local Nus-
selt number was on the second wave and its local maxima were located upstream of
the peak of each wave. The average Nusselt number of the wavy channel was always
higher than the one of the even channel for the Reynolds number range. Wang and
Vanka [70] studied the flow patterns and heat transfer characteristics in a sinusoidally
curved converging-diverging channel numerically. Implementing a curvilinear orthogo-
nal grid they got accurate numerical solutions. For Reynolds numbers below 180 steady
laminar flow was established, while beyond this value a transition to chaotic flow with
significant increase in heat transfer was observed. Laminar flow and heat transfer in a
sinusoidally curved converging-diverging channel were analysed by Garg and Maji [14]
using a finite difference method for solving the conservation equations. They show that
the flow separates mainly in the diverging region. However, for higher Reynolds num-
ber and amplitude-wavelength ratio the flow separates in the converging region as well.
The pressure distribution in the axial direction is a combination of a nearly linear and
11
2 State of the Art
periodic function. The pressure drop in the fully developed region is constant for each
cycle and increases with Reynolds number and amlitude-wavelength ratio. The Nusselt
number increases with Reynolds number as well, which is not the case in a parallel plate
channel. The analysis of the same physical problem, but inside a sinusoidal wavy tube
was conducted by Mahmud et al. [28] using the Finite-Volume method. They propose
a correlation for calculating the friction factor as a function of surface waviness and
Reynolds number. The surface waviness was formed to cause earlier flow separation,
larger pressure drop and higher heat transfer rate. Heat transfer and entropy generation
in a wavy channel under the influence of inclination was studied by Mahmud and Islam
[29] for the case of laminar free convection and the walls having different temperature.
The distributions of the average Nusselt number and of the average entropy generation
number have bin given. However, the angle of inclination is not expected to have any
essential effect on forced convection fluid flow.
Turbulent flow through the wavy wall channel and convective heat transfer were in-
vestigated by Dellil et al. [11]. As they reported, the two-layer model they implemented
(the standard k-εturbulence model combined with a one-equation model to resolve the
near-wall viscosity-affected regions) was successful in capturing most of the important
physical features of a turbulent flow with reasonable amount of memory storage and
computer time. The results of their simulations show that the minimum and the max-
imum Nusselt numbers are located near the separation and the reattachment points,
respectively, which supports the results of Russ and Beer [53]. The average Nusselt
number increases with the wave amplitude until a critical value, after which it remains
largely constant. However, the heat transfer enhancement is accompanied by an increase
in the pressure drop.
Figure 2.11: Distribution of the local Nusselt number for Re = 500 and Pr = 6.93,
[69].
A very comprehensive and detailed study of a laminar forced convection in a wavy
wall channel is given in a paper by Wang and Chen [69]. A simple coordinate trans-
formation method was employed to transform the wavy channel into a parallel plate
12
2.3 Two-dimensional wavy channel
Figure 2.12: Distribution of the average Nusselt number for Re = 500 and Pr = 6.93,
[69].
channel. Effects of the waviness, Reynolds number and Prandtl number on the local
friction coefficient and Nusselt number were investigated. Their results show that the
channel forms a highly complex flow pattern, which comprises a strong forward flow and
an oppositely directed recirculating flow with each wave. The harmonic distribution
of the Nusselt number has the same frequency as the wavy surface, Figure 2.11. The
maxima occur in the converging section of each wave near the minimal cross-section of
the channel and its positions are independent of the amplitude-wavelength ratio and the
Reynolds number, Figures 2.11 to 2.13. The minima of the local Nusselt number are
located downstream and upstream within a short distance of the largest cross-section
of each wave. With an increase in the Reynolds number, the maxima of the Nusselt
number and of the friction coefficient increase in each converging channel position, but
its minima remain almost constant in the diverging section. Similar distributions of the
Sherwood number obtained Nishimura et al. [47].
Figure 2.13: Distribution of the local Nusselt number for a/L = 0.2 and Pr = 6.93,
[69].
13
2 State of the Art
Detailed direct numerical simulations of laminar, transitional and chaotic flow in
converging - diverging wavy channels were performed by Guzm´an and Amon [18]. Mod-
ern dynamical system techniques such as time-delay reconstructions of pseudophase
spaces, autocorrelation functions, fractal dimensions and Eulerian-Lyapunov exponents
are used for the dynamical flow characterisation and stability analysis. Their simu-
lations indicate the presence of stable three-dimensional oscillatory flows. The value
of the largest Eulerian Lyapunov exponent verifies that the flow changes its behaviour
from a quasi-periodic predictable regime to an aperiodic, chaotic, unpredictable regime
at Reynolds numbers between 500 and 550. In Figure 2.14 are shown the periodically
fully developed streamlines obtained numerically by Guzm´an and Amon [18] for differ-
ent Reynolds numbers. An analysis of linear stability of two-dimensional steady flow in
Figure 2.14: Streamlines predicted numerically by Guzm´an and Amon [18] for different
Reynolds numbers.
wavy channels as in Figure 2.8a and 2.8d with ϕ= 180◦, was conducted by Cho et. al
[9]. Their numerical investigation reveals that the two-dimensional steady flow becomes
unstable depending on the waviness hight. Critical Reynolds numbers of the channel
as in Figure 2.2d with ϕ= 180◦are smaller for three-dimensional disturbances and
larger for two-dimensional disturbances in comparison to the symetric channel, Figure
2.2a, and in any case substantially lower than that for plane channel. Another linear
stability of the flow in wavy channels analysed was investigated by Blancher et. al [5],
who applied a Galerkin approach divergence-free Chebychev basis functions and trigono-
metric polynomials. They found the instability to set in as a Tollmien-Schlichting wave
at Re ≈90. It is also shown that the variations of the local heat transfer coefficient
are increased with increasing Prandtl number. Lahbabi and Chang [25] used a global
Galerkin method to solve the Navier-Stokes equation in wavy tubes. Their stationary
14
2.4 Three-dimensional wavy channel
solution is stable to axisymmetric disturbance at all Reynolds numbers, but destabilises
to azimuthal disturbance at a critical Reynolds number. A nonlinear analysis was carried
out to characterise the Hopf bifurcation.
An experimental study of flow and heat transfer in sinusoidal wavy passages with
different phase shift ϕbetween the sidewalls was conducted by Rush et al. [52], Figure
2.8d. The two walls are in-line with each other at ϕ= 0◦. At the other extreme, ϕ=
180◦, converging-diverging wavy channel is established, Figure 2.8a. The experimental
results confirm that flow instabilities cause a heat transfer enhancement. The fluid flow
is unstable and unsteady even in the laminar flow regime. For a given Reynolds number,
the fluid flow in the channels with ϕ= 0◦and ϕ= 90◦is unsteady farther upstream
than in case of ϕ= 180◦. However, the pressure drop data were not recorded in this
study.
As a further means of the heat and mass transfer enhancement, the influence of
the oscillatory flow in symmetric wavy walled channels was investigated experimentally
by Nishimura et al. [44], Nishimura et al [45] and Nishimura and Matsune [48] and
numerically by Nishimura and Matsune [49] and Lee et al. [26]. The influence of the
Reynolds number and Strouhal number on the heat and mass transfer has been found
to be very important.
Although the works mentioned above are performed with two-dimensional flows on
a macroscopic scale, they provide the fundamental background information that is basi-
cally also valid for the numerical experiments performed for channels with a pronounced
three-dimensional flow that will be presented in this work further below.
2.4 Three-dimensional wavy channel
Channels with three-dimensional wall corrugations are of much greater importance for
the practical applications than the two-dimesional wavy cannels, but only few papers
are dealing with fluid flow and heat transfer in such channels.
Witry et al. [73, 74] used the CFD code FLUENT and studied numerically the
fluid flow, pressure drop and heat transfer of both, inside and outside fluid of a dimple
thermoplate, Figure 2.15. Near the inlet and outlet of the thermoplate they analysed,
a number of dimples have been removed to allow the interior flow to re-distribute itself
without causing high pressure losses. The velocity and pressure distributions of the
inside fluid flow are given in Figures 2.16 and 2.17, respectively. Closer inspection of
the flow field around the dimples show the three-dimensionality of the problem. It can
also be seen, that the high inlet pressure is lost due to the direct hit against the dimple
facing the inlet and to the sudden enlargement in cross-sectional area, which indicates
the need for more dimples to be removed from the region near the inlet connection tube.
A rapid temperature drop, Figure 2.18, takes place in the vicinity of areas where the
flow slows down. The high fluid velocities, especially near the inlet, lead to a high heat
15
2 State of the Art
Figure 2.15: Dimple plate geometry numerically analysed by Witry et al. [73].
transfer coefficients. The values of the overall heat transfer coefficient are higher than
the expected ones for tubular heat exchangers. However, the authors did not make
any numerical experiment to show how the heat transfer depends on the thermoplate
geometry.
Sawyers et al. [56] combined analytical and numerical techniques to study the effect
of three-dimensional hydrodynamics on the enhancement of steady, laminar heat transfer
in corrugated channels with walls having constant, but different temperature. Since the
periodic boundary conditions have been adopted, the results refer to the developed flow.
The authors took the channel with two-dimensional sinusoidal corrugations, Figure 2.19,
as the base. It was found, that the heat transfer in such channels was higher than for
the flat plates due to the presence of recirculation zones. In the three-dimensional case,
Figure 2.20, the corrugations are sinusoidal in two orthogonal directions. A small mean
flow in the transverse direction allows particles to cross between the recirculation zones
and the main flow, which leads to an increase in the heat transfer.
From the literature review follows that the fluid flow with heat transfer in channels
with three-dimensional wall corrugations has benn studied only by few authors. A
dependency of heat transfer and pressure drop on geometrical parameter of such channels
was not investigated. A Nusselt number correlation has not been given. For that reason
in this thesis the fluid flow and heat transfer in the thermoplates with walls corrugated
sinusoidaly in two orthogonal directions that are spot-welded according to an appropriate
pattern were studied. The walls were kept at the same, constant temperature. The
dependency of the heat transfer and the pressure drop on thermoplate geometry was
investigated. The details will be given in Chapters 4 and 5. A part of these investigations
were published in papers by Mitrovi´c and Maleti´c in their papers [34, 35, 36].
16
2.4 Three-dimensional wavy channel
Figure 2.16: Internal flow velocity vectors in the plate radiator calculated by Witry et
al. [73].
Figure 2.17: Internal flow pressure distribution in the plate radiator calculated by
Witry et al. [73].
Figure 2.18: Internal flow temperature distribution in the plate radiator calculated by
Witry et al. [73].
17
3 Heat Transfer in a Parallel Plate
Channel
A parallel plate channel is a simple device that basically belongs to the large chan-
nel family. It is frequently used for comparison purposes when channels of a complex
geometry are to be used for numerical or laboratory experiments. For this reason, the
fluid flow and heat transfer in this channel is treated next under steady-state conditions,
assuming the plates to be at constant, but different temperatures.
3.1 Physical model and governing equations
The physical model of the so-called asymmetric Graetz problem is shown in Figure 3.1.
Two smooth parallel plates that confine the fluid space, a distance 2bapart, are kept at
constant but different temperatures ϑW1and ϑW2. Since the plates are infinite in the
y-direction (i.e. a >> 2b), this represents a two-dimensional model. The steady state
laminar flow is directed along the z-coordinate.
Figure 3.1: Illustration of the physical model.
As given in the monograph of Petukhov [50], the fluid velocity w(x, z) along the
z-axis, that takes into account the hydraulically developing region is calculated from an
expression deduced by Targ,
w(X, Z) = 3
2¯w(1 −X2)−2 ¯w
∞
X
n=1
1
ω2
n1−cos(ωnX)
cos(ωn)e−16ω2
nZ,(3.1)
19
3 Heat Transfer in a Parallel Plate Channel
where the dimensionless parameters are defined in the following manner
X=x
b, Z =1
Re
z
dh
, dh= 4b, Re =¯wdh
ν= 4 ¯wb
ν,(3.2)
ωnbeing the nth root of the equation tan(ω) = ω.
In case of hydrodynamically developed flow, Equation (3.1) reduces to the Poiseuille
parabolic velocity field
w(X) = 3
2¯w(1 −X2).(3.3)
For simplicity reasons, in many textbooks a constant velocity is adopted:
w(X) = ¯w. (3.4)
To describe the heat transfer between parallel plates mathematically , the energy
equation for stationary single phase forced two-dimensional convection is used
u∂ϑ
∂x +w∂ϑ
∂z =κ∂2ϑ
∂x2+∂2ϑ
∂z2,(3.5)
where κrepresents the thermal diffusivity; the viscous dissipation is neglected. Since
the distance between the plates is considered to be small in comparison to its length,
the convection in the x-direction and the diffusion in the z-direction can be neglected.
The energy Equation (3.5) becomes:
w∂ϑ
∂z =κ∂2ϑ
∂x2.(3.6)
According to the Figure 3.1, the temperature ϑmust satisfy the following boundary
conditions:
−b≤x≤+band z= 0 : ϑ(x, 0) = ϑIN ,(3.7)
x=−band z > 0 : ϑ(−b, z) = ϑW1,(3.8)
x= +band z > 0 : ϑ(+b, z) = ϑW2.(3.9)
3.2 Normalised temperature function and analytical
solution
Using the expression for the hydrodynamically developed flow (3.3) and introducing the
dimensionless quantities
Θ(X, Z∗) = ϑ(X, Z∗)−ϑW1
ϑIN −ϑW1
,(3.10)
Z∗=8
3
1
RePr
z
b, Pr =ν
κ,(3.11)
20
3.2 Normalised temperature function and analytical solution
the energy equation (3.6) becomes
(1 −X2)∂Θ
∂Z∗
=∂2Θ
∂X2.(3.12)
The boundary conditions (3.7)-(3.9) become correspondingly
−1≤X≤+1 and Z∗= 0 : Θ(X, 0) = 1,(3.13)
X=−1 and Z∗>0 : Θ(−1, Z∗) = ΘW1= 0,(3.14)
X= +1 and Z∗>0 : Θ(+1, Z∗) = ΘW2.(3.15)
To separate the variables in Xand Z∗directions for Equation (3.12) and to satisfy
the boundary conditions (3.13)-(3.15) the following ansatz was made
Θ(X, Z∗) = ΘW2
2(1 + X) +
∞
X
n=1
Anϕn(X)e−λ2
nZ∗,(3.16)
where the first term on the right-hand side represents the temperature profile of the
thermally developed flow, ϕn(X) are the eigenfunctions and λnthe eigenvalues of the
following eigenvalue problem:
ϕ00
n(X) + λ2
n(1 −X2)ϕn(X) = 0,(3.17)
ϕn(−1) = 0,(3.18)
ϕn(+1) = 0.(3.19)
Since all ϕnfunctions are orthogonal to each other, the scalar product R1
−1ϕnϕm
disappears, except for n=m. Multiplying (3.16) by ϕm(1 −X2), integrating from
X=−1 to X= +1 and setting back m=nthe expression for the Anis deduced:
An=(ΘW2−1)ϕ0
n(+1) + ϕ0
n(−1)
λnϕ0
n
∂ϕn
∂λnX=+1 −ϕ0
n
∂ϕn
∂λnX=−1.(3.20)
A solution of (3.17) can be obtained in form of series. However, the solution can
be written in more compact form using the Kummer confluent hypergeometric func-
tion 1F1((1 −λn)/4; 1/2; λnX2) and the Hermite polynomials H(λn−1)/2(λ1/2
nX), that are
already implemented in a number of commercial software packages (in this work Math-
ematica 5 was used) and that can be found, e.g., in books by Slater [60] or Abramowitz
and Stegun [1]. The solution of (3.17) has thus the following form
ϕn(X) = Dne−λnX2
21F11−λn
4;1
2;λnX2+EnHλn−1
2λ
1
2
nX,(3.21)
21
3 Heat Transfer in a Parallel Plate Channel
where Dnand Enare arbitrary constants. The constants Dncan be thought to be a
part of An, thus Dn= 1. Applying the boundary conditions (3.18) and (3.19) to the
solution (3.21) gives:
1F11−λn
4;1
2;λn+EnHλn−1
2−λ
1
2
n= 0,(3.22)
1F11−λn
4;1
2;λn+EnHλn−1
2λ
1
2
n= 0.(3.23)
These expressions suffice to obtain the constants Enand the eigenvalues λnas positive
roots of the characteristic equation, which follows from (3.22) and (3.22), and which can
be written as
1F11−λn
4;1
2;λnHλn−1
2λ
1
2
n−Hλn−1
2−λ
1
2
n= 0.(3.24)
It is useful and illustrative to rewrite the characteristic equation (3.24) as follows
1F11−λn
4;1
2;λn= 0,(3.25)
Hλn−1
2λ
1
2
n−Hλn−1
2−λ
1
2
n= 0,(3.26)
because it can be recognised that Equation (3.25) corresponds to the classical, symmet-
rical Graetz problem, whereas Equation (3.26) takes into account the asymmetry of the
thermal action on the fluid. The first 30 eigenvalues obtained from (3.25) and (3.26) are
listed in the Table 3.1.
Hatton and Turton [19] solved also this problem for eigenvalues, which are in a
satisfactory agreement with those from this work, but the solution for the eigenfunctions
ϕn(X) in their paper was not shown. This extension, Equation (3.21), was done by
Mitrovi´c, Maleti´c and Baˇcli´c [32]. For n= 3k+ 1, k= 0,1,2, ..., ∞, the eigenfunctions
are zero, ϕn(X) = 0, so that only the roots λnfor n= 3k+ 2 and n= 3k+ 3,
k= 0,1,2, ..., ∞have to be considered.
In the general solution (3.21) the eigenfunctions ϕncorresponding to the eigenvalues
λnwhere n= 3k+ 2 are called by Hutton and Turton [19] even eigenfunctions and they
have the following properties:
ϕn(−1) = ϕn(1) = 0, ϕn(0) = 1, ϕ0
n(−1) = −ϕ0
n(1), ϕ0
n(0) = 0,(3.27)
whereas for n= 3k+ 3 the eigenfunctions, called odd eigenfunctions by Hutton and
Turton [19], have the properties that are different from those given in Equation (3.27)
ϕn(−1) = ϕn(1) = ϕn(0) = 0, ϕ0
n(−1) = ϕ0
n(1), ϕ0
n(0) >0.(3.28)
22
3.2 Normalised temperature function and analytical solution
Table 3.1: Some eigenvalues λnof the asymmetrical Graetz problem
n λnn λnn λn
1 1.00000 11 13.668 21 27.667
2 1.68169 12 15.667 22 29.000
3 3.67229 13 17.000 23 29.667
4 5.00000 14 17.667 24 31.667
5 5.66996 15 19.667 25 33.000
6 7.66881 16 21.000 26 33.667
7 9.00000 17 21.667 27 35.667
8 9.66824 18 23.667 28 37.000
9 11.6679 19 25.000 29 37.667
10 13.0000 20 25.667 30 39.667
3.2.1 Characteristic temperature distributions
Before discussing the thermal asymmetry effect on the Nusselt number, i.e. heat transfer,
some typical temperature distributions are presented. Figure 3.2 shows, as an example,
the development of the temperature field when the fluid enters the channel at a tem-
perature below those of the plates, ϑIN < ϑW2< ϑW1. The surface Θ = Θ(X, Z∗),
Equation (3.16), bordering the thermally affected fluid region is attached to the plates
having the temperatures ΘW1= 0 and ΘW2= 0.4. As may be taken from the Figure,
the fluid temperature decreases along the coordinate Z∗, except for the immediate start
of the thermal action on the fluid near X= 0.
At fixed Z∗, the temperature profile shows a maximum, which moves toward the
upper plate (X= +1, ΘW2= 0.4), as Z∗increases, Figure 3.3. When the temperature
maximum reaches this plate, the plate becomes adiabatic with a zero heat flux at that
position. For the chosen temperature ΘW2, this occurs at Z∗= 0.6816, Figures 3.2 and
3.3.
The axial position Z∗with the zero heat flux is an important marker of the phys-
ical model. Regarding the heat transfer, the corresponding cross-section separates the
channel in two regions. In the upstream region, the fluid is heated by both plates, but
downstream, it is heated by the lower (X=−1, ΘW1= 0), and cooled by the upper plate
(X= +1, ΘW1= 1). This interplay is illustrated in Figure 3.4, where the slopes of the
dimensionless temperature at the plates are plotted along the coordinate Z∗according
to Equation (3.39). The curves correspond to the wall heat flux ˙qWgiven in equations
(3.30) and (3.31). For equal temperature slopes, (∂Θ/∂X)W1= (∂Θ/∂X)W2= 0.2, the
heat is transferred at the same rate from the lower plate to the upper one through the
fluid like through a solid wall, which corresponds to the thermally developed flow, see
the temperature profile at Z∗= 2 in Figure 3.3.
As a further example, Figure 3.5 gives the temperature development for ΘW2= 1,
23
3 Heat Transfer in a Parallel Plate Channel
Figure 3.2: Temperature field at specified boundary conditions, ΘW1= 0 and ΘW2=
0.4, with hydrodynamically developed flow, Equation (3.3).
Figure 3.3: Temperature profiles for some fixed values of Z∗. At Z∗= 0.6816 the upper
plate is adiabatic.
Figure 3.4: The temperature slope at the plates for ΘW1= 0 and ΘW2= 0.4.
24
3.2 Normalised temperature function and analytical solution
Figure 3.5: Temperature field at specified boundary conditions, ΘW1= 0 and ΘW2= 1,
with hydrodynamically developed flow, Equation (3.3).
Figure 3.6: Temperature profiles for some fixed values of Z∗for ΘW1= 0 and ΘW2= 1.
Figure 3.7: The temperature slope at the plates for ΘW1= 0 and ΘW2= 1.
25
3 Heat Transfer in a Parallel Plate Channel
Figure 3.8: Temperature field at specified boundary conditions, ΘW1= 0 and ΘW2= 2,
with hydrodynamically developed flow, Equation (3.3).
Figure 3.9: Temperature profiles for some fixed values of Z∗for ΘW1= 0 and ΘW2= 2.
Figure 3.10: The temperature slope at the plates for ΘW1= 0 and ΘW2= 1.
26
3.3 The Nusselt Number
that is, for the temperature of one plate equal to the fluid inlet temperature. The
corresponding temperature profiles and the temperature slopes are plotted in Figures
3.6 and 3.7.
Finally, Figure 3.8 shows the development of the temperature field for ΘW2= 2.
In this case, the fluid is heated by one and cooled by the other plate at the same
rate so that the mean fluid temperature remains unchanged along the flow direction,
¯
Θ(Z∗) = 1. Some temperature profiles and the temperature slopes at the plates are
displayed in Figures 3.9 and 3.10. Because of the total thermal asymmetry imposed, the
temperature slopes coincide for this case.
3.3 The Nusselt Number
The local Nusselt number Nu(z) is usually defined by using the local heat transfer
coefficient, h(z), the hydraulic diameter of the flow cross section, dh, and the thermal
conductivity of the fluid, k, as follows:
Nu(z) = h(z)dh
k.(3.29)
Since the plates are kept at different temperatures, the heat flux, the heat transfer
coefficients and the Nusselt numbers at the upper and at the lower plate are consequently
different.
At this place, it is to be noted, that some authors define an average Nusselt number
of the channel, that takes into account the heat transfer at both plates, see e.g. Nield [41].
However, such average Nusselt number has in case of thermal asymmetry no practical
relevance. For this reason, it is clearly distinguished between heat transfer characteristics
at the upper and at the lower plate in this work.
Having in mind the boundary conditions (3.8) and (3.9), i.e. the physical model in
Figure 3.1 and adopting the sign convention to give the heat flux direction according to
the temperature gradient (direction of decreasing temperature), the heat flux from the
lower and the upper plate, respectively, to the fluid is given by
˙qW1(z) = −k∂ϑ
∂xx=−b
=hW1(z)ϑW1−¯
ϑ(z),(3.30)
˙qW2(z) = −k∂ϑ
∂xx=+b
=hW2(z)¯
ϑ(z)−ϑW2.(3.31)
From these equations the heat transfer coefficients at the plates are immediately
27
3 Heat Transfer in a Parallel Plate Channel
calculated by
hW1(z) = −k
ϑW1−¯
ϑ(z)∂ϑ(x, z)
∂x x=−b
,(3.32)
hW2(z) = k
ϑW2−¯
ϑ(z)∂ϑ(x, z)
∂x x=+b
.(3.33)
The average (mixing-cup) fluid temperature is calculated from
¯
ϑ(z) = 1
2b¯wZ+b
−b
w(x, z)ϑ(x, z)dx. (3.34)
where the Poiseuille parabolic velocity profile given in Equation (3.3) was used.
With the dimensionless temperature defined in Equation (3.10), the heat transfer
coefficients (3.32) and (3.33) become
hW1(Z∗) = 1
¯
Θ(Z∗)
k
b∂Θ(X, Z∗)
∂X X=−1
,(3.35)
hW2(Z∗) = 1
ΘW2−¯
Θ(Z∗)
k
b∂Θ(X, Z∗)
∂X X=+1
,(3.36)
where the average dimensionless temperature and temperature slope are
¯
Θ(Z∗) = ¯
ϑ(Z∗)−ϑW1
ϑIN −ϑW1
=3
4Z+1
−1
(1 −X2)Θ(X, Z∗)dX, (3.37)
¯
Θ(Z∗) = ΘW2
2−3
4
∞
X
n=1
An
λ2
n
(ϕ0
n(+1) −ϕ0
n(−1)) e−λ2
nZ∗,(3.38)
∂Θ(X, Z∗)
∂X =ΘW2
2+
∞
X
n=1
Anϕ0
n(X)e−λ2
nZ∗.(3.39)
The corresponding Nusselt numbers are then
NuW1(Z∗) = 4
¯
Θ(Z∗)∂Θ(X, Z∗)
∂X X=−1
,(3.40)
NuW2(Z∗) = 4
ΘW2−¯
Θ(Z∗)∂Θ(X, Z∗)
∂X X=+1
.(3.41)
In the fully developed region (Z∗→ ∞) the summations in equations (3.38) and
(3.39) disappear and equations (3.40) and (3.41) immediately give
NuW1=NuW1=Nu∞= 4.(3.42)
28
3.3 The Nusselt Number
In comparison with the thermally developed flow, the local Nusselt numbers NuW1(Z∗)
and NuW2(Z∗), calculated from equations (3.40) and (3.41), can dramatically change in
the thermally developing region. Depending on the boundary conditions (3.13) to (3.15)
the average fluid temperature ¯
Θ(Z∗) can become equal to the wall temperature ΘW1= 0
or ΘW2at certain position Z∗. Due to this fact and the definition of the heat transfer
coefficient, the Nusselt Number NuW1or NuW2can become infinite, positive or nega-
tive, or, in other words, the Nusselt numbers can have a vertical asymptote, although
the temperature slope (3.39) changes always continuously.
Furthermore, the temperature slope defined in Equation (3.39) can become equal
to zero at some position Z∗, which characterises the adiabatic cross section at the upper
or at the lower plate.
Figure 3.11: Distribution of the local Nusselt numbers Nu(Z∗) along the flow direction
for ΘW1= 0 and ΘW2= 0.4.
Such an example is illustrated in Figure 3.11 for ΘW2= 0.4, which corresponds
to the temperature field and the temperature profiles shown in Figures 3.2 and 3.3.
While the Nusselt number NuW1(Z∗) of the lower plate (X=−1) is always positive and
continually decreasing along the flow direction, the situation is essentially different with
the Nusselt number NuW2(Z∗) of the upper plate (X= +1). The curve representing
NuW2(Z∗) shows a discontinuity with a vertical asymptote. Left to the asymptote,
NuW2(Z∗)>0, whereas on the right, NuW2(Z∗) is first negative, becoming zero and
then positive. With increasing Z∗, the Nusselt numbers take asymptotically the value in
the thermally developed flow, Nu∞= 4. Similar behaviour has already been reported
by Hatton and Turton [19].
Figure 3.12 gives a general idea of what happens when the dimensionless temper-
ature of the upper plate ΘW2increases from a negative value up to ΘW2= 0 first, and
than further up to ΘW2= 1. The situation where the fluid inlet temperature ϑIN is
greater than both wall temperatures, ϑW1and ϑW2, which corresponds to ΘW2<0, is
displayed on the left side of the Figure 3.12. The Nusselt number NuW2of the upper
plate is always positive and changes continuously, whereas the Nusselt number NuW1
of the lower plate experiences a discontinuity while passing the vertical asymptote. The
29
3 Heat Transfer in a Parallel Plate Channel
Figure 3.12: Effect of boundary conditions on Nusselt number for hydrodynamically
developed flow for ΘW1= 0 at different ΘW2.
30
3.3 The Nusselt Number
Figure 3.13: Effect of boundary conditions on Nusselt number for hydrodynamically
developed flow for ΘW1= 0 at different ΘW2.
31
3 Heat Transfer in a Parallel Plate Channel
position of this asymptote moves downstream as the dimensionless temperature of the
upper plate ΘW2rises in the interval −∞ <ΘW2<0. The situation is very similar if
ΘW2rises further in the interval 0 <ΘW2<1, right side of the Figure 3.12. This is
the case when the fluid inlet temperature ϑIN is less than both wall temperatures, ϑW1
and ϑW2. The upper and the lower plate change the roles from a heat transfer point
of view and now the Nusselt number of the lower plate NuW1is positive and changes
continuously, whereas NuW2shows the discontinuity.
Finally, in the region ΘW2≥1, which is the case if the fluid inlet temperature
ϑIN is between the plate temperatures, ϑW1> ϑIN > ϑW2or ϑW1< ϑIN < ϑW2, both
Nusselt numbers NuW1and NuW2are positive without discontinuities, Figure 3.13. At
ΘW2= 1, the NuW2curve is sandwiched between two horizontal asymptotes, Nu = 0
and Nu = 4, thereby showing an inflection point. For ΘW2→+∞, the behaviour is the
same as for ΘW2= 1, but this time, the Nusselt number NuW1is sandwiched between
the horizontal asymptotes.
For ΘW2= 0 both plates are kept at the same temperature, which corresponds to
the classical symmetric Graetz problem. In this case NuW1and NuW2are identical,
positive, change continuously, becoming Nu∞= 7.545 at Z∗→ ∞. If there are some
heat source terms in the energy equation, the average fluid temperature can exceed the
wall temperature, even if the heating is symmetric. Such an example in case of porous
channel is analysed as special case in the paper by Mitrovi´c and Maleti´c [37].
The asymmetric Graetz problem considered in this work is discontinuous in two
respects. One kind of the discontinuity represents the above discussed jump of the
Nusselt number from positive to negative values. The other kind of discontinuity is the
switch from one developed state to the other one, which means from the symmetric to
the asymmetric thermal conditions. In context with the latter discontinuity, the question
arises: How does the fluid decide at a small thermal asymmetry, ΘW2→0 or ϑW1≈ϑW2,
which one of the limiting values of the Nusselt numbers, Nu∞= 4.0 (asymmetric case) or
Nu∞= 7.545 (symmetric case), is to be followed? An answer to this question is not easy
to formulate, but Figure 3.14 provides a general idea of what happens when the thermal
asymmetry measured in terms of ΘW2is gradually decreased. The curves represent the
Nusselt numbers for different thermal asymmetries, whereas the two horizontal lines give
the limits of the Nusselt numbers in the fully developed heat transfer region for thermal
asymmetry (lower line) and for the thermal symmetry (upper line). At a sufficiently
large asymmetry (greater ΘW2, ΘW2>0 ) the fluid takes resolutely the ”right” way,
Figure 3.14, left side. On the contrary, at a very small asymmetry (ΘW2→0), one gets
the impression as if the fluid knew very little about its fate downstream the flow. The
Nusselt number NuW1first follows the thermal symmetry with Nu∞= 7.545 from below,
but it apparently remembers, further downstream, the actual boundary conditions and
switches down to the limit Nu∞= 4.0 at the thermal asymmetry. At the same time,
the Nusselt number NuW2is either smaller than 4.0 or larger than 7.545. The values
between these limits are inaccessible to NuW2. If ΘW2→0, but ΘW2<0, Figure 3.14,
right side, the Nusselt numbers behave in the same way, only with changed roles in
32
3.3 The Nusselt Number
comparison to the discussed situation, where ΘW2>0.
Figure 3.14: Transitions and asymptotical behaviour of the Nusselt numbers for chosen
ΘW2.
The thermal behaviour illustrated in Figures 3.11-3.14 has a strong impact on the
calculation of the heat transfer. The assumption of a thermally fully developed flow,
frequently adopted in practice, is unallowable in an asymmetric case. This is convincingly
demonstrated in Figure 3.14, illustrating at the same time the complexity expected when
trying to correlate the Nusselt numbers in terms of common quantities. The question
as how to calculate the heat transfer coefficient in such a case becomes important, and
it is by no means as trivial as it might appear. It is also obvious from Figure 3.14 that
the asymmetric Graet problem represents a typical example of physical models where
relatively small causes result in strong effects.
In the analysis presented so far, the hydrodynamically developed flow was adopted,
i.e. the Poiseuille parabolic velocity profile, expression (3.3), was used to find an an-
alytical solution to the energy equation (3.6). Adopting the plug flow (3.4) or the
hydrodynamically developing flow given in Equation (3.1), where the solution to the
energy equation (3.6) can be found only numerically, results in the same qualitative dis-
tributions of the Nusselt numbers as in Figures 3.12 and 3.13 for the hydrodynamically
developed flow. Figure 3.15 displays, as an example, the velocity effect on the Nusselt
number, where, for comparison reasons, all tree mentioned velocity distributions were
considered. In the case of a constant velocity (plug flow, Figure 3.15a), the position
of the vertical asymptote is shifted upstream in comparison to the hydrodynamically
developed flow (Figure 3.15b). The almost perfect agreement of the developing (Figure
3.15c) and developed flows is rooted in the relatively large Prandtl number Pr, which
affects the contribution of the sum in Equation (3.1) to the parabolic term essentially.
This becomes obvious if the coordinate Z∗is introduced in this equation by replacing
Zby (3/32)PrZ∗, as required by Equations (3.2) and (3.11). Even at a much smaller
33
3.3 The Nusselt Number
Prandtl number, the sum in Equation (3.1) is almost negligible for the conditions spec-
ified in Figure 3.15d. The results illustrated in Figure 3.15 do not qualitatively change
at other parameter combinations.
The two positions where the Nusselt number becomes infinite and zero and the
conditions for these phenomena to occur are of big importance, both from the practical
and from the academical point of view and how they depend on the boundary conditions
will by analysed in the next section.
3.3.1 Positions of the vertical asymptote and the adiabatic point
The local Nusselt number NuW1or NuW2can become infinite positive or negative, i.e.
a vertical asymptote can occur at some axial position where the mean fluid temperature
in the channel cross-section, Equation (3.38), becomes equal to the temperature of the
lower plate, ΘW1= 0, or of the upper plate, ΘW2. Independently of this, the temperature
slope, Equation (3.39), can become zero at one of the plates, at another axial position.
Consequently, the heat transfer coefficient αW1(Z∗) or αW2(Z∗) and the Nusselt number
NuW1or NuW2can become zero, which characterises the adiabatic point.
The positions that determine the distances between the channel inlet and the cross-
sections where the vertical asymptote of the Nusselt number and the adiabatic point
occur are denoted by ZA∗and Z0∗, respectively. The boundary conditions (3.7) to (3.9)
determine if these appearances are possible at all and where, at the lower plate or at the
upper plate, and at which distance from the channel inlet. In the dimensionless form
this situation is affected only by the dimensionless temperature of the upper plate, ΘW2,
since the dimensionless temperatures of the lower plate ΘW1= 0 and of the fluid at the
inlet ΘIN = 1 are permanently constant.
The vertical asymptote of the Nusselt number and the adiabatic point establish only
at the plate, the temperature of which is closer to the fluid inlet temperature ϑIN and
only if the fluid inlet temperature ϑIN is smaller or greater than the plate temperatures,
ϑW1and ϑW2. If ϑW1< ϑIN < ϑW2or ϑW1> ϑIN > ϑW2the vertical asymptote of
the Nusselt number and the adiabatic point do not occur at any plate. In terms of the
dimensionless temperature this means that there are three characteristical intervals of
ΘW2, namely:
•ΘW2<0; in this case the vertical asymptote and the adiabatic point establish at
the lower plate, X=−1, and for this reason the notation ΘW2(X=−1) is used in
this analysis,
•ΘW2= 0; this corresponds to the classical symmetrical Graetz problem, where the
Nusselt number curves change continuously and neither the vertical asymptote nor
the adiabatic point occur,
•0<ΘW2<1; the vertical asymptote and the adiabatic point establish at the upper
plate, X= +1, and this dimensionless temperature is noted here by ΘW2(X=+1),
35
3 Heat Transfer in a Parallel Plate Channel
•ΘW2>1; in this interval neither the vertical asymptote nor the adiabatic point
are possible at any plate.
As stated above, the vertical asymptote of the Nusselt number establishes at the
lower plate, X=−1, if
¯
Θ(ZA∗) = ΘW1= 0,(3.43)
or with (3.38)
ΘW2(X=−1) =3
2
∞
X
n=1
An
λ2
n
(ϕ0
n(+1) −ϕ0
n(−1)) e−λ2
nZA∗,(3.44)
which can be solved for a real ZA∗only if ΘW2<0.
For the upper plate, X= +1, similar condition is valid
¯
Θ(ZA∗) = ΘW2,(3.45)
or with (3.38)
ΘW2(X=+1) =−3
2
∞
X
n=1
An
λ2
n
(ϕ0
n(+1) −ϕ0
n(−1)) e−λ2
nZA∗.(3.46)
The Equation (3.46) can be solved for a real ZA∗if 0 <ΘW2<1, as stated above.
The conditions for the adiabatic point to occur are
∂Θ(X, Z∗)
∂X X=−1
,∂Θ(X, Z∗)
∂X X=+1
,(3.47)
or using Equation (3.39)
ΘW2+ 2
∞
X
n=1
Anϕ0
n(X)e−λ2
nZ0∗= 0,(3.48)
which can be solved for a real Z0∗if ΘW2<0, which corresponds to the lower plate,
X=−1, or if 0 <ΘW2<1, which corresponds to the upper plate, X= +1.
For ΘW2>1 there are no real solutions, neither for ZA∗nor for Z0∗.
Figure 3.16 shows the positions of the vertical Nusselt number asymptote and of the
adiabatic point (the reversal of the heat flux) at both plates. Plotted is the dimensionless
temperature ΘW2at zero heat transfer potential, ¯
Θ(Z∗) = ΘW1and ¯
Θ(Z∗) = ΘW2, and
at the zero temperature slope, (∂Θ/∂X)W1= 0 and (∂Θ/∂X)W2= 0, in dependence
of the coordinate Z∗. Since ΘW1= 0, the thermal symmetry is reached at ΘW2→0.
This means that the two positions, ZA∗and Z0∗, move downstream as ΘW2tends to
36
3.3 The Nusselt Number
Figure 3.16: Transitions and asymptotical behaviour of the Nusselt numbers for the
region ΘW2<1.
zero and the thermal asymmetry diminishes. Disregarding the limiting cases Z∗→0
and Z∗→ ∞, the inequality ZA∗< Z0∗is satisfied always.
From the physical reasoning it is irrelevant for the phenomena to occur, which of
the plates has a temperature closer to the fluid inlet temperature. Consequently, the
curves in Figure 3.16 are expected to be transferable into each other, and indeed, they
are reciprocal. Using the properties of the eigenfunctions, equations (3.27) and (3.28),
the following important relationship can be deduced
(1 −ΘW2(X=−1))(1 −ΘW2(X=+1)) = 1.(3.49)
This expression describes the dependency of the dimensionless temperatures of the upper
plate, ΘW2(X=−1) and ΘW2(X=+1), for the vertical Nusselt number asymptotes to occur
at the same axial position ZA∗, at the lower and at the upper plate, respectively. The
dependency (3.49) was first derived by Mitrovi´c, Maleti´c and Baˇcli´c [32].
37
3 Heat Transfer in a Parallel Plate Channel
38
4 Fluid Flow and Heat Transfer in
Thermoplates
4.1 Problem formulation
The objective of this numerical study is to obtain sets of the geometrical parameters that,
in interaction with process parameters, should pave the way for optimal heat transfer of
the inside fluid, which is assumed to pass the thermoplate as a single phase (coolant in
Figure 2.1).
The overall size of the thermoplate can vary in wide ranges and the connecting
tubes can be placed and shaped in a different way. The welding spots can be arranged
in a staggered or in a parallel manner and their shape and size can also vary. To
investigate the influence of all these parameters on the hydrodynamics and on the heat
transfer, a huge number of numerical experiments would be necessary. For each of these
experiments, a huge number of the numerical cells would be needed to simulate the fluid
flow and heat transfer in whole thermoplate. All this requires again a large computer
capacity, a very long computational time, especially in case of transient calculations,
and increases total costs of the numerical experiments.
Neglecting the influence of the connecting tubes and of the seam-welded thermoplate
edges, a few symmetry planes can be introduced, which reduces the size of the numerical
domain, i.e. the number of the numerical cells and the number of the geometrical
parameters that are to be varied.
The geometry and dimensions of all simulated three-dimensional models are shown
in Figure 4.1. The calculation domain is a strip of a thermoplate channel. The semi-
circles represent the welding spots. The strip is confined in spanwise direction by the
planes z= 0 and z=sTand in streamwise (x) direction by x= 0 and x=L.
Depending on the thermoplate type, Figure 2.2, and on the welding spot pattern
(staggered or parallel), the heat transfer surface area of the thermoplate, that is consid-
ered to be a three-dimensional wavy surface, can be approximated by employing different
mathematical functions. If the welding spots are arragned in a staggered manner, Fig-
ure 4.1a, this surface can be represented by trigonometrical functions in the following
manner:
y=±δ
41 + cos 2π
sL
xcos π
sT
z,(4.1)
39
4 Fluid Flow and Heat Transfer in Thermoplates
Figure 4.1: Geometry and dimensions of the three-dimensional simulated domain; a)
welding spots arrangend in a staggered manner, Equation 4.1, b) parallel welding spots,
Equation (4.2) with C=sL/2, c) parallel welding spots, Equation (4.2) with C=sL/4,
d) parallel welding spots, Equation (4.2) with C= 0, and e) thermoplate surface defined
using the exponential function (4.4).
40
4.1 Problem formulation
where the sign + is used for the upper, and −for the lower surface, respectively, and the
geometrical parameters δ,sLand sTcan be taken from Figure 4.1. For welding spots
arranged in parallel, Figure 4.1b, c and d, the following expression is applied
y=±δ
82−cos 2π
sL
(x−C)−cos 2π
sT
z.(4.2)
The cross-section area of the model with staggered welding spots, Figure 4.1a, is nearly
independent on x(in this case, it changes slightly only in the welding spots regions),
while at models with parallel welding spots it changes with x-coordinate noticeable. The
size of the inlet surface can be adjusted by constant C. In case of C= 0, Figure 4.1d,
the inlet surface is minimal and if C=sL/2 it is maximal, Figure 4.1b. For C=sL/4
or C= 3sL/4, the inlet surface of the model with parallel welding spots is in between
of those with C= 0 and C=sL/2 and is nearly same as the one of the model with
staggered arangement of the welding spots, Figure 4.1c.
The computational domain is determinated in each case by the surface y=f(x, z)
according to the Equation (4.1) or (4.2), the planes x= 0 and x=L, the symmetry
planes y= 0, z= 0 and z=sTor z=sT/2 for models with staggered and parallel
welding spots, respectively, and the cylindrical surfaces
(x−x0i)2+ (z−z0i)2=R2,(4.3)
where x0iand z0iare the centres of the welding spots obtained as the roots of the
equation y=f(x, z) = 0 representing the heat transfer wavy surface according to the
Equation (4.1) or (4.2).
However, the heat transfer surface area can be described also using other types of
mathematical functions, i.g.:
y=±δ
2 1−
N
X
i=1
exp −D(x−x0i)2+ (z−z0i)2!,(4.4)
where Nrepresents the total number of welding spots, Dis a positive constant, x0iand
z0iare the coordinates of the welding spot centres, which can be arranged freely in this
case, in a staggered, parallel or any other manner, and have to be defined previously.
Choosing an appropriate value Dthe curvature of the heat transfer surface can be
changed. Due to the asymptotic features of the exponential function, the hight of the
channel defined using Equation (4.4) is always slightly smaller than δ, which can be
neglected in the engineering practice, for sufficient high value of constant D. An example
of the thermoplate surface defined using the exponential function (4.4) is illustrated
schematically in Figure 4.1e.
In this context, it is to be noted that the mathematical contour of the thermoplate
according to Equations (4.1), (4.2) and (4.4) does not ideally represent the reality, but
the discrepancy as regards the thermohydraulic behaviour of the system is expected to
be of a marginal importance.
41
4 Fluid Flow and Heat Transfer in Thermoplates
The objective of this numerical analysis is to numerically find out the influence of the
thermoplate geometry on the heat transfer and pressure drop. Based on the thermoplate
that was used in the laboratory for the experimental research [38], whose geometry was
approximated by the Equation (4.1), only one of the geometrical parameter, sL,sT, or
δ, Figure 4.1a, and the Reynolds number have been varied firstly, whereas the other two
parameters have been kept constant, as it is specified in the Table 4.1.
Table 4.1: Parameter combinations for the numerical experiments.
sL[mm] sT[mm] δ[mm] L[mm] Re
28 ÷112 36 3.4 500 50 ÷4000
42 20 ÷50 3.4 500 50 ÷4000
42 36 2.5÷8.4 500 50 ÷4000
Although there are 21 different Models at 11 different Reynolds Numbers according
to the Table 4.1, which results in 231 numerical experiments, this parameter combination
does not suffice for a Nusselt number and pressure drop correlations to be stated, which
belongs to the objectives of this thesis as well.
For that reason, the parameter combination from the Table 4.1 has been supple-
mented to some extent, Table 4.2.
Table 4.2: Extended parameter combinations for the numerical experiments.
sL[mm] sT[mm] δ[mm] L[mm] Re
28 ÷112 20 ÷50 2.5÷8.4 500 50 ÷3000
Additionaly, a comparison between different thermoplate geometries and welding
spots arrangements, Figure 4.1 had to be made, Table 4.3.
Table 4.3: Different model geometries.
Model sL[mm] sT[mm] δ[mm] L[mm] Re
Equation (4.1), Figure 4.1a 40 35 4 300 50 ÷2000
Equation (4.2), Figure 4.1b 40 35 4 300 50 ÷2000
Equation (4.2), Figure 4.1c 40 35 4 300 50 ÷2000
Equation (4.2), Figure 4.1d 40 35 4 300 50 ÷2000
Equation (4.4), Figure 4.1e 40 35 4 300 50 ÷2000
42
4.2 The physical model
Further, the influence of time fluctuations of the fluid velocity at the inlet on the
heat transfer in thermoplate is to be investigated, as it is stated in Table 4.4.
Table 4.4: Inlet velocity for time dependent fluid flows.
uIN [m/s] sL[mm] sT[mm] δ[mm] L[mm] Re
0.1 40 35 4 300 450
0.1 (1 + 0.2 (sin(100t))) 40 35 4 300 450
0.1+0.04 P∞
n=0 H(t−nπ
100 ) 40 35 4 300 450
To find out the influence of the thermally asymmetric boundary conditions on the
heat transfer in thermoplates and to compare it to the heat transfer in a parallel plate
channel under the same boundary conditions, Chapter 3, the models defined in Table
4.5 are to be simulated.
Table 4.5: Thermally asymmetric boundary conditions.
TW1[K] TW2[K] sL[mm] sT[mm] δ[mm] L[mm] Re
333 293 40 35 4 300 50 ÷2000
333 303 40 35 4 300 50 ÷2000
333 306.33 40 35 4 300 50 ÷2000
333 317 40 35 4 300 50 ÷2000
333 319.67 40 35 4 300 50 ÷2000
333 325 40 35 4 300 50 ÷2000
333 343 40 35 4 300 50 ÷2000
333 353 40 35 4 300 50 ÷2000
4.2 The physical model
The fluid flow was considered to be laminar, single phase, incompressible and three-
dimensional. The velocity, pressure and temperature fields are governed by the well
43
4 Fluid Flow and Heat Transfer in Thermoplates
known equations of continuity, momentum and energy
∂ρ
∂t +∇·(ρw)= 0,(4.5)
ρDw
Dt =−∇p−∇·τ+ρf,(4.6)
ρcv
DT
Dt =−∇· ˙q−T∂p
∂T ρ
(∇·w)−(τ:∇w),(4.7)
where ρdenotes fluid density, cvthe specific heat capacity, wthe velocity vector, ˙qthe
heat flux vector, τthe stress tensor, Tthe temperature, pthe pressure and fthe body
force. By the derivation of the energy equation 4.7, the equation
i=i(v, T),(4.8)
was used, where irepresents the specific internal energy and vthe specific volume. In
case of compressible flows, an additional equation, the thermal equation of state
p=p(v, T),(4.9)
is needed to close the system of the equations.
In order to use the Equations (4.5) to (4.7), the stress ternsor must be expressed in
terms of velocity gradients and the heat flux vector in terms of temperature gradients [4].
For Newtonian fluids the components of the stress tensor τin the Cartesian coordinate
system are:
τxx =−2µ∂u
∂x +2
3µ(∇·w),(4.10)
τyy =−2µ∂v
∂y +2
3µ(∇·w),(4.11)
τzz =−2µ∂w
∂z +2
3µ(∇·w),(4.12)
τxy =τyx =−µ∂u
∂y +∂v
∂x,(4.13)
τyz =τzy =−µ∂v
∂z +∂w
∂y ,(4.14)
τzx =τxz =−µ∂w
∂x +∂u
∂z .(4.15)
Here are u,vand wcomponents of the velocity vector win x,yand zdirection respec-
tively and µthe dynamic viscosity.
44
4.2 The physical model
Using the Fourier’s low of heat conduction for constant k
˙q=−k∇T, (4.16)
where krepresents the thermal conducivity, and properties of the stress tensor (4.10 to
4.15) the governing equations (4.5) to (4.7) become:
∂ρ
∂t +∂(ρu)
∂x +∂(ρv)
∂y +∂(ρw)
∂z = 0,(4.17)
ρ∂u
∂t +u∂u
∂x +v∂u
∂y +w∂u
∂z =−∂p
∂x +µ∂2u
∂x2+∂2u
∂y2+∂2u
∂z2+ρfx,(4.18)
ρ∂v
∂t +u∂v
∂x +v∂v
∂y +w∂v
∂z =−∂p
∂y +µ∂2v
∂x2+∂2v
∂y2+∂2v
∂z2+ρfy,(4.19)
ρ∂w
∂t +u∂w
∂x +v∂w
∂y +w∂w
∂z =−∂p
∂z +µ∂2w
∂x2+∂2w
∂y2+∂2w
∂z2+ρfz,(4.20)
ρcv∂T
∂t +u∂T
∂x +v∂T
∂y +w∂T
∂z =k∂2T
∂x2+∂2T
∂y2+∂2T
∂z2+µΦ (4.21)
−T∂p
∂T ρ∂u
∂x +∂v
∂y +∂w
∂z ,(4.22)
where fx,fyand fzrepresent the components of the body force vector fand the dissi-
pation function Φ is defined as follows
Φ =2 ∂u
∂x2
+∂v
∂y2
+∂w
∂z 2!
+ ∂u
∂y +∂v
∂x2
+∂u
∂z +∂w
∂x 2
+∂v
∂z +∂w
∂y 2!
−2
3∂u
∂x +∂v
∂y +∂w
∂z 2
.
(4.23)
To complete the set of the equations, the model equations (4.17) to (4.21) are
45
4 Fluid Flow and Heat Transfer in Thermoplates
subject to the following boundary coditions:
x= 0 : w={uIN ,0,0}, T =TIN ;x=L:∂w
∂x = 0,∂T
∂x = 0,(4.24)
y= 0 : ∂w
∂y = 0,∂T
∂y = 0; y=f(x, z) : w= 0, T =TW,(4.25)
z= 0 : ∂w
∂z = 0,∂T
∂z = 0; z=sT:∂w
∂z = 0,∂T
∂z = 0.(4.26)
For the models with parallel welding spots, Figure 4.1, there is an additional sym-
metry plane at z=sT/2 , which is used to economise the calculations. In this case,
instead of the conditions (4.26), following conditions are used:
z= 0 : ∂w
∂z = 0,∂T
∂z = 0; z=sT
2:∂w
∂z = 0,∂T
∂z = 0.(4.27)
The small cylindrical surfaces representing the welding spots according to the Equa-
tion (4.3) are kept adiabatic
(x−x0i)2+ (z−z0i)2=R2:w= 0,∂T
∂n = 0,(4.28)
where nrepresents a coordinate normal to the cylindrical surfaces.
By the conditions (4.24), the momentum and heat transfer in the streamwise direc-
tion at x=Lare neglected. Clearly, this simplification is not met in reality, but the
errors thus introduced affect the corresponding fields only in the immediate vicinity of
the channel outlet. The pressure boundary at x=Lwas not used because of its negative
inpact on the calculating time.
In the practice the walls of the thermoplate wavy channel are usually at different
temperatures. The thermal asymmetry in parallel plate channel, which is used as a
reference for comparison, was analysed in the Chapter 3. If the lower and the upper wall
of the wavy cannel are kept at different temperatures, the symmetry plane at y= 0 and
the boundary conditions (4.25) are invalid. Instead, the following conditions would be
required:
y=−f(x, z) : w= 0, T =TW1;y=f(x, z) : w= 0, T =TW2,(4.29)
where TW1and TW2are the temperatures of the lower and of the upper wall, respectively.
Note that the lower and the upper thermoplate surface do not touch each other, because
the adiabatic cylindrical surface (4.3) is placed in between. That allows the lower and
the upper wall to have different temperatures, without any singularity problems.
The governing equations are solved using a commercial software StarCD by CD-
adapco. A brief overview of the numerical methods used by this programme is given in
Appendix A.
46
4.3 Process quantities
An example of the thermoplate model definition in StarCD using input files is
illustrated in Appendix B. In this program, the thermoplate surface is approximated
using Equation 4.1, without the symmetry plane at y= 0 and the wall temperatures cann
be different and specified by user, as well as other geometrical and process parameters.
Not only a thermoplate strip, but also the whole thermoplate with connecting tubes cann
be modelled by choosing appropriate values of the parameters. For other thermoplate
geometries or welding pattern, similar programmes have been used.
4.3 Process quantities
4.3.1 Symmetric heated channels
Once the governing equations are solved and the flow, the pressure and the temperature
fields are obtained, which makes possible the definition of other very important process
quantities.
The fluid bulk temperature (also called the mixing cup temperature), Tb(x), and
the local wall heat flux, ˙qW(x), averaged in spanwise (z) direction, can be calculated by
Tb(x) = RSuTdS
RSudS ,(4.30)
˙qW(x) = 1
ASW ZASW
(˙q·n)dASW ,(4.31)
where S=S(x) represents the cross-sectional flow area of the channel, ASW =ASW (x)
a narrow strip of the thermoplate surface y=f(x, z) bounded by the planes x−(ε/2)
and x+ (ε/2), εbeing the width of the numerical cell in the x-direction.
The heat flow rate, ˙
Q(x), and the average wall heat flux, ˙qmW (x), are obtained in
a similar way
˙
QW(x) = ZAW
(˙q·n)dAW,˙qmW (x) = ˙
QW(x)
AW
=1
AWZAW
(˙q·n)dAW,(4.32)
where AW=AW(x) denotes the heat transfer area according to the thermoplate surface
y=f(x, z), from the very beginning of the thermoplate, x= 0, to the current position
x+ε/2.
The local and the average heat transfer coefficients, hW(x) and hmW (x), are defined
by
hW(x) = ˙qW(x)
TW−Tb(x), hmW (x) = 1
xZx
0
hW(ξ)dξ, (4.33)
47
4 Fluid Flow and Heat Transfer in Thermoplates
with Tb(x) and ˙qW(x) according to Equations (4.30) and (4.31); TWrepresents the wall
temperature, that does not change with any coordinate in this analysis.
The average heat transfer coefficient, hmW (x), defined according to the Equation
(4.33) can be calculated if the distribution of the local heat transfer coefficient hW(x) is
known, which is the case in this numerical analysis. However, if these numerical results
are to be compared with the experimental ones, some difficulties could be experienced,
because, the average heat flux and the fluid temperature at the inlet and at the outlet
are measured in experiments rather than the temperature distribution of the flowing
fluid through the whole thermoplate and the heat flux distribution on its surface. This
would be either impossible or too expensive. For that reason, another definition of the
heat transfer coefficient is proposed:
hmW =˙qmW
∆Tm
,(4.34)
where hmW and ˙qmW denote the average heat transfer coefficient and the average heat flux
for the whole thermoplate, respectively, and ∆Tmrepresents the average temperature
difference between the thermoplate surface and the fluid. This temperature difference
can be calculated as follows
∆Tm=1
LZL
0
(TW−Tb(x)) dx =TW−Tbm,(4.35)
where Lrepresents the thermoplate length, Tb(x) the fluid temperature according to the
Equation (4.30) and Tbm the mean fluid temperature in whole thermoplate. However,
this average temperature difference (4.35) suffers from a lack of being inapplicable for
comparison with exmerimental results, for the same reasons as it is the case with the
average heat transfer coefficient, hmW (x), discussed above. An acceptable average tem-
perature difference can be deduced from the heat transfer balance, writing it first in
infinitesimal form, [38]:
d˙
QW= ˙mcpdT =hmW (TW−Tb(x))dA, (4.36)
and than in integral form
˙
QW= ˙mcp(TOUT −TIN ) = hmW ∆TmA, (4.37)
where ˙mand Arepresent the mass flow rate of the fluid through the thermoplate and
thermoplate surface, respectively. Integrating the Equation (4.36) from the inlet to the
outlet gives:
hmW =˙mcp
Aln TW−TIN
TW−TOUT
,(4.38)
where Tb(0) = TIN and Tb(L) = TOUT . Combining the Equations (4.38) and (4.37)
enables calculating of the average temperature difference, also called the logaritmic mean
temperature difference (LMTD), as follows:
∆Tm=TOUT −TIN
ln TW−TIN
TW−TOUT
.(4.39)
48
4.3 Process quantities
This logarithmic mean temperature difference depends only on the fluid themperatures
at the inlet and outlet, and the wall temperature. However, the temperature difference,
∆Tm, calculated according to Equation (4.39) does not take into account any changes
of fluid properties and heat transfer coefficient along the channel. For that reason a
generalised mean temperature difference, that would take into these changes, is often
recommended, e.g. by Utamura et al. [67]. To employ this generalisation method one
would need some other experimental data of the heat exchanger, and the generalised
mean temperature difference would have to be calculated numerically. This is, however,
beyond the scope of this thesis.
The local and the average Nusselt numbers, NuW(x) and NumW (x), are obtained
from
NuW(x) = hW(x)dh
k, NumW (x) = hmW (x)dh
k,(4.40)
with dhas the hydraulic diameter
dh=4S
U= 2d, (4.41)
Sbeing the cross-sectional flow area, Uits circumference, and dthe distance between
the parallel plates channel having the same volume Vas the actual thermoplate
d=V
sTL.(4.42)
Note that the hydraulic diameter dhdefined in this way does not change along the
channel length.
A further quantity involved in the presentation of the results is the Reynolds num-
ber,
Re =¯udh
ν,(4.43)
where ¯uis the average axial velocity component, and νthe kinematic viscosity.
In order to obtain a thermo-fluid characteristic of the thermoplate, the local pres-
sure, p(x), averaged in the plane y0zand the pumping power, P, are calculated as
follows
p(x) = 1
SZS
pdS, (4.44)
P(x) = (p(0) −p(x)) ˙
V , (4.45)
where p(0) = pIN and p(x) are the averaged inlet pressure and the average pressure at
distance xfrom the inlet, respectively; ˙
Vdenotes the volume flow rate (for x=Lwe
have pm(L) = pOUT ,Lbeing the channel length).
49
4 Fluid Flow and Heat Transfer in Thermoplates
4.3.2 Asymmetric heated channels
In the case of asymmetrical heating the symmetry plane at y= 0 and the boundary
conditions (4.25) are not applicable. Instead, the conditions (4.29) have to be employed.
Consequently, the quantities that are calculated at the upper and at the lower ther-
moplate surface, (4.31 to 4.40), will not have the same values, as in the case of the
symmetrical heating, and have to be defined separately.
The local heat flux at the lower wall, y=−f(x, z), and at the upper wall, y=
f(x, z), is defined:
˙qW1(x) = 1
ASW 1ZASW 1
(˙q·n)dASW 1,(4.46)
˙qW2(x) = 1
ASW 2ZASW 2
(˙q·n)dASW 2,(4.47)
where ASW 1and ASW 2are, similar to ASW , narrow strips of the thermoplate lower
surface, y=−f(x, z), and of the upper surface, y=f(x, z), respectively, bounded by
the planes x−(ε/2) and x+ (ε/2).
Similar, the average wall heat flux and the heat flow rate at the lower and at the
upper wall are:
˙qmW 1(x) = 1
AW1ZAW1
(˙q·n)dAW1,˙
QW1(x) = ZAW1
(˙q·n)dAW1,(4.48)
˙qmW 2(x) = 1
AW2ZAW2
(˙q·n)dAW2,˙
QW2(x) = ZAW2
(˙q·n)dAW2,(4.49)
where AW1and AW2are the lower and the upper wall surface, respectively, in range
between x= 0 and x+ε/2.
The local and the mean heat transfer coefficients at the lower and at the upper
thermoplate surface become
hW1(x) = ˙qW1(x)
TW1−Tb(x), hmW1(x) = 1
xZx
0
hW1(ξ)dξ, (4.50)
hW2(x) = ˙qW2(x)
TW2−Tb(x), hmW2(x) = 1
xZx
0
hW2(ξ)dξ. (4.51)
Since there are no experimental data available on asymmetric heated thermoplates,
no other definitions of the average heat transfer coefficient are considered.
50
4.3 Process quantities
Finaly, the Nusselt number at the lower and at theupper thermoplate surface are
defined as follows:
NuW1(x) = hW1(x)dh
k, NumW 1(x) = hmW 1(x)dh
k,(4.52)
NuW2(x) = hW2(x)dh
k, NumW 2(x) = hmW 2(x)dh
k.(4.53)
All the process quantities defined in this section are calculated in post-processing
part of the simulations with StarCD using a Fortran subroutine POSDAT.f, Appendix
C. Similar program have been used for transient calculations.
The dependency of the heat transfer coefficient on the size of the numerical mesh
is diskussed in Appendix D.
51
4 Fluid Flow and Heat Transfer in Thermoplates
52
5 Results
Because of the complex geometry of the computational domain, all the vector quan-
tities obtained from the governing equations are three-dimensional. This poses some
difficulties on the presentation of the results. For this reason, the interesting quantities
will be visualised in the following as general images and quantitatively illustrated in
characteristic planes or on surfaces, for one specific channel, i.e. for one combination
of the geometrical parameters δ= 3.4 mm, sL= 42 mm, sT= 36 mm, R= 5 mm and
L= 500 mm, Figure 4.1a. For other parameter combinations the similar flow behaviour
was observed.
5.1 Symmetric heated channels, steady state flows
5.1.1 The velocity field
The development of the velocity field in the symmetry plane y= 0 at different Reynolds
numbers is shown in Figure 5.1. At the smallest Reynolds number adopted in the
numerical experiments, Re = 50, the velocity field is largely smooth with relatively
small separation zones that are established behind the welding spots. Nevertheless, the
velocity fluctuations both in the transversal, z, and axial, x, direction are observed.
Increasing the Reynolds number the separation zones grow in the region between
the welding spots that are placed one behind another, finally reaching the next welding
spot. This is already the case at Re = 500, Figure 5.2. With the strictly laminar flow,
the fluid is recirculating and detained in these zones. The thermoplate surface area
corresponding to these recirculation zones is expected to be less effective for the heat
transfer.
At the Reynolds number Re = 4000 the fluid separation establishes between the
welding spots with reattachments at the contour of the neighbouring, downstream spots,
and a comparatively large portion of the channel is occupied by recirculation zones. In
the middle of the channel, there is a meandering fluid core that is bounded by the
recirculation zones and only touches the welding spots. Here is the velocity almost four
times greater than the inlet fluid velocity.
Similar flow structures have been reported by Wang and Vanka [70], on the basis of
the numerical experiments, and by Nishimura et al. [42] as well as by Rush et al. [52]
in laboratory experiments, all using two-dimensional channels.
53
5 Results
Figure 5.1: Velocity field in the symmetry plane y= 0 at different Reynolds numbers.
54
5.1 Symmetric heated channels, steady state flows
Figure 5.2: Velocity profiles in the welding spots cross-sections at Re = 4000.
55
5 Results
Figure 5.3: Pressure field in the symmetry plane y= 0 at Re = 4000 and its linkage
with the velocity field.
56
5.1 Symmetric heated channels, steady state flows
As an example, Figure 5.2 provides an idea about the velocity profile in y0zplanes
at different distance from the inlet (at the cross-sections where the welding spots are
placed) for Re = 4000. A very strong velocity gradient is observed in the immediate
vicinity of each welding spot, where the velocity jumps from zero on the welding spot
rims, to a value of about 1.5 m/s. Somewhere in the middle of the channel, in the kernel
flow region, the velocity reaches its maximum, and reverses on the opposite side to the
welding spot in respect to zaxis, where the fluid flow is dominated by an adverse pressure
gradient, Figure 5.3. Interestingly, in the recirculation zones, the velocity maximum does
not establish in the symmetry plane x0z, but is shifted towards the wall. The like is
observed in the kernel flow near the boundary of the recirculation zone, Figure 5.3. This
means that the velocity profiles in this case do not have a parabolic-like form with one
maximum, but they have two inflection points with two local maxima and one local
minimum, the latter being placed in the symmetry plane y= 0. Similar behaviour has
also been observed by Mahmud et al. [29], Figure 2.9.
The distribution of the total pressure in the symmetry plane, y= 0, at Re = 4000,
as well as its linkage to the velocity vector direction in the recirculation zones, is shown
in Figure 5.3. In the middle of the channel, the pressure decreases with increasing
distance from the channel inlet. However, in the recirculation zones the pressure is not
monotonically decreasing function of xand it even increases with xin the regions, where
the thermoplate walls converge, which is the case near the symmetry planes at z= 0
and z=sTnext to each welding spot in upstream direction. This piecewise adverse
pressure gradients cause fluid flow reversals, resulting in recirculation zones.
From Figure 5.2 one might infer that in the same cross-section turbulent flow may
coexist with the laminar flow (apart from the viscous boundary layer in the near-wall
regions). The extension of one flow mode occurs at the expenses of the other one, and
turbulence will occupy the larger portion of the channel cross-section, the larger the
Reynolds number. The curvature of the boundary is expected to support the generation
of G¨ortler’s vortices and hairpin vortex structures [16]. As a consequence, this would
result in a stronger pressure drop and an improvement of heat transfer in comparison
with those obtained from the purely laminar model. In this context it is to be noted
that a possible laminar-turbulent flow transitions are disregarded in these simulations.
5.1.2 The temperature field and the heat flux at the thermoplate
surface
Some insights into the development of the temperature field in the x0zsymmetry plane
of the channel at different Reynolds numbers are given in Figure 5.4. The fluid enters
the thermoplate at the inlet temperature TIN = 293 K and is heated by the thermoplate
walls, which are kept isothermal at TW= 333 K.
At the Reynolds number Re = 50 the fluid temperature increases very rapidly and
the thermal boundary layer penetrates the fluid bulk relatively fast, reaching the x0z
57
5 Results
Figure 5.4: Temperature field in the symmetry plane y= 0 at different Reynolds
numbers.
58
5.1 Symmetric heated channels, steady state flows
Figure 5.5: Wall heat flux at different Reynolds numbers.
59
5 Results
plane downstream the channel inlet at a distance smaller than sL/2. Near the welding
spots and in the recirculation zones, the fluid approaches the temperature of the walls
much faster than in the middle channel region.
Increasing the Reynolds number up to Re = 500, the temperature field changes
significantly. In the central channel region, z≈sT/2, the fluid remains thermally
unaffected by the heating plates, and the thermally developing region spreads up to
x≈2·sL. The recirculation zones are clearly recognised and the fluid reaches the
temperature of the walls primarily behind the welding spots in the downstream direction
and in the centres of the eddies, which is expected to reduce the heat transfer in these
regions.
Further increasing of the Reynolds number up to Re = 4000 does not change the
temperature field qualitatively. The thermally developing region occupies now a channel
length of about 3·sL. The fluid velocity is high, so that the fluid, even in the recirculation
zones is thermally active.
For better understanding of the heat transfer mechanism in thermoplates, the heat
flux should be analysed in detail. Its distribution as a vector intensity, |˙q|, at different
Reynolds numbers is illustrated in Figure 5.5. At all Reynolds numbers, the heat flux
is seen to decrease in the streamwise, x, direction, becoming smaller in the recirculation
zones, than in the central part of the strip. The heat flux fluctuates significantly not only
in spanwise, but also in streamwise direction. The local heat flux maxima are situated
upstream near the welding spots, since its local minima are placed behind the welding
spots in the flow direction.
Figure 5.6: Local, spanwise averaged, heat flux, ˙qW= ˙qW(x), and the mean heat flux,
˙qmW = ˙qmW (x), according to Equations (4.31) and (4.32) at different Reynolds numbers.
60
5.1 Symmetric heated channels, steady state flows
Figure 5.7: Local, spanwise averaged, heat transfer coefficient, hW=hW(x), and the
mean heat transfer coefficient, hmW =hmW (x), according to the Equation (4.33) at
different Reynolds numbers.
Figure 5.8: Heat transfer coefficients, hW=hW(x) and hmW =hmW (x), according to
the Equation (4.33), with fully developed region in a 500mm long channel at Re = 50.
61
5 Results
Further insights into the heat transfer fluctuations is provided by the local heat
flux, ˙qW= ˙qW(x), averaged over the width, sT, of the strip according to Equation
(4.31), the distribution of which is demonstrated in Figure 5.6 at different Reynolds
numbers. When analysing these heat flux distributions, the heat flux fields shown in
Figure 5.5 should be kept in mind. The quantity ˙qW(x) behaves damped oscillatory
along the channel, decaying at the same time both in magnitude and amplitude. The
positions of the maxima of ˙qW(x) nearly coincide (they are shifted in the upstream
direction) with the axial positions of the welding spots, which was already observed in
Figure 5.5. Averaging the heat flux ˙qWstreamwise as stated in Equation (4.32), thus
getting the quantity ˙qmW = ˙qmW (x), results in a similar distribution with oscillations
that are stronger damped.
The distributions of the local, spanwise over sTaveraged, heat transfer coefficient,
hW(x), along the channel, and of the heat transfer coefficient, hmW (x), averaged over the
current length, x, according to Equation (4.33), are illustrated in Figure 5.7. The curves
representing hW(x) and hmW (x) are quite similar to those for the heat flux in Figure 5.6.
The decrease of the heat transfer coefficients in the streamwise direction indicates that
the channel length of 150mm is too short for the thermal developing region to become
completed, even at the Reynolds number Re = 50, Figure 5.8. This heat transfer figure
changes when the channel is sufficiently long and the fully developed heat transfer region
is approached. Figure 5.7 illustrates this behaviour for the parameters noted in Figure
5.5.
The distributions of the local and average Nusselt number, NuW(x) and NumW (x)
defined in Equation 4.40 along the channel length are similar to the corresponding heat
transfer coefficients, hW(x) and hmW , respectively, and they are qualitatively comparable
with the corresponding ones reported by Garg and Maji [14], Russ and Beer [53], Saniei
and Dini [55] and Wang and Chen [69], Figures 2.11, 2.12 and 2.13, for two-dimensional
channels.
5.1.3 Influence of the geometrical parameters
In the numerical experiments, three geometrical parameters (sL,sTand δ) have been
varied, as well as the Reynolds number. The influence of varied geometrical parame-
ters on the process quantities will be shown in the following sections separately. For
this purpose, the geometrical parameter under consideration is varied together with the
Reynolds number, whereas the remaining two geometrical parameters are kept constant.
5.1.3.1 Welding spot spacing in the flow direction
Effects of the longitudinal welding spots pitch, sL, on the heat transfer coefficient for
500mm long channels at selected values of the Reynolds number and other geometrical
parameters are visualised in Figure 5.9. As may be taken from the diagram, a larger
62
5.1 Symmetric heated channels, steady state flows
Figure 5.9: Average heat transfer coefficient, hmW of 500 mm long channels at different
Re as a function of longitudinal welding spots pitch, sL;sL→ ∞: parallel plates.
Figure 5.10: Average Nusselt number, NumW of 500 mm long channels at different Re
as a function of longitudinal welding spots pitch, sL;sL/L → ∞: parallel plates.
63
5 Results
Figure 5.11: Average Nusselt number, NumW of 500 mm long channels as a function
of Re at different sL/L;sL/L → ∞: parallel plates.
Figure 5.12: Pumping power, P, Equation (4.45), needed for forced convection through
a 500 mm long channel as a function of the longitudinal welding spots pitch, sL, at
different Reynolds numbers, Re;sL→ ∞: parallel plates.
64
5.1 Symmetric heated channels, steady state flows
Reynolds number corresponds to a larger heat transfer coefficient. At the smallest
Reynolds number, Re = 50, the heat transfer is largely independent of sL, whereas at
other values of Re, certain dependency of hmW on sLis observed. At sL→ ∞,hmW
tends to the value for the channel formed by plane parallel plates. These values are
higher than the ones for the thermoplate at smaller Re. The situation reverses at larger
Re, where the heat transfer of the thermoplate becomes better. Certain maxima of hmW
are observed at sL≈80 mm at larger Reynolds numbers.
The Nusselt number, NumW depends on the longitudinal welding spots pitch, sL,
in the same way as the heat transfer coefficient hmW , since the hydraulic diameter and
the kinematic viscosity remain constant when the longitudinal welding spots pitch sLis
varied, Figure 5.10.
Figure 5.11 shows the effect of the Reynolds number, Re, on the heat transfer for
selected values of the geometrical parameters. The average Nusselt number NumW of the
thermoplate is larger than of the flat channel for Re ≥500. At Re ≤500 the secondary
flows are obviously too weak to outweigh the deterioration effect of recirculation zones
on the heat transfer, thus leading to lower NumW -values in comparison to those of the
flat channel.
The pumping power, P, calculated according to Equation (4.45), as a function of
the longitudinal welding spot pitch, sL, at different Reynolds numbers is visualised in
Figure 5.12. The pressure drop increases with increasing Reynolds number, resulting
in higher pumping power, needed for the forced convection through the channel. The
pumping power is a decreasing function of the welding spot spacing sL, at all Reynolds
numbers, having its minimum when sL→ ∞, which is the parallel plate channel.
5.1.3.2 Welding spot spacing in the cross direction
The distribution of the average heat transfer coefficient, hmW , as a function of transver-
sal welding spots pitch, sT, at different Reynolds numbers for 500 mm long channels is
shown in Figure 5.13. A larger Reynolds number corresponds to a larger heat trans-
fer coefficient, as it was the case when the longitudinal welding spots pitch was varied.
At all Reynolds numbers, Re, the heat transfer coefficient, hmW , behaves as a decreas-
ing function of the transversal welding spot pitch, sT. However, at a lower Reynolds
numbers, Re ≤500, this dependency is hardly recognisable.
Since the the hydraulic diameter and the thermal conductivity remain constant also
when the transversal welding spots pitch is varied, the dependency of the average Nusselt
number, NumW , on the spacing sTis very similar to the distribution of the average heat
transfer coefficient, hmW , Figure 5.14.
The same distribution of the average Nusselt number, NumW , as the one shown in
Figure 5.14 is visualised in Figure 5.15 as a function of the Reynolds number, Re,sL/L
being a parameter.
65
5 Results
Figure 5.13: Average heat transfer coefficient, hmW of 500 mm long channels at dif-
ferent Re as a function of transversal welding spots pitch, sT.
Figure 5.14: Average Nusselt number, NumW of 500 mm long channels at different Re
as a function of transversal welding spots pitch, sT.
66
5.1 Symmetric heated channels, steady state flows
Figure 5.15: Average Nusselt number, NumW of 500 mm long channels as a function
of Re at different sT/L.
Figure 5.16: Pumping power, P, Equation (4.45), needed for forced convection through
a 500 mm long channel as a function of the transversal welding spots pitch, sT, at
different Reynolds numbers, Re.
67
5 Results
Figure 5.16 visualises the dependency of the pumping power on the transversal
welding spots pitch, sT, calculated according to Equation (4.45). At lower Reynolds
numbers, Re ≤1000, the pumping power increases with increasing sT, but at larger
Reynolds numbers, it is first slightly decreasing, reaching a minimum, and than increas-
ing. In this context it is to be noted, that the pumping power is not only affected by
the pressure drop when the parameter sTis varied, as it was the case when the lon-
gitudinal welding spot spacing, sL, was varied. Increasing the parameter sT, the inlet
cross-section area of the channel is increasing as well, which indirectly increases the vol-
ume flow rate, when the Reynolds number is kept constant. In Figure 5.16 the influence
of both, pressure drop and volume flow rate is contained.
5.1.3.3 Maximal plate distance
The average heat transfer coefficient, hmW , as a function of the maximal plate distance,
δ, is illustrated in Figure 5.17. The heat transfer coefficient is a decreasing function of δ,
which is rooted in a corresponding reduction of the fluid velocity at the given Reynolds
number.
The corresponding average Nusselt number, NumW , as a function of the ratio δ/L
is shown in Figure 5.18. It is almost a linearly increasing function, but at the Reynolds
number Re = 50 the increase is hardly observable.
At first sight, the increase of NumW with increasing δ/L, Figure 5.18, could appear
contradictory, when comparing it with the distribution of the heat transfer coefficient,
Figure 5.17. The discrepancy between the heat transfer coefficient and Nusselt number
is associated with the increase of the hydraulic diameter, dh. The Nusselt number
calculated according to Equation (4.40) depends not only on the heat transfer coefficient,
but also on the hydraulic diameter, dh, which directly increases by rising the plate
distance δ.
The average Nusselt number, NumW , as a function of the Reynolds number, Re,
the ratio δ/L being a parameter, is illustrated in Figure 5.19.
Similar to the case, where the transversal welding spot pitch, sT, was varied, one
would expect that not only the pressure drop changes, when the the plate distance,
δ, is changed, but also the volume flow rate. However, the volume flow rate remains
unchanged when the plate distance, δ, i.e. the hydraulic diameter, dh, is varied, at a
constant Reynolds number. An increase of δ, i.e. dh, causes a decrease of the average
axial velocity, ¯u, which means, that the pumping power is in this case influenced only
by the pressure drop, as it was the case, when the parameter sLwas varied. Figure
5.20 provides an information about the dependency of the pumping power, P, on the
maximal plate distance, δ. It decreases rapidly with increasing plate distance, which is
connected with the velocity and pressure drop reduction.
68
5.1 Symmetric heated channels, steady state flows
Figure 5.17: Average heat transfer coefficient, hmW of 500 mm long channels at dif-
ferent Re as a function of maximal plate distance, δ.
Figure 5.18: Average Nusselt number, NumW of 500 mm long channels at different Re
as a function of maximal plate distance, δ.
69
5.1 Symmetric heated channels, steady state flows
5.1.3.4 Nusselt number correlation
The Nusselt number of the thermoplate is shown by the numerical experiments to depend
on the Reynolds number and the dimensionless geometrical parameters δ/L,sT/L and
sL/L.
Assuming a function of the form
NumW =C1sL
Ln1sT
Ln2δ
Ln3
Ren4Prn5,(5.1)
numerical results for the Nusselt number can be correlated in certain ranges of the
geometrical parameters. Since the Prandtl number was kept constant in this study
(Pr = 5.97), Equation 5.1 becomes
NumW =C2sL
Ln1sT
Ln2δ
Ln3
Ren4,(5.2)
where the constant C2takes into account the influence of the Prandtl number, C2=
C1Prn5. The constants are to be calculated using the least square method, which is
implemented in some mathematical softwares, such as Mathematica or Matlab.
Using the Equation (5.2) with the constants given in the Table 5.1, the average
Nusselt numbers of the models that dimensions are defined in Table 4.1 can be correlated
with the tolerance of ±10% for Re ≥500 and sL≤63.
Keeping the same form of the correlation equation as above, Equation (5.2), and
using the constants given in Table 5.2, the average Nusselt numbers of all the simulated
models, Table 4.2, can be approximated with the discrepancy of ±20%.
Table 5.1: Constants used in Equation (5.2).
C2n1n2n3n4
0.843 −0.03 −0.68 0.66 0.62
Table 5.2: Constants used in Equation (5.2).
C2n1n2n3n4
0.843 −0.249 −0.322 0.66 0.607
However, in order to achieve a better accuracy in all ranges of varied parameters
with, a complicated function for calculating the average Nusselt number, NumW , for
71
5 Results
example a five-dimensional polynomial
NumW =
i=I
X
i=0
j=J
X
j=0
k=K
X
k=0
l=L
X
l=0
m=M
X
m=0
ci,j,k,l,m δ
LisT
LjsL
LkRelPrm,(5.3)
would be necessary, where I,J,K,Land Mrepresent the polynomial degree for coordi-
nates δ/L,sT/L,sL/L,Re and Pr, respectively and ci,j,k,l,m are constants, that are also
in this case to be calculated through the least squares method, based on the simulation
results for the average Nusselt numbers, given in Appendix F, Table F.1.
All the numerical results for the Nusselt number can be approximated using Equa-
tion (5.3) with a deviation of ≤10% by setting I=J=K=L= 2, M= 0 and
m= 0, which results in 81 constants, that are to be calculated, Appendix E, Table E.
This would be, however, scarcely applicable in the engineering practice.
5.1.3.5 The thermo-fluid characteristic
The complex geometry of the thermoplate channel results in a larger heat transfer coef-
ficient in comparison to the channel with plane walls. This heat transfer improvement
is to be paid by a larger pressure drop, that is to say, by a larger pumping power at
given volume flow rate. To illustrate this interplay, one combines the heat flow rate,
˙
QW, and the pumping power, P, to a thermo-fluid characteristic of the heat transfer
device, see e.g. Stephan and Mitrovi´c [62], and Webb and Kimm [71]. For instance, the
ratio ˙
QW/P specifies the running costs of the heat transfer.
As an example, Figure 5.21 illustrates the quantity ˙
QW/P in dependence on the
plate geometry with Re as parameter. As is obvious from this Figure, the ration ˙
QW/P
increases with increasing sL/L, and the curves move up with rising Reynolds number.
However, at each Re, this ratio is larger for the channel having plane walls (sL/L → ∞).
This behaviour is emphasized in Figure 5.22 that shows the ratio of the thermo-fluid
characteristics ( ˙
QW/P)T P and ( ˙
QW/P)F C of the thermoplate and of the flat channel with
plane walls as function of the geometrical parameters for selected Reynolds numbers, Re.
The ratio is seen to increase with increasing sL/L, but it still remains smaller than unity,
the curve for Re = 50 at larger sL/L being an exception.
In general, at larger Reynolds numbers the thermoplate is more effective than at
smaller ones. Furthermore, the ratio ( ˙
QW/P)T P /(˙
QW/P)F C is smaller at smaller sL/L
and larger δ/L. Its dependency on sT/L is more complicated because of the maxima at
larger Reynolds numbers.
Summarising the results on thermo-fluid characteristics, one may conclude that the
flat channel formed by plane walls shows a better heat transfer characteristic, expressed
in terms of pumping power, than the corresponding thermoplate. However, a fabrication
of such a channel without welding spots seems scarcely possible in sizes required in
common practice and, from this point of view, the parallel plate channel serves merely
as a model for comparison purposes.
72
5.1 Symmetric heated channels, steady state flows
Figure 5.21: Thermo-fluid characteristic of a 500 mm long thermoplate channel as
function of geometry at selected Reynolds numbers; sL/L → ∞: parallel plates .
73
5 Results
Figure 5.22: Thermo-fluid characteristic of the thermoplate compared with the char-
acteristic of the flat channel at selected Reynolds numbers.
74
5.1 Symmetric heated channels, steady state flows
5.1.4 Welding spots pattern and thermoplate surface shape
The thermoplates with welding spots arranged in a staggered manner, whose surface is
described by employing of the trigonometrical functions, Equation (4.1), Figure 4.1a,
were analysed and described in previous chapters. However, thermoplates with different
welding spots arrangements and surface shapes have not been compared yet.
Figure 5.23: Distribution of the local heat transfer coefficient, Equation (4.33), for
different models, Figure 4.1a, c and e, at Re = 50.
Table 5.3: Volume flow rate comparison for different models, Figure 4.1.
Model
Re
50 100 500 1000 2000
uIN ˙
V uIN ˙
V uIN ˙
V uIN ˙
V uIN ˙
V
[m
s] [mm3
s] [m
s] [mm3
s] [m
s] [mm3
s] [m
s] [mm3
s] [m
s] [mm3
s]
1 0.011 781 0.022 1557 0.112 7786 0.224 15572 0.448 31144
2 0.011 1173 0.022 2346 0.112 11728 0.223 23456 0.447 46911
3 0.011 782 0.022 1563 0.112 7817 0.223 15634 0.447 31269
4 0.011 376 0.022 753 0.112 3764 0.224 7529 0.447 15057
5 0.006 706 0.012 1412 0.062 7060 0.124 14121 0.247 28242
Figures 5.23 and 5.24 show the local heat transfer coefficient, calculated according
to Equation (4.33) for models 1, 3 and 5, Figure 4.1a, c and e. These three models
are compared, because not only the overall geometrical parameters such as sL,sT,δ
75
5.1 Symmetric heated channels, steady state flows
Figure 5.26: Average Nusselt number at different Reynolds numbers for different mod-
els, Figure 4.1.
Figure 5.27: Pumping power number at different Reynolds numbers for different mod-
els, Figure 4.1.
and L, as well as the Reynolds numbers are the same, but also the volume flow rate
differs not significantly, Table 5.3. The channel with staggered welding spots (model
1, Figure 4.1a) seems to be the most effective one, whereas the channel with surface
approximated using the exponential function (4.4), model 5 in Figure 4.1e, gives the
lowest value for the heat transfer coefficient. The local heat transfer coefficient of the
channel with parallel welding spots, model 3 in Figure 4.1c, has the largest oscillation
amplitude, but its average value lies in between the other two, Maleti´c and Mitrovi´c
[31].
The models 2, 3 and 4 have basically the same geometry. They are only shifted in
77
5 Results
Figure 5.28: Thermo-fluid characteristic at different Reynolds numbers for different
models, Figure 4.1.
xdirection and have different inlet flow areas when other geometrical parameters (sL,
sT,δ,dhand L) are kept constant. For this reason, the volume flow rates for models 2,
3 and 4 are different at the same Reynolds number, Table 5.3, which causes differences
in heat transfer coefficient, Figure 5.25.
The average Nusselt number, Equation (4.40), and pumping power, Equation (4.45),
for all five simulated models (Figure 4.1) at different Reynolds numbers are given in
Figures 5.26 and 5.27, respectively. From Figures 5.26 and 5.27 it can be seen, that the
heat transfer enhancement has to be paid by higher pumping power.
Combining the heat transfer rate and pumping power into a thermo-fluid character-
istic, as it was done in the previous section, one can get some insights into the running
costs of the heat transfer in different models, Figure 5.28. Although the channel, whose
heat transfer surface is approximated by Equation (4.4) does not show the best heat
transfer characteristic in comparison to other models, the pressure drop in such a chan-
nel is much lower than in other channels. For that reason, the model 5, Figure 4.1e,
seems to be the most effective one.
At this place it has to be noted, that the thermo-fluid characteristic defined as
˙
Q/P could be misleading, when considered alone without the heat transfer rate ˙
Qand
pumping power P. Although the heat transfer rate ˙
Qand pumping power Phave the
same unit, the prices of these two arts of energy are often quite different and other costs,
like storage costs, were not taken into account in this analysis.
78
5.2 Symmetric heated channels, transient flows
5.2 Symmetric heated channels, transient flows
Models with unsteady fluid flow and heat transfer are considered in this section, as it
was stated in Table 4.4. Since the simulations of the time dependent flows are rather
time consuming, only three numerical experiments with the same thermoplate geometry,
but different definitions of the fluid velocity at the inlet were made.
To satisfy the condition that the Courant number is always smaller than unity,
Cr=¯u∆t
∆l≤1,(5.4)
which is required in numerical calculations of time dependent flows, the average inlet
velocity was chosen to be ¯uIN = 0.1m/s, the time step size ∆t= 0.001 s, which corre-
sponds to the Reynolds number Re = 450 and one thousand time steps for the duration
of the simulation of 1 s. ∆ldenotes the local cell size.
Development of the velocity field in the symmetry plane, y= 0, with time at
different time steps, 0 s < t < 1 s is visualised in Figure 5.29. At the very beginning of
the time interval, t= 0.05 s, one can notice the fluid separation behind of each welding
spot. There arise vortices, the centres of which move downstream with the time, until
the recirculation zones occupy whole regions between two neighbouring welding spots.
Tis happens already at t= 0.4 s and from this time step the fluid flow can be considered
as hydraulically fully stabilised, i.e. developed with the time.
The situation is much different, if the development of the temperature field in the
symmetry plane, y= 0 with the time is considered, Figure 5.30. Large portion of fluid
remains thermally unaffected until t= 0.2 s, with small regions behind the welding
spots being exceptions. This was expected, since the channel height and fluid velocity
the in these regions are low. The temperature field develops further, also after t= 0.4 s,
although the fluid flow was stabilised until this time step. Much the same conclusion
can be deduced after analysing of the heat flux time development at the thermoplate
surface, which is illustrated in Figure 5.31. The heat flux maximum is placed at the
inlet, x≈0, because of the very high temperature gradient. Apart from this one, the
local heat flux maxima are placed in the central section of the thermoplate,z≈sT/2,
i.e. between the two rows of the welding spots, where the fluid velocity is high, and
immediate before every welding spot, where the channel height converges and cold fluid
flows with high velocity, resulting in higher temperature gradients.
A comparison between the average heat transfer coefficients for steady and unsteady
flows is given in Figure 5.32. It can be seen that the average heat transfer coefficient
for the unsteady flow tends asymptotic toward the one of steady state flow. However,
the period of 2 s is not enough for the fluid flow to be completely stabilised from the
thermal point of view.
Figure 5.33 illustrates the time evolution of the local heat transfer coefficient. At
the beginning of the time interval, it has high values with very short developing region.
79
5 Results
Figure 5.29: Development of the velocity field in the symmetry plane y= 0 with time.
80
5.2 Symmetric heated channels, transient flows
Figure 5.30: Development of the temperature field in the symmetry plane y= 0 with
time.
81
5 Results
Figure 5.31: Development of the heat flux at the thermoplate wall with time.
82
5.2 Symmetric heated channels, transient flows
Figure 5.32: Heat transfer coefficient averaged in space for a 120 mm long model,
steady and unsteady flow.
As the time passes, the developing region grows and the heat transfer coefficient values
fall down with a tendency to reach the corresponding distribution of the steady state
flow.
Although the boundary layer theory is only conditionally adequate to study the
fluid flow and heat transfer in closed channels, under simplified conditions it can give
some insights into the difference in time periods needed for the fluid flow to be stabilised
hydraulically and thermally. Considering transient boundary layer along a flat plate, i.e.
Blasius and Rayleigh problems, see e.g. Schlichting and Gersten [58], the velocity and
temperature boundary layer thickness are described as
δw=c1√νt, (5.5)
δT=c2√κt, (5.6)
where ν,κand tdenote kinematic viscosity, thermal diffusivity ant time, respectively,
c1and c2are constants.
This means that the velocity and temperature boundary layers do not grow with the
same velocity. If the fluid flow between two parallel flat plates is considered, the velocity
and temperature boundary layers will appear at both plates. The velocity boundary
layers will meet in the middle of such channel at time t1, and the temperature layers
at time t2. Adopting c1≈c2, these two time periods, t1and t2, will be related by the
Prandtl number:
t2=c2
1
c2
2
ν
κt1≈Pr ·t1.(5.7)
Since the Prandtl number for liquids is Pr > 1, the temperature boundary layers obvi-
ously have to meet after the velocity boundary layers met, t2> t1. This means, that the
83
5 Results
Figure 5.33: Local heat transfer coefficient for a 120 mm long model, at different time
steps.
fluid flow will stabilise hydrodynamically first, and then thermally, which was observed
in the numerical experiments in this section.
Analysis of unsteady boundary layer over a wavy plate gives qualitative similar
results, see Saha and Debnath [54].
Figure 5.34: Different inlet fluid velocities.
As a further example of transient fluid flow, the fluid velocity at the inlet was varied
84
5.3 Asymmetric heated channels
from a constant value with time to oscillating functions, as it was stated in Table 4.4:
uIN1= 0.1 [m/s],(5.8)
uIN2= 0.1 (1 + 0.2 (sin(100t))) [m/s],(5.9)
uIN3= 0.1+0.04
∞
X
n=0
H(t−nπ
100) [m/s],(5.10)
where H(t) represents the Heaviside step function. These three velocity distributions
are visualised in Figure 5.34.
As a result, the heat transfer coefficient averaged in space, but not in time, hmW ,
was compared and this comparison is illustrated in Figure 5.35. The curves representing
the heat transfer distributions with oscillating fluid velocities at the inlet also oscillate
around the one with constant velocity. Since the difference in heat transfer coefficients
for transient flows averaged with time is smaller than 0.3%, it can be said, that the
fluctuations of the inlet fluid velocity defined by Equations (5.9) and (5.10), where the
frequency and amplitude of the velocity oscillations were kept constant, do not affect
the heat transfer performance of the thermoplate.
5.3 Asymmetric heated channels
In Chapter 3 it was explained in detail how the asymmetrical boundary conditions influ-
ence the heat transfer and the Nusselt number distribution in a parallel plate channel.
For that simple channel geometry the energy equation in the dimensionless form (3.12)
Figure 5.35: Average heat transfer coefficient (averaged in space) of a 300 mm long
model for different inlet fluid velocities.
85
5 Results
with the parabolic velocity profile could be solved analytically. Due to the definition
of the heat transfer coefficient and depending on the thermal asymmetry, the Nusselt
number could have a vertical asymptote or become zero.
Figure 5.36: Nusselt number distributions at different ΘW2; ΘW1= 0, Re = 100. The
thin lines represent the corresponding distributions for the parallel plate channel, Figure
3.13.
The similarities and differences in heat transfer in the parallel plate channel and
thermoplate are discussed in this section. One thermoplate with the geometrical para-
meters and thermally asymmetric boundary conditions specified in Table 4.5 was chosen
to compare the heat transfer character in thermoplate with the one in parallel plate
channel. The Nusselt number distributions for the lower and upper thermoplate wall,
NuW1and NuW2respectively, are illustrated in Figure 5.36. The thin lines represent
the corresponding distributions for the parallel plate channel, as it was shown in Figure
3.13, in Chapter 3.
The Figure 5.36 shows that the Nusselt number distributions, NuW1or NuW2, can
have a vertical asymptote and become zero also in case of a thermoplate, that is exposed
to thermally asymmetric boundary conditions, for certain combinations of the boundary
conditions TIN ,TW1and TW2, i.e. ΘW2, similar as in the case of the parallel plate
channel. The conditions for this to happen are the same as in the case of the parallel
plate channel, Section 3.3.1. Since the Nusselt number distribution in case of parallel
86
5.3 Asymmetric heated channels
Figure 5.37: Nusselt number distributions at different ΘW2; ΘW1= 0, Re = 100. The
thin lines represent the corresponding distributions for the parallel plate channel.
plates could become zero only once, in case of thermoplate it can become zero more than
once, due to its oscillating character.
Although the Equations (3.44), (3.46) and (3.48) can not be used to calculate the
positions where the Nusselt number experiences a vertical asymptote and zero in case
of thermoplate, it was observed that the dependency 3.49:
(1 −ΘW2(X=−1))(1 −ΘW2(X=+1)) = 1,(5.11)
holds in case of thermoplate, as well.
In Section 3.3 it has been observed, that in case of parallel plates the Nusselt
number distributions, NuW1and NuW2, falling from the infinity follow the limit of
Nu∞= 7.545 (symmetric case) in the region before the asymptote, i.e. 0 < Z∗< ZA∗,
first, and than the limit of Nu∞= 4.0, for Z∗→ ∞. The numerical simulations show
that the corresponding distributions of Nusselt number oscillate around the same value
before and after the vertical asymptote in case of heat transfer in thermoplate, Figure
5.36. This is an advantage of the thermoplate, since it seems to be not sensitive to the
thermal asymmetry, as it is the case with the parallel plate channel, and the thermal
asymmetry is rather the rule than the exception in the practice.
Furthermore, the Nusselt number distributions for the parallel plate channel does
not depend on the Reynolds number, because its influence is taken into account through
87
5 Results
Z∗, Equation (3.11). The diagram is only scaled in Z∗direction when the Reynolds
number is changed. The situation changes in chase of thermoplates, where the Nusselt
number distribution is influenced by the Reynolds number, not only through the scaling
in the Z∗direction. In case of thermoplate, the average Nusselt number value increases,
when the Reynolds number is raised, Figure 5.37.
5.4 Validation through Experiments
Despite various uses of thermoplates in heat exchangers, no relevant experimental data
on fluid flow and heat transfer in thermoplates have been reported in the literature, the
paper by Mitrovi´c and Peterson [38] being an exception.
Since the part I of this paper deals with single phase forced convection heat transfer
and pressure drop in a thermoplate, the data published there were convenient for the
comparison with the numerical results of this thesis. The thermoplate made of 0.8 mm
thick stainless steel sheets was installed vertically. It dimensions (width x height x
thickness) were 300 mm x 1000 mm x 5 mm, Figure 5.38. The test plate was switched
as a resistance heater into an electric circuit and the measured electrical resistance of the
plate was used to determine its wall temperature. The plate was electrically disconnected
from the rest of the plant by means of POM (Polyoxymethylene). For other details see
the paper by Mitrovi´c and Peterson, [38].
Due to the fact, that the experiments on heat transfer and pressure drop have been
made using only one thermoplate geometry, the calculated Nusselt number and pressure
drop values could be validated through experiments only as function of the Reynolds
number and mass flow rate, respectively, and the thermoplate geometry could not be
varied.
Marlotherm oil and distilled water were used as test fluids, because the electrical
resistance of these fluids is high enough to allow the application of the electrical resistance
heating method. The fluid properties of the Marlotherm oil and distilled water were
adopted to be temperature dependent.
From the variety of the experimental runnings, for the comparison with the numeri-
cal results were chosen only those, where the inlet fluid temperature and the thermoplate
wall temperature were kept approximately constant for all Reynolds numbers, Table 5.4,
(Marlotherm oil at Pr = 8.4 and water at Pr = 3.9 in paper by Mitrovi´c and Peterson
[38]).
The influence of the connecting tubes and thermoplate edges was not known in
advance. For that reason two sets of the numerical runnings have been made, one with
the numerical model that consists of a thermoplate strip and the second one was a
model describing the whole thermoplate inclusive the connecting tubes (by employing
of the symmetry planes at y= 0 and z= 0 it was enough to simulate 1/4 of the
thermoplate). The fact, that the maximal plate distance (interior channel hight) δis
88
5.4 Validation through Experiments
Figure 5.38: Geometry of the thermoplate used in experiments by Mitrovi´c and Pe-
terson, [38].
not constant in the reality, but a function of xand zcoordinates, was taken into account
by the numerical calculations. Still, a minor discrepancy between the real thermoplate
geometry and the numerical model could not be avoided. Further, the thermoplate wall
temperature, TW, that was estimated by the experimental runnings represents rather an
average temperature. This temperature was adopted for the wall boundary condition by
the corresponding numerical experiments, because the exact temperature distribution
over the thermoplate surface could not be measured by the available equipment. In
the runnings the temperature difference TW−TIN of only few degrees (less than 12K)
could be established, Table 5.4, which challenges the modelling additionally. Defining a
constant heat flux boundary condition at the thermoplate wall instead of the constant
temperature did not result in better agreement with the experimental results.
89
5 Results
Figure 5.39: Flow field and temperature field distribution across the thermoplate;
Marlotherm oil at Re = 1860, ˙mIN = 510kg/h, TIN = 328.2K, TW= 333.2K.
90
5.4 Validation through Experiments
Figure 5.40: Flow field and temperature field distribution across the thermoplate;
Marlotherm oil at Re = 5416, ˙mIN = 1500kg/h, TIN = 328.1K, TW= 337.2K.
91
5 Results
Table 5.4: Boundary conditions for the comparison with the experiments of Mitrovi´c
and Peterson.
Run
Marlotherm Oil Water
˙mIN [kg/h] TIN [K] TW[K] ˙mIN [kg/h] TIN [K] TW[K]
1 510 328.2 333.2 455 314.2 325.3
2 1000 327.5 335.8 670 313.4 321.2
3 1280 327.9 336.5 837 314.6 321.0
4 1500 328.1 337.2 1016 314.8 320.1
5 2020 328.3 337.4 1238 315.8 320.5
6 1485 315.4 319.5
7 1684 315.9 319.6
8 1960 315.6 318.6
9 2199 315.9 318.9
By calculations with the thermoplate strip model, the effects of the thermoplate
edges and connecting tubes could not be taken into account and had to be neglected. In
this case the simplification was made that the fluid flows in with the mean velocity uIN
and temperature TIN and passes homogeneously through the thermoplate that is kept
at constant temperature TW. Consequently, the recirculation zones that are expected to
occur in the thermoplate corners and the heat transfer reduction because of this could
not be captured with the thermoplate strip approach. The high temperature gradient in
the developing region (x≈0) over the whole thermoplate width results in locally higher
Nusselt number values (tend to the infinity), than it is expected to be in the reality.
This will have an impact on the overall heat transfer.
Modeling the whole thermoplate shows that there are recirculation zones in the
thermoplate corners near the inlet, but not at the other end, near the outlet, Figures 5.39
and 5.40. The fluid is heated up faster in the recirculation zones which causes reduction
of the heat transfer to some extent. But on the other hand, according to the continuity
equation the fluid is speeded up near the inlet, which increases the heat transfer (low
temperature regions in Figures 5.39 and 5.40). These effects may compensate themselves
if they are of the same order of magnitude. Depending on the Reynolds number, the
former or the latter effect can outweigh. Inspecting the flow fields shown in Figures 5.39
and 5.40) it can be seen, that the fluid hits the welding points with high flow velocity
near the inlet. Removing a few welding points from the near inlet region would reduce
the overall pressure drop, which is the same conclusion that Witry et al. [73, 74] deduced
from their numerical investigations. This could however allow the recirculation zones in
the thermoplate corners to be greater.
Figures 5.41 and 5.42 illustrate calculated Nusselt number compared to the the one
reported by Mitrovi´c and Peterson [38], for the thermoplate strip and whole thermoplate
model respectively. In principle, there is no big difference between these two diagrams.
92
5.4 Validation through Experiments
Figure 5.41: Nusselt number comparison for a thermoplate strip.
Figure 5.42: Nusselt number comparison for whole thermoplate.
The calculated Nusselt number values at lower Reynolds numbers for the thermoplate
strip, Figure 5.41, are somewhat greater than those of the whole thermoplate, Figure
5.42. This means, that there is an influence of the connecting tubes and recirculation
zones in the thermoplate corners only at lower Reynolds numbers, whereas at higher
Reynolds numbers the heat transfer reduction is compensated by the heat transfer in-
crease caused by the local speed up.
Further, Figures 5.41 and 5.42 show a good agreement between calculated and
measured Nusselt number values for the Marlotherm oil. In case of distilled water the
calculated Nusselt number is always clearly lower than the measured one, except for
the low Reynolds numbers. This indicates that the turbulent effects play an important
role at such regimes, which is not surprising. The differences between Nusselt number
93
5 Results
Figure 5.43: Pressure drop comparison for whole thermoplate.
behaviour for Marlotherm oil and water are associated with different fluid properties
of these two fluids. Calculations with the standard k-εmodel resulted in qualitative
similar flow distribution but with much higher Nusselt number values in comparison
to the measured ones. That was the reason for adopting the laminar model for all
calculations presented in this thesis.
A comparison between calculated and measured pressure drop over the thermoplate
as a function of the mass flow rate is given in Figure 5.43. Although the fluid properties
of Marlotherm oil and water are different, the calculated pressure drop values are almost
identical for these two fluids. The calculated values are lower than the measured ones
and the discrepance increases with increasing mass flow rate. This can be attributed to
the small inaccuracy in the geometry of the thermoplate connections with the inlet and
outlet tube.
As it could be seen in this section, a very good agreement between the calculated
and measured values for the Marlotherm oil could be achieved. For the distilled water a
satisfactory agreement could not be reached at all Reynolds numbers. This makes clear
that an appropriate turbulence model is anavoidable for the heat transfer quantification.
Nevertheless, for the comparison purposes and finding an optimal thermoplate geometry,
the methodology shown in this thesis is considered to be adequate.
94
6 Conclusions
This thesis deals with the forced convection and heat transfer in a parallel plate channel
and in thermoplates under thermally symmetrical and asymmetrical boundary condi-
tions.
In case of asymmetrically heated parallel plate channel the dimensionless energy
equation could be solved analytically in form of Hermite polynomials and hypergeometric
functions. Depending on the boundary conditions, the Nusselt number distribution at
the lower or at the upper plate could experience a vertical asymtote and a zero value. A
relation, that describes the dependency of the dimensionless temperatures of the upper
plate, ΘW2(X=−1), for the Nusselt number asymptotes to occur at the same axial position,
ZA∗, at the lower and at the upper plate, respectively, was deduced. Further, it was
shown, that a small thermal asymmetry reduces the Nusselt number from Nu∞= 7.545
(symmetric case) to Nu∞= 4.0 (asymmetric case).
The heat transfer can be increased up to four times in comparison to the parallel
plate channel when thermoplates, whose surface can be described through trigonometric
functions, are used. This enhancement has to be paid through higher pressure drop
and pumping power. The numerical experiments reveal that the staggered welding spot
pattern should be used to achieve better heat transfer. The welding spots pitch in the
flow direction, sL, does not play any important role in this case. Only at higher Reynolds
numbers an optimum of sL≈80mm occurs. The heat transfer is also improved chosing
lower values for the transversal welding spot pitch, sT, and the maximal plate distance,
δ.
If the pressure drop and the pumping power are very important by the implemen-
tation of the thermoplates, a thermoplate that has even surface for the most part, as it
is the case with the model 5, Figure 4.1e, is recommended. This thermoplate surface is
approximated using an exponential function. Using this kind of thermoplate surface will
not increase the heat transfer as much as in the case of using the model 1, Figure 4.1a,
in comparison to the parallel plate channel. However, the low values of the pressure
drop and pumping power result in a convenient thermo-fluid characteristic.
Contrary to the parallel plate channel, in the case of asymmetrical heating after
one of the thermoplate walls experienced a vertical asymptote, thermoplates remain
thermally active as in the region before the vertical asymptote. This distinguishes the
thermoplates, since the thermal asymmetry is rather the rule than an exception in the
practice. The higher Reynolds number results in higher Nusselt number when thermo-
plates are used, which is not the case with the prarallel plate channel.
95
6 Conclusions
A Nusselt number correlation as a function of the Reynolds number and three
geometrical parameters of the thermoplate could be found, where the constants are
calculated using the least squares method.
Some of the numerical results from this thesis could be compared to the laboratory
experiments. A very good agreement was observed when Marlotherm oil was used as
a test fluid. However, the comparison indicates, that the turbulence effects can not be
neglected particularly at higher Reynolds numbers.
96
Bibliography
[1] Abramowitz, M. and Stegun, I.A., 1972, Handbook of mathematical func-
tions,Dover Publications, New York.
[2] Anderson, J.D., 1995, Computational fluid dynamics, McGraw-Hill, Inc., New
York.
[3] Behrend, H.-J., 1993, Thermoblechw¨armetauscher und -apparate: Vielseit-
igkeit bewiesen, Schweizer Maschinenmarkt, 12, 95-97.
[4] Bird, R.B, Stewart, W.E. and Lightfoot, E.N., Transport phenomena (Revised
2nd Ed.), 2007, Wiley, New York.
[5] Blancher, S., Creff, R. and Quere, P.L., 1998, Effect of Tollmien Schlichting
wave on convective heat transfer in a wavy channel. Part I: Linear analysis,
International Journal of Heat and Fluid Flow, 19, 39-48.
[6] Burns, J.C, Parkers, T., 1967, Peristatic motion, Journal of Fluid Mechanics,
29, 731-743.
[7] CD adapco Group, 2004, Methodology StarCD version 3.24.
[8] Cherukat, P., Na, Y., Hanratty, T.J., and Mclaughlin, Y.B. 1998, Direct numer-
ical simulation of a fully developed turbulent flow over a wavy wall, Theoretical
and Computational Fluid Dynamics, 11, 109-134.
[9] Cho, K.J., Kim, M.U. and Shin, H.D., 1998, Linear stability of two-dimensional
steady flow in wavy-walled channels, Fluid Dynamics Research, 23, 349-370.
[10] Deiber, J.A. and Schowalter, W.R., 1979, Flow through tubes with sinusoidal
axial variations in diameter, AIChE Journal, 25, 638-645.
[11] Dellil, A.Z., Azzi, A. and Jubran, B.A., 2004, Turbulent flow and convective
heat transfer in a wavy wall channel, Heat and Mass Transfer, 40, 793-799.
[12] Fan, R., 2007, Numerische Untersuchungen zum W¨arme¨ubergang und Druck-
abfall in Thermoblechen, Bachelorarbeit, University of Paderborn.
97
Bibliography
[13] Ferziger, J.H. and Peri´c, M., 2002, Computational methods for fluid dynamics,
Springer, Berlin.
[14] Garg, V.K. and Maji, P.K., 1988, Laminar flow and heat transfer in a pe-
riodically converging-diverging channel, International Journal for Numerical
Methods in Fluids, 8, 579-597.
[15] Graetz, L., 1883, Ueber die W¨armeleitungsf¨ahigkeit von Fl¨ussigkeiten, Annalen
der Physik NF, 18, 79-94.
[16] Grass, A.J., Stewart, R.J. and Mansour-Tehrani, M.,1991, Vortical structures
and coherent motion in turbulent flow over smooth and rough boundaries,
Philosophical Transactions of the Royal Society of London, Series A, 336, 35-
65.
[17] Grimm, P., 2007, Numerische Untersuchung des W¨arme¨ubergangs an Ther-
moblechen, Studienarbeit, University of Paderborn.
[18] Guzm´an, A.M. and Amon, C.H., 1996, Dynamical flow characterization of
transitional and chaotic regimes in converging-diverging channels, Journal of
Fluid Mechanics, 321, 25-57.
[19] Hatton, A.P and Turton, J.S., 1962, Heat transfer in the thermal entry length
with laminar flow between parallel walls at unequal temperatures, International
Journal of Heat Mass Transfer, 5, 673-679.
[20] Heckel, W., 2007, Numerische Untersuchungen zur W¨arme¨ubertragung in
Thermoblechen, Studienarbeit, University of Paderborn.
[21] Hossain, M.Z. and Islam, A.K.M.S., 2004, Fully developed flow structures and
heat transfer in sine-shaped wavy channels, International Communications in
Heat Mass Transfer, 31, 887-896.
[22] Hossain, M.Z. and Islam, A.K.M.S., 2007, Numerical investigation of fluid flow
and heat transfer characteristics in sine, triangular, and arc-shaped channels,
Thermal Science, 11, 17-26.
[23] Kaka¸c, S., Shah, R.K. and Aung, W., 1987, Handbook of single-phase convec-
tive heat transfer,John Wiley & Sons, New York.
[24] Kim, S.K., 2001, An experimental study of flow in a wavy channel by PIV,
Proceedings of the 6th Asian Symposium on Visualisation, Busan, Korea.
[25] Lahbabi, A. and Chang, H.C., 1986, Flow in periodically constricted tubes:
Transition to inertial and nonsteady flows, Chemical Engineering Science, 41,
2487-2505.
98
Bibliography
[26] Lee, B.S., Kang, I.S. and Lim, H.C., 1999, Chaotic mixing and mass trans-
fer enhancement by pulsatile laminar flow in an axisimmetric wavy channel,
International Journal of Heat Mass Transfer, 42, 2571-2581.
[27] Li, Y., 2007, Numerische Untersuchungen zur Hydrodynamik und zum
W¨rme¨
bergang in Thermoblechen, Bachelorarbeit, University of Paderborn.
[28] Mahmud, S., Islam, A.K.M.S. and Mamun, M.A.H., 2002, Separation charac-
teristics of fluid flow insice two parallel plates with wavy surface, International
Journal of Engineering Science, 40, 1495-1509.
[29] Mahmud, S. and Islam, 2003, Laminar free convection and entropy generation
inside an inclined wavy enclosure, International Journal of Thermal Sciences,
42, 1003-1012.
[30] Mahmud, S., Islam, A.K.M.S. and Feroz, C.M., 2003, Flow and heat transfer
characteristic inside a wavy tube, Heat and Mass Transfer, 39, 387-393.
[31] Maleti´c, B. and Mitrovi´c, J., 2008, Influence of the thermoplate geometry on the
heat transfer, Proceedings of the 5th European Thermal-Sciences Conference
EUROTHERM 2008, HEX27, Eindhoven, The Netherlands.
[32] Mitrovi´c, J., Maleti´c, B. and Baˇcli´c, B.S., 2006, Some peculiarities of the asy-
metric Graetz problem, International Journal of Engineering Science, 44, 436-
455.
[33] Mitrovi´c, J. and Maleti´c, B. 2005, Effect of thermal assymetry on heat transfer
in a laminar annular flow, Chemical Engineering and Technology, 28, 1144-
1150.
[34] Mitrovi´c, J. and Maleti´c, B., 2006, Numerical simulation of fluid flow and heat
transfer in thermoplates, Proceedings of the 13th International Heat Transfer
Conference IHTC-13, HEX-10, Sydney, Australia.
[35] Mitrovi´c, J. and Maleti´c, B., 2007, Numerical simulation of fluid flow, heat
transfer and pressure drop in thermoplates, Proceedings of the 5th Interna-
tional Conference on Heat Transfer, Fluid Mechanics and Thermodynamics
HEFAT2007, MJ3, Sun City, South Africa.
[36] Mitrovi´c, J. and Maleti´c, B., 2007, Numerical simulation of fluid flow and heat
transfer in thermoplates, International Journal of Heat Exchangers, in Press.
[37] Mitrovi´c, J. and Maleti´c, B., 2007, Heat transfer with laminar forced convection
in a porous channel exposed to a thermal asymmetry, International Journal of
Heat and Mass Transfer, 50, 1106-1121.
[38] Mitrovi´c, J. and Peterson, R., 2007, Vapor condensation heat transfer in ther-
moplate heat exchanger, Chemical Engineering and Technology, 30, 907-919.
99
Bibliography
[39] M¨uhlthaler, W., 2002, Thermo-Blech W¨armeaustauscher-Systeme, Absolven-
ten Aktuell, 13, 2-3.
[40] Niˇceno, B. and Nobile, E., 2001, Numerical analysis of fluid flow and heat
transfer in periodic wavy channels, International Journal of Heat and Fluid
Flow, 22, 156-167.
[41] Nield, D.A., 2004, Forced convection in a parallel plate channel with asymmet-
ric heating, International Journal of Heat and Mass Transfer, 47, 5609-5612.
[42] Nishimura, T., Murakami, S., Arakawa, S. and Kawamura Y., 1990, Flow
observation and mass transfer in sinusoidal wavy-walled channel at moderate
Reynolds numbers for steady flow, International Journal of Heat and Mass
Transfer, 33, 835-845.
[43] Nishimura, T., Bian, Y.N., Matsumoto, K. and Kunitsugu K., 2003, Fluid
flow and mass transfer in a sinusoidal wavy-walled tube at moderate Reynolds
numbers for steady flow, Heat and Mass Transfer, 39, 239-248.
[44] Nishimura, T., Arakawa, S., Murakami, S. and Kawamura Y., 1989, Oscillatory
viscous flow in symmetric wavy-walled channels, Chemical Engineering Science,
44, 2137-2148.
[45] Nishimura, Murakami, S. and Kawamura Y., 1993, Mass transfer in a symmet-
ric sinusoidal wavy-walled channel for oscillatory flow, Chemical Engineering
Science, 48, 1793-1800.
[46] Nishimura, T., Ohori, Y. and Kawamura, Y.,1984, Flow characteristics in a
channel with symmetric wavy wall for steady flow, Journal of Chemical Engi-
neering of Japan, 17, 466-471.
[47] Nishimura, T., Ohori, Y., Kajimoto, Y. and Kawamura, Y.,1985, Mass transfer
characteristics in a channel with symmetric wavy wall for steady flow, Journal
of Chemical Engineering of Japan, 18, 550-555.
[48] Nishimura, T. and Matsune, S., 1996, Mass transfer enhancement in a sinu-
soidal wavy channel for pulsatile flow, Heat and Mass Transfer, 32, 65-72.
[49] Nishimura, T. and Matsune, S., 1998, Vortices and wall shear stresses in asym-
metric and symmetric channels with winusoidal wavy walls for pulsatile flow
at low Reynolds numbers, International Journal of Heat and Fluid Flow, 19,
583-593.
[50] Petukhov, B.S., 1967, Laminar fluid flow and heat transfer in ducts (in
Russian), Energiya, Moscow.
[51] Ralph, M.E., 1987, Steady flow structures and pressure drops in wavy-walled
tubes, Journal of Fluids Engineering, 109, 255-261.
100
Bibliography
[52] Rush, T.A., Newell, T.A. and Jacobi, A.M., 1999, An experimental study of
fluid flow and heat transfer in sinusoidal wavy passages, International Journal
of Heat and Mass Transfer, 42, 1541-1553.
[53] Russ, G. and Beer, H., 1997, Heat transfer and flow field in a pipe with sinu-
soidal wavy surface - I. Numerical investigation, II. Experimental investigation,
International Journal of Heat and Mass Transfer, 40, 1061-1081.
[54] Saha, L.M. and Debnath, L., 1973, Unsteady boundary layer flows induced by
a wavy plate, Pure and Applied Geophysics, 105, 802-809.
[55] Saniei, N. and Dini, S., 1993, Heat transfer characteristics in wavy-walled chan-
nel, ASME Journal of Heat Transfer, 115, 788-792.
[56] Sawyers, D.R., Sen, M. and Chang, H.-C., 1998, Heat transfer enhancement in
three-dimensional corrugated channel flow, International Journal of Heat and
Mass Transfer, 41, 3559-3573.
[57] Selvarajan, S., Tulapurkara, E.G. and Vasanta Ram, V., 1998, A numerical
study of flow through wavy-walled channels, International Journal for Numer-
ical Methods in Fluids, 26, 519-531.
[58] Schlichting, H. and Gersten, K., 1997, Grenzschicht-Theorie, Springer, Berlin.
[59] Shah, R.K. and London, A.L., 1978, Laminar flow forced convection in ducts,
supplement to advances in heat transfer, Academic Press, New York.
[60] Slater, L.J., 1960, Confluent hypergeometric functions, Cambridge University
Press, Cambridge.
[61] Sparrow, E.M. and Prata, A.T., 1983, Numerical solutions for laminar flow and
heat transfer in a periodically converging-diverging tube, with experimental
confirmation, Numerical Heat Transfer, Part A, 6, 441-461.
[62] Stephan, K. and Mitrovi´c, J., 1984, Maßnahmen zur Intensivierung des
W¨arme¨ubergangs, Chemie - Ingenieur - Technik, 56, 427-431.
[63] Stone, K. and Vanka, S.P., 1999, Numerical study of developing flow and heat
transfer in a wavy passage, Journal of Fluids Engineering, 121, 713-719.
[64] URL. DEG-Engineering GmbH. http://www.deg-engineering.de/
produktee.htm
[65] URL. GEA FlatPlate Inc.. http://www.heatexchangers.org/info/
heatexchangers/
[66] URL. OMEGA Heat Transfer Technology. http://www.heseco.com/
platecoil.htm
101
Bibliography
[67] Utamura, M., Nikitin, K. and Kato, Y., 2007, Generalization of logarithmic
mean temperature difference method for heat exchangers performance analysis,
Thermal Science and Engineering, 15, 163-173.
[68] Versteeg, H.K. and Malalasekera, W., 1995, An introduction to computational
fluid dynamics - The finite volume method, Prentice Hall, Harlow.
[69] Wang, C.C. and Chen, C.K., 2002, Forced convection in a wavy-wall channel,
International Journal of Heat and Mass Transfer, 45, 2587-2595.
[70] Wang, G. and Vanka S.P., 1995, Convective heat transfer in periodic wavy
passages, International Journal of Heat and Mass Transfer, 38, 3219-3230.
[71] Webb, R.L. and Kim, N.-H., Principles of enhanced heat transfer (2nd Ed.),
2005, Taylor & Francis, New York.
[72] Wirtz, R.A., Huang, F. and Greiner, M., 1999, Correlation of fully developed
heat transfer and pressure drop in a symmetrically grooved channel, Transac-
tions of the ASME, 121, 236-239.
[73] Witry, A., Al-Hajeri, M.H. and Bondok, A.A, 2005, Thermal performance of
automotive aluminium plate radiator, Applied Thermal Engineering, 25, 1207-
1218.
[74] Witry, A., Al-Hajeri, M.H. and Bondok, A.A, 2003, CFD Analyses of fluid
flow and heat transfer in patterned roll-bonded Aluminium plate radiators,
Proceedings of the 3rd International Conference on CFD in the Minerals and
Process Industries CSIRO, Melbourne, Australia.
102
A Numerical Method
As it is well-known from the literature, an analytical solution to the set of the governing
equations (4.17) to (4.21) can be obtained only for very simplified flow situations. How-
ever, for the engineering practice it is often of big importance to solve this set of partial
differential equations employing also complex bounadary conditions, which is only possi-
ble if a numerical method is used. For that reason, a number of numerical methods (the
finite difference method, the finite volume method, the finite element method, the vol-
ume of fluid and so on) and CFD1codes (FLUENT, StarCD, OpenFOAM, PHOENICS,
FLOW3D and others) have been developed. For the numerical calculations presented
in this thesis the commertial software StarCD that uses the finite volume method was
used. In the next sections a short overview of this method will be given; for more de-
tails the original literature sources such as the books by Versteeg and Malalasekera [68],
Anderson [2] and Ferziger and Peric [13] are recommended.
A.1 The finite volume method
If we consider the governing Equations (4.17) to (4.21) it is clear that there is a possibility
to write them in one common form introducing a general variable φ
∂(ρφ)
∂t +div(ρφw) = div(Γφgradφ) + Sφ.(A.1)
The governing equations (4.17) to (4.21) are obtained from (A.1) by setting φequal to
1, u,v,wor i(i.e. T) and selecting appropriate values for the diffusion coefficient Γ and
the source term Sφ.
In the finite volume method, the Equation (A.1) is used as a basic equation for
computationoal procedures. The crucial step of this method is the integration of the
Equation (A.1) over a three-dimensional control volume CV yielding
ZCV
∂(ρφ)
∂t dV +ZCV
div(ρφw)dV =ZCV
div(Γφgradφ)dV +ZCV
SφdV. (A.2)
Using the Gauss’ divergence theorem, the convective term and the diffusive term can
be rewritten as surface integral over the entire bounding surface of the control volume,
1Computational Fluid Dynamics
103
A Numerical Method
and for the steady state problems, the Equation (A.1) becomes
ZA
nρφwdA =ZA
nΓφgradφdA +ZCV
SφdV, (A.3)
whereas for the time dependent problems, an integration with respect to time tover a
small interval ∆tis necessary, which yields
Z∆t
∂
∂t ZCV
(ρφ)dV +Z∆tZA
nρφwdA =Z∆tZA
nΓφgradφdA +Z∆tZCV
SφdV, (A.4)
where the order of integration and differentiantion in the first term on the left hand
side has been changed, which is allowed, if the control volume is independent on time.
Steady state calculations will be considered first, and than time dependent ones.
A.1.1 Steady state calculations
The division of the domain into discrete control volumes, i.e. numerical cells, is the first
important step in the finite volume method. The usual notation convention of CFD
methods is shown in Figure A.1. A cell characterised by the node Pis surrounded by
Figure A.1: A three-dimensional cell with the node Pand neighbouring nodes W,E,
S,N,Band T, [68].
six neighbouring cells and their nodes are identified as west, east, south, north, bottom
and top (W,E,S,N,B,T). The six cell faces are denoted analogical by w,e,s,n,b
and t.
104
A.1 The finite volume method
The second important step of the finite volume method is the integration of the
governing Equation (A.2) over a control volume to get a discretised equation at its
nodal point P, which has the following general form:
aPφP=Xanbφnb +Su=aWφW+aEφE+aSφS+aNφN+aBφB+aTφT+Su,(A.5)
where the source term from the Equation (A.2) is taken to have a linear form:
ZCV
SφdV =¯
S∆V=Su+Spφp,(A.6)
and the coefficients aisatisfy the following relation
ap=Xanb + ∆F−Sp,(A.7)
where Frepresents the convective mass flux per unit area at cell face. The calculation
of the coefficients ap,anb and Fiand the system of the algebraic equations describing
the physical situation depend on the differencing scheme that is used.
A.1.1.1 Differencing schemes
A numerical solution to the transport equation is, in theory, indistinguishable from its
exact solution only if the number of computational cells is infinitely large. However, in
practical calculations only a finite number of cells can be used and the numerical results
will only be physically realistic when the discretisation scheme has certan fundamental
properties, and the most important ones are:
•Conservativeness. The flux of φleaving a control volume through a certain cell
face must be equal to the flux of φentering the adjacent control volume through
the same face, to ensure conservation of φ.
•Boundedness. The differencing scheme should produce such coefficients ai, such
that the matrix of the set of the algebraic equations is diagonally dominant. This
ensures that in absence of sources, the internal nodal values of the property φare
bounded by its boundary values.
•Transportiveness. It is very important that the differencing scheme is not influ-
enced by the direction and rate of the convective flow, which can be expressed in
terms of the dimensionless cell P´eclet number, that is defined as the ratio of the
convective mass flux per unit area and the diffusion conductance at cell faces. For
one-dimensional calculations P´eclet number becomes
Pe =F
D=ρu
Γ/∆x,(A.8)
where ∆xrepresents characteristic length, i.e. cell width in xdirection.This inter-
play is ilustrated in Figure A.2 as the dependency of a countour of constant φon
P´eclet number.
105
A Numerical Method
Figure A.2: Contours of constant φat different P´eclet numbers [68].
Depending on how the convective terms Fon the cell faces are calculated, there are
a few different differencing schemes:
•The central differencing scheme. To compute the cell face value the linear interpo-
lation between the two neighbouring cell nodes is used. This differencing scheme
has second-order accuracy in terms of Taylor series truncation error, but it is sta-
ble and accurate only if Pe < 2, which means only in diffusion-dominated low
Reynolds number flows or if the grid spacing is small, because it is not able to
identify the flow direction. This creates the need for other differencing schemes.
•The upwind differencing scheme. In convective flow, the cell face value is strongly
influenced by the value of property φat the upstream cell node. In this differencing
scheme it is taken that the value at the cell face is equal to the value at the upstream
node, so taking into account the flow direction, Figure A.3. This differencing
scheme has only first order accuracy and produces erroneous results when the flow
is not aligned with the grid lines. The error has a diffusion-like appearance and is
called false diffusion. However, because of its simplicity, the upwind differencing
scheme is widely used in CFD calculations.
•The hybrid differencing scheme. This differencing scheme represents a combination
of the central differencing scheme for Pe < 2 and the upwind scheme for Pe ≥2
•The power low differencing scheme is more accurate than the hybrid scheme. For
Pe ≥10 the diffusion is set to zero, and for 0 < Pe < 10 the flux at the cell face
is calculated according to a polynomial expression.
•The QUICK2differencing scheme uses a tree-point upstream-weighted quadratic
interpolation for cell face values. This scheme has greater accuracy than the central
differencing or hybrid schemes, but it can give some unterschoots and overshoots
and in complex flow calculations the stability problems may occur.
•Other higher order schemes were developed for calculating of convective terms, but
the implementation of the boundary conditions can be problematic.
2Quadratic Upstream Interpolation for Convective Kinetics
106
A.1 The finite volume method
Figure A.3: The upwind differencing scheme for one-dimensional flows [68].
A.1.1.2 The SIMPLE algorithm
In the previous section (A.1.1.1) it was assumed that the velocity field was known, but
in general it is not the case and the velocity field has to be computed. It is governed by
the momentum equations (4.18) to (4.20) and must also satisfy the continuity equation
(4.17). If the flow is compressible, the continuity equation is used as transport equation
for density and the pressure can be obtained employing the equation of state (4.9).
Otherwise, if the flow is incompressible, the density is constant and not linked to the
pressure.
The momentum equations are coupled, contain non-linear quantities and the pres-
sure gradient appears in all of them, without any additional equation for pressure. These
problems can be solved by adopting an iterative strategy such as the SIMPLE3algorithm
of Patankar and Spalding.
Here it is to be noted that if pressure and velocity components are calculated at
the same nodes, a highly irregular pressure field, i.g. checker-board pressure field, could
appear as constant [68]. One of the possibilities to solve this problem is to use a staggered
3Semi-Implicit Method for Pressure-Linked Equations
107
A Numerical Method
grid for the velocity components, Figure A.4. In that way, the pressure and other scalar
values are calculated at the node I, J, whereas the velocity components are calculated
at other cells, i, J and I, j, whose nodes coincide with the faces of the cell I, J. There
are, however, also other possibilities to solve this problem [13] and the StarCD software
does not use the staggereg grid.
In this approach, where velocity and pressure fields are not known, the pressure
gradient term of the momentum equations, which forms the main momentum source
term, is written as an additional source term of the discretised momentum equations.
In two-dimensional case the discretised momentum equations become
ai,J ui,J =Xanbun,b −pI,J −pI−1,J
∆xu
∆Vu+¯
S∆Vu,(A.9)
aI,jvI,j =Xanbvn,b −pI,J −pI,J−1
∆xv
∆Vv+¯
S∆Vv,(A.10)
where ∆Vuand ∆Vvdenote the volume of the uand v-cell, respectively.
The SIMPLE procedure is initiated after the velocity field ˜u∗and ˜v∗is guessed.
Now the coefficients ai,j and anb can be calculated using any of the differencing schemes
suitable for a concrete simulation. With the guessed pressure field p∗the discretised
momentum equations (A.9) and (A.10) can be solved for the new values u∗and v∗
ai,J u∗
i,J =Xanbu∗
n,b −p∗
I,J −p∗
I−1,J
∆xu
∆Vu+¯
S∆Vu,(A.11)
aI,jv∗
I,j =Xanbv∗
n,b −p∗
I,J −p∗
I,J−1
∆xv
∆Vv+¯
S∆Vv.(A.12)
Introducing the pressure and the velocity corrections p0,u0and v0the correct values
for the pressure and the velocity comopnents become
p=p∗+p0,(A.13)
u=u∗+u0,(A.14)
v=v∗+v0.(A.15)
Combining the equations (A.9) to (A.15), the expressions for the calculations of the
velocity corrections u0and v0are obtained as follows
ai,J u0
i,J =Xanbu0
n,b −p0
I,J −p0
I−1,J
∆xu
∆Vu+¯
S∆Vu,(A.16)
aI,jv0
I,j =Xanbv0
n,b −p0
I,J −p0
I,J−1
∆xv
∆Vv+¯
S∆Vv.(A.17)
108
A.1 The finite volume method
The main approximation of the SIMPLE algorithm is made by omitting the summa-
tion terms Panbu0
n,b and Panbu0
n,b in the Equations (A.16) and (A.17). The velocity
components uand vare than
ui,J =u∗
i,J +p0
I−1,J −p0
I,J
∆xuai,J
∆Vu+¯
S∆Vu,(A.18)
vI,j =v∗
I,j +p0
I,J−1−p0
I,J
∆xvaI,j
∆Vv+¯
S∆Vv.(A.19)
So far only the momentum equations have been considered, but the velocity com-
ponents should also satisfy the continuity equation (4.17), which, combined with the
eqations (A.16) and (A.17), results in the following relation
aI,J p0
I,J =Xanbp0
n,b +a0
I,J ,(A.20)
where ai,j = (ρdA)i,j,dA representing the cell face surface and the source term a0
i,j
denotes the imbalance caused by the incorrect velocity field u∗and v∗. The Equation
(A.20) can be solved for p0and is called the pressure correction equation.
Setting p∗=p, ˜u∗=uand ˜v∗=vthis procedure is repeated until the converged
solution is reached and the whole method is illustrated in Figure A.5.
The method can be improved by introducing of the unter-relaxation factors, since
the correction equations often diverge.
There are also other algorithms used for solving the convection-diffusion problems
with pressure linked equations, e.g. SIMPLER4, SIMPLEC5, PISO6and so on.
A.1.1.3 Convergence criteria
The StarCD solver provides the user with two types of output information, [7]. First one
is monitoring values of all the dependent variables at a cell, specified by user. In case
of steady state flows, necessary but not sufficient condition for achieving the solution
is that the variation in the monitoring values should be less than some small tolerance
specified by user. However, this information is not used for control purposes. Another
type of the output information, that is used both for monitoring and control purposes
is the residual. The residual at a particular cell Pand iteration kfor the variable φ,rk
φ,
is defined as the imbalance in the finite volume transport equation (A.5) arising from
incomplete solution, i.e.
rk
φ≡aPφk
P−Xaiφk
i−Su,(A.21)
4SIMPLE Revised
5SIMPLE Consistent
6Pressure Implicit with Splitting of Operators
109
A Numerical Method
Figure A.4: The two-dimensional, backward staggered grid [68].
where the summation is over all neighbouring cells to the specific cell P.
The normalised absolute residual sum, which is actually used to judge convergence,
is defined as follows
Rk
φ≡P|rk
φ|
Mφ
.(A.22)
The summation in the Equation (A.22) is done over all cells in the computational domain.
The scaling factor Mφused to normalise the residual sum is defined:
Mφ=X
p
apφk
p.(A.23)
At each iteration and for each variable, the values Rk
φare computed and printed out. To
decide when to terminate the calculations, the requirement
max(Rk
φ)< λ, (A.24)
stating that all residual have fallen below a tolerance prescribed by the user, must be
satisfied.
110
A Numerical Method
A.1.2 The finite volume method for time dependent flows
Before the Equation (A.4) can be integrated, i.e. discretised with respect to time, an
assumption about the time dependancy of the values φPand φnb with time has to be
made. To calculate the time integral from tto t+ ∆t, the values φat time tor t+ ∆t,
or a combination of these values can be used. This can be generalised by introducing a
weighting parameter αbetween 0 and 1 and the time integral becomes
Z∆t
φdt = [αφ + (1 −α)φ0]∆t, (A.25)
where φ0and φdenote the values at time level tand t+ ∆t, respectively.
Figure A.6: The transient SIMPLE algorithm, [68].
112
A.2 Discretised integrals for the process quantities
Now the Equation (A.4) can be discretised and in general form it becomes
aPφP=X[anb(αφnb + (1 −α)φ0
nb]+[a0
P−X(1 −α)anb]φ0
nb +Su,(A.26)
where the discretisation coefficients satisfy the following relation
aP=α(Xanb −SP) + a0
P.(A.27)
Setting α= 0 the value φP(time level t+∆t) depends on φ0
Pand φ0
nb (time level t),
which gives the explicit discretisation scheme. On the contrary, if αis set equal to 1, the
value φPat time level t+ ∆tdepends on values φnb at the same time level and on φ0
P,
which results in fully implicit discretisation scheme. A combination of these two schemes
is possible, and, e.g. for α= 1/2 we have so called Crank-Nicolson discretisation scheme,
which is implicit, based on central differencing and it is second order accurate in time.
The overall accuracy of the calculation depends, however, also on spatial discretisation,
so the Crank-Nicolson scheme is often combined with spatial central differencing scheme.
The SIMPLE algorithm can be extended to transient calculations. The discretised
time dependent continuity equation (4.17) will have the same form as (A.20), with the
only difference, that the source term a0
I,J is now the imbalance, caused not only by the
incorrect velocity values, but also includes the time change of the fluid dencity ρ.
In transient flow calculations with the implicit discretisation, the iterative procedure
described for steady state calculations such as SIMPLE is applied at each time step until
the convergence criteria is satisfied, Figure A.6.
A.2 Discretised integrals for the process quantities
Once the governing equations are solved numerically, the velocity, pressure and the
temperature field are known. Since these fields are not continuous, but rather discret
fields, the integrals in the Equations (4.30) to (4.33) and (4.46) to (4.51) have also to be
discretised.
Thus the the Equation (4.30) representing the mean fluid temperature at cross-
section at distance xfrom the inlet can be discretised using the volumes of all the cells
that are located between the current position xand x+ ∆x, Figure A.7
Tb(x) = PN
i=1 uiTi∆Vi
PN
i=1 ui∆Vi
,(A.28)
where Nmeans the number of cells placed in the region xand x+∆x, that is not constant
along the channel. The statement above using the cell volume Viis more suitable for
calculations where the staggered grid is not used, than the definition that would use the
cell face surfaces Si.
113
A Numerical Method
Figure A.7: Definitions of the elementary volume and surface.
Similar, the local heat flux at the wall is calculated as follows
˙q(x) = PM
i=1 qi∆ASi
PM
i=1 ∆ASi
,(A.29)
where qirepresents the heat flux magnitude calculated at the boundary cell face ∆ASi
and Mis the number of boundary cells that are placed between xand x+ ∆x.
For the average heat flux the equation of the same form as (A.29) is used, with the
only difference, that Min this case denote the number of boundary cells, that are placed
between x= 0 and x+ ∆x.
All other integrals defined in the section 4.3 are calculated in a similar way and all
the variables defined in the section 4.3 are calculated in a post-processing user subroutine
POSDAT.f that is used by the sofware StarCD and is given in Appendix B.
114
B Definition of the Model for StarCD
C*************************************************************************
PROGRAM GITTER
C*************************************************************************
C DEFINES WHOLE THERMOPLATE MODEL OR A SECTION *
C*************************************************************************
IMPLICIT NONE
INTEGER NCX,NCZ,JJ,PHI(110),PH(110),M,IT
INTEGER I,J,K,A,B,N1,N2,N3,N4,N5,N6,N7,N8,NV,NVX,NVZ,NVR
REAL L,UIN,TIN,TW1,TW2,TT,RHOIN,LX,LZ
DOUBLE PRECISION SL,SL1,ST,ST1,D,PI,VOL,VOL1
DOUBLE PRECISION R,R1,R2,R3,R4,R5,X1,X2,XM,XX,XVD,Y1A,Y1B,Y2
DOUBLE PRECISION F,FF,FF1,F1,F2,F2A,F3,F3A,F3B
DOUBLE PRECISION F4,F41,F42,F43,FE1,FE2,FA1,FA2
DOUBLE PRECISION F4A,F41A,F42A,F43A,DD
DOUBLE PRECISION DX,DZ,XV,YV,ZV
DOUBLE PRECISION XCS(110),YCS(110),ZCS(110)
DOUBLE PRECISION XC(110),YC(110),ZC(110)
C------------------------------------------------------------------------*
C OVERALL SIZE AND GEOMETRY PARAMETERS OF THE MODEL HAVE TO BE SPECIFIED *
C------------------------------------------------------------------------*
PARAMETER (LX=504)
PARAMETER (LZ=36)
PARAMETER (SL=42.)
PARAMETER (SL1=0.)
PARAMETER (ST=36.)
PARAMETER (ST1=0.)
PARAMETER (D=10.)
PARAMETER (R1=D/2)
PARAMETER (R2=15.)
PARAMETER (PI=3.14159)
PARAMETER (R3=12.)
PARAMETER (R4=10.)
PARAMETER (R5=8.1)
PARAMETER (R=8)
115
B Definition of the Model for StarCD
C------------------------------------------------------------------------*
C VALUES NEEDED FOR THE BOUNDARY CONDITIONS DEFINITION *
C------------------------------------------------------------------------*
PARAMETER(L=LX)
PARAMETER(UIN=0.3)
PARAMETER(TIN=293)
PARAMETER(TW1=333)
PARAMETER(TW2=333)
PARAMETER(TT=1.E-4)
PARAMETER(RHOIN=997.561)
PARAMETER(IT=2000)
C CELL NUMBERS IN X AND Z DIRECTION
PARAMETER (NCX=2*LX)
PARAMETER (NCZ=2*LZ)
C CELL SIZE IN X AND Z DIRECTION
PARAMETER (DX=LX/NCX)
PARAMETER (DZ=LZ/NCZ)
C NUMBER OV VERTEX LINES IN X AND Z DIRECTION
PARAMETER (NVX=NCX+1)
PARAMETER (NVZ=NCZ+1)
PARAMETER (NVR=NVX*NVZ)
C------------------------------------------------------------------------*
C INPUT FILES FOR PROSTAR *
C------------------------------------------------------------------------*
OPEN(UNIT=11,FILE=’VERT1.INP’,STATUS=’NEW’)
OPEN(UNIT=12,FILE=’VERT2.INP’,STATUS=’NEW’)
OPEN(UNIT=13,FILE=’VERT3.INP’,STATUS=’NEW’)
OPEN(UNIT=14,FILE=’VERT4.INP’,STATUS=’NEW’)
REWIND 11
REWIND 12
REWIND 13
REWIND 14
C------------------------------------------------------------------------*
C AT EVERY WELDING SPOT A NEV CYLINDRICAL COORDINATE SYSTEM IS DEFINED *
C------------------------------------------------------------------------*
DO J=1,55
XCS(2*J-1)=SL1+(J-1)*SL/2
XCS(2*J)=XCS(2*J-1)
ZCS(2*J-1)=(ST/2)*(1-(-1)**J)
ZCS(2*J)=ZCS(2*J-1)+2*ST
ENDDO
DO J=1,100
YCS(J)=0.
ENDDO
116
DO I=1,100
PHI(I)=-90
ENDDO
JJ=4
DO I=1,95
IF(XCS(I).LE.(LX-SL1).AND.ZCS(I).LE.(LZ-ST1)) THEN
JJ=JJ+1
XC(JJ)=XCS(I)
YC(JJ)=YCS(I)
ZC(JJ)=ZCS(I)
PH(JJ)=PHI(I)
WRITE(12,*) ’LOCAL’,JJ,’ CYLI’,XC(JJ),YC(JJ),ZC(JJ),0,PH(JJ),0
ENDIF
ENDDO
WRITE(12,*) ’CSYS 1’
C------------------------------------------------------------------------*
C VERTEX DEFINITION *
C------------------------------------------------------------------------*
DO K=0,1
DO J=1,NVZ
A=NVR*K+1+NVX*(J-1)
B=NVR*K+NVX*J
DO I=A,B
NV=I
XV=0+DX*(I-A)
ZV=DZ*(J-1)
FF=ST1*(SL1-XV)/SL1+3*ST
FF1=ST1*(SL1+XV-LX)/SL1+3*ST
C------------------------------------------------------------------------*
C FUNCTIONS FOR DESCRIBING OF THE THERMOPLATE SURFACE, EQUATION (5.1) *
C------------------------------------------------------------------------*
F1=D*(1+COS(PI*ZV/ST)*COS(2*PI*(XV-SL1)/SL))/4
C FUNCTIONDEFINING THE FIRST EDGE
F2=D*(4+COS(PI*ZV/ST)+(2-COS(PI*ZV/ST))*COS(PI*XV/SL1))
**SIN(PI*XV/(2*SL1))/8
C FUNCTION DEFINING THE LAST EDGE
F2A=D*(4+COS(PI*ZV/ST)+(2-COS(PI*ZV/ST))*COS(PI*(LX-XV)/SL1))
**SIN(PI*(LX-XV)/(2*SL1))/8
C FUNCTION DEFINING THE RIGHT EDGE
F3A=(1-COS(PI*(3*ST+ST1-ZV)/ST1))*COS(PI*(3*ST+ST1-ZV)/ST1)
**COS(2*PI*(XV-SL1)/SL)
F3B=-3*SIN(PI*(3*ST-5*ST1-ZV)/(2*ST1))
*-SIN(3*PI*(3*ST-ST1-ZV)/(2*ST1))
F3=D*(F3A+F3B)*(1-0.065*(1-COS(2*PI*(3*ST+ST1-ZV)/ST1)))/8
117
B Definition of the Model for StarCD
C FUNCTION DEFINING THE FIRST CORNER
F41=(1-0.1*(1-(COS(4*PI*(3*ST+ST1-ZV)/ST1)))**0.4)
F4=4*F2*F3/(1.22*D)*F41
C FUNCTION DEFINING THE LAST CORNER
F4A=4*F2A*F3/(1.22*D)*F41
C FUNCTION DEFINING THE CONNECTING TUBES
XM=2*SL1*ACOS(SqRT((COS(PI*ZV/ST)-1)/(COS(PI*ZV/ST)-2)))/PI
XX=0.8*(XM-SL1)+SL1
Y1A=R1
Y1B=SqRT((2*R1)**2-(ZV-R2)**2)/2
Y2=D*(4+COS(PI*ZV/ST)+(2-COS(PI*ZV/ST))*COS(PI*XX/SL1))
**SIN(PI*XX/(2*SL1))/8
FE1=XV*(Y2-Y1A)/XX+Y1A
FE2=XV*(Y2-Y1B)/XX+Y1B
FA1=(LX-XV)*(Y2-Y1A)/XX+Y1A
FA2=(LX-XV)*(Y2-Y1B)/XX+Y1B
C------------------------------------------------------------------------*
C THERMOPLATE REGIONS *
C------------------------------------------------------------------------*
IF((XV.GE.SL1).AND.(XV.LE.(LX-SL1))) THEN
IF(ZV.LE.(3*ST)) THEN
F=F1
ENDIF
ENDIF
IF(XV.LT.SL1.AND.ZV.LT.(3*ST+ST1/2).AND.ZV.LT.FF) THEN
F=F2
IF(XV.LT.XX.AND.ZV.LE.R2) THEN
IF(FE1.GT.F2)THEN
F=FE1
ENDIF
ENDIF
IF(XV.LT.XX.AND.ZV.GT.R2.AND.ZV.LE.(R2+2*R1)) THEN
IF(FE2.GT.F2) THEN
F=FE2
ENDIF
ENDIF
ENDIF
IF(XV.GT.(LX-SL1).AND.ZV.LT.(3*ST+ST1/2).AND.
*ZV.LT.FF1) THEN
F=F2A
IF(XV.GT.(LX-XX).AND.ZV.LE.R2) THEN
IF(FA1.GT.F2A)THEN
F=FA1
ENDIF
118
ENDIF
IF(XV.GT.(LX-XX).AND.ZV.GT.R2.AND.ZV.LE.(R2+2*R1)) THEN
IF(FA2.GT.F2A) THEN
F=FA2
ENDIF
ENDIF
ENDIF
IF(XV.GE.SL1/2.AND.XV.LE.(LX-SL1/2).AND.
*ZV.GT.3*ST.AND.ZV.GE.FF.AND.ZV.GE.FF1) THEN
F=F3
ENDIF
IF(XV.LT.SL1/2.AND.ZV.GE.(3*ST+ST1/2)) THEN
F=F4
ENDIF
IF(XV.GT.(LX-SL1/2).AND.ZV.GE.(3*ST+ST1/2)) THEN
F=F4A
ENDIF
YV=(2*K-1)*F
F=0
FF=0
WRITE(11,*) ’V ’,NV,’ ’,XV,’ ’,YV,’ ’,ZV
ENDDO
ENDDO
ENDDO
C------------------------------------------------------------------------*
C CELL DEFINITION *
C------------------------------------------------------------------------*
DO J=1,NVZ-1
A=1+NVX*(J-1)
B=NVX*J
DO I=A,B-1
N1=I
N2=I+1
N3=I+NVX+1
N4=I+NVX
N5=I+NVR
N6=I+NVR+1
N7=I+NVR+NVX+1
N8=I+NVR+NVX
WRITE(11,*) ’C ’,N1,’ ’,N4,’ ’,N3,’ ’,N2,’ ’,N5,’ ’,N8,’ ’,N7,’ ’,N6
ENDDO
ENDDO
C------------------------------------------------------------------------*
C CELLS TO BE REFINED IN X AND Z DIRECTION *
119
B Definition of the Model for StarCD
C------------------------------------------------------------------------*
WRITE(13,*) ’CSET ’,’NEWS ’,’GRAN ’,-R3,R3,-365,365,-10,10,5
DO I=6,JJ
WRITE(13,*) ’CSET ’,’ADD ’,’GRAN ’,-R3,R3,-365,365,-10,10,I
ENDDO
WRITE(13,*) ’CREFINE,2,2,1,CSET,0,0,MERG,NOCOUPLE’
WRITE(13,*) ’CSET,NEWS,GRAN’,-R4,R4,-365,365,-10,10,5
DO I=6,JJ
WRITE(13,*) ’CSET,ADD,GRAN’,-R4,R4,-365,365,-10,10,I
ENDDO
WRITE(13,*) ’CREFINE,2,2,1,CSET,0,0,MERG,NOCOUPLE’
C------------------------------------------------------------------------*
C CELLS TO BE DELETED *
C------------------------------------------------------------------------*
WRITE(13,*) ’CSET ’,’NEWS ’,’GRAN ’,-R5,R5,-365,365,-10,10,5
DO I=6,JJ
WRITE(13,*) ’CSET ’,’ADD ’,’GRAN ’,-R5,R5,-365,365,-10,10,I
ENDDO
WRITE(13,*) ’CDEL CSET’
WRITE(13,*) ’CSET ALL’
WRITE(13,*) ’CPLOT’
WRITE(13,*) ’ CSET ALL’
WRITE(13,*) ’ CPLOT’
C------------------------------------------------------------------------*
C CELLS TO BE REFINED IN Y DIRECTION *
C------------------------------------------------------------------------*
WRITE(13,*) ’ CREFINE,1,1,10,CSET,0,0,MERG,NOCOUPLE’
WRITE(13,*) ’ CSET ALL’
WRITE(13,*) ’ CPLOT’
C------------------------------------------------------------------------*
C COUPLES DEFINITION *
C------------------------------------------------------------------------*
DO I=1,20
X1=-1+(I-1)*50
X2=50+(I-1)*50
IF(X1.LE.LX) THEN
WRITE(13,*) ’ CSET NEWS GRAN’,X1,X2
WRITE(13,*) ’ CPCR CSET 1’
ENDIF
ENDDO
WRITE(13,*) ’CSET ALL’
WRITE(13,*) ’CPLOT’
C------------------------------------------------------------------------*
C CELLS AT THE WALL TO BE REFINED *
120
C------------------------------------------------------------------------*
WRITE(13,*) ’ VSET ALL’
WRITE(13,*) ’ VSET SUBSET GRAN 0.001’,
*LX-0.001,’ -50,50,0.001’,LZ-0.001
WRITE(13,*) ’ VSET SUBSET SURFACE’
WRITE(13,*) ’ CSET NEWS VSET ANY’
WRITE(13,*) ’CPLOT’
WRITE(13,*) ’CREFINE,2,2,3,CSET,0,0,MERG,NOCOUPLE’
WRITE(13,*) ’CSET ALL’
WRITE(13,*) ’CPLOT’
DO I=1,20
X1=-1+(I-1)*50
X2=50+(I-1)*50
IF(X1.LE.LX) THEN
WRITE(13,*) ’ CSET NEWS GRAN’,X1,X2
WRITE(13,*) ’ CPCR CSET 1’
ENDIF
ENDDO
WRITE(13,*) ’CSET ALL’
WRITE(13,*) ’CPLOT’
C------------------------------------------------------------------------*
C BOUNDARY CONDITIONS *
C------------------------------------------------------------------------*
WRITE(14,*) ’ RDEF,1,INLE $ $ $ $ $ $ $ $ $ $’
WRITE(14,*) ’ RNAM,1,INLET’
WRITE(14,*) ’ RDEF,2,OUTL $ $ $ $ $ $ $ $ $ $’
WRITE(14,*) ’ RNAM,2,OUTLET’
WRITE(14,*) ’ RDEF,3,WALL $ $ $ $ $ $ $ $ $ $’
WRITE(14,*) ’ RNAM,3,WALL’
WRITE(14,*) ’ RDEF,4,WALL $ $ $ $ $ $ $ $ $ $’
WRITE(14,*) ’ RNAM,4,WALL’
WRITE(14,*) ’ RDEF,5,SYMP $ $ $ $ $ $ $ $ $ $’
WRITE(14,*) ’ RNAM,5,SYMMETRY’
WRITE(14,*) ’ RDEF,6,SYMP $ $ $ $ $ $ $ $ $ $’
WRITE(14,*) ’ RNAM,6,SYMMETRY’
WRITE(14,*) ’ TEMP,ON’
WRITE(14,*) ’ RADI,OFF’
WRITE(14,*) ’ SOLAR,OFF’
WRITE(14,*) ’ CONJ,OFF’
WRITE(14,*) ’ ACCE,-1,0,0,1,0’
WRITE(14,*) ’ PMAT,1,FLUID,H2O_L,0’
WRITE(14,*) ’ DENS,CONSTANT,997.561’
WRITE(14,*) ’ LVIS,CONSTANT,0.000887’
WRITE(14,*) ’ SPEC,CONSTANT,4181.72’
121
B Definition of the Model for StarCD
WRITE(14,*) ’ COND, CONSTANT,0.620271’
WRITE(14,*) ’ MOLWT,18’
WRITE(14,*) ’ RELAX,0.7,0.3, ,0.95, ,1,1,1, , , , , ,’
WRITE(14,*) ’ SWEEP,10000,10000,10000,100000,,,100,,,’
WRITE(14,*) ’ RESID,0.1,0.1,0.1,0.05, , ,0.1, , , , ,’
WRITE(14,*) ’ WDATA,RESTART,10,0,STAR.PST’
WRITE(14,*) ’ PRFIELD,NONE,,,USER’
WRITE(14,*) ’ ITER’,IT,TT
WRITE(14,*) ’ MEMORY MAXCEL 15000000’
WRITE(14,*) ’ MEMORY MAXVRT 15000000’
WRITE(14,*) ’ MEMORY MAXCUT 15000000’
WRITE(14,*) ’ CSET,ALL’
WRITE(14,*) ’ VSET,ALL’
WRITE(14,*) ’ VMER,VSET,,0.0001,,,,,,,LOW,DELETE,C’
WRITE(14,*) ’ VCOMPRESS,ALL $Y’
WRITE(14,*) ’ CCOMPRESS $Y’
WRITE(14,*) ’ CSET,ALL’
WRITE(14,*) ’ VSET,ALL’
WRITE(14,*) ’ CPLOT’
WRITE(14,*) ’ BSET ALL’
WRITE(14,*) ’ BDEL BSET’
WRITE(14,*) ’ RDEF,1,INLET,STANDARD’
WRITE(14,*) UIN,0,0,1,0,TIN,RHOIN
WRITE(14,*) ’ RTURB,1,NONE’
WRITE(14,*) ’ RINLET,1,MASS,N’
WRITE(14,*) ’ RDEF,3,WALL,STANDARD’
WRITE(14,*) ’ NOSL,STAND,9,’
WRITE(14,*) ’ 0,0,0,1,0’
WRITE(14,*) ’ FIXED’,TW1,0
WRITE(14,*) ’ RDEF,4,WALL,STANDARD’
WRITE(14,*) ’ NOSL,STAND,9,’
WRITE(14,*) ’ 0,0,0,1,0’
WRITE(14,*) ’ FIXED’,TW2,0
C------------------------------------------------------------------------*
C BOUNDARIES LOCATION *
C------------------------------------------------------------------------*
C INLET
WRITE(14,*) ’ CSET,ALL’
WRITE(14,*) ’ CSET, NEWS, GRAN -2, 2,-10,10’,-10,R2+2*R1+0.01
WRITE(14,*) ’ CPLOT’
122
WRITE(14,*) ’ ZOOM OFF $REPLOT’
WRITE(14,*) ’ VIEW,-1,0,0 $REPLOT’
WRITE(14,*) ’ BZON,1,,’
WRITE(14,*) ’ BSET ALL’
WRITE(14,*) ’ BSET SUBSET GRAN’,-2*L,0.001
WRITE(14,*) ’ BSET INVERT’
WRITE(14,*) ’ BSET SUBSET REGION 1’
WRITE(14,*) ’ BDEL BSET’
C OUTLET
WRITE(14,*) ’ CSET NEWS GRAN’,L-2,2*L,-10,10,-10,R2+2*R1+0.01
WRITE(14,*) ’ CPLOT’
WRITE(14,*) ’ VIEW,1,0,0,$REPLOT’
WRITE(14,*) ’ ZOOM OFF $REPLOT’
WRITE(14,*) ’ BZON,2,,’
WRITE(14,*) ’ BSET ALL’
WRITE(14,*) ’ BSET SUBSET GRAN’,L-0.001,2*L
WRITE(14,*) ’ BSET INVERT’
WRITE(14,*) ’ BSET SUBSET REGION 2’
WRITE(14,*) ’ BDEL BSET’
C SYMPLANE LEFT
WRITE(14,*) ’ CSET ALL’
WRITE(14,*) ’ CSET NEWS GRAN’,-2*L,2*L,-50,50,-50,2
WRITE(14,*) ’ CPLOT’
WRITE(14,*) ’ VIEW,0,0,-1’
WRITE(14,*) ’ ZOOM OFF $REPLOT’
WRITE(14,*) ’ BZON,5,,’
WRITE(14,*) ’ BSET ALL’
WRITE(14,*) ’ BSET SUBSET GRAN’,-2*L,2*L,-50,50,-50,0.001
WRITE(14,*) ’ BSET INVERT’
WRITE(14,*) ’ BSET,SUBSET, REGION, 5’
WRITE(14,*) ’ BDEL BSET’
C SYMPLANE RIGHT
WRITE(14,*) ’ CSET ALL’
WRITE(14,*) ’ CSET NEWS GRAN’,-2*L,2*L,-50,50,ST-2,2*LZ
WRITE(14,*) ’ CPLOT’
WRITE(14,*) ’ VIEW,0,0,1’
WRITE(14,*) ’ ZOOM OFF $REPLOT’
WRITE(14,*) ’ BZON,6,,’
WRITE(14,*) ’ BSET ALL’
WRITE(14,*) ’ BSET SUBSET GRAN’,-2*L,2*L,-50,50,LZ-0.001,2*LZ
WRITE(14,*) ’ BSET INVERT’
WRITE(14,*) ’ BSET,SUBSET, REGION, 6’
WRITE(14,*) ’ BDEL BSET’
C LOWER WALL W1
123
B Definition of the Model for StarCD
WRITE(14,*) ’ CSET ALL’
WRITE(14,*) ’ ZOOM OFF’
WRITE(14,*) ’ CPLOT’
WRITE(14,*) ’ VIEW,0,-1,0,$REPLOT’
WRITE(14,*) ’ BZON,3,,’
WRITE(14,*) ’ BSET ALL’
WRITE(14,*) ’ BSET NEWS GRAN’,-2*LX,0.001
WRITE(14,*) ’ BSET ADD GRAN’, LX-0.001,2*LX
WRITE(14,*) ’ BSET ADD GRAN’, -2*LX,2*LX,-50,50,-LZ,0.001
WRITE(14,*) ’ BSET ADD GRAN’, -2*LX,2*LX,-50,50,LZ-0.001,2*LZ
WRITE(14,*) ’ BSET SUBSET REGION 3’
WRITE(14,*) ’ BDEL BSET’
DO I=5,JJ
WRITE(14,*) ’ VSET NEWS GRAN’,-R-0.001,R+0.001,-365,365,-50,50,I
WRITE(14,*) ’ BSET ADD VSET ALL’
ENDDO
WRITE(14,*) ’BMOD,BSET,7’
C UPPER WALL W2
WRITE(14,*) ’ CSET ALL’
WRITE(14,*) ’ VIEW,0,1,0,$REPLOT’
WRITE(14,*) ’ BZON,4,,’
WRITE(14,*) ’ BSET ALL’
WRITE(14,*) ’ BSET NEWS GRAN’,-2*LX,0.001
WRITE(14,*) ’ BSET ADD GRAN’, LX-0.001,2*LX
WRITE(14,*) ’ BSET ADD GRAN’, -2*LX,2*LX,-50,50,-LZ,0.001
WRITE(14,*) ’ BSET ADD GRAN’, -2*LX,2*LX,-50,50,LZ-0.001,2*LZ
WRITE(14,*) ’ BSET SUBSET REGION 4’
WRITE(14,*) ’ BDEL BSET’
CLOSE(11)
CLOSE(12)
CLOSE(13)
CLOSE(14)
END
124
C POSDAT.f
C*************************************************************************
SUBROUTINE POSDAT(KEY,VOL,U,TE,ED,T,P,VIST,DEN,CP,VISM,CON,
* F,ICLMAP,ICTID,RESOR,VF,FORCB,IRN,PREFM,LEVEL)
C POSDAT.f is a StarCD post-processing procedure
C*************************************************************************
C------------------------------------------------------------------------*
C STAR VERSION 3.24.000 *
C------------------------------------------------------------------------*
INCLUDE ’comdb.inc’
COMMON/USR001/INTFLG(100)
DIMENSION KEY(-NBMAXU:NCTMXU),VOL(NCTMXU),U(3,-NBMAXU:NCMAXU),
* TE(-NBMAXU:NCMAXU),ED(-NBMAXU:NCMAXU),T(-NBMAXU:NCTMXU,1+NSCU),
* P(-NBMAXU:NCMAXU),VIST(-NBMAXU:NCMAXU),DEN(-NBMAXU:NCTMXU),
* CP(-NBMAXU:NCTMXU),VISM(-NBMXVU:NCMXVU),CON(-NBMXCU:NCMXCU),
* F(3,-NBMAXU:NCMAXU),ICLMAP(NCTMXU),ICTID(NCTMXU),
* RESOR(89,-100:100),VF(NCDMXU),FORCB(3,NWLMX),IRN(NWLMX)
C------------------------------------------------------------------------*
C VARIABLES DEFINITION *
C------------------------------------------------------------------------*
DOUBLE PRECISION P1,P2,P3,P4,P5,P6,P7,P8,PIN,PTIN,PT
DOUBLE PRECISION T1,T2,T3,T4,T5,TT
DOUBLE PRECISION TW1,TW2,TWQ1,TWQ2,AREA1,AREA2
DOUBLE PRECISION HF1,HF2,HFQ1,HFQ2
DOUBLE PRECISION HF11,HF21,HFQ11,HFQ21
DOUBLE PRECISION HCF1,HCF2,HCFQ1,HCFQ2
DOUBLE PRECISION HCFQ11,HCFQ21
DOUBLE PRECISION TM1(2000,12),TM2(2000,14)
DOUBLE PRECISION A1,A2
DOUBLE PRECISION X,DX,G61,G62,G71,G72,G8,VV
DOUBLE PRECISION G5A,G5B
DOUBLE PRECISION Z,DZ,Z0,Z1,XMAX,AA,BB,A,B
INTEGER L,NPP,WR1,WR2,IMAX,IMAX1,N,N1,NPP1
INTEGER L1(5000)
INCLUDE ’usrdat.inc’
C DX is the cell size in X direction in mm
PARAMETER(DX=0.0005)
125
C POSDAT.f
C Residual tolerance
PARAMETER(TT=1.E-4)
C Number of the cross sections, steps in x direction
PARAMETER(NPP=1008)
C------------------------------------------------------------------------*
C ADDITIONAL COMMON BLOCKS *
C------------------------------------------------------------------------*
POINTER(P_IRTYP,IRTYP(0:NREGMX))
COMMON/C06004/P_IRTYP
POINTER(P_LX,LX(4))
COMMON /C06013/P_LX
POINTER(P_INDWL,INDWL(4))
COMMON/C06014/P_INDWL
POINTER(P_IBLMAP,IBLMAP(NBMAX))
COMMON/C1C003/P_IBLMAP
POINTER(P_IDC,IDC(4))
COMMON/C33001/P_IDC
POINTER(P_DRM,DRM(4))
COMMON/C33004/P_DRM
POINTER(P_LQ,LQ(6,NCTMAX))
COMMON /DCO5/P_LQ
POINTER(P_S,S(3,3,-NBMAX:NCTMAX))
COMMON /DCO14/P_S
POINTER(P_CX,CX(3,-NBMAX:NCTMAX))
COMMON/DCO21/P_CX
POINTER(P_HCB,HCB(NWSMX,1+NSC))
POINTER(P_HFB,HFB(NWSMX,1+NSC))
COMMON /DCO19A/P_HCB
COMMON /DCO19B/P_HFB
POINTER(P_SB,SB(3,NBN1MX:NCUTMX+1))
COMMON/DCO31R2/P_SB
COMMON /DROP1A/ NDR
COMMON /DROP3/ TWPHL
COMMON /INDMP1/ NDIR(-3:3),NDINAB(6)
COMMON /INDMP2/ NDIN(6)
COMMON /MPP430/ MASTER,SLAVE,HPCRUN,PARRUN
COMMON /MPP1/ NOPROC,MYID
COMMON /SLIPW/ NWSBT
LOGICAL MASTER,SLAVE,HPCRUN,PARRUN,TWPHL
C------------------------------------------------------------------------*
C OUTPUT FILES DEFINITION *
C------------------------------------------------------------------------*
IF(INTFLG(1).EQ.0) THEN
OPEN(85,file=’../results1.dat’)
126
OPEN(86,file=’../results2.dat’)
INTFLG(1) = 1
ENDIF
C------------------------------------------------------------------------*
C THE EVALUATION TAKE PLACE AT AFTER THE LAST ITERATION *
C------------------------------------------------------------------------*
IF(LEVEL.EQ.2) THEN
IMAX=ITERS+ITSTEP-1
IF((RESOR(1,1).LT.TT.AND.RESOR(2,1).LT.TT.AND.
*RESOR(3,1).LT.TT.AND.RESOR(4,1).LT.TT.AND.
*RESOR(5,1).LT.TT.AND.RESOR(6,1).LT.TT).OR.
*(ITER.EQ.IMAX)) THEN
C------------------------------------------------------------------------*
C INITIALIZING *
C------------------------------------------------------------------------*
T1=0.
T2=0.
T3=0.
T4=0.
T5=0.
P1=0.
P2=0.
P3=0.
P4=0.
P5=0.
P6=0.
P7=0.
PI=0.
PT=0.
PIN=0.
PTIN=0.
A=0.
B=DX
X=DX*1000/2
HF1=0.
HF2=0.
HF11=0.
HF21=0.
HFQ1=0.
HFQ2=0.
HFQ11=0.
HFQ21=0.
HCF1=0.
HCF2=0.
127
C POSDAT.f
HCFQ1=0.
HCFQ2=0.
HCFQ11=0.
HCFQ21=0.
POV1=0.
POV2=0.
POVQ1=0.
POVQ2=0.
DO 100 K=1,NPP
DO 110 J=1,6
TM(K,J)=0.
110 CONTINUE
100 CONTINUE
C------------------------------------------------------------------------*
C THE EVALUATION TAKES PLACE AT EACH STEP X *
C------------------------------------------------------------------------*
DO 200 K=1,NPP
X=DX/2+(K-1)*DX
A=(K-1)*DX
B=K*DX
C------------------------------------------------------------------------*
C VARIABLES EVALUATED AT CELL CENTRES *
C------------------------------------------------------------------------*
DO 210 I=1,NCELL
IF(CX(1,I).GE.A) THEN
IF(CX(1,I).LE.B) THEN
C VV is the velocity magnitude for the cell I
VV=SQRT(U(1,I)**2+U(2,I)**2+U(3,I)**2)
C is the total pressure for the cell I
PT=P(I)+DEN(I)*(VV**2)/2
C T1 is the volume integral rho*cp*u*T*dV
T1=T1+DEN(I)*CP(I)*U(1,I)*T(I,1)*VOL(I)
C T2 is the volume integral ro*cp*u*dV
T2=T2+DEN(I)*CP(I)*U(1,I)*VOL(I)
C T3 is the volume integral u*dV, volume flow rate*dx
T3=T3+U(1,I)*VOL(I)
C T4 is rho*u*dV, mass flow rate*dx
T4=T4+U(1,I)*DEN(I)*VOL(I)
C P1 is the volume integral p*dV
P1=P1+P(I)*VOL(I)
C P2 is the volume integral PT*dV
P2=P2+PT*VOL(I)
c P3 is the cross-section area
P3=P3+VOL(I)
128
C P4 is the mean density*dV
P4=P4+DEN(I)*VOL(I)
ENDIF
ENDIF
210 CONTINUE
C------------------------------------------------------------------------*
C VARIABLES EVALUATED AT WALL BOUNDARIES, THERMOPLATE SURFACE *
C------------------------------------------------------------------------*
DO 220 J=1,NWSBT
IB=INDWL(J)
IC=LX(IB)/7
IP=ICLMAP(IC)
INDX=MOD(LX(IB),7)
IBG=LQ(INDX,IC)
IBP=IBLMAP(-IBG)
IREG=IRN(IB)
IBTYP=IRTYP(IREG)
IF(CX(1,IBG).GE.A.AND.CX(1,IBG).LE.B) THEN
IF(ABS(HFB(IB,1)).GT.0.000001) THEN
C------------------------------------------------------------------------*
C LOWER WALL, W1 *
C------------------------------------------------------------------------*
IF(CX(2,IBG).LT.0) THEN
A1=SQRT((S(3,1,IC))**2+(S(3,2,IC))**2+(S(3,3,IC))**2)
C HF is the heat flux*surface, A<X<B
HF1=HF1+HFB(IB,1)*A1
C HFQ is the mean heat flux*surface
HFQ1=HFQ1+HFB(IB,1)*A1
TWQ1=TWQ1+T(IBG,1)*A1
C POV is the wall surface A<x<B
POV1=POV1+A1
C POVQ is the wall surface 0<x<B
POVQ1=POVQ1+A1
ENDIF
C------------------------------------------------------------------------*
C UPPER WALL, W2 *
C------------------------------------------------------------------------*
IF(CX(2,IBG).GT.0) THEN
A2=SQRT((S(3,1,IC))**2+(S(3,2,IC))**2+(S(3,3,IC))**2)
C HF is the heat flux*surface, A<X<B
HF2=HF2+HFB(IB,1)*A2
C HFQ is the mean heat flux*surface
HFQ2=HFQ2+HFB(IB,1)*A2
TWQ2=TWQ2+T(IBG,1)*A2
129
C POSDAT.f
C POV is the wall surface A<x<B
POV2=POV2+A2
C POVQ is the wall surface 0<x<B
POVQ1=POVQ1+A1
ENDIF
ENDIF
ENDIF
220 CONTINUE
T5=T1/T2
HF11=HF1/POV1
HF21=HF2/POV2
HFQ11=HFQ1/POVQ1
HFQ21=HFQ2/POVQ2
HCF1=HF11/(TW1-T5)
HCF2=HF21/(TW2-T5)
C------------------------------------------------------------------------*
C Introduced for parallel computing (this was needed for StarCD *
C version 3.24) *
if (nhpc.gt.1) then *
t1_sum = dgsum(t1) *
t2_sum = dgsum(t2) *
t3_sum = dgsum(t3) *
t4_sum = dgsum(t4) *
p1_sum = dgsum(p1) *
p2_sum = dgsum(p2) *
p3_sum = dgsum(p3) *
p4_sum = dgsum(p4) *
hf1_sum = dgsum(hf1) *
hf2_sum = dgsum(hf2) *
hfq1_sum = dgsum(hfq1) *
hfq2_sum = dgsum(hfq2) *
pov1_sum = dgsum(pov1) *
pov2_sum = dgsum(pov2) *
povq1_sum = dgsum(povq1) *
povq2_sum = dgsum(povq2) *
twq1_sum = dgsum(twq1) *
twq2_sum = dgsum(twq2) *
if (ihpc.ge.1) then *
t3 = t3_sum/dx *
t4 = t4_sum/dx *
t5 = t1_sum/t2_sum *
p5 = p1_sum/p3_sum *
p6 = p2_sum/p3_sum *
p7 = p4_sum/p3_sum *
130
hf11 = hf1_sum/pov1_sum *
hf21 = hf2_sum/pov2_sum *
hfq11 = hfq1_sum/povq1_sum *
hfq21 = hfq2_sum/povq2_sum *
tw1 = twq1_sum/pov1_sum *
tw2 = twq2_sum/pov2_sum *
hcf1 = hf11/(tw1-t5) *
hcf2 = hf21/(tw2-t5) *
hcfq1 = hcfq1+hcf1*dx *
hcfq2 = hcfq2+hcf2*dx *
hcfq11= hcfq1/b *
hcfq21= hcfq2/b *
IF(K.EQ.1) THEN *
pin = p1_sum/p3_sum *
ptin = p2_sum/p3_sum *
ENDIF *
g5a = t3*(pin-p5)/dx *
g5b = t3*(ptin-p6)/dx *
g61 = hfq1_sum *
g62 = hfq2_sum *
g71 = povq1_sum *
g72 = povq2_sum *
p8 = p8+p3_sum *
endif *
endif *
C ----- Introduced for parallel computing *
C------------------------------------------------------------------------*
TM1(K,1)=X*1000
TM1(K,2)=HF11
TM1(K,3)=HF21
TM1(K,4)=HFQ11
TM1(K,5)=HFQ21
TM1(K,6)=HCF1
TM1(K,7)=HCF2
TM1(K,8)=HCFQ11
TM1(K,9)=HCFQ21
TM1(K,10)=T5
TM1(K,11)=TW1
TM1(K,12)=TW2
TM2(K,1)=X*1000
TM2(K,2)=P5
TM2(K,3)=P6
TM2(K,4)=T3
TM2(K,5)=T4
131
C POSDAT.f
TM2(K,6)=G61
TM2(K,7)=G62
TM2(K,8)=G71
TM2(K,9)=G72
TM2(K,10)=P8
TM2(K,11)=G5A
TM2(K,12)=G5B
TM2(K,13)=(G61+G62)/G5A
TM2(K,14)=(G61+G62)/G5B
C------------------------------------------------------------------------*
C Data are written out at each step x *
C------------------------------------------------------------------------*
IF(K.eq.1) THEN
WRITE(85,*)’ x q1 q2 qm1 qm2 h1 h2 hm1 hm2 tb tw1 tw2’
WRITE(86,*)’ x pm ptm VS MSv Q1 Q2 S1 S2 V P Pt Q/P Q/Pt’
ENDIF
WRITE(85,*) (TM1(K,J),J=1,12)
WRITE(86,*) (TM2(K,J),J=1,14)
C------------------------------------------------------------------------*
C Some values have to be restored *
C------------------------------------------------------------------------*
T1=0.
T2=0.
T3=0.
T4=0.
T5=0.
P1=0.
P2=0.
P3=0.
P4=0.
P5=0.
P6=0.
P7=0.
PI=0.
PT=0.
HF1=0.
HF2=0.
HF11=0.
HF21=0.
HCF1=0.
HCF2=0.
HCF11=0.
HCF21=0.
POV1=0.
132
POV2=0.
TW1=0.
TW2=0.
TWQ1=0.
TWQ2=0.
t1_sum = 0
t2_sum = 0
t3_sum = 0
t4_sum = 0
p1_sum = 0
p2_sum = 0
p3_sum = 0
p4_sum = 0
hf1_sum = 0
hf2_sum = 0
pov1_sum = 0
pov2_sum = 0
twq1_sum = 0
twq2_sum = 0
200 CONTINUE
ENDIF
ENDIF
RETURN
END
133
C POSDAT.f
134
D Mesh Test and Other Numerical
Details
The complex geometry of the flow channel contains several areas with large gradients
of interesting quantities, which requires careful spatial resolutions. This is particularly
the case near the walls and in the separation and reattachment regions of the recir-
culation zones. This demands for an appropriate and expensive discretisation of the
computational domain.
To test the discretisation effect on heat transfer, a relatively coarse grid with equal
number of points in the y-direction at each x, z-position have been used, Figure 4.1,
and subsequently refined by increasing the number of grid points locally, or in the whole
domain. Figure D.1 illustrates the case of a grid refinement in the wall region. Figure D.2
shows the effect of grid refinement on the heat transfer coefficient for a larger Reynolds
number. For smaller Reynolds numbers, the deviations were still observed, but they
were considerable smaller. The curve a) in Figure D.2 lies considerably below the other
ones. This is due chiefly to the coarseness of the grid, which results in relatively large
errors, particularly in the wall region. The scatter of the curves b) to c) is much smaller
and is rooted, on one side, in the physical uncertainty due to the grid size and in the
numerical residua, on the other side. The reduction of the physical uncertainty by a
refinement of the grid is expected to result in a correspondingly larger error due to the
numerical residua, and an optimum regarding the reliability of the numerical results
is expected when these uncertainties are the same, or at least of the same order of
magnitude. A criterion for and objective judgement of the resulting uncertainty can
scarcely be formulated. For this reason, and with regard to the flood of the data to be
processed, the discretisation a) has been adopted in the numerical simulations.
The algorithm used was SIMPLE, a scalar solver type and conjugate gradient
method (see e.g. StarCD Methodology Handbook [7] and Versteeg and Malalasekera
[68]), at a maximal residual tolerance λ= 10−4.
135
D Mesh Test and Other Numerical Details
Figure D.1: Illustration of the grid used in the numerical simulations.
Figure D.2: Effect of mesh size on the numerical results: a) equidistant mesh, refined
at wall (actually used, Figure D.1), b) like a), but the number of cells doubled in all
three directions, c) like b), but the number of cells doubled in all three directions.
136
E Constants
Table E.1: Constants used in Equation (5.3) with I=J=K=L= 2 and M=m= 0.
c0,0,0,0,0−371.939 c1,0,0,0,04779.99 c2,0,0,0,0−14562.2
c0,0,0,1,02.7805 c1,0,0,1,0−32.5245 c2,0,0,1,090.3223
c0,0,0,2,0−0.001084 c1,0,0,2,00.0129844 c2,0,0,2,0−0.03698
c0,0,1,0,051611.7c1,0,1,0,0−645786 c2,0,1,0,02.018 ·106
c0,0,1,1,0−552.129 c1,0,1,1,06433.27 c2,0,1,1,0−17765.4
c0,0,1,2,00.22232 c1,0,1,2,0−2.65846 c2,0,1,2,07.55836
c0,0,2,0,0−1.316 ·106c1,0,2,0,02.017 ·107c2,0,2,0,0−7.1599 ·107
c0,0,2,1,025898.1c1,0,2,1,0−305430 c2,0,2,1,0850628
c0,0,2,2,0−10.026 c1,0,2,2,0120.587 c2,0,2,2,0−343.563
c0,1,0,0,012213.6c1,1,0,0,0−161076 c2,1,0,0,0512002
c0,1,0,1,0−99.9053 c1,1,0,1,01181.18 c2,1,0,1,0−3304.62
c0,1,0,2,00.03971 c1,1,0,2,0−0.480533 c2,1,0,2,01.37835
c0,1,1,0,0−1.79387 ·106c1,1,1,0,02.385 ·107c2,1,1,0,0−7.84 ·107
c0,1,1,1,020357.9c1,1,1,1,0−239478 c2,1,1,1,0666282
c0,1,1,2,0−8.23688 c1,1,1,2,099.6535 c2,1,1,2,0−285.594
c0,1,2,0,05.80779 ·107c1,1,2,0,0−9.229 ·108c2,1,2,0,03.326 ·109
c0,1,2,1,0−963007 c1,1,2,1,01.151 ·107c2,1,2,1,0−3.234 ·107
c0,1,2,2,0372.797 c1,1,2,2,0−4543.44 c2,1,2,2,013053.7
c0,2,0,0,0−121519 c1,2,0,0,01.675 ·106c2,2,0,0,0−5.502 ·106
c0,2,0,1,0891.516 c1,2,0,1,0−10735.7c2,2,0,1,030432.7
c0,2,0,2,0−0.35433 c1,2,0,2,04.35689 c2,2,0,2,0−12.6416
c0,2,1,0,02.0514 ·107c1,2,1,0,0−2.898 ·108c2,2,1,0,09.797 ·108
c0,2,1,1,0−185038 c1,2,1,1,02.223 ·106c2,2,1,1,0−6.279 ·106
c0,2,1,2,074.1139 c1,2,1,2,0−912.556 c2,2,1,2,02647.22
c0,2,2,0,0−7.988 ·108c1,2,2,0,01.27 ·1010 c2,2,2,0,0−4.527 ·1010
c0,2,2,1,08.7854 ·106c1,2,2,1,0−1.072 ·108c2,2,2,1,03.0568 ·108
c0,2,2,2,0−3358.9c1,2,2,2,041659.2c2,2,2,2,0−121107
137
E Constants
138
F Numerical Results
Table F.1: Simulation results for model parameters given in Tables 4.1 and 4.2.
sLsTδ Re Nu ∆p sLsTδ L Nu ∆p
[mm] [mm] [mm] [-] [-] [Pa] [mm] [mm] [mm] [mm] [-] [Pa]
28 36 3.4 50 4.45 67.4 63 25 4.2 500 10.53 258.6
28 36 3.4 200 5.58 335.8 63 25 4.2 1000 18.36 897.
28 36 3.4 500 8.95 1038.4 63 25 4.2 2000 30.53 2887.6
28 36 3.4 750 11.46 1797.0 63 25 4.2 3000 39.20 6200.7
28 36 3.4 1000 13.88 2755.4 63 25 6.3 50 5.78 4.4
28 36 3.4 1500 17.91 4961.2 63 25 6.3 500 14.03 102.8
28 36 3.4 2000 21.68 7725.4 63 25 6.3 1000 24.66 340.4
28 36 3.4 2500 25.49 10985.4 63 25 6.3 2000 38.48 1163.7
28 36 3.4 3000 29.27 14867.3 63 25 6.3 3000 49.16 2603.8
28 36 3.4 3500 32.66 19176.7 63 25 8.4 50 6.34 2.0
28 36 3.4 3800 36.23 22488.9 63 25 8.4 500 17.72 54.1
28 36 3.4 4000 36.76 24154.3 63 25 8.4 1000 29.14 176.0
35 30 3.4 50 4.96 49.7 63 25 8.4 2000 45.21 645.6
35 30 3.4 1000 15.95 2528.6 63 25 8.4 3000 57.35 1420.9
35 30 3.4 2000 34.67 2228.3 63 30 4.2 50 5.05 15.9
35 30 6.3 50 5.18 9.3 63 30 4.2 500 9.80 331.9
35 30 6.3 1000 24.11 557.3 63 30 4.2 1000 18.31 1100.9
35 30 6.3 2000 37.70 1679.4 63 30 4.2 2000 29.39 3965.3
35 35 3.4 50 4.42 55.1 63 30 4.2 3000 35.56 7796.5
35 35 3.4 1000 14.29 2495.3 63 30 6.3 50 5.58 5.0
35 35 3.4 2000 21.84 6940.1 63 30 6.3 500 14.08 128.6
35 35 6.3 50 4.71 10.2 63 30 6.3 1000 24.81 471.3
35 35 6.3 1000 21.36 534.3 63 30 6.3 2000 36.09 1559.1
35 35 6.3 2000 33.41 1587.4 63 30 6.3 3000 45.17 3189.1
35 36 3.4 50 4.18 57.2 63 30 8.4 50 6.07 2.2
35 36 3.4 200 5.33 295.9 63 30 8.4 500 18.30 69.0
35 36 3.4 500 9.24 970.9 63 30 8.4 1000 29.0 245.2
35 36 3.4 750 11.84 1683.3 63 30 8.4 2000 42.26 818.8
35 36 3.4 1000 14.34 2543.3 63 30 8.4 3000 52.20 1675.9
35 36 3.4 1500 18.63 4688.7 63 35 4.2 50 4.88 17.8
35 36 3.4 2000 21.98 7195.4 63 35 4.2 500 9.61 383.1
35 36 3.4 2500 25.69 10057.5 63 35 4.2 1000 16.99 1187.2
35 36 3.4 3000 28.40 13367.5 63 35 4.2 2000 26.34 3817.6
35 36 3.4 3500 31.39 17135.9 63 35 4.2 3000 32.06 7319.6
35 36 3.4 3800 33.71 18963.0 63 35 6.3 50 5.33 5.6
35 36 3.4 4000 34.00 21510.6 63 35 6.3 500 13.49 145.2
139
F Numerical Results
sLsTδ Re Nu ∆p sLsTδ L Nu ∆p
[mm] [mm] [mm] [-] [-] [Pa] [mm] [mm] [mm] [mm] [-] [Pa]
35 40 3.4 50 4.16 59.3 63 35 6.3 1000 22.85 474.2
35 40 3.4 1000 12.76 2467.6 63 35 6.3 2000 32.69 1475.4
35 40 3.4 2000 19.56 6281.6 63 35 6.3 3000 40.26 2832.0
35 40 6.3 50 4.38 10.6 63 35 8.4 50 5.76 2.5
35 40 6.3 1000 18.97 521.9 63 35 8.4 500 17.38 75.2
35 40 6.3 2000 29.91 1546.9 63 35 8.4 1000 26.95 245.2
40 30 3.4 50 2.74 28.8 63 35 8.4 2000 38.92 764.3
40 30 3.4 1000 6.23 1398.9 63 35 8.4 3000 48.32 1488.5
40 30 3.4 2000 7.52 4009.6 63 36 3.4 50 4.61 34.0
40 30 6.3 50 5.40 8.1 63 36 3.4 100 5.14 73.9
40 30 6.3 1000 24.01 544.2 63 36 3.4 200 6.29 184.1
40 30 6.3 2000 37.17 1644.2 63 36 3.4 300 6.97 371.5
40 35 3.4 50 4.57 48.8 63 36 3.4 400 7.67 519.7
40 35 3.4 1000 14.45 2433.7 63 36 3.4 500 8.43 729.4
40 35 3.4 2000 21.96 6795.9 63 36 3.4 750 11.82 1314.4
40 35 6.3 50 4.91 9.0 63 36 3.4 1000 14.29 2016.0
40 35 6.3 1000 21.47 523.2 63 36 3.4 1500 19.19 3873.3
40 35 6.3 2000 33.31 1550.5 63 36 3.4 2000 23.29 6219.0
40 40 3.4 50 4.31 53.4 63 36 3.4 2500 26.55 8902.7
40 40 3.4 1000 13.06 2409.2 63 36 3.4 3000 29.22 11867.2
40 40 3.4 2000 19.84 6667.7 63 36 3.4 3500 31.68 15203.2
40 40 6.3 50 4.54 9.7 63 36 3.4 3800 33.19 17381.8
40 40 6.3 1000 19.32 510.7 63 36 3.4 4000 34.03 18808.0
40 40 6.3 2000 29.98 1500.8 70 25 4.2 50 5.06 13.0
42 20 3.4 50 5.40 30.7 70 25 4.2 500 10.89 214.9
42 20 3.4 200 8.03 172.0 70 25 4.2 1000 15.9 675.8
42 20 3.4 500 12.02 699.3 70 25 4.2 2000 28.85 2049.3
42 20 3.4 750 17.60 1564.1 70 25 4.2 3000 37.42 4059.3
42 20 3.4 1000 21.95 2672.5 70 25 6.3 50 5.61 4.1
42 20 3.4 1500 29.41 5577.7 70 25 6.3 500 12.59 77.7
42 20 3.4 2000 35.45 9617.3 70 25 6.3 1000 22.25 258.4
42 20 3.4 2500 38.86 13485.5 70 25 6.3 2000 36.62 800.1
42 20 3.4 3000 42.25 18045.5 70 25 6.3 3000 47.14 1715.6
42 20 3.4 3500 46.17 23505.6 70 25 8.4 50 6.13 1.8
42 20 3.4 4000 49.13 29796.0 70 25 8.4 500 15.95 41.5
42 25 3.4 50 5.08 35.9 70 25 8.4 1000 27.12 129.4
42 25 3.4 200 7.03 223.4 70 25 8.4 2000 42.66 421.6
42 25 3.4 500 11.66 827.9 70 25 8.4 3000 54.33 895.8
42 25 3.4 750 15.59 1551.17 70 30 4.2 50 4.94 14.5
42 25 3.4 1000 19.38 2498.1 70 30 4.2 500 9.55 274.1
42 25 3.4 1500 23.55 4661.7 70 30 4.2 1000 16.58 896.6
42 25 3.4 2000 28.12 7377.0 70 30 4.2 2000 30.29 3195.0
42 25 3.4 2500 33.45 11008.2 70 30 4.2 3000 36.85 6574.6
42 25 3.4 3000 35.72 14256.4 70 30 6.3 50 5.44 4.5
42 25 3.4 3500 39.12 18622.0 70 30 6.3 500 12.77 105.6
42 25 3.4 4000 42.38 23465.0 70 30 6.3 1000 23.60 366.4
42 30 3.4 50 4.81 41.5 70 30 6.3 2000 36.36 1337.8
42 30 3.4 200 5.94 257.3 70 30 6.3 3000 45.64 2816.4
42 30 3.4 500 10.54 855.8 70 30 8.4 50 5.86 1.9
140
sLsTδ Re Nu ∆p sLsTδ L Nu ∆p
[mm] [mm] [mm] [-] [-] [Pa] [mm] [mm] [mm] [mm] [-] [Pa]
42 30 3.4 750 14.01 1565.5 70 30 8.4 500 15.97 51.2
42 30 3.4 1000 17.06 2439.3 70 30 8.4 1000 27.08 173.0
42 30 3.4 1500 20.47 4420.7 70 30 8.4 2000 41.23 635.8
42 30 3.4 2000 24.43 6953.1 70 30 8.4 3000 50.71 1388.0
42 30 3.4 2500 28.06 9851.9 70 35 4.2 50 4.81 16.2
42 30 3.4 3000 31.44 13033.4 70 35 4.2 500 8.87 328.0
42 30 3.4 3500 34.48 16748.9 70 35 4.2 1000 16.04 1037.3
42 30 3.4 4000 37.42 21357.0 70 35 4.2 2000 27.19 3719.1
42 36 2.5 50 4.37 116.4 70 35 4.2 3000 33.13 7268.8
42 36 2.5 200 4.99 668.2 70 35 6.3 50 5.26 5.1
42 36 2.5 500 7.48 2113.9 70 35 6.3 500 12.58 124.8
42 36 2.5 750 10.47 3715.1 70 35 6.3 1000 22.99 441.7
42 36 2.5 1000 12.01 5652.6 70 35 6.3 2000 33.39 1447.0
42 36 2.5 1500 15.36 9897.2 70 35 6.3 3000 41.23 2836.2
42 36 2.5 2000 18.05 15081.6 70 35 8.4 50 6.71 7.7
42 36 2.5 2500 20.79 21089.8 70 35 8.4 500 13.17 115.9
42 36 2.5 3000 23.30 27926.8 70 35 8.4 1000 16.53 323.3
42 36 2.5 3500 25.67 35598.7 70 35 8.4 2000 29.27 990.9
42 36 2.5 4000 27.88 43908.0 70 35 8.4 3000 43.49 2061.0
42 36 3.4 50 4.51 48.5 70 36 3.4 50 4.49 31.0
42 36 3.4 200 5.29 275.7 70 36 3.4 100 5.34 66.3
42 36 3.4 500 9.28 920.2 70 36 3.4 200 6.72 155.0
42 36 3.4 750 12.13 1624.4 70 36 3.4 300 7.24 276.5
42 36 3.4 1000 14.27 2453.6 70 36 3.4 400 7.47 424.1
42 36 3.4 1500 18.12 4448.6 70 36 3.4 500 7.85 586.9
42 36 3.4 2000 21.74 6855.0 70 36 3.4 750 10.61 1091.9
42 36 3.4 2500 25.03 9706.8 70 36 3.4 1000 13.19 1750.5
42 36 3.4 3000 28.03 12925.1 70 36 3.4 1500 18.57 3514.6
42 36 3.4 3500 30.92 16580.6 70 36 3.4 2000 23.37 5839.4
42 36 3.4 3800 32.4 18824.1 70 36 3.4 2500 27.61 8723.8
42 36 3.4 4000 33.55 20657.4 70 36 3.4 3000 30.58 11867.1
42 36 4.2 50 4.62 27.3 70 36 3.4 3500 32.96 15162.6
42 36 4.2 200 6.04 160.4 70 36 3.4 3800 34.36 17381.4
42 36 4.2 500 10.84 536.6 70 36 3.4 4000 35.15 18775.7
42 36 4.2 750 13.81 955.7 77 20 4.2 50 5.20 11.0
42 36 4.2 1000 16.36 1455.3 77 20 4.2 500 11.35 165.4
42 36 4.2 1500 20.95 2667.0 77 20 4.2 1000 15.28 412.8
42 36 4.2 2000 25.12 4163.1 77 20 4.2 2000 24.7 1145.8
42 36 4.2 2500 28.85 5929.9 77 20 4.2 3000 33.78 2203.7
42 36 4.2 3000 32.30 7891.1 77 20 6.3 50 5.71 3.4
42 36 4.2 3500 35.48 10128.1 77 20 6.3 500 13.19 56.0
42 36 4.2 4000 38.39 12629.5 77 20 6.3 1000 19.89 146.9
42 36 6.3 50 4.86 9.0 77 20 6.3 2000 33.23 436.2
42 36 6.3 200 8.26 53.1 77 20 6.3 3000 43.79 834.6
42 36 6.3 500 14.17 190.1 77 20 8.4 50 6.22 1.5
42 36 6.3 750 17.84 343.2 77 20 8.4 500 15.10 26.5
42 36 6.3 1000 21.30 528.4 77 20 8.4 1000 24.56 73.0
42 36 6.3 1500 27.47 995.6 77 20 8.4 2000 40.00 220.7
42 36 6.3 2000 32.85 1573.1 77 20 8.4 3000 51.59 428.4
141
F Numerical Results
sLsTδ Re Nu ∆p sLsTδ L Nu ∆p
[mm] [mm] [mm] [-] [-] [Pa] [mm] [mm] [mm] [mm] [-] [Pa]
42 36 6.3 2500 38.37 2286.2 77 20 8.4 50 5.14 12.3
42 36 6.3 3000 42.23 3086.6 77 25 8.4 500 10.96 191.8
42 36 6.3 3500 46.00 3940.54 77 25 8.4 1000 14.79 531.7
42 36 6.3 4000 49.79 4933.7 77 25 8.4 2000 26.51 1688.8
42 36 8.4 50 5.10 4.2 77 25 8.4 3000 35.92 3212.5
42 36 8.4 200 10.40 25.0 77 25 6.3 50 5.59 3.8
42 36 8.4 500 17.00 92.9 77 25 6.3 500 12.61 66.6
42 36 8.4 750 21.89 172.0 77 25 6.3 1000 20.81 214.5
42 36 8.4 1000 26.27 269.4 77 25 6.3 2000 35.11 635.5
42 36 8.4 1500 33.95 512.7 77 25 6.3 3000 45.99 1243.3
42 36 8.4 2000 40.53 818.1 77 25 8.4 50 6.07 1.7
42 36 8.4 2500 46.40 1185.2 77 25 8.4 500 14.75 32.8
42 36 8.4 3000 51.70 1644.5 77 25 8.4 1000 25.55 107.0
42 36 8.4 3500 56.48 2125.4 77 25 8.4 2000 41.68 326.3
42 36 8.4 4000 61.02 2631.3 77 25 8.4 3000 52.93 668.7
42 40 3.4 50 4.23 52.8 77 30 4.2 50 5.03 13.8
42 40 3.4 200 5.11 290.1 77 30 4.2 500 10.11 228.9
42 40 3.4 500 8.78 946.5 77 30 4.2 1000 15.27 733.5
42 40 3.4 750 11.49 1656.7 77 30 4.2 2000 27.78 2490.5
42 40 3.4 1000 13.36 2464.4 77 30 4.2 3000 36.69 5065.0
42 40 3.4 1500 16.96 4463.3 77 30 6.3 50 5.49 4.2
42 40 3.4 2000 20.39 6897.5 77 30 6.3 500 11.96 87.6
42 40 3.4 2500 23.24 9697.7 77 30 6.3 1000 21.71 303.1
42 40 3.4 3000 26.03 12883.2 77 30 6.3 2000 35.78 997.6
42 40 3.4 3500 28.62 16451.2 77 30 6.3 3000 45.69 2152.7
42 40 3.4 4000 31.08 20400.1 77 30 8.4 50 5.92 1.9
42 45 3.4 50 4.07 56.6 77 30 8.4 500 15.48 45.9
42 45 3.4 200 4.86 297.9 77 30 8.4 1000 26.26 157.4
42 45 3.4 500 7.96 948.7 77 30 8.4 2000 41.92 512.1
42 45 3.4 750 10.23 1631.8 77 30 8.4 3000 53.23 1134.5
42 45 3.4 1000 12.25 2460.0 77 36 3.4 50 4.87 29.0
42 45 3.4 1500 15.35 4384.2 77 36 3.4 500 8.45 507.4
42 45 3.4 2000 18.89 6887.7 77 36 3.4 750 9.89 952.7
42 45 3.4 2500 21.13 9536.7 77 36 3.4 1500 17.98 3098.1
42 45 3.4 3000 23.71 12657.5 77 36 3.4 2000 23.92 5543.6
42 45 3.4 3500 26.19 16187.4 77 36 3.4 2500 28.42 8312.2
42 45 3.4 4000 28.59 20058.7 77 36 3.4 3000 31.79 11692.7
42 50 3.4 50 3.88 59.9 77 36 3.4 3500 34.80 14917.1
42 50 3.4 200 4.68 305.9 77 36 3.4 3800 35.98 17215.9
42 50 3.4 500 7.48 959.1 77 36 3.4 4000 37.11 18637.6
42 50 3.4 750 9.33 1641.2 84 20 4.2 50 5.12 10.6
42 50 3.4 1000 11.05 2437.0 84 20 4.2 500 11.05 155.7
42 50 3.4 1500 14.12 4399.6 84 20 4.2 1000 14.76 383.4
42 50 3.4 2000 16.96 6754.1 84 20 4.2 2000 22.97 1019.9
42 50 3.4 2500 19.55 9448.4 84 20 4.2 3000 30.78 1905.0
42 50 3.4 3000 22.96 12916.3 84 20 6.3 50 5.60 3.3
42 50 3.4 3500 25.12 16297.7 84 20 6.3 500 12.85 52.3
42 50 3.4 4000 27.16 19977.1 84 20 6.3 1000 18.75 133.9
45 30 3.4 50 4.93 38.5 84 20 6.3 2000 30.77 380.1
142
sLsTδ Re Nu ∆p sLsTδ L Nu ∆p
[mm] [mm] [mm] [-] [-] [Pa] [mm] [mm] [mm] [mm] [-] [Pa]
45 30 3.4 1000 16.01 2325.4 84 20 6.3 3000 40.90 723.1
45 30 3.4 2000 24.35 6701.0 84 20 8.4 50 6.09 1.5
45 30 6.3 50 5.56 7.0 84 20 8.4 500 14.57 24.6
45 30 6.3 1000 24.13 527.4 84 20 8.4 1000 22.98 65.5
45 30 6.3 2000 36.87 1500.1 84 20 8.4 2000 37.29 191.3
45 35 3.4 50 4.66 43.6 84 20 8.4 3000 48.47 367.3
45 35 3.4 1000 14.516 2354.8 84 25 4.2 50 5.08 11.7
45 35 3.4 2000 21.96 6606.9 84 25 4.2 500 10.79 177.2
45 35 6.3 50 5.08 78.0 84 25 4.2 1000 14.28 456.4
45 35 6.3 1000 21.54 509.3 84 25 4.2 2000 23.99 1443.8
45 35 6.3 2000 33.12 1508.8 84 25 4.2 3000 33.06 2687.6
45 40 3.4 50 4.43 48.1 84 25 6.3 50 5.50 3.6
45 40 3.4 1000 13.23 2343.7 84 25 6.3 500 12.46 60.6
45 40 3.4 2000 20.05 6504.7 84 25 6.3 1000 19.06 174.2
45 40 6.3 50 4.70 8.8 84 25 6.3 2000 32.45 533.5
45 40 6.3 1000 19.54 498.4 84 25 6.3 3000 42.79 1014.1
45 40 6.3 2000 30.10 1459.7 84 25 8.4 50 5.95 1.6
49 25 4.2 50 5.44 17.2 84 25 8.4 500 14.15 28.9
49 25 4.2 500 12.16 401.5 84 25 8.4 1000 23.81 90.2
49 25 4.2 1000 21.61 1412.1 84 25 8.4 2000 39.24 268.7
49 25 4.2 2000 31.29 4679.8 84 25 8.4 3000 50.14 520.0
49 25 4.2 3000 38.55 9149.9 84 30 4.2 50 5.03 12.9
49 25 6.3 50 6.04 5.5 84 30 4.2 500 10.34 203.9
49 25 6.3 500 17.36 162.8 84 30 4.2 1000 14.03 608.9
49 25 6.3 1000 27.22 577.0 84 30 4.2 2000 25.77 2023.9
49 25 6.3 2000 39.41 1878.7 84 30 4.2 3000 34.61 3994.2
49 25 6.3 3000 49.23 3772.8 84 30 6.3 50 5.42 4.0
49 25 8.4 50 6.59 2.5 84 30 6.3 500 11.73 72.6
49 25 8.4 500 20.99 87.2 84 30 6.3 1000 20.00 247.7
49 25 8.4 1000 32.08 303.2 84 30 6.3 2000 33.81 766.9
49 25 8.4 2000 46.41 997.2 84 30 6.3 3000 43.87 1543.9
49 25 8.4 3000 56.96 1972.8 84 30 8.4 50 5.83 1.8
49 30 4.2 50 5.15 19.86 84 30 8.4 500 14.01 37.7
49 30 4.2 500 11.50 452.2 84 30 8.4 1000 24.72 123.4
49 30 4.2 1000 19.01 1378.1 84 30 8.4 2000 40.18 399.3
49 30 4.2 1000 19.01 1378.1 84 30 8.4 2000 40.18 399.3
49 30 4.2 2000 27.51 4191.8 84 30 8.4 3000 50.64 787.4
49 30 4.2 3000 34.12 8042.1 84 36 3.4 200 6.41 128.6
49 30 6.3 50 6.00 6.2 84 36 3.4 500 9.14 434.9
49 30 6.3 500 15.60 165.0 84 36 3.4 750 9.73 818.2
49 30 6.3 1000 23.86 512.1 84 36 3.4 1500 16.79 2691.2
49 30 6.3 2000 34.70 1571.9 84 36 3.4 2000 22.22 4727.5
49 30 6.3 3000 43.84 3073.3 84 36 3.4 2500 26.73 7045.7
49 30 8.4 50 7.29 9.1 84 36 3.4 3000 30.72 9790.8
49 30 8.4 500 12.74 163.7 84 36 3.4 3500 34.41 12975.2
49 30 8.4 1000 19.05 443.6 84 36 3.4 3800 36.36 15208.9
49 30 8.4 2000 35.85 1338.9 84 36 3.4 4000 37.74 16636.4
49 30 8.4 3000 48.61 2585.2 91 20 4.2 50 5.04 10.2
49 35 4.2 50 4.85 22.3 91 20 4.2 500 10.73 147.6
143
F Numerical Results
sLsTδ Re Nu ∆p sLsTδ L Nu ∆p
[mm] [mm] [mm] [-] [-] [Pa] [mm] [mm] [mm] [mm] [-] [Pa]
49 35 4.2 500 10.77 479.3 91 20 4.2 1000 14.29 358.6
49 35 4.2 1000 16.92 1363.5 91 20 4.2 2000 21.68 926.0
49 35 4.2 2000 24.62 3955.36 91 20 4.2 3000 28.87 1691.6
49 35 4.2 3000 30.78 7463.3 91 20 6.3 50 5.49 3.17
49 35 6.3 50 5.19 7.3 91 20 6.3 500 12.51 49.3
49 35 6.3 500 14.44 175.3 91 20 6.3 1000 17.82 123.3
49 35 6.3 1000 21.59 508.4 91 20 6.3 2000 28.94 338.2
49 35 6.3 2000 31.57 1503.8 91 20 6.3 3000 38.24 632.6
49 35 6.3 3000 40.14 2895.76 91 20 8.4 50 5.95 1.4
49 35 8.4 50 5.50 3.4 91 20 8.4 500 14.09 22.9
49 35 8.4 500 17.57 87.1 91 20 8.4 1000 21.76 59.4
49 35 8.4 1000 25.87 256.8 91 20 8.4 2000 35.05 168.5
49 35 8.4 2000 38.76 779.3 91 20 8.4 3000 45.67 319.12
49 35 8.4 3000 49.83 1530.0 91 25 4.2 50 5.00 11.25
49 36 3.4 50 4.66 41.9 91 25 4.2 500 10.56 165.1
49 36 3.4 100 5.26 98.84 91 25 4.2 1000 13.92 412.45
49 36 3.4 200 5.55 250.0 91 25 4.2 2000 22.00 1175.7
49 36 3.4 300 6.32 418.0 91 25 4.2 3000 30.48 2244.7
49 36 3.4 400 7.67 623.8 91 25 6.3 50 5.39 3.5
49 36 3.4 500 9.03 847.7 91 25 6.3 500 12.23 55.7
49 36 3.4 750 12.03 1538.9 91 25 6.3 1000 17.56 146.3
49 36 3.4 1000 14.43 2337.2 91 25 6.3 2000 30.06 445.4
49 36 3.4 1500 18.48 4284.5 91 25 6.3 3000 40.18 837.3
49 36 3.4 2000 21.89 6635.9 91 25 8.4 50 5.83 1.5
49 36 3.4 2500 25.07 9339.9 91 25 8.4 500 13.75 26.2
49 36 3.4 3000 28.09 12529.4 91 25 8.4 1000 21.99 73.1
49 36 3.4 3500 30.90 16000.6 91 25 8.4 2000 36.57 220.8
49 36 3.4 3800 32.50 18324.7 91 25 8.4 3000 47.57 425.3
49 36 3.4 4000 33.46 19816.8 91 30 4.2 50 4.96 12.3
50 30 3.4 50 4.94 35.3 91 30 4.2 500 10.26 185.4
50 30 3.4 1000 16.05 2209.0 91 30 4.2 1000 13.30 493.7
50 30 3.4 2000 24.99 6543.6 91 30 4.2 2000 23.47 1671.6
50 30 6.3 50 5.64 6.25 91 30 4.2 3000 32.54 3136.1
50 30 6.3 1000 24.51 508.55 91 30 6.3 50 5.33 3.78
50 30 6.3 2000 36.87 1544.6 91 30 6.3 500 11.76 63.8
50 35 3.4 50 4.72 39.7 91 30 6.3 1000 18.30 203.45
50 35 3.4 1000 14.50 2261.6 91 30 6.3 2000 31.69 609.5
50 35 3.4 2000 22.13 6398.0 91 30 6.3 3000 42.12 1212.3
50 35 6.3 50 5.23 7.14 91 30 8.4 50 5.65 1.7
50 35 6.3 1000 21.73 493.1 91 30 8.4 500 13.24 30.9
50 35 6.3 2000 33.13 1458.5 91 30 8.4 1000 23.03 102.3
50 40 3.4 50 4.53 43.9 91 30 8.4 2000 39.07 314.98
50 40 3.4 1000 13.32 2274.7 91 30 8.4 3000 39.07 314.98
50 40 3.4 2000 20.12 6314.5 91 36 3.4 50 4.74 25.8
50 40 6.3 50 4.87 7.9 91 36 3.4 200 6.25 119.4
50 40 6.3 1000 19.66 484.3 91 36 3.4 500 9.38 387.4
50 40 6.3 2000 30.16 1414.2 91 36 3.4 1500 15.42 2242.2
56 25 4.2 50 5.35 15.4 91 36 3.4 2500 24.92 5678.8
56 25 4.2 500 10.93 327.3 91 36 3.4 3000 28.65 7750.4
144
sLsTδ Re Nu ∆p sLsTδ L Nu ∆p
[mm] [mm] [mm] [-] [-] [Pa] [mm] [mm] [mm] [mm] [-] [Pa]
56 25 4.2 1000 20.66 1179.0 91 36 3.4 3500 32.23 10069.2
56 25 4.2 2000 31.71 4164.1 91 36 3.4 4000 35.51 12886.1
56 25 4.2 3000 39.45 8740.9 112 20 4.2 50 4.87 9.58
56 25 6.3 50 5.94 4.9 112 20 4.2 500 9.97 132.7
56 25 6.3 500 15.81 131.2 112 20 4.2 1000 13.29 316.2
56 25 6.3 1000 26.37 473.5 112 20 4.2 2000 19.06 784.4
56 25 6.3 2000 39.07 1711.2 112 20 4.2 3000 25.03 1385.4
56 25 6.3 3000 48.91 3573.5 112 20 6.3 50 5.26 2.9
56 25 8.4 50 6.78 2.2 112 20 6.3 500 11.68 43.8
56 25 8.4 500 19.16 68.1 112 20 6.3 1000 15.98 106.3
56 25 8.4 1000 30.00 238.9 112 20 6.3 2000 25.06 276.5
56 25 8.4 2000 45.14 867.2 112 20 6.3 3000 32.99 501.7
56 25 8.4 3000 57.51 1771.9 112 20 8.4 50 5.67 1.3
56 30 4.2 50 5.13 17.5 112 20 8.4 500 13.07 20.2
56 30 4.2 500 10.73 391.9 112 20 8.4 1000 18.94 50.1
56 30 4.2 1000 19.23 391.9 112 20 8.4 2000 30.41 134.9
56 30 4.2 2000 28.52 4147.4 112 20 8.4 3000 39.64 248.1
56 30 4.2 3000 34.73 7982.5 112 25 4.2 50 4.79 9.2
56 30 6.3 50 5.64 5.6 112 25 4.2 500 9.43 125.1
56 30 6.3 500 15.26 152.5 112 25 4.2 1000 12.53 300.8
56 30 6.3 1000 24.85 515.4 112 25 4.2 2000 17.60 757.6
56 30 6.3 2000 35.61 1618.9 112 25 4.2 3000 22.97 1350.7
56 30 6.3 3000 44.15 3167.25 112 25 6.3 50 5.09 2.8
56 30 8.4 50 6.13 2.5 112 25 6.3 500 11.06 41.37
56 30 8.4 500 19.25 80.5 112 25 6.3 1000 14.86 101.7
56 30 8.4 1000 29.31 270.1 112 25 6.3 2000 23.05 269.0
56 30 8.4 2000 42.68 862.1 112 25 6.3 3000 30.86 496.7
56 30 8.4 3000 53.29 1730.9 112 25 8.4 50 5.43 1.2
56 35 4.2 50 6.45 70.7 112 25 8.4 500 12.32 19.1
56 35 4.2 500 9.61 968.3 112 25 8.4 1000 17.45 48.0
56 35 4.2 1000 11.68 2532.4 112 25 8.4 2000 28.39 132.5
56 35 4.2 2000 17.31 6833.2 112 25 8.4 3000 37.73 247.5
56 35 4.2 3000 25.00 12674.4 112 30 4.2 50 4.77 9.86
56 35 6.3 50 5.31 6.3 112 30 4.2 500 9.32 135.3
56 35 6.3 500 14.16 162.5 112 30 4.2 1000 12.28 331.3
56 35 6.3 1000 22.22 493.7 112 30 4.2 2000 17.29 868.5
56 35 6.3 2000 31.92 1482.8 112 30 4.2 3000 23.27 1682.3
56 35 6.3 3000 39.80 2843.6 112 30 6.3 50 5.05 3.0
56 35 8.4 50 5.70 2.9 112 30 6.3 500 10.91 45.1
56 35 8.4 500 17.66 82.3 112 30 6.3 1000 14.61 113.9
56 35 8.4 1000 26.35 251.6 112 30 6.3 2000 23.26 329.4
56 35 8.4 2000 38.56 769.6 112 30 6.3 3000 32.04 627.6
56 35 8.4 3000 48.74 1510.37 112 30 8.4 50 5.37 1.3
56 36 3.4 50 4.63 37.5 112 30 8.4 500 12.13 20.9
56 36 3.4 100 5.14 83.9 112 30 8.4 1000 17.19 54.8
56 36 3.4 200 6.32 221.3 112 30 8.4 2000 29.239 166.5
56 36 3.4 300 6.61 393.7 112 30 8.4 3000 39.06 310.5
56 36 3.4 400 7.53 582.1 112 36 3.4 50 4.68 22.9
56 36 3.4 500 8.98 778.5 112 36 3.4 200 5.67 102.4
145
F Numerical Results
sLsTδ Re Nu ∆p sLsTδ L Nu ∆p
[mm] [mm] [mm] [-] [-] [Pa] [mm] [mm] [mm] [mm] [-] [Pa]
56 36 3.4 750 11.87 1423.6 112 36 3.4 500 8.86 312.8
56 36 3.4 1000 14.42 2205.5 112 36 3.4 750 10.59 532.8
56 36 3.4 1500 18.90 4144.6 112 36 3.4 1000 11.65 780.9
56 36 3.4 2000 22.44 6445.7 112 36 3.4 1500 13.64 1421.9
56 36 3.4 2500 25.48 9128.3 112 36 3.4 2000 15.88 2257.3
56 36 3.4 3000 28.34 12072.4 112 36 3.4 2500 18.68 3212.7
56 36 3.4 3500 31.02 15599.6 112 36 3.4 3000 22.23 4538.9
56 36 3.4 3800 32.58 17854.8 112 36 3.4 3500 25.65 5947.2
56 36 3.4 4000 33.53 19349.0 112 36 3.4 3800 28.28 7078.7
63 25 4.2 50 5.22 14.1 112 36 3.4 4000 29.23 7534.9
146
Acknowledgements
The present dissertation arose during my activity as a scientific coworker at the Insti-
tute of Thermal Process Engineering and Plant Technology, University of Paderborn,
Germany. Thereby I was supported by many persons and at this place I want to express
my gratitude to all of them.
First of all, I would like to thank my supervisor, Prof. Dr.-Ing Jovan Mitrovi´c,
for his comprehensive support, motivating discussions, encouragement and kindly help,
not only concerning the work at the Insitute, but also in everyday life. He was able to
recognise my professional abilities and helped its further developing.
Dr.rer.nat. Henning Raach iducted me into the world of Computational Fluid Dy-
namics and I would like to thank him for his generous support, kindness and friendliness.
I am grateful to the Ms. Rong Fan, Ms. Yan Li, Mr. Philipp Grimm and Mr.
Waldemar Heckel and Mr. Azer Nazdraji´c, who chose the Computational Fluid Dynam-
ics as a subject of their bachelor’s thesis and student research projects and me as their
tutor. A number numerical experiments that were used for creating the Nusselt number
correlation in this thesis were conducted by them.
Further, I gratefully acknowledge the help that I recieved from all my colleagues,
scientific coworkers, secretaries and technical staff from the Institute. Especially I would
like to thank Mr. Axel Keller from PC2Paderborn center for parallel computing.
Furthermore, I would like to mention all my schoolteachers and university professors
who taught me during my education, especially Prof. Dr.-Ing. Branislav B. Baˇcli´c, from
the University of Novi Sad, Serbia.
And last but not least I would like to thank my wife Sladjana and my sons Darko
und Oliver for their continuous support and encouragement, in spite of the long periods
of my absence they had to put up with.
Paderborn, July 2008
Boban Maleti´c
147
