Cascaded Lattice Boltzmann Method application in forced and natural convection

Lattice Boltzmann Method is applied to forced and natural convection heat transfer from the tube banks. Hot tubes are cooled by flowing or passive air. Two Reynolds numbers (Re = 80 and Re = 1600) and two Rayleigh numbers (Ra = 103 and Ra = 105) of the corresponding heat transfer regime are studied. The method itself is based on the recently derived cascaded collision operator not only for the fluid flow but also for the temperature field. Using this method and moderate space resolution of the lattice we were able to obtain stable and bounded simulations for the non-trivial geometry.


Introduction
Lattice Boltzmann Methods (LBM) are numerical methods based on the Boltzmann Transport Equation (BTE) [1]. They solve discrete version of BTE with finite velocity set in order to recover macroscopic variables (density, momentum, temperature) from the underlying distribution functions. The dynamic of the system is described by the collision operator which can be numerically realized in different ways. The most popular are SRT (Single Relaxation Time) and MRT (Multiple Relaxation Time) operators because of their simple form and easy implementation. Cascaded operator [2] can cure some of the numerical artefacts and offers enhanced Galilean invariance of the method [3]. This was the main reason why we derived such an operator for the one conservation law and examined its behavior in case of heat transfer in the recently published article [4]. In this paper we apply cascaded version of a collision operator for D 2 Q 5 lattice and one conservation law together with D 2 Q 9 cascaded collision operator for fluid flow [2] to forced and natural convection in complex arrays of cylindrical obstacles.

Cascaded Lattice Boltzmann method
In LBM, fluid is considered as an ensemble of fictitious particles described by the velocity distribution functions (DFs) denoted as f (x, ξ, t). The evolution of such DFs obeys Boltzmann's transport equation which reads where ξ is the microscopic velocity and Ω is the collision operator. Discretization of the velocity space in the above equation by finite set of velocities c i leads to the Lattice Boltzmann Equation (written in lattice units) where f i is the velocity distribution function related to the i th characteristic velocity, K is a collision matrix, k is a vector of moments of f i and F i is a force term. The scheme implementation in the natural convection is detailed in [5]. The macroscopic density ρ and macroscopic velocity u can be obtained from f i as In this work we use two population LBM i.e. we have D 2 Q 9 lattice model for fluid flow simulation and D 2 Q 5 lattice model for temperature field simulation presented in Fig. 1, where characteristic velocities c i = (c i,x , c i,y ) are defined for both models. The weight factors, w i for D 2 Q 9 and D 2 Q 5 are w 1 = 4/9, w 3,5,7,9 = 1/9, w 2,4,6,8 = 1/36 and w 1 = 1/3, w 2,3,4,5 = 1/6, respectively and sound speed c 2 s = 1/3 for both models. In order to solve equation for the internal energy (temperature) we add another set of DFs and their evolution equation where g i is the distribution function for temperature field, Q is a collision matrix and q is a vector of moments of g i . The temperature is obtained as a sum of all g i The two matrices used in the above evolution equations are constructed from the monomials of the characteristic velocities [4] and are given as follows where superscript T denotes the transpose and the vectors k and q are given by the following relations where components k 1 , k 2 , k 3 and q 1 are equal to zero due to the conservation of the density, momentum and internal energy. The implementation of the force terms F i are described in [5]. The term κ xxyy is a central moment defined by Finally, the parameters τ (called relaxation times) are related to hydrodynamic and nonhydrodynamic modes, those related to hydrodynamic modes are used to setup the transport parameters from the macroscopic partial differential equations and the rest can be used to tune stability and accuracy of the method, in present work we keep all non-hydrodynamic relaxation times equal to unity. In lattice units (superscript "+") we have the following relations for the relaxation times and transport parameters where ν and α denote kinematic viscosity and heat diffusivity (in lattice units). To impose initial and some of the boundary conditions we need equilibrium distribution functions for each characteristic velocity, these are given as follows The initial conditions are given by zero velocity and defined temperature (in lattice units) where N x and N y are lattice dimensions of the computational domain. Boundary conditions needed during the simulations are Dirichlet for inlet, extrapolation for outlet and no-slip conditions for the velocity at walls, Dirichlet, extrapolation and adiabatic conditions for the temperature at inlet, outlet and at walls. The no-slip and adiabatic conditions can be simulated by the bounce-back approach [1] fī(x bb , t the case of Dirichlet conditions (i.e. for wall with given temperature) we use "anti-bounce-back" conditions to setup the temperature of the walls For the inlet boundary condition we use the equilibrium DF where x in/d is the location of the inlet. The outlet is modeled by a first-order extrapolation where x out is location of the outlet, n is an outer normal and i is equal to directions that are pointing to the inside of the domain. Macroscopic equations solved by the aforementioned schemes are the following [3,4] ∂u with ν and α being the kinematic viscosity and heat diffusivity of the fluid. These material parameters can be set through the relaxation times of the LBM, which are located inside the expressions for the k and q vectors. The force term is present only when the natural convection flows are solved and is given by the Boussinesq approximation with g being the vector of gravitational acceleration and β is the linear isobaric thermal expansion coefficient defined by The accuracy and stability of the presented method was studied by the authors and results can be found in the recent articles [4] (forced convection) and [5] (natural convection).

Description of the problem
We solve forced and natural convection in the vicinity of the 5x5 hot tube banks. For the case of forced convection, the tubes are inserted into the adiabatic channel and in the case of natural convection the tubes are inside the adiabatic box. The geometry of the forced convection problem is presented in Fig. 2

Numerical simulations
Numerical setup for forced convection is following; the domain is divided into L x =768 and L y =512 lattice sites, the tube diameter in lattice units is D + =32, the Mach number is set to Ma=0.01 and the rest of the parameters (ν + and α + ) are computed from the dimensionless numbers. We simulate t + =25·10 4 dimensionless timesteps. Temperature and velocity magnitude fields in lattice units are presented in Figs. 4-7.
Numerical setup for natural convection is following; the domain is divided into L x =L y =512 lattice sites, the tube diameter in lattice units is D + =32, the Mach number is set to Ma=0.01 and the rest of the parameters (ν + and α + ) are computed from the dimensionless numbers. We simulate t + =25 · 10 4 dimensionless timesteps. Temperature and velocity magnitude fields in

Discussion
Comparing the forced convection for two different Re numbers, one can observe that for lower Re number big recirculation zone exist, right behind the bank (the flow is in the laminar regime but with unstable wake), limiting the heat transfer in tthe wo last columns of the bank. The  bypass flow above and under the bank confines the temperature wake into the band area which is nearly of the height of the bank. The flow in the bank is close to laminar and for the first two columns causes better heat transfer as can be inspected from the temperature field, where last two columns are very poorly cooled. For higher Re flow (which lies in the subcritical regime), the flow and temperature patterns are more complex. One can clearly see the vortex street instabilities, and also flow pattern inside the bank is unstable. In this case the heat transfer is more intense as the flow between tubes are transitional. One recirculation zone does not exist in this case and therefore heat is removed more efficiently from the last two columns.
During the natural convection, the situation is similar, for the low Ra number, flow is laminar and slow, which affects the heat transfer from the tubes. Two recirculation zones exist at the top of the cavity on both its sides and the temperature inside the bank is equal to the tubes temperature. The heat is removed from the sides and mainly from the top of the bank as can be seen from temperature fields. For higher Ra number, the flow become unstable and the heat is removed more efficiently from the bank. The unstable wakes behind the cylinders in the upper half of the bank is clearly visible now. As in the previous case, the first few rows transport heat more efficiently than the rest of the bank (the gradients of temperatures near cylinders are bigger).

Conclusion
Lattice Boltzmann Method with cascaded collision operators for both fluid flow and heat transfer was applied to thermal flows. We simulated forced and natural convection from the hot tube banks to the cold flowing and passive air. The method was able to deliver stable and bounded results. The cascaded operator for the case of one conservation law is still not well studied and we are currently working on its deep analysis in terms of accuracy and stability. The next step is the extension of the scheme to the three dimensions and exploration its ability to simulate industrial flows.