Lattice Boltzmann methods for complex micro-flows: applicability and limitations for practical applications

The non-slip wall boundary conditions used in the continuum LBM are based on perfect reflection, so the velocity and the temperature of a wall are not reflected in the distribution of the reflected particles. However, from a microscopic viewpoint, the wall boundary condition should include the physics on the wall because the fluid and the wall molecules interact with each other. Therefore, the incident particles are modeled to reflect the information of the Maxwell distribution function at the wall boundary. The modeled form (Ansumali and Karlin 2002) is written in the LBM frame as

$$\begin{eqnarray} &&\fl f_{{\alpha}} ({\bf{x}},t) = \frac{{\sum\nolimits_{{{\alpha}} '} {|({\boldsymbol{\xi }}'_{{\alpha}} - {\bf{u}}_{\rm{w}} ) \cdot {\bf{n}}|} f_{{{\alpha}} '} ({\bf{x}},t)}}{{\sum\nolimits_{{{\alpha}} '} {|({\boldsymbol{\xi }}'_{{\alpha}} - {\bf{u}}_{\rm{w}} ) \cdot {\bf{n}}|} f_{{{\alpha}} ',{\bf{w}}}^{\rm eq} ({\bf{x}},t)}}f_{{{\alpha}} ,{\bf{w}}}^{{\rm eq}} ({\bf{x}},t) = f_{{\alpha}} ^{\rm D} ({\bf{x}},t), \nonumber\\ &&{\rm if}\;\left[ {({\boldsymbol{\xi}}'_{{\alpha}} - {\bf{u}}_{\rm{w}} ) \cdot {\bf{n}} < 0;{\rm{ }}({\boldsymbol{\xi }}_{{\alpha}} - {\bf{u}}_{\bf{w}} ) \cdot {\bf{n}} > 0} \right], \end{eqnarray} \tag{ 20 }$$

where n is the unit wall normal vector, ξ '_α is the velocity of incident particles, f_{α ,w}^eq is the wall equilibrium distribution function, and the subscripts w,α ',α respectively denote the wall and the directions of the incident and reflected particles.

Succi (2002) proposed the combined bounce-back and specular-reflection boundary condition for slippage at a solid wall in a micro-flow. In this scheme, it is assumed that some of the particles hitting the solid wall bounce back but the others reflect specularly. It is expressed as

$$\begin{equation} f_{{\alpha}} ({\bf{x}},t + \delta t) = r_{\rm{s}} {\rm{ }}f_\beta ({\bf{x}},t) + (1 - r_{\rm{s}} ){\rm{ }}f_\gamma ({\bf{x}},t), \end{equation} \tag{ 21 }$$

where r_s is a bridge coefficient taking r_s = 0.0 - 1.0 and vectors ξ_β  = - ξ_α and ξ_γ  = ξ_α - 2(ξ_α ⋅ n)n are, respectively, the inverse and the specular symmetric velocity vectors of ξ_α .

In the context of the MRT model, by discussing the second-order slip boundary condition for the microscopic slip velocity, Guo et al (2008) derived the bridge coefficient as

$$\begin{equation} r_{\rm s} = \left[ {1 + \chi A_1 + \frac{{\tau '_{\rm s} (0)\delta _x }}{{8(\tilde\tau _{\rm s} (0))^2 }}} \right]^{ - 1} , \end{equation} \tag{ 22 }$$

where τ_s = 1/s_ν, $\tilde\tau _{\rm s} = \tau _{\rm s} - 0.5$, τ_s (0) = τ_s |_wall , τ '_s (0) = ∂ τ_s /∂n|_wall and n is the wall-normal direction. The coefficients are $\chi = \sqrt {\pi /6}$, $A_1 = \frac{{2 - \sigma _{\rm{a}} }}{{\sigma _{\rm{a}} }}(1 - 0.1817\sigma _{\rm{a}} )$ and σ_a is the accommodation coefficient.

The diffusive bounce-back boundary condition is the combination of the above diffuse-scattering and bounce-back boundary conditions for the wall boundary as

$$\begin{equation} f_{{\alpha}} ({\bf{x}},t + \delta t) = (1 - r){\rm{ }}f_\beta ({\bf{x}},t) + r{\rm{ }}f_{{\alpha}} ^{\rm{D}} ({\bf{x}},t + \delta t), \end{equation} \tag{ 23 }$$

where r is a probability coefficient taking r = 0.0 - 1.0 , the bounce-back vectors are ξ_β = - ξ_α and f_α ^D is the diffuse scattering model given by equation (20). This model turns into the diffuse scattering model when r = 1 , while it becomes the (non-slip) bounce-back model with r = 0 .

In the context of the MRT LBM with the D2Q9 model, the DBBC describes the velocity at near-wall lattice node 2 as

$$\begin{equation} U_2 = A_1 U_1 + A_2 \delta ta_x , \end{equation} \tag{ 24 }$$

where U_j is the velocity at lattice node j and a_x is the acceleration parallel to the wall (i.e. a = (a_x ,0) ). When lattice nodes 1 and 2 are located at δ x/2 and 3δ x/2 from the wall (the wall is situated halfway between nodes 0 and 1 as shown in figure 2), Chai et al (2010) obtained the parameters A₁ and A₂ as

$$\begin{equation} \fl A_1 = \frac{{3 - 2r + r\tau _\nu }}{{1 - r + r\tau _\nu }},\quad A_2 = \frac{{r(11/2 - 2\tau _q - 8\tau _\nu + 4\tau _q \tau _\nu ) + 4(\tau _q + \tau _\nu ) - 8\tau _q \tau _\nu - 5}}{{(2\tau _\nu - 1)(1 - r + r\tau _\nu )}}, \end{equation} \tag{ 25 }$$

where τ_ν = s_ν ^{- 1} and τ_q = s_q^{- 1} . In a fully developed 2D Poiseuille flow in the slip flow regime, by integrating the momentum equation with the second-order boundary condition,

$$\begin{equation} U|_{\rm wall} = C_1 \lambda \left. {\frac{{\partial u}}{{\partial y}}} \right|_{{\rm wall}} - C_2 \lambda ^2 \left. {\frac{{\partial ^2 u}}{{\partial y^2 }}} \right|_{{\rm wall}} . \end{equation} \tag{ 26 }$$

Zoom In Zoom Out Reset image size

For coefficients C₁ and C₂ , there are many proposals by the kinetic theory (e.g. Mitsuya 1993, Beskok et al 1996, Hadjiconstantinou 2003). When the velocity distribution is assumed to be parabolic, the velocity can be described as

$$\begin{equation} \frac{U}{{U_{\rm c} }} = 4\frac{y}{H}\left( {1 - \frac{y}{H}} \right) + \frac{{U_{\rm s} }}{{U_{\rm c}}}, \end{equation} \tag{ 27 }$$

where H is the channel height and U_c is the centerline velocity: U_c = a_x H² /(8ν ) . Thus, the slip velocity at the wall can be obtained as

$$\begin{equation} U_{\rm s} /U_{\rm c} = 4C_1 Kn + 8C_2 Kn^2 . \end{equation} \tag{ 28 }$$

By equations (24), (25) and (27), one can obtain

$$\begin{eqnarray} U_{\rm s }/U_{\rm c} &= & \frac{{4r}}{{2 - r}}\left( {\tau _\nu - \frac{1}{2}} \right)\frac{{\delta x}}{H} + \left( {\frac{{16(\tau _q - 1/2)(\tau _\nu - 1/2) - 3}}{3}} \right)\left( {\frac{{\delta x}}{H}} \right)^2 \nonumber\\ & = & \left( {\frac{{4r}}{{2 - r}}\sqrt {\frac{6}{\pi }} } \right)Kn + \left( {\frac{{32(\tau _q - 1/2)(\tau _\nu - 1/2) - 6}}{{(\tau _\nu - 1/2)^2 \pi }}} \right)Kn^2 . \end{eqnarray} \tag{ 29 }$$

Comparing equations (28) and (29), one can derive

$$\begin{equation} r = \frac{{2C_1 }}{{\sqrt {\frac{6}{\pi }} + C_1 }}, \end{equation} \tag{ 30 }$$

$$\begin{equation} \tau _q = \frac{1}{2} + \frac{{3 + \pi (2\tau _\nu - 1)^2 C_2 }}{{8(2\tau _\nu - 1)}}, \end{equation} \tag{ 31 }$$

as in Chai et al (2010).

Verhaeghe et al (2009) applied only the first-order terms of equations (28) and (29) and their expression for r was

$$\begin{equation} 1 - r = \frac{{3\mu - Kn\,Hc\,\rho _0 }}{{3\mu + Kn\,Hc\,\rho _0 }}, \end{equation} \tag{ 32 }$$

where c = δx/δt. (Equation (30) can be rewritten as equation (32) when c = C₁.) With this first-order slip velocity condition, Verhaeghe et al (2009) applied

$$\begin{equation} s_q = \frac{{8(2 - s_\nu )}}{{8 - s_\nu }}, \end{equation} \tag{ 33 }$$

following Ginzbourg and Adler (1994). When r = 0 for the pure bounce-back boundary condition, it is recognized that this condition of equation (33) guarantees U_s = 0 by equation (29).

Niu et al (2007) and Suga et al (2010b) confirmed that the D2Q21 LBGK LBM with the DSBC reasonably captured Knudsen number effects in channel flows. This can be confirmed in figure 4 where the velocity profiles of plane Poiseuille channel flows at Kn = 0.01–10 are compared with the data of the DSMC (Beskok and Karniadakis 1999). (Figures are re-plotted for this paper with re-calculated data.) These velocity profiles are normalized by the bulk mean velocity U_b . The 2D uniform Cartesian lattice whose wall normal node points are 100 (this density of the lattice was confirmed to be more than good enough) is used for the simulations. The correction function of the molecular mean free path Ψ 1 and the regularization procedure are coupled with the LBGK LBM. It can be seen that the less accurate discrete velocity model of the D2Q9 model also captures the increasing level of the slip velocity due to Kn and velocity profiles, although its results always deviate a little from those of the DSMC very near the walls, particularly when Kn ⩾ 0.1 . This is partly because in the D2Q9 model, only up to second-order terms are retained in the Hermite expansion of the Maxwellian distribution while the D2Q21 model retains terms of third-order in the Hermite expansion and uses a Gauss–Hermite quadrature of third-order to obtain the discrete velocity values. Note that the D2Q21 and D3Q39 models require three levels of particle velocity. For particles with two and three unit speeds, when the DSBC is applied, it is necessary to provide information of nodes inside the wall (pseudo nodes). In such cases, the wall surface values are simply taken for those pseudo nodes.

Zoom In Zoom Out Reset image size

Figure 5 compares the dependence of the flow characteristics on Kn. Figure 5(a) compares the slippage and centerline velocities: U_s and U_c . The D2Q21 model shows better agreement with the reference data (Karniadakis et al 2005). As shown in figure 5(b), as for the distribution of the normalized mass flow rate,

$$\begin{equation} Q = \frac{1}{{(\rho _0 a_x H^2 /c_{\rm s} )}}\int_0^H {\rho U\,{\rm{d}}y} , \end{equation} \tag{ 47 }$$

where a_x is the streamwise acceleration, both the D2Q21 and D2Q9 models predict reasonable results at least when Kn < 1. (Some comments on the slippage velocity by the LBGK LBM are in the Appendix.)

Zoom In Zoom Out Reset image size

In order to know the performance of the MRT LBMs with the DBBC to predict slippage velocities and flow rates, the same flow fields were discussed by Suga and Ito (2011). The compared LBMs were the MRT LBMs with the D2Q9 by the first-order and second-order slip-flow models. The first-order slip-flow model was the one proposed by Verhaeghe et al (2009) whereas one of the second-order models was Hadjiconstantinou's (2003) model, which applies C₁ = 1.11, C₂ = 0.61 , following the recommendation by Chai et al (2010). Another second-order model was Mitsuya's model, which was denoted the 1.5th-order model following the original naming by Mitsuya (1993). It applies C₁ = 1.0, C₂ = 2/9 . These coefficients are summarized in table 2.

Figure 6(a) compares the slip and centerline velocities: U_s and U_c . Although each agreement is not good enough at Kn > 0.1 (U_s is overpredicted up to 20% of U_b and U_c is underpredicted up to 10% of U_b ), the profiles of U_s and U_c by each model correspond well with the reference data at Kn ≈ 0.1. However, as shown in figure 6(b), for the mass flow rate the 1.5th-order slip-flow model leads to the best agreement among the slip-flow models. Interestingly, each profile shown in figure 6(b) corresponds well to the analytically obtained distributions in Mitsuya (1993). Only the different points in the model terms are the values of the probability coefficient r and the relaxation time τ_q for the energy flux. By equation (30), the probability coefficient is fixed as r = 0.84,0.84 and 0.89 by the first-order, 1.5-order and second-order slip-flow models, respectively. The relaxation time τ_q varies a lot depending on the Knudsen number. At Kn = 0.01, τ_q = 0.64,0.88 and 1.29 while it is 0.50, 239.26 and 655.96 at Kn = 10 by the first-order, 1.5-order and second-order slip-flow models, respectively, when the channel height is divided into 99 lattice spaces. Consequently, at Kn ≪ 0.1 , the difference in the values of r or τ_q is not significantly large and this results in nearly the same results of the flow characteristics. However, at Kn ≫ 1 , the difference in the values of τ_q becomes very large and accordingly it affects the predicted behavior such as the flow rate shown in figure 6(b).

Zoom In Zoom Out Reset image size

It is recognized that the regularization procedure filters the majority of such differences out. As shown in figures 7(a) and (b), with the regularization procedure, all the slip-flow models produce nearly the same flow characteristics even at Kn > 1. This means that with the regularization, the order of the slip-flow model may not be very important. A similar effect to that shown in figures 6(b) and 7(b) was reported by Zhang et al (2006b) for the same force-driven Poiseuille flows. They discussed the effect of the regularization on the flow rate tendency by the D2Q9 BGK model with the bounce-back boundary condition and found that by the regularization process, only the Navier–Stokes order effects were preserved in their result (which was much worse than in figure 7(b)). Their conclusion was that since the regularization enforced the system to be confined within the second-order Hermite moment space in the D2Q9, all the higher-order non-equilibrium contributions from both the numerical and physical processes responsible for the finite Knudsen phenomena were filtered out. Although their sub-models for the micro-flows such as for the relaxation time and the boundary condition are different from the results shown in figure 7, due to the similar mechanisms, it is considered that the difference in the slip-flow models is reduced. The effects of the regularization become more obvious in the complex flow geometry discussed later.

Zoom In Zoom Out Reset image size

Since the LBE does not converge to the BE, it is difficult to capture the Knudsen layer effects precisely. However, due to the models incorporated, it has often been reported that the LBM can capture them. Therefore, in this part of the section, the model performance is evaluated by the numerical solutions.

The near-wall velocities particularly with the correction Ψ of the mean free path are expected to show different profiles from the parabolic ones since Ψ was designed to take into account the characteristic Knudsen-layer effects. In order to investigate the Knudsen-layer correction by Sone (2007), analyses using the Couette flow simulations are performed as in Suga et al (2010b). In the Couette flows at a constant temperature, the Knudsen-layer correction of Sone reduces to

$$\begin{equation} \left[ \begin{array}{l} U_{{\rm{Gw}}} - U_{\rm{w}} \\ U_{\rm{K}} \\ \end{array} \right] = - \frac{{\sqrt \pi }}{2}Kn\frac{{{\rm{d}}U_{\rm{G}} }}{{{\rm{d}}y}}\left[ \begin{array}{l} k_0 \\ Y_0 (\eta _{{\rm{KL}}} ) \\ \end{array} \right], \end{equation} \tag{ 48 }$$

where U_Gw - U_w is the virtual slip velocity which is the difference between the wall velocity U_w and the velocity U_G at the wall. As illustrated in figure 8(a), U_G is the velocity extrapolated from the Navier–Stokes velocity profile in the bulk region. The Knudsen-layer correction of the velocity U_K is the difference between the velocity in the Knudsen layer and U_G : U_K = U- U_G . As shown in figure 8(b), both correction models Ψ 1 (equation (37)) and Ψ 2 (equation (39)) seem to reproduce well the solution of Sone (2007). Here the first-order slip flow DBBC and the regularization are applied to the MRT LBM by the D2Q9. However, Ψ 1 and Ψ 2 models predict the slip coefficient as k₀ = - 1.0664 and - 1.1705 at Kn = 0.11 while it is k₀ = - 1.01619 for the BGK model irrespective of Kn by Sone (2007). Moreover, as shown in figure 8(c), the predicted profiles of the Knudsen-layer function Y₀ (η_KL ) are too large or one order smaller than that obtained by Sone. Here, η_KL is the Knudsen-layer coordinate defined as $\eta _{{\rm{KL}}} = y/(Kn\sqrt \pi /2)$. (The profile of Ψ 1 by the D2Q9 model is almost invisible at the scale used.) In figure 8(c), the results by the LBGK LBMs with the DSBC are also compared. Generally, all the shown tendencies by Ψ 2 are the same and it is confirmed that the results are insensitive to the discrete velocity models except in the vicinity of the wall. Although the profiles are not shown here, at Kn = 1.1 even with the D2Q21, both the correction models predict the profiles that are not visible with the same scale of figure 8(c). These results are consistent with the discussions by Watari (2010), who proved that to reproduce the relaxation process in the Knudsen layer, the triple octagon model with 24 velocity directions, which was too elaborate to apply for complex geometry, was required.

Zoom In Zoom Out Reset image size

Therefore, it is confirmed that even though some correction model or wall function seems to work well to predict the Knudsen layer effects, we cannot generally expect much for the accuracy of the prediction of the Knudsen layer. This is consistent with the comments by Luo (2011). However, the correction of the mean free path is useful particularly for predictions of the slip velocity at Kn > 0.1 as shown in figure 9(a) whereas the flow rates are rather insensitive as in figure 9(b).

Zoom In Zoom Out Reset image size

Suga et al (2010b) discussed the flows in a nanochannel with a square cylinder placed at the center region, which is shown in figure 10 and called square cylinder flow hereafter. Periodic boundary conditions were applied to each direction and thus the flow regime was regarded as a part of an infinite cylinder array set in a nanochannel. The force was applied to the x-direction only and the simulations were carried out on the 2D grid whose size was 100 × 100 , which was confirmed to be good enough for grid-independent solutions.

Zoom In Zoom Out Reset image size

Figure 11 shows contour maps of the streamwise velocity magnitude calculated by the LBGK LBM with the DSBC. As is clearly seen, without the regularization procedure, the predicted flow field shows unphysical patterns around the cylinder (figure 11(b)). Although this tendency was not so strong, it was also confirmed in the lower Kn cases. When the projection of the LBE into the Hermite space is considered, it is understood that the collision step introduces an error because the distribution function does not automatically lie entirely within the Hermite space. In the case that the system is not far from equilibrium, such an error is small and can be ignored (Zhang et al 2006b). Latt and Chopard (2006) argued that since the number of distribution functions was normally greater than that of flow variables, the set of unknowns could be considered as overdetermined and thus the non-equilibrium part of the distribution included higher-order moments in principle. It is considered that this contribution of the non-equilibrium part of the distribution tends to be greater as Kn increases. However, the regularization procedure filters such higher-order non-equilibrium moments out, and the resultant flow field becomes dealiased. It can thus be considered that even though LBMs without the regularization work reasonably well at moderate Kn in wall parallel flows (the majority of the LBMs in the literature are in this category), their performance in complex micro-/nano-flows may be questionable.

Zoom In Zoom Out Reset image size

Figure 12 compares the normalized streamwise velocity U profiles by the LBGK model with the DSBC, the regularization and Ψ 1 at x/H = 0.0, 0.25, 0.5, 0.75 for the case of Kn = 0.11. The distribution profiles of the D2Q9 and D2Q21 models agree well with those of the MD simulation. Since it can be seen that the profiles of the D2Q9 model are also acceptable and the maximum difference from those of the D2Q21 model is less than 5% of the bulk mean velocity, the conventional D2Q9 model can also be applicable to simulate the velocity profiles of the complex micro-flow. Although Suga and Ito (2011) confirmed that slightly more satisfactory results were obtained by the MRT LBM with the D2Q9 model as shown in figure 13, the results of the LBGK models in figure 12 are also acceptable. Note that it is clear in figure 13 that with the regularization, the slip-flow models applied to the DBBC do not provide any difference as in the Poiseuille flows.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For discussing the flow characteristics by the wall boundary conditions in a developing nano-flow, flow in a nanochannel with an inserted flat plate was considered in Suga et al (2010a). As shown in figure 14, the mother channel with height H and domain length L had an inserted flat plate, which made another narrower channel region. The thickness and the length of the plate were 0.15H and L/2, respectively. The gap (the height of the finite narrower channel) between the plate and the channel bottom wall was 0.125H. The drive force and periodic boundary conditions were applied to the x-direction, and thus the flow regime was regarded as a fully developed force-driven nanochannel flow combined with developing narrower nanochannel flows inside. The considered Knudsen number was Kn = 0.14 based on the characteristic length H. The LBGK LBM simulations were carried out on a grid of 187 × 97 , which was confirmed to be good enough.

Zoom In Zoom Out Reset image size

Figure 15 compares the normalized streamwise velocity U profiles at x/L = 0.0 - 0.76 . Both the LBGK simulations were coupled with the regularization and Ψ 1 . The results with the DSBC agree with the profiles of the MD simulation with an error of up to 10% of U_b in the wider gap region. Although the slip and centerline velocities of Poiseuille flows were confirmed to be satisfactorily predicted by the BSBC model with b = 0.7 at Kn ⩽ 0.14 , the BSBC model produces less slippage velocities and greater centerline velocities (with up to 20% of U_b ) in figure 15. The initial flow development in the narrower channel at x/L = 0.22 does not satisfactorily agree with the MD profile although the DSBC provides better results. Hence, it can be recognized that the bridging coefficient for the BSBC model needs some dependence on local flow quantities.

Zoom In Zoom Out Reset image size

To discuss further the applicability to the complex flow geometry, as shown in figure 16, Suga and Ito (2011) considered a force-driven flow in a nanochannel with a triangular cylinder at the center region, which is called triangular cylinder flow hereafter. To describe the slanted faces of the triangular cylinder, an interpolation scheme was applied since Peng and Luo (2008) concluded that the accuracy of the interpolated bounce-back (IBB) boundary condition was better than that of the immersed-boundary method. Although they chose the quadratic IBB, there was no significant difference between the quadratic and linear IBB in Pan et al (2006). Thus, the linear IBB was applied.

Zoom In Zoom Out Reset image size

Figures 17(a)–(c) illustrate the procedures of the linear IBB boundary condition of Pan et al (2006) in a one-dimensional setting. As shown in figure 17(a), when the boundary x_W is located at the middle between nodes x_A and x_B (q = |x_A - x_W |/|x_A - x_B | = 1/2 ), the particle at node x_A travels and collides with the wall at x_W and reverses its momentum (the collision process completes instantly), then travels back to x_A . Thus, the incoming distribution function is simply equal to the corresponding outgoing one with the opposite momentum. This is the so-called halfway bounce-back (HWBB) condition and accurate in this case. The IBB condition thus generalizes the HWBB condition. When q < 1/2 as in figure 17(b), the particle at x_A should end up at x_C after the streaming–collision cycle. Likewise, the particle starting from x_C ends up at x_A after the streaming–collision cycle. This is accomplished by interpolating the distribution function for x_C before the streaming–collision process with the wall takes place. Similarly, when q > 1/2 as in figure 17(c), the incoming distribution function can be obtained by using the outgoing one located at x_C after the streaming–collision interaction with the wall takes place and the distribution function values at nearby nodes x_D and x_E . Thus, the linear IBB formulae for $f_{\rm{L}} (x_A ,t_{n + 1} ) = \hat f_{\rm{R}} (x_C ,t_n )$ may be written as

$$\begin{equation} f_{\rm{L}} (x_A ,t_{n + 1} ) = \left\{ \begin{array}{l} (1 - 2q)f_{\rm{R}} (x_A ,t_{n + 1} ) + 2q\hat f_{\rm{R}} (x_A ,t_n ),\quad q < 1/2, \\ \left( {1 - \frac{1}{{2q}}} \right)\hat f_{\rm{L}} (x_A ,t_n ) + \frac{1}{{2q}}\hat f_{\rm{R}} (x_A ,t_n ),\quad q \geqslant 1/2, \\ \end{array} \right. \end{equation} \tag{ 49 }$$

where f and $\hat f$ respectively denote the post- and pre-collision states of the distribution function. The subscripts L and R indicate left- and right-bound directions, respectively.

Zoom In Zoom Out Reset image size

In order to represent the curved surfaces, an interpolation scheme is also applied to the diffuse scattering boundary condition, which is called the interpolated diffuse scattering (IDS) boundary condition. As schematically shown in figure 17(d), the distribution functions can be interpolated as

$$\begin{equation} f_{\rm{R}} (x_{\rm W} ,t_{n + 1} ) = qf_{\rm{R}} (x_A ,t_n ) + (1 - q)f_{\rm{R}} (x_D ,t_n ), \end{equation} \tag{ 50a }$$

$$\begin{equation} f_{\rm{L}} (x_A ,t_{n + 1} ) = qf_{\rm{L}} (x_{\rm W} ,t_n ) + (1 - q)f_{\rm{L}} (x_{\rm W} ,t_{n + 1} ). \end{equation} \tag{ 50b }$$

The D2Q9 MRT LBM simulations with the IDBBC, the regularization and Ψ 1 were carried out on the 2D grids at the corresponding Kn = 0.14 to that of the MD result. The applied computational grid was 96 × 96 . Figure 18 compares the velocity profiles in the triangular-cylinder flow. As clearly seen in figure 18(a), without the regularization procedure, the distributions along the slant face at x/H = 0.43 and 0.5 are very different and show unphysical shapes, particularly by the 1.5th- and second-order slip-flow models. Indeed, with the higher-order slip-flow model, the slip velocity at the peak of the triangle tends to be oscillating as noticed in the profiles at x/H = 0.5 . As discussed before, such oscillations can be dealiased by introducing the regularization procedure. With the regularization process, the oscillations and discontinuity in the profiles disappear, as shown in figure 18(b). (As discussed so far, all the slip-flow models produce virtually the same velocity profiles as in figures 18(b) and (c).) It can also be observed that the discrepancy between the MRT LBM and the MD results is up to 20% of U_b at the most difficult position to predict.

Zoom In Zoom Out Reset image size

To confirm the tendencies observed in the 2D flows, 3D nano-mesh flows were discussed. Figure 19 illustrates a schematic view of the nano-mesh flow of Suga et al (2011). Since the periodic boundary condition was applied to each direction, this flow domain represented a kind of mesh geometry whose porosity was φ = 0.88 . This geometrical structure of the nano-mesh is similar to that of a model structure of cellular foams (e.g. Du Plessis and Masliyah 1988, Fourie and Du Plessis 2002). The force corresponding to a mean pressure gradient was applied to the x-direction. The characteristic Knudsen numbers from the MD simulation results were Kn = λ /d_h = 0.1 and 0.2 where d_h was the mean hydraulic diameter of the nano-mesh: d_h = 0.43H. The LBGK LBM simulations with the DSBC, the regularization and Ψ 1 were carried out on a regular uniform grid of 100 × 100 × 100 at the corresponding Kn. (After a series of grid dependency tests using grids of 40 × 40 × 40 , 100 × 100 × 100 and 140 × 140 × 140 , the results by the grid of 100 × 100 × 100 were confirmed to be enough for grid-independent solutions.)

Zoom In Zoom Out Reset image size

Figures 20(a)–(f) respectively compare the distribution profiles of the velocity U in the x -direction at sections A, B and C indicated in figure 19 at Kn = 0.1 and 0.2. (The figures are re-plotted with the additional data of the MRT LBM with the first-order slip-flow DBBC for this paper.) All the sections are in the symmetry plane of y/H = 0.5 as in figure 19. It is clearly shown that both the results of the LBGK simulations by the D3Q19 and the D3Q39 models agree with those of the MD simulation. Most of the LBM results are in the fluctuation range of the MD results although the LBM simulations predict peak velocities lower than those by the MD simulation at z/H = 0 and 1.0 in figures 20(a), (e) and (f). The discrepancies of the results are up to 20–30% of U_b . These discrepancies come from the predicted slip velocities at the wall surfaces. Since the LBM predicts slightly larger slip velocities at some locations as shown in figure 20(c) for example, the corresponding peak velocity distribution tends to be lower. It is also clear that the difference between the two discrete velocity models is small. It can also be seen that the difference between the results of the LBGK LBMs and the D3Q19 MRT LBM with the DBBC is marginal. Indeed, the results of the D3Q19 MRT LBM almost perfectly overlap with those of the D3Q39 LBGK LBM.

Zoom In Zoom Out Reset image size

Since it is confirmed that the profiles of the convenient D3Q19 model are comparable with those of the elaborate D3Q39 model, the superiority in the accuracy of the higher-order discrete velocity models may not be very important for simulating general flow characteristics in complex 3D flows as in the 2D flows discussed so far. We conclude that retaining terms of the third order in the Hermite expansion and using a Gauss–Hermite quadrature of the third order to obtain the discrete velocity values are not significantly important for predicting general flow profiles even in 3 D nano-porous media at least at 0.1 ⩽Kn ⩽ 0.2 .

Suga and Ito (2011) considered the channel flows whose walls were 3D bumpy shapes to discuss the performance of the LBM for the three-dimensionally curved surfaces. Figure 21 shows the computational domain and schematic views of the 3D bumpy channel flow considered. As shown in figures 21(a) and (b), the matrices of hemisphere bumps were placed on the top and bottom walls in a staggered fashion. The radius of the hemisphere was R = 0.4H where H is the channel height and the pitch of the bumps. Periodic boundary conditions were applied to the streamwise (x ) and spanwise (z ) directions. To drive the fluid, the force was applied to the x-direction only. The estimated Knudsen number in the MD simulation was Kn = 0.1. The corresponding MRT LBM simulations were carried out on the 3D grids using the linear IDBBC with the regularization and Ψ 1 by the D3Q19 model. The grid size used was 100 × 100 × 100 .

Zoom In Zoom Out Reset image size

Figure 22 compares the velocity profiles at the symmetry plane, z/H = 0.5 (figure 22(a)), and an off-center plane, z/H = 0.25 (figures 22(b)–(d)). As discussed in the 2D flow cases, with the regularization procedure, it is clearly seen that all the slip-flow models show the same predictive tendency. Generally, the agreement between the MRT LBM and the MD simulation results is very satisfactory. The maximum discrepancy between the results is indeed up to 15% of the bulk velocity U_b . However, unlike the 2D flows, even with the regularization, somehow kinky profiles still can be seen near walls at some sections such as x/H = 0.1 and 0.9 in the symmetry plane (figure 22(a)). Also weak oscillations appear near the bottom wall at x/H = 0.3 and 0.7 in the plane of z/H = 0.25 (figures 22(b)–(d)). Since the flow field and characteristics are much more complicated than those of the 2D test flows, it may not be easy to dampen such near-wall kinks or oscillations by the regularization procedure even at Kn = 0.1. It is worth re-considering the interpolation scheme for the diffusive bounce-back boundary condition used in the study because its accuracy of the curvature representation in the tangential direction is first order.

Zoom In Zoom Out Reset image size