International Journal of Heat and Mass Transfer: Faroogh Garoosi, Tew-Fik Mahdi

International Journal of Heat and Mass Transfer 172 (2021) 121163
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer

journal homepage: www.elsevier.com/locate/hmt
Presenting a novel higher-order bounded convection scheme for

simulation of multiphase flows and convection heat transfer
Faroogh Garoosi∗, Tew-Fik Mahdi
Department of Civil, Geological and Mining Engineering, Polytechnique Montreal, Montreal, Quebec, Canada
a r t i c l e i n f o a b s t r a c t
Article history: The primary aim of the current study is to enhance the stability and accuracy of the Volume-Of-Fluid
Received 2 December 2020 (VOF) method for modeling free-surface flows with large topological changes and high density ratio.
Revised 2 February 2021
For accurate capturing of fluid interfaces, a novel higher-order bounded convection scheme is first con-
Accepted 27 February 2021
structed based on the total variation diminishing (TVD) concept and is then employed for the discretiza-
Available online 24 March 2021
tion of convection terms in Navier-Stokes, energy and transport equations. In the second step, the clas-
Keywords: sical PISO algorithm is modified according to the two-step projection method (Chorin’s model) and the
Third-order TVD flux-limiter scheme combined model (PISOC) is then applied for the treatment of the pressure-velocity coupling. Moreover,
Improved PISO algorithm the second-order accurate piecewise-linear interface reconstruction technique (PLIC-ELVIRA) is used for
Multiphase flows determining the normal direction and curvature of the interface. The robustness and accuracy of the pro-
PLIC;VOF method posed models in handling multiphase flows with interface rupture and coalescence are verified against
ELVIRA scheme
several experimental and numerical benchmark solutions such as: dam break, Rayleigh-Taylor instability,
bubble rising, rotation of a slotted disk (Zalesak’s problem), deformation of a 2D disk and pure con-
vection of a step profile. The results show that, the proposed third-order TVD flux-limiter scheme can
considerably reduce the false-diffusion errors and ensure the boundedness of the volume fraction while
retaining the sharpness and shape of the interface. Furthermore, it is found that the proposed PISOC
algorithm has strong stability and convergence characteristics in strongly coupled multiphase problems
and is less susceptible to divergence when larger pressure under-relaxation factor is used. The perfor-
mance and effectiveness of the proposed modifications are further demonstrated by analyzing transient
entropy generation due to conjugate natural convection heat transfer in two different canonical test cases
(i.e. Differentially Heated Cavity and Rayleigh-Bénard Convection) and good agreements are found with
previously published works.
© 2021 Elsevier Ltd. All rights reserved.
1. Introduction cavities. In a pioneering work, Amine et al. [6] experimentally and

numerically investigated flow structure and heat transfer character-
The topic of conjugate conduction and natural convection heat istics of natural convection in a parallelepipedic cavity containing
transfer in closed enclosures has received considerable scholarly several solid blocks and concluded that the number and arrange-
attention in the last decade owing to its wide range of engineering ment of conductive bodies have a significant impact on the flow
applications, including building insulation, designing nuclear reac- intensity within the enclosure. Similar observations were reported
tor and solar energy collectors, heat exchangers, double pane win- by Zhang et al. [7], Ren et al. [8], Raji et al. [9] and Garoosi et al.
dows and microelectronic cooling devices [1–4]. It is well-known [10,11] who numerically analyzed the influence of the block’s frag-
that, the presence of internal conductive obstacles in the case of mentation upon the heat transfer rate and direction of the fluid
buoyancy driven flows causes significant alteration in the flow and motion inside the square cavity. They found that the subdivision
heat transfer characteristics of the thermal systems [5]. Over the of the conductive obstacle into the smaller segments tends to sup-
past decade, extensive effort has been devoted to analyzing and press the fluid motion and delay the onset of natural convection. In
optimizing conjugate natural convection heat transfer inside closed the same context, Hu et al. [12], Miroshnichenko et al. [13], Selime-
fendigil et al. [14], Wang et al. [15] and Zargartalebi et al. [16] sys-
tematically examined the influences of pertinent parameters such
∗
Corresponding author. as Prandtl number (Pr), Rayleigh number (Ra), inclination angles
E-mail addresses: [email protected] (F. Garoosi), of the enclosure, the relative thermal conductivity and thickness
tewfi[email protected] (T.-F. Mahdi).
https://doi.org/10.1016/j.ijheatmasstransfer.2021.121163
0017-9310/© 2021 Elsevier Ltd. All rights reserved.
F. Garoosi and T.-F. Mahdi International Journal of Heat and Mass Transfer 172 (2021) 121163
sure with differentially heated horizontal walls where temperature

Nomenclature gradient and gravity force are parallel but act in opposite direc-
tions (Rayleigh-Bénard convection) [18]. The study of convection
A surface area per unit depth A = 2(L + H ), m heat transfer in the semi concentric annulus such as C-shaped cav-
Be Bejan number ity is of great interest among researchers because it involves the
Cp specific heat, J kg-1 K-1 combination of two aforementioned groups. Kalidasan et al. [19],
Fb Body force Rahimi et al. [20], Makulati et al. [21] and Mliki et al. [22] nu-
g Gravitational acceleration, ms-2 merically investigated the hydrodynamic and thermal behavior of
h Smoothing length, m Rayleigh-Bénard convection in a two-dimensional C-shaped enclo-
H enclosure height, m sure. They proved that the onset of Rayleigh Bénard convection
k thermal conductivity,Wm-1 K-1 strongly depends upon the Ra and aspect ratio of the cold obsta-
L Length of the enclosure (m) cle. However, although the first-law analysis of thermodynamic can
Nu Average Nusselt number on the active walls provide useful information about flow structures and heat distribu-
p pressure, Nm-2 tion within the system, minimizing the total rate of entropy gener-
P dimensionless pressure ation is necessary for maximizing the overall energy efficiency of
Pr Prandtl number (= ν f /α f ) the thermal systems. The basic concept of Entropy Generation Min-
r distance imization (EGM) was first introduced by Bejan [23,24] who math-
Ra Rayleigh number (= gβ (Th − Tc )H 3 /αν ) ematically demonstrated that irreversibilities due to fluid friction
S Total entropy generation (W/m3 K) (viscous dissipation) and heat transfer (temperature gradients) are
SF Entropy generation due to viscous dissipation two main sources of entropy generation and exergy loss (destruc-
ST Entropy generation due to heat transfer tion of available work) in thermal systems. Al-Rashed et al. [25],
t time (s) Pordanjani et al. [26], Alkanhal et al. [27], Sivaraj et al. [28] and
t∗ dimensionless time, (t α /H 2 ) Magherbi et al. [29] applied second law of thermodynamics to in-
T temperature, K vestigate the effects of Ra on the entropy generation during the
T +T
T0 Bulk temperature, ( h 2 c ), K natural convection heat transfer inside the square cavity with and
u,v velocity components, ms-1 without internal heating. They found that as the Ra enhances, both
U, V Dimensionless velocity components average Nusselt number and entropy generation augment, signif-
x,y Cartesian coordinates, m icantly. A comprehensive review of latest efforts on the analysis
xij Distance between particle i and j in x-direction of entropy generation and natural convection heat transfer within
(xi j = x j − xi ) various enclosures can be found in works of Das et al. [30] and
yij Distance between particle i and j in y-direction Biswal et al. [31].
(yi j = y j − yi ) The principles of convection heat transfer and entropy gener-
X, Y Dimensionless Cartesian coordinates ation are reasonably well understood and documented, however,
modeling the hydrodynamics of multiphase flow with strong dis-
Greek symbols continuities is another challenging issue faced by the Computa-
α thermal diffusivity, m2 s-1 tional Fluid Dynamic (CFD) community. Since, multiphase flows
β Thermal expansion coefficient, K-1 with moving and deformable interfaces are ubiquitous in various
θ dimensionless temperature modern engineering technologies with different length scales such
μ dynamic viscosity, kg m-1 s-1 as coal gasification in molten slag [32], fluidized beds [33–35],
ν kinematic viscosity, m2 s-1 solidification-melt dynamics [36], combustion systems [37], atom-
ρ density, kg m-3 ization of liquid jets [38–40] and droplet impact [41,42], consider-

k curvature of interface
Y able research efforts have been devoted to developing robust nu-
ψ stream function(= − Y U ∂ Y + ψ (X, Y0 ))
0 merical approach for simulation of such flows. Generally, simula-
μT0 α 2
ϕ Irreversibility ratio ( ) or general depen- tion of multiphase flows can be accomplished through two differ-
kH 2 (Th −Tc )2
ent main frameworks, namely Eulerian and Lagrangian descriptions
dent variable
[43]. To date, most of the successful methods for investigations of
Subscripts multi-fluid and multiphase flows are designed based on Eulerian
c cold wall mesh-dependent techniques which in turn have their own pros
f fluid and cons. Among a large variety of numerical approaches in the
h hot wall class of Eulerian description, Volume-of-Fluid (VOF) [44], Level-
Total (tot) summation over the domain Set (LS) [45], Front Tracking model [46] and Two-Fluid approach
[47] are four top-ranking methods which have been widely used
for computation of multiphase flows with moving free surface. Re-
of the solid wall on the flow pattern, temperature distribution and cent applications of the VOF and the level set methods for sim-
heat transfer rate, for the case of conjugate natural convection in ulating bubble rising and dam-break flows involving several ob-
different geometries (square cavity and partitioned horizontal an- stacles can be found in works of Malgarinos et al. [48], Issakhov
nulus). They showed that, an increase in the wall-to-fluid thermal et al. [49,50] and Gu et al. [51]. Their results demonstrated that,
conductivity ratio leads to an intensification of convective circula- although two aforementioned methods can satisfactorily predict
tion and thereby heat transfer enhancement while thickness of the the pressure distribution of gas-liquid two phase flows, excessive
conductive wall has a minor impact on the average Nusselt num- growth in the thickness of the interface (or blurring of the free
ber especially at high Ra. surface) is still one of the main difficulties associated with these
Generally, from the applied thermal boundary conditions point models. To overcome this shortcoming, over the last two decades,
of view, natural convection mechanism in closed enclosure can be several researchers have focused on developing a robust and ac-
classified into two different groups: (a) enclosure with differen- curate Interface-Sharpening Algorithm for treatment of the mate-
tially heated vertical walls (DHC) [17] where direction of the buoy- rial interfaces in immiscible multiphase flows. Regardless of the
ancy force is perpendicular to temperature gradient, and (b) enclo- numerical approach used for simulation, the problem of preserv-
2
ing interface resolution, in general, can be classified into three ma- terms and transport equation. To enhance the convergence rate
jor groups, namely (Ⅰ) implementing high resolution differencing and consistency of the iterative solution procedure, the classical
schemes for discretization of the convective term [52], (Ⅱ) artificial PISO algorithm is combined with the two-step projection method
compression schemes to suppress the growth of interface thickness (Chorin’s model) and the proposed hybrid model (PISOC) is then
[53], and (Ⅲ) Adaptive Mesh Refinement (AMR) technique to han- implemented for the treatment of the pressure-velocity coupling.
dle large topology changes of the interface [54]. To preserve a high level of numerical accuracy, the ELVIRA scheme
The first class of methods attempts to define a set of rules is utilized for estimating the interface normal and curvature. The
with trade off precision-boundedness and its use in the simul- proposed modifications are validated against several benchmark
taneous reduction of numerical diffusion and unphysical oscilla- test cases including dam break, Rayleigh-Taylor instability, bubble
tions with a rapid, efficient, and cost-effective one-step compos- rising, rotation of a slotted disk (Zalesak’s problem), deformation
ite procedure. Pioneering work in this direction was carried out of a 2D disk, pure convection of a step profile and transient en-
by Leonard [52] through introducing the principle of the Nor- tropy generation due to conjugate natural convection heat transfer
malized Variable diagram (NVD). He constructed a new higher- in square enclosures.
order bounded convection scheme to mitigate false-diffusion er-
rors and suppress spurious oscillations (undershoots or overshoots 2. Problem statement and governing equations
that leads to non-physical negative value of density) near discon-
tinuities where strong gradients of variables exist. In the light of Schematic diagram of the physical models considered in this
the aforementioned study, Hill et al. [55], Denner et al. [56], Ling paper along with the corresponding boundary conditions are
et al. [57], Deng et al. [58] and Das et al. [59] applied differ- shown in Fig. 1 where cases 1 to 6 are employed for valida-
ent high order monotonicity-preserving schemes such as: CICSAM, tion and verification of the proposed models in dealing with
Hyper-C, Semi-Implicit Skewness Correction (SISC), MUSCL, HRIC the multiphase flow problems, whereas cases 7 and 8 are con-
and SMART to maintain the shape and sharpness of the interface sidered for analysis of transient entropy generation due to con-
in the complex convection dominated problems including those of jugate conduction-convection heat transfer in square enclosures
the rising bubble and melting with natural convection. Their re- filled with air (Pr=0.71). In all cases, the unsteady flow is assumed
sults showed that, implementing high-resolution bounded convec- to be two-dimensional and laminar. The thermo-physical proper-
tive scheme in transport equation plays a crucial role for deter- ties of the fluid are taken to be constant except for the density in
mining accurate topology features of the surface dynamics in the cases 7 and 8 which varies linearly with the temperature in accor-
coalescence and breakup processes. A comprehensive review and dance to the Oberbeck-Boussinesq approximation. The set of gov-
comparison of the behavior between eight widely-used oscillation- erning equations including continuity, momentum, energy and vol-
free convection schemes in a number of benchmark cases can be ume fraction for incompressible Newtonian fluids can be written in
found in works of Alves et al. [60] and Choi et al. [61]. Referring to the following dimensional forms [76,77]:
the second approach, Heyns et al. [53] introduced the concept of
Artificial Compression Velocity to remove false-diffusion errors and ∂ρ ∂ u ∂v
+ + = 0, (1)
restrict the thickness of the interface zone to the minimum pos- ρ ∂t ∂ x ∂ y
sible number of computational grids. Pozzetti et al. [62], Cifani
et al. [63], Nguyen et al. [64,65], So et al. [66] and Akhlaghi et al.
∂ ρm u ∂ ρm uu ∂ ρm vu ∂p ∂ ∂u ∂ ∂u
[67] extended former work and combined the anti-diffusion and + + =− + μm + μm + FST
∂t ∂x ∂y ∂x ∂x ∂x ∂y ∂y
artificial compression velocity techniques to control the interface
layer thickness between the lighter and heavier mediums. They (2)
solved an additional partial differential equation for volume frac-
tion field and concluded that the aforementioned technique can

considerably prevent blurring and wrinkling of moving interface, ∂ ρm v ∂ ρm uv ∂ ρm vv ∂p ∂ ∂v ∂ ∂v
+ + =− + μm + μm + Fg + FST
thereby reducing excessive numerical dissipation in that area. An- ∂t ∂x ∂y ∂y ∂x ∂x ∂y ∂y
other alternative approach to deal with the excessive smearing of (3)
interface is to use an Adaptive Mesh Refinement and coarsening
(AMR). A local AMR algorithm as a third-principle tool was orig-
inally proposed by Berger et al. [54] for solving hyperbolic sys- ∂T ∂T ∂T ∂ 2T ∂ 2T
ρ Cp +u +v =k + , (4)
tems of conservation laws. This technique was employed later by ∂t ∂x ∂y ∂ x2 ∂ y2
Hua et al. [68], Tan et al. [69], Ngo et al. [70], and Antepara et al.
[71] to further minimize non-physical smearing of the interface
∂ϕ ∂ϕ u ∂ϕv
in modeling of rising bubble and droplet dynamics. However, the + + + ∇ · (ϕ (1 − ϕ )uR ) = 0 (5)
results of Chen et al. [72], Schmidmayer et al. [73] and Li et al.
∂t ∂x ∂y
[74] showed that improper management of unnecessary elements The energy balance equation for pure solid region in cases 7
in the dynamic grid refinement/coarsening process can remarkably and 8 can be expressed as:
increase computational effort and even jeopardize convergence of 2
∂T ∂ T ∂ 2T
non-linear solvers. These findings were also supported by Das et al. = αr + (6)
[59] and Harvie et al. [75] who stated that embedding higher order ∂t ∂ x2 ∂ y2
differencing schemes (i.e. 5th order WENO) in local AMR technique
In the above equations, u and v are the mixture velocities in
without precaution can lead to unphysical velocities around the in-
x and y-directions and Cp , αr and k denote specific heat capac-
terface and substantial reduction of mass conservation accordingly.
ity, thermal diffusivity ratio (αr = αs /α f ) and thermal conductiv-
Based on the above literature review, the primary objective
ity, respectively. ρm and μm are the mixture density and viscos-
of the present study is to enhance the accuracy and stability of
ity of working fluids which can be computed via the volume-
Volume-Of-Fluid (VOF) model for simulation of multiphase flows
weighted average of the phases occupying each computational grid
with high density ratios. For this purpose, the third-order accu-
(0 < ϕ < 1) [78]:
rate flux-limiter function is first constructed based on the TVD
constraint and then employed for discretization of the convection ρm = ϕ ρ 1 + ( 1 − ϕ ) ρ 2 (7)
3
Fig. 1. Schematic diagram of eight different physical configurations under consideration with associated boundary conditions and coordinate system.
with n standing for the interface normal vector. The additional

μm = ϕ μ1 + ( 1 − ϕ ) μ2 (8) term ∇ · (ϕ (1 − ϕ ) uR ) in the phase advection equation is Artifi-
cial Compression [53] which helps in lessening the smearing of the
where subscripts 1 and 2 refer to the primary and secondary interface and retaining the sharpness of the interface region. Note
phases, respectively. The term Fg in Navier-Stokes equations is the that, the presence of ϕ (1 − ϕ ) term in Eq. (5) allows the Interface
body force which would be the gravity force (Fg = ρm g) in the Compression to be activated only in close vicinity of the interface
multiphase flow problems (cases 1 to 3) or buoyancy force (Fg = area. uR is the artificial compressive velocity which is given by:
ρ gβ (T − Tc )) in conjugate free convection problems (cases 7 and −
→ −
→ ∇ϕ
8) with Tc being reference temperature. FST is the volumetric sur- uR = Cϕ u
n = Cϕ u (12)
face tension force which acts on the interface region based on the |∇ϕ |
Brackbill’s model [79] as: Cϕ is a constant coefficient used to control the strength of the
= σ
compression velocity with the recommended value of unity [67].
FST k∇ϕ (9)
Once the temperature and velocity fields are determined, the total
where σ represents the interfacial tension coefficient between entropy generation (Stot ) per unit volume for cases 7 and 8 can be
phases and
k denotes the curvature of interface defined as: calculated via the second law of thermodynamics as follows [80]:

2
2
−
→ → k ∂T
∂T
· ∇ ) − −
→
1 n
n = −
k = −∇ · → ( −
→ n − (∇ · n ) (10) Stot = ST + SF = 2 +
n n T0 ∂x
∂y

2 2 2
−
→ μ ∂u ∂v ∂ u ∂v
n
n = −
, −
→
n = ∇ϕ (11) + 2
∂x
+2
∂y
+ +
∂y ∂x
(13)
→
n T0
4
where ST and SF are two main detrimental irreversibilities induced the pressure nodes are sited at the cell centroids while the ve-
by thermal dissipation and fluid friction, respectively. By defining locity components are located on the midpoints of control volume
the following pertinent parameters, Eq. (13) can be transformed faces, as depicted in Fig. 2 [83]. Motivated by works of Issa [117],
into the non-dimensional form as [81]: Chorin [118], Kim et al. [119], Nasri et al. [120] and Varsakelis et al.
[121], a novel hybrid PISO-Two-step projection algorithm (PISOC) is
2 , θ = T −T .
2
T −Tc
X= x
, Y = Hy , U = uH vH
α ,V = α , P
pH
= ρα
H h c
(14) first developed and then employed for handling pressure-velocity
tα gβ (Th −Tc )H 3 μT0 α 2
t∗ = , Ra = αν , Pr = αν . λ= 2 2, coupling between the mass and momentum equations. A more de-
H2 kH (Th −Tc )

2 2 tailed discussion and derivation procedure of the pressure-based
∂θ ∂θ coupled algorithm associated with the newly developed hybrid
ST = + (15) PISOC method are provided in appendix A (please see supplemen-
∂X ∂Y
tary material).

2 2 2 The classical second-order central difference scheme is used
∂U ∂V ∂U ∂V for discretization of the diffusion terms while the proposed third-
SF = λ 2 +2 + + (16)
∂X ∂Y ∂y ∂X order TVD bounded scheme is utilized for the treatment of the
convection terms in momentum, energy and volume fraction equa-
In Eq. (16), λ is the ratio between the viscous and thermal dis- tions. To further improve the performance and accuracy of the VOF
sipation which is set equal to 10-4 similar to works of Das et al. model, the ELVIRA technique pioneered by Pilliod et al. [84] is ap-
[81] and Ilis et al. [82]. The average total entropy generation (Stot ) plied for material interface reconstruction and its curvature calcu-
is then obtained via integrating the local entropy generation rates lation.
over the entire computational domain as:
3.1. A novel Third-order TVD Bounded scheme
1 1
ST = ST dX dY SF = SF dX dY Stot = ST + SF (17)
V V
In this subsection, a new third-order oscillation-free flux-limiter
The relative dominance of thermal dissipation to fluid friction function is developed based on the TVD constraint for convec-
irreversibility can be mathematically quantified using the local and tion discretization and the suppression of the interface smearing
average Bejan numbers defined as: in dealing with multifluid and multi-material flow problems. In or-

ST 1 der to simplify the derivation process of higher-order monotone
Be = Be = Be dX dY (18) convection scheme in the finite-volume description, the general
ST + SF V
convection-diffusion equation in Cartesian framework without the
Therefore, according to the above definition, one can conclude
pressure and transient terms is considered here as:
that for the values of 0 ≤ Be < 0.5 the fluid friction irreversibility
is dominant (ST < SF ) whereas 0.5 < Be ≤ 1 indicates dominance d d d dφ d dφ
of thermal dissipation (SF < ST ). For the particular case of Be =
(ρ uφ ) + (ρvφ ) = η + η (21)
dx dy dx dx dy dy
0.5, entropy generation due to both factors are identical (SF = ST ).
Meanwhile, the global heat transfer rate within the thermal sys- where the scalar (T, ϕ ) and vector quantities (u, v) in Eqs. (2) to
tems can be quantified through the average Nusselt number along (5) are replaced by the general transport variable φ . By integrating
the heated/cooled walls as follows: Eq. (21) over a two-dimensional control volume (CV) as sketched in
Fig. 2 and applying the Gauss’s divergence theorem in conjunction
hH ∂θ 1 H
∂θ with the classical central differencing scheme for diffusion term,
Nu = =− − Nu = dn (19)
k ∂→
n H 0 ∂n the following discretized equation can be derived:
As depicted in Fig. 1, the related initial and boundary conditions Fe φe − Fw φw + Fn φn − Fs φs = De (φE − φP ) − Dw (φP
for the present problems can be specified as follows:
− φW ) + Dn (φN − φP ) − Ds (φP − φS ) (22)
u = v=0 On the rigid walls of the enclosures in all cases where Fi = (ρ u )i Ai with A i = x = y and Di =
θ = 1 On the heated walls of the enclosures in cases 7 and 8 ηi y/x =ηi x/y stand for the convection flux and diffusion
θ = 0 On the cooled walls of the enclosures in cases 7 and 8 conductance cross the cell boundaries, respectively. By defining
∂θ the flux-limiter function (ψ (r )) and taking into account the flow
→ = 0 On the adiabatic walls of the enclosures in cases 7 and 8
∂−
n direction (u+ > 0 or u− < 0), the values of the property φ across
(20) the cell faces (φe , φw , φn and φs ) can be approximated as follows:
⎧
k f ( ∂θ ∂θ ) On the walls of the conductive obstacle in
→ ) f = ks ( − ⎪
1
φe = φP + ψ (re+ )(φE − φP ), re+ = (φP − φW )/(φE − φP )
∂−n ∂→n
s ⎪
⎪
⎪ 2
⎪
⎪
cases 7 and 8 ⎨φw = φW + 1 ψ (rw+ )(φP − φW ), rw+ = (φW − φWW )/(φP − φW )
∂p n +1
→ = 0 On the rigid walls of the enclosures in all cases
2 u+ > 0 (23)
∂− ⎪
⎪ 1
φn = φP + ψ (rn+ )(φN − φP ), rn+ = (φP − φS )/(φN − φP )
n ⎪
⎪
The initial conditions (at t = 0) inside the enclosures are based ⎪
⎪
2
⎩ 1
φs = φS + ψ (rs+ )(φP − φS ), rs+ = (φS − φSS )/(φP − φS )
on the assumptions that the cold working fluid is motionless (θ = 2
0, u = v = 0) under the zero-pressure filed ( p = 0). It is worth ⎧ 1

⎪
⎪φe = φP + ψ (re− )(φP − φE ), re− = (φE − φEE )/(φP − φE )
⎪
⎪ 2
mentioning that, although the governing equations are discretized ⎪
⎪
⎨φw = φP + 1 ψ (rw− )(φW − φP ), rw− = (φP − φE )/(φW − φP )
on the staggered-grid system but to close the Pressure Poisson 2 u− < 0 (24)
⎪
⎪ 1
φn = φN + ψ (rn− )(φP − φN ), rn− = (φN − φNN )/(φP − φN )
Equation (PPE) in the two-step projection part of the algorithm, the ⎪
⎪
−
→ ⎪
⎪
2
⎩
homogeneous Neumann condition for the pressure (∂ pn+1 /∂ n = 1
φs = φP + ψ (rs− )(φS − φP ), rs− = (φP − φN )/(φS − φP )
2
0) is imposed on the rigid walls during each time step.
where ri± stands for the ratio of upwind-side gradient to
3. Numerical methodology downwind-side gradient as described by Versteeg et al. [83]. More
+
precisely, for the positive flow direction, the gradient ratio of rw on
The system of governing partial differential equations together the West cell face is evaluated by ri+−1/2 = (φi−1 − φi−2 )/(φi − φi−1 )
with corresponding boundary conditions is discretized on a stag- (or ri+−1/2 = (φW − φW W )/(φP − φW )) whereas in the case of neg-
gered grid system by means of a control-volume method where ative flow direction, this ratio would be equal to ri−−1/2 =
5
Fig. 2. (a) Sweby diagram showing the proposed model along with the four other well-known convection schemes (the shaded areas represent the the TVD region). (b)
The 2D staggered grid arrangment used for discritisation of the governing equations where the scalar quantities (i.e. pressure) coincide with the cell faces of the u-control
volume.
(φi − φi+1 )/(φi−1 − φi ). In Eqs. (23) and (24), the term ψ (r ) is (monotonicity) of the solution and at least first-order accuracy of
known as a flux-limiter function which is designed by the com- the discretization scheme [85]. Moreover, similar to Van Leer lim-
bination of two piecewise non-linear functions as: iter function (ψ (r ) = (r + |r| )/(1 + r )) [89], the proposed model
⎧ passes through the points ψ (0 ) = 0 and ψ (1 ) = 1 which fulfills
⎨ (r + |r| ) × (r + 3 ) r≤1 consistency requirement for second order accuracy [88]. It is also
ψ (r ) = 3r 2 + 2r + 3 (25)
⎩ 2r + 2 r>1
easy to show that the slope of the limiter function at the intersec-
r+3 tion point (1,1) is the same as that of the classical QUICK scheme
To construct the above third-order accurate scheme, the neces- (ψ (r ) = (r + 3 )/4, dψdr(r ) |r=1 = 14 ) which in turn fulfills the desired
sary and sufficient criteria suggested by Sweby [85], Zijlema et al. third-order accuracy [52,86,90]. Furthermore, it can be seen from
[86], Leonard [52], Gao et al. [87], Waterson et al. [88] and Alves Fig. 2 that unlike the Minmod and Superbee models [91] (as the
et al. [60] are utilized to retain the monotonicity-preserving prop- lower and upper borders of the TVD region) which switch be-
erties of TVD scheme. It is evident from Fig. 2 that, the newly de- tween linear schemes with different slopes, the proposed model
veloped limiter is positive (ψ (r ) ≥ 0) and lies within the shaded is smooth, finite and continuous over the entire range of r. In
area in the TVD diagram which guarantees both the boundedness this context, it is relevant to mention that, both sub-functions of
6
the proposed limiter reach the point (1,1) with the same slope where γi is the constant parameter which depends on the flow di-
of magnitude ( dψdr(r ) |r=1 = 14 ) which in turn respects the convec- rection as:
tive stability criterion and fulfills Smooth Region Condition, accord-
γe = 1 if Fe > 0 , γe = 0 if Fe < 0
ingly. The effects of this crucial criterion on the consistency of
γw = 1 if Fw > 0 , γw = 0 if Fw < 0
the convection scheme were systematically investigated by Alves (30)
γn = 1 if Fn > 0 , γn = 0 if Fn < 0
et al. [60] and Gao et al. [87] who mathematically proved that,
γs = 1 if Fs > 0 , γs = 0 if Fs < 0
the sudden alteration in the slope of limiter function near the
critical point (1,1) can lead to spurious pressure oscillation near To provide deeper insight into the computer programming of
the material discontinuities in the multiphase flows. Before pro- proposed model, the FORTRAN code in term of subroutine is de-
ceeding further, it is worth to mention that in order to prevent veloped and documented as a supplementary material which can
or minimize excessive use of “IF” in computer programming pro- assist the user of a CFD code in constructing and implementing
cess and to avoid singularity of the denominator of subdomains such a TVD scheme.
(r ≤ 1, r > 1) in Eq. (25), special treatment should be applied to
each sub-domain of the limiter. For this purpose, by multiplying 3.2. Interface reconstruction scheme (ELVIRA algorithm)
the numerator and denominator of gradient ratio (r) with the value
of its denominator (e.g. r = A/B = A × B/B2 ), the following relation- The programming process of ELVIRA technique pioneered by
ships can be drawn: Pilliod et al. [84] is elaborated in this subsection. Generally, in the
⎧ PLIC-ELVIRA model, the interface is approximated by a straight line
⎪ (φ − φi−2 ) (φi−1 − φi−2 ) (φi − φi−1 )
⎪ri+−1/2 = rw+ = i−1
⎪ = × segment (y = mi x + d) where six possible gradients of the color
⎪
⎪ (φi − φi−1 ) (φi − φi−1 ) (φi − φi−1 )
⎪
⎪ ( φ − φ ) × ( φ − φ ) function (ϕ ) generated by the forward, central and backward finite-
⎪
⎨ =
i−1 i−2 i i−1
± (φi − φi−1 ) 2 difference schemes in the x and y-directions are used to determine
rw =
⎪ (φi − φi+1 ) (φi − φi+1 ) (φi−1 − φi ) the interface slope. As depicted in Fig. 3, the first three possible
⎪
⎪r −
= r −
= = ×
⎪
⎪
i−1/2 w
(φi−1 − φi ) (φi−1 − φi ) (φi−1 − φi ) angular coefficients of the straight line in the x-direction can be
⎪
⎪ (φ − φi+1 ) × (φi−1 − φi )
⎪
⎩ = i mathematically calculated by:
(φi−1 − φi )2
mx, f = ±|ϕEN + ϕE + ϕES − ϕN − ϕP − ϕS |
By applying the above technique on the West cell face and con- mx,c = ±|ϕEN + ϕE + ϕES − ϕW N − ϕW − ϕW S | × 0.5 (31)
sidering the flow direction, the sub-domains of the limiter can be mx,b = ±|ϕN + ϕP + ϕS − ϕW N − ϕW − ϕW S |
rewritten as follows:
(φi−1 −φi−2 )×(φi −φi−1 )
ri+−1/2 > 1,
(φi −φi−1 )2
> 1, (φi−1 − φi−2 ) × (φi − φi−1 ) > (φi − φi−1 )2
r>1 (φi −φi+1 )×(φi−1 −φi )
ri−−1/2 > 1, > 1, (φi − φi+1 ) × (φi−1 − φi ) > (φi−1 − φi )2
(φi−1 −φi )2
(φi−1 −φi−2 )×(φi −φi−1 ) (27)
ri+−1/2 ≤ 1,
(φi −φi−1 )2
≤ 1, (φi−1 − φi−2 ) × (φi − φi−1 ) ≤ (φi − φi−1 )2
r≤1 (φi −φi+1 )×(φi−1 −φi )
ri−−1/2 ≤ 1, ≤ 1, (φi − φi+1 ) × (φi−1 − φi ) ≤ (φi−1 − φi )2
(φi−1 −φi )2
Note that, the similar procedure should be applied to other The other three underlying interface slopes can be computed by
three faces of computational cell (re± , rn± and rs± ) by considering the rotating the gradients of row sums by 90 degrees given by:
related neighboring grid points. Finally, since the high-resolution
±1
schemes are typically prone to a numerical instability owing to the my, f = |ϕW N +ϕN +ϕEN −ϕW −ϕP −ϕE |
appearance of negative main coefficients, the proposed model can my,c = |ϕW N +ϕN +ϕEN±2
−ϕW S −ϕS −ϕES | (32)
be generalized and reformulated by incorporating the notation of my,b = |ϕW +ϕP +ϕE −±1
ϕW S −ϕS −ϕES |
upwind differencing scheme into the discretized equation as fol-
lows: It is important to note that in the above equations, the symbol
(±, plus-minus sign) should be interpreted based on the Table 1
aP φP = aW φW + aE φE + aS φS + aN φN + Su and Fig. 3. Take as an illustration, for cases A1, A2, C1 and C2
aW = Dw + max(Fw , 0 ) only the positive (+mi ) values of Eqs (31) and (32) should be used
aE = De + max(−Fe , 0 ) for interface reconstruction whereas for cases B1, B2, D1 and D2
(28)
aS = Ds + max(Fs , 0 ) the negative ones (−mi ) are regarded as potential slope candidates.
aN = Dn + max(−Fn , 0 ) Fig. 3 demonstrates various possible cases of interface orientation
aP = aW + aE + aS + aN + (Fe − Fw ) + (Fn − Fs ) with the corresponding fluid position within the uniform stencil of
3 × 3 cells which can be easily distinguished from each other by
One of the main advantages of the above format of discretiza-
defining the following two pertinent parameters:
tion is that the main coefficient (aP ) is always positive which can
effectively satisfy the needs for conservativeness and transportive- nx = (ϕEN + 2ϕE + ϕES ) − (ϕW N + 2ϕW + ϕW S ) (33)
ness. Su on the right-hand side of equations is called deferred cor-
rection which can be expressed as:
ny = (ϕW N + 2ϕN + ϕEN ) − (ϕW S + 2ϕS + ϕES ) (34)
1
Su = Fe (1 − γe )ψ (re− ) − γe ψ (re+ ) (φE − φP )
2 It is interesting to mention that, nx and ny in the above equa-
1 tions are very analogous to the Youngs’s model [92] who ini-
+ Fw −(1 − γw )ψ (rw −
) + γw ψ (rw+ ) (φP − φW ) tially introduced the concept of the PLIC technique for linear re-
2
1 construction of material interfaces. It is evident from Table 1 and
+ Fn (1 − γn )ψ (rn− ) − γn ψ (rn+ ) (φN − φP ) Fig. 3 that, depending on the signs of nx , ny and the absolute values
2
1 of the interface slope (|mi |), the twelve different interface config-
+ Fs −(1 − γs )ψ (rs− ) + γs ψ (rs+ ) (φP − φS ) (29) urations are likely to emerge within the main cell. Furthermore,
2
7
Fig. 3. Possible cases of interface orientation together with an example of interface reconstruction technique.
Table 1
Determination of the interface orientation and fluid position
variation rate of the color variation rate of the color Absolute value of Possible case of
function in the x-direction function in the y-direction the interface slop interface orientation
nx >0 ny < 0 |mi | > 1 A1

nx >0 ny < 0 |mi | ≤ 1 A2
nx >0 ny > 0 |mi | > 1 B1
nx >0 ny > 0 |mi | ≤ 1 B2
nx <0 ny > 0 |mi | > 1 C1
nx <0 ny > 0 |mi | ≤ 1 C2
nx <0 ny < 0 |mi | > 1 D1
nx <0 ny < 0 |mi | ≤ 1 D2
nx >0 ny ≈ 0 |mi | > 1020 V1
nx <0 ny ≈ 0 |mi | > 1020 V2
nx ≈0 ny < 0 |mi | < 10−20 H1
nx ≈0 ny > 0 |mi | < 10−20 H2
8
it is apparent that, each of the given position (except for hori- actual one which implies that the second iteration is needed to
zontal and vertical orientations) consists of three subgroups of in- reach a desired solution. Since, ϕTrial 1 < ϕActual , the lowest value
terface shapes (i.e. triangle, trapezoid or pentagon) which should of y-intercept is first replaced by the guessed one computed from
be determined through an iterative method until a desired re- the previous iteration, and then the averaging process (dT rial 2 =
sponse is achieved. This stage in the literature is commonly re- (1.4 + 1.05 )/2 = 1.225) is repeated. It should be mentioned that,
ferred to as Determination of Line Constant [93] where the line seg- in the case of ϕTrial 1 > ϕActual , the highest value of line constant
ment (y = mi x + d) is displaced iteratively within the stencil while should be replaced by the guessed value obtained from the previ-
keeping its alignment constant until the volume fraction enclosed ous trial. It can be seen from Fig. 3 that, the newly constructed in-
by the reconstructed interface line becomes equal to the actual terface line (yT rial 2 = 0.3x + 1.225) can successfully reproduce the
volume fraction predicted by transport equation. When the exact correct volume fraction. Having determined the true line constant,
location of the interface line is identified, the constructed line is the volume fractions in the neighboring nodes can be estimated as
then extended across the 3×3 block to determine the volume frac- depicted in Fig. 3. The above procedure should be applied to the
tion (ϕ ) in each of the surrounding 8 grids. In the last stage, any other five remaining slopes and the one that can minimize the fol-
candidate slopes that can minimize the following discrete error is lowing residual is the winner:
regarded as the solution:
mi = min Ei, j (mx, f ), Ei, j (mx,c ), Ei, j (mx,b ), Ei, j (my, f ), Ei, j (my,c ),

1
Ei, j (mi ) = ϕi+k, j+l (mi ) − ϕi+k, j+l (35) Ei, j ( my,b ) (36)
k,l=−1 Once the optimal interface line is determined, its normal vector
(ny y + nx x = d) can be calculated as follows:
Where Ei, j (mi ) is the total error between the actual volume
fractions obtained from the transport equation (ϕi+k, j+l ) and those mi y + x = mi d (37)
predicted by the interface reconstruction (ϕi+k, j+l (mi )). To shed
further light on the structure of computer programming of the where ny = mi and nx = 1 are the components of normal vector. Fi-
ELVIRA model, an analytical example is given here where the nally, similar to the previous sections, in order to provide deeper
linear function with formula y = 0.3x + 1.225 intersects the ver- insight into the computer programming of ELVIRA technique, the
tical sides of the nine-point cell stencil. The volume fractions FORTRAN code is developed and documented as supplementary
generated by this straight line is portrayed in Fig. 3 where the material.
upper half of the stencil is filled by the primary phase. As
a first step of the algorithm, the six possible interface slopes 4. Results and discussion
are calculated through Eqs. (31) and (32) namely: (mx, f = ±3/10,
mx,c = ±3/10, mx,b = ±3/10, my, f = ±480/947, my,c = ±384/571 In this section, eight different canonical test cases are employed
and my,b = ±960/961). In the second step, the interface orien- to demonstrate the accuracy and reliability of the proposed modi-
tation and position of the heavier fluid are determined using fications in handling convection heat transfer and non-hydrostatic
Eqs. (33) and (34). It can be seen from Table 1 that, according to free surface flow with complex interface morphologies and topo-
the values of nx = −1127/960 < 0 and ny = 763/192 > 0 (irrespec- logical changes. A classical 2D dam-break problem (case 1) is first
tive of |mi |), the existing cases for the test problem would be C1 adopted to verify the capability of the enhanced VOF model in
and C2 which in turn indicate that only positive values of inter- capturing moving free surfaces involving severe interfacial defor-
face slopes should be considered for the interface reconstruction mation, coalescence and breaking-up events associated with wa-
process namely: mx, f = +3/10, mx,c = +3/10 mx,b = +3/10, my, f = ter re-entry, impact pressure and splashing. To further assess the
+480/947, my,c = +384/571 and my,b = +960/961. For the first at- performance of the proposed third-order TVD bounded convection
tempt, mx, f = +3/10 is chosen as the angular coefficient of the scheme in suppressing false diffusion errors and spurious numer-
straight line (yInter f ace = 0.3x + d). Meanwhile, by taking into ac- ical smearing (i.e. smearing of the interface thickness), the non-
count the role of |mx, f | < 1 as a third condition along with nx < 0 linear development of 2D Rayleigh-Taylor instability is analyzed
and ny > 0 in Table 1, the possible interface mode is reduced to as the second benchmark test problem (case 2) where two im-
case C2. In the next stage, the shape of the fluid polygon among miscible fluids are evolved under the combined action of viscous
the three possible cases (i.e. triangle, trapezoid or a pentagon) and buoyancy forces. Simulation of a single gas bubble rising in a
should be identified. To accomplish this, the area of the biggest liquid column (case 3) is conducted to validate the effectiveness
possible fluid triangle and trapezoid enclosed in the central cell are of the ELVIRA technique in accurate estimation of surface tension
first calculated and then compared with the actual value of vol- force and its corresponding interface curvature. The classical Ro-
ume fraction. It can be seen from Fig. 3 that the volume fraction tation of a slotted disk (case 4), deformation of a 2D disk (case
obtained by the exact solution is positioned between the volume 5) and pure convection of a step profile (case 6) are chosen to
fraction of the biggest possible triangle and trapezoid (ϕTriangle < demonstrate the capability of the newly-developed TVD bounded
ϕActual < ϕTrapezoid ) which implies that the interface polygon has a convection scheme in capturing steep gradients. The robustness
shape of trapezoid. In the fourth stage, the linear function should and stability of the newly proposed semi-iterative PISOC-algorithm
be calibrated within the main grid to reproduce the correct volume in dealing with the velocity-pressure coupling problem is further
fraction (|ϕActual − ϕInter f ace | ≤ 10−10 ). As mentioned earlier, this demonstrated through analysis of transient entropy generation due
stage is equivalent of finding the line constant (d). From Fig. 3 one to conjugate natural convection heat transfer in square enclosures
can conclude that, according to ϕTriangle < ϕActual < ϕTrapezoid , the (cases 7 and 8).
line constant of interface should be located between y-intercepts
of the largest possible triangle and trapezoid polygons (0.7 ≤ d ≤ 4.1. Dam break (case 1)
1.4) created with the formulas yT riangle = 0.3x + 1.4 and yT rapezoid =
0.3x + 0.7, respectively. Hence, the arithmetic average of the high- The schematic diagram of the dam break model is given in
est and lowest y-intercepts (dT rial 1 = (1.4 + 0.7 )/2 = 1.05) is taken Fig. 1. In this numerical test, the rectangular water column (ρW =
as the first trial in the iterative process. However, it is evident 997 Kg/m3 and μW = 855 × 10−6 Pa.s) with width and height of
from Fig. 3 that the volume fraction produced by the new in- the W = 0.25 m and H = 0.5 m in hydrostatic equilibrium is ini-
terface line (yInter f ace = 0.3x + 1.05) is still lower than that of the tially installed on the LHS of the square container with dimensions
9
Fig. 4. Snapshots of the volume-fraction field and pressure contour with flow sequence identification simulated by the proposed model for case 1.
of [1m, 1m]. The rest of the enclosure is occupied by air (ρA = the water level height at the dam site (l1 = 0.1m) monotonically
1 Kg/m3 and μA = 185 × 10−7 Pa.s) as the outer medium (second lessens which implies that the conversion between potential and
phase) and the surface tension force (σ = 71.1 × 10−3 N/m) is taken kinetic energy takes place primarily at this stage. The second stage
into account. The simulation is conducted on 350 × 350 grid reso- involves the formation of the shock pressure due to impact of wa-
lution and time sequence of the collapsing water column in terms ter toe with the opposite vertical wall which in turn is accom-
of the volume fraction fields together with the pressure contours panied by the development of the ascending fluid jet and semi-
∗
plotted in Fig. 4 at different dimensionless time instants (t =
are hydrostatic pressure distribution on the bottom right corner of the
t g/H ). Generally, the dynamic evolution of the dam failure can enclosure (1.5 < t ∗ ≤ 2.26). A close inspection of Fig. 6 (c) shows
be classified into four distinct stages. The first stage encompasses that, the maximum impact pressure recorded by sensors 2 (h2 =
the destruction of the vertical gate and propagation of the wave 0.025m) at t ∗ = 1.59 s is approximately equal to p∗ = 2.057 which
front over the dry-deck until the surge front reaches the down- is in accordance with previously published measurements [94,95].
stream rigid wall (0 ≤ t ∗ ≤ 1.5). From Fig. 6 one can observe that, Furthermore, it can be seen that, the maximum run up of the liq-
10
Fig. 5. Comparisons of the interfaces at different time instants between present work and numerical results of Xu et al. [94] and Zheng et al. [95].
uid jet in the vicinity of the right wall exceeds twice the initial wa- travel towards the upstream left wall where the water splash-up
ter column height which is compatible with the numerical findings process and interface fragmentation/breaking are likely to occur on
of Colagrossi et al. [96]. However, as time progresses (2.26 < t ∗ ≤ the tip of the second plunging jet (t ∗ < 4.42).
3.42), due to the oncoming fluid flow and retarding effect of the Qualitative and quantitative comparison of the predicted results
gravity force, the fluid acceleration declines and subsequently the with numerical works of Xu et al. [94] and Zheng et al. [95] in Figs.
thickness of the upward-moving jet increases. This process corre- 5 and 6 show that, the trajectory of the tip of the plunging jet and
sponds to the continuous reduction in the pressure shock loading air entrapment phenomenon (closed cell structure) are well cap-
and establishment of the stagnation point (fluid trapping) on the tured and reproduced by the present model. The close-up snap-
corner of the enclosure. The third stage takes into account the ap- shots of the interface before and after wave breaking in Fig. 6 (a)
pearance of the plunging breaking wave where, due to the restor- clearly demonstrate that the sharpness (thickness) of the interface
ing action of gravity and the fluid-solid interaction adjacent to the is well controlled and confined to a maximum of 2-3 grid sten-
right wall, the rising water jet begins to overturn backwards and cil which in turn highlights the robustness of the proposed non-
eventually hits the underlying moving wetted deck, thereby form- linear flux limiter scheme (Eq. (25)) in eliminating interface smear-
ing the closed air-cushion structure and second shock pressure sce- ing and numerical diffusion errors. In order to better quantify the
nario (t ∗ = 4.3, p∗ = 0.704) within the container (3.42 < t ∗ ≤ 4.42). accuracy of the numerical model, the time variations of the water
Finally, in the last stage, the newly generated flying jet starts to front (X f ront ), non-dimensional values of water level height (h/H )
11
Fig. 6. Comparison of the predicted water column heights, fluid front and pressure variations on the verical walls between present study and numerical works of Xu et al.
[94] and Zheng et al. [95] for case 1.
at three different sections and pressure time histories recorded [1m, 2m] where a denser fluid with physical properties of ρH =
at six different sampling points (hi ) on the left and right vertical 1.8 Kg/m3 and μH = 0.018 Pa.s lies on top of a lighter fluid (ρH =
walls are presented in Fig. 6. It is evident that the predicted results 1 Kg/m3 , μH = 0.01 Pa.s) at t = 0 s with a locally perturbed inter-
are smooth without non-physical oscillations and are in excellent face of y = 1 − 0.15 × sin(2π x ). The instabilityis governed by the
agreement with those obtained by reference numerical solutions Reynolds and Atwood numbers of Re = ρH H Hg/μH = 420 and
[94,95]. However, there are some small discrepancy between the At = (ρR − 1 )/(ρR +1 ) = 2/7 where H = 1m denotes the width of
results which can be partially attributed to neglecting the surface the enclosure and Hg stands for the reference velocity with g
tension or due to the occurrence of the so-called tensile instability being the gravitational acceleration. Initially, the system is at rest
problem and its consequent particle clustering in the Lagrangian (u = 0 m/s) and no-slip boundary condition is imposed on all rigid
particle methods which trigger pressure oscillations in published walls. Similar to works of Rezavand et al. [97], Pahar et al. [98] and
results of [94,95]. Li et al. [99], owing to the relatively large radius of interface cur-
vature, the influence of the surface tension force is neglected. The
computations are conducted on a 200 × 400 grid resolution and
4.2. Rayleigh-Taylor instability (case 2)

predicted results are presented in Figs. 7 and 8 at various non-
dimensional time (T = t g/H ) instants. Generally, due to buoy-
As schematically shown in Fig. 1, the simulation is car-
ancy/gravity forces and initial distortion, the dynamics of the RTI
ried out in the rectangular enclosure with dimension of the
12
Fig. 7. Comparison of the predicted results with the numerical results of Rezavand et al. [97], Pahar et al. [98] and Li et al. [99] for the case of Rayleigh-Taylor instability
problem (case 2) at different time instants.
problem is characterized by falling of the denser fluid as a spike time goes on (2 < T ≤ 4), the heavy fluid penetrates further into
along the left wall and upward movement of the lighter fluid as a the less dense fluid and consequently the intensity of the trans-
bubble, which lead to the appearance of the mushroom-like struc- verse velocity augments and exceeds the critical value (see also
ture within the enclosure. Inspection of Figs. 7 and 8 reveals that, |umax | in Fig. 8 (d)). At this stage, the interface evolves into an
during the early stages of the RTI development (T ≤ 2), the con- asymmetric shape and undergoes severe morphological deforma-
tours of the velocity components and shapes of the bubble/spike tion which signifies that shear velocity and Kelvin-Helmholtz insta-
fronts remain entirely symmetrical with respect to the center of bility are coming into the picture [100]. This process is also accom-
the enclosure which implies that morphology and the growth of panied by the formation of the four secondary recirculating zones
the instability are primarily governed by the viscous force. As along the wings of the ascending and descending plumes which
13
Fig. 8. (a) and (b) Contours of the velocity in the x and y-directions at different time instants. (c) Maximum and minimum positions of the interface within the computational
domain, (d) time history of absolute values of maximum velocity components and (e) maximum positions of the interface at three different vertical sections.
causes the interface to roll up, resulting in the significant twisting with large deformations and topological changes. It is evident that
and distortion of the interface. Towards the end of the calculation the interface between two fluids is smooth and successfully cap-
(4 < T ≤ 6), due to the imposition of the no-slip and impermeable tured without any smearing. The time history of the bubble front
boundary conditions, more and more eddies are generated at the and spike tip is depicted in Fig. 8 (c) and good agreement is found
tails of the roll-ups and magnitude of the velocity field progres- with reference numerical measurements [98].
sively lessens which indicates that the RTI evolution enters into the
strongly nonlinear regime. In this circumstance, the major portion 4.3. Bubble rising (case 3)
of the nonlinear spike advances with the vortex shedding and cur-
vature effects of the ring eddies, leading to the development of the To further evaluate the performance of the ELVIRA technique
jet-like circulation and complex flow structure in the lower part of in calculating the interface-normal vector and its corresponding
the container. Qualitative comparison of obtained results with pre- surface-tension force (Eq. (10)), the hydrodynamic characteristics
viously published data [97–99] in Fig. 7 vividly confirms the ca- of a single rising bubble in quiescent liquid is examined in this
pability of the modified VOF model in handling multi-fluid flows subsection. As graphically shown in Fig. 1, the simulation is car-
14
Fig. 9. Comparison of the simulated time evolution of a single rising bubble (case 3) in two-dimensional rectangular enclosure with the numerical results of Zhao et al.
[101].
ried out in a rectangular duct with geometrical dimensions of evolution of the bubble trajectory with associated velocity contours
√ H, 6.48H] where the circular bubble (ρL = 1 Kg/m and μL =
[4 3 are presented in Figs. 9 and 10. It can be seen that, due to the
2/200 Pa.s) with nominal radius of H = 0.25m is initially im- existence of the density gradient between two phases, the buoy-
mersed
√ in a stationary heavy fluid (ρH = 100 Kg/m3 and μH = ancy force causes the air bubble to float up along the centerline
2/2 Pa.s) which occupies the entire enclosure. No-slip and im- of the enclosure while the heavy fluid sinks to the lower section
permeable boundary conditions are imposed on all rigid walls and at the vicinity of the vertical walls, creating a pair of counter-
calculations are performed on the 200 × 324 uniform mesh which rotating streamwise vortices on both sides of the bubble. How-
is fine enough to capture large deformation and breakup of the in- ever, during the floating process, owing to the downwelling mo-
to work of Zhao et al. [101], the Reynolds
terface. Similar num- tion of the denser fluid and shedding of two vortices from the
ber (Re = D 2Hg/υH ) based on the reference velocity ( 2Hg) and main body, the cores of the primary eddies are dragged towards
characteristic length (D = 2H) is set to Re = 100 where υH = υL de- the lower portion of the enclosure and subsequently the bubble
note the kinematic viscosity of the denser and lighter fluids. The shape begins to evolve into the horseshoe-like structure (2 ≤ t ∗ ≤
Weber number (W e = ρH g(2H ) /σ ) based on the properties of the
2
4.5). Fig. 10 reveals that this interface stretching/twisting process
surrounding liquid is chosen as 200, corresponding to the surface is also accompanied by significant enhancement in the flow in-
tension coefficient σ = 0.5 N/m. The components of the initial ve- tensity and the development of the wake structure in the rear of
locity vector and background pressure are set to zero ( p0 = 0 and the gas phase which causes the tails of the bubble to fold inward
u0 = 0 ms−1 ) and the obtained results in terms of the temporal and elongate in the y axis. This mechanism signifies that the ef-
15
Fig. 10. Contours of the velocity field together with the time history of the bubble front, maximum position of the interface at two different sections, and maximum absolute
values of velocity components in case 3.
fects of the shearing-velocity become more prominent as time ad- cending bubble. Finally, in the last stage (5 < t ∗ ≤ 6), due to both
vances. Note that this stage corresponds to maximum rate of bub- drift and lift forces generated by vortex shedding (wake-induced
ble deformation such that the bubble shape is reconfigured to the forces), two detached bubbles are directed towards the vertical
cap-like structure with terminal velocity (|Vmax | = 3.69). However, cross-section and begin to follow the convection path of the rolled-
Fig. 9 shows that, during the deceleration stage (4.5 < t ∗ ≤ 5), the up vortices. The corresponding velocity fields at this moment re-
bubble motion and its rate of deformation abruptly slow down semble approximately dipole-shaped pattern which in turn man-
and remain nearly constant, indicating that the fluid attached to ifest the formation of some secondary recirculating eddies inside
the vertical walls is likely to be in the regime of creeping flow the rising bubble.
[102,103]. Meanwhile, the shearing-stress in conjunction with sur- A comparison of calculated results with the previously pub-
face tension force starts to stretch the interface into thin layer. It lished data based on the VOF solution on unstructured trian-
is evident that during this sequence (4 < t ∗ ≤ 5), the bubble skirts gular grids with adaptive mesh refinement [101] in Figs. 9 and
become progressively squeezed and thinner, resulting in the occur- 10 demonstrates that, despite the negligible differences be-
rence of bubbles detachment and interface rupture behind the as- tween the two numerical outcomes, the hydrodynamic behavior
16
Fig. 11. Comparison between predicted results from the current work with the numerical data of Scardovelli et al. [93] and Owkes et al. [105] for case 4 (Zalesak’s disk).
of rising bubble including wake generation, interface filamenta- shown in Fig. 1, the calculation is performed in a square enclo-
tion/fragmentation together with the bubbles detachment process sure (H = 1 m) with 200 × 200 grid resolution where the notched
is well reproduced by the enhanced VOF method. disk with nominal radius of R = 0.15m is initially centered at
(0.5, 0.75 ). The width and the length of the slot are set equal
4.4. Rotation of a slotted disk (Zalesak’s problem) to W = 0.05m and L = 0.25m similar to works of Owkes et al.
[105] and Ming-Jian Li [106]. The notched disk is subjected to a
To further demonstrate the accuracy and robustness of the counter-clockwise rotating eddy given by:
proposed third-order bounded convection scheme in suppressing
false-diffusion error and maintaining the sharpness of the inter-
face, the solid-body rotation case of Zalesak [104] (case 4) is con- u = −2π (y − 0.5 )
(38)
sidered here as a fourth benchmark test problem. As schematically u = +2π (x − 0.5 )
17
Fig. 12. Transient evolution of deformation of a 2D disk (case 5) in a vortex velocity field predicted by the proposed model at different time instants.
The transient evolution of volume-fraction field at various time tial configuration and outperformed four other existing Eulerian
instants is depicted in Fig. 11. It is evident that, the general shape models tested by Cifani et al. [63] and Ming-Jian Li [106]. Finally,
of the slotted circle is successfully preserved after one full ro- from Table 2 one can deduce that, the predicted results in the cur-
tation. The qualitative and quantitative comparisons of the ob- rent work is more consistent with Lagrangian methods suggested
tained results with previously published data in terms of the fi- by Scardovelli et al. [93].
nal shape of the slotted disk and accumulated error (EZalesak =
N
N 4.5. Deformation of a 2D disk (case 5)
|ϕi, j − ϕ i, j |/ ϕ i, j ) are presented in Fig. 11 and Table 2.
i=1, j=1 i=1, j=1
The transient deformation of a weightless liquid circular disk
It can be seen that, despite a little smoothing on the sharp edges,
with nominal radius of R = 0.15m centered at (0.5, 0.75 ) within
the proposed improved VOF model can accurately maintain the ini-
a unit square domain (H = 1 m) is examined in this subsec-
18
Fig. 13. Qualitative comparison of obtained results for case 5 with numerical data of Owkes et al. [105] Li et al. [106].
Table 2
Zalesak’s problem (case 4): comparison of accumulated error.

N
N
Advection algorithms Error (E = |ϕi, j − ϕ i, j |/ ϕ i, j )
i=1, j=1 i=1, j=1
Present model/Eulerian 5.72 × 10−3

gVoFoam/Eulerian [63] 1.36 × 10−2
interFoam/Eulerian [63] 6.61 × 10−2
OpenFoam®/Eulerian [106] 7.03 × 10−2
Puckett/Stream/Eulerian [93] 1.00 × 10−2
Quadratic fit/Lagrangian [93] 5.47 × 10−3
Quadratic fit+continuity/Lagrangian [93] 4.16 × 10−3
tion as fifth canonical test case where due to the presence [105]:
of non-uniform vorticity field (shearing flow), the initially cir-
u = −2sin (π x ) sin(π y ) cos(π y ) cos( πTt )
2
cular disk begins to evolve into a spiral structure with max- (39)
imum distortion. The simulation is carried out on 256 × 256 v = +2sin2 (π y ) sin(π x ) cos(π x ) cos( πTt )
uniform grid system and periodic boundary condition is im- where u and v denote velocity components in x and y-directions.
posed on all enclosure walls. As sketched in Fig. 1, the com- The term t is time and T=8.0 represents the period of the veloc-
putational domain is subjected to a single vortex described by ity field. The sequential evolution of rotating disk at various time
19
Fig. 14. Simulation of the pure convecttion of a step profile (case 6) using three different convection schemes (a) the newely developed third-order bounded convection
scheme (Eq. (25)), (b) Minmod scheme and (c) the classical first-order Upwind scheme.
instants is portrayed in Fig. 12. It can be observed that contrary Table 3

Deformation of a 2D disk problem (case 5): comparison of accumulated
to previous test case (Zalesak’s disk), due to the imposition of the
error.
solenoidal velocity field, the circular fluid body is progressively
stretched and dragged towards the centre of the recirculating zone,
N
N
Advection algorithms Error (E = (ϕi, j − ϕ i, j )/ ϕ i, j )
giving rise to the development of the thin continuous spiral- i=1, j=1 i=1, j=1
shaped filament within the enclosure at t = T /2 = 4s. However, Present model/Eulerian 6.13 × 10−4
owing to the existence of cosinusoidal temporal factor (cos(π t /T )), gVoFoam/Eulerian [63] 2.95 × 10−3
interFoam/Eulerian [63] 5.58 × 10−2
the fluid flow starts moving backwards in a counterclockwise di-
rection (T /2 < t) until the elongated disk returns to its equilibrium
(initial) position at T=8.0. The comparison of the predicted results
with numerical works of Cifani et al. [63], Owkes et al. [105] and reported by Zhang et al. [107] who investigated the performance,
Ming-Jian Li [106] in Fig. 13 and Table 3 illustrates that, in spite accuracy and stability of nine different composite high-resolution
of severe topological change of rotating disk during its evolution, convective schemes in dealing with the multiphase flows with ma-
mass conservation and the thickness of the thin filament are ef- terial discontinuities.
fectively preserved by the proposed model. Furthermore, Fig. 13 il-
lustrates that, in contrast to the well-known Upwind and Minmod
4.6. Pure convection of a step profile (case 6)
schemes, the newly-proposed convection scheme (Eq. (25)) vividly
demonstrates excellent overall performance in terms of accuracy
The superiority and capability of the proposed non-oscillatory
and eliminating false diffusion errors. Similar observation was also
TVD flux-limiter scheme over the five classical monotone convec-
20
case 7, to induce the buoyancy force and Rayleigh-Bénard Convec-

tion (RBC) [108] in case 8, the horizontal walls of the enclosure
are kept at different constant temperatures (Tc < Th ) while the side
walls are maintained adiabatic. The spaces between the two con-
centric squares in both cavities are filled with air (Pr = 0.71) as a
working fluid and the thermal conductivity ratio of solid to fluid
is taken as Kr = 0.2 similar to work of Das et al. [109]. Based on
the grid independency test provided in Table 4, the simulations
are carried on 159 × 159 uniform mesh resolution and obtained re-
sults in the forms of streamlines, isotherms, local and average Be-
jan number, entropy generation due to heat transfer irreversibility
(ST ) and viscous dissipations (SF ) are presented in Figs. 16 to 21 at
various Rayleigh numbers (104 ≤ Ra ≤ 107 ).
Generally, due to the existence of temperature differences
Fig. 15. Comparison of ϕ profiles along a vertical line (x = 0.7m) of the domain
within the enclosures, the flow patterns and temperature distribu-
at steady state condition predicted by the proposed flux-limiter scheme with five
existing TVD convection schemes [107]. tions are characterized by upward movement of the heated fluid
and downward motion of the cooled one, which lead to the ther-
mal energy exchange between the two active walls. Fig. 16 re-
tion schemes in minimizing false-diffusion errors are further ver- veals that, in case 7, the buoyancy force induced by differentially
ified via simulation of the advective transport of a step profile in heated/cooled sidewalls causes the fluid to move upward along the
an oblique velocity field [107]. As shown in Fig. 1, the calculations left wall until it encounters the insulated top wall. After impinging
are performed in a square cavity where Dirichlet boundary condi- to the top wall, the expanded fluid changes its direction and travels
tions are applied to the left (ϕ = 1.0) and bottom (ϕ = 0) sides of horizontally towards the right wall as its internal energy decreases.
the domain while homogeneous Neumann boundary conditions are As the fluid approaches the downstream wall, because of its inter-
imposed on the right (∂ϕ /∂ x = 0) and top (∂ϕ /∂ y = 0) walls. The action with surrounding cold fluid, its temperature decreases and
computational domain is subjected to the uniform velocity field consequently the fluid becomes gradually colder and denser. The
(u = v = 1.0), oblique to the grid system at an angle of 45 degree relatively condensed fluid eventually reaches the boundary layer of
with respect to the horizontal direction. The governing conserva- the cold wall and starts to flow down and exchange its high level
tion equation of the problem is: of energy. Since the bottom wall is sealed, the heavy fluid moves
obliquely outward towards the hot area to complete its thermal
∂ϕ ∂ uϕ ∂vϕ
+ + + ∇ · (ϕ (1 − ϕ )uR ) = 0 (40) cycle. Hence, the clockwise-rotating eddy is established within the
∂t ∂x ∂y annulus. This process is repeated until no significant improvement
where ϕ represents the dependent variable and u and v are is observed in the rate of heat transfer (steady-state condition). It
the velocity components in the Cartesian framework (see also can be seen that at Ra = 104 , isotherms in case 7 are almost par-
Eq. (5)). To be consistent with literature [90,107], a uniformly allel to each other and evenly distributed which indicate that the
spaced grid size of 120 × 120 is utilized for the simulations and flow and temperature are in the weak coupling regime and heat
the first-order implicit Euler scheme is employed for the time transfer is governed by conduction mode. However, by increasing
discretization. The temporal evolution of volume fraction field the Rayleigh number up to 105 , the role of the buoyancy force be-
using three different schemes is displayed in Fig. 14. As ex- comes more pronounced and subsequently isotherms start to fol-
pected, the substantial improvement has been obtained by im- low the pattern of streamlines, signifying the onset of natural con-
plementing the proposed third-order TVD scheme whereas the vection at this stage. In addition, it is seen that due to the rela-
first-order bounded schemes (i.e. Minmod and Upwind) are char- tive dominance of convection mode, both isotherms and stream-
acterized by severe numerical diffusion which can be attributed lines are more clustered and compressed in the vicinity of the ac-
to the larger truncation errors generated by these schemes. This tive walls and two secondary vortices are appeared in both sides of
finding is also in accordance with the conclusions drawn from the conductive obstacle which manifests the presence of adverse
the previous benchmark problem and those reported by Dar- pressure gradients in those regions. With further increase in the
wish et al. [90]. The comparison of the obtained results with value of Rayleigh number (Ra = 106 ), the clustering structure be-
five other conventional TVD schemes in term of the computa- comes more prominent and distinct thermal boundary layers are
tional ϕ values along a vertical line of x = 0.7m at steady state established in close proximity of the vertical walls. In this circum-
condition in Fig. 15 once again confirms the superiority and ro- stance, due to substantial coupling of the heat and fluid flows,
bustness of the current non-oscillatory TVD flux-limiter scheme the isotherms and cores of the two secondary vortices get dis-
over the existing monotone convection schemes in capturing steep torted and shifted towards the sidewalls where high levels of ve-
gradients. locity gradients occur. The horizontally oriented isotherms across
the block in Fig. 16 confirm that, steep temperature gradients exist
4.7. Conjugate natural convection-conduction heat transfer (cases 7 in the vertical direction, signifying the high level of thermal mixing
and 8) within the enclosure. The absolute value of maximum stream func-
tion (|ψmax |) in Table 5 reveals that, further increment of Rayleigh
The transient entropy generation due to buoyancy-driven conju- number (Ra = 107 ) leads to remarkable intensification of convec-
gate thermal convection in two different horizontal annular ducts tive flow and heat exchange inside the enclosure. It can be seen
with a centered internal conducting obstacle (L = 0.3H) is ana- from Fig. 16 that this stage is dynamically characterized by the de-
lyzed in this subsection. As schematically shown in Fig. 1, the first velopment of small localized eddies on the corners of the cavity
benchmark problem (case 7) is strictly analogous to the classical and significant reduction in the thermal boundary layers thickness
Differentially Heated Cavity (DHC) [17] where gravity force is per- which support the dominance of the convective regime.
pendicular to the temperature gradient. In this case, the isother- However, by rotating the differentially heated vertical cavity
mal heating and cooling are imposed on the vertical walls whereas (case 7) by 90 degree in clockwise direction, the RBC configuration
the horizontal walls are thermally insulated. However, contrary to (case 8) will be established where gravity force becomes parallel
21
Fig. 16. Transient variation of isotherms and streamlines at different Rayleigh numbers for case 7. Pr = 0.71.
Table 4
Variations of Nu as a function of the grid size in cases 7 and 8 at low and high Rayleigh num-
bers.
Number of grids (Case 7)
Ra 39×39 59×59 79×79 99×99 119×119 139×139 159×159
104 2.1675 2.1832 2.1966 2.2051 2.2134 2.2168 2.2177

107 17.1438 17.1944 17.2132 17.2225 17.2286 17.2307 17.2319
Number of grids (Case 8)
Ra 39×39 59×59 79×79 99×99 119×119 139×139 159×159
104 2.2505 2.2601 2.2692 2.2784 2.2825 2.2846 2.2853

107 11.1563 11.3219 11.4134 11.4308 11.4392 11.4406 11.4417
22
Fig. 17. Transient variation of isotherms and streamlines at different Rayleigh numbers in case 8. Pr = 0.71.
to an imposed horizontal temperature gradient. By examining the RaCritical with RaCritical = 1708 being the critical Rayleigh number),
thermal structure and dynamics of RBC in Fig. 17 one can deduce (b) stable-convection-dominated regime where viscous force is
that, at the same nominal values of Rayleigh and Prandtl numbers, strong enough to overcome the oscillatory behavior of RBC insta-
RBC exhibits fundamentally different heat and momentum trans- bility where the CW eddy together with two secondary CCW vor-
port characteristics compared to DHC during its evolution. In gen- tices are developed within the enclosure (RaCritical < Ra ≤ 5 × 106 )
eral, according to linear stability theory [18], hydrodynamics and [108] and (c) unstable-convection-dominated regime where tran-
thermal behaviors of Rayleigh-Bénard Convection can be classified sition from laminar to turbulent flow occurs (5 × 106 < Ra ≤ 107 )
into the three distinct regimes, namely: (a) diffusion-dominated [110]. It should be noted that in the last stage, owing to the in-
regime where four symmetric CW and CCW eddies are formed herent unsteady nature of the RBC, the streamlines and isotherms
within the enclosure and isotherms are horizontally oriented (Ra ≤ constantly change their topology such that oscillatory and non-
23
Fig. 18. Transient variation of local Bejan number at different Rayleigh numbers for cases 7 and 8. Pr = 0.71.
stationary thermal plume structures are cyclically emerged and ductive obstacle. However, since the Ra is greater than the critical
disappeared inside the enclosure. Based on the above discussions value, the vertical symmetry feature ultimately dies out and the
and information from the literature, the thermal hydraulic behav- clockwise recirculating eddy is created inside the enclosure. More-
ior of RBC portrayed in Fig. 17 can be qualitatively analyzed as fol- over, isotherms are vertically stratified and weak thermal gradient
lows. is formed inside the obstacle which in turn supports the incep-
Starting from a cold isothermal initial state at Ra = 104 , it can tion of the convection at this Ra. As shown in Fig. 17, the early
be seen that the process consists of the development of the as- stages of flow evolution at Ra=105 are qualitatively similar to the
cending thermal plume in the close proximity to the hot bound- transient variations of temperature and velocity fields at Ra=104
ary which causes the fluid to move upward along the centerline where flow bifurcation occurs. As time elapses, the buoyancy force
of the enclosure and eventually impact on the bottom of the con- overweighs the retarding effects of the viscous force, leading to the
24
Fig. 19. Transient variation of local entropy generation due to heat transfer irriversibility (ST ) at different Rayleigh numbers for cases 7 and 8. Pr = 0.71.
formation of the strong CW rotating cell and two distinct thermal of the enclosure and (b) two secondary counter-clockwise eddies
boundary layers near the horizontal walls. Similar to the previous which are spatially confined to the smaller regions near the top-
case, the significant twisting of the isotherms and streamlines are left and bottom-right corners of the enclosure where isotherms are
interpreted to reflect improvements in the convection heat trans- largely clustered and characterized by the heat entrapment. It is
fer within the annulus. With increasing Rayleigh number (Ra=106 ), worth recalling that, similar to the work of Ouertatani et al. [108],
due to high buoyancy-to-inertia ratio, a pair of ascending thermal the flow structure and temperature pattern together with the two
plume is generated in the left and right halves of the obstacle. minor convective cells are completely symmetric with respect to
The formation of these swaying and swirling structures of ther- the center of the cavity. It is worth noting that, such a change
mal plume is also responsible for the development of some sec- of flow structure is known as the symmetry-breaking bifurcation
ondary vortices within the enclosure. As time progressed over the which typically occurs in the transient natural convection between
period of 0.0187 ≤ t ≤ 0.0525, two separate plumes begin to coa- two parallel plates with a heated horizontal base and cooled upper
lesce, creating two distinct recirculation zones, namely: (a) the pri- walls [111,112]. It is evident that as the time increases (0.0525 <
mary clockwise vortex (Bénard cell) which occupies the major part t), the primary cell grows in size and squeezes two CCW vor-
25
Fig. 20. Transient variations of average Nusselt number (Nu), total entropy generation (Stot ), average Bejan number (Be), entropy generation due to thermal (ST ) and viscous
(SF ) dissipations as function of the dimentional time at different Rayleigh numbers (104 ≤ Ra ≤ 107 ) for cases 7 and 8. (Pr = 0.71, L = 0.3H, Kr = 0.2).
tices towards the corners of the enclosure until the steady-state tions and structural alterations. It is evident that, at the very be-
condition is achieved. As expected, with further increase of the ginning of the transient evolution (0 ≤ t ≤ 0.0059), three distinct
Rayleigh number (Ra=107 ), the strength of the convective flow thermal plumes and unstable multi-cellular structure are formed
augments and subsequently RBC enters into the unsteady convec- within the enclosure. As time goes on (0.0059 < t ≤ 0.0637), the
tion regime where flow and thermal fields undergo severe distor- buoyancy force begins to disrupt the symmetrical pattern of the
26
Fig. 21. Zoomed-in view of average Nusselt number with associated snapshots of isotherms and streamlines at Ra=107 for case 8 (conjugate Rayleigh-Bénard convection).
vortices, leading to the initiation of convective instability and the convection effects, the local Bejan number values in most parts
occurrence of vortex breakdown, accordingly. of the enclosure in both cavities are close to unity which implies
The time-evolution of the local Bejan number and entropy gen- that most of the exergy destruction takes place as a results of the
eration due to heat transfer irreversibility associated with cases 7 heat transfer irreversibility (SF < ST ). The corresponded local en-
and 8 at various Rayleigh number are displayed in Figs. 18 and tropy generation due to thermal dissipation in Fig. 19 shows that
19. A detailed examination of Fig. 20 reveals that, at low Rayleigh the active zones of ST in case 7 are more confined to the small
number (Ra=104 ) due to diffusion-dominated regime and weak regions near the bottom-left and top-right corners of the enclo-
27
Fig. 22. Comparison of predicted results with the numerical data of Das et al. [109] for the case of Differentially Heated Cavity with an internal conductive body. Kr =0.2,
Pr=0.71 and L=0.5H.
sure where maximum heat transfer rate occurs. In contrast, ST tropy generation due to both factors (ST , SF ) augment whereas a
is almost trivial in middle portions of the adiabatic walls where reverse trend is found in the values of average Bejan number. This
isotherms are smooth curves and uniformly scattered due to in- behavior can be attributed to the fact that, by increasing the Ra,
significant heat flow in those areas. It is interesting to note that, the thickness of the thermal boundary layer decreases and inten-
this course of the event also takes place in case 8 but in re- sity of the convection flow and recirculation patterns augment. The
vers manner such that due to considerable amount of heat flow reduction in the thermal boundary layer thickness and augmen-
along the horizontal walls, heat transfer irreversibility (ST ) is found tation in the compression of temperature profiles are interpreted
to be maxima in the bottom-left and top-right junctions of the to reflect the intensification of temperature gradient, heat transfer
isothermal-adiabatic walls whereas vertical walls have negligible rate and ST . On the other hand, an increase in the buoyancy force
contributions to ST . Table 5 illustrates that, as the Ra enhances, en- promotes the fluid motion and elevates velocity gradients, leading
28
Table 5
The effects of the Rayleigh number on the flow intensity, average Nusselt number and entropy generation rate for cases 7 and 8
(steady-state condition). Note that, case 8 never reaches a steady state condition at Ra=107 .
Case 7
Ra Nu |Umax | |ψmax | Be Stot Stot,max ST ST,max SF SF,max

4
10 2.217 15.217 3.459 0.618 3.240 25.463 2.185 14.503 1.054 14.914
105 4.624 40.902 8.849 0.299 23.061 629.852 4.523 63.850 18.538 590.402
106 9.023 120.314 16.537 0.188 332.456 16362.64 8.687 325.676 323.768 16227.41
107 17.231 350.577 30.429 0.163 5465.67 373316.035 16.027 1553.121 5449.64 372930.16
Case 8
Ra Nu |Umax | |ψmax | Be Stot Stot,max ST ST,max SF SF,max

4
10 2.285 19.369 4.276 0.579 3.543 17.587 2.257 11.743 1.286 16.695
105 3.853 77.987 17.062 0.265 30.091 471.496 3.740 34.435 26.351 471.378
106 6.595 273.305 52.611 0.106 473.863 15582.646 6.027 147.585 467.836 15582.543
107 11.441 785.986 105.237 0.088 6826.085 329040.76 9.196 398.455 6816.88 329040.51
to significant enhancement in the fluid friction irreversibility ac- For getting further insight into the effects of Rayleigh number
cording to Eq. (16). However, since, entropy generation due to fluid on the performance of thermal systems in cases 7 and 8, time his-
friction irreversibility (SF ) grows much faster compared to the ther- tories of the average Nusselt number (Nu), total entropy genera-
mal dissipation (ST ), the local and average Bejan number decline tion (Stot ) due to both factors (ST , SF ) and average Bejan number
rapidly and become less than half (Be < 0.5) which indicates that (Be) are depicted in Fig. 20. As expected for all cavities, both heat
the dominant source of irreversibility is shifted toward the vis- transfer rate (Nu) and overall entropy generation augment with in-
cous dissipation. This trend is continued or even accentuated as creasing Rayleigh number while a reverse trend is observed in the
the Ra increases. This assertion is well borne out by Fig. 18 where values of Be. As stated before, an increase of Ra improves buoy-
due to the existence of flow separation and weak fluid circulation ancy force and consequently the heat flow effects, hence, the val-
in the left and right portions of the conductive obstacle in case ues of both ST and SF increase. However, since fluid fiction irre-
7, the fluid movement is restrained and thereby the local Bejan versibility is more sensitive than the thermal dissipation to any al-
number becomes greater than half in those areas. However, the teration in Ra, the average Bejan number drops immensely. It can
situation gets worst in case 8. As deduced from the magnitude be seen that, at Ra = 104 , the value of the Be is greater than half,
of streamlines in Table 5 (|ψmax |), because of high intense fluid which indicates that total entropy generation (Stot ) is more influ-
circulation cells, the entire enclosure in Rayleigh-Bénard Convec- enced by ST rather than viscous dissipation (Stot ≈ ST ). This conjec-
tion (RBC) acts as an active site of the entropy generation due to ture is well reflected in Fig. 20 where the variations of Stot have
friction irreversibility (compare also the values of |ψmax | between the same pattern of alteration with ST . As Rayleigh number ap-
cases 7 and 8 in Table 5). Further inspection of Table 5 also re- proaches the flow-transition value (Ra = 105 ), the fluid flow be-
veals that the optimal situation based on the highest heat trans- comes stronger and consequently the contribution of SF to Stot is
fer rate and lowest entropy generation also takes place in case 7. invigorated dramatically so that total entropy generation begins
Take as example, the ratios of average Nusselt number over that to follow the transient behavior of the fluid friction irreversibil-
of total entropy generation (Nu/Stot ) in case 7 at Ra = 104 , 105 , 106 ity (Stot ≈ SF ). This enhancement mechanism also manifests itself
and 107 are equal to 0.68, 0.2, 0.027 and 0.0031 while these ratios in form of the decreasing trend in the profiles of average Bejan
for case 8 are approximately 0.64, 0.12, 0.014 and 0.0016, respec- number in both cavities. It can be seen that the degradation in the
tively. In fact, these outcomes convey that Differentially Heated values of Be is accelerated and intensified as the flow enters into
Cavity (DHC) outperformed the RBC from the viewpoints of opti- the convection-dominated regime (Ra ≥ 106 ).
mal design, energy saving and minimization of entropy generation. However, before closing this section, there are some other in-
Similar observations were also reported by Anandalakshmi et al. teresting points that can be drawn from Fig. 20. It is apparent
[113] who examined the effects of the Ra on the entropy gener- that, at low Rayleigh number (Ra = 104 ), the aforementioned per-
ation during free convection heat transfer in rhombic enclosures tinent parameters tend to reach the constant value with the small-
with various inclination angles. Furthermore, it is apparent from amplitude fluctuation around the state of equilibrium while, at
Figs. 18 and 19 that, similar to velocity and temperature fields, high Rayleigh numbers (particularly beyond the critical value of
the contours plots of local Be and ST in both cases are symmet- RaCritical = 5200), the oscillatory behavior starts to emerge in their
rical with respect to the center of the enclosure which can be re- profiles [29]. This physical behavior can be attributed to the strong
garded as an extra validation from the standpoints of energy bal- interaction between hot and cold fluids at high Ra which causes
ance and mass conservation. By comparing Figs. 18 and 19 one can the induction of high-frequency internal waves and hydraulic jump
conclude that, although the maximum entropy generation due to in the velocity and temperature fields. These observations are in
heat transfer (ST,max ) still occurs near the heated surfaces where accordance with numerical findings of Magherbi et al. [29] and
the isotherms are largely compressed, but due to imposition of Schladow [114] who discussed that for RaCritical < Ra, when the ini-
the no-slip boundary conditions and enhanced convection effect tial condition is somewhat far from the equilibrium state, fluid
at Ra=106 , the significant velocity gradients are also generated in flow and heat transfer are likely to undergo nonlinear oscillation
close vicinity of the rigid walls (ST SF ), leading to the notable around the steady-state condition. Another interesting feature is
reduction in the values of local Bejan number in those areas (Be < that in both cavities due to assigning zero-temperature initial con-
0.5). Finally, at Ra=107 due to high velocity gradients generated dition (θ = 0), the average Nusselt number along the hot and cold
by the buoyancy force, viscous dissipation increases immensely walls decreases and increases asymptotically and reaches the iden-
and covers the more part of the domain which implies that at tical value which implies that the amount of the heat released
this stage the substantial amount of available work (exergy) might from the hot surface is totally absorbed by the cold wall which
be utilized to overcome the energy losses induced by frictional in turn demonstrates the accuracy of the proposed model from the
irreversibility. first law of thermodynamics viewpoint. Finally, as pointed up ear-
29
lier, by increasing the Rayleigh number from 106 to 107 , the flow • It was shown that, the enhanced VOF model can robustly han-
pattern in the Rayleigh-Bénard Convection (RBC) is evolved directly dle large topological changes of complex multiphase flows in-
from the unicellular steady regime to the multi-cellular oscillatory volving interface rupture and coalescence.
solution. This mechanism is also accompanied by the formation of • It was mathematically shown that, the ELVIRA technique can
a symmetry/asymmetry Hopf bifurcation structure (Pitchfork bifur- exactly reproduce any arbitrary linear interfaces.
cation) [115,116] and vortex breakdown which leads to the striking • The results revealed that, the improved VOF model outper-
fluctuations in the heat transfer rate within the enclosure. To bring formed the existent conventional VOF method available in the
this point out more clearly, the enlarged view of time-distributions OpenFoam® platform.
of Nu with associated streamlines and isotherms snapshots are re- • The results revealed that, from the Entropy Generation Mini-
plotted in Fig. 21. The weakly turbulent flow at Ra=107 manifests mization (EGM) and optimal design viewpoints, the Differen-
itself through the emergence of semi-periodic thermal plumes and tially Heated Cavity (DHC) outperformed the Rayleigh-Bénard
unsteady minor vortices. It is evident that during the advection Convection (RBC) in terms of the highest heat transfer rate and
process (0 ≤ t ≤ 0.8), the direction and position of the flow cir- lowest rate of entropy generation.
culation and travelling wave are constantly changing (and inter-
nally displaced towards the downstream and upstream end walls) Declaration of Competing Interest
within the enclosure, meaning that the internal waves and pertur-
bations generated by the initial conditions cannot be suppressed The authors declare that they have no known competing finan-
by the boundary layer. Figure 22 shows the comparison of the pre- cial interests or personal relationships that could have appeared to
dicted results with numerical outcomes of Das et al. [109], for the influence the work reported in this paper.
case of Differentially Heated Cavity with internal conductive obsta- The authors declare the following financial interests/personal
cle (Kr = 0.2 and L = 0.5H). It is apparent that the obtained results relationships which may be considered as potential competing in-
in terms of streamlines, isotherms and average Nusselt number are terests:
in excellent agreement with previously published data.
Acknowledgments
5. Conclusion This work is financially supported by the Merit scholarship pro-

gram for foreign students (PBEEE) of the Ministère de l’Éducation
In this study, a modified version of the Volume-Of-Fluid model et de l’Enseignement supérieur du Québec (MEES) (Grant No.
was introduced for simulation of multiphase and free-surface flows 263222), (Programme V1-Groupe: G2-20431).
with high density contrast. To circumvent the difficulties regarding
the physical discontinuities at the material interface, a novel third-
Supplementary materials
order bounded flux-limiter function was first constructed based on
the TVD constraints and then employed for discretization of the
Supplementary material associated with this article can be
convection terms in the momentum, energy and volume fraction
found, in the online version, at doi:10.1016/j.ijheatmasstransfer.
equations. To further enhance the performance of the proposed
2021.121163.
model, the classical PISO algorithm was first modified based on
the Two-step projection (PISOC) method and then was used for the
References
treatment of the velocity-pressure coupling. To deal with the sur-
face tension at the interface, the procedure for implementing the [1] N. Zhao, J. Yang, H. Li, Z. Zhang, S. Li, Numerical investigations of laminar heat
ELVIRA technique was first described in detail and then adopted transfer and flow performance of Al2O3-water nanofluids in a flat tube, Int.
to calculate the interface normal vector and its corresponding cur- J. Heat Mass Transf. 92 (2016) 268–282, doi:10.1016/j.ijheatmasstransfer.2015.
08.098.
vature. The accuracy and robustness of the improved VOF model [2] M. Chakkingal, S. Kenjereš, I. Ataei-Dadavi, M.J. Tummers, C.R. Kleijn, Numer-
were verified against eight different canonical test cases including: ical analysis of natural convection with conjugate heat transfer in coarse–
dam break, Rayleigh-Taylor instability, bubble rising, rotation of a grained porous media, Int. J. Heat Fluid Flow. 77 (2019) 48–60.
[3] A.I. Alsabery, M.A. Sheremet, A.J. Chamkha, I. Hashim, Conjugate natural con-
slotted disk (Zalesak’s problem), deformation of a 2D disk, pure
vection of Al2O3–water nanofluid in a square cavity with a concentric solid
convection of a step profile and transient entropy generation due insert using Buongiorno’s two-phase model, Int. J. Mech. Sci. 136 (2018)
to conjugate natural convection heat transfer in square enclosures. 200–219.
[4] M.J. Uddin, O.A. Bég, M.N. Uddin, Energy conversion under conjugate conduc-
Based on the obtained results, the main conclusions of present
tion, magneto-convection, diffusion and nonlinear radiation over a non-lin-
study can be summarized as follows: early stretching sheet with slip and multiple convective boundary conditions,
Energy 115 (2016) 1119–1129.
• It was found that, the newly proposed bounded convection [5] S. Pandey, Y.G. Park, M.Y. Ha, An exhaustive review of studies on natural con-
vection in enclosures with and without internal bodies of various shapes, Int.
scheme can successfully suppress false diffusion error and J. Heat Mass Transf. 138 (2019) 762–795.
maintain sharpness of the interface. [6] A. Amine, J.K. Platten, M. Hasnaoui, Thermal convection around obstacles: the
• The results revealed that, the proposed non-linear flux-limiter case of Sierpinski carpets, Exp. Fluids. 36 (2004) 717–727.
[7] D.-D. Zhang, J.-H. Zhang, D. Liu, F.-Y. Zhao, H.-Q. Wang, X.-H. Li, Inverse conju-
function can successfully control and maintain the thickness of gate heat conduction and natural convection inside an enclosure with multi-
the interface to a maximum of 2-3 grid stencils. ple unknown wall heating fluxes, Int. J. Heat Mass Transf. 96 (2016) 312–329.
• It was found that, the combination of the two classical PISO and [8] Q. Ren, C.L. Chan, Natural convection with an array of solid obstacles in an
enclosure by lattice Boltzmann method on a CUDA computation platform, Int.
Two-step projection method (PISOC) ensures that the mixed ve-
J. Heat Mass Transf. 93 (2016) 273–285.
locity in the vicinity of the material discontinuity is free of spu- [9] A. Raji, M. Hasnaoui, M. Naïmi, K. Slimani, M.T. Ouazzani, Effect of the sub-
rious spatial oscillations. division of an obstacle on the natural convection heat transfer in a square
cavity, Comput. Fluids. 68 (2012) 1–15.
• It was found that, the proposed PISOC algorithm can effectively
[10] F. Garoosi, F. Talebi, Numerical simulation of conjugate conduction and natu-
accelerate the convergence rate of numerical solution without ral convection heat transfer of nanofluid inside a square enclosure containing
introducing unphysical pressure fluctuations. a conductive partition and several disconnected conducting solid blocks using
• It was found that the newly developed semi-iterative hybrid the Buongiorno’s two phase model, Powder Technol 317 (2017) 48–71.
[11] F. Garoosi, M.M. Rashidi, Two phase flow simulation of conjugate natural con-
PISOC algorithm is robust and consistent with superior conver- vection of the nanofluid in a partitioned heat exchanger containing several
gence characteristics to the standard iterative PISO algorithm. conducting obstacles, Int. J. Mech. Sci. 130 (2017) 282–306.
30
[12] J.-T. Hu, X.-H. Ren, D. Liu, F.-Y. Zhao, H.-Q. Wang, Conjugate natural convection [41] Z. Wang, R. Chen, X. Zhu, Q. Liao, D. Ye, B. Zhang, X. He, L. Jiao, Dynamic be-
inside a vertical enclosure with solid obstacles of unique volume and multiple haviors of the coalescence between two droplets with different temperatures
morphologies, Int. J. Heat Mass Transf. 95 (2016) 1096–1114. simulated by the VOF method, Appl. Therm. Eng. 131 (2018) 132–140.
[13] I.V Miroshnichenko, M.A. Sheremet, Radiation effect on conjugate turbulent [42] F. Yeganehdoust, R. Attarzadeh, I. Karimfazli, A. Dolatabadi, A numerical anal-
natural convection in a cavity with a discrete heater, Appl. Math. Comput. ysis of air entrapment during droplet impact on an immiscible liquid film,
321 (2018) 358–371. Int. J. Multiph. Flow. 124 (2020) 103175.
[14] F. Selimefendigil, H.F. Öztop, Conjugate natural convection in a nanofluid [43] G. Tryggvason, R. Scardovelli, S. Zaleski, Direct numerical simulations of
filled partitioned horizontal annulus formed by two isothermal cylinder sur- gas–liquid multiphase flows, Cambridge University Press, 2011.
faces under magnetic field, Int. J. Heat Mass Transf. 108 (2017) 156–171. [44] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free
[15] L. Wang, Y. Zhao, X. Yang, B. Shi, Z. Chai, A lattice Boltzmann analysis of the boundaries, J. Comput. Phys. 39 (1981) 201–225, doi:10.1016/0021-9991(81)
conjugate natural convection in a square enclosure with a circular cylinder, 90145-5.
Appl. Math. Model. 71 (2019) 31–44. [45] M. Sussman, P. Smereka, S. Osher, A level set approach for computing
[16] H. Zargartalebi, M. Ghalambaz, K. Khanafer, I. Pop, Unsteady conjugate natural solutions to incompressible two-phase flow, J. Comput. Phys. 114 (1994)
convection in a porous cavity boarded by two vertical finite thickness walls, 146–159.
Int. Commun. Heat Mass Transf. 81 (2017) 218–228. [46] G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber,
[17] G. De Vahl Davis, Natural convection of air in a square cavity: A bench mark J. Han, S. Nas, Y.-J. Jan, A front-tracking method for the computations of mul-
numerical solution, Int. J. Numer. Methods Fluids. 3 (1983) 249–264, doi:10. tiphase flow, J. Comput. Phys. 169 (2001) 708–759.
10 02/fld.1650 030305. [47] M. Ishii, T. Hibiki, Thermo-fluid dynamics of two-phase flow, Springer Science
[18] P.H. Kao, R.J. Yang, Simulating oscillatory flows in Rayleigh-Bénard convection & Business Media, 2010.
using the lattice Boltzmann method, Int. J. Heat Mass Transf. 50 (2007) 3315– [48] I. Malgarinos, N. Nikolopoulos, M. Gavaises, Coupling a local adaptive grid re-
3328, doi:10.1016/j.ijheatmasstransfer.2007.01.035. finement technique with an interface sharpening scheme for the simulation
[19] K. Kalidasan, R. Velkennedy, P.R. Kanna, Laminar natural convection of Cop- of two-phase flow and free-surface flows using VOF methodology, J. Comput.
per-Titania/Water hybrid nanofluid in an open ended C-shaped enclosure Phys. 300 (2015) 732–753.
with an isothermal block, J. Mol. Liq. 246 (2017) 251–258. [49] A. Issakhov, Y. Zhandaulet, A. Nogaeva, Numerical simulation of dam break
[20] A. Rahimi, A.D. Saee, A. Baghban, A. Kasaeipoor, H. Ashrafi, E.H. Malek- flow for various forms of the obstacle by VOF method, Int. J. Multiph. Flow.
shah, Double-MRT lattice Boltzmann simulation of natural convection in a 109 (2018) 191–206.
C-shaped heat exchanger, Powder Technol 336 (2018) 465–480. [50] A. Issakhov, M. Imanberdiyeva, Numerical simulation of the movement of wa-
[21] N. Makulati, A. Kasaeipoor, M.M. Rashidi, Numerical study of natural convec- ter surface of dam break flow by VOF methods for various obstacles, Int. J.
tion of a water–alumina nanofluid in inclined C-shaped enclosures under the Heat Mass Transf. 136 (2019) 1030–1051.
effect of magnetic field, Adv. Powder Technol. 27 (2016) 661–672. [51] Z.H. Gu, H.L. Wen, C.H. Yu, T.W.H. Sheu, Interface-preserving level set method
[22] B. Mliki, M.A. Abbassi, A. Omri, Z. Belkacem, Lattice Boltzmann analysis of for simulating dam-break flows, J. Comput. Phys. 374 (2018) 249–280.
MHD natural convection of CuO-water nanofluid in inclined C-shaped enclo- [52] B.P. Leonard, Simple high-accuracy resolution program for convective
sures under the effect of nanoparticles Brownian motion, Powder Technol 308 modelling of discontinuities, Int. J. Numer. Methods Fluids. 8 (1988)
(2017) 70–83. 1291–1318.
[23] A. Bejan, Entropy generation through heat and fluid flow, Wiley, 1982. [53] J.A. Heyns, A.G. Malan, T.M. Harms, O.F. Oxtoby, Development of a compres-
[24] A. Bejan, Entropy generation minimization: The new thermodynamics of sive surface capturing formulation for modelling free-surface flow by us-
finite-size devices and finite-time processes, J. Appl. Phys. 79 (1996) ing the volume-of-fluid approach, Int. J. Numer. Methods Fluids. 71 (2013)
1191–1218. 788–804.
[25] A.A.A.A. Al-Rashed, Investigating the effect of alumina nanoparticles on heat [54] M.J. Berger, J. Oliger, Adaptive mesh refinement for hyperbolic partial differ-
transfer and entropy generation inside a square enclosure equipped with two ential equations, J. Comput. Phys. 53 (1984) 484–512.
inclined blades under magnetic field, Int. J. Mech. Sci. 152 (2019) 312–328. [55] S. Hill, D. Deising, T. Acher, H. Klein, D. Bothe, H. Marschall, Boundedness-p-
[26] A.H. Pordanjani, S. Aghakhani, A.A. Alnaqi, M. Afrand, Effect of alumina reserving implicit correction of mesh-induced errors for VOF based heat and
nano-powder on the convection and the entropy generation of water inside mass transfer, J. Comput. Phys. 352 (2018) 285–300.
an inclined square cavity subjected to a magnetic field: Uniform and non-u- [56] F. Denner, B.G.M. van Wachem, Compressive VOF method with skewness cor-
niform temperature boundary conditions, Int. J. Mech. Sci. 152 (2019) 99–117. rection to capture sharp interfaces on arbitrary meshes, J. Comput. Phys. 279
[27] T.A. Alkanhal, M. Sheikholeslami, A. Arabkoohsar, R. ul Haq, A. Shafee, Z. Li, (2014) 127–144.
I. Tlili, Simulation of convection heat transfer of magnetic nanoparticles in- [57] K. Ling, W.-Q. Tao, A sharp-interface model coupling VOSET and IBM for sim-
cluding entropy generation using CVFEM, Int. J. Heat Mass Transf. 136 (2019) ulations on melting and solidification, Comput. Fluids. 178 (2019) 113–131.
146–156, doi:10.1016/j.ijheatmasstransfer.2019.02.095. [58] X. Deng, S. Inaba, B. Xie, K.-M. Shyue, F. Xiao, High fidelity discontinuity-re-
[28] C. Sivaraj, M.A. Sheremet, MHD natural convection and entropy generation of solving reconstruction for compressible multiphase flows with moving inter-
ferrofluids in a cavity with a non-uniformly heated horizontal plate, Int. J. faces, J. Comput. Phys. 371 (2018) 945–966.
Mech. Sci. 149 (2018) 326–337. [59] P. Das, H.S. Udaykumar, A sharp-interface method for the simulation of shock-
[29] M. Magherbi, H. Abbassi, A.Ben Brahim, Entropy generation at the onset of -induced vaporization of droplets, J. Comput. Phys. 405 (2020) 109005.
natural convection, Int. J. Heat Mass Transf. 46 (2003) 3441–3450. [60] M.A. Alves, P.J. Oliveira, F.T. Pinho, A convergent and universally bounded in-
[30] D. Das, M. Roy, T. Basak, Studies on natural convection within enclosures of terpolation scheme for the treatment of advection, Int. J. Numer. Methods
various (non-square) shapes–A review, Int. J. Heat Mass Transf. 106 (2017) Fluids. 41 (2003) 47–75.
356–406. [61] S.K. Choi, H.Y. Nam, M. Cho, A comparison of higher-order bounded convec-
[31] P. Biswal, T. Basak, Entropy generation vs energy efficiency for natural convection schemes, Comput. Methods Appl. Mech. Eng. 121 (1995) 281–301.
tion based energy flow in enclosures and various applications: A review, Re- [62] G. Pozzetti, B. Peters, A multiscale DEM-VOF method for the simulation of
new. Sustain. Energy Rev. 80 (2017) 1412–1457, doi:10.1016/j.rser.2017.04.070. three-phase flows, Int. J. Multiph. Flow. 99 (2018) 186–204.
[32] W. Duan, Y. Gao, Q. Yu, T. Wu, Z. Wang, Numerical simulation of coal gasifi- [63] P. Cifani, W.R. Michalek, G.J.M. Priems, J.G.M. Kuerten, C.W.M. van der Geld,
cation in molten slag: Gas-liquid interaction characteristic, Energy 183 (2019) B.J. Geurts, A comparison between the surface compression method and an
1233–1243. interface reconstruction method for the VOF approach, Comput. Fluids. 136
[33] X. Li, M. Liu, Y. Ma, T. Dong, D. Yao, Experiments and meso-scale modeling of (2016) 421–435.
phase holdups and bubble behavior in gas-liquid-solid mini-fluidized beds, [64] V.-T. Nguyen, W.-G. Park, A volume-of-fluid (VOF) interface-sharpening
Chem. Eng. Sci. 192 (2018) 725–738. method for two-phase incompressible flows, Comput. Fluids. 152 (2017)
[34] A. Tschöpe, M. Wyrwoll, M. Schneider, K. Mandel, M. Franzreb, A Magneti- 104–119.
cally Induced Fluidized-bed Reactor for Intensification of Electrochemical Re- [65] V.-T. Nguyen, V.-D. Thang, W.-G. Park, A novel sharp interface capturing
actions, Chem. Eng. J. (2019) 123845. method for two-and three-phase incompressible flows, Comput. Fluids. 172
[35] T. Wang, T. Tang, Q. Gao, Z. Yuan, Y. He, Experimental and numerical inves- (2018) 147–161.
tigations on the particle behaviours in a bubbling fluidized bed with binary [66] K.K. So, X.Y. Hu, N.A. Adams, Anti-diffusion method for interface steepening
solids, Powder Technol 362 (2020) 436–449. in two-phase incompressible flow, J. Comput. Phys. 230 (2011) 5155–5177.
[36] J. Yan, W. Yan, S. Lin, G.J. Wagner, A fully coupled finite element formulation [67] M. Akhlaghi, V. Mohammadi, N.M. Nouri, M. Taherkhani, M. Karimi, Multi-
for liquid–solid–gas thermo-fluid flow with melting and solidification, Com- -Fluid VoF model assessment to simulate the horizontal air–water intermit-
put. Methods Appl. Mech. Eng. 336 (2018) 444–470. tent flow, Chem. Eng. Res. Des. 152 (2019) 48–59.
[37] J. Zhao, W. Liu, J. Zhao, L. Grekhov, Numerical investigation of gas/liquid [68] J. Hua, D. Mortensen, A front tracking method for simulation of two-phase
two-phase flow in nozzle holes considering the fuel compressibility, Int. J. interfacial flows on adaptive unstructured meshes for complex geometries,
Heat Mass Transf. 147 (2020) 118991. Int. J. Multiph. Flow. 119 (2019) 166–179.
[38] H. Ström, S. Sasic, O. Holm-Christensen, L.J. Shah, Atomizing industrial [69] H. Tan, An adaptive mesh refinement based flow simulation for free-surfaces
gas-liquid flows–development of an efficient hybrid vof-lpt numerical frame- in thermal inkjet technology, Int. J. Multiph. Flow. 82 (2016) 1–16.
work, Int. J. Heat Fluid Flow. 62 (2016) 104–113. [70] L.C. Ngo, H.G. Choi, A multi-level adaptive mesh refinement for an integrated
[39] F. Giussani, F. Piscaglia, G. Saez-Mischlich, J. Hèlie, A three-phase VOF solver finite element/level set formulation to simulate multiphase flows with surface
for the simulation of in-nozzle cavitation effects on liquid atomization, J. tension, Comput. Math. with Appl. (2019).
Comput. Phys. (2019) 109068. [71] O. Antepara, N. Balcázar, J. Rigola, A. Oliva, Numerical study of rising bubbles
[40] J. Huang, X. Zhao, Numerical simulations of atomization and evaporation in with path instability using conservative level-set and adaptive mesh refine-
liquid jet flows, Int. J. Multiph. Flow. 119 (2019) 180–193. ment, Comput. Fluids. 187 (2019) 83–97.
31
[72] X. Chen, V. Yang, Thickness-based adaptive mesh refinement methods for [97] M. Rezavand, C. Zhang, X. Hu, A weakly compressible SPH method for vio-
multi-phase flow simulations with thin regions, J. Comput. Phys. 269 (2014) lent multi-phase flows with high density ratio, J. Comput. Phys. 402 (2020)
22–39. 109092.
[73] K. Schmidmayer, F. Petitpas, E. Daniel, Adaptive Mesh Refinement algorithm [98] G. Pahar, A. Dhar, Mixed miscible-immiscible fluid flow modelling with in-
based on dual trees for cells and faces for multiphase compressible flows, J. compressible SPH framework, Eng. Anal. Bound. Elem. 73 (2016) 50–60.
Comput. Phys. 388 (2019) 252–278. [99] Y.L. Li, L. Wan, Y.H. Wang, C.P. Ma, L. Ren, Numerical investigation of interface
[74] H. Li, W.T. Leung, M.F. Wheeler, Sequential local mesh refinement solver with capturing method by the Rayleigh-Taylor instability, dambreak and solitary
separate temporal and spatial adaptivity for non-linear two-phase flow prob- wave problems, Ocean Eng 194 (2019) 106583.
lems, J. Comput. Phys. 403 (2020) 109074, doi:10.1016/j.jcp.2019.109074. [100] X. He, S. Chen, R. Zhang, A lattice Boltzmann scheme for incompressible mul-
[75] D.J.E. Harvie, M.R. Davidson, M. Rudman, An analysis of parasitic current tiphase flow and its application in simulation of Rayleigh–Taylor instability, J.
generation in volume of fluid simulations, Appl. Math. Model. 30 (2006) Comput. Phys. 152 (1999) 642–663.
1056–1066. [101] Y. Zhao, H.H. Tan, B. Zhang, A high-resolution characteristics-based implicit
[76] E.G. Puckett, A.S. Almgren, J.B. Bell, D.L. Marcus, W.J. Rider, A high-order pro- dual time-stepping VOF method for free surface flow simulation on unstruc-
jection method for tracking fluid interfaces in variable density incompressible tured grids, J. Comput. Phys. 183 (2002) 233–273.
flows, J. Comput. Phys. 130 (1997) 269–282. [102] F. Denner, Wall collision of deformable bubbles in the creeping flow regime,
[77] O. Ubbink, R.I. Issa, A method for capturing sharp fluid interfaces on arbitrary Eur. J. Mech. 70 (2018) 36–45.
meshes, J. Comput. Phys. 153 (1999) 26–50. [103] M.M. Larimi, A. Ramiar, Two-dimensional bubble rising through quiescent
[78] G.Y. Soh, G.H. Yeoh, V. Timchenko, An algorithm to calculate interfacial area and non-quiescent fluid: Influence on heat transfer and flow behavior, Int.
for multiphase mass transfer through the volume-of-fluid method, Int. J. Heat J. Therm. Sci. 131 (2018) 58–71.
Mass Transf. 100 (2016) 573–581. [104] S.T. Zalesak, Fully multidimensional flux-corrected transport algorithms for
[79] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling sur- fluids, J. Comput. Phys. 31 (1979) 335–362.
face tension, J. Comput. Phys. 100 (1992) 335–354. [105] M. Owkes, O. Desjardins, A computational framework for conservative,
[80] A. Bejan, Convection heat transfer, John wiley & sons, 2013. three-dimensional, unsplit, geometric transport with application to the vol-
[81] D. Das, L. Lukose, T. Basak, Role of multiple discrete heaters to minimize en- ume-of-fluid (VOF) method, J. Comput. Phys. 270 (2014) 587–612.
tropy generation during natural convection in fluid filled square and triangu- [106] M.-J. Li, Interaction between free surface flow and moving bodies with a dy-
lar enclosures, Int. J. Heat Mass Transf. 127 (2018) 1290–1312. namic mesh and interface geometric reconstruction approach, Comput. Math.
[82] G.G. Ilis, M. Mobedi, B. Sunden, Effect of aspect ratio on entropy generation with Appl. 81 (2021) 649–663.
in a rectangular cavity with differentially heated vertical walls, Int. Commun. [107] D. Zhang, C. Jiang, D. Liang, L. Cheng, A review on TVD schemes and a
Heat Mass Transf. 35 (2008) 696–703. refined flux-limiter for steady-state calculations, J. Comput. Phys. 302 (2015)
[83] H.K. Versteeg, W. Malalasekera, An introduction to computational fluid dy- 114–154.
namics: the finite volume method, Pearson Education, 2007. [108] N. Ouertatani, N. Ben Cheikh, B. Ben Beya, T. Lili, Numerical simulation of
[84] J.E. Pilliod Jr, E.G. Puckett, Second-order accurate volume-of-fluid algorithms two-dimensional Rayleigh–Bénard convection in an enclosure, Comptes Ren-
for tracking material interfaces, J. Comput. Phys. 199 (2004) 465–502. dus Mécanique 336 (2008) 464–470.
[85] P.K. Sweby, High resolution schemes using flux limiters for hyperbolic con- [109] M.K. Das, K.S.K. Reddy, Conjugate natural convection heat transfer in an in-
servation laws, SIAM J. Numer. Anal. 21 (1984) 995–1011. clined square cavity containing a conducting block, Int. J. Heat Mass Transf.
[86] M. Zijlema, On the construction of a third-order accurate monotone convec- 49 (2006) 4987–50 0 0.
tion scheme with application to turbulent flows in general domains, Int. J. [110] J. Zhao, W. Cai, Y. Jiang, Study on frequency patterns of 2D square
Numer. Methods Fluids. 22 (1996) 619–641. Rayleigh–Bénard convection filled with air, Eur. J. Mech. 74 (2019) 280–290.
[87] W. Gao, H. Li, Y. Liu, Y.-J. Jian, An oscillation-free high order TVD/CBC-based [111] K.T. Yang, Transitions and bifurcations in laminar buoyant flows in confined
upwind scheme for convection discretization, Numer. Algorithms. 59 (2012) enclosures, 1988.
29–50. [112] D. Mishra, K. Muralidhar, P. Munshi, Experimental study of Rayleigh–Benard
[88] N.P. Waterson, H. Deconinck, Design principles for bounded higher-order con- convection at intermediate Rayleigh numbers using interferometric tomogra-
vection schemes–a unified approach, J. Comput. Phys. 224 (2007) 182–207. phy, Fluid Dyn. Res. 25 (1999) 231.
[89] B. Van Leer, Towards the ultimate conservative difference scheme. II. Mono- [113] R. Anandalakshmi, T. Basak, Natural convection in rhombic enclosures with
tonicity and conservation combined in a second-order scheme, J. Comput. isothermally heated side or bottom wall: entropy generation analysis, Eur. J.
Phys. 14 (1974) 361–370. Mech. 54 (2015) 27–44.
[90] M.S. Darwish, F.H. Moukalled, Normalized variable and space formulation [114] S.G. Schladow, Oscillatory motion in a side-heated cavity, J. Fluid Mech. 213
methodology for high-resolution schemes, Numer. Heat Transf. 26 (1994) (1990) 589–610.
79–96. [115] S.C. Saha, Y. Gu, Natural convection in a triangular enclosure heated from be-
[91] S. Bidadi, S.L. Rani, Quantification of numerical diffusivity due to TVD low and non-uniformly cooled from top, Int. J. Heat Mass Transf. 80 (2015)
schemes in the advection equation, J. Comput. Phys. 261 (2014) 65–82. 529–538.
[92] D.L. Youngs, Time-dependent multi-material flow with large fluid distortion, [116] Z. Alloui, P. Vasseur, M. Reggio, Natural convection of nanofluids in a shallow
Numer. Methods Fluid Dyn. (1982). cavity heated from below, Int. J. Therm. Sci. 50 (2011) 385–393.
[93] R. Scardovelli, S. Zaleski, Interface reconstruction with least-square fit and [117] R.I. Issa, Solution of the implicitly discretised fluid flow equations by opera-
split Eulerian–Lagrangian advection, Int. J. Numer. Methods Fluids. 41 (2003) tor-splitting, J. Comput. Phys. 62 (1986) 40–65.
251–274. [118] A.J. Chorin, Numerical solution of the Navier-Stokes equations, Math. Comput.
[94] X. Xu, X.-L. Deng, An improved weakly compressible SPH method for simulat- 22 (1968) 745–762.
ing free surface flows of viscous and viscoelastic fluids, Comput. Phys. Com- [119] S.-W. Kim, T.J. Benson, Comparison of the SMAC, PISO and iterative time-ad-
mun. 201 (2016) 43–62. vancing schemes for unsteady flows, Comput. Fluids. 21 (1992) 435–454.
[95] X. Zheng, Q.W. Ma, W.Y. Duan, Incompressible SPH method based on Rankine [120] Z. Nasri, A.H. Laatar, J. Balti, Natural convection enhancement in an
source solution for violent water wave simulation, J. Comput. Phys. 276 (2014) asymmetrically heated channel-chimney system, Int. J. Therm. Sci. 90 (2015)
291–314. 122–134.
[96] A. Colagrossi, M. Landrini, Numerical simulation of interfacial flows by [121] C. Varsakelis, M.V Papalexandris, A numerical method for two-phase flows of
smoothed particle hydrodynamics, J. Comput. Phys. 191 (2003) 448–475. dense granular mixtures, J. Comput. Phys. 257 (2014) 737–756.
32

International Journal of Heat and Mass Transfer: Faroogh Garoosi, Tew-Fik Mahdi

Uploaded by

Copyright:

Available Formats

International Journal of Heat and Mass Transfer: Faroogh Garoosi, Tew-Fik Mahdi

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

International Journal of Heat and Mass Transfer: Faroogh Garoosi, Tew-Fik Mahdi

Uploaded by

Copyright:

Available Formats

International Journal of Heat and Mass Transfer 172 (2021) 121163

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer

Presenting a novel higher-order bounded convection scheme for

1. Introduction cavities. In a pioneering work, Amine et al. [6] experimentally and

sure with differentially heated horizontal walls where temperature

with  n standing for the interface normal vector. The additional

0, u = v = 0) under the zero-pressure ﬁled ( p = 0). It is worth ⎧ 1

nx >0 ny < 0 |mi | > 1 A1

Present model/Eulerian 5.72 × 10−3

instants is portrayed in Fig. 12. It can be observed that contrary Table 3

case 7, to induce the buoyancy force and Rayleigh-Bénard Convec-

Number of grids (Case 7)

Ra 39×39 59×59 79×79 99×99 119×119 139×139 159×159

104 2.1675 2.1832 2.1966 2.2051 2.2134 2.2168 2.2177

Ra 39×39 59×59 79×79 99×99 119×119 139×139 159×159

104 2.2505 2.2601 2.2692 2.2784 2.2825 2.2846 2.2853

Ra Nu |Umax | |ψmax | Be Stot Stot,max ST ST,max SF SF,max

Ra Nu |Umax | |ψmax | Be Stot Stot,max ST ST,max SF SF,max

5. Conclusion This work is ﬁnancially supported by the Merit scholarship pro-

You might also like

with n standing for the interface normal vector. The additional