12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
I UNSTABLE NATURAL CONVECTION IN

A VERTICAL CHANNEL WITH HOT-COLD

WALL CONFIGURATION
ChungGang Li
Department of Computational Science
Kobe University
1-1 Rokkodai, Nada-ku, Kobe 657-8501, Japan
cgli@aquamarine.kobe-u.ac.jp

Makoto Tsubokura
Department of Computational Science
Kobe University
1-1 Rokkodai, Nada-ku, Kobe 657-8501, Japantsubo@tiger.kobe-u.ac.jp

ABSTRACTThe unstable phenomena induced by the natural convection
in an open-ended vertical channel with a hot-cold wall
configuration are investigated. The compressible solver
combining absorbing and non-reﬂecting boundary conditions
for extremely low Mach numbers is applied to eliminate theproblem of requiring a priori knowledge of the ﬂow rate. It isfound that the existing empirical experience for the naturalconvection in a vertical channel cannot be applied to thecurrent configuration. In order to better reflect the underlying
physics, the turbulent behavior near the wall is investigated andthe result shows that Grossmann-Lohse (GL) theory, originallyproposed of horizontal (Rayleigh-Benard) convection, can beapplied to describe the current flow behavior. From thecorrelation between the Rayleigh number and Nusselt number,it can be estimated that this kind of unstable phenomenon canenhance the heat transfer by around 10%.

applied to investigate the flow and temperature fields to better
understand the physics. Moreover, to further understand themechanism of unstable phenomena, the effect of the Rayleighnumber on the heat transfer is also investigated. From the result,it can be known that, even at a low Rayleigh number, once theinteraction between hot and cold walls is strong enough,instead of the laminar flow field, the similar Rayleigh-Benard
convection is shown and this kind of convection can enhance
the heat transfer.

Physical model and Governing Euqation
The open-ended finite length channel with extra regions isshown in Fig. 1. The streamwise, vertical and spanwise
directions are denoted x1, x2 and x3, and their corresponding
velocities are u1, u2 and u3, respectively. Gravity is in thenegative x1 direction. For simplicity, only the x1 x2 plane isshown here as the boundary condition in the spanwise direction,
x3, is periodic. The channel has width l2 and the origin is at
point c. The whole computational domain can be separated tothree parts, the real physical domain abdc, the lower artificial
buffer zone a0b0ba, and the upper artificial buffer zone cdd0c0.
The related parameters used in this study are listed in Table 1.
The real physical domain abdc consists of the heat wall
ac with an isothermal condition Th and the cold wall bd with anisothermal condition Tc. The governing equations here are theoriginal Navier–Stokes equations,
U F1 F2 F3



S.
(1)
t x1 x2 x3

Introduction
The natural convection in a vertical channel with hot-coldwalls is very common to see in our daily life such as a wallmounted-panel heater. Differently from one of the most wellknown cases, the natural convection in an enclosure, how towell conditionally assign the boundary condition at inlet andoutlet is a challenging topic. Especially, due to the different
direction of the buoyancy force in a hot-cold wall configuration,
the opposite velocities located on the same boundary for both
upper and lower apertures in the channel increases thedifficulty of performing the simulation. Besides the issue of theboundary condition, to the authors’ best knowledge, the current
configuration has not been numerically studied before, so how
to appropriately analyze the results and describe the flow
characteristics using the known flow patterns to understand thenew physical insights should be also taken into consideration.
Therefore, this study aims to investigate the unstablephenomena induced by the natural convection in a verticalchannel with asymmetric hot-cold walls at moderate Rayleighnumbers from 5.2×105 to 7.6×106. To the authors’ best
knowledge, the physical model proposed here has not beenstudied yet due to the opposite directions of the velocity at thesame boundary condition. To overcome this problem, the nonreflecting boundary condition with the absorbing boundary
condition developed in [1] is applied to the current physicalmodel. Besides, the unstable phenomena near the wall generatesmall flow structures, which is different from past cognitions
about the flow pattern at the same Rayleigh number. To study
the scale of the fluid property, following Ng et al. [2], the GL
theory, originally proposed for the horizontal convection, isi.e. compressible DNS. The quantities included in U and Fiare  ui

  




  uiu1  P i1   Ai1  u1  ,U    u 2  , Fi    uiu2  P i 2   Ai 2  uiu3  P i 3   Ai 3  u3  e 
(
e

P
)
u

A
u

k

T
/

x




i
ij j
i

0





(



)
g
0


,
S 
0


0


 (    ) gu 
0
1


(2)

where Aij  u j / xi  u j / xi  2 / 3(iu ) ij and the ideal gas
equation is written as p   RT1

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022[1], for the artificial buffer zones a0b0ba, 1 and 1 areexpressed as1  Vout  c  [( x1  ac ) / aa0 ]3 ,(7)

The viscosity and thermal conductivity of the fluid are based
upon Sutherland’s law:
T 3 T0  110,
T0 T  110

 (T )  0 ( ) 2(3) 1   out  c  [( x1  ac ) / aa0 ]3 , (T ) R
k (T ) ,(  1) PrWhere 0  1.1842kg / m3

 0  1.85  105 N  s / m 2 ,,g  9.8m / s 2

T0  298.0592 K,(8)

where Vout and  out are 1.25 and c is calculated from Eq. (9).

(4)

c  0.5  u 2 (  1)2  4c 2 ,,(9)where   100(U 0 / c ) and c is the speed of sound. Similarly,2  1.4 ,for the artificial buffer zones cdd0c0, 1 and 1 are expressedR  287 J / kg / K and Pr  0.72 .

as1  Vin  c  ( x1 / cc0 )3 ,(10)

 1   out  c  ( x1 / cc0 ) ,

(11)

3

however,  1 is set as 0.035 to not to excessively damp out thefluctuation when the flow enters the real physical domain.
In addition to the absorbing boundary condition, toprevent thermal diffusion from the edge of the arstificial bufferzone back into the domain, thereby polluting the results, themodified local one-dimensional inviscid (LODI) relations isapplied at the edge of the artificial buffer zone as the nonreﬂecting boundary condition. For more detail information
about the absorbing and non-reﬂecting boundary conditions
mentioned above, Li et al. [1] can be referred.

Numerical Method
The fluid speed induced by the natural convection is nearly the
incompressible flow so the original compressible solver should
be modified for this study. The original Roe scheme with apreconditioning method and dual time stepping, which has
been proved to be a powerful tool for the transient naturalconvection simulation in [3,4,5] is applied to resolve Eq. (1).
Hence, the new governing equations employed to thecurvilinear coordinates ( , ,  ) to better resolve the near-wallphysics isU p U F1 F2 F3





S ,
(9)

t
  where  is the preconditioning matrix derived by Weiss andSmith [6], U p is the primitive form of [ P , u1 , u2 , u3 , T ] / J in

Figure 1. Physical Model

Owing to the buoyancy force generated by both the heatwall ac and cold wall bd that drives the fluid flow, it can be
expected that the opposite direction of the velocity will besimultaneously existed on the ab and cd. This means that theedge a0b0 and c0do should be as far as possible from the domain
of interest, abcd, to mimic the infinite environment. Therefore,
the artificial buffer zones, a0b0ba and cdd0c0, are added to thereal physical domain to reduce the computational requirementsof mimicking the infinite environment.
U F1  F2 F3  S ,(5)tx1 x1 x2 x3
where    target 
  




  u1  (  u1 ) target   u1   1   u 2  , and    1   u2  (  u2 ) target  ,
(6)




u

(
u
)


 3
3 target   u3  e   e  (  e)target which J is the Jacobian transformation matrix,  and t are
the artificial and physical times, respectively, and U is theconservative form of (  ,  u1 ,  u2 ,  u3 ,  e) / J . For S , in thereal physical domain, abdc, S  S / J ; in the artificial bufferzone, a0b0ba and cdd0c0, S  ( / x1   ) / J  S / J .
The discretized form of Eq. (9) is
U  U pk 3U k 1  4U n  U n 1


2t

1
1

( F1k 11  F1k 11 ) ( F k 1  F2 k 1 1 ) ,i , j,ki, j  ,k 2i , j  12 , k
 i  2 , j , k

2
2
k 1
p

where t arg et is the target value at infinity. So, target  1.1842kg / m3 , ui , target  0m / s and ptarget  101300Pa .



For the artificial buffer zones, the absorbing boundary
condition is applied. Comparing with Eq. (1), the governing
equations in artificial buffer zones have two additional terms;
the artificial convection term,  / x1 , to drive the convection(10)

1( F3k 1 1  F3k 1 1 )  S
i, j,k  i , j , k  2
2

The terms U p /  and U / t are approximated by afirst-order forward difference and a second-order backward
difference, respectively, and the terms F1 /  , F2 /  , andvelocity to the desired velocity at a0b0 and c0do, and thedamping term,  , to eliminate disturbances. Following Li et al.

F3 / 2

are approximated by a central difference. The

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

superscripts k and n indicate the iteration numbers of artificial
time and the proceeding step of real time, respectively. Whenthe artificial time U p /  converges to  ( 102 ) , theFor the boundary condition, all walls are assumed to be
isothermal and no-slip conditions. For a0a, b0b, cc0 and dd0 in
the lateral direction, Neumann condition is adopted. In the x3
direction, the periodic condition is also implemented to reducethe computational effort.

magnitude of the ( k  1)th iteration of the artificial time term isapproximately equivalent to the magnitude of the ( n  1)thtime step of the real time, reducing Eq. (10) to the originalNavier–Stokes equation including the transient term.
Finally, Eq. (10) can be rearranged asI3[  1M


2t
,
(11)
 1 ( Apk   B pk    C pk )]U p   1R k

Result and Discussion
The natural convection in a vertical channel as shown in
Fig. 1 is investigated and the Rayleigh number based on thewidth of the channel is 5.2×105. Two different cases are
conducted. One is that the wall temperatures in the left andright are Tc  T0  20 K and Th  T0  40 K , respectively(CaseI). The other is that the wall temperatures in the left andright are Tc  T0 and Th  T0  60 K , respectively (CaseII). The
detail computational parameters are listed in Table. 1. Please
note here. The temperature differences between two walls forthese two cases are the same so that the Rayleigh numbers arealso identical.

where M  U p /  ,   is the central-difference operator,
Ap  F1k / U p

is

R  S  (3U  4U  U
k

k

n

the
n 1

flux

Jacobianand) / 2t  (  F   F2   F3 )
k
1

k

k

.

Eq. (11) is solved using the Lower-Upper Symmetric-GaussSeidel (LUSGS) implicit method.
In the computation of Rk on the right-hand side of Eq. (11),the terms of Fi in Eq. (2) based on Cartesian coordinates canbe divided into two parts. One is the inviscid term,  ui  uiu1  P i1 Finviscid    uiu2  P i 2  ,  uiu3  P i3  (  e  P)u i

Table. 1 Computational parameters(12)

the other is the viscous term,




0


  Ai1Fviscous     Ai 2
(13)
,
A

 i3


T 
  Aiju j  

xi 

The Roe scheme is employed in the discretization of theFinviscid term and, when the preconditioning matrix is added, isexpressed as
1
1
F( FR  FL )  1 Ap U P ,(14)1 
inviscid ,i 222
For the reconstruction of FR and FL, the fifth-order monotonic
upstream-centred scheme (MUSCL) without a limiter function
to prevent turbulent fluctuations from attenuating, is adopted
and expressed as
U L ,i 1/ 2  1 / 60 ,
(15)
(2U i  2  13U i 1  47U i  27U i 1  3U i  2 )( 3U i  2  27U i 1  47U i  13U i 1  2U i  2 ),2.0  103 s

WallTemperature

CaseI: Tc  T0  20 K ; Th  T0  40K

T  Th  Tc

60 K

H

0.045 m

L

0.4 m

Ra

5.2×105

U T

CaseII: Tc  T0 ; Th  T0  60K

g  TH  0.3

Domain size

0.65  0.045  0.02 m3

x1  x2  x3

1  0.1  2 mm 3

N1  N 2  N 3

400×200×100 (Physical domain)

Fig 2. shows the contour of the instantaneous temperatureand the instantaneous velocity for both CaseI and CaseII. For
the CaseII, except the region near the outlet, the flow field andthe temperature field are almost stationary. On the other hand,
CaseI shows unstable phenomena. Generally speaking, at thismoderate Rayleigh, the natural convection in a vertical channel
should be stationary such as shown in CaseII. Therefore, it canbe known that the existing empirical experience for the verticalchannel flow can’t be applied to Case I.

U R ,i 1/ 2  1 / 60 

t

(16)

Aside from the inviscid term, the derivative terms in Aij arecomputed by a second-order central difference. For theartificial convection term,  / x1 , in Eq. (13), the first-orderbackward difference is used, which is expressed as i  i 1
U  U i 1
 1 i
,
x1
x1
x1

(17)

3

Figure 2. The instantaneous temperature and velocity
contour12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

to the mean temperature,  T  k (d T / d x2 ) and due to the

To investigate the unstable phenomena in CaseI, thestatistical data of the natural convection in a vertical channel
but with periodic conditions in both spanwise and streamwise
directions at the similar Ra number, 5.4×105, in [2] using DNS
is compared in Fig. 3. Please note here. The discrepancy
especially near the center of the channel height, x2 / H  0.5 ,is caused by the different conditions of the temperature on thewall. In [2], the condition of TH  T0  1 / 2  T andtemperature fluctuations,  T   k (T  / x j ) 2 are shown. Basedon these results, inside the boundary layer, the physicalproperties of the mean dominate the flow field. On the other
hand, outside the boundary layer, the physical properties of thefluctuation dominate the flow field. The flow field shows bothlaminar and turbulent behaviors, which is consistent with the
GL theory, therefore, can be explained as follows: because
away from the wall, the turbulence due to the naturalconvection is not strong enough to generate the boundary layer,
the kinetic dissipation is eventually decomposed into thelaminar behavior inside the boundary layer and turbulent
behavior outside the boundary layer.

Tc  T0  1 / 2  T is assigned. On the other hand, in CaseI, thecondition of TH  T0  2 / 3  T and Tc  T0  1 / 3  T is
assigned. The asymmetric distribution indicates that although
the results are qualitatively good agreement with the DNS
results, the current boundary condition can count the effectfrom the ambient condition that the periodic condition adopted
in [2] couldn’t. Besides, it can be also confirmed that theunstable phenomena here can be maintained to generate theturbulent structures near the wall.

Figure 4. The distribution of the kinematic dissipation and thethermal dissipation
Figure 3. Statistical data

The correlation between the Rayleigh number and Nusselt
number at the half-length of the channel, x1=0.2m, is shown in
Fig. 5. The distribution shows a discontinuous trend between
Ra=4.0×106 and Ra=6.0×106. Different from the rest of theother distribution, with increasing the Rayleigh number, theincrement of the Nusselt number can’t be observed during thisrange. Based on the results of Ra=4.19×106 and Ra=5.74×106
shown in Fig. 5, the reason that the Nu is not increased can beblamed on the lack of these unstable phenomena. At
Ra=4.19×106, the interaction of the heat transfer between twowalls is disappeared so the effect of enhancing the heat transferfrom this instability does not exist anymore. Hence, the trend of
the increment of the Nusselt numbers can’t be maintained when
Ra is larger than 4.19×106. Moreover, by extending the
distribution of Ra<4.19×106, the differences of the distribution
between Ra<4.19×106 and Ra>5.74×106 can be seen as theeffect of the unstable phenomena, which can be estimatedaround 10%. In other words, the unstable phenomena will
increase the performance of the heat transfer by around 10%.

To study the flow and thermal fields near the wall, theboundary-layer behavior is investigated. Following [2], for the
hot wall, the kinetic boundary-layer thickness  u is calculated
by intercepting u  x2  du / dx2 and u  umax . Similarly, thethermal boundary-layer thickness T

is calculated byintercepting T  Th  x2  dT / dx2 and T  Th  2 / 3  T .
Through these calculations, the effects on the flow can beseparated. One is from the boundary layer due to the velocity
gradient, du / dx2 (Inside the BL) and the other is from thebulk region caused by the bulk velocity, u  umax (Outside theBL). These distributions are compared with the kinetic
dissipation and the thermal dissipation in Fig. 4. to better
understand the physical insight. The kinetic dissipation due tothe mean velocity,  u   (d u 2 / d x2 ) and due to the velocity
fluctuations  u    (ui / x j ) 2 , and the thermal dissipation due

412th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 5. The correlation between Ra and Nu

Conclusion
The natural convection in an open-ended vertical channel
with a hot-cold wall configuration has been studied for the firsttime. The compressible solver with the hybrid boundary
condition is utilized to address the problem of the unknown
flow rate at the inlet and outlet. Based on the statistical data, it
can be known that the existing empirical experience for thevertical channel flow cannot be applied to the present study.
Instead, it is found that GL theory, originally proposed forhorizontal convection, can well describe the unstablephenomena here. From the result of the Rayleigh-Nusselt
number correlation, it can be known that the unstablephenomena due to the interaction between two walls canenhance the heat transfer by around 10%.

REFERENCES

[1]

[2]
[3]

[4]

[5]

[6]

C.-G. Li, M. Tsubokura, W.-S. Fu, N. Jansson, W.-H.
Wang, Compressible direct numerical simulation with ahybrid boundary condition of transitional phenomena innatural convection, Int. J. Heat Mass Transf. 90 (2015).
C.S. Ng, A. Ooi, D. Lohse, D. Chung, Vertical naturalconvection: Application of the unifying theory of thermal
convection, J. Fluid Mech. 764 (2015) 349–361.
W.-S. Fu, C.-G. Li, C.-C. Tseng, An investigation of adual-reflection phenomenon of a natural convection in athree dimensional horizontal channel without Boussinesq
assumption, Int. J. Heat Mass Transf. 53 (2010) 1575–
1585.
C.G. Li, M. Tsubokura, W.S. Fu, N. Jansson, W.H. Wang,
Compressible direct numerical simulation with a hybrid
boundary condition of transitional phenomena in naturalconvection, Int. J. Heat Mass Transf. 90 (2015) 654–664.
C.G. Li, A compressible solver for the laminar–turbulent
transition in natural convection with high temperature
differences using implicit large eddy simulation, Int.
Commun. Heat Mass Transf. 117 (2020).
J.M. Weiss, W.A. Smith, Preconditioning applied tovariable and constant density flows, AIAA J. 33

5