Numerical Investigation of Conjugate Heat Transfer and Natural Convection Using the Lattice-Boltzmann Method for Realistic Thermophysical Properties

Abstract

To enable the lattice-Boltzmann method (LBM) to account for temporally constant but spatially varying thermophysical properties, modifications must be made. Recently, many methods have emerged that can account for conjugate heat transfer (CHT). However, there still is a lack of information on the possible physical property range regarding realistic properties. Therefore, two test cases were investigated to gain further insight. First, a differentially heated cavity filled with blocks was used to investigate the influence of CHT on the error and stability of the LBM simulations. Reference finite volume method (FVM) simulations were carried out to estimate the error. It was found that a range between $0.5$ to $1.5$ is recommended for the fluid relaxation time to balance computational effort, stability, and accuracy. In addition, realistic thermophysical properties of fluids and solids were selected to test whether the lattice-Boltzmann method is suitable for simulating relevant industry-related applications. For a stable simulation, a mesh with 64 times more lattices was needed for the most extreme test case. The second test case was an insulated cavity with a heating pad as the local heat source, which was investigated in terms of the accuracy of a transient simulation and compared to a FVM simulation. It was found that the fluid-phase relaxation time mainly determines the error and that large thermal relaxation times for the solid improve accuracy. Observed deviations from the FVM reference simulations ranged from approximately 20% to below 1%, depending on collision operator and combination of relaxation times. For processes with a large temperature spread, the temporally constant thermophysical properties of the LBM are the primary constraint.

1. Introduction

Coupled convective (forced or natural) heat transport in fluids and heat conduction in solids is known as conjugate heat transfer (CHT) [1] and is present in almost all technical processes. The applications are almost countless and cover a wide range of industries. For example, the cooling of electronic components [2,3,4,5], which must provide more and more power in ever-smaller spaces, and efficiency improvements cannot fully compensate for the increasing heat generation. In addition to forced convection, natural convection plays a major role, for example, in devices intended for domestic use because fans should be avoided due to their undesirable noise emission and power consumption.

CHT affects the entire energy sector, both conventional and renewable, because it is the governing physical process in all heat-exchanger designs, such as pillow-plate [6], fin-and-tube [7], cross-flow [8], double-pipe [9], flat-plate [10], tube-bundle [11], and shell-and-tube [12]. It also plays an important role when investigating porous media [13,14].

Another important industry in which CHT is of great interest is transportation, whether to improve the cooling of blades in aero engines [15,16], or in the automotive industry either for conventional internal combustion engines [17,18] or, for example, cooling of battery packs for electrical vehicles [19].

Fast and accurate prediction of CHT is highly desired in industry. Building and testing a prototype is costly and takes much time. Measuring temperatures, especially in proximity to or on surfaces, is challenging. In addition, knowledge of the flow field allows more advanced optimizations than would be possible with prototypes. An established alternative is, therefore, computational fluid dynamics (CFD), which can save enormous costs and significantly reduce the time to market.

The geometries of the applications, as mentioned earlier, are often very complex, and accurate prediction of heat transfer demands a high-quality mesh. Generating such a mesh for standard CFD methods such as the finite volume method (FVM) usually requires weeks of work by highly skilled personnel. For porous media, creating a suitable mesh may even be impossible. Here, the lattice-Boltzmann method (LBM) offers a clear advantage, as meshing can be highly automated and usually takes only a few seconds.

Turbulent flow significantly affects heat transfer by enhancing it tremendously. Therefore, the fast and accurate simulation of turbulent flow is crucial. The additional computations required for a large eddy simulation (LES) subgrid model have a relatively small impact on simulation time because LBM is typically limited by memory bandwidth.

LBM has made significant progress in accuracy and stability in recent years, and all needed types of boundary conditions are readily available. In addition, the algorithm can be massively parallelized relatively effortlessly due to its local formulation. With advanced compilers and parallel algorithms from the C++ standard library, it is even possible to develop it for cross-platform architectures, as recently demonstrated by Latt et al. [20]. All of this makes it a compelling alternative to more established methods, such as FVM, for studying CHT problems in industry.

However, the problem arises that the standard LBM can only be used for athermal, incompressible flow problems. The most popular approach to incorporating the energy equation when using the LBM is to include a second set of equations, which recover an advection-diffusion equation (ADE) for temperature. Unfortunately, this type of energy equation does not automatically satisfy the conjugate conditions (temperature continuity and continuity of normal heat flux) at walls and interfaces.

In recent years, many researchers have developed methods for the LBM to solve this problem. Korba et al. [21] investigated the accuracy of CHT schemes for the LBM and divided the methods into several groups, for example:

(a)

Interface-considering schemes: For example, Li et al. [22] developed such a method, where the basic equations of LBM remain unchanged, but in the neighboring cells of the interface, the populations are corrected so that the heat flux density across the interface remains constant.

(b)

Additional source term: One example is the study of Karani and Huber [23], where by adding an additional source term to the standard LBM the conjugate conditions at interfaces can be enforced.

(c)

Modified equilibrium distribution: Lu et al. [24] presented a unification of methods in which the enthalpy is introduced into the ADE as the transported quantity instead of the temperature, together with incorporating an additionally needed source term into the equilibrium distribution of the resting population.

A more recent CHT model, which is, therefore, not included in [21], is the work of Kiani-Oshtorjani et al. [25], in which both a source term and a modified equilibrium distribution were used. In this work, however, results are shown only for conductive heat transfer, although the model should also be able to handle moving fluids. Although Korba et al. [21] presented a synopsis of the currently available methods together with an investigation of the accuracy of these methods, it is still not clear which range of thermophysical properties relevant for typical industrial applications can be simulated stably and accurately together with realistic boundary conditions.

The application of the LBM has recently gained popularity in studying conjugate heat transfer on a pore scale. For example, Paknahad et al. [26] studied metal foam at pore scale using a regularized collision operator for the thermal solver, but only steady-state conditions were considered. In another study, Korba and Li [27] investigated thermal convection and conjugate heat transfer in porous media and resolved the turbulent flow through direct numerical simulation.

This study aims to provide a deeper understanding of the limitations of currently available CHT methods for LBM and to propose an applicable methodology for thermal flows oriented towards industrial-relevant problems with realistic thermophysical properties. Therefore, it was first reviewed which of the currently available methods for CHT for the LBM are suitable for simple and complex geometries together with realistic thermophysical properties, particularly at high temperatures. To gain further insight into the stability and accuracy of LBM for CHT, a differentially heated cavity filled with blocks was investigated by extensive parameter studies. The dimensions of the geometry were taken from Merrikh and Lage [28]. As a first step, the range of relaxation times for the LBM in which fast and accurate results can be expected was determined based on steady-state results. The individual relaxation times in the LBM are directly linked to the thermophysical properties of the fluids and solids. Furthermore, the influence of the solid-to-fluid ratio of the thermophysical properties on the error when using the CHT method was examined by performing a parameter study on the relaxation times of the energy equation. Finally, the test case was evaluated for selected cases of realistic thermophysical properties and boundary conditions, which have not been investigated before, to the authors’ knowledge.

The analysis is then carried out considering an insulated cavity with an electric heating pad as the heat source regarding the transient behavior of the LBM in the presence of a heat source term and a more complex convective boundary condition on the outer domain boundaries. As the basis for the error estimations, the deviations from reference FVM simulations were used for both test cases.

The remaining paper is divided into four additional sections. Section 2 presents the governing equations and numerical methods, including all boundary conditions used in this work. Section 3 explains the investigated test cases and lists the used thermophysical properties. Section 4 presents and discusses the results. Section 5 concludes the paper and gives a brief outlook on future work.

2. Numerical Methods

The numerical methods section is divided into two subsections. In the first subsection, the governing equations of CHT are introduced. The second subsection presents the LBM, which can be used to solve these equations and also discusses in detail the modifications that must be made to the standard LBM.

2.1. Governing Equations

The conjugate heat transfer (CHT) problems considered in this work can be mathematically described using the following equations: the continuum Equation (1), the momentum Equation (2), and a simplified form of the energy Equation (3):

$$\sum _{j=1}^{N_D}\frac{{\partial }^{}{u}_j}{\partial {x_j}^{}}=0$$

(1)

$$\frac{{\partial }^{}{u}_i}{\partial t^{}}+\sum _{j=1}^{N_D}{u}_i\frac{{\partial }^{}{u}_i}{\partial {x_j}^{}}=-\frac{1}{\rho }\frac{{\partial }^{}p}{\partial {x_i}^{}}+\sum _{j=1}^{N_D}\frac{{\partial }^{}}{\partial {x_j}^{}}\left(\nu \frac{{\partial }^{}{u}_i}{\partial {x_j}^{}}\right)+S_i$$

(2)

$$\rho c_p\frac{\partial T}{\partial t}+\sum _{j=1}^{N_D}\frac{\partial }{\partial x_j}\left[u_j\rho c_{p}T\right]-\sum _{j=1}^{N_D}\frac{\partial }{\partial x_j}\left[k\frac{\partial T}{\partial x_j}\right]=0$$

(3)

In Equations (1)–(3), $\rho$ is the density, $u_i$ and $u_j$ are the velocity vector in index notation, p is the pressure, T is the temperature, $x_j$ is the position vector, t is the time, $N_D$ is the number of dimensions, $\nu$ is the kinematic viscosity, k is the thermal conductivity, and $c_p$ is the specific heat capacity (for constant pressure).

2.2. Lattice-Boltzmann Method

The LBM evolution equation reads

$$f_{\alpha }(\mathit{x}+{\mathit{c}}_{\alpha }\delta t,t+\delta t)=f_{\alpha }(\mathit{x},t)+{\Omega }_{\alpha }+S_{\alpha },$$

(4)

where $f_{\alpha }$ are velocity probability density functions for particle clusters, also called populations. $\mathit{x}$ is the position vector, t is the time, and $\delta t$ is the time increment. ${\mathit{c}}_{\alpha }$ is the discrete velocity of the population $f_{\alpha }$, ${\Omega }_{\alpha }$ is the collision operator, and S is an arbitrary source term that could represent, for example, a body force.

The set of discrete velocities depends on the lattice configuration, typically denoted by $D{N}_{D}Q{N}_{\alpha }$, where $N_D$ is the number of dimensions and $N_{\alpha }$ is the number of discrete velocities. For example, the $D2Q9$ lattice is the most commonly used lattice for two-dimensional fluid simulations, while the $D2Q5$ is predominantly used for advection-diffusion type problems, where the reduced set of discrete velocities is sufficient for recovering an ADE on the macroscopic level.

Figure 1a shows the $D2Q9$ lattice configuration, where eight velocities are defined by directly connecting the current lattice and its neighbors. The ninth velocity is zero and assigned to the rest population.In Figure 1b, it can be seen that the $D2Q5$ lattice can be interpreted as a pruned version of the $D2Q9$.

Density $\rho$ and the momentum vector $\rho \mathit{u}$ can be recovered as zero- and first-order moments of the velocity probability density functions, respectively, as

$$\begin{array}{cc}\hfill \rho (\mathit{x},t)& =\sum _{\alpha =1}^{N_{\alpha }}{f}_{\alpha }(\mathit{x},t)\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \rho \mathit{u}(\mathit{x},t)& =\sum _{\alpha =1}^{N_{\alpha }}{\mathit{c}}_{\alpha }{f}_{\alpha }(\mathit{x},t)\hfill \end{array}$$

(6)

From Equations (5) and (6), the velocity vector $\mathit{u}$ can be readily calculated. Several collision operators are available. The most commonly used are the single-relaxation-time (SRT), often referred to as the Bhatnagar–Gross–Krook (BGK) [29] model, and the multiple-relaxation-time (MRT) [30] collision operators. The SRT is defined by a single relaxation time, while the MRT schemes possess as many relaxation times as lattice velocities. A special case of MRT is the two-relaxation-time (TRT) [31] scheme. In addition to the standard MRT, there are more sophisticated variants, such as the cascaded [32] and cumulant [33] collision operators, which will not be considered in this study. The simplest variant is the BGK [29], given by

$${\Omega }_{\alpha }=\frac{\delta t}{{\tau }_f}\left(f_{\alpha }^{\mathrm{eq}}-f_{\alpha }\right),$$

(7)

where ${\tau }_f$ is the relaxation time, linked to the viscosity of the fluid by

$${\tau }_f=\frac{\nu }{c_s^2}+\frac{\delta t}{2}.$$

(8)

and $c_s$ is the lattice speed of sound, which depends on the lattice configuration (e.g., D2Q9: $c_s=\sqrt{1/3}$). It is important to mention that a LBM simulation is typically performed on a uniform grid (i.e., $\delta x=\delta y=\delta z=\mathrm{const}.$) and in lattice units defined by a set of conversion factors resulting from a conveniently rescaled system, where the time increment $\delta t$, the spatial step size $\delta x$ and the reference density of the fluid ${\rho }_{f,0}$ in lattice units are unity (i.e., $\delta t=1$, $\delta x=1$, ${\rho }_{f,0}=1$). The system must be rescaled in such a way that the non-dimensional numbers determining the flow behavior are equal to the original unscaled system.

The TRT [31] uses two different relaxation times. One for the symmetric (even) part of the populations $f_{\alpha }^s=\left(f_{\alpha }+f_{\overline{\alpha }}\right)/2$ and one for the anti-symmetric (odd) part $f_{\alpha }^a=\left(f_{\alpha }-f_{\overline{\alpha }}\right)/2$, where $f_{\overline{\alpha }}$ is the population associated with the discrete velocity ${\mathit{c}}_{\overline{\alpha }}$ in opposite direction to ${\mathit{c}}_{\alpha }$ and the corresponding symmetric and anti-symmetric equilibria are calculated similarly:

$${\Omega }_{\alpha }=\frac{\delta t}{{\tau }_f^s}\left(f_{\alpha }^{\mathrm{eq},\mathrm{s}}-f_{\alpha }^s\right)+\frac{\delta t}{{\tau }_f^a}\left(f_{\alpha }^{\mathrm{eq},\mathrm{a}}-f_{\alpha }^a\right)$$

(9)

The symmetric relaxation time ${\tau }_f^s$ is related to the fluid viscosity $\nu$ (similar to Equation (8)), while the anti-symmetric relaxation time ${\tau }_f^a$, given by

$${\tau }_f^a=\frac{\Lambda +0.5({\tau }_f^s-0.5)}{{\tau }_f^s-0.5}$$

(10)

can be used to fine tune the behavior of the collision operator by choosing an appropriate value for $\Lambda$ [34]. Since $\Lambda$ has no direct physical meaning, it is often called magic operator. In the present work, where not otherwise stated, the TRT was optimized for stability by setting $\Lambda =1/4$ [34].

Theoretically, it would be possible to relax every population $f_{\alpha }$ towards its equilibrium $f_{\alpha }^{\mathrm{eq}}$ with its own individual relaxation time ${\tau }_f$. By transforming the populations into momentum space and performing the collision there, it is possible to relax the moments directly rather than populations, simplifying the task of finding individual relaxation times to some extent. However, based on the choice of moments, there are still free parameters that must be fine-tuned for stability and accuracy. For example, the considered moments $m_{\alpha }$, following from the Gram–Schmidt procedure, for the D2Q9 lattice configuration are density $\rho$, kinetic energy e, energy squared $ϵ$, x-momentum $j_x$, x-energy-flux $q_x$, y-momentum $j_y$, y-energy-flux $q_y$, normal stress $p_{xx}$, and shear stress $p_{xy}$.

The MRT [30] collision operator in vector notation is

$$\begin{array}{cc}\hfill \Omega & =\delta t\left({\mathit{M}}^{-1}\mathit{S}\mathit{M}\left({\mathit{f}}^{\mathrm{eq}}-\mathit{f}\right)\right)\hfill \\ \hfill & =\delta t\left({\mathit{M}}^{-1}\mathit{S}\left({\mathit{m}}^{\mathrm{eq}}-\mathit{m}\right)\right),\hfill \end{array}$$

(11)

where $\mathit{M}$ is a transformation matrix from population into momentum space, and ${\mathit{M}}^{-1}$ is the inverse Matrix accomplishing the reverse operation. $\mathit{S}$ is a matrix containing the relaxation times for the individual moments $m_{\alpha }$. Another advantage of the MRT is that the equilibrium moments $m_{\alpha }^{\mathrm{eq}}$ can be computed analytically and do not need to be approximated like their counterparts in population space $f_{\alpha }^{\mathrm{eq}}$.

Using the Chapman–Enskog theory [35], it can be shown that the LBM, using the most common lattice configurations, recovers the incompressible athermal Navier–Stokes equations ((Equations (1) and (2)) as introduced in Section 2.1.

The most commonly used method for incorporating the energy equation is to introduce a second population $g_{\alpha }$ to simulate an ADE for the temperature as a passive scalar. It should be noted that conservation of energy is not necessarily guaranteed in this way since viscous dissipation, pressure work, and kinetic energy are neglected. However, this approximation is often acceptable.

The LBM evolution equation for the ADE for the temperature T reads formally similar to its Navier–Stokes counterpart (Equation (4)):

$$g_{\alpha }(\mathit{x}+{\mathit{c}}_{\alpha }^g\delta t,t+\delta t)=g_{\alpha }(\mathit{x},t)+{\Omega }_{\alpha }^g+S_{\alpha }^g$$

(12)

The main difference is in the definition of the equilibrium distributions $g_{\alpha }^{\mathrm{eq}}$ and that, in general, fewer discrete velocities are needed to recover an ADE compared to the incompressible, athermal Navier–Stokes equations. Again, using Chapman–Enskog expansion, it could be shown that the LBM for the ADE recovers an energy equation in the form of

$$\frac{\partial T}{\partial t}+\sum _{j=1}^{N_D}\frac{\partial }{\partial x_j}\left[u_{j}T\right]-\sum _{j=1}^{N_D}\frac{\partial }{\partial x_j}\left[a\frac{\partial T}{\partial x_j}\right]=0,$$

(13)

where $u_j$ is the velocity vector $\mathit{u}$ in index notation and a is the thermal diffusivity.

Even though the temperature T is connected to the flow-field by the advective term in Equation (13), it can only be seen as a one-way coupling since the density in LBM is not a function of the temperature. For relatively small temperature changes, a popular method to restore a two-way coupling between the Navier–Stokes equations and the energy equation is the Boussinesq approximation, which was also used in the present work.

There are several methods to include the buoyancy source term resulting from the Boussinesq approximation into the LBM. In the present work, the exact difference method (EDM) of Kupershtokh [36] was used. In accordance with the Boussinesq approximation, a force density vector field is first calculated from

$$\mathit{F}(\mathit{x},t)={\mathit{g}}^a{\rho }_0\beta \left(T-T_0\right),$$

(14)

assuming this force field is constant during a time step and the time step is small enough in order for the forces to be considered pulse-like, the change in momentum can be approximated as

$$\Delta \left(\rho \mathit{u}\right)=\mathit{F}(\mathit{x},t)\Delta t.$$

(15)

If, in addition, it is assumed that the density remains constant, the change in velocity can be calculated from

$$\Delta \mathit{u}={\mathit{u}}_s-\mathit{u}=\mathit{F}(\mathit{x},t)\frac{\Delta t}{\rho }.$$

(16)

Hence, the shifted macroscopic velocities are given by

$${\mathit{u}}_s=\mathit{u}+\mathit{F}\frac{\Delta t}{\rho }$$

(17)

and are then used to calculate a shifted equilibrium distribution $f_s^{\mathrm{eq}}=f^{\mathrm{eq}}\left({\mathit{u}}_s\right)$. Finally, the difference in $f^{\mathrm{eq}}$ and $f_s^{\mathrm{eq}}$ is added as a source term to the LBM during the collision step as

$$S=f^{\mathrm{eq}}-f_s^{\mathrm{eq}}.$$

(18)

The presented LBM would suffice for all further investigations if only fluid flow were of interest. The problem, however, is that in the context of CHT, in order to account for a spatially varying (but temporally constant) density $\rho$, specific isobaric heat capacity $c_p$, and thermal conductivity k within the domain, at least an energy equation in the form of Equation (3) instead of Equation (13) needs to be considered.

Following the review paper of Korba and Li [21], the methods to solve this problem can be divided into several groups, for example:

(a)

Interface-considering schemes: Li et al. [22] developed such a method, where the basic equations of LBM remain unchanged, but in the neighboring cells of the interface, the populations $g_{\alpha }$ are corrected so that the heat flux density $\dot{q}$ across the interface remains constant. Although this method would have advantages for cases where the wall is not located exactly between two lattices, it was not used in the present work because it is not clear how to account for the discontinuity of the wall-normal vectors at sharp corners.

(b)

Additional source term: By comparing Equation (3) with Equation (13), a source term $S_c$ can be identified, which, when added to the LBM using an appropriate method, yields Equation (3) on a macroscopic level. Karani and Huber [23] developed such a method. However, as also pointed out by Korba and Li [21], it is problematic that this method uses a finite difference approximation of the term $1/\left(\rho c_p\right)$. The derivative of a jump is not well-defined mathematically. Therefore, this method was not used in the present work.

(c)

Modified equilibrium distribution: This approach was proposed by Lu et al. [24], where instead of the temperature T, the enthalpy $H=\rho c_{p}T$ is introduced into the ADE as the transported quantity. The additionally required source term arising from this change in variables is included in the modified equilibrium distribution of the rest population $f_{N_{\alpha }}^{\mathrm{eq}}$. This method was chosen in the present study and will be briefly introduced in the next section.

2.2.1. Conjugate Heat Transfer Method

The basic LBM evolution equation remains unchanged and is recalled from the previous section as

$$g_{\alpha }(\mathit{x}+{\mathit{c}}_{\alpha }^g\delta t,t+\delta t)=g_{\alpha }(\mathit{x},t)+{\Omega }_{\alpha }^g+S_{\alpha }^g.$$

(19)

Instead of the temperature T, the enthalpy H is recovered as zero-order moment $g_{\alpha }$ of the ADE

$$H=\rho c_{p}T=\sum _{\alpha =1}^{N_{\alpha }}{g}_{\alpha }$$

(20)

and the modified equilibrium distribution reads in population space

$$\begin{array}{cc}\hfill g_{\alpha =N_{\alpha }}^{\mathrm{eq}}& =T\left[\rho c_p+{\left(\rho c_p\right)}_0(w_{\alpha }^g-1)\right]\hfill \\ \hfill g_{\alpha }^{\mathrm{eq}}& =w_{\alpha }^{g}T\left[{\left(\rho c_p\right)}_0+\rho c_p{\mathit{c}}_{\alpha }^g·\mathit{u}\right],\hfill \end{array}$$

where $w_{\alpha }^g$ are the lattice weights for the corresponding lattice velocities $c_{\alpha }^g$ and ${\left(\rho c_p\right)}_0$ is a reference-specific volumetric heat capacity at constant pressure, usually defined as the minimal value of the problem. In the present work, the fluid term is always smaller than the solid term, which is usually valid if the fluid is a gas. This equilibrium distribution can be used by the BGK and TRT collision operators. The equilibrium vector in momentum space [21] is for the standard ADE in terms of the temperature T, given by

$${\mathit{m}}_T^{\mathrm{eq}}={\left[0,T,uT,vT,2/3T\right]}^T$$

(21)

and for the CHT method, the modified momentum equilibrium vector reads, according to [21],

$${\mathit{m}}_H^{\mathrm{eq}}={\left[0,\rho c_{p}T,\rho c_{p}uT,\rho c_{p}vT,4\rho c_{p}T-10/3{\left(\rho c_p\right)}_{0}T\right]}^T.$$

(22)

2.2.2. Boundary Conditions

Figure 2 schematically depicts the lattice configuration near the wall as used in the present work, where ${\mathit{x}}_s$ is the location of a solid lattice, ${\mathit{x}}_b$ is the location of the first fluid lattice, and ${\mathit{x}}_w$ is the location of the wall, which is assumed to be approximately between the two adjacent lattices.

For the fluid populations at the resting walls, the half-way bounce-back [37] was used:

$$f_{\alpha }({\mathit{x}}_b,t+\delta t)=f_{\overline{\alpha }}^{*}({\mathit{x}}_b,t),$$

(23)

where the wall-normal population $f_{\alpha }$, at the cell location next to the wall ${\mathit{x}}_b$, at the new time level $t+\delta t$, is set to the post-collision population of the opposite direction $f_{\overline{\alpha }}^{*}$ from the current time step.

Three types of boundary conditions were required for the thermal populations. A Dirichlet-type boundary condition was implemented using the anti-bounce-back method of Ginzburg et al. [38]. For adiabatic boundary conditions, the zero-gradient method of Malaspinas [39] was used, which estimates a wall temperature based upon a Taylor series expansion in proximity to the wall. This temperature was then enforced by using the anti-bounce-back method of Ginzburg et al. [38]. Additionally, the convective heat flux boundary condition of Huang et al. [40] was used, which allows setting a specific value for the heat-transfer coefficient h.

The anti-bounce-back boundary condition [38] reads

$$g_{\alpha }({\mathit{x}}_b,t+\delta t)=-g_{\overline{\alpha }}^{*}({\mathit{x}}_b,t)+2w_{\alpha }{g}_{\alpha }^{\mathrm{eq}}({\mathit{x}}_w,t+\delta t),$$

(24)

where $g_{-\alpha }^{*}({\mathit{x}}_b,t)$ is the population into the opposite direction of $g_{\alpha }({\mathit{x}}_b,t+\delta t)$ of the old time step at the boundary node position ${\mathit{x}}_b$, $w_{\alpha }$ is the corresponding lattice weight and $g_{\alpha }^{\mathrm{eq}}$ is the equilibrium distribution, which depends on the wall temperature $T_w$ and the wall velocity vector ${\mathbf{v}}_w$.

The adiabatic boundary condition follows the approach proposed by Malaspinas [39] but differs in the lattice configuration near the wall. While [39] assumes the lattice-center coincident with the wall interface, in this study, the interface is located approximately between two adjacent lattices. Therefore, the evaluation points of the Taylor series expansion are different.

Evaluating the Taylor series expansion one and two lattices away from the wall at indices $i+1$ and $i+2$ (see Figure 2), neglecting third- and higher order terms, and setting the temperature gradient at the wall to zero $\left(\partial T/\partial x\right){|}_{x_w}=0$, the wall temperature readily follows as

$$T_w=\frac{25T_{i+1}-9T_{i+2}}{16}.$$

(25)

The computed wall temperature $T_w$ is then plugged into the anti-bounce-back formulation (Equation (24)) to fulfill the adiabatic boundary requirement automatically.

The convective heat flux boundary condition was implemented as proposed by Huang et al. [40] and reads

$$g_{\alpha }=\frac{g_{\alpha }^n+g_{\overline{\alpha }}^n+\left(w_{\alpha }+w_{\overline{\alpha }}\right)\frac{h}{k}\delta x{T}_{\mathrm{amb}}}{1+\frac{h}{k}\delta x}-g_{\overline{\alpha }},$$

(26)

where superscript n indicates adjacent populations at the interface (e.g., $g_{\alpha }^n$), h is the heat transfer coefficient, k is the heat conductivity, and $T_{\mathrm{amb}}$ is the ambient temperature on the other side of the interface.

3. Test Cases

In the present work, two test cases were investigated. First, a differentially heated cavity filled with blocks was used to gain insight into the relationship between relaxation times and errors. The influence of the CHT method on the error was studied in detail for a wide range of fictive combinations of thermophysical properties. During the simulations, the thermophysical properties are not modeled as temperature-dependent. To the best of the authors’ knowledge, this has not been demonstrated for LBM and is an intrinsic limitation of the standard method. Finally, three realistic combinations of thermophysical properties were tested. For estimating the error of a LBM simulation, the relative root mean squared deviation (RMSD) of the temperature gradient at the hot wall (see Figure 3) from a corresponding FVM simulation was calculated as:

$$\mathrm{rel}.\mathit{RMSD}=\sqrt{\frac{{\sum }_{n=1}^{N_y}{\left(\left({\nabla }_x{T}_w^{LBM}-{\nabla }_x{T}_w^{FVM}\right)/{\nabla }_x{T}_w^{FVM}\right)}^2}{N_y}},$$

(27)

where $N_y$ is the number of cells (or lattices) in the y-coordinate direction along the wall, and ${\nabla }_x{T}_w$ is the wall-normal temperature gradient at the hot wall approximated by the forward difference quotient

$${\nabla }_x{T}_w=\frac{{\partial }^{}{T}_w}{\partial x^{}}\approx \frac{T_f-T_w}{x_f-x_w},$$

(28)

where $T_w=T_h$ is the hot-wall temperature, $x_w$ is the position of the wall, and $x_f$ is the location of the first fluid cell adjacent to the hot wall.

The second test case, an insulated cavity with an electric heating pad as a heat source, was investigated to scrutinize the transient behavior of the LBM compared to a FVM simulation. Compared to the first test case, a more complex convective boundary condition was used at the outer domain boundaries.

All FVM reference simulations were performed using the commercial software package ANSYS Fluent 19.2 [41]. To allow direct comparison, the FVM simulations were set up as similar as possible to the LBM simulations, using a uniform Cartesian mesh and (temporally) constant thermophysical properties. As in the LBM simulations, a buoyancy term based on the Boussinesq approximation was included. Pressure-velocity coupling was achieved by using a coupled algorithm, and the pressure interpolation was carried out by the body-force-weighted scheme. Both the momentum and energy equations were discretized using the second-order upwind method. A first-order implicit method was chosen for the transient formulation. Implementation details can be found in the Ansys Fluent Theory Guide [42]. Mesh refinement studies were conducted to ensure that the mesh resolution did not influence the reference solutions.

3.1. Cavity Filled with Solid Blocks

The geometry is defined by Merrikh and Lage [28], who used FVM to evaluate various configurations. Though the LBM was already compared to FVM in studies such as Karani and Huber [23], and Lu et al. [24], a thorough investigation of the limitations of the CHT methods and their accuracy for realistic physical properties was still lacking.

In the present work, $N=16$ solid blocks were used, and the solid to empty cavity volume ratio is $\gamma =0.36$. This results in the following dimensions: $L=1.0$, $B=0.15$, and $D=0.1$. Figure 3 shows the schematics of the simulation domain.

The test case was studied under the influence of a gravitational field $\mathit{g}$. The left wall is fixed at a constant high temperature $T_h$, while a constant cold Temperature $T_c$ is considered at the right wall. The top and bottom walls are assumed to be adiabatic. The no-slip boundary condition is selected for all walls.

The test case was investigated for a Rayleigh number of

$$\mathit{Ra}=\mathit{Gr}\mathit{Pr}=1_{\mathrm{E}}-5.$$

(29)

$\mathit{Gr}$ is the Grashof number, defined as

$$\mathit{Gr}=\frac{\left|\mathit{g}\right|\beta \Delta T{L}^3}{{\nu }^2},$$

(30)

$\nu$ is the dynamic viscosity, $\beta$ is the (isobaric) thermal expansion coefficient, $\left|\mathit{g}\right|$ is the magnitude of the gravitational field $\mathit{g}$, $\Delta T=T_h-T_c$ is the temperature difference between hot and cold wall, and L is the length and height of the square cavity. $\mathit{Pr}$ is Prandtls number, given by

$$\mathit{Pr}=\frac{\nu }{a}=\frac{\mu c_p}{\lambda },$$

(31)

where a is the thermal diffusivity, $\lambda$ is the thermal conductivity, $\mu$ is the dynamic viscosity, and $c_p$ is the specific isobaric heat capacity.

For the FVM simulations, a fictitious fluid was defined, corresponding to the Rayleigh number Ra = $1\times {10}^{-5}$. The temperature difference ($\Delta T=T_h-T_c$) should not be too small to avoid round-off errors. Finally, the product of the isobaric thermal expansion coefficient and temperature difference must be small (i.e., $\beta \Delta T\ll 1$). The magnitude of the gravitational acceleration was conveniently chosen as $\left|\mathit{g}\right|=10$ to yield pleasing numerical values for the thermophysical properties. These considerations led to the thermophysical properties and boundary conditions listed in

Table 1. Fictitious thermophysical properties and boundary conditions for the first test case.

$\mathit{\rho }$	1	$\mathrm{k}\mathrm{g}$/${\mathrm{m}}^3$	$\mathit{\mu }$	$0.001$	$\mathrm{k}\mathrm{g}$/$\mathrm{m} \mathrm{s}$	L	1	$\mathrm{m}$
$c_p$	100	$\mathrm{J}$/$\mathrm{k}\mathrm{g}$ $\mathrm{K}$	$\beta$	$0.0005$	1/$\mathrm{K}$	$T_c$	$563.15$	$\mathrm{K}$
$\lambda$	$0.1$	$\mathrm{W}$/$\mathrm{m}$ $\mathrm{K}$	$\left|\mathit{g}\right|$	10	$\mathrm{m}$/${\mathrm{s}}^2$	$T_h$	$583.15$	$\mathrm{K}$

. The LBM simulations were performed in lattice units. Thus, the fictitious thermophysical properties were only needed to convert the system into a physical dimensional system to allow direct comparison with results from the FVM.

For the parametric study, the thermophysical properties of the solids ($k_s$ and ${\left(\rho c\right)}_s$) were calculated from the solid/fluid ratios of the thermal conductivities $r_k=k_s/k_f$, and the product of density and heat capacity $r_c={\left(\rho c\right)}_f/{\left(\rho c\right)}_s$. In the case where the product of density and heat capacity is smaller for the fluid than for the solid, which is typically valid if the fluid is a gas, i.e., ${\left(\rho c\right)}_f<{\left(\rho c\right)}_s\to {\left(\rho c\right)}_0={\left(\rho c\right)}_f$, the relationship between the fluid relaxation time ${\tau }_f$ and the ratio of the thermal conductivities $r_k$ is as follows:

$$\begin{array}{cc}\hfill {\tau }_{gf}& =\frac{1}{\mathit{Pr}}{\tau }_f+\frac{\mathit{Pr}-1}{2\mathit{Pr}}\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill {\tau }_{gs}& =\frac{r_k}{\mathit{Pr}}{\tau }_f+\frac{\mathit{Pr}-r_k}{2\mathit{Pr}}\hfill \end{array}$$

(33)

Note that $r_c$ does not appear in Equations (32) and (33).

Three combinations of realistic thermophysical properties representative of high-temperature applications were selected for fluids and solids (see

Table 2. Physical properties for the verification test cases.

	Fluid	$\mathit{\rho }$	${\mathit{c}}_{\mathit{p}}$	$\mathit{\lambda }$	$\mathit{\mu }$	$\mathit{\beta }$
Case	Solid	$\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$	$\mathrm{J}/\mathrm{k}\mathrm{g} \mathrm{K}$	$\mathrm{W}/\mathrm{m} \mathrm{K}$	$\mathrm{k}\mathrm{g}/\mathrm{m} \mathrm{s}$	$1/\mathrm{K}$
A	Air	0.5663	1056.6	0.04737	3.15787 × 10${}^{-5}$	0.001605
	Wool	145	840	0.07699
B	Air	1.0928	1007.44	0.02808	1.96342 × 10${}^{-5}$	0.003102
	Steel	7800	458.9	49.1
C	Flue gas	0.2842	1360	0.0454	8.2204 × 10${}^{-5}$	0.00125
	Steel	7800	599.3	27.28

). The properties for air were obtained at the respective mean temperature levels from the open-source thermophysical property library CoolProp [43,44]. In contrast, the properties for flue gas, used in case C, are the mass-averaged values from a FVM simulation of a walking hearth type reheating furnace carried out by Prieler et al. [45]. The properties of wool used in case A are from internal communications with a mineral wool producer. The steel properties in cases B and C are for low-alloyed steel and were obtained from the graph given in [45].

All simulations were started with a resting fluid, and the temperature field was initialized with the mean temperature $T_0=(T_c+T_h)/2$. For the parametric study on relaxation times, a mesh study showed that a grid resolution of 280 by 280 lattices is sufficient for the LBM simulations. When conducting a mesh independence study, it is important to choose a representative and sensitive quantity. In this work, the average Nusselt number along the hot wall, as defined by Merrikh and Lage [28], was selected:

$$N{u}_{avg}=\frac{h_{avg}L}{{\lambda }_f}=\frac{L}{T_h-T_c}\left[-{\int }_0^L\frac{{\partial }^{}T}{\partial x^{}}{|}_{x=0}dy\right]$$

(34)

Figure 4 shows the results of the mesh independence study for the TRT collision operator. It can be seen that for all relaxation times $\tau$, the grid resolution of 280 by 280 lattices is within the asymptotic range. Furthermore, it could be argued that even the resolution of 240 by 240 lattices is within 1% of the average Nusselt number obtained with a grid resolution of 560 by 560 lattices, with the exception of the largest relaxation time of $\tau =1.5$, where the deviation would still be less than 2%, however. For the investigation of realistic thermophysical properties, higher mesh resolutions up to 2240 by 2240 lattices were needed for stability reasons. The actual mesh resolutions used per case are given in the discussion of the results.

3.2. Insulated Cavity with a Heating Pad

As a second test case for a realistic combination of physical properties, an insulated cavity filled with air and heated by a heating pad was selected. Figure 5 shows the geometry and dimensions in mm, where the gray area consists of mineral wool that serves as insulation. The cavity in the center is filled with air, and the red section represents the heating pad. The thermophysical properties of air and wool are listed in

Table 3. Physical properties of air and wool for the insulated cavity with heating-pad test case.

	Air	Wool
$\rho$	0.4551	145	$\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$
$c_p$	1090	840	$\mathrm{J}/\mathrm{k}\mathrm{g} \mathrm{K}$
$\lambda$	0.05536	0.11343	$\mathrm{W}/\mathrm{m} \mathrm{K}$
$\mu$	3.6281 × 10${}^{-5}$		$\mathrm{k}\mathrm{g}/\mathrm{m} \mathrm{s}$
$\beta$	0.0013075		$1/\mathrm{K}$

.

The simulation was initialized at time $t=0\mathrm{s}$ with a constant temperature of $T=23{}^{\circ }$C and the fluid at rest. The heating pad provides a constant electrical heating power of $P_{\mathrm{el}}=205\mathrm{W}$. Assuming natural convection outside of the domain, a convective boundary condition with an ambient temperature of $T_{\mathrm{amb}}=23{}^{\circ }$C was used together with a heat transfer coefficient of $h=4\mathrm{W}/{\mathrm{m}}^2$ for the outer domain boundary. The case setup and boundary conditions are given in

Table 4. Mesh resolution, boundary conditions, and relaxation times for the insulated cavity with heating pad test case.

$N_x$	335	−	${\mathrm{P}}_{\mathrm{el}}$	205	$\mathrm{W}$	${\mathit{\tau }}_{\mathit{f}}$	0.55	−
$N_y$	330	−	$\left|\mathit{g}\right|$	9.81	$\mathrm{m}/{\mathrm{s}}^2$	${\tau }_{gf}$	0.5711	−
$\Delta x$	1	$\mathrm{m}\mathrm{m}$	$T_{\mathrm{amb}}$	23	${}^{\circ }\mathrm{C}$	${\tau }_{gs}$	0.5615	−

.

4. Results and Discussion

4.1. Cavity Filled with Solid Blocks

The investigations for this test case are divided into four parts. First, correct implementation of the CHT method is demonstrated. Second, a parametric study for the relaxation times was performed, with all relaxation times set equal (${\tau }_f={\tau }_g$). This served to determine a reasonable range of relaxation times to be used in the third part of the study, where the influence of the CHT method was examined in more detail. In addition, the magnitude of the deviation without the influence of the CHT method could be determined. For each LBM simulation, a corresponding FVM simulation was performed and used as the reference for estimating the error by calculating the relative root mean squared deviation (RMSD) of the wall-normal temperature gradient at the hot wall. In the fourth part of this test case, realistic combinations of physical properties were investigated.

Figure 6 shows, as an example, the predictions for the temperature field (Figure 6a) together with predictions for the velocity field and streamlines (Figure 6b) obtained with the LBM using a relaxation time of ${\tau }_f={\tau }_g=0.9$ and the TRT collision operator. Since ${\tau }_f={\tau }_g$ implies that the thermal diffusivities of fluid and solid are the same, the influence of the advective heat transport can easily be seen in the temperature field. Without the presence of a gravitational field, the temperature distribution between the hot and cold walls would be linear and quasi-one-dimensional. Figure 6b shows that the fluid predominantly flows at high velocity clockwise between the outer layer of solid blocks and the four middle blocks. Note that the dimensional values were obtained from the simulation in lattice units using the fictitious medium defined by the thermophysical properties and boundary conditions given in Table 1.

4.1.1. Implementation Validation

To ensure correct implementation of the CHT model [24], qualitative and quantitative comparison to previously published work by Merrikh and Lage [28], Raji et al. [46], Karani and Huber [23], Lu et al. [47], and Lu et al. [24] is presented. Figure 7 shows the steady-state predictions for the temperature field (Figure 7a) and the streamlines of the velocity field (Figure 7b) as obtained by Lu et al. [24]. In the same figure, the temperature field (Figure 7c) is shown together with the velocity field and streamlines (Figure 7d) obtained using the solver implemented in this work. Both the temperature field and the streamlines for the velocity field agree very well in proximity to the hot and cold walls. A slight deviation can be observed in the central part of the temperature field.

To gain further confidence in the implementation, the results for the average Nusselt number (Equation (34)) were compared with those presented in several other publications, considering variations in the ratio of heat conductivities. As the ratio of the thermal diffusivities was chosen to be unity, the ratio of the heat capacities was scaled accordingly. The Prandtl number for these selected test cases was set to unity ($Pr=1$) and the Rayleigh number was $Ra=1\times {10}^5$. In addition, the results obtained in the present work were also compared to reference FVM simulations. As can be seen in

Table 5. Predictions for $N{u}_{avg}$ along the hot wall (Equation (34)) for various ratios of fluid and solid heat conductivities obtained with the TRT collision operator together with FVM reference simulations, comparison to previously published results from other authors [23,24,28,46,47].

${\mathit{\lambda }}_{\mathit{s}}/{\mathit{\lambda }}_{\mathit{f}}$	[28]	[46]	[23]	[47]	[24]	FVM	LBM
0.1	0.813	0.785	0.7969	0.8143	0.8215	0.8119	0.8178
1	1.233	1.193	1.185	1.2491	1.2425	1.2318	1.2367
10	2.030	2.066	2.031	2.0566	2.0518	2.0262	2.0347
100	2.313	2.394	2.4506	2.3351	2.3355	2.3071	2.3106

, the obtained predictions agree very well with results previously presented by other authors. Additionally, they are in very good accordance with the reference FVM simulations. Furthermore, the confidence in adequate implementation is strengthened by not only the averaged Nusselt number along the hot wall agreeing very well with the carried out FVM simulations, but also the distribution of the local Nusselt number along the hot wall, as can be seen in Figure 8. They agree closely and are essentially indistinguishable; only for the most extreme case of ${\lambda }_s/{\lambda }_f=100$ can a slight deviation between the LBM and FVM predictions be observed.

Overall, the validation studies provide a high level of confidence in the adequacy and accuracy of the implemented numerical method.

In addition, simulations were performed for a different number of solid blocks and compared with the results of Merrikh and Lage [28]. For these simulations, a mesh of the same quality as in the mesh study was created by removing or adding lattices in such a way that the number of lattices per channel and per solid block remained constant. This results in mesh resolutions of 210 by 210 for the 9-block test case, 280 by 280 for the 16-block base case, and 420 by 420 cells for the 36-block test case. Predictions for $N{u}_{avg}$ for various solid-to-fluid heat conductivity ratios and different number of blocks are listed in

Table 6. Predictions for $N{u}_{avg}$ for various solid-to-fluid heat conductivity ratios and different number of blocks. On the left-hand side of each pipe (|) are the results obtained by Merrikh and Lage [28] reported and on the right-hand side, results obtained with the TRT collision operator of the present work (PW).

Blocks	Source	${\mathit{\lambda }}_{\mathit{s}}/{\mathit{\lambda }}_{\mathit{f}}=0.1$	${\mathit{\lambda }}_{\mathit{s}}/{\mathit{\lambda }}_{\mathit{f}}=1$	${\mathit{\lambda }}_{\mathit{s}}/{\mathit{\lambda }}_{\mathit{f}}=10$	${\mathit{\lambda }}_{\mathit{s}}/{\mathit{\lambda }}_{\mathit{f}}=100$
9	[28] | PW	1.053 | 1.048	1.383 | 1.376	2.140 | 2.134	2.418 | 2.405
16	[28] | PW	0.813 | 0.818	1.233 | 1.237	2.030 | 2.035	2.313 | 2.311
36	[28] | PW	0.676 | 0.676	1.098 | 1.098	1.922 | 1.923	2.211 | 2.221

. On the left-hand side of each pipe (|) are the results obtained by Merrikh and Lage [28], and on the right-hand side are the results obtained using the TRT collision operator of the present work (PW). The excellent agreement in each case demonstrates that the current implementation is reliable and accurate for various block configurations.

4.1.2. Base Error

Figure 9 shows the relative root mean squared deviation (RMSD) in % from the reference FVM for the different collision operators, where the TRT collision operator was simulated twice. Once with a magic number of $\Lambda =1/4=TRT-S$ and once with $\Lambda =3/16=TRT-A$. S indicates a fine-tuned anti-symmetric relaxation time ${\tau }_f^a$ for stability and A for accuracy [34], respectively.

The relative simulation time is defined as ${\Delta t/\Delta t|}_{\tau =1.5}$ because the computational effort scales linearly with the time-step size, which depends on the relaxation time given by the diffusive scaling as

$$\Delta t_p=\Delta x_p^2\frac{{\nu }_{LBM}}{{\nu }_P},$$

(35)

where the mesh resolution determines $\Delta x_p$. Assuming the mesh resolution is kept constant, the time step only depends on the viscosity in lattice units ${\nu }_{LBM}$ which is given by ${\nu }_{LBM}=({\tau }_f-0.5)/3$.

From Figure 9, a viable relaxation time in the range of $0.6\le {\tau }_f\le 1.5$ was derived. The upper bound is determined by a deviation of approximately $5\%$ from the reference simulation, and the lower bound follows from simulation time considerations where choosing $\tau =0.6$ instead of $\tau =1.5$ leads to an order of magnitude longer simulation time.

As the limit $\tau =0.5$ is approached, the computational cost increases exponentially, and the simulation becomes less stable. The BGK operator fails, for this test case, at a lower limit of $0.56$. Interestingly, the BGK operator achieves a smaller deviation from the FVM predictions for values of $\tau <1$; however, the deviation increases rapidly for relaxation times above $\tau =1$. It should also be mentioned that it is inherently less stable than TRT and MRT. Second, the influence of the magic operator $\Lambda$ on the overall prediction is small. The lines for both types of TRT simulations in Figure 9 almost coincide. The MRT yields very similar results for $\tau <1$ but behaves better for larger values of $\tau$. For the present test case, there was no significant advantage in terms of stability to using the MRT. It is, therefore suggested that the TRT collision operator is used, as it is easier to implement, can be fine-tuned with a single parameter, and is about 20–30% faster overall.

For this test case and a Rayleigh number of Ra = $1\times {10}^5$, the stability limits during the parametric study were observed as:

$$\begin{array}{cc}\hfill BGK& 0.58\le \tau \le 4.0\hfill \\ \hfill TRT-A& 0.56\le \tau \le 5.9\hfill \\ \hfill TRT-S& 0.56\le \tau \le 5.9\hfill \\ \hfill MRT& 0.56\le \tau \le 7.0\hfill \end{array}$$

As expected, the BGK operator is slightly less stable than the TRT and MRT. Surprisingly, the $TRT-S$ operator, which is fine-tuned for stability, demonstrated no clear advantage compared to the $TRT-A$ (fine-tuned for accuracy). The MRT collision operator offers similar stability to TRT at the lower end but is more stable for larger relaxation times. However, it does not matter since the error is already intolerably high for $\tau =5.9$.

4.1.3. Influence of the CHT Method on the Deviation

In the second part of the investigation, the relaxation time for the energy equation of the solids ${\tau }_{gs}$ was varied. For the sake of simplicity, the fluid’s relaxation times for the Navier–Stokes equation and the energy equation were set equal, i.e., ${\tau }_f={\tau }_{gf}$. This results in an implicit variation in the ratios of the physical properties $r_{\lambda }={\lambda }_s/{\lambda }_f$, and $r_c={\left(\rho c\right)}_s/{\left(\rho c\right)}_f$.

In the case where the product of density and heat capacity is smaller for the fluid than for the solid, which is typically valid if the fluid is a gas, i.e., ${\left(\rho c\right)}_f<{\left(\rho c\right)}_s\to {\left(\rho c\right)}_0={\left(\rho c\right)}_f$, the relationships between the thermal relaxation times, the fluid relaxation time ${\tau }_f$, and the ratio of the thermal conductivities $r_{\lambda }$ are as follows:

$$\begin{array}{cc}\hfill {\tau }_{gf}& =\frac{1}{\mathit{Pr}}{\tau }_f+\frac{\mathit{Pr}-1}{2\mathit{Pr}}\hfill \\ \hfill {\tau }_{gs}& =\frac{r_{\lambda }}{\mathit{Pr}}{\tau }_f+\frac{\mathit{Pr}-r_{\lambda }}{2\mathit{Pr}}\hfill \end{array}$$

(36)

Note that ${\tau }_f={\tau }_{gf}$ implies $\mathit{Pr}=1$ and that the ratio $r_c$ does not appear in Equation (36). From this variation, it is possible to determine a range of possible material combinations for a stable and accurate LBM simulation.

Based on the findings from the last section (Section 4.1.2), four different relaxation times for the fluid were selected: ${\tau }_f=1.5,1.2,0.9,0.6$. The upper and lower limits are based on two assumptions: First, the error should remain small ($\approx <5\%$). Second, the physical time step should not be too small.

Figure 10 shows the relative RMSD from the corresponding FVM simulation for the BGK, TRT, and MRT collision operators using the enthalpy-based CHT method of Lu et al. [24]. During a sweep of the relaxation time of the solid ${\tau }_{gs}$, the fluid relaxation times ${\tau }_f$ and ${\tau }_{gf}$ were kept constant.

Two conclusions can be drawn from Figure 10. First, the base error is determined by the fluid relaxation time ${\tau }_f$. Second, starting from this base deviation, the RMSD of the simulation grows rapidly when approaching the limit of ${\tau }_{gs}=0.5$ and continuously decreases as the relaxation time ${\tau }_{gs}$ becomes larger. This is surprising since in the standard LBM (without an additional CHT method), the error becomes smaller approaching ${\tau }_f=0.5$, but at the cost of a much smaller time step and stability issues.

Interestingly, the BGK operator performs very well for relaxation times less than one for the fluid part (${\tau }_f<1.0$). However, the deviation in the BGK depends much more strongly on the fluid relaxation time than TRT and MRT, which is an unpleasant behavior. This is most likely due to the viscosity dependence of the location of the wall when using the half-way bounce-back boundary condition, which directly affects the approximation of the wall-normal temperature gradient at the hot wall. The MRT collision operator provides a slightly better solution but at a higher computational cost of approximately 20–30%.

4.1.4. Realistic Physical Properties

Finally, turning to realistic thermophysical properties reveals an issue that is less linked to the CHT method but more to the LBM itself. An investigation of the realistic properties defined in Table 2 at a Rayleigh number of $Ra=1\times {10}^5$ would make no sense since this would result in very small temperature differences $\Delta T$ with no practical relevance (see

Table 7. Rayleigh numbers and temperature differences for the fictitious fluid and the realistic material pairings air/wool, air/steel, and flue gas/steel.

	Fictitious	Air/Wool	Air/Steel	Flue Gas/Steel
g	10	9.81	9.81	9.81	$\mathrm{m}/{\mathrm{s}}^2$
$Ra$	1 $\times {10}^5$	1 $\times {10}^5$	1 $\times {10}^5$	1 $\times {10}^5$	−
$\Delta T$	20	0.02804	0.001116	0.2771	$\mathrm{K}$
$\Delta T$	20	20	20	20	$\mathrm{K}$
$Ra$	1 $\times {10}^5$	7.13 $\times {10}^7$	1.33 $\times {10}^9$	2.35 $\times {10}^8$	−
$\beta \Delta T$	0.01	0.03	0.06	0.03	−

).

Therefore, for the assessment of the LBM for realistic parameters, the temperature difference was fixed to $\Delta T=20\mathrm{K}$. This leads to substantially higher Rayleigh numbers, which are listed in Table 7. Further increasing the temperature differences would lead to even higher Rayleigh numbers, thus requiring a turbulence model. As turbulence is an inherently three-dimensional phenomenon, it is widely agreed that two-dimensional LES is non-physical. Therefore, the temperature spread was not further increased as the selected boundary conditions already pushed the boundaries of what was demonstrated in the open literature before. However, it is planned to extend the solver into three dimensions in a future work where a LES sub-grid model can be implemented with relative ease.

Table 8. Stability matrix (✔ stable, x unstable) for realistic thermophysical property test cases for different cell counts and relaxation times ${\tau }_f$. Note that ${\tau }_{gf}$ and ${\tau }_{gs}$ are calculated accordingly.

	Air/Wool			Air/Steel			Flue Gas/Steel
	$\mathit{Ra}=\mathbf{7}.\mathbf{13}\times {\mathbf{10}}^{\mathbf{7}}$			$\mathit{Ra}=\mathbf{1}.\mathbf{33}\times {\mathbf{10}}^{\mathbf{9}}$			$\mathit{Ra}=\mathbf{2}.\mathbf{35}\times {\mathbf{10}}^{\mathbf{8}}$
$N_x\times N_y/{\tau }_f$	0.6	0.9	1.5	0.6	0.9	1.5	0.6	0.75	0.9	1.5
280 × 280	✔	x	x	x	x	x	x	x	x	x
560 × 560	✔	✔	x	x	x	x	x	x	x	x
1120 × 1120	✔	✔	x	x	x	x	x	✔	x	x
2240 × 2240	✔	✔	✔	✔	x	x	x	✔	✔	✔

shows the stability of the test cases with realistic physical properties. First, it must be said that only the air/wool case could be simulated with the coarse grid of 280 by 280 lattices used for all results presented so far. Increasing the resolution to 560 by 560 allowed a stable simulation of the air/wool test case with a relaxation time of ${\tau }_f=0.9$. For a fluid relaxation time of ${\tau }_f=1.5$, the mesh needs to be refined further. A resolution of 2240 by 2240 yielded a stable result. However, because of the extremely high computational cost associated with such a fine mesh, such a high relaxation time is not advisable.

For the flue gas/steel test case, the lowest relaxation time ${\tau }_f$ for which a stable simulation could be obtained was ${\tau }_f=0.75$. This is because, for the flue gas/steel test case, the Prandtls number is greater than unity. From Equation (36), it is known that the relaxation time for the energy equation is smaller than for Navier–Stokes (${\tau }_{gf}<{\tau }_f$). Choosing ${\tau }_f=0.75$ results in ${\tau }_{gf}=0.6$ for the flue gas/steel test case. From that, it can be followed that the smallest relaxation time of the bunch (${\tau }_f$, ${\tau }_{gf}$, and ${\tau }_{gs}$) defines the stability of the simulation.

For the air/steel test case, the mesh had to be refined to a resolution of 2240 by 2240 latices, which allowed for a stable simulation. With a Rayleigh number of $Ra=1.33\times {10}^9$, the flow is probably in the transitional region and could benefit from a turbulence model, but this is beyond the scope of the present work.

4.2. Cavity with a Heating Pad

The simulation was run using the TRT collision operator and a relaxation time ${\tau }_f=0.55$ until a simulation time of t = 25,000$\mathrm{s}\approx 6.94\mathrm{h}$ was reached. The relaxation time was slightly reduced below our recommendation in the last section because, otherwise, the local Mach number for this test case would have been too large for accurate simulations. Figure 11 shows the final temperature field for the entire domain, and the velocity field zoomed into the cavity along with the streamlines. The flow and temperature field inside the cavity is quite similar to Rayleigh–Bernard convection and exhibits the counter-rotating vortices characteristic of natural convection, where the hot air rises in the center and, as it cools on the walls at the top of the cavity, falls back down to the heating pad.

However, the transient behavior is the most significant regime for this type of problem. Figure 12a shows a comparison of the temperature profile at the centerline (x = const.) obtained with the LBM and FVM, respectively, after a simulation time of $t=7200\mathrm{s}$. The results show that the temperature profiles agree very well, but the LBM simulation is slightly hotter. This can also be seen in Figure 12b, which shows the time trace of a temperature probe at position x = center, $y=200\mathrm{m}\mathrm{m}$, indicating that the LBM simulation heats up marginally quicker, but both methods agree closely once the system reaches steady-state conditions. The difference is most likely due to the implementation of the heating-pad source term in the energy equation of the LBM. Further research is needed to improve the implementation. However, the results are very promising and support our conclusion that LBM equipped with proper CHT methods can be an alternative to standard CFD techniques for this kind of application.

5. Conclusions and Outlook

Two different test cases were investigated. First, a differentially heated cavity filled with solid blocks was used to explore the stability and accuracy of the unmodified double-distribution LBM. Without modification, the thermal relaxation times for the solid and fluid are the same (${\tau }_{gs}={\tau }_{gf}$), resulting in equal thermal diffusivities $a_f=a_s$. For simplicity, the fluid-flow relaxation time ${\tau }_f$ was chosen equal to the thermal relaxation times, implicitly setting the Prandtl number to unity $\left(\mathit{Pr}=1\right)$. It was confirmed that the accuracy of the LBM comes as a compromise between simulation time and stability. When the relaxation time is decreased, the simulation usually becomes more accurate, but at the expense of stability and computational effort. Above a certain value for the relaxation time, the results become implausible. According to our research, in similar cases, it would be advisable to use a relaxation time in the range of $0.6\le {\tau }_f\le 1.5$. It was also confirmed that the TRT collision operator represents the best compromise regarding accuracy, stability, and computation time compared to BGK and MRT.

Next, a variation in the solid relaxation time ${\tau }_{gs}$ was performed at four different fluid relaxation times ${\tau }_f$. It was found that the principal source of error concerns the fluid solver and depends on the relaxation time ${\tau }_f$. Besides that, increasing the relaxation time of the solid ${\tau }_{gs}$ leads to a more accurate result, which is surprising because it is the other way around for standard LBM without CHT. Further investigations are required to understand this seemingly contradicting behavior better.

If the temperature difference is increased to practical values, the Rayleigh number increases drastically, which requires refining the mesh significantly. This is a fact often overlooked by previous LBM research, where the restriction to small Ra numbers (Ra = $1\times {10}^5$) would lead to very small temperature differences, well below 1 K, for realistic thermophysical properties. Subsequently, and for test cases of greater practical importance, a turbulence model would have to be used, which is beyond the scope of the present work. However, the second test case demonstrated that the LBM provides accurate and fast predictions, especially if the transient solution is of interest. Still, the main limitations, especially for processes with large temperature differences, are the temporal constant thermophysical properties.

In future work, it is planned to extend the solver to three dimensions, allowing the inclusion of a turbulence model in the form of a LES sub-grid model with relative ease. This will allow consideration of test cases with higher Rayleigh numbers and investigation if a turbulence model can further stabilize the CHT model.