International Journal
on Marine Navigation
and Safety of Sea Transportation
Volume 3
Number 2
June 2009
173
NOMENCLATURE
cp : specific heat at a constant pressure J/(kg-K)
D : height of cavity, m
g : gravitational acceleration, m
2
/s
h : convection heat transfer coefficient, W/(m
2
-K)
k : thermal conductivity, W/(m-K)
p : pressure, N/m
2
q
′′
: heat flux, W/m
2
t : time, s
T : temperature, K
u : horizontal velocity, m/s
U : dimensionless horizontal velocity
v : vertical velocity, m/s
V : dimensionless vertical velocity
x,y : coordinates, m
X : dimensionless horizontal coordinate
Y : dimensionless vertical coordinate
ΔT : temperature difference (T
H
-T
C
)
α : thermal diffusivity, m
2
/s
β : thermal expansion coefficient, 1/K
θ : dimensionless temperature
µ : dynamic viscosity, kg/(m-s)
ν : kinematics viscosity, m
2
/s
ρ : density, kg/m
3
σ : surface tension, N/m
σT : temperature coefficient of surface tension, N/
(m-K)
τ : dimensionless time
ψ : stream function, m
2
/s
Ψ : dimensionless stream function
ω : vorticity, 1/s
Ω : dimensionless vorticity
Dimensionless Numbers
Pr : Prandtl number
Ra : Rayleigh number
Ma : Marangoni number
Nu : Nusselt number
Nu
: Average Nusselt number
Pe : Peclet number
Subscripts
H : hot
C : cold
l : local
1 INTRODUCTION
From its modest origins, maritime transportation has
always been the dominant support of global trade.
Therewithal, the importance of maritime transporta-
tion increases in parallel with technologic evolu-
tions, technical improvements and economic devel-
opment of countries. High level of economic growth,
industrialization, technological evolutions and ur-
banization for developed countries result in an in-
crease in energy demand. It was determined that in
2005, 86% of primary energy demand in the world
supplied from petroleum and derivatives. Therefore,
maritime transportation is the most important trans-
portation mode for transferring of petroleum and its
derivatives to procure energy demand. With the cur-
rent data illustrated that transportation of fossil fuels
by seaway reached approximately twenty seven bil-
lion tons in the year 2007 (UNCTAD, 2007).
A Numerical Study of Combined Natural and
Marangoni Convection in a Square Cavity
K. Cicek & A. Cihat Baytas
Istanbul Technical University, Tuzla, Istanbul, Turkey
ABSTRACT: Through the aim of this study, the effects of combined buoyancy-driven flows and thermo ca-
pillary flows, which are emerged from temperature differences, on fluid flow and heat transfer numerically
investigated with differentially heated side walls in a free surface square cavity. The study has been accom-
plished with three milestones to achieve the right solutions. For every milestone Navier-Stokes, continuity
and energy equations are discretized by using finite volume method and grids with 52 x 52 control volumes.
Results are presented Pr=1, Pr=7 and Pr=100. The effect of positive and negative Marangoni number on fluid
flow and heat transfer at different Rayleigh number are considered and discussed.
174
Increase in the volume of carriage of petroleum
and its derivatives with maritime transportation have
increased the risk of the maritime accidents such as
oil spill, collision, grounding of ships that threaten to
environment, ecosystems, and aquatic life. For this
reason many academic researches and studies are di-
rectly focused out on prevention of pollution of ma-
rine environment with exploring fluid mechanics
and heat transfer events in cargo tanks of ships such
as Grau et al. 2004, Oro et al. 2006, Pallares et al.
2004, Segerre-Perez et al. 2007.
The main aim of present study is to shed light on
fluid flow and heat transfer in cargo tank of ships.
The present work is a numerical study natural con-
vection due to the temperature difference between
left and right wall and Marangoni convection due to
the free surface effect in a two-dimensional square
cavity. Natural convection is induced by the differ-
ence in temperature between vertical walls, and it is
represented by the Rayleigh number (Ra). Maran-
goni convection flow directly related to the surface
tension gradient with respect to temperature which
acts as a force applied to the free surface of the cavi-
ty, and it is represented by the Marangoni number
(Ma). The presence of free surface can not only alter
the flow field and heat transfer but also prove to
have an impact on the process because of surface
tension variations (Behnia, Stella, Guj 1995; Berg-
man & Ramadhyani 1986; Smith & Davis 1983).
The study is conducted numerically under the as-
sumption of steady laminar flow with three mile-
stones. In the first milestone, natural convection in a
two-dimensional square cavity has been investigated
with the left vertical wall is a constant temperature
T
h
, the right wall is constant T
c
and all other walls
are assumed adiabatic. In the second milestone,
combined natural and Marangoni convection in a
two-dimensional square cavity with a top free sur-
face has been investigated with the left vertical wall
is a constant temperature T
h
, the right wall is con-
stant T
c
, and all other walls are adiabatic. In the third
and last milestone combined natural and Marangoni
convection in a two-dimensional square cavity with
a top free surface has been investigated with the left
vertical wall is a constant temperature T
c
, the right
and bottom wall is constant T
c
, whilst top surface is
adiabatic. For the second and third milestones, the
top surface deformation and interactions with the
gaseous phase are neglected. In every three steps,
the Navier-Stokes, continuity and energy equations
are solved using finite volume method and grids
with 52 x 52 control volumes. Results are presented
Pr=7 for first and second milestones, Pr=100 for
third milestone. The effect of positive and negative
Marangoni number on the fluid flow and heat trans-
fer at different Rayleigh numbers are considered and
discussed.
2 MATHEMATICAL MODEL
A schematic diagram of every milestone shows in
Figure 1, Figure 2 and Figure 3 respectively.
Figure 1. Schematic diagram of the physical model and coordi-
nate system of first milestone.
Figure 2. Schematic diagram of the physical model and coordi-
nate system of second milestone.
Figure 3. Schematic diagram of the physical model and coordi-
nate system of third milestone.
In Figure 1, it is assumed that the left vertical
wall of the cavity is a constant value T
H
. The right
vertical wall is held at a constant temperature T
C
,
while the horizontal walls are adiabatic. In Figure 2,
it is assumed that the left vertical wall of the cavity
is a constant value T
H
. The right vertical wall is held
at a constant temperature T
C,
top horizontal wall is a
175
free surface and adiabatic and bottom horizontal
wall is adiabatic. In Figure 3, cargo tank is modeled
to investigate fluid motion and heat transfer in tank.
In Figure 3, it is assumed that the left vertical wall of
the cavity is a constant value T
H
. The right vertical
wall and bottom horizontal walls are held at a con-
stant temperature T
C
,
while
the top horizontal wall is
a free surface and adiabatic.
To model the liquid motion in cavity, we use the
conservation equation for mass, momentum and en-
ergy for two-dimensional, steady and laminar flow.
All the physical properties of fluid, μ, k and c
p
are
considered constant except density, in buoyancy
term, which obeys Boussinesq approximation. In the
energy conservation, we neglect the effect of com-
pressibility and viscous dissipation. With these as-
sumptions the continuity, momentum and energy
equations can be written as;
0
uv
xy
∂∂
+=
∂∂
(1)
2
dv
gP v
dt
ρρ µ
= −∇ + ∇
(2)
22
22
T T T TT
uv
t x y xy
α
∂ ∂ ∂ ∂∂
++= +
∂ ∂ ∂ ∂∂
(3)
Introducing the vorticity w as;
wv=∇×
(4)
to get vorticity transport equation by taking the curl
of momentum equation to eliminates the
P∇
term,
the momentum equation can be rewritten in terms of
the vorticity defined above as;
2
T
uv g
t xy x
ωωω
νω β
∂∂∂ ∂
++=∇+
∂∂∂ ∂
(5)
The stream function (ψ) for two dimensional prob-
lem is defined such that;
v
ψ
=∇×
(6)
the vorticity transport equation can be obtained from
Eqs. (5) and (6), which further gives;
2
ωψ
= −∇
(7)
The governing equations are converted into the non-
dimensional form by defining the non-dimensional
variables:
2
3
2
, , , V ,
, , , , Pr
s
x y uD vD D
XYU
DD
TT
t g TD
Ra
TD
ω
αα α
ψ αβ ν
θτ
α να α
= = = = Ω=
−
∆
Ψ= = = = =
∆
(8)
Based on these non-dimensional variables, the gov-
erning equations are obtained as follow;
0
UV
XY
∂∂
+=
∂∂
(9)
22
22
Pr PrU V Ra
X Y XY X
θ
τ
∂Ω ∂Ω ∂Ω ∂ Ω ∂ Ω ∂
+⋅ +⋅ = ⋅ + + ⋅⋅
∂ ∂ ∂ ∂∂ ∂
(10)
2
UV
XY
θθθ
θ
τ
∂∂∂
+⋅ +⋅ =∇
∂∂∂
(11)
2
−Ω = ∇ Ψ
(12)
, V=U
YX
∂Ψ ∂Ψ
= −
∂∂
(13)
Other physical quantities of interest in the present
study are the average and local Nusselt numbers for
the hot and cold walls; these variables are defined
respectively as;
() ( )
0
,
l hot l cold
x xD
Nu Nu
XX
θθ
= =
∂∂
=−=−
∂∂
(14)
1
1
0
0
l
Nu dY Nu dY
X
θ
∂
=−=−
∂
∫∫
(15)
3 INITIAL & BOUNDARY CONDITIONS
In order to obtain results of the conservation equa-
tions we define initial and boundary conditions. For
each milestone initial conditions are at
0
τ
=
;
0UV
θ
= = =Ψ=Ω=
But it is necessary to define boundary conditions of
each milestone separately. In first milestone the
boundary conditions are at
0
τ
≥
;
X = 0,
01Y≤≤
,
0 & 1UV
θ
= =Ψ= =
(16a)
X = 1 ve
01Y≤≤
de
0 & 0UV
θ
= =Ψ= =
(16b)
Y = 0,
01X≤≤
,
0 & 0UV
Y
θ
∂
= =Ψ= =
∂
(16c)
Y = 1 ve
01X≤≤
de
0 & 0UV
Y
θ
∂
= =Ψ= =
∂
(16d)
Also it is necessary to define non-dimensional vorti-
city boundary conditions and with the help of Eqs.
(4) the boundary conditions can be obtained as for
vertical walls;
2
2
X
∂Ψ
Ω=−
∂
(16e)
and for horizontal walls;
2
2
Y
∂Ψ
Ω=−
∂
(16f)
176
In the second and third milestones free surface
condition must be taken in consideration. The pres-
ence of this type of boundary condition can not only
alter the flow field and heat transfer characteristics
but may also prove to have an impact on the process
because of surface tension variations. (Smith & Da-
vis 1983; Bergman & Ramadhyani 1986). In gen-
eral, for most liquids there is a variation of surface
tension with temperature. As a result, the interface
between air and liquid which is subjected to a tem-
perature gradient can initiate a bulk flow due to sur-
face tension variations. This flow which is due to a
temperature gradient applied normally to the free
surface is known as Marangoni convection (Behnia
et al 1995). After that explanation, the boundary
condition at a free surface can be written as,
1
yD
yD
uT
y Tx
σ
µ
=
=
∂ ∂∂
=⋅⋅
∂ ∂∂
(17)
Based on non-dimensional variables expressed in
Eqs (8), Eqs 17 are obtained as follow;
U
Ma
YX
θ
∂∂
=
∂∂
(18)
,
T
d TD d
Ma
dT dT
σσ
σ
µα
∆
=−=
Eqs. (18) can be rearranged with the help of
, V=U
YX
∂Ψ ∂Ψ
= −
∂∂
and obtained as follow;
2
2
Ma
YX
θ
∂Ψ ∂
=
∂∂
(19)
This boundary condition applied to a vorticity
boundary condition for top horizontal wall as follow;
2
2
Ma
YX
θ
∂Ψ ∂
Ω=− =−
∂∂
(20)
In the second milestone initial and boundary con-
ditions are same with the first milestone except top
horizontal wall (free surface condition).
In the third milestone initial and boundary condi-
tions are same with the second milestone except bot-
tom horizontal wall. In this milestone the bottom
wall is held to a constant temperature T
c
. For this
reason the boundary condition for bottom wall is de-
fined as follow,
Y = 0,
01X≤≤
,
0 & =0UV
θ
= =Ψ=
(21)
4 SOLUTION PROCEDURE
The differential equations, represented by equations
(9) to (11), together with respect boundary condi-
tions for every milestone, equations are (16), (20)
and (21), are solved using the fine volume method
described in Patankar (1980). In this method solu-
tion domain is divided into small finite control vol-
umes. The differential equations are integrated into
each of those control volumes. From this integration
there were algebraic equations which, when solved
simultaneously or separately, supplied velocity,
stream function, vorticity and temperature compo-
nents. A power-law scheme is adopted for the con-
vection-diffusion formulation (Patankar 1980).
The discretization equations are solved iterative-
ly, using the line by line method known as Thomas
algorithm or TDMA (tridiagonal matrix algorithm).
An over relaxation parameter of 1.85 was used in
order to obtain stable convergence for the solution of
vorticity transport and energy equations.
In order to assess the accuracy of our numerical
procedure for every milestone, we have tested our
algorithm based on the grid size (42 X 42) for the
first milestone with the work of Davis (1983a), Da-
vis (1983b) and Hortman and Peric (1990), for the
second milestone we have tested our algorithm
based on the grid sizes (52 X 52) with the work of
Behnia et al (1995).
Grid independence tests were conducted for all
the configurations studied in this work. Three differ-
ent grid size (42 X 42, 52 X 52, 62 X 62) were used
and average Nusselt numbers are compared for dif-
ferent grid size. Because of the small differences be-
tween 52 X 52 and 62 X 62 grids, 52 X 52 grid was
chosen for second and third milestones. For the first
milestone 42 X 42 grid was chosen because of
achieve effective benchmark with the work of Davis
(1983a), Davis (1983b) and Hortman and Peric
(1990). The numerical solution is considered to be
converged when the maximum convergence criteria
was smaller than 10
-7
5 NUMERICAL RESULTS AND DISCUSSION
In Figure 4 illustrate the stream function and iso-
therm contours of first milestone numerical results
for various Ra=10
3
and 10
5
and Pr=1. In general, flu-
id circulation is strongly dependent on Rayleigh
number as we have seen in Figure 4a and 4b.
Change in the values of Ra has its influence on the
stream lines and isotherms. The centre of circulation
pattern moved towards to hot (left) wall of the cavi-
ty. Streamlines pattern shows that there is a strong
upward flow near the hot wall side and downward
flow near the cold wall.
177
(a)
(b)
Figure 4. Stream lines and isotherm contours for (a) Ra=103,
(b) Ra=105
In Figure 5-7 illustrate the stream function and
isotherm contours of second milestone numerical re-
sults for Ra=0 and 10
4
, Ma=
±
10
3
and Pr=7.
Figure 5. Stream lines and isotherm contours for Ra=0,
Ma=10
3
Figure 6. Stream lines and isotherm contours for Ra=10
4
,
Ma=10
3
Figure 7. Stream lines and isotherm contours for Ra=10
4
,
Ma=-10
3
In general, fluid circulation is strongly dependent
on both Rayleigh number and Marangoni number as
we have seen in Figure 5-7. Change in the values of
Ra and Ma have influence on the stream lines and
isotherms. In Figure 5 only the effects of surface
tension are present, the Rayleigh number being set to
zero. The flow is driven due to effect of velocity
gradients on the top surface. In Figure 6 it is seen
that the effect of buoyancy on the flow field is added
to the positive Marangoni one. So the flow structure
is still constituted by a single large main circulation.
Due to the effect of buoyancy forces that are distrib-
uted in the bulk of the cavity the main vortex is
stronger than in Figure 5 and encompasses almost
completely the entire flow domain. For positive Ma
there is a strong thermal boundary layer located
close to the upper-right corner due to the effect of
Marangoni convection in the vicinity of the free sur-
face. In Figure 7, due to the negative value of Ma-
rangoni number, there is an opposition of effects be-
tween buoyancy and surface tension. For this reason,
it is possible to clearly distinguish the effect of the
two different contributions. In this case, it is clearly
evident how the flow is strongly divided into two re-
gions; the upper one where the motion is induced by
surface tension and the lower one driven by buoyan-
cy. Also in Figure 7 isotherm contours shows that
the negative Ma causes a strong penetration of hot
fluid in the centre of the cavity due to the combined
effect of the two counter-rotating recirculation.
In Figure 8-10 illustrate the stream function and
isotherm contours of third and last milestone numer-
ical results for Ra=0 and 10
4
, Ma=
±
10
3
and Pr=100.
Figure 8. Stream lines and isotherm contours for Ra=0,
Ma=10
3
Figure 9. Stream lines and isotherm contours for Ra=10
4
,
Ma=10
3
0.00 0.10 0. 20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
0.00 0.10 0 .20 0.30 0. 40 0.50 0. 60 0.70 0.80 0. 90 1.00
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
0.00 0.10 0 .20 0.30 0. 40 0.50 0. 60 0.70 0.80 0. 90 1.00
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
0.00 0.10 0.20 0.30 0. 40 0.50 0.60 0. 70 0.80 0.90 1. 00
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
0.00 0.10 0 .20 0.30 0. 40 0.50 0. 60 0.70 0.80 0. 90 1.00
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
178
Figure 10. Stream lines and isotherm contours for Ra=10
4
,
Ma=-10
3
The bottom wall which is held at a constant tem-
perature T
C
doesn’t make any changes on stream
line contours but have influence on the isotherm
contours. As a result there is a strong boundary layer
located close to lower-left corner. In Figure 8-10
stream lines contours don’t show any difference but
isotherm contours with the influence of bottom wall
show significance varieties.
6 CONCLUSION
This study try to shed light on fluid flow and heat
transfer in cargo tank of ships with investigating of
combined natural and Marangoni convection in a
liquid.
For the case of positive Marangoni number the
fluid structure is mainly constituted by a dominant
core vortex structure which becomes stronger as ei-
ther of the Rayleigh and Marangoni numbers are in-
creased.
More complex and interesting is the flow struc-
ture for negative values of Marangoni numbers. In
this case, the fluid is stratified and flow filed is
clearly divided into two separate regions where the
motion is driven by Marangoni (upper region) and
buoyancy (lower region) respectively.
Also for the bottom wall at a constant tempera-
ture T
c
, it is shown that there is a strong thermal
boundary layer is formed at the lower-left corner. It
is clearly evident of lower-left corner is a strong heat
loss point.
Consequently, Marangoni number is as important
parameter as Rayleigh number for fluid motion and
dispersion of temperature in cavity with free surface
and for positive and negative values it changes both
stream lines and isotherms contours dispersion in the
cavity.
REFERENCES
Behnia, M., Stella, F., & Guj, G. 1995. A numerical study of
three dimensional combined buoyancy and thermocapillary
convection, International Journal of Multiphase Flow,
21(3): 529-542.
Bergman, T.L. & Ramadhyani, S. 1986. Combined buoyancy
and thermocapillary driven convection in open square cavi-
ties. Journal of Numerical Heat Transfer, 9: 441-451.
Bergman, T.L. and Keller, J.R., 1988. Combined Buoyancy,
Surface Tension Flow in Liquid Metals, Numerical Heat
Transfer, 13, 49-63.
Davis, G.D.V., 1983a, Natural Convection of Air in a Square
Cavity: A Bench Mark Numerical Solution, International
Journal for Numerical Methods in Fluids, 3, 249-264.
Davis, G.D.V., 1983b, Natural Convection of Air in a Square
Cavity: A Comprasion Exercise, International Journal for
Numerical Methods in Fluids, 3, 227-248.
Faghri, A. and Zhang, Y., 2006. Transport phenomena in mul-
tiphase systems. Academic Press, London.
Grau, F.X., Valencia, L., Fabregat, A., Pallares, J. and Cuesta
I., 2005. Modelization and simulation of the fluid dynamics
of the fuel in sunken tankers and of the dispersion of the
fuel spill, Symposium on Marine Accidental Oil Spills, Vi-
go, Spain, July 13-16.
Hortman, M. and Peric, M., 1990. Finite volume multigrid pre-
diction of laminar narutal convection: bench-mark solution,
International Journal for Numerical Methods in Fluids, 11,
189-207.
Oro, J.M.F., Morros, C.S. and Diaz, K.M.A., 2006. Numerical
simulation of the fuel oil cooling process in a wrecked ship,
Journal of Fluid Engineering, 128, 1390-1393.
Pallares J., Cuesta I. and Grau F.X., 2004. Numerical simula-
tion of the fuel oil cooling in the sunken prestige tanker.
The ASME-ZSIS International Thermal Science Seminar
II, Bled, Slovenia, June 13-16, 439-445.
Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow,
Hemissphere Publishing Corporation, Washington D.C.
Segerra-Perez, C.D., Olivia, A., Trias, X., Lehmkuhl, O. and
Capdevila, A., 2007. Numerical simulation of thermal and
fluid dynamic behavior of fuel oil in sunken ships, Sympo-
sium on Marine Accidental Oil Spills, Vigo, Spain, June 5-
8, 31.
Smith, M.K. & Davis, S.H. 1983. Instabilities of dynamic
thermocapillary liquid layers, part 1. Convective instabili-
ties. Journal of Fluid Mechanics, 132: 119-144.
0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00
0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
