Numerical investigation of laminar compressible natural convection flow in asymmetrically and isothermally heated open-ended inclined channel

International Journal of Heat and Mass Transfer 130 (2019) 83–97

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Numerical investigation of laminar compressible natural convection flow in asymmetrically and isothermally heated open-ended inclined channel Deboprasad Talukdar a,⇑, Chung-Gang Li b, Makoto Tsubokura b a b

Graduate School of System Informatics, Kobe University, Kobe 657 8501, Japan Faculty of Graduate School of System Informatics, Kobe University, Kobe 657 8501, Japan

a r t i c l e

i n f o

Article history: Received 8 June 2018 Received in revised form 11 September 2018 Accepted 18 October 2018

Keywords: Numerical simulation Laminar natural convection Compressible flow Inclined channel Asymmetrical and isothermal heating Composite correlation

a b s t r a c t This study presents a numerical investigation of laminar compressible natural convection flow induced under high-temperature difference in an open-ended inclined parallel-walled channel heated isothermally and asymmetrically at the top wall (hot plate facing downwards) while the lower plate is considered adiabatic. The investigation is carried out with air ðPr ¼ 0:72Þ as the working fluid for a temperature difference of 110 °C from the ambient, over of range of modified Rayleigh number 6:2  101 to 1:6  104 and inclination of channel varying between 10 and 90° with respect to horizontal position. For the simulation, compressible governing equations without Boussinesq approximation is solved by employing new modified all-speed Roe scheme with matching preconditioning method, dual time-stepping technique and utilizing the modified Local One-dimensional Inviscid (LODI) relations suitable for compressible natural convection flow as non-reflecting boundary conditions at the inlet and outlet of the channel. Average Nusselt number based on inter-plate spacing increases with the increase of modified Rayleigh number and also with the increase of inclination from the horizontal position. This increase in average Nusselt number can be described by the variation of the level of thermal saturation inside the channel. Significant reduction in average Nusselt number is observed in the combined zone of lower angle of inclination and lower modified Rayleigh number. Visualization of fluid flow (temperature contour and velocity magnitude contour) indicates the existence of both fully developed flow and developing flow regime within the span of investigated modified Rayleigh number. Based on the investigation, a composite correlation between average Nusselt number and product of modified Rayleigh number and the sine of the angle of inclination considering as a single parameter is presented which can be conveniently used for designing of engineering application dealing with natural convection heat transfer under hightemperature difference. Ó 2018 Elsevier Ltd. All rights reserved.

1. Introduction Passive systems based natural convection heat transfer in an open-ended inclined channel are of considerable interest because of its application in cooling of electronic cabinets, solar energy collection system and building energy management through natural ventilation. Many past experimental and numerical investigations have been undertaken to understand the impact of channel aspect ratio and angle of inclination in determining the overall thermal performance in the channel.

⇑ Corresponding author. E-mail addresses: [email protected] (D. Talukdar), cgli@aquamarine. kobe-u.ac.jp (C.-G. Li), [email protected] (M. Tsubokura). https://doi.org/10.1016/j.ijheatmasstransfer.2018.10.076 0017-9310/Ó 2018 Elsevier Ltd. All rights reserved.

Azevedo and Sparrow [1] conducted experimental investigations on an inclined isothermal channel at three different inclinations ðh ¼ 0 ; 30 and 45 from vertical position) for three different heating modes – both walls heated, only top wall heated and only bottom wall heated with water as working fluid. All the average Nusselt number results obtained for only top wall heated channel was correlated by using a single parameter in form of ðS=HÞRaS cos h where h is the angle of inclination from the vertical position. Moreover, a global correlation based on the heat transfer results for all investigations including the heating modes, inclinations of the channel, inter-plate spacing and modified Rayleigh number was suggested acceptable within 10 percent deviation. Onur et al. [2] and Onur and Aktas [3] performed an experimental study on asymmetrically and isothermally heated inclined channel (inclined at 0°, 30° and 45° with the vertical) with air. Experiments

84

D. Talukdar et al. / International Journal of Heat and Mass Transfer 130 (2019) 83–97

Nomenclature T0 Tw DT P P0 Pr RaL RaS ðS=LÞRaS NuS mx L0 L1 L S £ Nx Ny x; y u; v R

Ambient Fluid Temperature (K) Hot Wall Temperature (K) Temperature difference (K) pressure ({Pa}) Ambient Fluid Pressure ({Pa}) Prandtl number Rayleigh number based on length of heated surface Rayleigh number based on length of width of channel modified Rayleigh number average Nusselt number mass flow rate (kg/s) length of channel section upstream to-heated surface (m) length of channel section downstream to-heated surface (m) length of heated surface (m) channel inter-plate spacing (m) angle of inclination number of grids in x coordinate number of grids in y coordinate Cartesian coordinates (m) velocities in x and y direction ({m/s}) gas constant (J/kg/K)

were carried out with various temperature differences and interplate spacing. The representation of the results in form of the effectiveness of the opposite wall showed that the heat transfer depends on both inter-plate spacing and channel inclination. The Nusselt number for each particular angle of inclination in [3] is correlated with the modified Rayleigh number in a form similar to that of Azevedo and Sparrow [1]. Baskaya et al. [4] conducted three-dimensional laminar numerical investigations for top wall heated (maintained at temperature T H Þ and bottom wall passively heated (maintained at temperature T C ) channel inclined at 0 ; 30 and 45 with the vertical position and for plate separation ranging from 2 mm to 33 mm. The temperature difference in the simulation was considered close to 50 °C. So, in accordance with the study by Gray and Giorgini [5], compressible governing equations without Boussinesq approximation was solved. The average Nusselt number presented in form of NuS ¼ a½ðS=HÞRaS b had good agreement with the experimental findings of Onur and Aktas [3] and concluded that overall heat transfer depends both on plate separation and channel inclination. Said et al. [6] produced numerical results for incompressible natural convection in a twodimensional symmetrically heated inclined parallel-walled channel. Results were obtained for a range of modified Rayleigh number 0–1000 and range of angle of inclination 0–90° from the vertical position for a single aspect ratio of the channel. Their results indicated a reduction of average Nusselt number with the increase of inclination angle and especially significant reductions were observed at the very high angle of channel inclination. The results for the average Nusselt number was graphically represented versus single parameter ðS=HÞRaS cos h similar to [1] but no correlation has been suggested. Straatman et al. [7] undertook a numerical investigation to understand the effects of channel inclination (0–30° with vertical) on local heat flux distribution and flow behaviour near the inlet of a channel heated isothermally and symmetrically. Manca et al. [8] performed experimental analysis of laminar natural convection in channel inclined in range 60—90 from the vertical position for both symmetric and axisymmetric (the configuration similar to [1]) heating with uniform heat flux and air as working fluid. A glo-

g M ref M local

acceleration due to gravity (m=s2 Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Mach number ( ðgbDTLÞ=ðcRT 0 Þ) local Mach number

Greek symbols b volumetric expansion coefficient ð1=KÞ l dynamic viscosity ðNs=m2 Þ m kinematic viscosity, l=q (m2 =sÞ c specific heat ratio q density ðkg=m3 Þ q0 ambient Fluid density ðkg=m3 Þ a Grid Density Constant Superscripts k artificial time stepping iteration count n physical time stepping iteration count Subscripts i grid index in stream-wise direction j grid index in wall normal direction L; R left and right cell face value

bal correlation for the Nusselt number with that of the angle of inclination and Rayleigh number for all modes of heating and all configurations was presented. Manca and Nardini [9] provided a composite relation between average Nusselt number based on the average wall-temperature, modified Rayleigh number and angle of inclination for symmetric and asymmetric heating. Mittelman et al. [10] reported a composite (fully developed flow and developing flow) relationship between average Nusselt number based on wall temperature, channel modified Rayleigh number and channel inclination angle (between 30 and 90° with respect to horizontal) through analytical approach under uniform heat flux boundary conditions at the wall. It is understood from literature review that composite correlations for average Nusselt number related to modified Rayleigh number have been suggested for natural convection flow in an open-ended inclined channel with uniform heat flux and air as working fluid. Whereas for isothermal boundary condition, correlations for average Nusselt number mostly pertaining to developing flow regime has been suggested. In addition, the effect of inter-plate spacing and modified Rayleigh number on the overall heat transfer has also been discussed in past literatures. However, numerical investigation of natural convection flows induced under very high-temperature difference (>30 °C) where consideration of compressibility of fluid along with governing equation without Boussinesq approximation becomes important has not been addressed much except in [4]. But in [4], at the channel inlet, mass inflow is considered proportional to pressure difference along with stagnation pressure set equal to zero. At the outlet boundary of the channel, static pressure was set equal to the pressure of the ambient. Numerical investigation of laminar incompressible natural convection flow in an asymmetrically heated inclined channel by Brangeon et al. [11] showed that the boundary conditions set in [4] are a kind of global boundary condition. The study in [11] highlights the importance of choice of local conditions for pressure as the inlet-outlet boundary condition in natural convection flow similar to a channel-chimney system. It was concluded in [11] that local pressure boundary condition provides improved results in comparison to the adoption of global conditions at the inlet and

85

D. Talukdar et al. / International Journal of Heat and Mass Transfer 130 (2019) 83–97

outlet of the channel. Fu et al. [12] proposed Modified Local Onedimensional Inviscid (LODI) relations suitable to be adopted as boundary conditions at the inlet and outlet of the open-ended channel for compressible natural convection flow. This boundary conditions have been demonstrated to depend on local values and do not require prior knowledge of the velocity and pressure values at the inlet and outlet. Hence simulation of compressible natural convection flow induced under high-temperature difference of 110 °C in top-wall heated open-ended inclined channel adopting the modified LODI relations differentiates the present investigation from past works of literature. The present study also aims to propose a composite correlation for average Nusselt number suitable for compressible natural convection in asymmetrically and isothermal heated channel (top wall heated). Quere et al. [13] highlighted that natural convection flow induced under high-temperature difference similar to consideration in the current study, becomes compressible and under such condition numerical investigation with compressible flow solver is essential. Numerical investigations of compressible natural convection were conducted by utilizing the Roe scheme with matching preconditioning technique in [12,14,15]. But the accuracy of the preconditioned Roe scheme depends on the accurate prediction of the global cut-off Mach number which is required to treat the pseudo sound speed terms both in numerator and denominator of the dissipation term. Global cut-off Mach number is usually considered either based on the inlet Mach number or maximum Mach number in the flow. In case of the natural convection since no information about the flow is known prior to the simulation, accurate prediction of global cut-off Mach number becomes difficult. Many versions of low-Mach fix for the shock capturing scheme have been discussed in [16] and also in [17,18]. Li and Gu [16] proposed a new modified preconditioned All-speed Roe scheme where the pseudo sound speed terms in the numerator and denominator of dissipation are treated differently based on the local Mach number and global cut-off Mach number respectively. In [16], this scheme was shown to achieve better accuracy in comparison to the preconditioned Roe scheme in investigating flow with Mach number 0.001. In [19], this new modified preconditioned Allspeed Roe scheme was employed to investigate compressible laminar natural convection flow induced under high-temperature difference and was shown to produce accurate and competent results for such applications. So similar to [19], new preconditioned All-speed Roe scheme is employed in the current investigation. The present numerical investigation focuses on the effect of angle of inclination and channel inter-plate spacing on overall heat transfer in a parallel-walled inclined channel. The channel is heated asymmetrically and isothermally at the top wall while the lower wall is maintained adiabatic. The investigation is conducted with 1

4

a wide range of modified Rayleigh number 6:2  10 to 1:6  10 and a wide range of angle of inclination 10—90 with respect to horizontal position. The overall heat transfer is presented in terms of average Nusselt number based on channel inter-plate spacing. Mathematical representation relating the average Nusselt number with modified Rayleigh number and angle of inclination is presented which can be assumed to be valid for the investigated range of modified Rayleigh number and angle of inclination.

Fig. 1. The physical model of open-ended inclined channel.

RaL ¼

source length ðLÞ in the channel is calculated based on RaL 1  107 where

m2

ð1Þ

Pr

where g is acceleration due to gravity, b is volumetric expansion coefficient, m is the kinematic viscosity of the fluid, L is the length of the heat source and Pr is the Prandtl number of fluid. All fluid properties used in Eq. (1) are calculated at ambient temperature ðT 0 Þ. The dimension of the heat source and channel considered herein the investigation is provided in Table 1. L0 , L and L1 are considered as the global constraint in the current investigation. Modified Rayleigh number ½ðS=LÞRaS  is varied by only varying the channel inter-plate spacing ðSÞ. The study is conducted with air ðPr0:72Þ as working fluid and a temperature difference ðDTÞ of 110 °C is imposed of the heated wall ðT w Þ with respect to the ambient temperature ðT 0 Þ. Gravity is considered downward and ambient temperature (T 0 ) and ambient pressure (P 0 ) are 298:0592 K and 101300 Pa respectively. Air is drawn inside the channel at T 0 due to the pressure difference induced by buoyancy force generated. For the current investigation, compressible governing equations without Boussinesq approximation is considered. Such considerations lead to the following compressible governing equations in conservation form for a calorically and thermally perfect gas with constant Prandtl number and constant ratio of specific heat capacities:

@U @F @G þ þ ¼S @t @x @y

ð2Þ

where

2

3

q

6 qu 7 7 6 U¼6 7 4 qv 5

ð3Þ

qE 2 6 6 F¼6 4

6 6 G¼6 6 4

3

qu

7 qu2 þ p  sxx 7 7 5 quv  sxy @T qEu þ pu  k @x  usxx  v sxy

2

2. Mathematical formulation and numerical method The physical model considered herein the investigation is a two-dimensional open-ended parallel-walled inclined channel of finite length as shown in Fig. 1. The angle of inclination ð£Þ for the channel is considered from the horizontal. The isothermal heat

gbDTL3

ð4Þ

3

qv

7 quv  syx 7 7 2 7 qv þ p  syy 5 @T qEv þ pv  k @y  usyx  v syy

ð5Þ

Table 1 Dimension of heated surface and channel. RaL

LðmÞ

L1 ðmÞ

L0 ðmÞ

Nx  Ny

Dxmin  Dymin ðmmÞ

107

0:0978

0:0978

0:0489

500  40

0:489  0:14

86

D. Talukdar et al. / International Journal of Heat and Mass Transfer 130 (2019) 83–97

2 6 6 S¼6 4

0 qgSinð£Þ

qgCosð£Þ

qðugSinð£Þ þ v gCosð£ÞÞ

3

The flux terms consist of inviscid and viscous terms.

7 7 7 5

ð6Þ

p

qðc  1Þ

1 þ ðu2 þ v 2 Þ 2

ð7Þ

  @u @ v @u þ þ 2l @x @y @x

ð8Þ

2 3

  @u @ v @v þ 2l þ @x @y @y

ð9Þ

syy ¼  l sxy ¼ syx

G ¼ Ginv scid þ Gv iscous

ð15BÞ

F iþ1;j ¼ F c þ F d

For the convective part, values of the variables at the cell interface are obtained by interpolating values of the variables initialized at the cell centre using fifth order Upstream-centered Scheme for Conservation Laws (MUSCL) without limiters as proposed by Abalakin et al. [23] in the following way:

ð10Þ

2

U Riþ1 2

The equation of state is considered as follows:

P ¼ qRT

ð11Þ

In this study, temperature dependent dynamic viscosity and thermal conductivity are given by Sutherland’s law as follows:

lðTÞ ¼ l0 kðTÞ ¼

T T0

32

T 0 þ 110 T þ 110

ð12Þ

lðTÞcR ðc  1ÞPr

ð13Þ

Ns 105 m 2

Reference dynamic Viscosity ðl0 Þ ¼ 1:85  Specific heat ratio (cÞ ¼ 1:4, Acceleration due to gravity ðgÞ ¼ 9:81 m=s2 , Gas Constant ðRÞ ¼ 287 J=Kg=K, Reference density ðq0 Þ ¼ 1:1842 kg=m3 at T 0

It is to be noted that the order of interpolation was reduced nearer to the boundary. Dissipation term under Low-Mach number limit is considered in the following manner similar to description in [16],

8 2 2 2 39 0 > Dq 3 q 3 > > > > > > 7 7 7> 6 6 6 > = 6 DðquÞ 7 6 qu 7 6 n x 7> 1< 7 7 6 6 6 7 Fd ¼  jU n j6 7 þ dU 6 7 þ dp6 7 6 Dðqv Þ 7 6 qv 7 6 ny 7> 2> > > > 5 5 4 4 4 5> > > > > : ; DðqEÞ qH 0 U n ¼ nx u þ ny v

at T 0

 dU ¼ c0 

The physical geometry (Fig. 1) considered in the current investigation is similar to study in [12]. Uniform grid distribution is considered in stream-wise direction with smallest grid spacing (Dxmin Þ kept similar to consideration in [12]. For all the investigations, a non-uniform grid is generated in the wall normal direction. The minimum grid spacing in stream-wise and wall-normal direction is indicated in Table 1. To generate non-uniform grid distribution in wall-normal direction, the governing equations are converted into curvilinear coordinates utilizing Jacobian matrices which are constructed using the following algebraic equation as provided in [14]

   S 1 2ðj  0:5Þ 1 Þ tanh að1  2 tanhðaÞ Ny

ð17Þ

ð18Þ

where

where

yj ¼

1 ½3U iþ2 þ 27U iþ1 þ 47U i  13U i1 þ 2U i2  60 1 ½2U iþ3  13U iþ2 þ 47U iþ1 þ 27U i  3U i1  ¼ 60

U Liþ1 ¼

  @u @ v þ ¼l @y @x



ð16Þ

2

2 3

sxx ¼  l

ð15AÞ

The inviscid term can be further divided into convective ðF c Þ and dissipative terms ðF d Þ,

where

E¼

F ¼ F inv scid þ F v iscous

ð14Þ

where S is the channel inter-plate spacing, Ny is the total number of grids in wall-normal direction, a is the grid density constant and index j varies from 1 to Ny . With a higher value of a, more grids are located near to the wall. Roe scheme in form of curvilinear coordinates has been implemented within the numerical code setup following the detailed report provided by Kermani et al. [20,21]. Procedure to implement the modified All-speed Roe scheme with matching preconditioning method of Weiss and Smith [22] has been discussed in [16]. The implementation of the modified preconditioned All-speed Roe scheme in curvilinear coordinates to develop the current numerical setup has been conducted by the following procedures mentioned in [16,20,21].

dp ¼

ð19Þ 

1  h0 U0 Dp U0 þ DU n U n  h0jU n j 2 ~c ~c 2 qhc

U0 1  h0 U0 Dp þ ½c0  jU n j þ U n qDU n ~c ~c 2

1 h0 ¼ min½M2local ; 1; U 0 ¼ ð1 þ h0 ÞU n 2 ffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2 0 0 4c2 h þ ð1  h0 Þ U 2n and c ¼ 2 h ¼ min½maxðM2ref ; M 2local Þ; 1 and ~c ¼ min½maxðM local ; M ref Þ; 1c

ð20Þ

ð21Þ

ð22Þ

ð23Þ

All variables used to calculate the dissipation part is obtained using the Roe averages,

pffiffiffiffiffi pffiffiffiffiffiffi uL qL þ uR qR pffiffiffiffiffi pffiffiffiffiffiffi ; qL þ qR pffiffiffiffiffi pffiffiffiffiffiffi pffiffiffiffiffi pffiffiffiffiffiffi vL q þ vR q HL q þ HR q ¼ pffiffiffiffiffiL pffiffiffiffiffiffi R ; H ¼ pffiffiffiffiffiL pffiffiffiffiffiffi R qL þ qR qL þ qR pffiffiffiffiffipffiffiffiffiffiffi

q ¼ qL qR ; u ¼

v

ð24Þ

All the primitive variables - q; u; v ; T; P are initially considered at the center of a grid cell. So the second order central differencing method as shown in Eq. (25) is used to discretize and solve the viscous terms at the center ðxi ; yj Þ of grid cell,

  ðuiþ1;j  2ui;j þ ui1;j Þ @ @u ¼ ; @x @x ðDxÞ2   @ @u ðuiþ1;jþ1  uiþ1;j1  ui1;jþ1 þ ui1;j1 Þ ¼ @x @y 4Dx Dy

ð25Þ

D. Talukdar et al. / International Journal of Heat and Mass Transfer 130 (2019) 83–97

The preconditioning matrix demonstrated by Weiss and Smith [22] suitable for compressible Low Mach number flow is utilized in solving the governing Eq. (2). Additionally, dual time stepping technique is adopted for speedy convergence of the solution. First-order forward-difference and central-difference are used for discretizing the temporal and spatial terms in Eq. (2) respectively. Second order Runge-Kutta method is implemented for time integration of the solution.

U 1 ¼ U k þ DsLðU k Þ 1 1 1 U kþ1 ¼ U k þ U 1 þ DsLðU 1 Þ 2 2 2

ð26Þ

where

"  1 # 3Ds ; LðU k Þ ¼ Rk M C þ M 2Dt !  " #

1

 3U k  4U n þ U n1 1 k Rk ¼   F iþ1;j  F ki1;j þ Gki;jþ1  Gki:j1 þ Sk 2 2 2 2 2Dt Dx Dy

where k is inner artificial time stepping iteration count, n is outer physical time stepping iteration count and C is preconditioning matrix The convergence criteria for the inner artificial time stepping is controlled by minimum value among the two parameters: the iteration number set equal to 100 and the following equation,

PN

1 jT

Ds

kþ1

PN

 Tkj

k 1T

< 105

ð27Þ

where N is the total number of grid points in the computational domain and k is the artificial time step iteration count. When the term of artificial time stepping converges then the magnitude of the ðk þ 1Þth term is equal to the magnitude of the ðn þ 1Þth term of the physical time stepping. In [22], Weiss and Smith highlighted that Ds shall be calculated by the following equation for maintaining the stability of multistage explicit time integration scheme and Dt can be chosen as per the desired accuracy.

  Dx Dx2 Ds ¼ Min CFL 0 ; VNM 0 u þc m

ð28Þ

where

k is inner artificial time stepping count L1 ¼ u

@T ð1  cÞ @p þ @x qc @x

L2 ¼ u

@v @x 

 @p @u  qðu0  c0  uÞ @x @x



 @p @u  qðu0 þ c0  uÞ @x @x

L3 ¼ ðu0 þ c0 Þ

L4 ¼ ðu0  c0 Þ u0 ¼

1 1 ð1 þ hÞu; c0 ¼ 2 2

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4hc2 þ ð1  hÞ2 u2

and h ¼ min½maxðM2ref ; M 2local Þ; 1 The gradient terms in L are calculated using first-order forwarddifference at the inlet and first-order backward difference at the outlet boundaries. Additionally, similar to description in [12], the magnitude of L4 at the outlet-boundary and magnitude of L3 at the inlet-boundary is assigned a value equal to zero to avoid the reflection of acoustic waves into the computational domain. No-slip boundary condition is considered on all wall surfaces. An adiabatic boundary condition is considered on all wall surfaces except the heated surface segment. The wall boundary conditions are set in the following manner as shown in Eqs. (30) and (31) using one ghost cell grid outside the wall indicated with 0 and first grid from the wall indicated with 1 as shown in Fig. 2 below, Boundary conditions for the adiabatic wall are set using index ði; 0Þ and ði; 1Þ at the center of grid cell as follows,

uði; 0Þ ¼ uði; 1Þ; Pði; 0Þ ¼ Pði; 1Þ;

v ði; 0Þ ¼ v ði; 1Þ; Tði; 0Þ ¼ Tði; 1Þ

ð30Þ

Boundary conditions for the isothermal wall are set using index ði; 0Þ and ði; 1Þ at the center of grid cell as follows,

where CFL is the Courant-Friedrichs-Lewy condition, VNM is Von Neumann number, u0 is the preconditioned velocity, c0 is the preconditioned speed of sound and m is kinematic viscosity In the present investigation, Ds is calculated by Eq. (28). However, the CFL number has to be kept very low for maintaining the stability of the scheme. Dt equal to 0:0005 s is considered in the present investigation. Fu et al. [12] proposed modified Local One-dimensional Inviscid (LODI) relations and demonstrated its suitability as inlet and outlet boundary conditions especially for compressible natural convection flow. The advantage of such boundary conditions is in its implementation without having any prior knowledge of the flow. Similar inlet and outlet boundary conditions are considered for the current investigation and as described below,

Ds ½L3 ðu0 þ c0  uÞ  L4 ðu0  c0  uÞ 2c0 Ds ½L3  L4  ukþ1 ¼ uk  2qc0

pkþ1 ¼ pk 

v kþ1 ¼ v k þ DsL2

  1 ðc  1Þ 1 T kþ1 ¼ T k  Ds L1 þ þ 0 fL3 ðu0 þ c0  uÞ  L4 ðu0  c0  uÞg 2c q c ð29Þ

87

Fig. 2. Numerical grid distribution on the wall.

88

D. Talukdar et al. / International Journal of Heat and Mass Transfer 130 (2019) 83–97

uði; 0Þ ¼ uði; 1Þ; Pði; 0Þ ¼ Pði; 1Þ;

v ði; 0Þ ¼ v ði; 1Þ; Tði; 0Þ ¼ 2T w  Tði; 1Þ

ð31Þ

The investigation is conducted with a wide range of modified Rayleigh number 6:2  101 to 1:6  104 by changing the channel inter-plate spacing and a wide range of angle of inclination 10  90 with respect to horizontal position.

3. Results and discussion 3.1. Grid independence test The total height of the channel, the width of the channel corresponding to the maximum value of modified Rayleigh number 1:6  104 in the current investigation is considered similar to [12]. In [12], uniform grid in stream-wise direction and uniform grid spacing in the wall normal direction was considered for investigating natural convection flow at Rayleigh number (calculated using Eq. (1)) range 104  105 . In the current investigation, similar uniform grid spacing in stream-wise direction is considered. Also, Rayleigh number equal to 107 is used in calculating the length of the heated surface. This results in a thin boundary layer forming close to the wall. So, non-uniform grid distribution is considered in the wall-normal direction in the current investigation by using a value of a ¼ 1:0 in Eq. (14) and hence controlling the grid density near to the wall. Tests are conducted in the vertical channel 𣠼 90 Þ using 500  40 and 250  40 grids orientation in stream-wise and wall-normal direction to assess the effect of grid spacing on the solution. Comparison of the stream-wise velocity,

wall-normal velocity, non-dimensional temperature ððT  T 0Þ =DTÞ measured at the center of the heated surface and local distribution of the Nusselt number along the length of the heated surface is shown in Fig. 3. Fig. 3 show very small difference between the profiles obtained using the two tested grid orientations. Based on the comparison, 500  40 grid orientation is considered in the current study. Additionally, a test conducted with 500  40 and value of a ¼ 1:5 in Eq. (14) resulted in the average Nusselt number being improved by about 1:4%. Also, values of the variables q; u; v ; T; P at the center of the grid located nearest to the wall achieved a value nearer to its respective value assumed at the wall through boundary condition. So, it can be assumed that with a larger number of grids placed nearer to the wall or larger number of grids used in the wall-normal direction, improvement of approximately 1– 1.5% can be achieved from the current presented results. 3.2. Validation of numerical Set-up The validation of the set-up of the numerical scheme in the present investigation is conducted for two cases (i) horizontal channel 𣠼 0 Þ asymmetrically heated only on the lower surface and (ii) asymmetrically heated vertical channel 𣠼 90 Þ. Velocity vector plot, velocity magnitude contour and temperature contour obtained for the horizontal channel is shown in Fig. 4. The velocity vector plot shows that the hot air rises from the bottom surface and exits from both channel openings after flowing parallel to the upper wall. In turn, the cold air enters from both channel openings and flow parallel to the lower walls. This phenomenon creates two circulation zone forming a mirror image by having a demarcation at the middle section of the heated surface

Fig. 3. Comparison of the two test grid - 500  40 and 250  40 for (a) Stream-wise velocity, (b) wall-normal velocity, (c) non-dimensional temperature at the center of the heated surface and (d) local Nusselt number distribution along the heated surface.

D. Talukdar et al. / International Journal of Heat and Mass Transfer 130 (2019) 83–97

89

Fig. 4. (a) Velocity vector profile, (b) velocity magnitude contour and (c) temperature contour for horizontal channel asymmetrically heated on the lower surface.

as observed from the velocity magnitude contour. Temperature contour shows the rise of the thermal plume from the heated surface towards the upper wall and moving towards the exit of the channel. The local Nusselt number variation ðNux ¼ h i L Þ along the heated surface and velocity profile kðTÞ @T k0 ðT w T 0 Þ @y across the channel spacing at location 0:122 m in the axial direction

from inlet is shown in Fig. 5. Maximum velocity is found to occur near to the two wall and minimum velocity (approximately stagnant region) occurring at the channel middle section in the axial direction. The minimum velocity region in the middle section of the channel is also visualized in velocity vector plot shown in Fig. 4. This flow phenomenon captured in the current simulation is found to be similar as described in [6] for the horizontal channel.

Fig. 5. (a) Local Nusselt number variation along the heated lower surface of the channel, (b) velocity profile across the channel spacing at location 0.122 m measured from channel inlet in axial direction.

90

D. Talukdar et al. / International Journal of Heat and Mass Transfer 130 (2019) 83–97

The insensitiveness of the changes in heat transfer in terms of average Nusselt number with increase in inter-plate separation at a given Rayleigh number for the vertical and inclined channel was demonstrated in [3]. For Rayleigh number RaL 3  107 and

Table 2 Comparison of the average Nusselt number obtained in present simulation and correlation in [24] Angle of inclination from horizontal

Average NuL : Present simulation

Average NuL : Correlation in Chen et al. [24]

90°

29.067

28.88

angle of inclination 0° from the vertical (in present investigation £ ¼ 90 Þ, it was shown that opposite wall effectiveness on heat transfer diminishes for aspect ratio near to 0:2. In the present investigation, the average Nusselt number based on the length of the heated surface for RaL 1  107 and aspect ratio equal to 0:2 can be considered approximately equal to heat transfer along a vertical plate. Hence the average Nusselt number obtained is validated against a well-known correlation reported by Chen et al. [24] and shown in Table 2. Since good agreement of the fluid flow dynamics in the horizontal channel and average Nusselt number for the vertical channel has been reached, the numerical setup can be considered adequate for current numerical investigation.

Fig. 6. (A) Velocity vector plot, (B) distribution of pressure distribution ðPaÞ inside the channel and (C) distribution of stream-wise velocity across the non-dimensional width of the channel at section (a) entrance of channel, (b) center of heat surface and (c) exit of channel for channel inclined at 60° and at modified Rayleigh number 1:6  104 .

D. Talukdar et al. / International Journal of Heat and Mass Transfer 130 (2019) 83–97

91

3.3. Inclined channel

all heat transfer in terms of the average Nusselt number can be represented in the form of

Literature surveys indicate that for any angle of inclination, the effect of inter-plate spacing can be represented by modified Rayleigh ½ðS=LÞRaS  which is a combination of channel aspect ratio (Ratio of channel spacing to the length of the heat source) and Rayleigh number based on channel spacing. Thereby, the requirement of separate consideration of inter-plate spacing and Rayleigh number is mitigated. Additionally, the sine of the angle of inclination takes into consideration the effect of the inclination of the channel with respect to gravity. Hence based on the literature survey, over-

NuS ¼ f ððS=LÞRaS ; sin £Þ

ð32Þ

where

RaS ¼ NuS ¼

gbDTS3

m2 hS k

Pr

ð33Þ ð34Þ

Fig. 7. Temperature Contour at angle of inclination (a) 75 , (b) 45 and (c) 15 in sequential manner for modified Rayleigh number (A) 1:52  102 , (B) 1  103 , (C) 5  103 and (D) 1:6  104 .

92

D. Talukdar et al. / International Journal of Heat and Mass Transfer 130 (2019) 83–97

In Eq. (32), ðS=LÞRas and sin £ act as independent variables and in Eq. (34), h is the heat transfer coefficient and k is the thermal conductivity of the fluid. The investigation is performed at 8 different angles of inclination (10 ; 15 ; 20 ; 30 ; 45 ; 60 ; 75 &90 Þ with respect to the horizontal position. For each inclination, dimensionless spacing ðS=LÞ ranging from 0:05 to 0:2 is utilized varying the modified Rayleigh number in the range of 6:2  101  1:6  104 . The flow dynamics inside the inclined channel is discussed through the visualization of the velocity vector profile, pressure distribution profile and velocity profile across the width of the channel at different sections of the channel shown in Fig. 6 for the angle of inclination 60 and at modified Rayleigh number 1:6  104 . It is seen that air is drawn into the channel from the inlet where the velocity profile is almost uniform and accelerated upwards towards the exit of the channel with the maximum being near to the center of the heated surface. In the pressure profile, it is observed that the pressure decreases towards the center of heat surface helping in air being drawn into the channel and then increases further downstream to facilitate the flow of air from inside of channel to outside. Since the channel is heated at the top wall, the flow velocity at the center of the heated surface achieves a greater velocity nearer to the top wall. At the exit section of the channel the flow trends to remain towards the top wall. In Fig. 7, visualization of the temperature contour at a steady state condition for representative four different modified Rayleigh number at the angle of inclination 75 ; 45 and 15 is presented. It is observed that temperature field is fully developed (fluid attains a temperature near to wall temperature very near to the hot surface) for lower modified Rayleigh number and developing for higher modified Rayleigh number for all the angle of inclination. Similar to the description in [25], in Table 3, the average exit temperature expressed in terms of ðT exit av g  T 0 Þ=ðT w  T 0 Þ for varied modified Rayleigh number and varied angle of inclination in investigated. This value near to 1 indicates a high level of thermal saturation inside the channel. It is observed that for the modified Rayleigh

towards the heated surface. It can be assumed that these phenomena influence the reversal of flow from the exit section towards the inside of the channel since starting from this location, the reversed flow is observed to form a V-shaped pattern adjacent to the unheated bottom surface. With decrease in the angle of inclination, the flow reversal phenomenon is observed to get diminished and height and width of V-shaped pattern is observed to get reduced. No flow reversal is observed at the angle of inclination at 15 . At a particular angle of inclination, the flow reversal phenomenon is observed to occur in the higher range of modified Rayleigh number. The V-shaped pattern reduces in height and width with decreasing modified Rayleigh number. It can be assumed that the flow tends to become fully developed with the decrease in modified Rayleigh number which results in the disappearance of the flow reversal phenomenon in the lower range of modified Rayleigh number for all angle of inclinations. Natural convection flows are considered to be laminar at modified Rayleigh number less than 105 [33]. In

number 1:5  102 the average exit temperature is very much near to 1 for the angle of inclination 75 ; 45 and 15 indicating a high level of thermal saturation. The average exit temperature moves away from 1 with the increase of the modified Raleigh number for the angle of inclination 75 ; 45 and 15 : So from this table, it can be considered that within the investigated range of the modified Rayleigh number both thermally developed and developing flow is experienced. The high level of thermal saturation obtained in the channel with very small inter-plate spacing results in lesser heat transfer inside the channel which is reflected in the very low value of average Nusselt number shown later. Velocity magnitude contour shown in Fig. 8 shows the flow variation inside the channel for different combinations of angle of inclination and modified Rayleigh number. At a fixed modified Rayleigh number and at the higher angle of inclination ð75 ; 60 ; 45 Þ of the channel, partial reversal of flow from the exit of the channel is observed. The flow after entering the channel is seen to be sucked Table 3 Effect of average exit temperature verifying the level of thermal saturation inside the channel. Modified Rayleigh number

ðT exit av g  T 0 Þ=ðT w  T 0 Þ Angle of Inclination 15

45

75

0:99

0:87

0:80

1  103

0:79

0:6

0:55

5  103

0:6

0:48

0:45

1:6  104

0:54

0:43

0:4

1:52  10

2

Fig. 8. (A) Velocity Magnitude Contour at RaS ðS=LÞ ¼ 1:6  104 for indicated angle of inclinations of the channel, (B) and (C) Velocity Magnitude Contour for angle of inclination of channel 75 and 45 respectively at RaS ðS=LÞa) 1:6  104 , b) 5  103 , c) 1  103 and d) 1:5  102 .

D. Talukdar et al. / International Journal of Heat and Mass Transfer 130 (2019) 83–97

the current study, the maximum modified Rayleigh number considered is  104 which is less than that of criteria usually considered for laminar flow. Since the flow is well within the laminar regime, it can be assumed that only two-dimensional features shall exist in the flow. The observation from the visualization of velocity magnitude contour in Fig. 8 ensures that the flow hold two-dimensional features even under the presence of strong heat transfer induced

93

under high temperature difference. Additionally, it can be concluded from Fig. 8 that both developing and fully developed flow exist within the investigated range of modified Rayleigh number depending on the angle of inclination. Comparison of velocity profile at the entrance of channel, the center of heat surface and at the exit section of the channel for four representatives modified Rayleigh number at the angle of inclina-

Fig. 9. Velocity profile at (a) entrance of channel, (b) center of heat surface and (c) exit section of channel arranged vertically in serial manner at modified Rayleigh number (A) 1:5  102 , (B) 1  103 , (C) 5:0  103 and (D) 1:6  104 for angle of inclination .

94

D. Talukdar et al. / International Journal of Heat and Mass Transfer 130 (2019) 83–97

tion 75 ; 45 and 15 is presented in Fig. 9. A fully developed velocity profile is observed for lower modified Rayleigh number and developing flow profile for higher modified Rayleigh number for all the angle of inclination. The developed flow regime is seen to spanning across a wider range of modified Rayleigh number for angle of inclination 15 , in comparison to that of an angle of inclination 75 . It is observed that with the change of the inclination of the channel from the vertical position towards the horizontal position, the flow at exit section tends to remain nearer to the top wall.

The observation made in the previous section show that combined developed and developing flow regime is obtained for the range of the modified Rayleigh number investigated. These flow phenomena directly influence the overall heat transfer which is reflected in terms of average Nusselt number. In Fig. 10, the average Nusselt number ðNus Þ is plotted versus the modified Rayleigh number ðRa ¼ ðS=LÞRas Þ for all angle of inclination of channel investigated. It is observed that for a particular angle of inclination, the average Nus increases with the increase of Ra . Idea from

Fig. 10. variation of average Nusselt number and the respective correlation for angle of inclination (a) 90 , (b) 75 , (c) 60 , (d) 45 , (e) 30 and (f) 15 .

D. Talukdar et al. / International Journal of Heat and Mass Transfer 130 (2019) 83–97

literature survey and flow visualization in present investigation suggests that composite correlation (consisting of relation valid for combined developing and developed flow) in the form similar to Sparrow and Azevedo [25] and Bar-Cohen and Roshennow [26] for isothermal channel is suitable to correlate the average Nusselt number (Nus Þ with Ras ðS=LÞ for each angle of inclination. The form of correlation considered is shown in Eq. (35). Hence, the variation of average Nusselt number with modified Rayleigh number and correlation calculated using regression analysis is presented in Fig. 10 for each particular angle of inclination.

" NuS ¼

A ðRas ðS=LÞÞ2

þ

B ðRas ðS=LÞÞ0:5

#0:5 ð35Þ

Additionally, the correlation in the form of NuS ¼ a½RaS ðS=LÞb , reported in experimental investigation by Onnur and Akstas [3] and numerical investigation Baskaya et al. [4] for natural convection in channel heated isothermally at the top wall for angle of inclination 90 ; 60 and 45 has been plotted in the corresponding figures for comparison with past literatures in Fig. 10. It is observed that the average Nusselt number obtained from present investigation increasingly deviates from compared literature with the

95

decrease of the modified Rayleigh number. This is due to the reduction of heat transfer inside the channel with increasing level of thermal saturation. The average Nusselt number obtained for all angle of inclination is presented in one graphical representation versus Ras ðS=LÞ in Fig. 11a. It is observed that with the decrease of the angle of inclination from 90 to 10 the average Nusselt number decreases considerably with the decrease of the modified Rayleigh number, especially beyond 30 . Fig. 11b shows the rate of reduction of the average Nusselt number with the variation of the modified Rayleigh number for particular angle of inclination. The rate of reduction at a particular angle of inclination is calculated with respect to the corresponding higher consecutive Nusselt number value. Major reduction in average Nusselt number is observed in combination region of the angle of inclination less than 30 and lower modified Rayleigh number. Similar to data presented in the literature by Azevedo and Sparrow [1] and SAID et al. [5], the plot of the average Nusselt number (Nus Þ versus ðS=LÞRas Sin£ as a single parameter is shown in Fig. 12. It is noted here that since substantial reduction in value of average Nusselt number is obtained for angle of inclination less than 30 , the data plotting is divided into two sections (i) 90–30° and (ii) less

Fig. 11. (a) Variation of average Nusselt Number with Modified Rayleigh number for each of angle of inclination of channel and (b) percentage of reduction of average Nusselt Number with respect to modified Rayleigh number for each angle of inclination of channel.

Fig. 12. Average Nusselt Number versus the product of modified Rayleigh number and sine of angle of inclination as single parameter for (a) 90–30° and (b) less than 30°.

96

D. Talukdar et al. / International Journal of Heat and Mass Transfer 130 (2019) 83–97

than 30 , This is also done to achieve a better correlation that can be applied to estimate the heat transfer and might offer better design control with minimum deviation from the proposed correlation. In Fig. 12a and b, it is observed that the curves plotted in Fig. 11a at angle of inclination (i) 90–30° and (ii) less than 30 , collapse to a single curve without much observable scattering. The correlation equation obtained through regression analysis is indicated in Fig. 12. Maximum deviation for the average Nusselt number obtained from the present investigation with respect to the predicted average Nusselt number which can be obtained from the correlation is approximately within 8% for angle of inclination between 90° and 30° and approximately 18% for angle of inclination less than 30 . Low Mach number non-Boussinesq approximation was used in [27–31] to investigate the stability phenomena in natural convection and mixed convection flow induced under large temperature difference. Density was considered to be non-linearly and inversely

varying with temperature while the pressure was assumed constant in the calculation of density since the Mach number in the flow is low. Fig. 13a shows the variation of pressure, density and temperature normalized by their respective ambient values at the center of the channel in the axial direction for modified Rayleigh number 1:6  104 and channel inclined at 45 . It is observed that the pressure remains approximately constant along the whole length of the channel while the density is observed to vary in close correlation (inversely proportional) to the temperature distribution. So, the observations in current investigation are found to be similar to the assumption in [27–31], indicating that the variation of density happens primarily due to spatial variation of temperature induced under large temperature difference. Effect of compressibility of flow generated under such a large temperature difference was found to

be small in [27–31] at overall temperature difference e ¼ 12 DT 0T of 0:2 considered in the current investigation. Based on this conclusion,

Fig. 13. (a) Variation of normalized pressure, density and temperature at center of channel inclined at angle of inlcination 45 in the axial direction for modified Rayleigh number 1:6  104 and (b) contour of Mach number variation in the channel at modified Rayleigh number 1:6  104 for angle of inlcination 75–15°.

D. Talukdar et al. / International Journal of Heat and Mass Transfer 130 (2019) 83–97

it is important to understand the effect of compressibility of flow on the results obtained in the current investigation. Mach number contour at modified Rayleigh number 1:6  104 and for the different angle of inclination of the channel is shown in Fig. 13b. For all range of angle of inclinations, the absolute magnitude of the variation of

[3]

[4]

4

Mach number is found to be small, in order of 10 . Since, the compressibility of the flow is proportional to the square of the variation of Mach number of the flow, in the current investigation, the effect of compressibility can be assumed to be small. The flow can be considered to be weakly compressible attributed to the density difference generated due to spatial variation of temperature. Similar to the conclusion drawn by Chenoweth and Paolucci [32], it can be assumed that although the flow is weakly compressible under such high overall temperature difference but the velocity and temperature fields which are observed to be controlling the flow either being developing or developed are still strongly influenced by the large overall temperature difference.

[5] [6]

[7]

[8] [9] [10]

[11]

4. Conclusion The present study reports the numerical solution of compressible natural convection induced under very high-temperature difference in an inclined parallel-walled channel. Simulation is conducted over a range of modified Rayleigh number 6:2  101  1:6  104 and angle of inclination varied from 10° to 90° with respect to the horizontal position. Results show that average Nusselt number increases with the increase of modified Rayleigh number and also increases with the increase of angle of inclination starting from 10 . Significant reduction in average Nusselt number is observed in the combined region of the lower angle of inclination (<30 Þ and lower modified Rayleigh number. This is attributed to the increase of the level of thermal saturation inside the channel. Fully developed flow is obtained at lower modified Rayleigh number and developing flow is observed at higher modified Rayleigh number. The average Nusselt number is correlated with the RaS ðS=LÞSinð£Þ using a composite relation valid for isothermal boundary condition. Based on the outcome of the investigation, correlation is proposed for the angle of inclination (i) 90–30° and (ii) less than 30 . This correlation shall provide a good and reliable estimate of the overall heat transfer inside inclined channel for natural convection flow with a high-temperature difference which can be used in many engineering applications ranging from the cooling of electronic equipment to building energy management.

[12]

[13]

[14]

[15] [16] [17] [18] [19]

[20] [21] [22] [23] [24]

Conflict of interest [25]

We wish to confirm that there are no known conflicts of interest associated with this publication. Acknowledgements The authors wish to acknowledge the financial support received from the Ministry of Education, Culture, Sports Science and Technology (MEXT), Japan in form of scholarship (Scholarship no. 91606000758). A part of this research was supported by MEXT as ‘‘Priority Issue on post-K computer” (Development of Innovative Design and Production Processes that would lead the way for the manufacturing industry in near future. References [1] L.F.A. Azevedo, E.M. Sparrow, Natural convection in open-ended inclined channels, J. Heat Transfer 107 (1985) 893–901. [2] N. Onur, M. Sivrioglu, M.K. Aktas, An experimental study on the natural convection heat transfer between inclined plates (Lower plate isothermally

[26] [27]

[28] [29] [30]

[31]

[32]

[33]

97

heated and upper plate thermally insulated as well as unheated), Heat Mass Transfer 32 (1997) 471–476. M. Onur, M.K. Aktas, An experimental study on the effect of opposing wall on natural convection along an inclined hot plate facing downward, Int. Comm. Heat Mass Transfer 25 (3) (1998) 389–397. S. Baskaya, M.K. Aktas, N. Onur, Numerical simulation of the effects of plate separation and inclination on heat transfer in buoyancy driven open channels, Heat Mass Transfer 35 (1999) 273–280. D.D. Gray, A. Giorgini, The validity of the Boussinesq approximation for liquids and gases, Int. J. Heat Mass Transfer 19 (1976) 545–551. S.A.M. Said, M.A. Habib, H.M. Badr, S. Anwar, Numerical investigation of natural convection inside an inclined parallel-walled channel, Int. J. Numer. Meth. Fluids 49 (2005) 569–582. A.G. Straatman, D. Naylor, J.M. Floryan, J.D. Tarasuk, A study of natural convection between inclined isothermal plates, J. Heat Transfer 116 (1994) 243–245. O. Manca, S. Nardini, V. Naso, 11th Australian Fluid Mechanics Conference, 1992, pp. 131–134. O. Manca, S. Nardini, Composite relations for air natural convection in tilted channels, Heat Transfer Eng. 20 (1999) 64–72. G. Mittelman, A. Alshare, J.H. Davidson, Composite relation for laminar free convection in inclined channels with uniform heat flux boundaries, Int. J. Heat Mass Transfer 52 (2009) 4689–4694. B. Brangeon, P. Joubert, A. Bastide, Numerical investigation of natural convection in an asymmetrically heated inclined channel-chimney systems importance of the choice of artificial inlet-outlet boundary conditions, in: Proceedings of BS2013, 13th Conf. of Int. Building Performance Simulation Association, pp. 542–549. W.S. Fu, C.G. Li, C.P. Huang, J.C. Hunag, An investigation of a high temperature difference natural convection in a finite length channel without Boussinesq assumption, Int. J. Heat Mass Transfer 52 (2009) 2571–2580. P. Le Quere, C. Weisman, H. Paillere, J. Vierendeels, E. Dick, R. Becker, M. Braack, J. Locke, Modelling of natural convection flows with large temperature differences: a benchmark problem for low Mach number solvers. Part 1. Reference solutions, ESAIM: Math. Model. Numer. Anal. 39 (2005) 609–616. 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 natural convection, Int. J. Heat Mass Transfer 90 (2015) 654–664. Y. Shen, G. Zha, Simulation of Flows at All Speeds with Implicit High-Order WENO Schemes AIAA, 2009, Paper – 1312. X.S. Li, C.W. Gu, Mechanism of Roe-type schemes for all-speed flows and its application, Comput. Fluids 86 (2013) 56–70. X.S. Li, Uniform algorithm for all-speed shock-capturing schemes, Int. J. Comput. Fluid Dyn. 28 (2014) 329–338. X.S. Li, X.L. Li, All-speed Roe scheme for the large eddy simulation of homogeneous decaying turbulence, Int. J. Comput. Fluid Dyn. 30 (2016) 69–78. D. Talukdar, C.G. Li, M. Tsubokura, Numerical investigation of buoyancy-driven compressible laminar flow using new method preconditioned all-speed roe scheme, Int. Commun. Heat Mass Transfer 98 (2018) 74–84. M.J. Kermani, E.G. Plett, Roe scheme in Generalized Coordinates: Part I Formulations, AIAA, 2001, Paper-0086. M.J. Kermani, E.G. Plett, Roe scheme in Generalized Coordinates: Part IIApplication to Inviscid Viscous flows, AIAA, 2001, Paper-0087. J.M. Weiss, W.A. Smith, Preconditioning applied to variable and constant density flows, AIAA J. 33 (1995) 2050–2056. I. Abalakin, A. Dervieux, T. Atiana Kozubksaya, A vertex centered high order MUSCL scheme applying to linearized Euler acoustics, HAL, RR-4459, 2006. T.S. Chen, H.C. Tien, B.F. Armaly, Natural convection on horizontal, inclined and vertical plates with variable surface temperature or heat flux, Int. J Heat Mass Transfer 29 (1986) 1465–1478. E.M. Sparrow, L.F.A. Azevedo, Vertical channel natural convection spanning between the fully-developed limit and the single-plate boundary-layer limit, Int. J Heat Mass Transfer 28 (1985) 1847–1857. A. Bar-Cohen, W.M. Rosenhow, Thermally optimum spacing of vertical, natural convection cooled, parallel plates, J. Heat Transfer 106 (1984) 116–123. S.A. Suslov, S. Paolucci, Stability of natural convection flow in a tall vertical enclosure under non-Boussinesq conditions, Int. J. Heat Mass Transfer 38 (1995) 2143–2157. S.A. Suslov, S. Paolucci, Stability of mixed-convection flow in a tall vertical channel under non-Boussinesq conditions, J. Fluid Mech. 302 (1995) 91–115. S.A. Suslov, S. Paolucci, Nonlinear analysis of convection flow in a tall vertical enclosure under non-Boussinesq conditions, J. Fluid Mech. 344 (1997) 1–41. S.A. Suslov, S. Paolucci, Nonlinear stability of mixed convection flow under non-Boussinesq conditions. Part 1. Analysis and bifurcations, J. Fluid Mech. 398 (1999) 61–85. S.A. Suslov, S. Paolucci, Nonlinear stability of mixed convection flow under non-Boussinesq conditions. Part 2. Mean flow characteristics, J. Fluid Mech. 398 (1999) 87–108. D.R. Chenoweth, S. Paolucci, Natural convection in an enclosed vertical air layer with large horizontal temperature differences, J. Fluid Mech. 169 (1986) 173–210. A.S. Kaiser, B. Zamora, A. Viedma, Numerical correlation for natural convective flows in isothermal heated, inclined and convergent channels, for high Rayleigh numbers, Comput. Fluids 38 (2009) 1–15.