Interface Treatment for Conjugate Conditions in the Lattice Boltzmann Method for the Convection Diffusion Equation

Abstract

The lattice Boltzmann method (LBM) has emerged as an attractive numerical method for fluid flows and thermal and mass transport. For LBM modeling of transport between different phases or materials of distinct properties, effective treatment for the conjugate conditions at the interface is required. Recognizing the benefit of satisfying the conjugate conditions in each time step without iterative computations using LBM, various interface schemes have been proposed in the last decade. This chapter provides a review of those interface schemes, with a focus on the comparison of numerical accuracy and convergence orders. It is shown that in order to preserve the second-order accuracy in LBM, the local interface geometry must be considered; and the modified geometry-ignored interface schemes result in degraded convergence orders and/or much higher error magnitude. It is also verified that with appropriate interface schemes, interfacial transport with scalar and flux jumps can be effectively modeled.

1. Introduction

Heat and mass transfer between multi-phases or different materials with interfacial conjugate conditions is frequently encountered in fundamental sciences and numerous engineering applications involving fluid dynamics, thermal transport, materials sciences, and chemical reactions. Examples are cooling of turbine blades, heat exchangers and electronic devices, thermal insulation on heat pipes and chemical reactors, heat conduction in composite materials, and heat and mass transfer between solid particles and their surrounding fluids [1, 2, 3, 4, 5, 6, 7, 8], to name a few. The most well-known conjugate conditions include the continuity of both the temperature (concentration) and the heat (mass) flux at the interface. Other conjugate conditions, such as with temperature (concentration) jumps and/or flux discontinuities [9], and Henry’s law relationship [10], are also noticed at fluid-solid interfaces or interfaces of two solids or fluids of different thermal (mass diffusion) properties.

The non-smoothness or discontinuities/jumps in the physical or transport properties, and consequently in the distribution of the temperature (concentration) field across the interface, pose a great challenge to any numerical method applied to solve the interface problems. Development of accurate and efficient numerical schemes to treat the interface conditions has attracted much attention in the literature, such as the immersed boundary method (IBM) [11, 12], the immersed interface method (IIM) [13, 14], the ghost fluid method (GFM) [15, 16], the sharp interface Cartesian grid method [17, 18], and the matched interface and boundary (MIB) method [19]. Most of these methods are formulated in the finite-difference, finite-volume or finite-element frameworks.

When applying those traditional numerical methods, a popular approach to implement the conjugate conditions is to employ iterative schemes, in which a Dirichlet interface condition is imposed for one phase or material and a Neumann interface condition for the other. The heat and mass transfer in each phase is separately solved, and the continuity or prescribed jump condition at the interface could be satisfied after multiple iterations. For conjugate transport with complex interface geometry, the iterative schemes would become difficult to implement and they normally necessitate a considerable amount of computational effort.

The lattice Boltzmann method (LBM), which has emerged as an attractive alternative numerical method for modeling fluid flows and heat mass transfer (see [20, 21, 22] and Refs. therein), has been demonstrated to be an effective and efficient numerical approach for conjugate interface conditions in tandem with the convection diffusion equation (CDE) [9, 23]. In this chapter, we present a critical review of the various interface schemes proposed in the literature, with a focus on the comparison of numerical accuracy.

The well-known features of the LBM method include its explicit algorithm, ease in implementation, capability to treat complex geometry, and compatibility with parallel computing [20, 21]. Boundary condition treatment is essential to the integrity of LBM since the kinetic theory-based method deals directly with the microscopic distribution functions (DFs) rather than the macroscopic conservation equations. Earlier LB models treat the collision effects with a single-relaxation-time (SRT) approximation, commonly referred to as the Bhatnagar-Gross-Krook (BGK) model [24, 25, 26]. However, the SRT model is limited such that it can only describe isotropic diffusion [20]. In recent years, models such as the two-relaxation-time (TRT) [27, 28] and multiple-relaxation-time (MRT) [20, 29, 30] LB models have been proposed that can handle anisotropic diffusion. Representative LB models proposed in the literature include the general BGK model by Shi and Guo [31] for the nonlinear CDE, the D3Q7/D2Q5 MRT models by Yoshida and Nagaoka [20] for the general convection anisotropic diffusion equation, and the D1Q3/D2Q9/D3Q19 MRT models by Chai and Zhao [30] for the general nonlinear convection anisotropic diffusion equation, to name a few. The MRT models have improved numerical accuracy and stability compared to the SRT models [20, 28, 30]. The D3Q7/D2Q5 model proposed in [20] is used for this review, as it preserves second order spatial accuracy when recovering the general CDE following an asymptotic analysis. Based on the D3Q7/D2Q5 LB models, Li et al. [21] proposed second-order accurate boundary treatments for both the Dirichlet and Neumann conditions; they have also established a general framework for heat and mass transfer simulations with direct extension to curved boundary situations. In their framework, explicit analytical expressions were developed to relate the macroscopic quantities, such as boundary temperature (concentration) and their fluxes, and interior temperature (concentration) gradients, to the microscopic DFs in the LB model.

The first work that explicitly addressed the fluid-solid interface condition in LBM was conducted by Wang et al. [6]. They proposed a simple “half lattice division” (HLD) treatment in which no special treatment is required and the temperature and flux continuity condition at the interface was automatically satisfied for steady cases. Improvement of this HLD based scheme was conducted by several groups (see a review in Ref. [23]) for unsteady cases. Importantly, with the interface treated as a shared boundary between the adjacent domains, the boundary conditions by Li et al. [21] were applied to interface conditions and particular interface schemes were proposed and verified in [23] for standard conjugate conditions, and in [9] for conjugate heat and mass transfer with interfacial jump conditions. This idea of developing analytical relationships for the DFs to satisfy the conjugate conditions was also extended to handle general interface conditions in [32, 33]. In all the previous schemes in [9, 23, 32, 33], the local geometry was taken into account and the second-order accuracy of the LBM solution can be preserved.

There is another category of interface schemes that has attracted interest in the LBM community. In those schemes, additional source terms [34], alternative LBE formulations [35, 36], or modified equilibrium DFs [37, 38], were proposed to handle the conjugate conditions. The main motivation for those schemes is to avoid the consideration of the interface geometry or topology, which can be a challenge in complex systems such as porous media. As pointed out in [33], however, these schemes usually suffer from degraded numerical accuracy and/or convergence orders. This perspective will be illustrated in detail in this chapter.

The rest of this chapter is organized as follows. In Section 2, the different types of conjugate conditions are presented. The LB model for the general CDE is briefly described in Section 3. In Section 4, the representative interface schemes are summarized. Discussion on the numerical accuracy of the selected interface schemes is provided in Section 5 with representative numerical examples. Conclusions and outlook are given in Section 6.

2. Conjugate conditions in heat and mass transfer

In order to define the conjugate conditions, we begin by defining two domains 1 and 2, as shown in Figure 1. The conjugate interface conditions, from a heat and mass transfer perspective, can be defined as:

Figure 1.

Illustration of the local geometry of an interface in the lattice (filled circles: lattice nodes in domain 1, filled squares: interface nodes, and open circles: lattice nodes in domain 2). With permission from [9].

$${\Phi }_{\mathit{df}}={\Phi }_{\mathit{ds}}+{\varphi }_{\text{jump}}$$E1a

$$\mathbf{n}·{\left(k\nabla \varphi +\rho c_p\mathbf{u}\varphi \right)}_f=\mathbf{n}·{\left(k\nabla \varphi +\rho c_p\mathbf{u}\varphi \right)}_s+\mathbf{n}·{\mathbf{q}}_{\text{jump}}^T\text{in heat transfer},$$E1b

$$\mathbf{n}·{\left(D_m\nabla \varphi +\mathbf{u}\varphi \right)}_f=\mathbf{n}·{\left(D_m\nabla \varphi +\mathbf{u}\varphi \right)}_s+\mathbf{n}·{\mathbf{q}}_{\text{jump}}^C\text{in mass transfer},$$E1c

where n represents the normal direction, $\varphi$ the macroscopic scalar variable of interest (temperature or concentration), k the thermal conductivity, $\rho$ the density, $c_p$ the heat capacity, $\mathbf{u}$ the velocity vector, D_m the mass diffusivity, and ${\varphi }_{\text{jump}}$ and ${\mathbf{q}}_{\text{jump}}$ the jump conditions at the interface.

Some examples of jump conditions can be found in cases such as concentration jumps (Henry’s law [10]) or temperature jumps at the interface [9]. Eqs. (1a)–(1c) reduce to the standard conjugate conditions in [23] with no jumps and zero normal velocity; they can also be extended to yield two general relationships between interfacial scalar values and their fluxes as in [32, 33].

3. Lattice Boltzmann model for the general CDE

The governing heat and mass transfer equation within each domain can be written as a general convection diffusion equation as:

$$\frac{\partial \varphi }{\partial t}+\frac{\partial }{\partial x_j}\left({\upsilon }_j\varphi \right)=\frac{\partial }{\partial x_j}\left(D_{\mathit{ij}}\frac{\partial \varphi }{\partial x_j}\right)+G$$E2

where $D_{\mathit{ij}}$ represents the diffusion coefficient, and $G$ is the general source term.

For fluid flow simulations, the D2Q9/D3Q19 LB models are the most popular selections due to their accuracy and robustness [39]. While for the scalar CDE (2), the D2Q5/D3Q7 LB models are most widely used [20, 40]. To recover the CDE to second-order accuracy, the evolution equation follows

$$g_{\alpha }\left(\mathbf{x}+{\mathbf{e}}_{\alpha }\mathit{\delta t},t+\mathit{\delta t}\right)-g_{\alpha }\left(\mathbf{x},t\right)={\left[\mathbf{L}·\left(\mathbf{g}-{\mathbf{g}}^{\mathrm{eq}}\right)\left(\mathbf{x},t\right)\right]}_{\alpha }+{\omega }_{\alpha }G\left(\mathbf{x},t\right)\mathit{\delta t}$$E3

where the microscopic distribution function, g_α(x, t) $\equiv$ g(x, ξ_α, t), is defined in the discrete velocity space, ξ is the particle velocity vector that is discretized to a small set of discrete velocities {ξ_α|α = 0, 1, …, m $-$ 1}, e_α is the αth discrete velocity vector, δt is the time step, L is the collision operator, $g_{\alpha }^{\mathrm{eq}}\left(\mathbf{x},t\right)$ is the equilibrium distribution function, and ω_α is the weight coefficient. The macroscopic scalar variable is obtained from

$$\varphi \left(\mathbf{x},t\right)=\sum _{\alpha =0}^{m-1}{g}_{\alpha }\left(\mathbf{x},t\right)$$E4

The equilibrium distribution function can be defined as [20, 40]

$$g_{\alpha }^{\mathrm{eq}}={\omega }_{\alpha }\varphi \left(1+\frac{{\mathbf{e}}_{\alpha }·\mathbf{u}}{c_s^2}\right)$$E5

where c_s is the speed of sound with $c_s=c/\sqrt{3}=\left(\mathit{\delta x}/\mathit{\delta t}\right)/\sqrt{3}=1/\sqrt{3}$.

When using the multiple-relaxation-time (MRT) collision operator, and the collision-streaming process for efficient computations, Eq. (3) is split into two steps:

Collision step

$${\widehat{g}}_{\alpha }\left(\mathbf{x},t\right)=g_{\alpha }\left(\mathbf{x},t\right)-{\left[{\mathbf{M}}^{‐1}\mathbf{S}\left(\mathbf{m}\left(\mathbf{x},t\right)-{\mathbf{m}}^{\left(\mathrm{eq}\right)}\left(\mathbf{x},t\right)\right)\right]}_{\alpha }+{\omega }_{\alpha }G\left(\mathbf{x},t\right)\mathit{\delta t},\text{and}$$E6

Streaming step

$$g_{\alpha }\left(\mathbf{x}+{\mathbf{e}}_{\alpha }\mathit{\delta t},t+\mathit{\delta t}\right)={\widehat{g}}_{\alpha }\left(\mathbf{x},t\right).$$E7

In the above, M is a matrix to transform the distribution functions ${\mathbf{g}}^{\left(\mathrm{eq}\right)}$ to their moments ${\mathbf{m}}^{\left(\mathrm{eq}\right)}$ by the relation ${\mathbf{m}}^{\left(\mathrm{eq}\right)}=\mathbf{M}{\mathbf{g}}^{\left(\mathrm{eq}\right)}$. S is a matrix of relaxation coefficients ${\tau }_{\mathit{ij}}$. In the D3Q7 model, the equilibrium moments of the distribution functions are ${\mathbf{m}}^{\left(\mathrm{eq}\right)}={\left(\varphi ,\mathit{u\varphi },\mathit{v\varphi },\mathit{w\varphi },\mathit{a\varphi },0,0\right)}^T,$where u, v, and w are the macroscopic velocity components, and $a$ is a constant related to the weight coefficients. For details about the matrices and the constants in the LB models the reader can refer to [20, 32].

4. Interface schemes for conjugate conditions

One unique feature of the LBM method is that both the Dirichlet-type boundary value and the Neumann-type boundary flux, i.e., temperature/concentration gradient, can be obtained from a simple moment of the distribution functions with appropriate boundary schemes [20, 21]. It eliminates finite-difference type approximation schemes for the flux. This idea can also be applied to construct interface schemes by treating the interface as a shared boundary between the two adjacent domains [9, 23, 32, 33]. We consider the basic situation with zero convective flux (u_n = 0) and the normal of the interface parallel to the discrete velocity vector, i.e., parallel straight interface; the interface scheme for more general situations such as curved geometry can be similarly constructed as shown in [9, 23, 32, 33]. The conjugate conditions in Eqs. (1b) and (1c) thus reduce to (see Figure 1)

$${\mathrm{\Phi }}_{\mathit{nf}}=-D_f\frac{\partial {\varphi }_f}{\partial n_f}=\sigma D_s\frac{\partial {\varphi }_s}{\partial n_s}+q_{\text{jump}}=-\sigma {\mathrm{\Phi }}_{\mathit{ns}}+q_{\text{jump}},$$E8

with σ = 1, $q_{\text{jump}}=q_{\text{jump}}^{\mathrm{C}}$ in mass transfer, and $\sigma =\frac{{\left(\rho c_p\right)}_s}{{\left(\rho c_p\right)}_f}$, $q_{\text{jump}}=q_{\text{jump}}^{\mathrm{T}}/{\left(\rho c_p\right)}_f$ in heat transfer.

Depending on whether the local interface geometry is considered, the interface schemes fall into two different categories, as discussed in Sections 4.1 and 4.2.

4.1 Interpolation-based interface schemes

According to the second-order interpolation-based boundary schemes developed in [21], the relationships between the distribution functions and the interfacial ϕ values and their fluxes for each domains can be obtained: for the Dirichlet condition treatment,

$$g_{\overline{\alpha }}\left({\mathbf{x}}_f,t+\mathit{\delta t}\right)=c_{d1}{\widehat{g}}_{\alpha }\left({\mathbf{x}}_f,t\right)+c_{d2}{\widehat{g}}_{\alpha }\left({\mathbf{x}}_{\mathit{ff}},t\right)+c_{d3}{\widehat{g}}_{\overline{\alpha }}\left({\mathbf{x}}_f,t\right)+c_{d4}{\epsilon }_D{\mathrm{\Phi }}_{\mathit{df}},$$E9a

$$g_{\alpha }\left({\mathbf{x}}_s,t+\mathit{\delta t}\right)=c_{d1}^{\ast }{\widehat{g}}_{\overline{\alpha }}\left({\mathbf{x}}_s,t\right)+c_{d2}^{\ast }{\widehat{g}}_{\overline{\alpha }}\left({\mathbf{x}}_{\mathit{ss}},t\right)+c_{d3}^{\ast }{\widehat{g}}_{\alpha }\left({\mathbf{x}}_s,t\right)+c_{d4}^{\ast }{\epsilon }_D{\mathrm{\Phi }}_{\mathit{ds}}.$$E9b

Similarly, for the Neumann condition treatment,

$$g_{\overline{\alpha }}\left({\mathbf{x}}_f,t+\mathit{\delta t}\right)=c_{n1}{\widehat{g}}_{\alpha }\left({\mathbf{x}}_f,t\right)+c_{n2}{\widehat{g}}_{\alpha }\left({\mathbf{x}}_{\mathit{ff}},t\right)+c_{n3}{\widehat{g}}_{\overline{\alpha }}\left({\mathbf{x}}_f,t\right)+c_{n4}\left(\mathit{\delta t}/\mathit{\delta x}\right){\mathrm{\Phi }}_{n\overline{\alpha }},$$E10a

$$g_{\alpha }\left({\mathbf{x}}_s,t+\mathit{\delta t}\right)=c_{n1}^{\ast }{\widehat{g}}_{\overline{\alpha }}\left({\mathbf{x}}_s,t\right)+c_{n2}^{\ast }{\widehat{g}}_{\overline{\alpha }}\left({\mathbf{x}}_{\mathit{ss}},t\right)+c_{n3}^{\ast }{\widehat{g}}_{\alpha }\left({\mathbf{x}}_s,t\right)+c_{n4}^{\ast }\left(\mathit{\delta t}/\mathit{\delta x}\right){\mathrm{\Phi }}_{\mathit{n\alpha }},$$E10b

where ${\mathrm{\Phi }}_{n\overline{\alpha }}$ and ${\mathrm{\Phi }}_{\mathit{n\alpha }}$ are the respective interfacial fluxes along the discrete lattice velocity directions ${\mathbf{e}}_{\overline{\alpha }}$ and ${\mathbf{e}}_{\alpha }$, x_f and x_ff are the first and second interior lattice nodes along ${\mathbf{e}}_{\overline{\alpha }}$ direction in Domain 1 (x_ff = x_f + ${\mathbf{e}}_{\overline{\alpha }}$ δt), and x_s and x_ss are the lattice nodes along ${\mathbf{e}}_{\alpha }$ direction in Domain 2 (x_ss = x_s + ${\mathbf{e}}_{\alpha }$ δt = x_s − ${\mathbf{e}}_{\overline{\alpha }}$ δt), respectively (see Figure 1). All the coefficients in Eqs. (9) and (10) are only related to the local geometry as denoted by the link intersection fraction, Δ, at the interface [21, 23].

When ${\mathbf{e}}_{\overline{\alpha }}$ is aligned with the interface normal directions, ${\mathrm{\Phi }}_{n\overline{\alpha }}={\mathrm{\Phi }}_{\mathit{nf}}$ and ${\mathrm{\Phi }}_{\mathit{n\alpha }}={\mathrm{\Phi }}_{\mathit{ns}}$ are readily noticed. Hence, Eqs. (1), (8), (9a), (9b), (10a) and (10b) constitute a linear system of six equations, and the six unknowns $g_{\overline{\alpha }}\left({\mathbf{x}}_f,t+\mathit{\delta t}\right)$, $g_{\alpha }\left({\mathbf{x}}_s,t+\mathit{\delta t}\right)$, ${\mathrm{\Phi }}_{\mathit{df}}$, ${\mathrm{\Phi }}_{\mathit{ds}}$, ${\mathrm{\Phi }}_{\mathit{nf}}$, and ${\mathrm{\Phi }}_{\mathit{ns}}$ can be analytically solved. The interface scheme thus becomes [9]:

$$\begin{array}{l}{g}_{\overline{\alpha }}\left({\mathbf{x}}_f,t+\mathit{\delta t}\right)=A_1^f{\widehat{g}}_{\alpha }\left({\mathbf{x}}_f,t\right)+A_2^f{\widehat{g}}_{\alpha }\left({\mathbf{x}}_{\mathit{ff}},t\right)+A_3^f{\widehat{g}}_{\overline{\alpha }}\left({\mathbf{x}}_f,t\right)\\ +B_1^f{\widehat{g}}_{\overline{\alpha }}\left({\mathbf{x}}_s,t\right)+B_2^f{\widehat{g}}_{\overline{\alpha }}\left({\mathbf{x}}_{\mathit{ss}},t\right)+B_3^f{\widehat{g}}_{\alpha }\left({\mathbf{x}}_s,t\right)\\ +{\chi }_q{q}_{\text{jump}}+{\chi }_t{\varphi }_{\text{jump}},\end{array},$$E11a

$$\begin{array}{l}{g}_{\alpha }\left({\mathbf{x}}_s,t+\mathit{\delta t}\right)=A_1^s{\widehat{g}}_{\overline{\alpha }}\left({\mathbf{x}}_s,t\right)+A_2^s{\widehat{g}}_{\overline{\alpha }}\left({\mathbf{x}}_{\mathit{ss}},t\right)+A_3^s{\widehat{g}}_{\alpha }\left({\mathbf{x}}_s,t\right)\\ +B_1^s{\widehat{g}}_{\alpha }\left({\mathbf{x}}_f,t\right)+B_2^s{\widehat{g}}_{\alpha }\left({\mathbf{x}}_{\mathit{ff}},t\right)+B_3^s{\widehat{g}}_{\overline{\alpha }}\left({\mathbf{x}}_f,t\right)\\ +{\gamma }_q{q}_{\text{jump}}+{\gamma }_t{\varphi }_{\text{jump}}.\end{array}.$$E11b

The coefficients in Eqs. (11a) and (11b) are now determined by the geometry fraction Δ and the property ratio σ. It is worth noting that there is an adjustable parameter in those coefficients since the second-order Dirichlet boundary scheme allows one adjustable parameter, as shown in [21], where three particular Dirichlet schemes were also presented:

$$\text{Scheme}1:c_{d1}=\left\{\begin{array}{l}-2\mathrm{\Delta },\left(0\le \mathrm{\Delta }\le 0.5\right),\\ -\frac{1}{2\mathrm{\Delta }},\left(\mathrm{\Delta }>0.5\right),\end{array}\right.$$E12a

$$\text{Scheme}2:c_{d1}=-2\left(1-\mathrm{\Delta }\right),\text{and}$$E12b

$$\text{Scheme}3:c_{d1}=-1.$$E12c

The corresponding interface schemes, can thus also be obtained. Those schemes will be numerically verified in Section 5.2 for a test case including interfacial jump conditions.

When the straight interface is located “halfway” between the lattice nodes (∆ = 0.5), the unknown DFs can be calculated by only knowing two single-node post-collision DFs, i.e., without interpolation:

$$g_{\overline{\alpha }}\left({\mathbf{x}}_f,t+\mathit{\delta t}\right)=\left(\frac{1-\sigma }{1+\sigma }\right){\widehat{g}}_{\alpha }\left({\mathbf{x}}_f,t\right)+\left(\frac{2\sigma }{1+\sigma }\right){\widehat{g}}_{\overline{\alpha }}\left({\mathbf{x}}_s,t\right)+\frac{q_{\text{jump}}}{1+\sigma }+\frac{{\epsilon }_D{\mathit{\sigma \varphi }}_{\text{jump}}}{1+\sigma },$$E13a

$$g_{\alpha }\left({\mathbf{x}}_s,t+\mathit{\delta t}\right)=-\left(\frac{1-\sigma }{1+\sigma }\right){\widehat{g}}_{\overline{\alpha }}\left({\mathbf{x}}_s,t\right)+\left(\frac{2}{1+\sigma }\right){\widehat{g}}_{\alpha }\left({\mathbf{x}}_f,t\right)+\frac{q_{\text{jump}}}{1+\sigma }-\frac{{\epsilon }_D{\varphi }_{\text{jump}}}{1+\sigma }.$$E13b

Furthermore, for the most simplified case of ∆ = 0.5 and σ = 1, Eqs. (13a) and (13b) reduce to:

$$g_{\overline{\alpha }}\left({\mathbf{x}}_f,t+\mathit{\delta t}\right)={\widehat{g}}_{\overline{\alpha }}\left({\mathbf{x}}_s,t\right)+q_{\text{jump}}/2+{\epsilon }_D{\varphi }_{\text{jump}}/2,$$E14a

$$g_{\alpha }\left({\mathbf{x}}_s,t+\mathit{\delta t}\right)={\widehat{g}}_{\alpha }\left({\mathbf{x}}_f,t\right)+q_{\text{jump}}/2-{\epsilon }_D{\varphi }_{\text{jump}}/2.$$E14b

For straight interfaces where Δ ≠ 0.5 and for curved interfaces, the complete interface conditions can be found in [9, 23].

It should be noted that second-order accurate boundary schemes can also be obtained using only the single lattice node next to the boundary, as demonstrated in [32], instead of using interpolations in Eqs. (9) and (10). However, such boundary schemes were constructed with complex coefficients that are related not only to geometry-related Δ, but also to the LB model-related relaxation coefficient. Interested readers are referred to [32] for details and such interface schemes are not discussed in this chapter.

4.2 Modified geometry-ignored interface schemes

The second-order interpolation based interface scheme developed in [23] has attracted much interest. In the last 5 years, there have been various modified interface schemes proposed with the objective of simplifying the original scheme, as it becomes complex and computationally expensive when applied to curved or irregular interfaces. The applicability of those modified schemes was demonstrated in those publications while their accuracy and convergence order have not been fully investigated. In this section, we present three groups of those modified schemes that do not account for the local interface geometry. Most of those schemes were formulated for conjugate heat transfer problems, and conjugate mass transfer can be similarly handled. Comparison of their numerical accuracy with the original interpolation based scheme will be presented in Section 5.1.

4.2.1 Group 1: sourcing term addition

In the first group, additional source terms were introduced to the lattice nodes in the domains next to the interface. For example, the following additional source was given in [34]:

$$S_{\mathit{\text{conj}}}=\frac{\partial }{\partial x_j}\left(\frac{1}{\rho c_p}\right)·\left(-k\frac{\partial \varphi }{\partial x_j}+\rho c_p{u}_j\varphi \right)$$E15

In LBM, the total flux in the second bracket can be conveniently obtained from the moment of the nonequilibrium DFs [20, 21, 34]. In [34], the gradient of the heat capacity-related term in Eq. (15) was computed from a first-order one-sided finite-difference (FD) scheme: $\frac{\partial }{\partial x_j}{\left(\frac{1}{\rho c_p}\right)}_k=\left[{\left(\frac{1}{\rho c_p}\right)}_k-{\left(\frac{1}{\rho c_p}\right)}_{\mathrm{avg}}\right]/0.5\mathit{\delta x}$ with ${\left(\rho c_p\right)}_{\mathrm{avg}}=\left[{\left(\rho c_p\right)}_k+{\left(\rho c_p\right)}_{k+1}\right]/2$.

It should be noted that the above introduction of the source term and the calculation of the heat capacity-related gradient cause two issues: first, for adjacent domains with distinct $\rho c_p$ values, such as fluid-solid interfaces, a discontinuity shows up and the gradient term cannot be resolved with FD schemes; second, with the simple first-order FD approximation, the LBM solution would preserve, at most, up to first-order accuracy. The first-order accuracy was demonstrated in [34], and it will be further discussed in Section 5 with a numerical test.

4.2.2 Group 2: enthalpy-based formulation

Another area of interest is found in enthalpic formulations for the LBE. With the definition of an “enthalpic term” $h^{\ast }={\left(\rho c_p\right)}_0\varphi$, where ${\left(\rho c_p\right)}_0$ is a reference heat capacity, the governing CDE (2) becomes [27, 28]:

$$\frac{\partial h^{\ast }}{\partial t}+\nabla ·\left(\mathbf{u}{h}^{\ast }\right)=\nabla ·\left(\frac{k}{\rho c_p}\nabla h^{\ast }\right)-\frac{k}{{\left(\rho c_p\right)}_0}\nabla h^{\ast }·\nabla \left(\frac{{\left(\rho c_p\right)}_0}{\rho c_p}\right)-\frac{h^{\ast }{\left(\rho c_p\right)}_0}{\rho c_p}\mathbf{u}·\nabla \frac{\rho c_p}{{\left(\rho c_p\right)}_0}$$E16

Comparing Eq. (16) with Eq. (2), the last two terms in Eq. (16) need to be included in the source term implementation in LBM simulations. Clearly, those additional terms also have the heat capacity-related gradients and thus the discontinuity effect. In [36], the gradient was approximated from ${\nabla }_{j}f=\frac{1}{{c_s}^2\mathit{\delta t}}\sum _{\alpha }{w}_{\alpha }f\left(\mathbf{x}+\left({\mathbf{e}}_{\alpha }\mathbf{·}\mathbf{x}\right)\mathit{\delta t}\right)e_{\mathit{\alpha j}}$, which reduces to a central FD schemes in the Cartesian grid.

4.2.3 Group 3: modified equilibrium distribution functions

In this group, modified equilibrium DFs were introduced. The heat capacity is typically involved in the modified equilibrium DFs, such as [37, 38]

$$g_{\alpha }^{\mathrm{eq}}=\left\{\begin{array}{l}\varphi \left(\eta c_p-c_{p0}\right)+{\omega }_{\alpha }\varphi c_p\left[c_{p0}/c_p+\left({\mathbf{e}}_{\alpha }\cdot \mathbf{u}\right)/c_s^2\right],\alpha \ne 0\\ {\omega }_{\alpha }\varphi c_p\left[c_{p0}/c_p+\left({\mathbf{e}}_{\alpha }\cdot \mathbf{u}\right)/c_s^2\right],\alpha =0\end{array}\right.$$E17

When using the MRT D2Q5 model the equilibrium moments are calculated to be

$${\mathbf{m}}^{\mathrm{eq}}={\left[c_p\varphi ,{\mathit{uc}}_p\varphi ,{\mathit{vc}}_p\varphi ,4c_p\varphi -\frac{10}{3}{c}_{p0}\varphi ,0\right]}^T.$$E18

The temperature is solved from $\varphi =\sum _{\alpha }{g}_{\alpha }/\eta c_p$. A key note should be made regarding the relationship between transport properties and the relaxation coefficient, τ, in LBM for Group 3. In all previous models, the thermal diffusivity is related to τ as $D=\left(\tau -0.5\right)c_s^2\mathit{\delta t}$. However, since a virtual heat capacity correction is employed in this modified equilibrium DF group, the thermal conductivity, rather than the diffusivity, is related to τ as $k=c_{p0}\left(\tau -0.5\right)c_s^2\mathit{\delta t}$.

5. Numerical accuracy of interface schemes

In order to verify the applicability of the different interface schemes to simulate conjugate heat and mass transfer and compare their accuracy, we consider two benchmark cases with analytical solutions: (i) 2D convection-diffusion in a channel with two-layered fluids, and (ii) 2D diffusion within a circular domain of two solids with interfacial jump conditions. The computational domains are depicted in Figures 2 and 9, respectively. Those two cases have been widely employed in [9, 23, 32, 33, 34] to verify the various interface schemes.

Figure 2.

Schematic layout of the lattice on a 2D channel containing two fluids in domain 1 (0 ≤ y ≤ h) and domain 2 (h ≤ y ≤ H). With permission from [23].

For straight interfaces, the following relative L2 norm errors are defined to check the numerical accuracy and convergence orders following [23]:

$$E_2={\left[\sum _{x,y}{\left({\varphi }_{\mathrm{LBE}}-{\varphi }_{\mathrm{ex}}\right)}^2/\sum _{x,y}{\varphi }_{\mathrm{ex}}^2\right]}^{1/2}$$E19

$$E_{2,\varphi int}=\left[\sum _{x,y=h}{\left({\varphi }_{1,2}\left|{}_{\mathrm{LBE}}\right)-{\varphi }_{1,2}\left|{}_{\mathrm{ex}}\right)\right)}^2/\right.{\left.\sum _{x,y=h}{\left({\varphi }_{1,2}\left|{}_{\mathrm{ex}}\right)\right)}^2\right]}^{1/2}$$E20

$$E_{2,\text{qint}}={\left[\sum _{x,y=h}{\left(D_{1,2}\frac{\partial {\varphi }_{1,2}}{\partial y}\left|{}_{\mathrm{LBE}}\right)-D_{1,2}\frac{\partial {\varphi }_{1,2}}{\partial y}\left|{}_{\mathrm{ex}}\right)\right)}^2/\sum _{x,y=h}{\left(D_{1,2}\frac{\partial {\varphi }_{1,2}}{\partial y}\left|{}_{\mathrm{ex}}\right)\right)}^2\right]}^{1/2}$$E21

where E₂ evaluates the overall error in all the interior lattice nodes in the two domains, E_{2, ϕint} and E_{2, qint} evaluate the respective relative errors of the interfacial ϕ value and its flux. For circular interfaces, E_{2, ϕint} and E_{2, qint} are evaluated at the interface nodes along the curved geometry. The computation of the interfacial quantities follows that in [21, 23].

5.1 Two-D convection-diffusion in a channel with two fluids

The computational domain is depicted in Figure 2. The two fluids are assumed immiscible and both have the same velocity u = (U, 0). The characteristic Péclet number, Pe, is defined as Pe = UH/D₁.

When considering isotropic diffusion, the governing CDE can be expressed as

$$\frac{\partial {\varphi }_{1,2}}{\partial t}+U\frac{\partial {\varphi }_{1,2}}{\partial x}=D_{1,2}\left(\frac{{\partial }^2{\varphi }_{1,2}}{\partial x^2}+\frac{{\partial }^2{\varphi }_{1,2}}{\partial y^2}\right)$$E22

We consider only the steady case with sinusoidal boundary conditions on the horizontal walls ${\varphi }_1\left(x,y=0\right)={\varphi }_2\left(x,y=H\right)=cos\left(2\mathit{\pi x}/L\right)$, and periodic conditions in the x-direction. Taking into account the standard conjugate conditions: ϕ₁ = ϕ₂, k₁∂ϕ₁/∂y = k₂∂ϕ₂/∂y at y = h, the analytical solution to Eq. (22) can be solved (see [23] for details).

The interpolation-based interface scheme in Section 4.1 and the modified schemes in Groups 1–3 in Section 4.2 are implemented with the D2Q5 MRT-LB model. For all cases presented, Pe = 20 is used with H = 2 h.

Figures 3a and b show the scalar value and flux profiles at selected locations along the x-direction using a diffusivity ratio D₂₁ = D₂/D₁ = 2 and a thermal conductivity ratio k₂₁ = k₂/k₁ = 3 (consequently $\sigma ={\left(\rho c_p\right)}_2/{\left(\rho c_p\right)}_1$ = 1.5). Note that Δ = 0.5 is used for all schemes and thus the geometry effect is not included. Other simulation parameters are: τ₁ = 0.65, H = 64. Good agreement between LBM and analytical solutions is observed in Figures 3a and b for all interface schemes. As a further step, Figures 4a and b compares the interfacial scalar and flux values using the same parameters as for Figure 3. Noticeable discrepancy between numerical and analytical solutions is observed in Figure 4a and b when using the interface schemes of Group 1 and Group 2; while both Group 3 and the original interface scheme in [23] have good agreement. This is mainly due to the presence of the discontinuity in approximating the heat capacity gradient in Groups 1 and 2 when $\sigma ={\left(\rho c_p\right)}_2/{\left(\rho c_p\right)}_1\ne 1$. This will be further demonstrated in the L₂ norm errors.

Figure 3.

Profiles of the (a) scalar variable, and (b) flux values at selected vertical lines in the channel.

Figure 4.

Comparison of (a) interfacial scalar value, and (b) interfacial fluxes.

Figures 5–7 show the L₂ errors defined in Eqs. (19)–(21), respectively. The physical and geometric parameters are D₂₁ = 10, k₂₁ = 100, σ = 10, and Δ = 0.5. Simulation parameters in LBM include τ₁ = 0.55, τ₂ = 1.0 for the original scheme [23] and Groups 1 and 2 with (τ₂₋0.5)/(τ₁₋0.5) = D₂₁, and τ₁ = 0.5025, τ₂ = 0.75 for Group 3 with (τ₂₋0.5)/(τ₁₋0.5) = k₂₁.

Figure 5.

Relative L2 norm error, E2, for the interior scalar versus the grid resolution, 1/H, for steady convection-diffusion in the channel at Δ = 0.5.

Figure 6.

Relative L2 norm error, E2, ϕint, for the interfacial scalar versus the grid resolution, 1/H, for steady convection-diffusion in the channel at Δ = 0.5.

Figure 7.

Relative L2 norm error, E2, qint, for the interfacial flux versus the grid resolution, 1/H, for steady convection-diffusion in the channel at Δ = 0.5.

Figure 5 clearly shows that the original scheme in [23] and Group 3 are able to preserve the second-order accuracy in LBM. However, Groups 1 and 2 show only a linear convergence at low resolution, and it reduces further towards zeroth-order convergence at high resolution. Similar observations can be found in the errors for the interfacial scalar values in Figure 6. For the interfacial flux errors shown in Figure 7, Group 2 always exhibits zeroth-order accuracy, while Group 1 exhibits linear convergence in one domain and zeroth-order in the other. As previously mentioned, a discontinuity approximation is present in the development of interface schemes in Groups 1 and 2. This is the direct cause for the degradation of the order of accuracy and the much higher error magnitude in these two groups.

To further demonstrate the necessity of taking into account the local geometry to preserve second-order accuracy, Figure 8 presents the L₂ norm errors for the interior field at different Δ values. The parameters used are D₂₁ = 5, k₂₁ = 50, σ = 10, (τ₁, τ₂) = (0.515, 0.575) for the interpolation-based scheme [23] and Groups 1 and 2, and (τ₁, τ₂) = (0.515, 1.25) for Group 3. For the interpolation-based scheme with an adjustable variable, c_d1 = −1 is used.

Figure 8.

Relative L2 norm error, E2, for the interior scalar versus the grid resolution, 1/H, for steady convection-diffusion in the channel at different Δ values.

It is clear in Figure 8 that only the interpolation-based scheme considering the local geometry is able to preserve second-order accuracy at different Δ values. The error behavior for Groups 1 and 2 is similar to that in Figure 5, the near identical errors at high resolution for different Δ values confirm the error caused by the same discontinuity; moreover, while Group 3 presents second-order convergence for Δ = 0.5, the convergence order drops to first-order for Δ ≠ 0.5. This is similar to the behavior of the “half-lattice division” scheme discussed in [23]. It is evident that when the local interface geometry is not considered for cases Δ ≠ 0.5, the inherent second-order accuracy in LBM computation can be lost.

5.2 Two-D diffusion in a circular domain with jump conditions

This test is applied to demonstrate the applicability and accuracy of the interpolation-based interface scheme to simulate transport in complex geometry and with interfacial jump conditions. Figure 9 schematically presents the lattice layout and computational domain for 2D diffusion within a circular domain. On the outer circle with radius, R₂, a Dirichlet boundary condition is applied as ϕ₂(r = R₂) = cosφ. With the following conjugate conditions at the interface:

Figure 9.

Schematic layout of the computational domain for circular diffusion. With permission from [9].

$${\varphi }_1={\varphi }_2+{\varphi }_{\mathit{\text{jump}}}={\varphi }_2+{\varphi }_{0}cos\phi \mathrm{at}r=R_1$$E23

and

$$-D_1\frac{{\partial \varphi }_1}{\partial r}=-D_2\frac{{\partial \varphi }_2}{\partial r}+q_{\mathit{\text{jump}}}=-D_2\frac{{\partial \varphi }_2}{\partial r}+q_{0}cos\phi \mathrm{at}r=R_1$$E24

the analytical solution can be found [9]. For the numerical LBM computation, the parameters are R₂/R₁ = 2, D₂/D₁ = 10, ${\tau }_1=0.525$, and ${\tau }_2=0.75$. As previously mentioned, three schemes were presented in [21] for adjustable parameters of c_d1, those three schemes in Eqs. (12a)–(12c) are also applied here. For the flux jump case, the jump conditions are set as ϕ₀ = 0 and q₀ = D₂/R₁; and for the temperature jump case, those are ϕ₀ = 0.5 and q₀ = 0.

Figures 10 and 11a, b provide the respective relative L₂ norm errors for the interior temperature and the interfacial temperature and flux with the conjugate schemes implemented.

Figure 10.

Relative L2 norm error, E2, for the interior temperature versus the grid resolution for 2D diffusion in a circular domain. With permission from [9].

Figure 11.

Relative L2 norm errors (a) E2_tint for the interfacial temperature, and (b) E2_qint for the interfacial flux, versus the grid resolution for 2D diffusion in a circular domain. With permission from [9].

The results in Figures 10 and 11 demonstrate the first-order accuracy at high resolution for all cases. It should also be noted that the temperature and flux jump conditions at the interface do not affect the order of convergence. The decrease of the convergence order from second-order in Section 5.1 to first-order in this test is due to the implementation of the Cartesian decomposition method [21] that was used to convert the normal fluxes into those in the discrete velocity directions. It is expected that the modified geometry-ignored interface schemes in Groups 1–3 would result in much higher error magnitude for curved interfaces, and Groups 1 and 2 would yield only zeroth-order convergence for all the three quantities of interest. This will be presented in future publications.

6. Conclusions

The chapter presents a brief review of the interface schemes within the scope of the lattice Boltzmann method (LBM) for conjugate transport between multiphases or different materials. Compared to the interface schemes developed to satisfy the macroscopic conjugate conditions using traditional CFD methods, the LBM method deals with the microscopic distribution functions (DFs); the physical conjugate conditions can be converted to those for the DFs, and they are satisfied in each time-marching step without iterations. In the last decade, a number of interface schemes have been proposed. The interpolation-based schemes [9, 23, 33] taking into account the local interfacial geometry are able to preserve the second-order accuracy in LBM for straight interfaces; while those “modified” geometry-ignored schemes [34, 35, 36, 37, 38] have at most first-order accuracy in general, and with the introduction of heat capacity-related discontinuity in those schemes (e.g., Groups 1 and 2), the order of accuracy becomes essentially zeroth order.

Furthermore, it is verified that when using the interpolation-based schemes, the interfacial jump conditions can be conveniently modeled with no effect on the order of accuracy of the LBM solutions.

Curved geometry also has a substantial effect and it reduces the order of accuracy of LBM solutions in modeling conjugate heat and mass transfer problems. In addition, the interpolation-based schemes would demand for a higher computational cost than those modified schemes. The readers are thus recommended to take into account both numerical accuracy and computational cost when selecting effective interface schemes for curved geometries.

Acknowledgments

LL acknowledges the support and start-up fund from the Department of Mechanical Engineering at Mississippi State University.