Lattice Boltzmann Simulation of Fluid Flow Characteristics in a Rock Micro-Fracture Based on the Pseudo-Potential Model

Abstract

Slip boundary has an important influence on fluid flow, which is non-negligible in rock micro-fractures. In this paper, an improved pseudo-potential multi-relaxation-time (MRT) lattice Boltzmann method (LBM), which can achieve a large density ratio, is introduced to simulate the fluid flow in a micro-fracture. The model is tested to satisfy thermodynamic consistency and simulate Poiseuille flow in the case of large liquid-gas density ratio. The slip length is used as an index for evaluating the flow characteristics, and the effects of wall wettability, micro-fracture width, driving pressure and liquid-gas density ratio on the slip length are discussed. The results demonstrate that the slip length increases significantly with the increase of the wall contact angle in rock micro-fracture. And the liquid-gas density ratio has an important impact on the slip length, especially for the hydrophobic wall. Moreover, under the laminar flow regime the driving pressure and the micro-fracture width has little effect on the slip length.

1. Introduction

The fluid flow in rock micro-fractures is a topic of great importance for a wide range of scientific problems, such as water conservancy, oil recovery in low-permeability oilfields, nuclear waste treatment and others [1,2,3,4]. Classical Navier-Stokes hydrodynamics with non-slip boundaries, in which the fluid velocity at a fluid-solid interface equals the solid velocity, is successfully applied to seepage research at the macro-scale. However, a series of experimental and numerical results indicate that the non-slip boundary notion is invalid at the micro-scale [5,6]. Besides, there are a large number of micro-fractures in the rock with the width < 1 mm [7]. The slip boundaries should be taken into account during the fluid flow in these micro-fractures [8,9]. Meanwhile, the surfaces of micro-fractures show different wettabilities due to its diverse mineral composition [10], and it has a strong influence on the wall slip. Therefore, great attention should be paid to studying the slip conditions in micro-fractures with different wettabilities.

Slip length is an important dynamic parameter for quantifying the slip condition of a liquid flowing over the solid surface, and it’s difficult to measure directly through experiments due to the objective limits of the current techniques [11,12]. Thus, slip length studies in micro-fluidics have focused on numerical simulations. The lattice Boltzmann method (LBM) offers great potential for micro-scale flow simulations owing to its mesoscopic particle background [13]. Many scholars have studied micro-scale fluid flow with the LBM and proposed several typical multiphase models, such as the color model, the free-energy model, the pseudo-potential model, etc. Among these, the pseudo-potential model realizes the phase separation by implementing micro-molecular interactions. It has received considerable attention due to the fact the interface between different phases does not need to be tracked or captured. Based on this model, Zhang et al. [14] found that the non-slip boundary conditions could produce an apparent slip. Chen et al. [15] simulated the Couette flow and found that there was a direct relationship between the magnitude of slip length and the strength of solid-fluid interaction. Zhang et al. [16] investigated the droplet flow velocity in the micro-channels, and discussed the wall slip caused by the roughness. Kunert [17] studied the influence of fluid viscosity on the slip length. However, almost all the studies do not consider the effect of liquid-gas density ratio on the wall slip due to the difficulty of achieving a large liquid-gas density ratio with the original pseudo-potential model [18].

Based on the LBM, an improved pseudo-potential model is introduced to simulate the liquid flow in a rock micro-fracture, which could achieve a large liquid-gas density ratio. The validity of the proposed model is tested on two benchmark cases: the thermodynamic consistency and Poiseuille flow. The slip length is discussed considering the effects of wall wettability, micro-fracture width, driving pressure and liquid-gas density ratio.

2. Slip Condition

Bernoulli first proposed the non-slip boundary model (Figure 1a) where the fluid velocity was zero relative to the physical wall. Subsequently, some scholars have supposed that the fluid layer near the wall has a certain thickness (Figure 1b), in which no relative movement occurs between the fluid and the wall. In addition, some experimental results indicated that the fluid velocity was not zero at the wall (Figure 1c). The above three slip conditions have been confirmed, and they are deemed to be related to the wall wettability. When the fluid flow over a strong hydrophilic wall, it presents non-slip or negative slip, and a fluid usually shows a positive slip condition at a strong hydrophobic wall. As shown in Figure 1c, based on the assumption that the distance to the virtual wall, at which the extrapolated fluid velocity at a constant shear would be equal to the wall velocity, Navier defined the slip length $l$ as:

$$l=v_{\mathrm{s}}/{\tau }_{\mathrm{wall}},$$

(1)

where $v_{\mathrm{s}}$ is the slip velocity at the physical wall; ${\tau }_{\mathrm{wall}}$ is the local shear near the physical wall.

3. Coexistence Densities and Maxwell Construction Rule

In nature, the fluids in rock micro-fractures are mostly in a liquid-gas two phase situation. A two phase fluid could be described by its equation of state, which is an equation relating the fluid state variables in thermodynamics, such as pressure $P$, molar volume $V$, and temperature $T$. Figure 2 shows a series of $P$ − $V$ curves of an equation of state with different temperatures, where the bulk (the region not close to any interface) pressure $P\left(V,T\right)$ is a function of the molar volume $V$ and the temperature $T$. For different temperatures, the fluid may show subcritical, critical, and supercritical behaviors. At high temperature (above the critical temperature $T_c$) the fluid shows supercritical behavior and no distinct liquid and gas phases can be discerned. Below the critical temperature, phase separation into liquid and gas is possible, and different molar volumes $V_l$ and $V_g$ of the certain substance may coexist at a single pressure $P_0$ at equilibrium, where $V_l$ and $V_g$ are the liquid and gas phases molar volumes, respectively. At the critical temperature, the coexistence of different molar volumes is also impossible, nevertheless the first and second derivative of $P$ at $T_c$ should be zero [19]:

$$\{\begin{array}{l}\frac{\partial P\left(T_{\mathrm{c}}\right)}{\partial V}=0\\ \frac{{\partial }^{2}P\left(T_{\mathrm{c}}\right)}{\partial V^2}=0\end{array},$$

(2)

The Carnahan-Starling equation of state (C-S EOS) is one of classic equations of state for non-ideal fluid, which is given by [19]:

$$P\left(\rho ,T\right)=\rho RT\frac{1+b\rho /4+{\left(b\rho /4\right)}^2-{\left(b\rho /4\right)}^3}{{\left(1-b\rho /4\right)}^3}-a{\rho }^2,$$

(3)

where $a$ is the attraction parameter; $b$ is the repulsion parameter;$R$ is the gas constant and $\rho$ is the fluid density. Rearranging and taking the derivatives of Equation (3), $a$ and $b$ can be presented as function of the critical pressure $P_c$ and the $T_c$:

$$\{\begin{array}{l}a=0.4963R^2{T}_c^2/P_c\\ b=0.18727R{T}_c/P_c\end{array},$$

(4)

Then, $T_c$ can be given as:

$$T_c=0.3773a/\left(bR\right),$$

(5)

In this paper, the parameters in Equation (3) are set as $a=1$, $b=4$, $R=1$, respectively [20]. Then $T_c=0.0943$. The reduced temperature coefficient $t_r$ ($0<t_r<1$) is introduced to ensure the temperature $T=t_r{T}_c$ below the critical temperature, which makes the system at the subcritical behavior. The coexistence of different molar volumes at a single pressure could be found according to Maxwell construction rule [13]:

$$\{\begin{array}{l}{\displaystyle {\int }_{V_g}^{V_l}\left[P_0-P\left(V,T\right)\right]}\mathrm{d}V=0\\ P_0=P\left(V_{\mathrm{g}},T\right)=P\left(V_l,T\right)\end{array},$$

(6)

Replacing the molar volume $V$ as the reciprocal of density $V=1/\rho$, Maxwell construction rule Equation (6) can be rewritten as:

$$\{\begin{array}{l}{\displaystyle {\int }_{{\rho }_g}^{{\rho }_l}\left[P_0-P\left(\rho ,T\right)\right]}\frac{1}{{\rho }^2}\mathrm{d}\rho =0\\ P_0=P\left({\rho }_{\mathrm{g}},T\right)=P\left({\rho }_l,T\right)\end{array},$$

(7)

where ${\rho }_l$ and ${\rho }_{\mathrm{g}}$ are the liquid and gas phase coexistence densities, respectively. A series of coexistence densities with different temperatures could be directly obtained by numerical method according to Maxwell construction rule.

4. Numerical Model

4.1. The MRT-LBM

Using the Guo’s force scheme [21], the evolution equation of LBM can be given as follows:

$$\mathrm{Collision}\mathrm{step}:f_{\alpha }(r,t+{\delta }_t)-f_{\alpha }(r,t)={\mathsf{\Omega }}_{\alpha }+{\delta }_t{F}_{\alpha },$$

(8)

$$\mathrm{Streaming}\mathrm{step}:f_{\alpha }(r+{\mathit{e}}_{\alpha }{\delta }_t,t+{\delta }_t)=f_{\alpha }(r,t+{\delta }_t),$$

(9)

where $f_{\alpha }$ is the density distribution function in the $\alpha$ direction; $r$ represents the coordinate of the node;${\mathit{e}}_{\alpha }$ is the discrete velocity;${\delta }_t$ is the time step; ${\mathsf{\Omega }}_{\alpha }$ is the collision operator; $F_{\alpha }$ is the forcing term.

For the D2Q9 model:

$${\mathit{e}}_{\alpha }=\{\begin{array}{ll}\left(0,0\right)& \alpha =0\\ \left(\cos \left[\left(\alpha -1\right)\pi /2\right],sin\left[\left(\alpha -1\right)\pi /2\right]\right)& \alpha =1,2,3,4\\ \left(\sqrt{2}\cos \left[\left(\alpha -5\right)\pi /2+\pi /4\right],\sqrt{2}sin\left[\left(\alpha -5\right)\pi /2+\pi /4\right]\right)& \alpha =5,6,7,8\end{array},$$

(10)

The single-relaxation-time lattice Boltzmann method (SRT-LBM) can simplify the collision operator ${\mathsf{\Omega }}_{\alpha }$ of the Equation (8). However, it comes at a cost of reduced accuracy and stability. Compared with SRT-LBM, the multi-relaxation-time lattice Boltzmann method (MRT-LBM), which offers a series of relaxation-time can overcome these problems [22]. So the MRT-LBM is applied in the present work. The collision operator ${\mathsf{\Omega }}_{\alpha }$ and the Guo’s forcing term $F_{\alpha }$ are respectively given by [23]:

$${\mathsf{\Omega }}_{\alpha }=-{\left(M^{-1}\Lambda M\right)}_{\alpha \beta }\left(f_{\beta }(r,t)-f_{\beta }^{\mathrm{eq}}(r,t)\right),$$

(11)

$$F_{\alpha }=S_{\alpha }-0.5{\left(M^{-1}\Lambda M\right)}_{\alpha \beta }{S}_{\beta },$$

(12)

where $M$ is an orthogonal transformation matrix; $\Lambda$ is a diagonal matrix; $f^{\mathrm{eq}}$ is the equilibrium density distribution function;$S$ is the forcing term in the moment space; $\alpha$ and $\beta$ is the row and column number of a matrix, respectively. For the D2Q9 model, $M$ and $\Lambda$ can be written as:

$$M=\left[\begin{array}{rrrrrrrrr}1& 1& 1& 1& 1& 1& 1& 1& 1\\ -4& -1& -1& -1& -1& 2& 2& 2& 2\\ 4& -2& -2& -2& -2& 1& 1& 1& 1\\ 0& 1& 0& -1& 0& 1& -1& -1& 1\\ 0& -2& 0& 2& 0& 1& -1& -1& 1\\ 0& 0& 1& 0& -1& 1& 1& -1& -1\\ 0& 0& -2& 0& 2& 1& 1& -1& -1\\ 0& 1& -1& 1& -1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& -1& 1& -1\end{array}\right],$$

(13)

$$\Lambda =diag\left({\tau }_{\rho }^{-1},{\tau }_e^{-1},{\tau }_{\zeta }^{-1},{\tau }_j^{-1},{\tau }_q^{-1},{\tau }_j^{-1},{\tau }_q^{-1},{\tau }_{\upsilon }^{-1},{\tau }_{\upsilon }^{-1}\right),$$

(14)

where $M$ contains all the relaxation rates; it is set as $\Lambda =diag\left(1.0,1.1,1.1,1.0,1.1,1.0,1.1,0.8,0.8\right)$ [20] in this paper. The kinematic viscosity $\upsilon$ can be derived from one of the relaxation coefficients ${\tau }_{\upsilon }$: $\upsilon =c_s^2\left({\tau }_{\upsilon }-0.5\right){\delta }_t$, where $c_s=c/\sqrt{3}$ is lattice sound speed, and $c$ is the lattice constant, which is set as 1 for the present.

The moments $m$ can be obtained by $m=Mf$; the equilibrium moments $m^{\mathrm{eq}}$ can be straightforwardly computed by $m^{\mathrm{eq}}=M{f}^{\mathrm{eq}}$. Alternatively, it can be constructed more precisely and efficiently from the density $\rho$ and the velocity $v$ by [20]:

$$m^{eq}=\rho {\left[1,\left(-2+3{\left|v\right|}^2\right),\left(1-3{\left|v\right|}^2\right),v_x,-v_x,v_y,-v_y,\left(v_x^2-v_y^2\right),v_{\mathrm{x}}{v}_y\right]}^{\mathrm{T}},$$

(15)

where $v=\left(v_x,v_y\right)$ is the macroscopic velocity given by $v={\displaystyle \sum _{\alpha =1}^N{e}_{\alpha }{f}_{\alpha }/\rho }+F/2\rho$; and $\rho ={\displaystyle \sum _{\alpha =0}^N{f}_{\alpha }}$ is the macroscopic density; $F=\left(F_x,F_y\right)$ is the total force.

Through Equations (11) and (12), the Equation (8) can be transformed as follows:

$$m^{*}=m-\Lambda \left(m-m^{eq}\right)+{\delta }_t\left(I-\Lambda /2\right)\overline{S},$$

(16)

where $m^{*}$ are the moments after the collision step; $I$ is the unit tensor; $\overline{S}$ is the forcing term in the moment space which is given by [20]:

$$\overline{S}={\left[0,6\left(v_{\mathrm{x}}{F}_x+v_y{F}_y\right),-6\left(v_{\mathrm{x}}{F}_x+v_y{F}_y\right),F_x,-F_x,F_y,-F_y,2\left(v_{\mathrm{x}}{F}_x-v_y{F}_y\right),\left(v_{\mathrm{x}}{F}_y+v_y{F}_x\right)\right]}^{\mathrm{T}},$$

(17)

After the collision step, the moments $m^{*}$ should be transformed to the density distribution function by $f^{*}=M^{-1}{m}^{*}$, then the $f^{*}$ can be streamed to the next node through Equation (9).

4.2. The Original Pseudo-Potential Model

The total force in Equation (17) is given by $F=F_{\mathrm{ex}}+F_{\mathrm{sc}}$, where $F_{\mathrm{ex}}$ is the external force such as gravitational force and buoyancy force. The inter-molecular force $F_{\mathrm{sc}}=F_{\mathrm{int}}+F_{\mathrm{w}}$ can be obtained by the pseudo-potential $\psi \left(\rho \right)$ presented in the following, which is a function to mimic the molecular interactions that cause phase separation. $F_{\mathrm{int}}$ is the inter-molecular force within fluids; $F_{\mathrm{w}}$ is the adhesive force between fluid particles and solid wall. Both $F_{\mathrm{int}}$ and $F_{\mathrm{w}}$ can be expressed as [24,25,26]:

$$F_{\mathrm{int}}=-G\psi \left(\rho \left(r\right)\right){\displaystyle \sum _{\alpha =1}^{N}w\left({\left|e_{\alpha }\right|}^2\right)}\psi \left(\rho \left(r+e_{\alpha }\right)\right)e_{\alpha },$$

(18)

$$F_w=-GS\left(r\right){\displaystyle \sum _{\alpha =1}^{N}w\left({\left|e_{\alpha }\right|}^2\right)}\psi \left(\rho \left(r+e_{\alpha }\right)\right)e_{\alpha },$$

(19)

$$S\left(r\right)=\{\begin{array}{ll}0& \left(r\ne r_w\right)\\ \psi \left({\rho }_w\right)& \left(r=r_w\right)\end{array},$$

(20)

where $G$ is the interaction strength, which is set as 6.0 in the present paper; $w\left({\left|e_{\alpha }\right|}^2\right)$ is the weight, $w\left(1\right)=1/3$ and $w\left(2\right)=1/12$; $r_w$ represents the solid wall node; ${\rho }_w$ is the wall density. It should be noticed that ${\rho }_w$ is not a real physical density. It is merely an artificial parameter to acquire different wall affinities.

Using the Chapman-Enskog analysis, the following macroscopic equation can be derived from Equations (11), (15)–(17) [20]:

$$\nabla \cdot P=\nabla \cdot \left(\rho c_s^{2}I\right)-F,$$

(21)

Through the Taylor expansion, the Equation (18) can be rewritten as [26,27]:

$$\begin{array}{ll}F& =-G{c}^2\left[\psi \nabla \psi +\frac{c^2}{6}\psi \nabla \left({\nabla }^2\psi \right)\right]+\cdots \\ & =-\frac{G{c}^2}{2}\nabla {\psi }^2-\frac{G{c}^4}{6}\left[\nabla \left(\psi {\nabla }^2\psi \right)-{\nabla }^2\psi \nabla \psi \right]+\cdots \\ & =-\frac{G{c}^2}{2}\nabla {\psi }^2-\frac{G{c}^4}{6}\nabla \left(\psi {\nabla }^2\psi \right)+\frac{G{c}^4}{6}\left[\nabla \cdot \left(\nabla \psi \nabla \psi \right)-\frac{1}{2}\nabla |\nabla {\psi |}^2\right]+\cdots \end{array}$$

(22)

Combining Equations (21) and (22), the normal pressure tensor $P_{\mathrm{n}}$ of the interface can be given by [26]:

$$P_{\mathrm{n}}=\rho c_s^2+\frac{G{c}^2}{2}{\psi }^2+\frac{G{c}^4}{12}\left[\alpha {\left(\frac{\mathrm{d}\psi }{\mathrm{d}n}\right)}^2+\beta \psi \frac{{\mathrm{d}}^2\psi }{\mathrm{d}{n}^2}\right],$$

(23)

where the subscript n denotes the normal direction of the interface; $\alpha$ and $\beta$ are constants for the D2Q9 model. According to Equation (23) and the physical requirement that the interface pressure should be equal to the bulk (the region not close to any interface) pressure at equilibrium, the following mechanical stability condition can be written as:

$${\displaystyle {\int }_{{\rho }_g}^{{\rho }_l}\left(P_0-\rho c_s^2-\frac{G{c}^2}{2}{\psi }^2\right)}\frac{{\psi }^{\prime }}{{\psi }^{1+\epsilon }}\mathrm{d}\rho =0,$$

(24)

where ${\psi }^{\prime }=d\psi /d\rho$, $\epsilon =-2\alpha /\beta$, and $P_0=P\left({\rho }_l\right)=P\left({\rho }_g\right)$. In the pseudo-potential LB model, the coexistence densities (${\rho }_l$ and ${\rho }_g$) are determined by the mechanical stability condition Equation (24). Maxwell construction rule, which determines the thermodynamic coexistence, is built in the requirement of Equation (7). To satisfy the Equation (7), Sbragaglia and Shan [28] have proposed the following interaction potential $\psi \left(\rho \right)$:

$$\psi \left(\rho \right)=\{\begin{array}{ll}\exp \left(-1/\rho \right)& ,\epsilon =0\\ {\left(\frac{\rho }{\epsilon +\rho }\right)}^{1/\epsilon }& ,\epsilon \ne 0\end{array},$$

(25)

which gives ${\psi }^{\prime }/{\psi }^{1+\epsilon }=1/{\rho }^2$. With the above potential function, Equations (7) and (22) will be same, thereby the thermodynamic consistency (Maxwell construction rule) will be obtained. But, to be consistent with the equation of state in the thermodynamic theory, the potential $\psi \left(\rho \right)$ must be chosen as [18,29]:

$$\psi \left(\rho \right)=\sqrt{\frac{2\left(P-\rho {\mathrm{c}}_s^2\right)}{G{c}^2}},$$

(26)

Obviously, Equations (25) and (26) cannot be obtained at the same time, so this original pseudo-potential model cannot content the thermodynamic consistency.

4.3. The Improved Pseudo-Potential Model

To resolve the problem of thermodynamic inconsistency, the parameter $\epsilon$ in Equation (24) can be tuned to approximately achieve the thermodynamic consistency when $\psi \left(\rho \right)$ is defined by Equation (26). Through changing the original forcing term, Equation (17) is rewritten as follows [20]:

$${\overline{S}}^{*}=\left[\begin{array}{c}0\\ 6\left(v_{\mathrm{x}}{F}_x+v_y{F}_y\right)+\frac{12\sigma {\left|F\right|}^2}{{\psi }^2{\delta }_{\mathrm{t}}\left({\tau }_e-0.5\right)}\\ -6\left(v_{\mathrm{x}}{F}_x+v_y{F}_y\right)-\frac{12\sigma {\left|F\right|}^2}{{\psi }^2{\delta }_{\mathrm{t}}\left({\tau }_{\varsigma }-0.5\right)}\\ F_x\\ -F_x\\ F_y\\ -F_y\\ 2\left(v_{\mathrm{x}}{F}_x-v_y{F}_y\right)\\ \left(v_{\mathrm{x}}{F}_y-v_y{F}_x\right)\end{array}\right],$$

(27)

Meanwhile, Equation (23) could be given as:

$$P_{\mathrm{n}}=\rho c_s^2+\frac{G{c}^2}{2}{\psi }^2+\frac{G{c}^4}{12}\left[\left(\alpha +24G\sigma \right){\left(\frac{d\psi }{dn}\right)}^2+\beta \psi \frac{d^2\psi }{d{n}^2}\right],$$

(28)

which leads to $\epsilon =-2\left(\alpha +24G\sigma \right)/\beta$, where $\sigma$ is a correction parameter which could tune the mechanical stability condition Equation (24) to approximately get the thermodynamic consistency.

4.4. Flowchart of the MRT-LBM

Based on the improved pseudo-potential model, the MRT-LBM has been proposed to simulate the fluid flow in the rock micro-fracture. The flowchart is shown in Figure 3.

4.5. Verification

Two reference problems are simulated to verify the improved pseudo-potential MRT-LBM in this section. One case is a bubble in gas, which is applied to test the thermodynamic consistency with large liquid-gas density ratio. The other is Poiseuille flow, in which the slip boundary with the slip length $l$ close to zero is applied to simulate the non-slip behavior. It should be noticed that the units in this paper are all lattice units.

4.5.1. The Thermodynamic Consistency

A bubble in gas is simulated in a 100 × 100 lattice computational domain. The periodic boundary is set for its four sides, and the bubble is placed in the middle. The liquid-gas density coexistence curves are obtained by changing the temperature $T=t_r{T}_c$, which are compared with the theoretical curve given by Maxwell construction rule [30]. As shown in Figure 4, the results of the original forcing scheme ($\sigma =0$), and the improved model with $\sigma =0.09,0.12$ are in good agreement with the theoretical results in the liquid phase, but there are large deviation between them in the gas phase. The numerical solutions of the improved forcing scheme with $\sigma =0.11$ agrees well with the theoretical ones in both phases, which indicate that the proposed model could satisfy the thermodynamic consistency in the case of large liquid-gas density ratio, so the correction parameter $\sigma$ is set as 0.11 in the following work.

4.5.2. Poiseuille Flow

Poiseuille flow is simulated with the slip length $l$ close to zero ($l=-0.02$). A 800 × 100 uniform lattice system is chosen in this case, the top and bottom boundaries are solid walls with the bounce-back boundaries, and the left and right sides are set as the periodic boundaries. To analyze the influence of different discretization features, we simulate the Poiseuille flow with the different mesh sizes (400 × 50, 800 × 100, and 1600 × 200). In this test, the liquid kinematic viscosity $u$ is 0.1, the liquid is driven by the pressure of 3 × 10⁻⁴ along the horizontal direction. The wall density ${\rho }_w$ is set as 0.12 to ensure the slip length $l$ close to zero ($l=-0.02$). The reduced temperature coefficient $t_r$ is 0.55. The correction parameter $\sigma$ is 0.11. Other parameters are set as follows: lattice spaces are equal in horizontal and vertical direction, ${\delta }_x={\delta }_y=1$, and the time step is ${\delta }_t=1$. Figure 5 shows the velocity profile obtained by the proposed model, which are compared with the theoretical solutions of Poiseuille flow (non-slip behavior) [13]. It is obvious that the fluid velocity given by numerical method shows a good match with the theoretical results and the grids number have little influence on the computed result. The maximum error is less than 3.24%. Hence, the improved pseudo-potential MRT-LBM can be used to simulate the fluid flow in the narrow channel successfully.

5. Results and Discussion

In the present study, the slip condition is discussed in a rock micro-fracture with different wettabilities. The wall wettability is the macroscopic expression of the interaction between fluid particles and the solid wall, so the relationship between the contact angle and the wall density is tested based on the improved the improved pseudo-potential model and the slip length, as an important parameter for judging the slip condition, is investigated considering the effects of contact angle, driving pressure, and liquid-gas density ratio.

5.1. Contact Angle Test

To obtain the relationship between the contact angle and the wall density, a 200 × 200 lattice domain is used, the top and bottom boundaries are solid walls and the left and right sides are set as the periodic boundaries. A liquid droplet is placed at the solid surface. The liquid-gas density ratio is set as 293:1 in this test. Different wall densities are applied to describe the changing wall affinities, Figure 6 presents different contact angles, and the relationship between the contact angle and the wall density is listed in

Table 1. Relationship between the contact angle and the wall density.

Wall Density (ρw)	0.02	0.06	0.10	0.14	0.18
Contact Angle (θ)	148.7°	115.9°	92.6°	76.0°	61.9°

.

5.2. Discussion

Based on the improved pseudo-potential MRT-LBM, a size of 800 × 100 lattice domain is applied to simulate the rock micro-fracture with different wettabilities. As shown in Figure 7, the left and right sides of the micro-fracture are set as the periodic boundaries; the top and bottom of the model are the solid walls. The volume ratio of liquid to gas is 9:2 in the initial state, and the liquid phase is placed in the middle of the micro-fracture. To ensure a laminar flow regime, the driving pressure of 3 × 10⁻⁴ is applied in the horizontal direction after the liquid-gas phases separated, in which the calculated Reynolds number is 1.49. The slip length, as an important parameter for judging the slip degree, is used to discuss the effects of contact angle, driving pressure, and liquid-gas density ratio on the fluid flow state in the micro-fracture.

5.2.1. Effect of the Contact Angle

In order to study the effect of the contact angle on the slip length, the liquid-gas density ratio is set as 293:1, and the micro-fracture widths are set as 25, 50, and 100, respectively. Figure 8 shows the slip length against the contact angle.

It can be observed that the slip length increases from −0.27 to 4.9 as the contact angle increasing from 61.9° to 148.7°. The larger the contact angle, the more significant the slip length changes. This finding is in agreement with that by Tsu-Hsu Yen [31] and Bladimir Ramos-Alvarado [32]. Moreover, in the case of different fracture widths, the variation trends of slip length with the contact angle are very close to each other, which indicate that the micro-fracture width has little effect on the slip length.

To analyze the physical mechanism of wall slip, the liquid density profiles (the fracture width is 100) for different contact angles are shown in Figure 9. As the contact angle increases, the liquid density near the wall decreases and the slip length increases, which confirms the existing slip theory [33,34]. The wall slip is induced by the changes of the liquid density near the solid surface, which is attributed to the varying interaction force between the liquid and the solid wall. With the increases of interaction force, the liquid density near the wall increases. Thereby it shows non-slip or even negative slip behaviors. As the interaction force decreases, the liquid density close to the wall surface decreases and the slip length increases. Even when the liquid density is less than a certain value, the liquid density near the wall surface is reduced to the gas density, so that a “gas layer” is formed at the solid surface.

5.2.2. Effect of the Driving Pressure

To investigate the effect of the driving pressure on the slip length, the micro-fracture width W is set as 100, the liquid-gas density ratio is 293:1, and the contact angles are 61.9°, 76.0°, 92.6°, 115.9° and 148.7°, respectively. When the driving pressure ranges from 3 × 10⁻⁴ to 1.5 × 10⁻³, the maximum Reynolds number is 7.48, which shows that the liquid flow is in a laminar flow regime. Figure 10 presents the relationship between the driving pressure and the average flow velocity, which indicates that the flow velocity linearly increases with the driving pressure. In Figure 11, the relationship between the driving pressure and the slip length shows that the driving pressure has almost no influence on the slip length. Xiang [35] and Huang [36] have also yielded the similar results by other numerical methods.

5.2.3. Effect of the Liquid-Gas Density Ratio

To study the influence of liquid-gas density ratio on the slip length, the liquid-gas density ratios are chosen in the range from 22:1 to 293:1. In this case, the micro-fracture width is 100, the driving pressure is 3 × 10⁻⁴, and the contact angles are 61.9°, 76.0°, 92.6°, 115.9°, and 148.7°, respectively. It should be noticed that changing liquid-gas density ratio will affect the contact angle, so the wall density ${\rho }_{\mathrm{w}}$ should be lightly tuned to ensure the contact angle unchanged. Numerical results of the slip length against the density ratio are shown in Figure 12. In the case of $\theta =148.7°$, the slip length increases from 0.18 to 4.54 with the liquid-gas density ratio from 22:1 to 293:1. The relationship between the slip length and the density ratio is close to linear. The above phenomena is induced by the strong changes of the liquid density near the wall and the formation of the “gas layer” at the solid surface. When the contact angles $\theta$ are 115.9°, 92.6°, 76.0°, and 61.9°, the slip lengths increase from 0.09 to 0.77, −0.09 to 0.12, −0.21 to −0.13, and −0.30 to −0.27, respectively, the growth slope is relatively smaller. Thus, the liquid-gas density ratio has an important influence on the fluid flow characteristics, especially for the hydrophobic wall.

6. Conclusions

Based on the improved pseudo-potential model, the MRT-LBM is proposed to investigate the fluid flow in the rock micro-fracture, and it is verified according to two benchmark cases. The slip length is used to evaluate the flow characteristics in the micro-fracture, and the effects of the contact angles, the driving pressure, and the liquid-gas density ratio on the slip length are discussed, the corresponding conclusions are summarized as follows:

(1)

With increasing contact angle, the slip length increases at the wall, and the larger the contact angle, the more obvious the slip length changes.

(2)

Under the laminar flow regime, the fluid flow velocity is proportional to the driving pressure, but there is almost no change in the slip length with the driving pressure increasing, so the driving pressure has almost no impact on the slip length.

(3)

The slip length increases with the increasing of the liquid-gas density ratio, and the larger the wall contact angle, the more remarkable it shows. The liquid-gas density ratio has an important influence on the fluid flow characteristics, especially for the hydrophobic wall.

(4)

The wall slip is induced by the changing of liquid density near the solid surface, which is attributed to the varying interaction force between fluid particles and the solid wall. With the decrease of interaction force, the liquid density near the wall decreases and the slip length increases.

And as the increase of interaction force, it shows no slip or even negative slip at the solid wall.