Open access peer-reviewed chapter

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

Written By

David Korba and Like Li

Submitted: March 31st, 2019 Reviewed: April 9th, 2019 Published: May 14th, 2019

DOI: 10.5772/intechopen.86252

Chapter metrics overview

770 Chapter Downloads

View Full Metrics


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.


  • conjugate conditions
  • boundary conditions
  • heat and mass transfer
  • lattice Boltzmann
  • numerical accuracy

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].

n·kϕ+ρcpuϕf=n·kϕ+ρcpuϕs+n·qjumpTin heat transfer,E1b
n·Dmϕ+uϕf=n·Dmϕ+uϕs+n·qjumpCin mass transfer,E1c

where n represents the normal direction, ϕ the macroscopic scalar variable of interest (temperature or concentration), k the thermal conductivity, ρ the density, cp the heat capacity, u the velocity vector, Dm the mass diffusivity, and ϕjump and qjump 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:


where Dij 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


where the microscopic distribution function, gα(x, t 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αeqxt is the equilibrium distribution function, and ωα is the weight coefficient. The macroscopic scalar variable is obtained from


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


where cs is the speed of sound with cs=c/3=δx/δt/3=1/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


Streaming step


In the above, M is a matrix to transform the distribution functions geq to their moments meq by the relation meq=Mgeq. S is a matrix of relaxation coefficients τij. In the D3Q7 model, the equilibrium moments of the distribution functions are meq=ϕ00T,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 (un = 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)


with σ = 1, qjump=qjumpC in mass transfer, and σ=ρcpsρcpf, qjump=qjumpT/ρcpf 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,


Similarly, for the Neumann condition treatment,


where Φnα¯ and Φ are the respective interfacial fluxes along the discrete lattice velocity directions eα¯ and eα, xf and xff are the first and second interior lattice nodes along eα¯ direction in Domain 1 (xff = xf + eα¯ δt), and xs and xss are the lattice nodes along eα direction in Domain 2 (xss = xs + eα δt = xseα¯ δ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 eα¯ is aligned with the interface normal directions, Φnα¯=Φnf and Φ=Φ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α¯xft+δt, gαxst+δt, Φdf, Φds, Φnf, and Φns can be analytically solved. The interface scheme thus becomes [9]:


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:


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αxst+δt=1σ1+σĝα¯xst+21+σĝαxft+qjump1+σεDϕjump1+σ. E13b

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


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]:


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: xj1ρcpk=1ρcpk1ρcpavg/0.5δx with ρcpavg=ρcpk+ρcpk+1/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 ρcp 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=ρcp0ϕ, where ρcp0 is a reference heat capacity, the governing CDE (2) becomes [27, 28]:


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 jf=1cs2δtαwαfx+eα·xδteα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]


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


The temperature is solved from ϕ=αgα/ηcp. 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=τ0.5cs2δ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=cp0τ0.5cs2δ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]:


where E2 evaluates the overall error in all the interior lattice nodes in the two domains, E2, ϕint and E2, qint evaluate the respective relative errors of the interfacial ϕ value and its flux. For circular interfaces, E2, ϕint and E2, 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/D1.

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


We consider only the steady case with sinusoidal boundary conditions on the horizontal walls ϕ1xy=0=ϕ2xy=H=cos2πx/L, and periodic conditions in the x-direction. Taking into account the standard conjugate conditions: ϕ1 = ϕ2, k1∂ϕ1/∂y = k2∂ϕ2/∂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 D21 = D2/D1 = 2 and a thermal conductivity ratio k21 = k2/k1 = 3 (consequently σ=ρcp2/ρcp1 = 1.5). Note that Δ = 0.5 is used for all schemes and thus the geometry effect is not included. Other simulation parameters are: τ1 = 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 σ=ρcp2/ρcp11. This will be further demonstrated in the L2 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 57 show the L2 errors defined in Eqs. (19)(21), respectively. The physical and geometric parameters are D21 = 10, k21 = 100, σ = 10, and Δ = 0.5. Simulation parameters in LBM include τ1 = 0.55, τ2 = 1.0 for the original scheme [23] and Groups 1 and 2 with (τ2−0.5)/(τ1−0.5) = D21, and τ1 = 0.5025, τ2 = 0.75 for Group 3 with (τ2−0.5)/(τ1−0.5) = k21.

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 L2 norm errors for the interior field at different Δ values. The parameters used are D21 = 5, k21 = 50, σ = 10, (τ1, τ2) = (0.515, 0.575) for the interpolation-based scheme [23] and Groups 1 and 2, and (τ1, τ2) = (0.515, 1.25) for Group 3. For the interpolation-based scheme with an adjustable variable, cd1 = −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, R2, a Dirichlet boundary condition is applied as ϕ2(r = R2) = cosφ. With the following conjugate conditions at the interface:

Figure 9.

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




the analytical solution can be found [9]. For the numerical LBM computation, the parameters are R2/R1 = 2, D2/D1 = 10, τ1=0.525, and τ2=0.75. As previously mentioned, three schemes were presented in [21] for adjustable parameters of cd1, those three schemes in Eqs. (12a)(12c) are also applied here. For the flux jump case, the jump conditions are set as ϕ0 = 0 and q0 = D2/R1; and for the temperature jump case, those are ϕ0 = 0.5 and q0 = 0.

Figures 10 and 11a, b provide the respective relative L2 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.



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


  1. 1. Roe B, Jaiman R, Haselbacher A, Geubelle PH. Combined interface boundary condition method for coupled thermal simulations. International Journal for Numerical Methods in Fluids. 2008;57:329
  2. 2. Dorfman A, Renner Z. Conjugate problems in convective heat transfer: Review. Mathematical Problems in Engineering. 2009;2009:927350
  3. 3. Errera MP, Chemin S. Optimal solutions of numerical interface conditions in fluid-structure thermal analysis. Journal of Computational Physics. 2013;245:431
  4. 4. Li L, Mei R, Klausner JF, Hahn DW. Heat transfer between colliding surfaces and particles. ASME Journal of Heat Transfer. 2012;134:011301
  5. 5. Zhang J, Yang C, Mao ZS. Unsteady conjugate mass transfer from a spherical drop in a simple extensional creeping flow. Chemical Engineering Science. 2012;79:29
  6. 6. Wang J, Wang M, Li Z. A lattice Boltzmann algorithm for fluid-solid conjugate heat transfer. International Journal of Thermal Sciences. 2007;46:228-234
  7. 7. Wang M, Wang J, Pan N, Chen S. Mesoscopic predictions of the effective thermal conductivity for microscale random porous media. Physical Review E. 2007;75:036702
  8. 8. Chen L, Kang Q, He YL, Tao WQ. Pore-scale simulation of coupled multiple physicochemical thermal processes in micro reactor for hydrogen production using lattice Boltzmann method. International Journal of Hydrogen Energy. 2012;37:13943
  9. 9. Guo K, Li L, Xiao G, AuYeung N, Mei R. Lattice Boltzmann method for conjugate heat and mass transfer with interfacial jump conditions. International Journal of Heat and Mass Transfer. 2015;88:306-322
  10. 10. Lu JH, Lei HY, Dai CS. Analysis of Henry’s law and a unified lattice Boltzmann equation for conjugate mass transfer problem. Chemical Engineering Science. 2019;199:319-331
  11. 11. Peskin CS. The immersed boundary method. Acta Numerica. 2002;11:479-517
  12. 12. Mittal R, Iaccarino G. Immersed boundary methods. Annual Review of Fluid Mechanics. 2005;37:239-261
  13. 13. LeVeque RJ, Li ZL. The immersed interface method for elliptic equations with discontinuous coefficients and singular sources. SIAM Journal on Numerical Analysis. 1994;31:1019-1044
  14. 14. Deng SZ, Ito K, Li ZL. Three-dimensional elliptic solvers for interface problems and applications. Journal of Computational Physics. 2003;184:215-243
  15. 15. Fedkiw RP, Aslam T, Merriman B, Osher S. A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method). Journal of Computational Physics. 1999;152:457-492
  16. 16. Liu XD, Fedkiw RP, Kang M. A boundary condition capturing method for Poisson’s equation on irregular domains. Journal of Computational Physics. 2000;160:151-178
  17. 17. Udaykumar HS, Mittal R, Rampunggoon P, Khanna A. A sharp interface Cartesian grid method for simulating flows with complex moving boundaries. Journal of Computational Physics. 2001;174:345-380
  18. 18. Marella S, Krishnan S, Liu H, Udaykumar HS. Sharp interface Cartesian grid method I: An easily implemented technique for 3D moving boundary computations. Journal of Computational Physics. 2005;210:1-31
  19. 19. Zhou YC, Zhao S, Feig M, Wei GW. High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources. Journal of Computational Physics. 2006;213:1-30
  20. 20. Yoshida H, Nagaoka M. Multiple-relaxation-time lattice Boltzmann model for the convection and anisotropic diffusion equation. Journal of Computational Physics. 2010;229:7774
  21. 21. Li L, Mei R, Klausner JF. Boundary conditions for thermal lattice Boltzmann equation method. Journal of Computational Physics. 2013;237:366
  22. 22. Li L, Mei R, Klausner JF. Heat transfer evaluation on curved boundaries in thermal lattice Boltzmann equation method. ASME Journal of Heat Transfer. 2014;136:012403
  23. 23. Li L, Chen C, Mei R, Klausner JF. Conjugate heat and mass transfer in the lattice Boltzmann equation method. Physical Review E. 2014;89:043308
  24. 24. Krivovichev GV. Linear Bhatnagar-Gross-Krook equations for simulation of linear diffusion equation by lattice Boltzmann method. Applied Mathematics and Computation. 2018;325:102-119
  25. 25. Perko J. A single-relaxation-time lattice Boltzmann model for anisotropic advection-diffusion equation based on the diffusion velocity flux formulation. Computational Geosciences. 2018;22:1423-1432
  26. 26. Zhang XX, Bengough AG, Crawford JW, Young IM. A lattice BGK model for advection and anisotropic dispersion equation. Advances in Water Resources. 2002;25:1-8
  27. 27. Ginzburg I. Equilibrium-type and link-type lattice Boltzmann models for generic advection and anisotropic-dispersion equation. Advances in Water Resources. 2005;28:1171-1195
  28. 28. Ginzburg I, d’Humières D, Kuzmin A. Optimal stability of advection-diffusion lattice Boltzmann models with two relaxation times for positive/negative equilibrium. Journal of Statistical Physics. 2010;139:1090-1143
  29. 29. Chai Z, Shi B, Guo Z. A multiple-relaxation-time lattice Boltzmann model for general nonlinear anisotropic convection-diffusion equations. Journal of Scientific Computing. 2016;69:355-390
  30. 30. Chai Z, Zhao TS. Lattice Boltzmann model for the convection-diffusion equation. Physical Review E. 2013;87:063309
  31. 31. Shi B, Guo Z. Lattice Boltzmann model for nonlinear convection-diffusion equations. Physical Review E. 2009;79:016701
  32. 32. Hu Z, Huang J, Yong WA. Lattice Boltzmann method for convection-diffusion equations with general interfacial conditions. Physical Review E. 2016;93:043320
  33. 33. Mu Y-T, Gu Z-L, He P, Tao W-Q. Lattice Boltzmann method for conjugated heat and mass transfer with general interfacial conditions. Physical Review E. 2018;98:043309
  34. 34. Karani H, Huber C. Lattice Boltzmann formulation for conjugate heat transfer in heterogeneous media. Physical Review E. 2015;91:023304
  35. 35. Rihab H, Moudhaffar N, Sassi BN, Patrick P. Enthalpic lattice Boltzmann formulation for unsteady heat conduction in heterogeneous media. International Journal of Heat and Mass Transfer. 2016;100:728
  36. 36. Chen S, Yan YY, Gong W. A simple lattice Boltzmann model for conjugate heat transfer research. International Journal of Heat and Mass Transfer. 2017;107:862
  37. 37. Gao D, Chen Z, Chen L, Zhang D. A modified lattice Boltzmann model for conjugate heat transfer in porous media. International Journal of Heat and Mass Transfer. 2017;105:673-683
  38. 38. Chen S. Simulation of conjugate heat transfer between fluid-saturated porous media and solid wall. International Journal of Thermal Sciences. 2018;124:477-483
  39. 39. Yu D, Mei R, Luo L-S, Shyy W. Viscous flow computations with the method of lattice Boltzmann equation. Progress in Aerospace Science. 2003;39:329-367
  40. 40. Li L, Mei R, Klausner JF. Lattice Boltzmann models for the convection-diffusion equation: D2Q5 vs. D2Q9. International Journal of Heat and Mass Transfer. 2017;108:41-62

Written By

David Korba and Like Li

Submitted: March 31st, 2019 Reviewed: April 9th, 2019 Published: May 14th, 2019