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
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 , and Henry’s law relationship , 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 . 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 . 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  for the nonlinear CDE, the D3Q7/D2Q5 MRT models by Yoshida and Nagaoka  for the general convection anisotropic diffusion equation, and the D1Q3/D2Q9/D3Q19 MRT models by Chai and Zhao  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  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.  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. . 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. ) for unsteady cases. Importantly, with the interface treated as a shared boundary between the adjacent domains, the boundary conditions by Li et al.  were applied to interface conditions and particular interface schemes were proposed and verified in  for standard conjugate conditions, and in  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 , 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 , 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:
Some examples of jump conditions can be found in cases such as concentration jumps (Henry’s law ) or temperature jumps at the interface . Eqs. (1a)–(1c) reduce to the standard conjugate conditions in  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 represents the diffusion coefficient, and 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 . 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,
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,
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 (
with σ = 1, in mass transfer, and , 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 , the relationships between the distribution functions and the interfacial
Similarly, for the Neumann condition treatment,
where and are the respective interfacial fluxes along the discrete lattice velocity directions and ,
When is aligned with the interface normal directions, and are readily noticed. Hence, Eqs. (1), (8), (9a), (9b), (10a) and (10b) constitute a linear system of six equations, and the six unknowns , , , , , and can be analytically solved. The interface scheme thus becomes :
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 , 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:
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 , 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  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  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 :
In LBM, the total flux in the second bracket can be conveniently obtained from the moment of the nonequilibrium DFs [20, 21, 34]. In , the gradient of the heat capacity-related term in Eq. (15) was computed from a first-order one-sided finite-difference (FD) scheme: with .
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 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 , and it will be further discussed in Section 5 with a numerical test.
4.2.2 Group 2: enthalpy-based formulation
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 , the gradient was approximated from , which reduces to a central FD schemes in the Cartesian grid.
4.2.3 Group 3: modified equilibrium distribution functions
When using the MRT D2Q5 model the equilibrium moments are calculated to be
The temperature is solved from . A key note should be made regarding the relationship between transport properties and the relaxation coefficient,
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.
For straight interfaces, the following relative L2 norm errors are defined to check the numerical accuracy and convergence orders following :
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
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 , and periodic conditions in the
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,
Figures 3a and b show the scalar value and flux profiles at selected locations along the
Figures 5–7 show the
Figure 5 clearly shows that the original scheme in  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
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 . 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,
the analytical solution can be found . For the numerical LBM computation, the parameters are
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  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.
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.