Hybrid lattice Boltzmann––Finite difference formulation for combined heat transfer problems by 3D natural convection and surface thermal radiation

Journal Pre-proof

Hybrid lattice Boltzmann – finite difference formulation for combined heat transfer problems by 3D natural convection and surface thermal radiation A. Nee PII: DOI: Reference:

S0020-7403(19)34681-8 https://doi.org/10.1016/j.ijmecsci.2020.105447 MS 105447

To appear in:

International Journal of Mechanical Sciences

Received date: Revised date: Accepted date:

9 December 2019 11 January 2020 13 January 2020

Please cite this article as: A. Nee , Hybrid lattice Boltzmann – finite difference formulation for combined heat transfer problems by 3D natural convection and surface thermal radiation, International Journal of Mechanical Sciences (2020), doi: https://doi.org/10.1016/j.ijmecsci.2020.105447

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Elsevier Ltd. All rights reserved.

Highlights  Hybrid lattice Boltzmann method is proposed.  Developed hybrid model is successfully validated against experimental and up to date numerical data.  Interaction of 3D natural convection and surface thermal radiation is studied.  Computation speed of hybrid lattice Boltzmann is more than 19 times higher than pure finite difference approach.

Hybrid lattice Boltzmann – finite difference formulation for combined heat transfer problems by 3D natural convection and surface thermal radiation A. Nee Research and Educational Center of I.N. Butakov, National Research Tomsk Polytechnic University, Russia, Tomsk E-mail: [email protected] Abstract. In this study, a hybrid model for 3D natural convection combined with surface thermal radiation in a closed differentially heated cube was developed. Within this model, numerical procedure for fluid flow was considered in terms of the lattice Boltzmann method under Bhatnagar-Gross-Krook approximation with D3Q19 scheme. On the other hand, unsteady three-dimensional energy equation was solved by means of the finite difference technique. Developed hybrid approach was validated against experimental and up to date numerical data. Mathematical modelling was conducted for two-dimensional and threedimensional problem formulations under variation of Rayleigh number, conduction-radiation number and surface emissivity. It was found that discrepancy between 2D and 3D results was increased with an enhancement in the radiation heat transfer rate. Computational performance of hybrid lattice Boltzmann method was several times higher than conventional CFD approach. Keywords: hybrid lattice Boltzmann, finite difference method, 3D natural convection, surface thermal radiation, BGK approximation. 1. Introduction The last few decades combined heat transfer problems via natural convection and radiation received great attention. This phenomenon is encountered in many engineering applications such as photovoltaics, solar collectors, microelectronic cooling etc. Numerical studies are predominantly conducted for two-dimensional problem formulation [1–11]. Along with that computational performance of central processing units (CPU) is increased year by year. As a result, heat transfer and fluid flow regularities are extensively analyzed for 3D cavities last decade. However, it should be noted that this kind of problems is still a challenge due to significant complexity in numerical solution of unsteady three-dimensional mass, momentum and energy equations. Moreover, complexity is drastically increased when taking into account radiation. The most relevant studies in the field of 3D natural convection combined with thermal radiation are listed below. Interaction of 3D natural convection and radiation is predominantly studied by means of the differentially heated cubic cavity [12–18]. Borjini et. al. [12]

analyzed the effect of volumetric radiation on 3D laminar flow. The governing equations were formulated in terms of the vorticity – vector potential – temperature variables. It was found that radiation significantly altered the flow pattern. Wang et. al. [13] considered a similar problem formulation. However, the authors used the moving least square meshless method to examine the effect of optical thickness, albedo and heat source on heat transfer and fluid flow regularities. Cherifi et. al. [14] studied the interaction of 3D double-diffusive natural convection with thermal radiation in a cubic cavity filled by a non-gray fluid mixture. The radiation transport equation was solved by means of the discrete ordinate method. It was found that 2D and 3D results obtained in the mid-plane of the cavity represented the same trends of isotherms and isoconcentrations. Wu et. al. [15] examined the effect of surface radiation on natural convection in a 3D flat-box type cooling channel of photovoltaic thermal system with tilt angle of 30 . While the papers [12–15] analyzed the laminar flow, the studies [16–18] are devoted to turbulent natural convection with volumetric radiation at the Rayleigh number up to 3  109 . One side open 3D cylindrical cavity was considered in [19]. Wu et. al. performed numerical study of the combined natural convection and radiation heat losses. The effect of tilt angle, aspect ratio and surface emissivity on local and mean heat transfer and fluid flow characteristics was analyzed. It was found that radiation weakened the natural convection in the cavity. The effect of thermal radiation on 3D magnetohydrodynamic flow is examined in [20–24]. Al-Rashed et. al. [20] numerically studied the interaction of magnetic field and volumetric radiation in terms of the low electrically conductive dielectric oxide melt. It was found that magnetic field intensified the 3D character of the flow. When neglecting the radiation, the flow was quasi-dimensional in the area near the core of the cavity. The similar problem formulation was considered by Zhang et. al. [21]. On the other hand, papers [22–24] analyzed the effect of magnetic field and thermal radiation on 3D natural convective flow of nanofluid over a vertical stretching sheet. Conjugate 3D natural convection with radiation was considered in [25, 26]. Hachem et. al. [25] presented an immersed volume method within a finite element approach for heat transfer by conduction in a hat-shaped disk, natural convection and volumetric radiation. Martyushev and Sheremet [26] performed a numerical study of natural convection and surface thermal radiation in a cubic cavity bounded by heat-conducting finite thickness walls with a heat source. The governing

equations were formulated in terms of the vorticity – vector potential – temperature dimensionless variables. Study of combined heat transfer by 3D natural convection and thermal radiation in terms of the real engineering applications such as melting, solar collectors etc. was performed in [27–29]. Mousavi et. al. [27] examined heat transfer from a vertical finned heat sink. In order to find an optimum configuration, ten various configurations of interrupted, staggered and capped finned heat sink were studied. Heat dissipation rate by natural convection and radiation were separately analyzed for each case. Modelling of unsteady heating and melting of PWR core after SBO was performed by Zhang et. al. [28]. The effect of heat generation from a zirconium–water reaction, radiation and the adoption of an IVR strategy on heating and melting of the reactor core was analyzed. A 3D inclined rectangular solar channel was considered by Kumar and Premachandran [29]. The aim of the work was to understand how atmospheric wind affected the heat transfer in natural convection based solar air heaters. The above studies were performed in terms of the conventional computational fluid dynamics (CFD) techniques. On the other hand, alternative tool for CFD problems called the lattice Boltzmann method (LBM) is extensively developing by now. The key advantages of LBM over traditional numerical techniques based on Navier-Stokes equations solution are numerical simplicity and natural parallelism. The most commonly used thermal LBM model [30–32] is double distribution function (DDF). A little attention is paid to another thermal LBM model called hybrid approach [33] where energy equation is solved by means of the conventional numerical techniques. This approach has numerous advantages over DDF such as higher numerical stability [34], less memory consumption etc. To the best of my knowledge, very little papers are published where combined heat transfer is analyzed in terms of the hybrid LBM. Thus, the motivation of this study is to build a numerically efficient hybrid LB-FDM model for combined heat transfer problems via 3D natural convection and thermal surface radiation. Within this model, fluid dynamics will be solved by means of the LBM with D3Q19 scheme and Bhatnagar-Gross-Krook approximation and energy equation will be solved in terms of the finite difference technique. 2. Problem formulation 2.1. Geometrical and physical models

Fig. 1. Differentially heated cubic cavity filled with radiatively non-participating medium A closed cube filled with radiatively non-participating medium is considered. Heat transfer via radiation was performed only between grey diffusive walls of the cavity. It was assumed that thermophysical properties of the fluid were temperature-independent. Unsteady laminar 3D flow was considered. The boundary condition of the first kind was applied at the opposite vertical walls. Other boundaries assumed to be heat insulated. No-slip condition for flow was set at the solid walls. It was assumed that the effects of viscous energy dissipation were negligible. 2.2. Mathematical model The discretized lattice Boltzmann equation for incompressible low Mach number fluid flow with the Bhatnagar-Gross-Krook (BGK) approximation is as follows:

f k  x  ck  t , t  t   f k  x , t   t     f keq  x , t   f k  x , t   t  F .

(1)

Fig. 2. The lattice velocities for three-dimensional nineteen velocity (D3Q19) scheme. The lattice velocities for D3Q19 model (Fig. 2.) are given by: 0, k  1,  ck   1,0,0  ,  0,0, 1 ,  0, 1,0  , k  2..7,  k  8..19.  1,0, 1 ,  1, 1,0  ,  0, 1, 1 ,

The equilibrium distribution function is as follows:

f

eq k

 3   ck  u  9   ck  u 2 3   u  u    wk    1    , 2 2 2 c 2  c 2  c  

(2)

where wk is weighting factor defined as:

1  3 ; k  1,  1 wk   ; k  2..7, 18 1  36 ; k  8..19.  Macroscopic density and velocity are recovered as: 19

   fk .

(3)

k 1

u

1



19

  ck  f k . k 1

Relaxation frequency is calculated as:

(4)



1 , 3   0.5

(5)

Force term in equation (4) under Boussinesq approximation can be determined as follows: T T     g   T  h c   cz   c  u  2   F  f eq , R T

(6)

Transient energy equation in dimensionless variables for natural convection under Boussinesq approximation is as follows:   2   2   2      1 U V W   2  2 .  2 τ X Y Z Y Z  Ra  P r  X

(7)

In order to discretize the equation (7), Euler scheme for temporal term, monotonic scheme of Samarskii for convective terms and central differences for diffusive terms were used. After approximation the following equations can be obtained:



1 3 i , j ,k n



n i , j ,k

τ

 U in, j ,k 



1 3 i 1, j ,k n



1 3 i 1, j , k n

2  hx

 U in, j ,k 



1 3 i 1, j , k n

 2

1 3 i , j ,k n



1 3 i 1, j , k n

2  hx

 (8)

1 1 1 n n n   1 3 3 3    2     1 Ra  Pr  h  i 1, j , k i , j ,k i 1, j , k  n x  1  U    i , j , k ,  2 2 Ra  Pr  hx   



2 3 i , j ,k n



1 3 i , j ,k n

τ

 Vi ,nj ,k 



2 3 i , j 1,k n



2 3 i , j 1,k n

2  hy

 Vi ,nj ,k 



2 3 i , j 1, k n

 2

2 3 i , j ,k n



2 3 i , j 1, k n

2  hy

 (9)

2 2 2 1 n n n   3 3 3     2    Ra  Pr  h 1  i , j 1,k i , j ,k i , j 1,k  y   1  Vi ,nj ,k   2 .  2 h Ra  Pr  y     



n 1 i , j ,k

 τ

2 3 i , j ,k n

W

n i , j ,k



in,j1,k 1  in,j1,k 1 2  hz

 1 Ra  Pr  hz   1  Wi ,nj ,k   2 Ra  Pr  

W

n i , j ,k

1

   



in,j1,k 1  2  in,j1,k  in,j1,k 1

n 1 i , j ,k 1

2  hz n 1 i , j ,k 2 z

 2 h



n 1  i , j ,k 1

.  

 (10)

The equations (8) – (10) were reduced to the standard tridiagonal form and solved by means of the sweep method. In order to reduce equation (7) to its dimensionless form, the following relations were used: 

T  T0 t u v w x y z ; X  ;Y  ; Z  ; U  ; V ; W ; ; t0 Vnc Vnc Vnc Th  Tc L L L Vnc  g    Th  Tc   L ; t0 

L . Vnc

The net flux method was used to compute the surface thermal radiation. Within this method, surfaces of the cavity were divided into isothermal areas, where the radiation characteristics would be constant. Then, radiosity and irradiation were calculated. Dimensionless radiosity is determined as follows: N

J i   i  i 4  1   i   i  j  J j ,

(11)

j 1

where J i is non-dimensional radiosity from the i element,  i  4 is self-radiation of N

the i element,

1   i   i j  J j is

reflected radiant energy coming from

j 1

surrounding j elements that can be seen by the i element,  is view factor. The irradiation is given by: N

Qr ,i  J i  i  j  J j .

(12)

j 1

The following relations were used to reduce radiosity and irradiation to its non-dimensional form: q' q '' J ; Q  .   Th4 r   Th4 View factors for cubic cavity were determined by means of the analytical relations and graphs presented in [35]. Radiosity equation (11) was solved in terms of the Gauss elimination method. Equations (11), (12) were recomputed at every time step since transient problem formulation was considered. 2.3. Initial and boundary conditions Initial conditions for equations (1), (7) are as follows: f  X ,Y,Z ,0   f eq  X ,Y,Z ,0   0,

  X ,Y,Z ,0  U  X ,Y,Z ,0   V  X ,Y,Z ,0   W  X ,Y,Z ,0   0

  X ,Y,Z ,0   1. Boundary conditions for equations (7) are as follows:

X  0, X  1, 0  Y  1, 0  Z  1:

  N r  Qr  0, X

Z  0, Z  1, 0  X  1, 0  Y  1:

  N r  Qr  0, Z

Y  0, 0  X  1, 0  Z  1:  0.5, Y  1, 0  X  1, 0  Z  1:  0.5, In order to determine distribution function at the solid walls, a bounce back condition [36, 37] was used. 2.4. Grid study and validation The following convergence criteria were used to obtain steady-state solution:

U

2

V 2 W 2   n

U

2

V 2 W 2 

n 1

 107 ,

n  n1  107. Before the main numerical procedure, grid refinement test was conducted. It should be noted that numerical stability was achieved with the grid size of 273 when Ra=106 , Nr =150 and  =0.5 . However, only three cases will be presented in order to avoid too dense graph. Figure 3 shows results of grid study.

b a c Fig. 3. Steady-state temperature (a), velocity component along Y-axis (b) and Zaxis (с) in the planes of X=0.5 (section of Y=0.5) and Y=0.5 (section of Z=0.5) when Ra=106 , Nr =150 and  =0.5 . It was found that the deviation in temperature and velocity components did not exceed 1% with a number of 613 nodes. Thus, mathematical modelling was performed on a uniform 613 mesh. Three-dimensional hybrid LB-FDM model was validated against available experimental and up to date numerical data. Temperature along X-axis in the mid-

plane was used as a characteristic criterion. It should be noted that numerical data of other researchers were obtained with 2D models. Figure 4 presents temperature profiles when Ra=3.55 105 and   0.995 .

Fig. 4. Comparison of temperature profiles in the mid-plane obtained numerically with 2D and 3D simulations and experimentally by Bajorek and Lloyd [38] when Ra=3.55 105 and   0.995 . When comparing numerical results with benchmark experimental data, a reasonable conclusion could be made that both 2D and 3D simulation in a satisfactory agreement with experimental study of Bajorek and Lloyd [38] in terms of the temperature. But still, it was clearly seen that 3D modelling produced more accurate results. Higher values of temperature in the center of the cavity for 2D simulation were associated with the higher radiation heat transfer rate in this case of simulation. 3. Results and discussion Combined heat transfer by natural convection and surface thermal radiation was numerically studied in terms of the cube filled with air. The Prandtl number was set to 0.71. Variation range of the Rayleigh number was 104  Ra  106 . Surface emissivity was varied in a range of 0    1. Conduction-radiation number was changed from 50 to 150. The results of mathematical modelling are presented in terms of the isosurfaces, streamslices, temperature contours, temperature and velocity profiles. The main attention was given to comparison of 2D and 3D

problem formulations, computation time of developed hybrid approach and pure finite difference technique. 3.1. Temperature and flow Steady-state isosurfaces and streamslices with variation of the Rayleigh number are presented in Fig. 5.

a

b

c

f e d Fig. 5. Steady-state isosurfaces (a, b, c) and streamslices (d, e, f) with variation of Rayleigh number and conduction-radiation number when  =0.5 : a, d) Ra=104 , N r =50 ; b, e) Ra=105 , N r =90 ; c, f) Ra=106 , Nr =150 When Ra=104 (Fig. 5 d), a uniform flow pattern in the mid-plane was formed in the cube. However, the cores of convective cells near the left and right walls were oppositely displaced from the center to the hot and cold surfaces. This flow behavior was probably associated with the additional walls heating by surface thermal radiation. Moreover, this factor led to isosurfaces bending (Fig. 5 a) near the walls. With an increment in the buoyancy force to 105 the fluid velocity was increased. As a result, inner flows (Fig. 5 e) inside a large-scale convective cell were formed in the mid-plane. Along with that the cores of vortices near the left and right walls were advanced to the corners of the cube. Thermal stratification (Fig. 5 b) could be noted in the center of the cavity. Further increase in the

Rayleigh number to 106 led to an enhancement of natural convective flows. A complex flow pattern (Fig. 5 f) was formed in the mid-plane and the cores of convective cells near the left and right walls achieved the corners of the cavity. Moreover, an oblique stratification was observed in the cube due to surface thermal radiation. Comparative study of 2D and 3D results is presented below. Vorticity – stream function – temperature formulation was used for two-dimensional simulation. Numerical procedure is described in detail in previous works [41, 42]. Figure 6 shows temperature contours when varying the Rayleigh number.

a

b

c

f e d Fig. 6. Steady-state temperature contours with variation of Rayleigh number and conduction-radiation number when  =0.5 for 2D (a, b, c) and 3D (d, e, f) problem formulations in the mid-plane: a, d) Ra=104 , N r =50 ; b, e) Ra=105 , N r =90 ; c, f)

Ra=106 , Nr =150 It was clearly seen that deviations of 2D and 3D temperature in the midplane were increased with an increment in the buoyancy force and conductionradiation number. In the conditions under study, variation of the Rayleigh number was performed in terms of the characteristic size of the cavity. Thus, an increase in the cavity length also led to an increment in the conduction-radiation number. Surface thermal radiation was enhanced. When Ra=104 and N r =50 (Fig. 6 a, d),

2D and 3D temperature contours in the mid-plane were hardly distinguishable. Probably, this factor was associated with a relative low radiation heat transfer rate. An increase in the buoyancy force to 105 and conduction-radiation number to 90 led to a discrepancy in isotherms shape in the center of the cavity (Fig. 6 b, e). The angle of thermal stratification was steeper in the case of the 2D problem formulation. This regularity was due to higher intensity of surface radiation since view factors for 2D formulation was greater. Along with that enhancement of surface radiation weakened natural convective flow along the top and bottom walls. Hence, temperature was higher near the top wall and lower near the bottom wall in the case of 3D formulation. The same regularity was observed when

Ra=106 and Nr =150 . Moreover, discrepancy in the 2D and 3D was more pronounced in this case of simulation (Fig. 6 c, f). To sum up, errors in heat transfer and fluid flow characteristics were increased with an enhancement of the surface thermal radiation since view factors were different for 2D and 3D problem formulations. Probably, the contribution of buoyancy force in total discrepancy of 2D and 3D results was minor. Figure 7 shows temperature and velocity profiles with and without surface thermal radiation in the mid-plane and section of Y=0.5.

a

b

c

d

Fig. 7. Temperature (a, b) and velocity component (c, d) profiles with variation of surface emissivity for 2D and 3D problem formulations in the mid-plane and section of Y=0.5 when Ra=106 and Nr =150 : a, c)  =0 ; b, d)  =0.2 It was interesting to note that temperature for 2D and 3D simulations was in a good agreement in the case of no radiation (Fig. 7 a). However, a discrepancy was found for velocity component (Fig. 7 c). This discrepancy was probably associated with the flow evolution simultaneously along three coordinates whereas one of the coordinate was excluded in the case of 2D modelling. On the other hand, the similar velocity component trends were observed for 2D and 3D results (Fig. 7 c). When taking into account the surface thermal radiation, error in dimensionless temperature near the walls (Fig. 7 b) approximately achieved 25 %. Moreover, a significant deviation in velocity component (Fig. 7 d) was found. Summarizing the above, 2D and 3D temperature satisfactory correlated whereas velocity component profiles had different trends with a relatively weak surface thermal radiation rate. 3.2. Heat transfer rate Integral analysis of heat transfer was performed in terms of the mean Nusselt numbers. Convective and radiative heat transfer rate was computed for the vertical isothermal hot wall. The mean convective, radiative and total Nusselt numbers are given by:  dX  dZ ,  Y Y 0 0 0

1 1

NuC   1 1

NuR    N r  Qr  Y 0 dX  dZ , 0 0

NuT  NuC  NuR .

Figure 8 presents variation of mean Nusselt numbers with surface emissivity and Rayleigh number obtained with 2D and 3D simulation.

a

b

c Fig. 8. Variation of convective (a), radiative (b) and total (c) Nusselt numbers with surface emissivity and Rayleigh number obtained with 2D and 3D models As could be assumed, buoyancy driven flow (Fig. 8 a) was enhanced with an increment in the Rayleigh number. On the other hand, it was found that an increase in the surface emissivity led to a decrease in the convective Nusselt numbers for both 2D and 3D problem formulations. However, more significant reduction in convective mechanism of heat transfer was observed in the case of twodimensional simulation. This regularity was probably associated with the higher values of view factors for 2D model. Thus, the heating rate of the top wall was higher. As a result, velocity of an up-ward flow near the hot wall was reduced more significantly in the case of 2D simulation. Therefore, mean convective Nusselt numbers obtained with 3D model had the higher values in comparison with 2D problem formulation when taking into account surface thermal radiation. On the contrary, radiative Nusselt number (Fig. 8 b) was increased with an increment in the surface emissivity and higher values of NuR were observed in the case of 3D simulation. This regularity associated with the enhancement of radiative heat transfer when increasing the surface emissivity. Two-dimensional model produced lower values of NuR since radiation heat transfer was performed between boundaries. On the other hand, radiative energy was transferred between planes when taking account 3D problem formulation. As a result, irradiation for the hot wall was higher in the case of 3D simulation. However, it should be noted that radiative Nusselt number for the top and bottom walls was higher for 2D model since view factors had the higher values in comparison with 3D model. When analyzing overall heat transfer rate (Fig. 8 c), a positive trend was observed with variation of surface emissivity. Weakening of convective mechanism of heat transfer with an increment in the  was counterbalanced by enhancement of radiation heat transfer.

3.3 Computational performance In this sub-section, hybrid lattice Boltzmann method will be compared with a conventional CFD approach. Computational performance will be discussed in terms of the running time. Vorticity – vector potential – temperature formulation [44] was used for pure finite difference method. In order to approximate the governing equations, locally one-dimensional scheme of Samarskiiy [45] was applied. Convective terms were discretized by the monotonic scheme of Samarskii [45] whereas diffusive terms were approximated by the central differences. In general, numerical procedure is quite similar to [41, 42]. Intel core i5-4400 with frequency of 3,1 GHz was used for all computations. Figure 9 shows computation time of one time step versus grid size.

Fig. 9. Variation of computation time of one time step with grid size in terms of the hybrid Lattice Boltzmann and pure finite difference technique when Ra=106 , Nr =150 and   0.995 . It was found that computational performance of hybrid LB-FDM model was significantly better than pure finite difference technique. In the case of pure FDM, computation time of one time step was exponentially increased with an increment in the grid size. On the other hand, hybrid LB-FDM approach had an approximately linear trend. It was interesting to note that the more was mesh points the clearer was the difference in computational performance of hybrid lattice Boltzmann method versus pure FDM. The running times necessary to obtain a steady-state solution for various Rayleigh numbers are presented in Table 1.

Convergence criteria for temperature and velocity were set to 107 for both hybrid LB-FDM and pure FDM approaches. Numerical codes were developed in the MatLab 2017a. The running times were estimated with the tic-toc function. Table 1. Variation of running times with the Rayleigh number Ra Running time, s Pure FDM Hybrid LB-FDM 3 34383.8 891 10 49036.32 1675.39 104 5 72108.31 3715.41 10 It was clearly seen that hybrid lattice Boltzmann method consumed much less time to obtain steady-state solution in comparison with pure finite difference technique in terms of the 3D natural convection combined with surface thermal radiation. In general, this scenario was predictable since only one partial differential equation for flow was solved in LBM instead of six equations in the case of FDM. As could be assumed, an increment in the Rayleigh number led to an increase in the running time. This regularity was associated with an enhancement of natural convective mechanism of heat transfer. Fluid flow was intensified. As a result, more time was required to achieve the convergence criteria. To sum up, hybrid lattice Boltzmann method surpassed conventional CFD approach approximately more than 19 times in terms of the running time. Conclusion A hybrid lattice Boltzmann – finite difference model was developed for 3D natural convection combined with surface thermal radiation. Within this model, fluid flow was solved by means of the LB BGK approximation with D3Q19 scheme whereas heat transfer was considered in terms of the finite difference technique. The main findings are as follows:  When taking into account the surface thermal radiation, an oblique stratification was formed in the cubic cavity. Moreover, opposite flows were observed near the left and right walls.  In the case of no radiation, temperature and velocity component obtained with 2D and 3D simulation were in a satisfactory agreement. On the other hand, the discrepancy between 2D and 3D results was increased with an enhancement in the radiation heat transfer rate. Error in temperature near the walls approximately achieved 25 % when surface emissivity was equal to 0.2. It should be noted that the contribution of buoyancy force in total discrepancy of 2D and 3D results was probably minor.

 Radiation reduced mean convective Nusselt numbers. However, overall heat transfer rate was higher since weakening of natural convection with an increment in the surface emissivity was counterbalanced by enhancement of radiation heat transfer.  Hybrid lattice Boltzmann method had an approximately linear trend of computation time of one time step with variation of grid size. Moreover, LBFDM model was required much less time to obtain steady-state solution in comparison with pure finite difference technique. When Ra=103 , Ra=104 and Ra=105 , computational performance of hybrid lattice Boltzmann method was correspondingly 38, 29 and 19 times higher than conventional CFD approach. Acknowledgments The research is carried out at Tomsk Polytechnic University within the framework of Tomsk Polytechnic University Competitiveness Enhancement Program. In addition, the author would like to thank the reviewers for their valuable comments and suggestions.

Nomenclature a : Thermal diffusivity, m2 /s ; c : Lattice speed, m/s ck : Particle speeds, m/s

x, y, z; T : Temperature, K; T0 : Initial temperature, K;

f: Distribution function, kg/m3 ;

Th : Temperature of the hot wall, K;

g : Gravitational acceleration, m/s 2 ;

Tc : Temperature of the cold wall, K;

J : Dimensionless radiosity; L : Length of the cavity; N r : Conduction-radiation number=

wk: Weighting factor; Greek symbols  : Coefficient of thermal expansion,

σ  T

K -1 ;  : Solid-fluid interfaces emissivity;  : View factor;

4

h

X ,Y , Z : Dimensionless analogues of

 L  /  λ  Th  Tc   ; Mean

convective

Nusselt

number; Mean NuR :

radiative

Nusselt

NuC :

number; NuT : Mean total Nusselt number;

Pr: Prandtl number= / a ; '

q : Radiosity, W/m ;

 : Stefan–Boltzmann

q '' : Irradiation, W/m2 ; Qr : Dimensionless irradiation;

Rayleigh

 g    (T

 : Kinematic viscosity, m/s 2 ;

 : Density, kg/m3 ;

2

Ra:

 : Dimensionless temperature;  : Heat conductivity coefficient, W/(m  K) ;

number=

 Tc )  L3  /   a  ;

constant,

W/  m2  K 4  ;

 : Dimensionless time;  : Dimensional time step in energy

t : Time, s;

equation;  : Relaxation frequency, 1/s;

t0 : Time scale, s;

 : Collision operator, kg/(s  m3 ) ;

h

t : Time step in lattice Boltzmann equation; u, v, w : Velocities on the x and y axes, respectively, m/s; U ,V ,W : Dimensionless analogues of u, v, w ;

Vnc : Velocity scale, m/s;

x, y, z : Cartesian coordinates, m;

Every part of the work entitled ―Hybrid lattice Boltzmann – finite difference formulation for combined heat transfer problems by 3D natural convection and surface thermal radiation‖ was solely performed by Alexander Nee. Declaration of interests

xThe author declares that he has no known competing financial interests or

personal relationships that could have appeared to influence the work reported in this paper.

References 1. Akiyama, M., and Chong, Q.P., 1997, ―Numerical analysis of natural convection with surface radiation in a square enclosure,‖ Numer. Heat Tr. A-Appl., 32, pp. 419–433. http://dx.doi.org/10.1080/10407789708913899 2. Ridouane, E.H., Hasnaoui, M., and Campo, A., 2006, ―Effects of surface radiation on natural convection in a rayleigh-benard square enclosure: steady and unsteady conditions,‖ Heat Mass Transfer, 42, pp. 214 – 225. https://doi.org/10.1007/s00231-005-0012-7 3. Vivek, V., Sharma, A.K., and Balaji, C., 2012, ―Interaction effects between laminar natural convection and surface radiation in tilted square and shallow enclosures,‖ Int. J. Therm. Sci., 60, pp. 70–84. https://doi.org/10.1016/j.ijthermalsci.2012.04.021 4. Bouafia, M., Hamimid, S., and Guellal, M., 2015, ―Non-Boussinesq convection in a square cavity with surface thermal radiation,‖ Int. J. Therm. Sci., 96, pp. 236– 247. https://doi.org/10.1016/j.ijthermalsci.2015.04.017 5. Ibrahim, A., Saury, D., and Lemonnier, D., 2013, ―Coupling of turbulent natural convection with radiation in an air-filled differentially-heated cavity at Ra = 1.5×109,‖ Comput. Fluids., 88, pp. 115–125. https://doi.org/10.1016/j.compfluid.2013.09.006 6. Saravanan, S., and Sivaraj, C., 2015 ―Combined natural convection and thermal radiation in a square cavity with a nonuniformly heated plate,‖ Comput. Fluids., 117, pp. 125–138. https://doi.org/10.1016/j.compfluid.2015.05.005 7. Ashish Gad, M., and Balaji, C., 2010, ―Effect of surface radiation on RBC in cavities heated from below,‖ Int. Commun. Heat. Mass., 37, pp. 1459–1464. https://doi.org/10.1016/j.icheatmasstransfer.2010.08.003 8. Sharma, A.K., Velusamy, K., Balaji, C., and Venkateshan, S.P., 2007, ―Conjugate turbulent natural convection with surface radiation in air filled

rectangular enclosures,‖ Int. J. Heat. Mass. Tran., 50, pp. 625–639. https://doi.org/10.1016/j.ijheatmasstransfer.2006.07.022 9. Xaman, J.P., Hinojosa, J.F., Flores, J.J., and Cabanillas, R.E., 2008, ―Effect of the surface thermal radiation on turbulent natural convection in tall cavities of facade elements,‖ Heat Mass Transfer, 45, pp. 177–185. https://doi.org/10.1007/s00231-008-0393-5 10. Lari, K., Baneshi, M., Gandjalikhan Nassab, S.A., Komiya, A., and Maruyama, S., 2011, ―Combined heat transfer of radiation and natural convection in a square cavity containing participating gases,‖ Int. J. Heat. Mass. Tran., 54, pp. 5087–5099. https://doi.org/10.1016/j.ijheatmasstransfer.2011.07.026 11. Mahian, O., Kolsi, L., Amani, M., Estellé, P., Ahmadi, G., Kleinstreuer, C., Marshall, J.S., Taylor, R.A., Abu-Nada, E., Rashidi, S., Niazmand, H., Wongwises, S., Hayat, T., Kasaeian, A., Pop, I., 2019, ―Recent advances in modeling and simulation of nanofluid flows—Part II: Applications,‖ Physics Reports, 791, pp. 1–59. https://doi.org/10.1016/j.physrep.2018.11.003 12. Naceur Borjinia, M., Ben Aissia, H., Halouani, and K., Zeghmati, B., 2008, ―Effect of radiative heat transfer on the three-dimensional Buoyancy flow in cubic enclosure heated from the side,‖ Int. J. Heat. Fluid. Fl., 29, pp. 107–118. https://doi.org/10.1016/j.ijheatfluidflow.2007.07.012 13. Wang, C.-A., Sadat, H., and Tan, J.Yu., 2017, ―Simulations of coupled radiative and convective heat transfer in 2D and 3D semitransparent medium by the moving least square method,‖ Numer. Heat Tr. B-Fund., 72, pp. 33–54. http://dx.doi.org/10.1080/10407790.2017.1338083 14. Cherifi, M., Laouar-Meftah, S., Benbrik, A., Lemonnier, D., and Saury, D., 2016, ―Interaction of radiation with double-diffusive natural convection in a threedimensional cubic cavity filled with a non-gray gas mixture in cooperating cases,‖ Numer. Heat Tr. A-Appl., 69, pp. 479–496. https://doi.org/10.1080/10407782.2015.1090233 15. Wu, S.-Y., Wu, Y.-Y., and Xiao, L., 2014, ―Effect of thermal radiation on convection heat transfer in cooling channel of photovoltaic thermal system,‖ Journal of Mechanical Science and Technology, 28, pp. 3353–3360. https://doi.org/10.1007/s12206-014-0746-6 16. Wu, T., and Lei, C., 2015, ―On numerical modelling of conjugate turbulent natural convection and radiation in a differentially heated cavity,‖ Int. J. Heat. Mass. Tran., 91, pp. 454–466. https://doi.org/10.1016/j.ijheatmasstransfer.2015.07.113 17. Soucasse, L, Rivièrea, Ph., and Soufiani, A., 2014, ― Subgrid-scale model for radiative transfer in turbulent participating media,‖ J. Comput. Phys., 257, pp. 442–459. http://dx.doi.org/10.1016/j.jcp.2013.10.006 18. Soucasse, L, Rivièrea, Ph., and Soufiani, A., 2016, ―Natural convection in a differentially heated cubical cavity under the effects of wall and molecular gas radiation at Rayleigh numbers up to 3 ×10 9,‖ Int. J. Heat. Fluid. Fl., 61, pp. 510– 530. http://dx.doi.org/10.1016/j.ijheatfluidflow.2016.06.012 19. Wu, S.-Y., Guo, F.-H., and Xiao, L., 2014, ―Numerical investigation on combined natural convection and radiation heat losses in one side open cylindrical

cavity with constant heat flux,‖ ,‖ Int. J. Heat. Mass. Tran., 71, pp. 573–584. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2013.12.064 20. Al-Rashed, A.A.A.A., Kolsi, L., Tashkandi, M.A., Hasani Malekshah, E., Chamkha, A.J., and Borjini, M.N., 2019, ―Three-dimensional combined radiationmagnetoconvection of low electrically conductive dielectric oxide melt,‖ Int. J. Numer. Method H., 29, pp. 3611–3637. https://doi.org/10.1108/HFF-06-2018-0263 21. Zhang, J.-K, Li, B.-W., Dong, H., Luo, X.-H, and Lin H., 2017, ―Analysis of magnetohydrodynamics (MHD) natural convection in 2D cavity and 3D cavity with thermal radiation effects,‖ Int. J. Heat. Mass. Tran., 112, pp. 216–223. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.04.105 22. Ghadikolaei, S.S., and Gholinia, M., 2019, ―Terrific effect of H2 on 3D free convection MHD flow of C2H6O2–H2O hybrid base fluid to dissolve Cu nanoparticles in a porous space considering the thermal radiation and nanoparticle shapes effects,‖ Int. J. Hydrogen Energ., 44, pp. 17072–17083. https://doi.org/10.1016/j.ijhydene.2019.04.171 23. Nayak, M.K., Akbar, N.S., Pandey, V.S., Khan, Z.H., Tripathi, D., 2017, ―3D free convective MHD flow of nanofluid over permeable linear stretching sheet with thermal radiation,‖ Powder Technol., 315, pp. 205–215. https://doi.org/10.1016/j.powtec.2017.04.017 24. Ghadikolaei, S.S., and Gholinia, M., Hoseini, M.E., and Ganji, D.D., 2019, ―Natural convection MHD flow due to MoS 2 –Ag nanoparticles suspended in C2H6O2–H2O hybrid base fluid with thermal radiation,‖ J. Taiwan Inst. Chem. E., 97, pp. 12–23. https://doi.org/10.1016/j.jtice.2019.01.028 25. Hachem, E., Digonnet, H., Massoni, E., and Coupez, T., 2012, ―Immersed volume method for solving natural convection, conduction and radiation of a hatshaped disk inside a 3D enclosure,‖ Int. J. Numer. Method H., 22, pp. 718–741. https://doi.org/10.1108/09615531211244871 26. Martyushev, S.G., and Sheremet, M.A., 2014, ―Conjugate natural convection combined with surface thermal radiation in a three-dimensional enclosure with a heat source,‖ Int. J. Heat. Mass. Tran., 73, pp. 340–353. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2014.02.009 27. Mousavi, H., Rabienataj Darzi, A.A., Farhadi, M., and Omidi, M., 2018, ―A novel heat sink design with interrupted, staggered and capped fins,‖ Int. J. Therm. Sci., 127, pp. 312–320. https://doi.org/10.1016/j.ijthermalsci.2018.02.003 28. Zhang, X., Chen, H., Liu, F., and Zhan, D., 2018, ―Transient 3D simulation for heating and melting process of PWR core after SBO,‖ Ann. Nucl. Energy, 115, pp. 256–267. https://doi.org/10.1016/j.anucene.2018.01.019 29. Kumar, D., and Premachandran, B., 2019, ―Effect of atmospheric wind on natural convection based solar air heaters,‖ Int. J. Therm. Sci., 138, pp. 263–275. https://doi.org/10.1016/j.ijthermalsci.2018.12.010 30. Sharma, K.V., Straka, R., Tavares, F.W., 2017, ―New Cascaded Thermal Lattice Boltzmann Method for simulations of advection-diffusion and convective heat transfer,‖ Int. J. Therm. Sci., 118, pp. 259–277. https://doi.org/10.1016/j.ijthermalsci.2017.04.020

31. Sajjadi, H., Amiri Delouei, A., Sheikholeslami, M., Atashafrooz, M., Succi, S., 2019, ―Simulation of three dimensional MHD natural convection using double MRT Lattice Boltzmann method,‖ Physica A, 515, pp. 474–496. https://doi.org/10.1016/j.physa.2018.09.164 32. Frapolli, N., Chikatamarla, S.S., Karlin, I.V., 2018, ―Entropic lattice Boltzmann simulation of thermal convective turbulence,‖ Comp. Fluids, 175, pp. 2–19. https://doi.org/10.1016/j.compfluid.2018.08.021 33. Bettaibi, S., Kuznik, F., Sediki, E., 2014, ―Hybrid lattice Boltzmann finite difference simulation of mixed convection flows in a lid-driven square cavity,‖ Phys. Lett. A, 378, pp. 2429–2435. https://doi.org/10.1016/j.physleta.2014.06.032 34. Lallemand, P., Lou, L.-S., 2003, ―Hybrid finite-difference thermal lattice Boltzmann equation,‖ Int. J. Modern Physics B., 17, pp. 41–47. 35. Sparrow, E. M., and Cess, R. D., 1970, Radiation Heat Transfer, Brook/Cole Publishing Company, Belmont, California. 36. Peng, Y., Shu, C., Chew, Y.T., 2003, ―A 3D incompressible thermal lattice Boltzmann model and its application to simulate natural convection in a cubic cavity‖, J. Comput. Phys., 193, pp. 260–274. https://doi.org/10.1016/j.jcp.2003.08.008 37. Mezrhab, A., Moussaoui, M.A., Jami, M., Naji, H., Bouzidi, M., 2010, ―Double MRT thermal lattice Boltzmann method for simulating convective flows,‖ Phys. Lett. A., 374, pp. 3499–3507. https://doi.org/10.1016/j.physleta.2010.06.059 38. Bajorek, S.M., Lloyd, J.R., 1982, ―Experimental investigation of natural convection in partitioned enclosures,‖ J. Heat Transf., 104, pp. 527–532. https://doi.org/ 10.1115/1.3245125 39. Patil, S, Sharma, A. K., Velusamy, K., 2016, ―Conjugate laminar natural convection and surface radiation in enclosures: Effects of protrusion shape and position,‖ Int. Commun. Heat Mass., 76, pp. 139–146. http://dx.doi.org/10.1016/j.icheatmasstransfer.2016.05.006 40. Nia, M. F., Nassab, S. A. G., Ansari, A. B., 2018, ―Transient combined natural convection and radiation in a double space cavity with conducting walls,‖ Int. J. Therm. Sci., 128, pp. 94–104. https://doi.org/10.1016/j.ijthermalsci.2018.01.021 41. Kuznetsov, G.V., Kurilenko, N.I., and Nee, A.E., 2018, ―Mathematical modelling of conjugate heat transfer and fluid flow inside a domain with a radiant heating system,‖ Int. J. Therm. Sci., 131, pp. 27–39. https://doi.org/10.1016/j.ijthermalsci.2018.05.010 42. Kuznetsov, G.V., and Nee, A.E., 2018, ―Conduction, convection, and radiation in a closed cavity with a local radiant heater,‖ Frontiers in Heat and Mass Transfer, 10, pp. 1–9. https://doi.org/10.5098/hmt.10.26 43. Wang, H., Xin, S., and Le Quere, P., 2006, ―Etude numerique du couplage de la convection naturelle avec le rayonnement de surfaces en cavite carree remplie d'air,‖ Comptes Rendus Mecanique, 334, pp. 48–57. https://doi.org/10.1016/j.crme.2005.10.011

44. N. S. Gibanov, M. A. Sheremet, Natural convection in a cubical cavity with different heat source configurations, Thermal Science and Engineering Progress 7 (2018) 138–145. https://doi.org/10.1016/j.tsep.2018.06.004 45. Samarskii, A.A., 1977, The theory of difference schemes [in Russian], Nauka, Moscow.

Graphical Abstract