Three-Dimensional Simulations of Anisotropic Slip Microflows Using the Discrete Unified Gas Kinetic Scheme

Abstract

The specific objective of the present work study is to propose an anisotropic slip boundary condition for three-dimensional (3D) simulations with adjustable streamwise and spanwise slip length by the discrete unified gas kinetic scheme (DUGKS). The present boundary condition is proposed based on the assumption of nonlinear velocity profiles near the wall instead of linear velocity profiles in a unidirectional steady flow. Moreover, a 3D corner boundary condition is introduced to the DUGKS to reduce the singularities. Numerical tests validate the effectiveness of the present method, which is more accurate than the bounce-back and specular reflection slip boundary condition in the lattice Boltzmann method. It is of significance to study the lid-driven cavity flow due to its applications and its capability in exhibiting important phenomena. Then, the present work explores, for the first time, the effects of anisotropic slip on the two-sided orthogonal oscillating micro-lid-driven cavity flow by adopting the present method. This work will generate fresh insight into the effects of anisotropic slip on the 3D flow in a two-sided orthogonal oscillating micro-lid-driven cavity. Some findings are obtained: The oscillating velocity of the wall has a weaker influence on the normal velocity component than on the tangential velocity component. In most cases, large slip length has a more significant influence on velocity profiles than small slip length. Compared with pure slip in both top and bottom walls, anisotropic slip on the top wall has a greater influence on flow, increasing the 3D mixing of flow. In short, the influence of slip on the flow field depends not only on slip length but also on the relative direction of the wall motion and the slip velocity. The findings can help in better understanding the anisotropic slip effect on the unsteady microflow and the design of microdevices.

1. Introduction

Surface characteristics play a critical role in designing and fabricating microfluidic devices. Superhydrophobicity is an important aspect of surface characteristics, which can significantly control flow and reduce drag [1,2,3,4,5,6,7,8,9,10]. Based on the knowledge of fluid mechanics, the no-slip boundary condition is valid at the solid–liquid interface. However, the liquid slip velocity is observed on the superhydrophobic surface. Unlike gas slip caused by Knudsen effects, liquid slip is modelled with two strategies [11]: introducing the force that repels water in the multiphase system; and modelling the slip boundary condition. For the former strategy, the forces are not well understood and determined. Existing research recognizes the critical role played by the latter strategy; therefore, it is emphasized in the present work.

The appropriate boundary condition is an important area in the simulation of fluid flows [12,13]. Recently, researchers have shown an increased interest in slip boundary conditions. For example, Min and Kim [14,15] directly modelled the hydrophobic wall with Navier’s slip boundary condition by direct numerical simulation (DNS) based on the macroscopic Navier–Stokes equations. However, the Navier’s slip boundary condition cannot be introduced to the lattice-Boltzmann-based methods, because the lattice-Boltzmann-based methods track the particle distribution function instead of the macroscopic variables. The lattice Boltzmann method (LBM) is appropriate for simulating mesoscopic physics that are hard to describe macroscopically. It has the advantages of simple algorithms, natural parallelization, and saving computing costs for simulating microflow [16], which proves to be a promising tool. Slip boundary conditions have been proposed and employed in the LBM, such as bounce-back and specular reflection (BSR) [11,17], discrete full diffusive and specular reflection (DSR) [18,19], discrete full diffusive and bounce-back (DBB) [20,21], and tangential momentum accommodation coefficient (TMAC) scheme [22]. Coupled with Navier’s slip model [23], liquid slip can be characterized and adjusted by slip length. With assumption of the 2D unidirectional flow at a lower Reynolds number, the relation between the slip length and combination parameter of the coupled schemes [24,25,26] is determined. However, it is not valid in three-dimensional (3D) flows and turbulent flows. Moreover, the slip should be considered in streamwise and spanwise directions [1]. Existing studies about rice leaves have found that the anisotropic groove-like microstructures can lead to the anisotropic slip behavior of water droplets on the surfaces [27,28,29]. More recent studies have considered superhydrophobic surfaces with randomly distributed textures in streamwise and spanwise directions [8,30,31,32,33,34,35], indicating that the slip is anisotropic. For example, both spanwise and streamwise slip lengths have been measured on a randomly textured superhydrophobic surface [35]. Therefore, the present study considers the anisotropic slip in streamwise and spanwise directions for the 3D system based on the nonlinear velocity profiles near the wall, which is close to the actual situation.

Moreover, Guo et al. proposed a new finite-volume scheme named discrete unified gas kinetic scheme (DUGKS) based on the lattice Boltzmann equation [36]. The DUGKS has also been proved to be a promising numerical tool to simulate fluid flow [37,38,39,40,41,42], such as 3D turbulent channel flow. It is found that the DUGKS, even with a coarse nonuniform mesh, is overall better than the LBM [37]. So, the numerical simulation will be performed by the DUGKS in this work. Up to now, no attention has been paid to the 3D slip boundary condition in the DUGKS. Therefore, in the present study, the new slip boundary condition is proposed for the DUGKS instead of the LBM. The DUGKS is performed to simulate 3D flows with the proposed slip boundary condition. It is hoped that we provide a superior method to describe and characterize anisotropic slip on superhydrophobic surfaces.

To study the effect of anisotropic slip on flows in benchmark geometries, the present method is applied to the lid-driven cavity flow. It is an important problem in the field of fluid mechanics due to its applications, such as cooling of electronic gadgets, oil extraction, design of heat exchangers, solar ponds, acoustic liner, float glass productions, insulation materials, multiscreen gadgets for nuclear reactors, coating, food processing, crystal growth, etc. [43,44,45,46,47,48,49]. Moreover, it has capability in exhibiting important phenomena such as eddies, secondary flows, instabilities, transition, and turbulence [50]. The liquid slip flow in two-sided, orthogonal, oscillating, micro-lid-driven cavities has largely been oversighted. Recently, two-sided motion [51,52,53,54] and oscillatory flows in the cavity have also caught the necessary attention, except single-sided steady flows. The purposes of the literature cornering oscillatory flows in lid-driven cavities include:

1. To test and validate numerical solvers, such as least-squares finite element methods [55], Taylor-series-expansion- and least-squares-based lattice Boltzmann methods [56] and conservative level-set methods [57]. 2. To understand industrial applications, such as surface viscometer [58] and optimization of fluid mixing [59,60]. 3. To understand the flow characteristics, such as the single-sided oscillatory rarefied gas flows inside two- and three-dimensional cavities [61,62,63], and two-sided oscillating flows in two-dimensional lid-driven cavities [64].

The studies mentioned above have been solely restricted to the no-slip flow; the liquid slip flow in lid-driven cavities have largely been oversighted. Slip should be carefully considered in the design of micro-devices with moving parts. This work will generate fresh insight into the effects of anisotropic slip on the 3D flow in a two-sided oscillating micro-lid-driven cavity. In this work, for the first time, there are mainly two types of slip distribution: (1) pure streamwise slip emerges on both top and bottom wall surfaces; (2) both streamwise and spanwise slips emerge on the top wall surface. In a 3D global coordinate system, the unit velocity vector of the moving top and bottom walls is (1,0,0) and (0,1,0), respectively. The motion direction of the top and bottom walls is orthogonal. To the best of the authors’ knowledge, no such study has been conducted before.

It is hoped that the present study will provide a superior method and contribute to a deeper understanding of the anisotropic slip. The paper is organized as follows. In Section 2, the D3Q19 lattice model and DUGKS are introduced, and the new slip boundary condition and corner boundary condition for the D3Q19 lattice model are derived theoretically. In Section 3, numerical validation is performed by simulating the 3D microchannel flow. Numerical results of 3D flow in a two-sided, oscillating, lid-driven cavity are discussed in Section 4. Section 5 gives the conclusions.

2. Numerical Methods

2.1. D3Q19 Lattice Model

D3Q19 lattice model [65] is adopted in this work. As shown in Figure 1, there are 19 discrete velocities in the D3Q19 lattice model, including one rest velocity (α = 0) and 18 non-rest velocities (α = 1, ..., 18).

As shown in

Table 1. The velocity set for the D3Q19 lattice model.

$\mathbf{Velocities} {\mathit{\xi }}_{\mathit{\alpha }}$	Number	$\mathbf{Length} |{\mathit{\xi }}_{\mathit{\alpha }}|$	$\mathbf{Weight} {\mathit{W}}_{\mathit{\alpha }}$
(0,0,0)c	1	0	1/3
(±1, 0, 0)c, (0, ±1, 0)c, (0, 0, ±1)c	6	1	1/18
(±1, ±1, 0)c, (±1, 0, ±1)c, (0, ±1, ±1)c	12	$\sqrt{2}$	1/36

, the velocity set includes the velocities $\left\{{\xi }_{\alpha }\right\}$ and the corresponding weights $\left\{W_{\alpha }\right\}$. The speed of the sound is $c=\sqrt{3RT}=1$.

2.2. Discrete Unified gas Kinetic Scheme

The discrete Boltzmann equation with the Bhatnagar–Gross–Krook (BGK) collision model [66] is the governing equation of the DUGKS:

$$\frac{\partial f_{\alpha }}{\partial t}+{\mathit{\xi }}_{\alpha }\cdot \nabla f_{\alpha }=\mathsf{\Omega }\equiv \frac{f_{\alpha }{}^{\mathbf{e}q}-f_{\alpha }}{\tau }$$

(1)

It is assumed that fluid particles move with the velocity ${\mathit{\xi }}_{\alpha }$ at position $\mathbf{x}$ and time t, so the velocity distribution function is $f_{\alpha }=f_{\alpha }\left(\mathbf{x},{\mathit{\xi }}_{\alpha },t\right)$. Ω and τ represent the collision term and relaxation time, respectively.

$f_{\alpha }{}^{\mathbf{e}q}$ represents the Maxwellian equilibrium distribution function, which is approximated by Taylor expansion around zero particle velocity at low Mach number:

$$f_{\alpha }{}^{\mathbf{e}q}=W_{\alpha }\rho \left[1+\frac{{\mathit{\xi }}_{\alpha }\cdot u}{c_s^2}+\frac{{\left({\mathit{\xi }}_{\alpha }\cdot u\right)}^2}{2c_s^4}-\frac{{|u|}^2}{2c_s^2}\right], c_s^2=RT$$

(2)

The velocities $\left\{{\mathit{\xi }}_{\alpha }\right\}$ and the corresponding weights $\left\{W_{\alpha }\right\}$ are presented in Table 1.

The computational domain is divided into cuboid cells V_i,j,k = Δx_iΔy_jΔz_k centered at x_i,j,k = (x_i, y_j, z_k) in the DUGKS. As a new finite volume scheme, the volumes-averaged values of the distribution function and collision term need to be computed. For example, the volumes-averaged value of the distribution function $f_{\alpha }^n({\mathbf{x}}_{i,j,k})$ is computed as,

$$f_{\alpha }^n({\mathbf{x}}_{i,j,k})=\frac{1}{|V_{i,j,k}|}{\displaystyle \underset{V_{i,j,k}}{\int }{f}_{\alpha }\left(\mathbf{x},t^n\right)}dV$$

(3)

In the DUGKS, the governing equation needs to be integrated on each cell, and the time step Δt is assumed to be constant. Equation (1) is integrated on a cell V_i,j,k centered at x_i,j,k from time $t^n=n\mathsf{\Delta }t$ to time $t^{n+1}=\left(n+1\right)\mathsf{\Delta }t$, and the evolution of the distribution is advanced from tⁿ to tⁿ⁺¹ as,

$$f_{\alpha }^{n+1}-f_{\alpha }^n=-\frac{\mathsf{\Delta }t}{|V_{i,j,k}|}{ℱ}_{\alpha }^{n+1/2}+\frac{\mathsf{\Delta }t}{2}\left[{\mathsf{\Omega }}_{\alpha }^n+{\mathsf{\Omega }}_{\alpha }^{n+1}\right]$$

(4)

The scalar variable ${ℱ}_{\alpha }^{n+1/2}$ represents the flux across the cell interface,

$$\begin{array}{l}{ℱ}_{\alpha }^{n+1/2}\left({\mathbf{x}}_{i,j,k}\right)={\displaystyle {\int }_{\partial V_{i,j,k}}\left({\mathit{\xi }}_{\alpha }\cdot \mathit{n}\right)}{f}_{\alpha }^{n+1/2}({\mathbf{x}}_b)dS=\\ =[f_{\alpha }^{n+1/2}({\mathbf{x}}_{i,j,k}+0.5\mathsf{\Delta }{x}_i{\mathbf{e}}_x)-f_{\alpha }^{n+1/2}({\mathbf{x}}_{i,j,k}-0.5\mathsf{\Delta }{x}_i{\mathbf{e}}_x)]{\xi }_{\alpha ,x}\mathsf{\Delta }{y}_j\mathsf{\Delta }{z}_k\\ +[f_{\alpha }^{n+1/2}({\mathbf{x}}_{i,j,k}+0.5\mathsf{\Delta }{y}_j{\mathbf{e}}_y)-f_{\alpha }^{n+1/2}({\mathbf{x}}_{i,j,k}-0.5\mathsf{\Delta }{y}_j{\mathbf{e}}_y)]{\xi }_{\alpha ,y}\mathsf{\Delta }{x}_i\mathsf{\Delta }{z}_k\\ +[f_{\alpha }^{n+1/2}({\mathbf{x}}_{i,j,k}+0.5\mathsf{\Delta }{z}_k{\mathbf{e}}_z)-f_{\alpha }^{n+1/2}({\mathbf{x}}_{i,j,k}-0.5\mathsf{\Delta }{z}_k{\mathbf{e}}_z)]{\xi }_{\alpha ,z}\mathsf{\Delta }{x}_i\mathsf{\Delta }{y}_j\end{array}$$

(5)

where $f_{\alpha }^{n+1/2}({\mathbf{x}}_b)$ represents the distribution at the cell interface x_b at the time t^n+1/2 = tⁿ + h (h = Δt/2), and e_x, e_y, and e_z are unit vectors in x, y, and z directions, respectively.

For clarity, new distribution functions are introduced:

$${\tilde{f}}_{\alpha }^n\equiv f_{\alpha }^n-\frac{\mathsf{\Delta }t}{2}\mathsf{\Omega }\left(f_{\alpha }^n\right),{\tilde{f}}_{\alpha }^{+,n}\equiv f_{\alpha }^n+\frac{\mathsf{\Delta }t}{2}\mathsf{\Omega }\left(f_{\alpha }^n\right)$$

(6)

The collision term can be expanded in the BGK collision model, and Equation (6) can be converted to the following equations:

$$f_{\alpha ,j}^n=\frac{2\tau }{2\tau +\mathsf{\Delta }t}{\tilde{f}}_{\alpha ,j}^n+\frac{\mathsf{\Delta }t}{2\tau +\mathsf{\Delta }t}{f}_{\alpha ,j}^{\mathbf{e}q,n} , {\tilde{f}}_{\alpha ,j}^{+,n}=\frac{2\tau -\mathsf{\Delta }t}{2\tau +\mathsf{\Delta }t}{\tilde{f}}_{\alpha ,j}^n+\frac{2\mathsf{\Delta }t}{2\tau +\mathsf{\Delta }t}{f}_{\alpha ,j}^{\mathbf{e}q,n} .$$

(7)

The evolution equation of DUGKS from tⁿ to tⁿ⁺¹ is simplified as:

$${\tilde{f}}_{\alpha ,j}^{n+1}={\tilde{f}}_{\alpha ,j}^{+,n}-\frac{\mathsf{\Delta }t}{|V_{i,j,k}|}{ℱ}_{\alpha ,j}^{n+1/2}$$

(8)

Based on the conservation of mass, momentum, the density ρ, and velocity u can be computed from ${\tilde{f}}_{\alpha }$:

$${\rho }^n={\displaystyle \sum _{\alpha }{\tilde{f}}_{\alpha }^n},{\rho }^n{\mathbf{u}}^n={\displaystyle \sum _{\alpha }{\mathit{\xi }}_{\alpha }{\tilde{f}}_{\alpha }^n}$$

(9)

All other forms of the distribution function can be expressed in terms of ${\tilde{f}}_{\alpha }$ and $f_{\alpha }{}^{\mathbf{e}q}$. So, the distribution function ${\tilde{f}}_{\alpha }$ is mainly computed in the code.

The critical step is to evaluate the interface flux ${ℱ}_{\alpha ,j}^{n+1/2}$. The midpoint integral formula is employed to evaluate ${ℱ}_{\alpha ,j}^{n+1/2}$, due to its easy implementation and fast computation. For DUGKS with higher-order accuracy, more intermediate time steps need to be selected; for example, the flux at the cell interface at t * = t_n + ∆t/6 and t * = t_n+ 3∆t/4 need calculating.

In the present study, Equation (1) is integrated within a half time step (h = Δt/2) along the characteristic line with the endpoint (x_b) located at the cell interface, and the following formula is obtained:

$$f_{\alpha }^{n+1/2}\left({\mathbf{x}}_b\right)-f_{\alpha }^n\left({\mathbf{x}}_b-h{\mathit{\xi }}_{\alpha }\right)=\frac{h}{2}\left[\mathsf{\Omega }\left(f_{\alpha }^{n+1/2}\left({\mathbf{x}}_b\right)\right)+\mathsf{\Omega }\left(f_{\alpha }^n\left({\mathbf{x}}_b-h{\mathit{\xi }}_{\alpha }\right)\right)\right]$$

(10)

Similarly, new distribution functions are introduced and can be computed by expanding the collision term:

$${\overline{f}}_{\alpha }^{n+1/2}\left({\mathbf{x}}_b\right)\equiv f_{\alpha }^{n+1/2}\left({\mathbf{x}}_b\right)-\frac{h}{2}\mathsf{\Omega }\left(f_{\alpha }^{n+1/2}\left({\mathbf{x}}_b\right)\right)=\frac{2\tau +h}{2\tau }{f}_{\alpha }^{n+1/2}\left({\mathbf{x}}_b\right)-f_{\alpha }^{\mathbf{e}q,n+1/2}\left({\mathbf{x}}_b\right)$$

(11)

$$\begin{array}{c}{\overline{f}}_{\alpha }^{+,n}\left({\mathbf{x}}_b-h{\mathit{\xi }}_{\alpha }\right)\equiv f_{\alpha }^n\left({\mathbf{x}}_b-h{\mathit{\xi }}_{\alpha }\right)+\frac{h}{2}\mathsf{\Omega }\left(f_{\alpha }^n\left({\mathbf{x}}_b-h{\mathit{\xi }}_{\alpha }\right)\right),\\ {\overline{f}}_{\alpha }^{+,n}=\frac{2\tau -h}{2\tau +h}{\overline{f}}_{\alpha }^n+\frac{2h}{2\tau +h}{f}_{\alpha }^{\mathbf{e}q,n}.\end{array}$$

(12)

With new distribution functions, Equation (10) is converted to the following equation:

$${\overline{f}}_{\alpha }^{n+1/2}\left({\mathbf{x}}_b\right)={\overline{f}}_{\alpha }^{+,n}\left({\mathbf{x}}_b-h{\mathit{\xi }}_{\alpha }\right)$$

(13)

With the Taylor expansion around the endpoint (x_b) located at the cell interface, the right term of Equation (13) can be approximated as:

$${\overline{f}}_{\alpha }^{+,n}\left({\mathbf{x}}_b-h{\mathit{\xi }}_{\alpha }\right)={\overline{f}}_{\alpha }^{+,n}\left({\mathbf{x}}_b\right)-h{\mathit{\xi }}_{\alpha }\cdot \nabla {\overline{f}}_{\alpha }^{+,n}\left({\mathbf{x}}_b\right)$$

(14)

where ${\overline{f}}_{\alpha }^{+,n}\left({\mathbf{x}}_b\right)$ and its gradients $\nabla {\overline{f}}_{\alpha }^{+,n}\left({\mathbf{x}}_b\right)$ can be approximated by the linear interpolations. In x-direction, e.g.,

$$\begin{array}{l}\frac{\partial {\overline{f}}_{\alpha }^{+,n}\left({\mathbf{x}}_{i,j,k}+0.5\mathsf{\Delta }{x}_i{\mathbf{e}}_x\right)}{\partial x}\approx \frac{{\overline{f}}_{\alpha }^{+,n}\left({\mathbf{x}}_{i+1,j,k}\right)-{\overline{f}}_{\alpha }^{+,n}\left({\mathbf{x}}_{i,j,k}\right)}{(\mathsf{\Delta }{x}_i+\mathsf{\Delta }{x}_{i+1})/2},\\ {\overline{f}}_{\alpha }^{+,n}\left({\mathbf{x}}_{i,j,k}+0.5\mathsf{\Delta }{x}_i{\mathbf{e}}_x\right)\approx {\overline{f}}_{\alpha }^{+,n}\left({\mathbf{x}}_{i,j,k}\right)+\frac{\partial {\overline{f}}_{\alpha }^{+,n}\left({\mathbf{x}}_{i,j,k}+0.5\mathsf{\Delta }{x}_i{\mathbf{e}}_x\right)}{\partial x}\frac{\mathsf{\Delta }{x}_i}{2},\end{array}$$

(15)

The distribution function ${\overline{f}}_{\alpha }^{+,n}$ can be computed from ${\tilde{f}}_{\alpha }$, as follows,

$${\overline{f}}_{\alpha }^{+,n}=\frac{2\tau -h}{2\tau +\mathsf{\Delta }t}{\tilde{f}}_{\alpha }^n+\frac{3h}{2\tau +\mathsf{\Delta }t}{f}_{\alpha }^{\mathbf{e}q,n}$$

(16)

Then, we obtain the function ${\overline{f}}_{\alpha }^{n+1/2}\left({\mathbf{x}}_b\right)$. The density and velocity at the cell interface can also be evaluated, which can be used for the equilibrium distribution function $f_{\alpha }^{\mathbf{e}q,n+1/2}\left({\mathbf{x}}_b\right)$,

$${{\rho }^{n+1/2}|}_{{\mathbf{x}}_b}={\displaystyle \sum _{\alpha }{\overline{f}}_{\alpha }^{n+1/2}}\left({\mathbf{x}}_b\right),{{\left(\rho \mathbf{u}\right)}^{{}^{n+1/2}}|}_{{\mathbf{x}}_b}={\displaystyle \sum _{\alpha }{\mathit{\xi }}_{\alpha }{\overline{f}}_{\alpha }^{n+1/2}}\left({\mathbf{x}}_b\right)$$

(17)

Finally, the flux ${ℱ}_{\alpha ,j}^{n+1/2}$ is evaluated according to Equation (5) after the distribution function $f_{\alpha }^{n+1/2}$ at the cell interface is determined by Equation (11). The tracked distribution function ${\tilde{f}}_{\alpha }$ can be updated to the next time step after the flux is obtained.

Particularly, the following equation will be used in the DUGKS,

$${\tilde{f}}_{\alpha }^{+,n}=\frac{4}{3}{\overline{f}}_{\alpha }^{+,n}-\frac{1}{3}{\tilde{f}}_{\alpha }^n$$

(18)

For the present DUGKS, the relaxation time τ is computed from τ = μ/p, where μ is the dynamic viscosity coefficient. p (p = ρRT) is the pressure, where R is the specific gas constant. In the following simulations, the temperature T is constant in the isothermal flow with $c_s{}^2=RT=1/3$. The time step Δt is determined by the Courant–Friedrichs–Lewy (CFL) condition [67], which is independent of the relaxation time τ for all flow regimes.

2.3. The Present Slip Boundary Condition

It can be seen intuitively that, considering the actual case with anisotropic slip, slip boundary conditions derived by adopting two-dimensional unidirectional flow are not valid. Therefore, a new anisotropic slip boundary condition is proposed in 3D flows. In this work, the anisotropic slip boundary condition is characterized and constructed by adjusting the relative magnitude of the streamwise and spanwise slip lengths. It is noted that x, y, and z denote the streamwise, spanwise, and wall-normal directions, respectively. It is noted that the DUGKS tracks the distribution function ${\tilde{f}}_{\alpha }$, unlike the LBM.

Considering the impermeable wall boundary (U_Wz = 0), the unknown distributions are obtained by the specular reflection (${\tilde{f}}_{\alpha }^{sr}$) and the stress exerted by the wall (${\tilde{f}}_{\alpha }^w$).

$${\tilde{f}}_{\alpha }={\tilde{f}}_{\alpha }^{sr}+{\tilde{f}}_{\alpha }^w({\mathit{\xi }}_{\alpha }\cdot \mathbf{n}>0)$$

(19)

As shown in Figure 2, the unknown distributions are ${\tilde{f}}_4,{\tilde{f}}_8,{\tilde{f}}_9,{\tilde{f}}_{12},{\tilde{f}}_{14}$.

The specular reflection ${\tilde{f}}_{\alpha }^{sr}$ can be obtained by:

$${\tilde{f}}_4^{sr}={\tilde{f}}_3,{\tilde{f}}_8^{sr}={\tilde{f}}_{10},{\tilde{f}}_9^{sr}={\tilde{f}}_7,{\tilde{f}}_{12}^{sr}={\tilde{f}}_{13},{\tilde{f}}_{14}^{sr}={\tilde{f}}_{11}$$

(20)

With ${\tilde{f}}_{\alpha }^{sr}$ determined, $U_{SR}$ can be obtained by:

$$\rho ={\displaystyle \sum _{\alpha }{\tilde{f}}_{\alpha }}$$

(21)

$$\rho U_{SR}={\displaystyle \sum _{{\mathit{\xi }}_{\alpha }\cdot \mathbf{n}\le 0}{\left({\tilde{f}}_{\alpha }{\mathit{\xi }}_{\alpha }\right)}_{||}+}{\displaystyle \sum _{{\mathit{\xi }}_{\alpha }\cdot \mathbf{n}>0}{\left({\tilde{f}}_{\alpha }^{sr}{\mathit{\xi }}_{\alpha }\right)}_{||}}$$

(22)

The stress ${\tilde{f}}_{\alpha }^w$ contributes to the change in tangential momentum caused by the wall surface, as shown in the following equations:

$${\sigma }^{\prime }\rho \left(U_W-U_{SR}\right)={\displaystyle \sum _{{\mathit{\xi }}_{\alpha }\cdot \mathbf{n}>0}{\left({\tilde{f}}_{\alpha }^w{\mathit{\xi }}_{\alpha }\right)}_{||}}$$

(23)

$$0={\displaystyle \sum _{{\mathit{\xi }}_{\alpha }\cdot \mathbf{n}>0}{\left({\tilde{f}}_{\alpha }^w{\mathit{\xi }}_{\alpha }\right)}_{\perp }}$$

(24)

$$0={\tilde{f}}_{\alpha }^w(\mathrm{for} \mathrm{normal} \mathrm{direction})$$

(25)

$U_W$ and $U_{SR}$ are the tangential velocity of the wall and the average tangential velocity under the specular reflection by the impermeable boundary, respectively. ${\sigma }^{\prime }$ represents a modified tangential momentum accommodation coefficient. $n$ is the normal direction of the wall toward the fluid field. The subscripts “$||$” and “$\perp$” represent the tangential and normal directions, respectively. The sum of normal parts is zero, ensuring that the function ${\tilde{f}}_{\alpha }^w$ only changes the tangential momentum. It is also shown that the calculation of density is not determined by ${\displaystyle \sum {\tilde{f}}_{\alpha }^w}$.

The following relations can be obtained for the case in Figure 2.

$$\mathrm{x} \text{-} \mathrm{direction}: \rho U_{SR\mathrm{x}}={\tilde{f}}_1-{\tilde{f}}_2+2({\tilde{f}}_7-{\tilde{f}}_{10})$$

(26)

$$\mathrm{y} \text{-} \mathrm{direction}: \rho U_{SR\mathrm{y}}={\tilde{f}}_8-{\tilde{f}}_5+2({\tilde{f}}_{13}-{\tilde{f}}_{11})$$

(27)

The momentum can change due to the shear stress imposed on the wall:

$$\mathrm{x} \text{-} \mathrm{direction}:{{\sigma }^{\prime }}_{\mathrm{x}}\rho \left(U_{W\mathrm{x}}-U_{SR\mathrm{x}}\right)={\displaystyle \sum _{{\mathit{\xi }}_{\alpha }\cdot \mathbf{n}>0}{\left({\tilde{f}}_{\alpha }^w{\mathit{\xi }}_{\alpha }\right)}_{\mathrm{x}}}={\tilde{f}}_9^w-{\tilde{f}}_8^w$$

(28)

$$\mathrm{y} \text{-} \mathrm{direction}:{{\sigma }^{\prime }}_{\mathrm{y}}\rho \left(U_{W\mathrm{y}}-U_{SR\mathrm{y}}\right)={\displaystyle \sum _{{\mathit{\xi }}_{\alpha }\cdot \mathbf{n}>0}{\left({\tilde{f}}_{\alpha }^w{\mathit{\xi }}_{\alpha }\right)}_{\mathrm{y}}}={\tilde{f}}_{12}^w-{\tilde{f}}_{14}^w$$

(29)

For positive normal direction of the wall (i.e., +z direction):

$$0={\displaystyle \sum _{{\mathit{\xi }}_{\alpha }\cdot \mathbf{n}>0}{({\tilde{f}}_{\alpha }^w{\mathit{\xi }}_{\alpha })}_z}={\tilde{f}}_4^w+{\tilde{f}}_8^w+{\tilde{f}}_9^w+{\tilde{f}}_{12}^w+{\tilde{f}}_{14}^w$$

(30)

$$0={\tilde{f}}_4^w$$

(31)

Based on the Maxwell equilibrium distribution function and Ref. [22], the density can be calculated,

$$\rho ={\tilde{f}}_0+{\tilde{f}}_1+{\tilde{f}}_2+{\tilde{f}}_5+{\tilde{f}}_6+{\tilde{f}}_{15}+{\tilde{f}}_{16}+{\tilde{f}}_{17}+{\tilde{f}}_{18}+2({\tilde{f}}_3+{\tilde{f}}_7+{\tilde{f}}_{10}+{\tilde{f}}_{11}+{\tilde{f}}_{13})$$

(32)

Then, ${\tilde{f}}_8^w, {\tilde{f}}_9^w, {\tilde{f}}_{12}^w, {\tilde{f}}_{14}^w$ and $U_{SR\mathrm{x}}, U_{SR\mathrm{y}}$ are still unknown. ${{\sigma }^{\prime }}_{\mathrm{x}},{{\sigma }^{\prime }}_{\mathrm{y}}$ are the manually adjustable parameters, which are related to the streamwise and spanwise slip lengths, respectively.

There are five known equations in the system:

$$\{\begin{array}{c}\begin{array}{c}\rho U_{SR\mathrm{x}}={\tilde{f}}_1-{\tilde{f}}_2+2({\tilde{f}}_7-{\tilde{f}}_{10})\\ \rho U_{SR\mathrm{y}}={\tilde{f}}_8-{\tilde{f}}_5+2({\tilde{f}}_{13}-{\tilde{f}}_{11})\\ {{\sigma }^{\prime }}_{\mathrm{x}}\rho \left(U_{W\mathrm{x}}-U_{SR\mathrm{x}}\right)={\displaystyle \sum _{{\mathit{\xi }}_{\alpha }\cdot \mathbf{n}>0}{\left({\tilde{f}}_i^w{\mathit{\xi }}_{\alpha }\right)}_{\mathrm{x}}}={\tilde{f}}_9^w-{\tilde{f}}_8^w\\ {{\sigma }^{\prime }}_{\mathrm{y}}\rho \left(U_{W\mathrm{y}}-U_{SR\mathrm{y}}\right)={\displaystyle \sum _{{\mathit{\xi }}_{\alpha }\cdot \mathbf{n}>0}{\left({\tilde{f}}_i^w{\mathit{\xi }}_{\alpha }\right)}_{\mathrm{y}}}={\tilde{f}}_{12}^w-{\tilde{f}}_{14}^w\end{array}\\ 0={\tilde{f}}_8^w+{\tilde{f}}_9^w+{\tilde{f}}_{12}^w+{\tilde{f}}_{14}^w\end{array}$$

(33)

To make the system closed, the hypothesis is proposed:

$${\tilde{f}}_8^w=-{\tilde{f}}_9^w or {\tilde{f}}_{12}^w=-{\tilde{f}}_{14}^w$$

(34)

In summary, a new slip boundary condition is proposed for the upper horizontal wall boundary in 3D:

$$\begin{array}{l}{\tilde{f}}_4={\tilde{f}}_3\\ {\tilde{f}}_8={\tilde{f}}_{10}-\frac{1}{2}{{\sigma }^{\prime }}_{\mathrm{x}}\left\{\rho U_{W\mathrm{x}}-({\tilde{f}}_1-{\tilde{f}}_2+2({\tilde{f}}_7-{\tilde{f}}_{10}))\right\}\\ {\tilde{f}}_9={\tilde{f}}_7+\frac{1}{2}{{\sigma }^{\prime }}_{\mathrm{x}}\left\{\rho U_{W\mathrm{x}}-({\tilde{f}}_1-{\tilde{f}}_2+2({\tilde{f}}_7-{\tilde{f}}_{10}))\right\}\\ {\tilde{f}}_{12}={\tilde{f}}_{13}+\frac{1}{2}{{\sigma }^{\prime }}_{\mathrm{y}}\left\{\rho U_{W\mathrm{y}}-({\tilde{f}}_6-{\tilde{f}}_5+2({\tilde{f}}_{13}-{\tilde{f}}_{11}))\right\}\\ {\tilde{f}}_{14}={\tilde{f}}_{11}-\frac{1}{2}{{\sigma }^{\prime }}_{\mathrm{y}}\left\{\rho U_{W\mathrm{y}}-({\tilde{f}}_6-{\tilde{f}}_5+2({\tilde{f}}_{13}-{\tilde{f}}_{11}))\right\}\end{array}$$

(35)

With the external forcing term, the local velocities $u$ are computed by,

$$u=\frac{1}{\rho }{\displaystyle \sum _i{\mathit{\xi }}_{\alpha }{\tilde{f}}_{\alpha }}+\frac{\mathsf{\Delta }t\cdot \overrightarrow{a}}{2}$$

(36)

To eliminate the numerical slip due to force tangential to the wall, the external forcing term is introduced to the new slip boundary condition:

$$\begin{array}{l}{\tilde{f}}_4={\tilde{f}}_3\\ {\tilde{f}}_8={\tilde{f}}_{10}-\frac{1}{2}{{\sigma }^{\prime }}_{\mathrm{x}}\left\{\rho (U_{W\mathrm{x}}-0.5\mathsf{\Delta }t{a}_x)-({\tilde{f}}_1-{\tilde{f}}_2+2({\tilde{f}}_7-{\tilde{f}}_{10}))\right\}\\ {\tilde{f}}_9={\tilde{f}}_7+\frac{1}{2}{{\sigma }^{\prime }}_{\mathrm{x}}\left\{\rho (U_{W\mathrm{x}}-0.5\mathsf{\Delta }t{a}_x)-({\tilde{f}}_1-{\tilde{f}}_2+2({\tilde{f}}_7-{\tilde{f}}_{10}))\right\}\\ {\tilde{f}}_{12}={\tilde{f}}_{13}+\frac{1}{2}{{\sigma }^{\prime }}_{\mathrm{y}}\left\{\rho (U_{W\mathrm{y}}-0.5\mathsf{\Delta }t{a}_{\mathrm{y}})-({\tilde{f}}_6-{\tilde{f}}_5+2({\tilde{f}}_{13}-{\tilde{f}}_{11}))\right\}\\ {\tilde{f}}_{14}={\tilde{f}}_{11}-\frac{1}{2}{{\sigma }^{\prime }}_{\mathrm{y}}\left\{\rho (U_{W\mathrm{y}}-0.5\mathsf{\Delta }t{a}_{\mathrm{y}})-({\tilde{f}}_6-{\tilde{f}}_5+2({\tilde{f}}_{13}-{\tilde{f}}_{11}))\right\}\end{array}$$

(37)

Similar manipulations can be applied to the lower wall and side walls boundary.

2.4. Relation between the Combination Parameters and Slip Lengths

Then, the relation between combination parameters (${{\sigma }^{\prime }}_{\mathrm{x}},{{\sigma }^{\prime }}_{\mathrm{y}}$) and slip lengths (b_x, b_y) is deduced to implement the new slip boundary condition. Previous research on the derivation of the relation is studied by taking the two-dimensional unidirectional steady flow as an example, which is expressed as [24,25,26]:

$$\rho =const,u_y=0,a_y=0,\frac{\partial \varphi }{\partial x}=0,\frac{\partial \varphi }{\partial t}=0$$

(38)

where $\varphi$ denotes flow variable, such as the velocity or density.

In this work, it is assumed that the anisotropic slip includes two components in streamwise and spanwise directions. Take the liquid slip on a horizontal plane (perpendicular to the z-axis) as an example. The slip length includes two components, b_x, b_y in the x and y directions, respectively. With the assumption of anisotropic slip, the simulation and characterization of the slip effect on a superhydrophobic surface can match the actual situation.

The upper wall in 3D is employed to derive the relationship between the parameters ${{\sigma }^{\prime }}_{\mathrm{x}},{{\sigma }^{\prime }}_{\mathrm{y}}$ and the slip lengths b_x, b_y.

The flow near the wall satisfies the continuity equation:

$$\frac{\partial \rho }{\partial t}+\frac{\partial \rho u_x}{\partial x}+\frac{\partial \rho u_y}{\partial y}+\frac{\partial \rho u_z}{\partial z}=0$$

(39)

With the assumption of no density change in the incompressible flow, the continuity equation can be written as,

$$\frac{\partial u_x}{\partial x}+\frac{\partial u_y}{\partial y}+\frac{\partial u_z}{\partial z}=0$$

(40)

In the 3D directional flows, the velocity distribution near the wall is assumed to be satisfied the following equations:

$$\frac{\partial u_x}{\partial y}=0, \frac{\partial u_y}{\partial x}=0, u_z=0$$

(41)

As shown in Figure 2, the local coordinate system $(\overline{x},\overline{y},\overline{z})$ in 3D is established in the lattice unit, where a node on the wall is served as the origin. The $\overline{x}-\overline{y}$ plane is located on the wall. The $\overline{z}$ direction is normal along the wall. In the local coordinate system, the boundary of the upper wall is located at the plane $\overline{z}=1$, and the solid is located in the region ($\overline{z}<1$), e.g., the plane $\overline{z}=0$.

As shown in Ref. [68], the measured velocity profiles across the channel with a parabolic fit are observed and recorded. Therefore, in this work, it is assumed that the nonlinear velocity profiles near the wall are parabolic in the $\overline{z}$ direction, conforming to the quadratic term fitting. The velocity profiles can be linear when the quadratic coefficient is 0. Then, the function of velocity profiles near the wall will be simplified as follows:

$$u_x(\overline{x},\overline{z})=u_x(\overline{x})+{\alpha }_1{\overline{z}}^2+{\beta }_1\overline{z}+e_1$$

(42)

$$u_y(\overline{y},\overline{z})=u_y(\overline{y})+{\alpha }_2{\overline{z}}^2+{\beta }_2\overline{z}+e_2$$

(43)

The acceleration can be approximated by relations as follows [69]:

$$a_x\approx -\nu \frac{{\partial }^2{u}_x(\overline{x},\overline{z})}{\partial {\overline{z}}^2},a_y\approx -\nu \frac{{\partial }^2{u}_y(\overline{y},\overline{z})}{\partial {\overline{z}}^2}.$$

(44)

Then, the coefficients can be obtained by the acceleration, α₁ = −0.5 a_x/ν, α₂ = −0.5 a_y/ν.

The slip velocity can be expressed as:

$$u_{sx}=u_x(\overline{x},\overline{z}){|}_{\overline{z}=1}-U_{Wx},u_{sy}=u_y(\overline{y},\overline{z}){|}_{\overline{z}=1}-U_{Wy}.$$

(45)

The slip velocity on a wall is characterized in the form of a Navier slip boundary condition in both the streamwise direction and the spanwise direction [35]:

$$u_{sx}=b_x\frac{\partial u_x}{\partial z}{|}_{\mathrm{wall}} , u_{sy}=b_y\frac{\partial u_y}{\partial z}{|}_{\mathrm{wall}}$$

(46)

Considering Equations (42) and (43), (46) can be written as:

$$u_{sx}=b_x\frac{\partial u_x}{\partial \overline{z}}{|}_{\overline{z}=0}=b_x{\beta }_1,u_{sy}=b_y\frac{\partial u_y}{\partial \overline{z}}{|}_{\overline{z}=0}=b_y{\beta }_2.$$

(47)

With the Taylor expansion around the $\overline{z}$ = 1 in the local coordinate system, $u_x(\overline{z})$ and $u_y(\overline{z})$ can be approximated,

$$u_x(\overline{z}){|}_{\overline{z}=2}=u_x(\overline{z}){|}_{\overline{z}=1}+\mathsf{\Delta }\overline{z}\frac{\partial u_x(\overline{z})}{\partial \overline{z}}{|}_{\overline{z}=1}+\frac{\mathsf{\Delta }{\overline{z}}^2}{2}(\frac{{\partial }^2{u}_x(\overline{z})}{\partial {\overline{z}}^2}){|}_{\overline{z}=1}+\mathrm{O}(\mathsf{\Delta }{\overline{z}}^3)$$

(48)

$$u_y(\overline{z}){|}_{\overline{z}=2}=u_y(\overline{z}){|}_{\overline{z}=1}+\mathsf{\Delta }\overline{z}\frac{\partial u_y(\overline{z})}{\partial \overline{z}}{|}_{\overline{z}=1}+\frac{\mathsf{\Delta }{\overline{z}}^2}{2}(\frac{{\partial }^2{u}_y(\overline{z})}{\partial {\overline{z}}^2}){|}_{\overline{z}=1}+\mathrm{O}(\mathsf{\Delta }{\overline{z}}^3)$$

(49)

With the assumption of parabolic velocity profiles near the wall, substitute the above equations, and the following equations can be obtained,

$$\begin{array}{l}{u}_x(\overline{x},\overline{z}){|}_{\overline{z}=2}-u_x(\overline{x},\overline{z}){|}_{\overline{z}=1}=(\mathsf{\Delta }{\overline{z}}^2+2\mathsf{\Delta }\overline{z}){\alpha }_1+\mathsf{\Delta }\overline{z}{\beta }_1,\\ u_y(\overline{y},\overline{z}){|}_{\overline{z}=2}-u_x(\overline{y},\overline{z}){|}_{\overline{z}=1}=(\mathsf{\Delta }{\overline{z}}^2+2\mathsf{\Delta }\overline{z}){\alpha }_2+\mathsf{\Delta }\overline{z}{\beta }_2.\end{array}$$

(50)

Then, the relations between the coefficients are given as,

$${\beta }_1=\frac{(3-\mathsf{\Delta }{\overline{z}}^2-2\mathsf{\Delta }\overline{z})}{\mathsf{\Delta }\overline{z}-1+\mathrm{O}(\mathsf{\Delta }\overline{z})}{\alpha }_1$$

(51)

$${\beta }_2=\frac{(3-\mathsf{\Delta }{\overline{z}}^2-2\mathsf{\Delta }\overline{z})}{\mathsf{\Delta }\overline{z}-1+\mathrm{O}(\mathsf{\Delta }\overline{z})}{\alpha }_2$$

(52)

For D3Q19 lattice model, with $\rho u={\displaystyle \sum _{\alpha =0}^{18}{\mathit{\xi }}_{\alpha }{f}_{\alpha }}$ known, the local velocities $u_x,u_y$ at $\overline{z}$ = 1 and 2 can be calculated as,

$$\begin{array}{l}\rho u_x{|}_{\overline{z}=1}={\tilde{f}}_1^1-{\tilde{f}}_2^1+{\tilde{f}}_7^1-{\tilde{f}}_{10}^1+{\tilde{f}}_9^1-{\tilde{f}}_8^1,\\ \rho u_x{|}_{\overline{z}=2}={\tilde{f}}_1^2-{\tilde{f}}_2^2+{\tilde{f}}_7^2-{\tilde{f}}_{10}^2+{\tilde{f}}_9^2-{\tilde{f}}_8^2,\end{array}$$

(53)

$$\begin{array}{l}\rho u_y{|}_{\overline{z}=1}={\tilde{f}}_6^1-{\tilde{f}}_5^1+{\tilde{f}}_{13}^1-{\tilde{f}}_{11}^1+{\tilde{f}}_{12}^1-{\tilde{f}}_{14}^1,\\ \rho u_y{|}_{\overline{z}=2}={\tilde{f}}_6^2-{\tilde{f}}_5^2+{\tilde{f}}_{13}^2-{\tilde{f}}_{11}^2+{\tilde{f}}_{12}^2-{\tilde{f}}_{14}^2.\end{array}$$

(54)

Combined with the proposed slip boundary condition, the following relations can be obtained:

$$\begin{array}{l}{\tilde{f}}_9-{\tilde{f}}_8={\tilde{f}}_7-{\tilde{f}}_{10}+{{\sigma }^{\prime }}_{\mathrm{x}}\left\{\rho U_{W\mathrm{x}}-[{\tilde{f}}_1-{\tilde{f}}_2+2({\tilde{f}}_7-{\tilde{f}}_{10})]\right\},\\ {\tilde{f}}_{12}-{\tilde{f}}_{14}={\tilde{f}}_{13}-{\tilde{f}}_{11}+{{\sigma }^{\prime }}_{\mathrm{y}}\left\{\rho U_{W\mathrm{y}}-[{\tilde{f}}_6-{\tilde{f}}_5+2({\tilde{f}}_{13}-{\tilde{f}}_{11})]\right\}.\end{array}$$

(55)

Then,

$$\begin{array}{l}\rho u_x{|}_{\overline{z}=1}=(1-{{\sigma }^{\prime }}_{\mathrm{x}})[{\tilde{f}}_{{}_1}^1-{\tilde{f}}_{{}_2}^1+2({\tilde{f}}_{{}_7}^1-{\tilde{f}}_{{}_{10}}^1)]+{{\sigma }^{\prime }}_{\mathrm{x}}\rho U_{W\mathrm{x}},\\ \rho u_x{|}_{\overline{z}=2}=(1-{{\sigma }^{\prime }}_{\mathrm{x}})[{\tilde{f}}_{{}_1}^2-{\tilde{f}}_{{}_2}^2+2({\tilde{f}}_{{}_7}^2-{\tilde{f}}_{{}_{10}}^2)]+{{\sigma }^{\prime }}_{\mathrm{x}}\rho U_{W\mathrm{x}},\\ \rho u_y{|}_{\overline{z}=1}=(1-{{\sigma }^{\prime }}_{\mathrm{y}})[{\tilde{f}}_{{}_6}^1-{\tilde{f}}_{{}_5}^1+2({\tilde{f}}_{{}_{13}}^1-{\tilde{f}}_{{}_{11}}^1)]+{{\sigma }^{\prime }}_{\mathrm{y}}\rho U_{W\mathrm{y}},\\ \rho u_y{|}_{\overline{z}=2}=(1-{{\sigma }^{\prime }}_{\mathrm{y}})[{\tilde{f}}_{{}_6}^2-{\tilde{f}}_{{}_5}^2+2({\tilde{f}}_{{}_{13}}^2-{\tilde{f}}_{{}_{11}}^2)]+{{\sigma }^{\prime }}_{\mathrm{y}}\rho U_{W\mathrm{y}}.\end{array}$$

(56)

The difference value between velocity at $\overline{z}$ = 1 and $\overline{z}$ = 2 can be written as:

$$\begin{array}{l}{u}_x{|}_{\overline{z}=2}-u_x{|}_{\overline{z}=1}=\frac{1}{\rho }(1-{{\sigma }^{\prime }}_{\mathrm{x}})\{({\tilde{f}}_{{}_1}^2-{\tilde{f}}_{{}_2}^2)-({\tilde{f}}_{{}_1}^1-{\tilde{f}}_{{}_2}^1)+2[({\tilde{f}}_{{}_7}^2-{\tilde{f}}_{{}_{10}}^2)-({\tilde{f}}_{{}_7}^1-{\tilde{f}}_{{}_{10}}^1)]\},\\ u_y{|}_{\overline{z}=2}-u_y{|}_{\overline{z}=1}=\frac{1}{\rho }(1-{{\sigma }^{\prime }}_{\mathrm{y}})\{({\tilde{f}}_{{}_6}^2-{\tilde{f}}_{{}_5}^2)-({\tilde{f}}_{{}_6}^1-{\tilde{f}}_{{}_5}^1)+2[({\tilde{f}}_{{}_{13}}^2-{\tilde{f}}_{{}_{11}}^2)+({\tilde{f}}_{{}_{13}}^1-{\tilde{f}}_{{}_{11}}^1)]\}.\end{array}$$

(57)

Inspired by Guo et al. [69], the collision and streaming rule in the LBM is adopted to establish the relationship between velocities near the wall. Considering the collision and streaming rule of LBE with BGK operator [70], the following relations can be obtained:

$$\mathrm{x} \text{-} \mathrm{direction}:{\tilde{f}}_7^1-{\tilde{f}}_{10}^1={\tilde{\overline{f}}}_7^2-{\tilde{\overline{f}}}_{10}^2,{\tilde{f}}_9^2-{\tilde{f}}_8^2={\tilde{\overline{f}}}_9^1-{\tilde{\overline{f}}}_8^1$$

(58)

$$\mathrm{y} \text{-} \mathrm{direction}:{\tilde{f}}_{13}^1-{\tilde{f}}_{11}^1={\tilde{\overline{f}}}_{13}^2-{\tilde{\overline{f}}}_{11}^2,{\tilde{f}}_{12}^2-{\tilde{f}}_{14}^2={\tilde{\overline{f}}}_{12}^1-{\tilde{\overline{f}}}_{14}^1$$

(59)

where $\tilde{\overline{f}}$ denotes the tracked distribution function in the collision.

Then, Equation (57) can be simplified as follows:

$$u_x{|}_{\overline{z}=2}-u_x{|}_{\overline{z}=1}=\frac{3{{\sigma }^{\prime }}_{\mathrm{x}}}{\tau (1-{{\sigma }^{\prime }}_{\mathrm{x}})}(u_x{|}_{\overline{z}=1}-U_{W\mathrm{x}})$$

(60)

$$u_y{|}_{\overline{z}=2}-u_y{|}_{\overline{z}=1}=\frac{3{{\sigma }^{\prime }}_{\mathrm{y}}}{\tau (1-{{\sigma }^{\prime }}_{\mathrm{y}})}(u_y{|}_{\overline{z}=1}-U_{W\mathrm{y}})$$

(61)

The slip velocity could be calculated as:

$$u_{sx}=b_x\frac{\partial u_x}{\partial \overline{z}}{|}_{\overline{z}=0}=b_x{\beta }_1=u_x{|}_{\overline{z}=1}-U_{W\mathrm{x}}=\frac{\tau (1-{{\sigma }^{\prime }}_{\mathrm{x}})}{3{{\sigma }^{\prime }}_{\mathrm{x}}}(u_x{|}_{\overline{z}=2}-u_x{|}_{\overline{z}=1})$$

(62)

$$u_{sy}=b_x\frac{\partial u_y}{\partial \overline{z}}{|}_{\overline{z}=0}=b_y{\beta }_1=u_y{|}_{\overline{z}=1}-U_{W\mathrm{y}}=\frac{\tau (1-{{\sigma }^{\prime }}_{\mathrm{y}})}{3{{\sigma }^{\prime }}_{\mathrm{y}}}(u_y{|}_{\overline{z}=2}-u_y{|}_{\overline{z}=1})$$

(63)

Finally, the relations between the slip lengths and parameters are obtained:

$$b_x=\frac{\tau (1-{{\sigma }^{\prime }}_{\mathrm{x}})}{3{{\sigma }^{\prime }}_{\mathrm{x}}}(u_x{|}_{\overline{z}=2}-u_x{|}_{\overline{z}=1})/{\beta }_1=(\mathsf{\Delta }{\overline{z}}^2+2\mathsf{\Delta }\overline{z})\frac{{\alpha }_1}{{\beta }_1}+\mathsf{\Delta }\overline{z}$$

(64)

$$b_y=\frac{\tau (1-{{\sigma }^{\prime }}_{\mathrm{y}})}{3{{\sigma }^{\prime }}_{\mathrm{y}}}(u_y{|}_{\overline{z}=2}-u_y{|}_{\overline{z}=1})/{\beta }_1=(\mathsf{\Delta }{\overline{z}}^2+2\mathsf{\Delta }\overline{z})\frac{{\alpha }_2}{{\beta }_2}+\mathsf{\Delta }\overline{z}$$

(65)

With the known values of ${\alpha }_1/{\beta }_1,{\alpha }_2/{\beta }_2$, Equations (64) and (65) can be simplified,

$${{\sigma }^{\prime }}_{\mathrm{x}}=\frac{1}{1+\frac{3b_x}{A\tau }},{{\sigma }^{\prime }}_{\mathrm{y}}=\frac{1}{1+\frac{3b_y}{A\tau }}.$$

(66)

where $A$ denotes the correction coefficient and is determined by:

$$A=(\mathsf{\Delta }{\overline{z}}^2+\mathsf{\Delta }\overline{z})\frac{\mathsf{\Delta }\overline{z}-1+O(\mathsf{\Delta }{\overline{z}}^3)}{3-\mathsf{\Delta }{\overline{z}}^2-2\mathsf{\Delta }\overline{z}}+\mathsf{\Delta }\overline{z}$$

(67)

where $\mathsf{\Delta }\overline{z}$ is the lattice grid spacing.

For the upper horizontal wall boundary in 3D, a new slip boundary condition is significantly determined by Equations (35), (66), and (67). Similar derivation and operation can be applied to the lower wall and side walls.

2.5. Corner Boundary Condition

The above discussion on boundary conditions focuses on straight surfaces. Considering the singularity, the treatment of corners should not be ignored in numerical simulations of the flow, such as the lid-driven cavity flow. Although corners account for only a few nodes, these corners should not be underestimated because the particle can stream in the fluid domain, which has influences on the performance of the numerical simulation. One single point may contaminate the numerical solution everywhere [71]. One of the earliest systematic approaches to treating corners in DUGKS was proposed by Guo et al. [72]. However, this approach is limited to 2D implementations. In this work, to reduce the singularities and improve the performance of numerical simulation, an approach to treating the corner boundary condition is proposed for the DUGKS based on the D3Q19 model with 19 independent moments [73].

$$\begin{array}{l}{0}^{\mathrm{th}}:\rho ={\displaystyle \sum _i{f}_i}; 1^{\mathrm{st}}:\rho u_{\alpha }={\displaystyle \sum _i{f}_i}{\xi }_{i\alpha }; 2^{\mathrm{nd}}:{\Pi }_{\alpha \beta }={\displaystyle \sum _i{f}_i}{\xi }_{i\alpha }{\xi }_{i\beta };\\ 3^{\mathrm{rd}}:Q_{\alpha \beta \gamma }={\displaystyle \sum _i{f}_i}{\xi }_{i\alpha }{\xi }_{i\beta }{\xi }_{i\gamma }; { 4}^{\mathrm{th}}:S_{\alpha \beta \gamma \delta }={\displaystyle \sum _i{f}_i}{\xi }_{i\alpha }{\xi }_{i\beta }{\xi }_{i\gamma }{\xi }_{i\delta }.\end{array}$$

(68)

The 0th moment has 1 equation, the 1st moment has 3 equations, the 2nd moment has 6 equations, the 3rd moment has 6 equations, and the 4th moment has 3 equations. In the 3D domain, there are 12 unknowns at every corner. So, 12 linearly independent moments are required. For the D3Q19 model, as shown in Figure 1b, the unknown functions are $f_1,f_3,f_6,f_7,f_9,f_{10},f_{11},f_{12},f_{13},f_{15},f_{16},f_{17}$, considering the low-south-west corner. We select the momenta $\rho u_x,\rho u_y,\rho u_z$, the momentum fluxes and shear stresses ${\Pi }_{xx},{\Pi }_{yy},{\Pi }_{zz},{\Pi }_{xy},{\Pi }_{xz},{\Pi }_{yz}$, and three higher-order moments $Q_{xxy},Q_{yyz},Q_{xzz}$ as 12 linearly independent moments.

$$\mathrm{rank}\left[\begin{array}{cc}{\tilde{f}}_1+{\tilde{f}}_9-{\tilde{f}}_{10}+{\tilde{f}}_{15}-{\tilde{f}}_{16}+{\tilde{f}}_{17}& \rho u_x\\ {\tilde{f}}_6-{\tilde{f}}_{11}+{\tilde{f}}_{12}+{\tilde{f}}_{13}-{\tilde{f}}_{15}+{\tilde{f}}_{16}+{\tilde{f}}_{17}& \rho u_y\\ {\tilde{f}}_3+{\tilde{f}}_7-{\tilde{f}}_9+{\tilde{f}}_{10}+{\tilde{f}}_{11}-{\tilde{f}}_{12}+{\tilde{f}}_{13}& \rho u_z\\ {\tilde{f}}_1+{\tilde{f}}_9+{\tilde{f}}_{10}+{\tilde{f}}_{15}+{\tilde{f}}_{16}+{\tilde{f}}_{17}& {\Pi }_{xx}\\ {\tilde{f}}_6+{\tilde{f}}_{11}+{\tilde{f}}_{12}+{\tilde{f}}_{13}+{\tilde{f}}_{15}+{\tilde{f}}_{16}+{\tilde{f}}_{17}& {\Pi }_{yy}\\ {\tilde{f}}_3+{\tilde{f}}_7+{\tilde{f}}_9+{\tilde{f}}_{10}+{\tilde{f}}_{11}+{\tilde{f}}_{12}+{\tilde{f}}_{13}& {\Pi }_{zz}\\ {\tilde{f}}_{17}-{\tilde{f}}_{15}-{\tilde{f}}_{16}& {\Pi }_{xy}\\ {\tilde{f}}_7-{\tilde{f}}_9-{\tilde{f}}_{10}& {\Pi }_{xz}\\ {\tilde{f}}_{13}-{\tilde{f}}_{11}-{\tilde{f}}_{12}& {\Pi }_{yz}\\ {\tilde{f}}_{17}-{\tilde{f}}_{15}+{\tilde{f}}_{16}& Q_{xxy}\\ {\tilde{f}}_7+{\tilde{f}}_9-{\tilde{f}}_{10}& Q_{xzz}\\ {\tilde{f}}_{13}+{\tilde{f}}_{11}-{\tilde{f}}_{12}& Q_{yyz}\end{array}\right]=12$$

(69)

The moments can be approximated by the Chapman–Enskog expansion as follows:

$${\Pi }_{\alpha \beta }={\Pi }_{\alpha \beta }^{(0)}+\epsilon {\Pi }_{\alpha \beta }^{(1)}+O({\epsilon }^2)\approx {\Pi }_{\alpha \beta }^{(0)},Q_{\alpha \beta }=Q_{\alpha \beta }^{(0)}+\epsilon Q_{\alpha \beta }^{(1)}+O({\epsilon }^2)\approx Q_{\alpha \beta }^{(0)}$$

(70)

where the equilibrium part of the momentum flux tensor (${\Pi }_{\alpha \beta }^{(0)}$) and the third-order moment ($Q_{\alpha \beta }^{(0)}$) can be expressed as:

$$\begin{array}{l}{\Pi }_{\alpha \beta }^{(0)}={\displaystyle \sum _i{f}_i^{(0)}}{\xi }_{i\alpha }{\xi }_{i\beta }=\rho u_{\alpha }{u}_{\beta }+\rho c_s^2{\delta }_{\alpha \beta },\\ Q_{\alpha \beta \gamma }^{(0)}={\displaystyle \sum _i{f}_i^{(0)}}{\xi }_{i\alpha }{\xi }_{i\beta }{\xi }_{i\gamma }=\rho c_s^2(u_{\alpha }{\delta }_{\beta \gamma }+u_{\beta }{\delta }_{\alpha \gamma }+u_{\gamma }{\delta }_{\alpha \beta }).\end{array}$$

(71)

The velocity is set to 0 at the corner, and some terms are assumed to be negligible. The momentum fluxes and shear stresses ${\Pi }_{xx},{\Pi }_{yy},{\Pi }_{zz},{\Pi }_{xy},{\Pi }_{xz},{\Pi }_{yz}$, and three higher-order moments $Q_{xxy},Q_{yyz},Q_{xzz}$ are written as follows:

$$\begin{array}{l}{\Pi }_{xx}={\Pi }_{yy}={\Pi }_{zz}=\rho c_s^2=\rho /3,\\ {\Pi }_{xy}=\rho u_x{u}_y=0,{\Pi }_{xz}=\rho u_x{u}_z=0,{\Pi }_{yz}=\rho u_y{u}_z=0,\\ Q_{xxy}=\frac{\rho }{3}(u_x{\delta }_{xy}+u_x{\delta }_{xy}+u_y{\delta }_{xx})=\frac{\rho }{3}{u}_y=0,\\ Q_{yyz}=\frac{\rho }{3}(u_y{\delta }_{yz}+u_y{\delta }_{yz}+u_z{\delta }_{yy})=\frac{\rho }{3}{u}_z=0,\\ Q_{xzz}=\frac{\rho }{3}(u_x{\delta }_{zz}+u_x{\delta }_{xz}+u_z{\delta }_{xz})=\frac{\rho }{3}{u}_x=0.\end{array}$$

(72)

The unknown functions are calculated as:

$$\begin{array}{l}{\tilde{f}}_1=-\frac{\rho }{3}+{\tilde{f}}_2+2{\tilde{f}}_5+4{\tilde{f}}_{14}+4{\tilde{f}}_{18}\\ {\tilde{f}}_3=-\frac{\rho }{3}+{\tilde{f}}_4+2{\tilde{f}}_2+4{\tilde{f}}_8+4{\tilde{f}}_{18}\\ {\tilde{f}}_6=-\frac{\rho }{3}+{\tilde{f}}_5+2{\tilde{f}}_4+4{\tilde{f}}_8+4{\tilde{f}}_{14}\\ {\tilde{f}}_7=\frac{\rho }{6}-{\tilde{f}}_2-{\tilde{f}}_8-2{\tilde{f}}_{18}\\ {\tilde{f}}_9={\tilde{f}}_8\\ {\tilde{f}}_{10}=\frac{\rho }{6}-{\tilde{f}}_2-{\tilde{f}}_8-2{\tilde{f}}_{18}\\ {\tilde{f}}_{11}={\tilde{f}}_{14}\\ {\tilde{f}}_{12}=\frac{\rho }{6}-{\tilde{f}}_4-{\tilde{f}}_{14}-2{\tilde{f}}_8\\ {\tilde{f}}_{13}=\frac{\rho }{6}-{\tilde{f}}_4-{\tilde{f}}_{14}-2{\tilde{f}}_8\\ {\tilde{f}}_{15}=\frac{\rho }{6}-{\tilde{f}}_5-{\tilde{f}}_{18}-2{\tilde{f}}_{14}\\ {\tilde{f}}_{16}={\tilde{f}}_{18}\\ {\tilde{f}}_{17}=\frac{\rho }{6}-{\tilde{f}}_5-{\tilde{f}}_{18}-2{\tilde{f}}_{14}\end{array}$$

(73)

The density at the low-south-west corner is calculated as:

$$\rho ={\tilde{f}}_0+2({\tilde{f}}_2+{\tilde{f}}_4+{\tilde{f}}_5)+4({\tilde{f}}_8+{\tilde{f}}_{14}+{\tilde{f}}_{18}).$$

(74)

Similar manipulations can be applied to other corner nodes.

2.6. Algorithm

The updating of ${\tilde{f}}_{\alpha }$ from time $t^n=n\mathsf{\Delta }t$ to time $t^{n+1}=\left(n+1\right)\mathsf{\Delta }t$ in the DUGKS can be performed as the following brief algorithm.

• Initialize the density, velocity, and viscosity. Obtain the values of $f_{\alpha }^{\mathbf{e}q,n}$ and ${\tilde{f}}_{\alpha }^n$ at time t = 0.
• Compute the distribution functions ${\overline{f}}_{\alpha }^{+,n}$ and ${\tilde{f}}_{\alpha }^{+,n}$ using Equations (16) and (18).
• Compute the value of $\nabla {\overline{f}}_{\alpha }^{+,n}\left({\mathbf{x}}_b\right)$ and ${\overline{f}}_{\alpha }^{+,n}\left({\mathbf{x}}_b\right)$ using Equation (15).
• Compute the distribution function ${\overline{f}}_{\alpha }^{n+1/2}\left({\mathbf{x}}_b\right)$ using Equations (14) and (13).
• Get the macro values of density and velocity using Equation (17). Compute the equilibrium distribution function $f_{\alpha }^{\mathbf{e}q,n+1/2}\left({\mathbf{x}}_b\right)$.
• Compute the distribution function $f_{\alpha }^{n+1/2}\left({\mathbf{x}}_b\right)$ using Equation (13). Obtain the flux ${ℱ}_{\alpha }^{n+1/2}$ by Equation (5).
• For the unknown distribution functions at the boundary or corner, the boundary conditions are employed, such as Equations (35) or (73).
• Update the distribution function ${\tilde{f}}_{\alpha }^{n+1}$ using Equation (8). Obtain the macro values of density and velocity.
• Repeat steps (2)–(8) until the convergence criterion is satisfied.

In C++ DUGKS computer code, the convergence criterion for attaining the steady-state solution is ${\displaystyle \sum |u\left(t\right)-u\left(t-1000\mathsf{\Delta }t\right)|}/{\displaystyle \sum |u\left(t\right)|}<{10}^{-6}$, where $u\left(t\right)$ represents the velocity in the flow field.

3. Numerical Validation

The flow in a 3D channel is a fundamental case in science and engineering. Only a few references on the anisotropic slip boundary condition are available for comparison, so the flow in the 3D channel is selected as the case for numerical validation.

3.1. Comparison with Single-Component Lattice Boltzmann Simulation

In Ref. [11], the slip boundary condition is modelled by combining the bounce-back and specular reflection (BSR) scheme using the single-component lattice Boltzmann method. The relevant parameters in the present simulation remain the same as those in Ref. [11]. As shown in Figure 3, the microchannel’s length, width, and height are 600 μm, 300 μm, and 30 μm, respectively. With the grid convergence study, the spatial discretization with resolution 400 × 200 × 20 (X, Y, and Z directions, respectively) is selected for the subsequent numerical simulations. The inlet and outlet along the X-direction adopt the periodic boundary condition. The remaining four walls adopt the no-slip/slip boundary conditions. For the no-slip case, the bounce-back scheme in LBM is used without treating the corner in Ref. [11], and the bounce-back scheme in DUGKS is used with the corner boundary condition in the present work. For the slip case, the BSR scheme is employed in Ref. [11], and the present method is employed in this work.

In the no-slip case, the present results agree well with the exact solution [74] and experimental result [9], as shown in Figure 4a. It is observed that the present results agree a little better with the exact solution than the BSR scheme in Ref. [11], which shows that the 3D corner boundary condition improves the accuracy.

In the slip case, the present results are closer to the experimental results [9] than those in Ref. [11], as shown in Figure 4b. The present method is more accurate than the BSR scheme in Ref. [11], which may be partly explained by the case that the BSR scheme can generate numerical slip, but the present method with the external force term can eliminate the numerical slip.

3.2. Comparison with Direct Numerical Simulation

In Ref. [14], the value of the streamwise slip length is set to equal the spanwise slip length. Furthermore, the case where the value of streamwise slip length is not equal to the spanwise slip length should be considered. With the different values of streamwise and spanwise slip lengths, the effect of anisotropic slip on velocity profiles and drag has been addressed using direct numerical simulations (DNS) of a turbulent channel flow [75]. In Ref. [75], the Navier slip length boundary condition adopts a linear slip length model. In this work, the present boundary condition is related to the second partial derivative of the velocity, with the assumption of nonlinear velocity profiles near the wall.

To test the accuracy of predicting the drag and velocity, the present method is performed in six different cases at Re_τ = u_τ0δ/ν = 180: (1) Case 1: b_x⁺⁰ = 0.1, b_y⁺⁰ = 1; (2) Case 2: b_x⁺⁰ = 0.316, b_y⁺⁰ = 1; (3) Case 3: b_x⁺⁰ = 3.16, b_y⁺⁰ = 1; (4) Case 4: b_x⁺⁰ = 10, b_y⁺⁰ = 1; (5) Case 5: b_x⁺⁰ = 0.631, b_y⁺⁰ = 1; (6) Case 6: b_x⁺⁰ = 2.51, b_y⁺⁰ = 10. δ and ν denote the channel half-height and kinematic viscosity, respectively. u_τ0 denotes the wall shear (friction) velocity in channel flow with no-slip walls. It is noted that the superscript +0 indicates slip length scales given in units of the viscous length scale ν/u_τ0 in the respective no-slip reference case [76].

The numerical results of the present method are compared to the data in Ref. [75]. The mean streamwise velocity profile is shown in Figure 5, and the root-mean-square (rms) velocity fluctuations are shown in Figure 6. As shown in Figure 5 and Figure 6, the present method is also accurate in predicting the velocity profiles in a turbulent channel flow with an anisotropic slip wall. Similar conclusions to those reported by A. Busse and N. D. Sandham [75] are obtained: the streamwise slip length is mainly responsible for determining mean velocity profiles. Streamwise slip length always has a reducing effect on the intensity of the turbulent fluctuations, and the reducing effect will increase with increasing slip length. Finite streamwise slip length can limit the turbulence-intensifying effects of infinite spanwise slip, thereby limiting the adverse effects of spanwise slip.

To investigate the influence of an anisotropic slip on drag, the DUGKS simulations are conducted by adjusting the streamwise and spanwise slip lengths with the present slip boundary condition. The investigated slip length values are selected randomly. The present results are compared with those in Ref. [75].

The percentage change in drag (DD) is defined by DD = (dp/dx–dp/dx|₀)/(dp/dx|₀), where dp/dx and dp/dx|₀ represent the mean streamwise pressure gradient in the present and reference case, respectively. If DD < 0, the drag is reduced. The DD values are obtained from numerical results in the case of Re_τ0 = 180 based on friction velocity u_τ0 in the reference case.

As shown in Figure 7, the dots match well with the lines, indicating that the present method is also accurate in predicting the change in drag. The same trends reported by Min and Kim [14] are recovered: drag is reduced in cases with pure streamwise slip and isotropic slip, but drag is increased in cases with pure spanwise slip.

4. Application to the Two-Sided Orthogonal Oscillating Micro-Lid-Driven Cavity Flow

4.1. Problem Description

The problem is the micro-lid-driven cavity flow with two moving walls, as shown in Figure 8. For two-sided oscillating wall motion, the top and bottom walls move with oscillating velocity U = U₀cos(ωt), where U₀ = 1.0 m/s, oscillating frequency ω = 1875 π/s. The directions of the top and bottom walls are positive X-direction and Y-direction, respectively. The size of the cavity is 400 μm × 400 μm × 400 μm. The cavity is filled with water. The density and kinematic viscosity of water are $\rho =1000 \mathrm{kg}/{\mathrm{m}}^3$ and $\upsilon =1.0\times {10}^{-6}{\mathrm{m}}^2/\mathrm{s}$, respectively. The Reynolds number is calculated as Re = U₀ × L/ν = 1.0 × 0.0004/10⁻⁶ = 400. For simplicity, ω has been dimensionalized as ω’ = ωL/U₀ = 0.75π, and St = ωL²/ν = ω’Re = 300π.

4.2. Convergence Validation Study

To choose an optimal lattice size utilizing fewer computational resources, lattice size convergence is studied. Numerical simulations with two-sided uniformly moving wall motions are performed at Re = 400 using three lattice sizes: 80³ (coarse), 120³ (medium), and 160³ (fine). Figure 9 shows the negligible improvement in the velocity profile on increasing the lattice size from 120³ to 160³.

So, the spatial discretization with resolution 120 × 120 × 120 is used for performing subsequent numerical simulations with two-sided orthogonal oscillating wall motions. To keep Re = 400, U_0lattice = 1.0/15 and the kinematic viscosity is set to ν_lattice = 0.02. The present slip boundary condition is applied to the top and bottom wall, and the corner boundary condition is applied to four corner nodes in the cavity. The four side walls remain at rest with the no-slip boundary condition.

4.3. Results and Discussion

There are mainly two types of slip distribution: pure streamwise slip emerges on both the top and bottom wall surfaces; and both streamwise and spanwise slips emerge on the top wall surface. For comparison, fourteen cases are considered: (a) both top and bottom walls: b_x = b_y = 0; (b) top wall: b_x = 0.1, b_y = 0, bottom wall: b_x = b_y = 0; (c) top wall: b_x = b_y = 0, bottom wall: b_x = 0,b_y = 0.1; (d) top wall: b_x = 0,b_y = 0.1, bottom wall: b_x = b_y = 0; (e) top wall: b_x = 0.1, b_y = 0, bottom wall: b_x = 0, b_y = 0.1; (f) top wall: b_x = 0.1,b_y = 0.1, bottom wall: b_x = b_y = 0; (g) top wall: b_x = 0.1, b_y = 0, bottom wall: b_x = 0,b_y = 0.05; (h) top wall: b_x = 0.1, b_y = 0.05, bottom wall: b_x = b_y = 0; (i) top wall: b_x = 0.05, b_y = 0, bottom wall: b_x = 0, b_y = 0.1; (j) top wall: b_x = 0.05, b_y = 0.1, bottom wall: b_x = b_y = 0; (k) top wall: b_x = 0.1, b_y = 0, bottom wall: b_x = 0, b_y = 0.2; (l) top wall: b_x = 0.1, b_y = 0.2, bottom wall: b_x = b_y = 0; (m) top wall: b_x = 0.2, b_y = 0, bottom wall: b_x = 0, b_y = 0.1; (n) top wall: b_x = 0.2, b_y = 0.1, bottom wall: b_x = b_y = 0. In Ref. [9], their work yields a slip length of approximately 1 μm at the wall coated with hydrophobic octadecyltrichlorosilane for water flow. In the present work, considering actual value of slip length, the values of 0.05, 0.1, and 0.2 in the lattice unit correspond to 0.25 μm, 0.5 μm, and 1 μm, respectively. The symbols b_x and b_y represent the slip length in the X and Y directions, respectively.

The velocity components and vorticity are Important and common parameters to describe the flow. The present work performs a comprehensive parametric study discussing flow velocity components and vorticity. It is noted that U, V, and W are used to denote the velocity component in the X, Y, and Z directions, respectively. The vorticity magnitude is calculated as √{(∂W/∂Y-∂V/∂Z)² + (∂U/∂Z-∂W/∂X)² + (∂V/∂X-∂U/∂Y)²}.

The contours for velocity U, V, and W and the vorticity magnitude of cases (a-n) at t = T, 0.25 T and 0.5 T are shown in Supplementary Material. It is found that nonphysical phenomena and numerical singularity do not exist, which shows that the present method is effective and the present results are credible. Furthermore, the centerline velocity profiles in the X, Y, and Z directions at t = T, 0.25 T and 0.5 T are shown in Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18, Figure 19, Figure 20, Figure 21 and Figure 22.

Figure 10 shows U (the velocity component in the X direction) along the centerline Z-axis at t = T and its magnified view. As shown in Figure 10, the 14 curves are divided into four groups according to the level of b_x (b_x = 0, 0.05, 0.1, and 0.2): group 1: cases (a, c, and d); group 2: cases (i and j); cases (b, e, f, g, h, k, and l); cases (m, and n). The slip length b_x has greater influence on U than b_y. It is found that velocity profiles in each group are very close. Therefore, for b_x at the same level, the existence of b_y on the top or bottom wall has almost no influence on the change in U along the centerline Z-axis. The greater b_x is, the greater the influence it has on the change in U along the centerline Z-axis for z/L < 0.9 (U < 0). For z/L > 0.9, U increases rapidly, and the closer the location is to the top wall, the faster U increases, and the greater the velocity gradient. When the curves intersect, the relative magnitude of the curves changes. The distribution of the intersection points is chaotic, as shown in the red circle in Figure 10. Therefore, when b_x = 0.1, the promotion effect of b_y on increasing U will change with the change in of z.

Figure 11 shows U (the velocity component in the X direction) along the centerline Y-axis at t = T and its magnified view. As shown in Figure 11, U is negative along the centerline Y-axis in cases (c-n), which indicates that the existence of b_y on the top wall or bottom wall results in the negative U along the centerline Y-axis. The anisotropic slip on the top wall with a larger slip length has a greater influence on the negative U along the centerline Y-axis, such as case (l) and case (n). For y/L < 0.52018, the absolute value of U in case (n) is less than that in case (l). For y/L > 0.52018, the absolute value of U in case (l) is less than that in case (n). Maybe there are more intersecting points near y/L = 0.5 and y/L = 0.915 because of the motion of the top and bottom walls and the interaction of the sidewalls and fluid.

Figure 12 shows V (the velocity component in the Y direction) along the centerline Z-axis at t = T and its magnified view. As shown in Figure 12, the results of different slip combinations are relatively close with a slight difference, indicating that anisotropic slip has a minor influence on V along the centerline Z-axis. However, it can still be seen that the slip combination in case (l) has the greatest influence on V along the centerline Z-axis, where the absolute value of negative V is the largest, as shown in the magnified view. Near the top wall, the anisotropic slip in case (l) results in the maximum positive V. For z/L > 0.06554, the intersection points of the curves are mostly evenly distributed along the centerline Z-axis, indicating that the strong or weak influence of slip distribution on V will change at most positions along the centerline Z-axis. It may be caused by the time-dependent motion of both top and bottom walls. Maybe, in this condition, the velocity V is mainly influenced by the disordered flow driven by the walls, and the influence of slip on V is negligible. So, the effect of slip may be greatly reduced in the disordered flow.

Figure 13 shows V (the velocity component in the Y direction) along the centerline X-axis at t = T. As shown in Figure 13, all curves have two troughs. For the no-slip condition, two troughs are located at x/L ≈ 0.07 and x/L ≈ 0.93. The slip condition makes the troughs closer to the center than the no-slip condition. Under the action of anisotropic slip, two peaks appear in the curve. The fluctuation range is large at large slip lengths, such as case (l) and case (n). Compared with pure slip in both top and bottom walls, anisotropic slip on the top wall results in stronger fluctuation disturbance and increases the 3D mixing of flow. The results of different slip combinations are similar near the side walls, and a great difference occurs in the cavity (0.2 < x/L < 0.8). The influence of slip on flow hardly propagates to the side walls, but mostly to the cavity.

Figure 14 shows W (the velocity component in the Z direction) along the centerline X-axis at t = T and its magnified view. As shown in Figure 14, 14 curves were divided into four groups according to the level of b_x: b_x= 0,0.05, 0.1, 0.2. The larger b_x is, the greater the influence it has on the change in W along the centerline X-axis, and the larger the fluctuation range is, the closer the position of the peak or trough is to the center. For the same level of b_x, the existence of b_y on the top or bottom wall has little effect on the change in W along the centerline X-axis. For x/L ≈ 0.5, the direction of W reverses under the interaction of the top and bottom walls. This can be explained by the fact that the liquid inside the cavity flows clockwise.

Figure 15 shows W (the velocity component in the Z direction) along the centerline Y-axis at t = T. As shown in Figure 15, the large slip length will enhance the oscillation phenomenon of W along the centerline Y-axis and promote the increase in fluctuation amplitude, such as case (m) and case (n). This can be explained by the fact that the large slip length significantly increases the moving velocity of the wall. The existence of b_y on the top wall has a stronger effect on enhancing the amplitude of the left peak, and the existence of b_y on the bottom wall has a stronger effect on enhancing the amplitude of the right peak. The direction of b_y intersects with the motion direction of the top wall vertically, enhancing the disturbance of W near the left side wall; the direction of b_y is the same as the motion direction of the bottom wall, enhancing the disturbance of W near the right sidewall. So, the influence of slip on flow also depends on the slip direction and wall motion direction.

Figure 16 shows U (the velocity component in the X direction) along the centerline Z-axis at t = 0.25 T and its magnified view. As shown in Figure 16, the 14 curves are divided into four groups according to the level of b_x (b_x = 0, 0.05, 0.1, 0.2): group 1: cases (d, a, and c); group 2: cases (j and i); cases (l, f, h, b, g, e, and k); cases (n and m). It is found that velocity profiles for U in each group are very close. Therefore, for b_x at the same level, the existence of b_y on the top or bottom wall has almost no influence on the change in U along the centerline Z-axis. So, the direction of slip is also an important consideration. The larger b_x results in a larger peak near the top wall, which has a greater influence on the change in U along the centerline Z-axis. When b_x and b_y are fixed, the anisotropic slip on the top wall has a greater effect on the positive U than the pure slip on the top and bottom walls. This may be explained by the directional inconsistency. When U > 0, there is no intersection point in 14 curves, and with fixed b_x = 1, larger b_y on the top wall results in larger U. With the existence of b_y on the top wall, the increase in velocity is facilitated.

Figure 17 shows U (the velocity component in the X direction) along the centerline Y-axis at t = 0.25 T. As shown in Figure 17, the pure slip on the top and bottom walls makes the curve symmetrical, and the anisotropic slip on the top wall makes the trough closer to the right-side wall. The asymmetry can be caused by the direction of slip normal to the wall motion direction.

Figure 18 shows V (the velocity component in the Y direction) along the centerline Z-axis at t = 0.25 T and its magnified view. With the existence of b_y, the slip velocity component in the Y direction appears on the top wall, and the larger slip length makes the non-negative value of V larger, such as with the case (l). In the magnified view, the order of peak value is k, c, i, e, m, g, a, d, b, j, h, f, l, and n. With the existence of b_y on the bottom wall, the larger b_y results in a larger peak value. The existence of b_y can contribute to influencing the flow. With fixed b_y = 0.1, the larger b_x results in a smaller peak value.

At t = 0.25 T, the top and bottom walls move with oscillating velocity U = U₀cos(ωt) = 0, but the slip drives the flow in the cavity. The trough near the bottom wall is more obvious in Figure 19 than that in Figure 13, owing to the different oscillating velocity. The trend in curves in Figure 20 is similar to that in Figure 14. The order of peak value near the right side wall in Figure 21 is consistent with that in Figure 15, which shows that the oscillating velocity of the top and bottom wall has a weaker influence on W than on U and V. It is found that W is positive along the centerline Y-axis in case (n) at t = 0.25 T. Anisotropic slip with large slip length can result in the disruptive change. Maybe, in this condition, the flow is dominated by anisotropic slip.

As shown in Figure 22, the profiles of U along the centerline Y-axis at t = 0.5 T show no qualitative similarity with those at t = T and t = 0.25 T. The other five types of profiles at t = 0.5 T show some approximative mirror symmetry with those at t = T. Because the top and bottom walls move with oscillating velocity U = U₀cos(ωt), the direction of velocity at t = 0.5 t is opposite to that at t = T.

As shown in Figures S1, S5, S6, S10, S11 and S15 in the Supplementary Materials, the contour of vorticity magnitude is concentrated on the top and bottom walls, owing to shear stress affected by the motion of the top and bottom walls. As shown in

Table 2. Maximum vorticity magnitude at t = T, 0.25 T, 0.5 T.

	t	a	b	c	d	e	f	g	h	i	j	k	l	m	n
Case
T		11.3139	11.3139	11.4836	11.3139	11.4836	11.3139	11.3988	11.3139	11.4836	11.3139	11.6533	11.3139	11.4836	11.3139
0.25 T		1.08977	1.20004	1.0898	1.09042	1.20004	1.20603	1.20004	1.20154	1.07988	1.08963	1.20004	1.22381	2.40009	2.40309
0.5 T		11.3139	11.3139	11.1442	11.3139	11.1442	11.3139	11.2291	11.3139	11.1442	11.3139	10.9745	11.3139	11.1442	11.3139

, the maximum vorticity magnitude at t = T and 0.5 T is about an order of magnitude larger than that at t = 0.25 T, owing to the time-dependent oscillating velocity of the top and bottom walls. Compared to the no-slip case, the maximum vorticity magnitude in slip cases changes very little at t = T and 0.5 T, and the maximum percentage change is 3% in case (k). Compared to the no-slip case, case (m) and case (n) obtain about a 120% increase in the percentage of the maximum vorticity magnitude at t = 0.25 T. It is found that the maximum vorticity magnitude makes no change at t = T and 0.5 T when the anisotropic slip exists on the top wall.

5. Conclusions

The present method is validated by simulating the microchannel flow in 3D. Compared with the reference data, the present method is more accurate than the bounce-back and specular reflection slip boundary condition in LBM in Ref. [11]. The effect of anisotropic slip boundary conditions on turbulent flow is investigated by considering different slip lengths in streamwise and spanwise directions. Good agreement with DNS results shows that the present method is also accurate and stable to simulate fluid slip on 3D hydrophobic microchannel walls in a turbulent flow. The present method is effectively accurate and stable to capture velocity profiles and predict drag changes to study the effect of anisotropic slip. Then, the present method is applied to the two-sided, orthogonal, oscillating, micro-lid-driven cavity flow. Some findings are obtained from the simulation results, which can help in better understanding the anisotropic slip effect on the unsteady microflow and the design of microdevices:

The oscillating velocity of the wall has a weaker influence on W than on U and V. In most cases, large slip length has a more significant influence on velocity profiles than small slip length. However, for V along the centerline Z-axis at t = 0.25 T, the larger streamwise slip length on the top wall results in a smaller peak value with a fixed spanwise slip length. Compared with pure slip in both top and bottom walls, anisotropic slip on the top wall has a greater influence on flow, increasing the 3D mixing of flow. In short, the influence of slip on the flow field depends not only on slip length but also on the relative direction of the wall motion and the slip velocity.