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

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.


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

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 includes the velocities {ξ α } and the corresponding weights {W α }. The speed of the sound is c = √ 3RT = 1.

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: It is assumed that fluid particles move with the velocity ξ α at position x and time t, so the velocity distribution function is f α = f α (x, ξ α , t). Ω and τ represent the collision term and relaxation time, respectively. f α eq represents the Maxwellian equilibrium distribution function, which is approximated by Taylor expansion around zero particle velocity at low Mach number: The velocities {ξ α } and the corresponding weights {W α } 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 n α (x i,j,k ) is computed as, 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∆t to time t n+1 = (n + 1)∆t, and the evolution of the distribution is advanced from t n to t n+1 as, The scalar variable F n+1/2 α represents the flux across the cell interface, α (x i,j,k − 0.5∆x i e x )]ξ α,x ∆y j ∆z k +[ f n+1/2 α (x i,j,k + 0.5∆y j e y ) − f n+1/2 α (x i,j,k − 0.5∆y j e y )]ξ α,y ∆x i ∆z k +[ f n+1/2 α (x i,j,k + 0.5∆z k e z ) − f n+1/2 α (x i,j,k − 0.5∆z k e z )]ξ α,z ∆x i ∆y j (5) where f n+1/2 α (x b ) represents the distribution at the cell interface x b at the time t n+1/2 = t n + 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: The collision term can be expanded in the BGK collision model, and Equation (6) can be converted to the following equations: The evolution equation of DUGKS from t n to t n+1 is simplified as: Based on the conservation of mass, momentum, the density ρ, and velocity u can be computed from f α : , 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: Similarly, new distribution functions are introduced and can be computed by expanding the collision term: With new distribution functions, Equation (10) is converted to the following equation: With the Taylor expansion around the endpoint (x b ) located at the cell interface, the right term of Equation (13) can be approximated as: where f +,n α (x b ) and its gradients ∇ f +,n α (x b ) can be approximated by the linear interpolations. In x-direction, e.g., The distribution function f +,n α can be computed from f α , as follows, Then, we obtain the function f n+1/2 α (x b ). The density and velocity at the cell interface can also be evaluated, which can be used for the equilibrium distribution function f eq,n+1/2 α Finally, the flux F n+1/2 α,j is evaluated according to Equation (5) after the distribution function f n+1/2 α at the cell interface is determined by Equation (11). The tracked distribution function f α can be updated to the next time step after the flux is obtained. Particularly, the following equation will be used in the DUGKS, 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.

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 f α , unlike the LBM.
Considering the impermeable wall boundary (U Wz = 0), the unknown distributions are obtained by the specular reflection ( f sr α ) and the stress exerted by the wall ( f w α ).
As shown in Figure 2, the unknown distributions are f 4 , f 8 , f 9 , f 12 , f 14 .  The specular reflection sr f  can be obtained by: 11  14  13  12  7  9  10  8  3   4~ , , , , can be obtained by: The stress w f  contributes to the change in tangential momentum caused by the wall surface, as shown in the following equations: The specular reflection f sr α can be obtained by: With f sr α determined, U SR can be obtained by: The stress f w α contributes to the change in tangential momentum caused by the wall surface, as shown in the following equations: 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. σ represents a modified tangential momentum accommodation coefficient. n is the normal direction of the wall toward the fluid field. The subscripts "||" and "⊥" represent the tangential and normal directions, respectively. The sum of normal parts is zero, ensuring that the function f w α only changes the tangential momentum. It is also shown that the calculation of density is not determined by ∑ f w α . The following relations can be obtained for the case in Figure 2.
The momentum can change due to the shear stress imposed on the wall: For positive normal direction of the wall (i.e., +z direction): Based on the Maxwell equilibrium distribution function and Ref. [22], the density can be calculated, Then, f w 8 , f w 9 , f w 12 , f w 14 and U SRx , U SRy are still unknown. σ x , σ 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: To make the system closed, the hypothesis is proposed: In summary, a new slip boundary condition is proposed for the upper horizontal wall boundary in 3D: With the external forcing term, the local velocities u are computed by, To eliminate the numerical slip due to force tangential to the wall, the external forcing term is introduced to the new slip boundary condition: Similar manipulations can be applied to the lower wall and side walls boundary.

Relation between the Combination Parameters and Slip Lengths
Then, the relation between combination parameters (σ x , σ 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]: where φ 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 σ x , σ y and the slip lengths b x , b y .
The flow near the wall satisfies the continuity equation: ∂ρ ∂t In the 3D directional flows, the velocity distribution near the wall is assumed to be satisfied the following equations: As shown in Figure 2, the local coordinate system (x, y, z) in 3D is established in the lattice unit, where a node on the wall is served as the origin. The x − y plane is located on the wall. The z direction is normal along the wall. In the local coordinate system, the boundary of the upper wall is located at the plane z = 1, and the solid is located in the region (z < 1), e.g., the plane 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 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: The acceleration can be approximated by relations as follows [69]: Then, the coefficients can be obtained by the acceleration, α 1 = −0.5 a x /ν, α 2 = −0.5 a y /ν. The slip velocity can be expressed as: 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]: Considering Equations (42) and (43), (46) can be written as: With the Taylor expansion around the z = 1 in the local coordinate system, u x (z) and u y (z) can be approximated, With the assumption of parabolic velocity profiles near the wall, substitute the above equations, and the following equations can be obtained, Then, the relations between the coefficients are given as, For D3Q19 lattice model, with ρu = 18 ∑ α=0 ξ α f α known, the local velocities u x , u y at z = 1 and 2 can be calculated as, (54) Combined with the proposed slip boundary condition, the following relations can be obtained: Then, The difference value between velocity at z = 1 and z = 2 can be written as: 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: where f denotes the tracked distribution function in the collision. Then, Equation (57) can be simplified as follows: The slip velocity could be calculated as: Finally, the relations between the slip lengths and parameters are obtained: With the known values of α 1 /β 1 , α 2 /β 2 , Equations (64) and (65) can be simplified, where A denotes the correction coefficient and is determined by: where ∆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.

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].
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 ρu x , ρu y , ρu z , the momentum fluxes and shear stresses Π xx , Π yy , Π zz , Π xy , Π xz , Π yz , and three higher-order moments Q xxy , Q yyz , Q xzz as 12 linearly independent moments.
The moments can be approximated by the Chapman-Enskog expansion as follows: where the equilibrium part of the momentum flux tensor (Π (0) αβ ) and the third-order moment (Q (0) αβ ) can be expressed as: The velocity is set to 0 at the corner, and some terms are assumed to be negligible. The momentum fluxes and shear stresses Π xx , Π yy , Π zz , Π xy , Π xz , Π yz , and three higher-order moments Q xxy , Q yyz , Q xzz are written as follows: The unknown functions are calculated as: The density at the low-south-west corner is calculated as: Similar manipulations can be applied to other corner nodes.

Algorithm
The updating of f α from time t n = n∆t to time t n+1 = (n + 1)∆t in the DUGKS can be performed as the following brief algorithm.
Obtain the values of f eq,n α and f n α at time t = 0.

3.
Compute the value of ∇ f Compute the distribution function f (14) and (13).
For the unknown distribution functions at the boundary or corner, the boundary conditions are employed, such as Equations (35) or (73). 8.
Update the distribution function f n+1 α using Equation (8). Obtain the macro values of density and velocity. 9.

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.

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. periodic boundary condition. The remaining four walls adopt the no-slip/slip bound conditions. For the no-slip case, the bounce-back scheme in LBM is used without trea the corner in Ref. [11], and the bounce-back scheme in DUGKS is used with the co boundary condition in the present work. For the slip case, the BSR scheme is employe 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] experimental result [9], as shown in Figure 4a. It is observed that the present results a a little better with the exact solution than the BSR scheme in Ref. [11], which shows 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.  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.

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.

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: . δ 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 Figures 5 and 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.
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 Figures 6. As shown in Figures 5 and 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]. 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|0)/(dp/dx|0), where dp/dx and dp/dx|0 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. The percentage change in drag (DD) is defined by DD = (dp/dx-dp/dx| 0 )/(dp/dx| 0 ), where dp/dx and dp/dx| 0 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.

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 = U0cos(ωt), where U0 = 1.0 m/s, oscillating frequency ω = 1875 π/s. The

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 0 cos(ωt), where U 0 = 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 ρ = 1000 kg/m 3 and υ = 1.0 × 10 −6 m 2 /s, respectively. The Reynolds number is calculated as Re = U 0 × L/ν = 1.0 × 0.0004/10 −6 = 400. For simplicity, ω has been dimensionalized as ω' = ωL/U 0 = 0.75π, and St = ωL 2 /ν = ω'Re = 300π.
The problem is the micro-lid-driven cavity flow

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 3 (coarse), 120 3 (medium), and 160 3 (fine). Figure 9 shows the negligible improvement in the velocity profile on increasing the lattice size from 120 3 to 160 3 .
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.

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 3 (coarse), 120 3 (medium), and 160 3 (fine). Figure 9 shows the negligible improvement in the velocity profile on increasing the lattice size from 120 3 to 160 3 . 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, U0lattice = 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.

Results and Discussion
There  [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) 2 + (∂U/∂Z-∂W/∂X) 2 + (∂V/∂X-∂U/∂Y) 2 }. 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 Figures 10-22.
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) 2 + (∂U/∂Z-∂W/∂X) 2 + (∂V/∂X-∂U/∂Y) 2 }.
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 Figures 10-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 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 0 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 0 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, 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.  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 bx (bx = 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 bx has greater influence on U than by. It is found that velocity profiles in each group are very close. Therefore, for bx at the same level, the existence of by on the top or bottom wall has almost no influence on the change in U along the centerline Z-axis. The greater bx 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

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, microlid-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.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/e24070907/s1. Figure S1: Contours for Velocity U, V, W and Vorticity magnitude at t = T in case (a) and case (b). Figure S2: Contours for Velocity U at t = T in cases (c-n). Figure S3: Contours for Velocity V at t = T in cases (c-n). Figure S4: Contours for Velocity W at t = T in cases (c-n). Figure S5: Contours for Vorticity magnitude at t = T in cases (c-n). Figure S6: Contours for Velocity U, V, W and Vorticity magnitude at t = 0.25 T in case (a) and case (b). Figure S7: Contours for Velocity U at t = 0.25 T in cases (c-n). Figure S8: Contours for Velocity V at t = 0.25 T in cases (c-n). Figure S9: Contours for Velocity W at t = 0.25 T in cases (c-n). Figure S10: Contours for Vorticity magnitude at t = 0.25 T in cases (c-n). Figure S11: Con-tours for Velocity U, V, W and Vorticity magnitude at t = 0.5 T in case (a) and case (b). Figure S12: Contours for Veloc-ity U at t = 0.5 T in cases (c-n). Figure S13: Contours for Velocity V at t = 0.5 T in cases (c-n). Figure S14: Contours for Velocity W at t = 0.5 T in cases (c-n). Figure S15: Contours for Vorticity magnitude at t = 0.5 T in cases (c-n). Figure S16: Contours for Velocity U, V, W and Vorticity magnitude at t = 0.75 T in case (a) and case (b). Figure S17: Contours for Velocity U at t = 0.75 T in cases (c-n). Figure S18: Contours for Velocity V at t = 0.75 T in cases (c-n). Figure S19: Contours for Velocity W at t = 0.75 T in cases (c-n). Figure S20: Contours for Vorticity magnitude at t = 0.75 T in cases (c-n).

Data Availability Statement:
The data that support the findings of this study are available within the article.

Conflicts of Interest:
The authors declare no conflict of interest.