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

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

$$\frac{\partial \phi}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_j} \left( \nu\_j \phi \right) = \frac{\partial}{\partial \mathbf{x}\_j} \left( D\_{i\bar{j}} \frac{\partial \phi}{\partial \mathbf{x}\_j} \right) + G \tag{2}$$

Streaming step

*DOI: http://dx.doi.org/10.5772/intechopen.86252*

refer to [20, 32].

to (see **Figure 1**)

heat transfer.

condition treatment,

*<sup>g</sup><sup>α</sup>* **<sup>x</sup>***<sup>s</sup>* ð Þ¼ *; <sup>t</sup>* <sup>þ</sup> *<sup>δ</sup><sup>t</sup> <sup>c</sup>* <sup>∗</sup>

**69**

*<sup>g</sup><sup>α</sup>* **<sup>x</sup>***<sup>s</sup>* ð Þ¼ *; <sup>t</sup>* <sup>þ</sup> *<sup>δ</sup><sup>t</sup> <sup>c</sup>* <sup>∗</sup>

with <sup>σ</sup> = 1, *<sup>q</sup>*jump <sup>¼</sup> *<sup>q</sup>*<sup>C</sup>

*gα*ð**x** þ **e***αδt; t* þ *δt*Þ ¼ ^*gα*ð Þ **x***; t :* (7)

In the above, **M** is a matrix to transform the distribution functions **g**ð Þ eq to their moments **<sup>m</sup>**ð Þ eq by the relation **<sup>m</sup>**ð Þ eq <sup>¼</sup> **Mg**ð Þ eq . **<sup>S</sup>** is a matrix of relaxation coefficients *τij*. In the D3Q7 model, the equilibrium moments of the distribution functions are **<sup>m</sup>**ð Þ eq <sup>¼</sup> ð Þ *<sup>ϕ</sup>; <sup>u</sup>ϕ; <sup>v</sup>ϕ; <sup>w</sup>ϕ; <sup>a</sup>ϕ;* <sup>0</sup>*;* <sup>0</sup> *T,* where *<sup>u</sup>*, *<sup>v</sup>*, and *<sup>w</sup>* 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

*Interface Treatment for Conjugate Conditions in the Lattice Boltzmann Method…*

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 finitedifference 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

**4. Interface schemes for conjugate conditions**

Φ*nf* ¼ �*Df*

**4.1 Interpolation-based interface schemes**

*∂ϕf ∂nf*

*<sup>d</sup>*<sup>1</sup> ^*g<sup>α</sup>* **<sup>x</sup>***<sup>s</sup>* ð Þþ *; <sup>t</sup> <sup>c</sup>* <sup>∗</sup>

Similarly, for the Neumann condition treatment,

*<sup>n</sup>*<sup>1</sup> ^*g<sup>α</sup>* **<sup>x</sup>***<sup>s</sup>* ð Þþ *; <sup>t</sup> <sup>c</sup>* <sup>∗</sup>

¼ *σDs*

*∂ϕs ∂ns*

jump in mass transfer, and *<sup>σ</sup>* <sup>¼</sup> ð Þ *<sup>ρ</sup>cp <sup>s</sup>*

schemes fall into two different categories, as discussed in Sections 4.1 and 4.2.

Depending on whether the local interface geometry is considered, the interface

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

*<sup>g</sup><sup>α</sup>* **<sup>x</sup>***<sup>f</sup> ; <sup>t</sup>* <sup>þ</sup> *<sup>δ</sup><sup>t</sup>* <sup>¼</sup> *cd*<sup>1</sup> ^*g<sup>α</sup>* **<sup>x</sup>***<sup>f</sup> ; <sup>t</sup>* <sup>þ</sup> *cd*<sup>2</sup> ^*g<sup>α</sup>* **<sup>x</sup>***ff ; <sup>t</sup>* <sup>þ</sup> *cd*<sup>3</sup> ^*g<sup>α</sup>* **<sup>x</sup>***<sup>f</sup> ; <sup>t</sup>* <sup>þ</sup> *cd*4*εD*Φ*df ,* (9a)

*<sup>d</sup>*<sup>2</sup> ^*g<sup>α</sup>* **<sup>x</sup>***ss* ð Þþ *; <sup>t</sup> <sup>c</sup>* <sup>∗</sup>

*<sup>g</sup><sup>α</sup>* **<sup>x</sup>***<sup>f</sup> ; <sup>t</sup>* <sup>þ</sup> *<sup>δ</sup><sup>t</sup>* <sup>¼</sup> *cn*<sup>1</sup> ^*g<sup>α</sup>* **<sup>x</sup>***<sup>f</sup> ; <sup>t</sup>* <sup>þ</sup> *cn*<sup>2</sup> ^*g<sup>α</sup>* **<sup>x</sup>***ff ; <sup>t</sup>* <sup>þ</sup> *cn*<sup>3</sup> ^*g<sup>α</sup>* **<sup>x</sup>***<sup>f</sup> ; <sup>t</sup>* <sup>þ</sup> *cn*4ð Þ *<sup>δ</sup>t=δ<sup>x</sup>* <sup>Φ</sup>*n<sup>α</sup>,* (10a)

*<sup>n</sup>*<sup>2</sup> ^*g<sup>α</sup>* **<sup>x</sup>***ss* ð Þþ *; <sup>t</sup> <sup>c</sup>* <sup>∗</sup>

þ *q*jump ¼ �*σ*Φ*ns* þ *q*jump*,* (8)

, *<sup>q</sup>*jump <sup>¼</sup> *<sup>q</sup>*<sup>T</sup>

jump*= ρcp <sup>f</sup>* in

*<sup>d</sup>*4*εD*Φ*ds:* (9b)

*<sup>n</sup>*4ð Þ *δt=δx* Φ*n<sup>α</sup>,* (10b)

ð Þ *<sup>ρ</sup>cp <sup>f</sup>*

*<sup>d</sup>*<sup>3</sup> ^*g<sup>α</sup>* **<sup>x</sup>***<sup>s</sup>* ð Þþ *; <sup>t</sup> <sup>c</sup>* <sup>∗</sup>

*<sup>n</sup>*<sup>3</sup> ^*g<sup>α</sup>* **<sup>x</sup>***<sup>s</sup>* ð Þþ *; <sup>t</sup> <sup>c</sup>* <sup>∗</sup>

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

$$\mathbf{g}\_a(\mathbf{x} + \mathbf{e}\_a \delta t, t + \delta t) - \mathbf{g}\_a(\mathbf{x}, t) = \left[\mathbf{L} \cdot (\mathbf{g} - \mathbf{g}^{\text{eq}})(\mathbf{x}, t)\right]\_a + \alpha\_a G(\mathbf{x}, t) \delta t \tag{3}$$

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***<sup>α</sup>* is the *α*th discrete velocity vector, *δt* is the time step, **L** is the collision operator, *g*eq *<sup>α</sup>* ð Þ **x***; t* is the equilibrium distribution function, and *ωα* is the weight coefficient. The macroscopic scalar variable is obtained from

$$\phi(\mathbf{x},t) = \sum\_{a=0}^{m-1} \mathbf{g}\_a(\mathbf{x},t) \tag{4}$$

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

$$\mathbf{g}\_a^{\text{eq}} = a \boldsymbol{\alpha} \boldsymbol{\phi} \left( \mathbf{1} + \frac{\mathbf{e}\_a \cdot \mathbf{u}}{c\_s^2} \right) \tag{5}$$

where *cs* is the speed of sound with *cs* <sup>¼</sup> *<sup>c</sup><sup>=</sup>* ffiffiffi <sup>3</sup> <sup>p</sup> <sup>¼</sup> ð Þ *<sup>δ</sup>x=δ<sup>t</sup> <sup>=</sup>* ffiffiffi <sup>3</sup> <sup>p</sup> <sup>¼</sup> <sup>1</sup>*<sup>=</sup>* ffiffiffi <sup>3</sup> <sup>p</sup> .

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

$$\hat{\mathbf{g}}\_a(\mathbf{x},t) = \mathbf{g}\_a(\mathbf{x},t) - \left[\mathbf{M}^1\mathbf{S}\left(\mathbf{m}(\mathbf{x},t) - \mathbf{m}^{\text{(eq)}}(\mathbf{x},t)\right)\right]\_a + \omega\_a G(\mathbf{x},t)\delta t,\text{ and}\tag{6}$$

*Interface Treatment for Conjugate Conditions in the Lattice Boltzmann Method… DOI: http://dx.doi.org/10.5772/intechopen.86252*

Streaming step

$$
\mathcal{g}\_a(\mathbf{x} + \mathbf{e}\_a \delta t, t + \delta t) = \hat{\mathcal{g}}\_a(\mathbf{x}, t). \tag{7}
$$

In the above, **M** is a matrix to transform the distribution functions **g**ð Þ eq to their moments **<sup>m</sup>**ð Þ eq by the relation **<sup>m</sup>**ð Þ eq <sup>¼</sup> **Mg**ð Þ eq . **<sup>S</sup>** is a matrix of relaxation coefficients *τij*. In the D3Q7 model, the equilibrium moments of the distribution functions are **<sup>m</sup>**ð Þ eq <sup>¼</sup> ð Þ *<sup>ϕ</sup>; <sup>u</sup>ϕ; <sup>v</sup>ϕ; <sup>w</sup>ϕ; <sup>a</sup>ϕ;* <sup>0</sup>*;* <sup>0</sup> *T,* where *<sup>u</sup>*, *<sup>v</sup>*, and *<sup>w</sup>* 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 finitedifference 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**)

$$\Phi\_{\mathfrak{n}f} = -D\_f \frac{\partial \phi\_f}{\partial \mathfrak{n}\_f} = \sigma D\_s \frac{\partial \phi\_s}{\partial \mathfrak{n}\_s} + q\_{\text{jump}} = -\sigma \Phi\_{\mathfrak{n}s} + q\_{\text{jump}},\tag{8}$$

with <sup>σ</sup> = 1, *<sup>q</sup>*jump <sup>¼</sup> *<sup>q</sup>*<sup>C</sup> jump in mass transfer, and *<sup>σ</sup>* <sup>¼</sup> ð Þ *<sup>ρ</sup>cp <sup>s</sup>* ð Þ *<sup>ρ</sup>cp <sup>f</sup>* , *<sup>q</sup>*jump <sup>¼</sup> *<sup>q</sup>*<sup>T</sup> jump*= ρcp <sup>f</sup>* 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,

$$\mathcal{g}\_{\overline{a}}(\mathbf{x}\_{\mathcal{f}},t+\delta t) = c\_{d1}\hat{\mathbf{g}}\_{a}(\mathbf{x}\_{\mathcal{f}},t) + c\_{d2}\hat{\mathbf{g}}\_{a}(\mathbf{x}\_{\mathcal{f}},t) + c\_{d3}\hat{\mathbf{g}}\_{\overline{a}}(\mathbf{x}\_{\mathcal{f}},t) + c\_{d4}\varepsilon\_{D}\Phi\_{d\mathcal{f}},\tag{9a}$$

$$\mathbf{g}\_a(\mathbf{x}\_t, t + \delta t) = c\_{d1}^\* \hat{\mathbf{g}}\_{\overline{a}}(\mathbf{x}\_t, t) + c\_{d2}^\* \hat{\mathbf{g}}\_{\overline{a}}(\mathbf{x}\_t, t) + c\_{d3}^\* \hat{\mathbf{g}}\_a(\mathbf{x}\_t, t) + c\_{d4}^\* \varepsilon\_D \Phi\_{ds}.\tag{9b}$$

Similarly, for the Neumann condition treatment,

$$\mathbf{g\_{\overline{a}}}(\mathbf{x\_{\{f,t}}}t+\delta t) = c\_{n1}\hat{\mathbf{g\_{\underline{a}}}}(\mathbf{x\_{\{f,t}}}t) + c\_{n2}\hat{\mathbf{g\_{\underline{a}}}}(\mathbf{x\_{\{f,t}}}t) + c\_{n3}\hat{\mathbf{g\_{\underline{a}}}}(\mathbf{x\_{\{f,t}}}t) + c\_{n4}(\delta t/\delta \mathbf{x})\Phi\_{\overline{a}\overline{a}} \tag{10a}$$

$$\mathbf{g}\_a(\mathbf{x}\_t, t + \delta t) = c\_{n1}^\* \hat{\mathbf{g}}\_{\overline{a}}(\mathbf{x}\_t, t) + c\_{n2}^\* \hat{\mathbf{g}}\_{\overline{a}}(\mathbf{x}\_{tt}, t) + c\_{n3}^\* \hat{\mathbf{g}}\_a(\mathbf{x}\_t, t) + c\_{n4}^\*(\delta t / \delta \mathbf{x}) \Phi\_{n\alpha} \tag{10b}$$

where Φ*n<sup>α</sup>* and Φ*n<sup>α</sup>* are the respective interfacial fluxes along the discrete lattice velocity directions **e** *<sup>α</sup>* and **e***α*, **x***<sup>f</sup>* and **x***ff* are the first and second interior lattice nodes along **e** *<sup>α</sup>* direction in Domain 1 (**x***ff* = **x***<sup>f</sup>* + **e** *<sup>α</sup> δt*), and **x***<sup>s</sup>* and **x***ss* are the lattice nodes along **e***<sup>α</sup>* direction in Domain 2 (**x***ss* = **x***<sup>s</sup>* + **e***<sup>α</sup> δt* = **x***<sup>s</sup>* � **e** *<sup>α</sup> δ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].

*<sup>g</sup><sup>α</sup>* **<sup>x</sup>***<sup>s</sup>* ð Þ¼ *; <sup>t</sup>* <sup>þ</sup> *<sup>δ</sup><sup>t</sup>* ^*g<sup>α</sup>* **<sup>x</sup>***<sup>f</sup> ; <sup>t</sup>* � � <sup>þ</sup> *<sup>q</sup>*jump*=*<sup>2</sup> � *<sup>ε</sup>Dϕ*jump*=*2*:* (14b)

For straight interfaces where Δ 6¼ 0.5 and for curved interfaces, the complete

*Interface Treatment for Conjugate Conditions in the Lattice Boltzmann Method…*

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

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

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

> � �*<sup>k</sup> <sup>∂</sup><sup>ϕ</sup> ∂xj*

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-

þ *ρcpujϕ*

� � (15)

0*:*5*δx* with

1 *ρcp* � �

*<sup>k</sup>* <sup>¼</sup> <sup>1</sup> *ρcp* � � *k* � <sup>1</sup> *ρcp* � � avg � �,

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

interface conditions can be found in [9, 23].

*DOI: http://dx.doi.org/10.5772/intechopen.86252*

**4.2 Modified geometry-ignored interface schemes**

interpolation based scheme will be presented in Section 5.1.

*Sconj* <sup>¼</sup> *<sup>∂</sup> ∂xj*

> *∂xj* 1 *ρcp* � �

it will be further discussed in Section 5 with a numerical test.

*<sup>k</sup>* þ *ρcp* � � *k*þ1

<sup>h</sup> i.2.

*4.2.1 Group 1: sourcing term addition*

given in [34]:

*ρcp* � �

**71**

difference (FD) scheme: *<sup>∂</sup>*

avg ¼ *ρcp* � �

discussed in this chapter.

When **e** *<sup>α</sup>* is aligned with the interface normal directions, Φ*n<sup>α</sup>* ¼ Φ*nf* and Φ*n<sup>α</sup>* ¼ Φ*ns* are readily noticed. Hence, Eqs. (1), (8), (9a), (9b), (10a) and (10b) constitute a linear system of six equations, and the six unknowns *<sup>g</sup><sup>α</sup>* **<sup>x</sup>***<sup>f</sup> ; <sup>t</sup>* <sup>þ</sup> *<sup>δ</sup><sup>t</sup>* � �, *g<sup>α</sup>* **x***<sup>s</sup>* ð Þ *; t* þ *δt* , Φ*df* , Φ*ds*, Φ*nf* , and Φ*ns* can be analytically solved. The interface scheme thus becomes [9]:

$$\begin{split} \mathbf{g}\_{\overline{\alpha}}(\mathbf{x}\_{\mathcal{f}},t+\delta t) &= \mathbf{A}\_{1}^{f}\hat{\mathbf{g}}\_{a}\left(\mathbf{x}\_{\mathcal{f}},t\right) + \mathbf{A}\_{2}^{f}\hat{\mathbf{g}}\_{a}\left(\mathbf{x}\_{\mathcal{f}},t\right) + \mathbf{A}\_{3}^{f}\hat{\mathbf{g}}\_{\overline{\alpha}}(\mathbf{x}\_{\mathcal{f}},t) \\ &+ \mathbf{B}\_{1}^{f}\hat{\mathbf{g}}\_{\overline{\alpha}}(\mathbf{x}\_{\star},t) + \mathbf{B}\_{2}^{f}\hat{\mathbf{g}}\_{\overline{\alpha}}(\mathbf{x}\_{\mathcal{s}},t) + \mathbf{B}\_{3}^{f}\hat{\mathbf{g}}\_{a}(\mathbf{x}\_{\star},t) \,. \end{split} \tag{11a)}$$
 
$$\begin{split} \mathbf{g}\_{a}(\mathbf{x}\_{\mathcal{s}},t+\delta t) &= \mathbf{A}\_{1}^{s}\hat{\mathbf{g}}\_{\overline{\alpha}}(\mathbf{x}\_{\mathcal{t}},t) + \mathbf{A}\_{2}^{s}\hat{\mathbf{g}}\_{\overline{\alpha}}(\mathbf{x}\_{\mathcal{s}},t) + \mathbf{A}\_{3}^{s}\hat{\mathbf{g}}\_{a}(\mathbf{x}\_{\mathcal{t}},t) \\ &+ \mathbf{B}\_{1}^{s}\hat{\mathbf{g}}\_{a}(\mathbf{x}\_{\mathcal{f}},t) + \mathbf{B}\_{2}^{s}\hat{\mathbf{g}}\_{a}(\mathbf{x}\_{\mathcal{f}},t) + \mathbf{B}\_{3}^{s}\hat{\mathbf{g}}\_{\overline{\alpha}}(\mathbf{x}\_{\mathcal{f}},t) \,. \end{split} \tag{11b)}$$

$$+\gamma\_q q\_{\rm jump} + \gamma\_t \phi\_{\rm jump}.$$

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} = \begin{cases} -2\Delta, & (0 \le \Delta \le 0.5), \\ -\frac{1}{2\Delta}, & (\Delta > 0.5), \end{cases} \tag{12a}$$

Scheme 2 : *cd*<sup>1</sup> ¼ �2 1ð Þ � Δ *,* and (12b)

$$\text{Scheme } \mathfrak{Z}: \mathfrak{c}\_{d1} = -\mathbf{1}.\tag{12c}$$

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:

$$\mathbf{g}\_{\overline{a}}(\mathbf{x}\_{\mathcal{f}},t+\delta t) = \left(\frac{\mathbf{1}-\sigma}{\mathbf{1}+\sigma}\right)\hat{\mathbf{g}}\_{a}(\mathbf{x}\_{\mathcal{f}},t) + \left(\frac{2\sigma}{\mathbf{1}+\sigma}\right)\hat{\mathbf{g}}\_{\overline{a}}(\mathbf{x}\_{\iota},t) + \frac{q\_{\text{jump}}}{\mathbf{1}+\sigma} + \frac{\varepsilon\_{D}\sigma\phi\_{\text{jump}}}{\mathbf{1}+\sigma},\tag{13a}$$

$$\mathcal{g}\_a(\mathbf{x}, t + \delta t) = -\left(\frac{\mathbf{1} - \sigma}{\mathbf{1} + \sigma}\right) \hat{\mathcal{g}}\_{\overline{a}}(\mathbf{x}, t) + \left(\frac{\mathbf{2}}{\mathbf{1} + \sigma}\right) \hat{\mathcal{g}}\_a(\mathbf{x}\_f, t) + \frac{q\_{\text{jump}}}{\mathbf{1} + \sigma} - \frac{\varepsilon\_D \phi\_{\text{jump}}}{\mathbf{1} + \sigma}. \tag{13b}$$

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

$$\mathbf{g}\_{\overline{a}}(\mathbf{x}\_{\delta}, t + \delta t) = \hat{\mathbf{g}}\_{\overline{a}}(\mathbf{x}\_{t}, t) + q\_{\text{jump}}/2 + \varepsilon\_{D} \phi\_{\text{jump}}/2,\tag{14a}$$

*Interface Treatment for Conjugate Conditions in the Lattice Boltzmann Method… DOI: http://dx.doi.org/10.5772/intechopen.86252*

$$\mathbf{g}\_a(\mathbf{x}\_\dagger, \mathbf{t} + \delta \mathbf{t}) = \hat{\mathbf{g}}\_a(\mathbf{x}\_\dagger, \mathbf{t}) + q\_{\text{jump}}/2 - \varepsilon\_D \phi\_{\text{jump}}/2. \tag{14b}$$

For straight interfaces where Δ 6¼ 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\_{conj} = \frac{\partial}{\partial \mathbf{x}\_j} \left(\frac{\mathbf{1}}{\rho c\_p}\right) \cdot \left(-k \frac{\partial \phi}{\partial \mathbf{x}\_j} + \rho c\_p u\_j \phi\right) \tag{15}$$

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-

$$\begin{split} \text{difference (FD) scheme:} & \frac{\partial}{\partial \mathbf{x}\_{j}} \Big( \frac{1}{\rho \mathbf{c}\_{p}} \Big)\_{k} = \left[ \left( \frac{1}{\rho \mathbf{c}\_{p}} \right)\_{k} - \left( \frac{1}{\rho \mathbf{c}\_{p}} \right)\_{\text{avg}} \right] \Big/ \mathbf{0.5} \text{\& with } \\ \left( \rho \mathbf{c}\_{p} \right)\_{\text{avg}} &= \left[ \left( \rho \mathbf{c}\_{p} \right)\_{k} + \left( \rho \mathbf{c}\_{p} \right)\_{k+1} \right] / 2. \end{split}$$

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" *<sup>h</sup>*<sup>∗</sup> <sup>¼</sup> *<sup>ρ</sup>cp* � � <sup>0</sup>*ϕ*, where *ρcp* � � <sup>0</sup> is a reference heat capacity, the governing CDE (2) becomes [27, 28]:

$$\frac{\partial h^{\*}}{\partial t} + \nabla \cdot (\mathbf{u}h^{\*}) = \nabla \cdot \left(\frac{k}{\rho c\_{p}} \nabla h^{\*}\right) - \frac{k}{\left(\rho c\_{p}\right)\_{0}} \nabla h^{\*} \cdot \nabla \left(\frac{\left(\rho c\_{p}\right)\_{0}}{\rho c\_{p}}\right) - \frac{h^{\*} \left(\rho c\_{p}\right)\_{0}}{\rho c\_{p}} \mathbf{u} \cdot \nabla \frac{\rho c\_{p}}{\left(\rho c\_{p}\right)\_{0}} \tag{16}$$

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 <sup>∇</sup>*<sup>j</sup> <sup>f</sup>* <sup>¼</sup> <sup>1</sup> *cs* <sup>2</sup>*δ<sup>t</sup>* ∑*<sup>α</sup> w<sup>α</sup> f*ð Þ **x** þ ð Þ **e***α*�**x** *δt e<sup>α</sup>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]

$$\mathbf{g}\_a^{\text{eq}} = \begin{cases} \phi \left( \eta \mathbf{c}\_p - \mathbf{c}\_{p0} \right) + \alpha\_a \phi c\_p \left[ c\_{p0} / c\_p + (\mathbf{e}\_a \cdot \mathbf{u}) / c\_s^2 \right], & a \neq \mathbf{0} \\\\ \alpha\_a \phi c\_p \left[ c\_{p0} / c\_p + (\mathbf{e}\_a \cdot \mathbf{u}) / c\_s^2 \right], & a = \mathbf{0} \end{cases} \tag{17}$$

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

$$\mathbf{m}^{\text{eq}} = \begin{bmatrix} c\_p \phi, & \omega c\_p \phi, & \nu c\_p \phi, & 4c\_p \phi - \frac{\mathbf{10}}{3} c\_{p0} \phi, & \mathbf{0} \end{bmatrix}^T. \tag{18}$$

*E*<sup>2</sup> ¼ ∑ *x, y*

*x, <sup>y</sup>*¼*<sup>h</sup>*

*<sup>∂</sup><sup>y</sup>* <sup>j</sup>LBE � *<sup>D</sup>*1*,* <sup>2</sup>

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

*<sup>∂</sup><sup>y</sup>* <sup>j</sup>ex � �<sup>2</sup>

*∂ϕ*1*,*<sup>2</sup>

*E*2*, <sup>ϕ</sup>*int ¼ ∑

*D*1*,* <sup>2</sup>

*E*2*,* qint ¼ ∑

**Figure 2.**

*x, y*¼*h*

*(h* ≤ *y* ≤ *H). With permission from [23].*

quantities follows that in [21, 23].

number, *Pe*, is defined as *Pe = UH/D*1.

(see [23] for details).

**73**

*∂ϕ*1*,*<sup>2</sup> *∂t*

þ *U*

*∂ϕ*1*,*<sup>2</sup>

*<sup>∂</sup><sup>x</sup>* <sup>¼</sup> *<sup>D</sup>*1*,* <sup>2</sup>

in the *x*-direction. Taking into account the standard conjugate conditions:

*<sup>ϕ</sup>*LBE � *<sup>ϕ</sup>*ex ð Þ<sup>2</sup>

*Schematic layout of the lattice on a 2D channel containing two fluids in domain 1 (0* ≤ *y* ≤ *h) and domain 2*

*Interface Treatment for Conjugate Conditions in the Lattice Boltzmann Method…*

*DOI: http://dx.doi.org/10.5772/intechopen.86252*

*<sup>ϕ</sup>*1*,* <sup>2</sup>jLBE � *<sup>ϕ</sup>*1*,* <sup>2</sup>jex � �<sup>2</sup>

*∂ϕ*1*,* <sup>2</sup>

where *E*<sup>2</sup> evaluates the overall error in all the interior lattice nodes in the two domains, *E*2, *<sup>ϕ</sup>*int and *E*2, qint evaluate the respective relative errors of the interfacial *ϕ* value and its flux. For circular interfaces, *E*2, *<sup>ϕ</sup>*int and *E*2, qint are evaluated at the interface nodes along the curved geometry. The computation of the interfacial

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

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 *ϕ*1ð Þ¼ *x; y* ¼ 0 *ϕ*2ð Þ¼ *x; y* ¼ *H* cos 2ð Þ *πx=L* , and periodic conditions

*ϕ*<sup>1</sup> = *ϕ*2, *k*1*∂ϕ*1*/∂y* = *k*2*∂ϕ*2*/∂y* at *y* = *h*, the analytical solution to Eq. (22) can be solved

*∂*2 *ϕ*1*,* <sup>2</sup> *∂x*<sup>2</sup> þ *∂*2 *ϕ*1*,*<sup>2</sup> *∂y*<sup>2</sup>

� � (22)

*<sup>∂</sup><sup>y</sup>* <sup>j</sup>ex � �<sup>2</sup> " #<sup>1</sup>*=*<sup>2</sup>

ex " #<sup>1</sup>*=*<sup>2</sup>

#<sup>1</sup>*=*<sup>2</sup> 2

*=* ∑ *x, y ϕ*2

> *=* ∑ *x, y*¼*h*

*=* ∑ *x, y*¼*h*

4 (20)

*<sup>ϕ</sup>*1*,* <sup>2</sup>jex � �<sup>2</sup>

*D*1*,* <sup>2</sup>

*∂ϕ*1*,* <sup>2</sup>

(19)

(21)

The temperature is solved from *ϕ* ¼ ∑*<sup>α</sup> 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 *<sup>τ</sup>* as *<sup>D</sup>* <sup>¼</sup> ð Þ *<sup>τ</sup>* � <sup>0</sup>*:*<sup>5</sup> *<sup>c</sup>*<sup>2</sup> *<sup>s</sup> δ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 *<sup>τ</sup>* as *<sup>k</sup>* <sup>¼</sup> *cp*0ð Þ *<sup>τ</sup>* � <sup>0</sup>*:*<sup>5</sup> *<sup>c</sup>*<sup>2</sup> *<sup>s</sup> δ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–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 [23]:

*Interface Treatment for Conjugate Conditions in the Lattice Boltzmann Method… DOI: http://dx.doi.org/10.5772/intechopen.86252*

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

$$E\_2 = \left[\sum\_{\mathbf{x}, \mathbf{y}} (\phi\_{\text{LBE}} - \phi\_{\text{ex}})^2 / \sum\_{\mathbf{x}, \mathbf{y}} \phi\_{\text{ex}}^2\right]^{1/2} \tag{19}$$

$$E\_{2,\ \phi\text{int}} = \left[ \sum\_{\mathbf{x}, \mathbf{y} = \mathbf{h}} \left( \phi\_{\mathbf{1}, 2} |\_{\text{LEE}} - \phi\_{\mathbf{1}, 2} |\_{\text{ex}} \right)^2 / \sum\_{\mathbf{x}, \mathbf{y} = \mathbf{h}} \left( \phi\_{\mathbf{1}, 2} |\_{\text{ex}} \right)^2 \right]^{1/2} \tag{20}$$

$$E\_{2,\text{ quint}} = \left[ \sum\_{\mathbf{x}, \mathbf{y} = \mathbf{h}} \left( D\_{\mathbf{1},2} \frac{\partial \phi\_{\mathbf{1},2}}{\partial \mathbf{y}} |\_{\text{LEE}} - D\_{\mathbf{1},2} \frac{\partial \phi\_{\mathbf{1},2}}{\partial \mathbf{y}} |\_{\text{ex}} \right)^2 \Big/ \sum\_{\mathbf{x}, \mathbf{y} = \mathbf{h}} \left( D\_{\mathbf{1},2} \frac{\partial \phi\_{\mathbf{1},2}}{\partial \mathbf{y}} |\_{\text{ex}} \right)^2 \right]^{1/2} \tag{21}$$

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

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

$$\frac{\partial \phi\_{1,2}}{\partial t} + U \frac{\partial \phi\_{1,2}}{\partial \mathbf{x}} = D\_{1,2} \left( \frac{\partial^2 \phi\_{1,2}}{\partial \mathbf{x}^2} + \frac{\partial^2 \phi\_{1,2}}{\partial \mathbf{y}^2} \right) \tag{22}$$

We consider only the steady case with sinusoidal boundary conditions on the horizontal walls *ϕ*1ð Þ¼ *x; y* ¼ 0 *ϕ*2ð Þ¼ *x; y* ¼ *H* cos 2ð Þ *πx=L* , and periodic conditions in the *x*-direction. Taking into account the standard conjugate conditions: *ϕ*<sup>1</sup> = *ϕ*2, *k*1*∂ϕ*1*/∂y* = *k*2*∂ϕ*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*.

discontinuity in approximating the heat capacity gradient in Groups 1 and 2 when

*Interface Treatment for Conjugate Conditions in the Lattice Boltzmann Method…*

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

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

0.006 0.02 0.04 0.08

*Relative L2 norm error, E2, for the interior scalar versus the grid resolution, 1/H, for steady convection-*

0.006 0.02 0.04 0.08

*Relative L2 norm error, E2, ϕint, for the interfacial scalar versus the grid resolution, 1/H, for steady convection-*

Group 1 Group 2 Group 3 Ref. [23]

> Group 1, int,1 Group 1, int,2 Group 2, int,1 Group 2, int,2 Group 3, int,1 Group 3, int,2 Ref. [23], int,1 Ref. [23], int,2

<sup>1</sup> 6¼ 1. This will be further demonstrated in the *L*<sup>2</sup> norm errors.

*σ* ¼ *ρcp* 

<sup>2</sup>*= ρcp* 

*DOI: http://dx.doi.org/10.5772/intechopen.86252*

3 with (*τ*2�0.5)/(*τ*1�0.5) = *k*21.

10-4

10-4

*diffusion in the channel at Δ = 0.5.*

10-3

10-2

10-1

10<sup>0</sup>

*diffusion in the channel at Δ = 0.5.*

**Figure 5.**

**Figure 6.**

**75**

10-3

10-2

10-1

**Figures 3a** and **b** show the scalar value and flux profiles at selected locations along the *x*-direction using a diffusivity ratio *D*<sup>21</sup> = *D*2/*D*<sup>1</sup> = 2 and a thermal conductivity ratio *k*<sup>21</sup> = *k*2/*k*<sup>1</sup> = 3 (consequently *σ* ¼ *ρcp* <sup>2</sup>*= ρcp* <sup>1</sup> = 1.5). Note that Δ = 0.5 is used for all schemes and thus the geometry effect is not included. Other simulation parameters are: *τ*<sup>1</sup> = 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

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

*Interface Treatment for Conjugate Conditions in the Lattice Boltzmann Method… DOI: http://dx.doi.org/10.5772/intechopen.86252*

discontinuity in approximating the heat capacity gradient in Groups 1 and 2 when *σ* ¼ *ρcp* <sup>2</sup>*= ρcp* <sup>1</sup> 6¼ 1. This will be further demonstrated in the *L*<sup>2</sup> norm errors.

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

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

#### **Figure 5.**

*Relative L2 norm error, E2, for the interior scalar versus the grid resolution, 1/H, for steady convectiondiffusion 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 convectiondiffusion 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 convectiondiffusion in the channel at Δ = 0.5.*

To further demonstrate the necessity of taking into account the local geometry to preserve second-order accuracy, **Figure 8** presents the *L*<sup>2</sup> norm errors for the interior field at different Δ values. The parameters used are *D*<sup>21</sup> = 5, *k*<sup>21</sup> = 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, *cd*<sup>1</sup> = �1 is used.

layout and computational domain for 2D diffusion within a circular domain. On the

*ϕ*<sup>1</sup> ¼ *ϕ*<sup>2</sup> þ *ϕjump* ¼ *ϕ*<sup>2</sup> þ *ϕ*<sup>0</sup> cos *φ* at *r* ¼ *R*<sup>1</sup> (23)

þ *q*<sup>0</sup> cos *φ* at *r* ¼ *R*<sup>1</sup> (24)

*∂ϕ*<sup>2</sup> *∂r*

the analytical solution can be found [9]. For the numerical LBM computation, the parameters are *R*2*/R*<sup>1</sup> = 2, *D*2*/D*<sup>1</sup> = 10, *τ*<sup>1</sup> ¼ 0*:*525, and *τ*<sup>2</sup> ¼ 0*:*75. As previously mentioned, three schemes were presented in [21] for adjustable parameters of *cd*1,

*Relative L2 norm error, E2, for the interior temperature versus the grid resolution for 2D diffusion in a circular*

outer circle with radius, *R*2, a Dirichlet boundary condition is applied as *ϕ*2(*r* = *R*2) = cos*φ*. With the following conjugate conditions at the interface:

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

*Interface Treatment for Conjugate Conditions in the Lattice Boltzmann Method…*

*DOI: http://dx.doi.org/10.5772/intechopen.86252*

þ *qjump* ¼ �*D*<sup>2</sup>

and

**Figure 10.**

**77**

*domain. With permission from [9].*

**Figure 9.**

�*D*<sup>1</sup>

*∂ϕ*<sup>1</sup>

*<sup>∂</sup><sup>r</sup>* ¼ �*D*<sup>2</sup>

*∂ϕ*<sup>2</sup> *∂r*

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 Δ 6¼ 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 Δ 6¼ 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

#### **Figure 8.**

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

*Interface Treatment for Conjugate Conditions in the Lattice Boltzmann Method… DOI: http://dx.doi.org/10.5772/intechopen.86252*

#### **Figure 9.**

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

layout and computational domain for 2D diffusion within a circular domain. On the outer circle with radius, *R*2, a Dirichlet boundary condition is applied as *ϕ*2(*r* = *R*2) = cos*φ*. With the following conjugate conditions at the interface:

$$
\phi\_1 = \phi\_2 + \phi\_{jump} = \phi\_2 + \phi\_0 \cos \rho \text{ at } r = R\_1 \tag{23}
$$

and

$$-D\_1 \frac{\partial \phi\_1}{\partial r} = -D\_2 \frac{\partial \phi\_2}{\partial r} + q\_{jump} = -D\_2 \frac{\partial \phi\_2}{\partial r} + q\_0 \cos \rho \text{ at } r = R\_1 \tag{24}$$

the analytical solution can be found [9]. For the numerical LBM computation, the parameters are *R*2*/R*<sup>1</sup> = 2, *D*2*/D*<sup>1</sup> = 10, *τ*<sup>1</sup> ¼ 0*:*525, and *τ*<sup>2</sup> ¼ 0*:*75. As previously mentioned, three schemes were presented in [21] for adjustable parameters of *cd*1,

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

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

*Interface Treatment for Conjugate Conditions in the Lattice Boltzmann Method…*

LL acknowledges the support and start-up fund from the Department of

Department of Mechanical Engineering, Mississippi State University, Mississippi,

© 2019 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium,

\*Address all correspondence to: likeli@me.msstate.edu

provided the original work is properly cited.

effective interface schemes for curved geometries.

*DOI: http://dx.doi.org/10.5772/intechopen.86252*

Mechanical Engineering at Mississippi State University.

**Acknowledgements**

**Author details**

USA

**79**

David Korba and Like Li\*

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

those three schemes in Eqs. (12a)–(12c) are also applied here. For the flux jump case, the jump conditions are set as *ϕ*<sup>0</sup> = 0 and *q*<sup>0</sup> *= D*2*/R*1; and for the temperature jump case, those are *ϕ*<sup>0</sup> = 0.5 and *q*<sup>0</sup> = 0.

**Figures 10** and **11a, b** provide the respective relative *L*<sup>2</sup> norm errors for the interior temperature and the interfacial temperature and flux with the conjugate schemes implemented.

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

*Interface Treatment for Conjugate Conditions in the Lattice Boltzmann Method… DOI: http://dx.doi.org/10.5772/intechopen.86252*

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.
