A Dirichlet boundary condition for the thermal lattice Boltzmann method

Highlights

• A boundary condition for thermal lattice Boltzmann simulations is introduced.

• The method is second order accurate.

• The method is well suited for particles containing systems.

Abstract

In this work we introduce a boundary condition for thermal lattice Boltzmann simulations that contain a Dirichlet boundary condition by bouncing back the non-equilibrium distribution of the energy distribution function. To this end the thermal lattice Boltzmann equation is modified by introducing an additional collision term that takes into account the thermal diffusivity and local solid volume fraction of a lattice (partially) covered by the solid phase. Asymptotic analysis of the boundary condition confirms that it is of second order accuracy. The method is validated using (i) an analytical solution for the Nusselt number correlation of a single sphere in an unbounded stationary fluid and (ii) direct numerical simulations of the heat transfer between a fluid and individual particles.

Introduction

Particulate two phase, non-isothermal flows are encountered frequently in industrial applications, e.g. gas-fluidized bed reactors. The understanding of the complicated interaction between fluid mechanics and thermal effects is essential to improve further the design of reactors that are commonly applied in the chemical industry and, hence, to improve their efficiency. In the past decades, computational fluid dynamics (CFD) has been applied extensively to compute the fluid flow, heat transfer and chemical reactions in complex systems, e.g. packed or fluidized reactors. However, the flow structure in these systems spans a wide range of length scales, e.g. from (sub-) millimeters (particle level) to meters (large bubbles or slugs). For example, Freund et al. (2003) simulated single phase, reacting flows in a randomly packed bed using a lattice Boltzmann method and found that integral quantities such as pressure drop were influenced significantly by the local packing structure, a feature which is often neglected in volume-averaged CFD simulations. In contrast, in direct numerical simulations (DNS), the flow around individual particles is fully resolved and interactions between the fluid and particles are taken into account by imposing the appropriate boundary conditions at the particles’ surface. Additionally, DNS has the potential to complement experiments as DNS provides a versatile tool to extract data that is otherwise extremely difficult, if not impossible, to obtain experimentally. Among the different DNS methods available, the lattice Boltzmann method has several advantages when compared to “conventional” methods such as the finite volume method (Chen and Doolen, 1998), viz. (1) the convective and collisional operators are linear; (2) the fluid pressure can be calculated simply by an equation of state; (3) complex boundaries are relatively easy to implement and (4) ease of parallelization as computations are local.

Recently the extension of LBM to thermal flows has received significant attention (Yoshida and Nagaoka, 2010; Li et al., 2013; Pareschi et al., 2016; Wang et al., 2016). In general, three different approaches have been proposed: (i) a multispeed model (Alexander et al., 1993), (ii) a passive-scalar model (Bartoloni et al., 1993; Shan, 1997) and (iii) a double population model (He et al., 1998; Guo et al., 2007).

Alexander et al. (1993) developed the multispeed approach which requires only the density distribution function. To recover the macroscopic energy equation, additional lattice speeds and equilibrium density distribution functions including higher order velocity terms are required. A disadvantage of this approach is that it suffers from numerical instabilities for high Rayleigh number flows and that the Prandtl number cannot be varied independently since a single relaxation time is used for the collision operation (McNamara et al., 1995). In the passive-scalar model (Bartoloni et al., 1993; Shan, 1997), the temperature is described as a passive scalar and simulated using a separate distribution function which is independent of the density distribution function. This approach enhances the numerical stability when compared to the multispeed model (Guo et al., 2002). The main disadvantage of the passive scalar model is that the viscous dissipation and compression work done by the pressure is not taken into account. He et al. (1998) proposed a double-population model at the incompressible limit. In the double population model, the internal energy distribution function is derived directly from the second order moment of the continuous Boltzmann equation. The shortcoming of this model is that the energy evolution equation includes temporal and spatial derivatives of the macroscopic fluid velocity, which increases the complexity of its implementation. More recently, Guo et al. (2007) proposed a double population model using a total energy distribution function. Similar to the model proposed by He et al. (1998), a new variable (functions of the energy and density distributions) is introduced to avoid the implicitness of the scheme since the original energy evolution equation does not compute explicitly the energy distribution functions. However, the simple bounce back condition has to be applied to the “old” variables, i.e. the energy distribution function. Hence, a clear disadvantage of this approach is that the implementation of the boundary conditions becomes cumbersome.

To overcome some of the shortcomings of the double-population model, Peng et al. (2003) proposed a simplified thermal lattice Boltzmann model for incompressible thermal flow. The work reported here is based on this model. Eqs. (1) and (2) are the governing equations for the density and energy distribution functions, respectively (as used by Peng et al. (2003)):

$$f_i(\overrightarrow{x}+{\overrightarrow{c}}_i\Delta t,t+\Delta t)=f_i(\overrightarrow{x},t)-\frac{1}{{\tau }_v}[f_i(\overrightarrow{x},t)-f_i^{eq}(\overrightarrow{x},t)]$$

$$g_i(x+{\overrightarrow{c}}_i\Delta t,t+\Delta t)=g_i(x,t)-\frac{1}{{\tau }_c}(g_i(x,t)-{g_i}^{eq}(x,t))$$

where f_i and g_i are, respectively, the density and energy distribution functions, ${\overrightarrow{c}}_i$ is the discrete lattice velocity and $\overrightarrow{x}$ and t are position and time, respectively. The superscript ‘eq’ denotes the equilibrium state, τ_v and τ_c are the single relaxation times related to the viscosity and thermal diffusivity, respectively and Δt is the time step applied. The model of Peng et al. (2003) is based on the assumption that the viscous dissipation and compression work done by the pressure can be neglected in incompressible flows. The numerical simulation of natural convection in a square cavity confirmed that the simplified thermal model of Peng et al. (2003) was of similar accuracy as the model proposed by He et al. (1998) which included viscous dissipation and compression work done by pressure. The D3Q19 and D3Q7 lattices (Fig. 1) for the density and energy distributions are given as follows:

$${\overrightarrow{c}}_i=c{\overrightarrow{e}}_i=\{\begin{array}{c}(0,0,0),i=0\hfill \\ (\pm 1,0,0)c,(0,\pm 1,0)c,(0,0,\pm 1)c,i=1...6\hfill \\ (\pm 1,\pm 1,0)c,(\pm 1,0,\pm 1)c,(0,\pm 1,\pm 1)c,i=7...18\hfill \end{array}$$

with c being the ratio of Δx to Δt. The discrete equilibrium density distribution function is expressed as

$$\begin{array}{c}{f}_i^{eq}=w_{i,v}\rho [1+3\frac{{\overrightarrow{c}}_i·\overrightarrow{u}}{c^2}+\frac{9}{2}\frac{{({\overrightarrow{c}}_i·\overrightarrow{u})}^2}{c^4}-\frac{3}{2}\frac{{\overrightarrow{u}}^2}{c^2}]\hfill \\ w_{i,v}=\{\begin{array}{c}\frac{1}{3},i=0\hfill \\ \frac{1}{18},i=1...6\hfill \\ \frac{1}{36},i=7...18\hfill \end{array}\hfill \end{array}$$

where $\overrightarrow{u}(\overrightarrow{x},t)$ is the velocity of the fluid. The discrete equilibrium energy distribution is:

$$\begin{array}{c}{g}_i^{eq}=w_{i,c}\varphi \left(1+4·\frac{{\overrightarrow{c}}_i·\overrightarrow{u}}{c^2}\right)\hfill \\ w_{i,c}=\{\begin{array}{c}\frac{1}{4},i=0\hfill \\ \frac{1}{8},i=1...6\hfill \end{array}\hfill \end{array}$$

The macroscopic properties such as density, momentum and energy can be obtained through the following summations:

$$\rho ={\sum }_i{f}_i,\rho \overrightarrow{u}={\sum }_i{f}_i{\overrightarrow{e}}_i,\varphi ={\sum }_i{g}_i,$$

The viscosity and thermal diffusivity are calculated as:

$$\begin{array}{c}\nu =\frac{1}{3}{c}^2\left({\tau }_v-\frac{\Delta t}{2}\right)\hfill \\ D=\frac{1}{4}{c}^2\left({\tau }_c-\frac{\Delta t}{2}\right)\hfill \end{array}$$

A key aspect in every numerical scheme is the accurate description of the boundary condition, in particular for multi-boundary systems such as packed or fluidized beds. For the energy equation, the method that is easiest to implement in thermal LBM is the bounce back scheme.

He et al. (1998) extended the bounce back rule introduced by Zou and He (1997) for thermal boundary conditions. Ginzburg (2005) proposed a multi-reflection approach for the Dirichlet boundary condition. When compared to the bounce back or linear interpolation method, this approach yields smoother macroscopic fluid properties (e.g. velocity and pressure) in the vicinity of the boundary (Ginzburg and d'Humieres, 2003). In the multi-reflection method, the equilibrium energy distribution is decomposed into symmetric and anti-symmetric parts. The symmetric part is tuned to achieve a second or third order accurate Dirichlet boundary condition by expanding the unknown energy distribution function at the boundary. With regards to the multiple relaxation time-thermal lattice Boltzmann method, Li et al. (2013) developed a bounce-back-based interpolation method for the Dirichlet boundary condition using three post-collision energy distributions. Li et al. (2013) demonstrated that the Dirichlet boundary condition was second order accurate and that the scheme cannot be determined uniquely since the moments of the leading-order solution of the energy distribution is sufficient to satisfy the Dirichlet boundary condition. The second order Dirichlet boundary condition proposed by Ginzburg (2005) can be considered as one solution of Li's boundary condition.

In addition to bounce back schemes, Guo et al. (2002) proposed a non-equilibrium extrapolation method for the double population model by decomposing the unknown energy distribution into an equilibrium and a non-equilibrium part at the boundary nodes. In this method, the equilibrium part of the energy distribution at the boundary node is calculated using the known boundary temperature, whereas the non-equilibrium part is extrapolated from neighboring fluid nodes. To improve the accuracy of the boundary condition, cut cell or sub-grid information has been introduced (Ginzburg, 2005; Li et al., 2013). However, the computational cost and the selection of appropriate interpolation schemes have to be assessed carefully when addressing complex boundaries.

In this paper, the hydrodynamic boundary condition proposed by Noble and Torczynski (1998) for isothermal flows is extended to the thermal lattice Boltzmann method. The original method for the hydrodynamic boundary condition can handle with great ease boundary dominated problems (e.g. porous media) in which the solid-fluid interface does not conform to the grid. The thermal boundary condition proposed here takes into account the fluid diffusivity and sub-grid effects such as fluid cells being partially covered by particles. Unlike the boundary conditions mentioned above, the new thermal boundary condition is of higher accuracy than the half way bounce back scheme, and less expensive than the multi-reflection and bounce-back-based interpolation methods. Moreover, the newly proposed boundary condition is purely local, whereas the multi-reflection and bounce-back-based interpolation methods need at least the information of two fluid nodes to achieve the desired accuracy at the boundary.