# The Seepage Model Considering Liquid/Solid Interaction in Confined Nanoscale Pores

Authors: Xiaona Cui; Erlong Yang; Kaoping Song; Yuming Wang
Journal: Geofluids
 (2018)

---

## Abstract

Different from conventional reservoirs, nanoscale pores and fractures are dominant in tight or shale reservoirs. The flow behaviors of hydrocarbons in nanopores (called “confined space”) are more complex than that of bulk spaces. The interaction between liquid hydrocarbons and solid pore wall cannot be neglected. The viscosity formula which is varied with the pore diameter and interaction coefficient of liquids and solids in confined nanopores has been introduced in this paper to describe the interaction effects of hydrocarbons and pore walls. Based on the Navier-Stokes equation, the governing equation considered liquid/solid effect in two dimensions has been established, and approximate theoretical solutions to the governing equations have been achieved after mathematic simplification. By introducing the vortex equation, the complex numerical seepage model has been discretized and solved. Numerical results show that the radial velocity distribution near the solid wall has an obvious change when considering the liquid/solid interaction. The results consist well with that approximate mathematical solution. And when the capillary radius is smaller, the liquid and solid interaction coefficientn is greater. The liquid and solid interaction obviously cannot be neglected in the seepage model if the capillary radius is small than 50 nm when n>0.1. The numerical model has also been further validated by two types of nanopore flow tests: from pore to throat and inversely from throat to pore. There is no big difference in flow regularity of throat to pore model considering when liquid/solid interaction or not, whereas the liquid/solid interaction of pore to throat model totally cannot be overlooked.


---

## Body

## 1. Introduction

Unconventional oil and gas are playing an increasingly important role in global energy supply. Unlike conventional reservoirs, tight or shale reservoirs possess unique characteristics, such as ultralow permeability and strong heterogeneity [1]. The length scales associated with transport phenomena in tight and shale formations are rich, from nanoscale phenomena to field-scale applications [2–4]. Nanoscale pores and nanoscale fractures are still dominant in tight or shale reservoirs. The extremely small size of tight or shale pores, leads to strong intermolecular forces between fluid particles and solid surface [5–7]. The flow behaviors of hydrocarbons in nanopores (called “confined space”) are more complex than those in bulk spaces [8]. The deviations in the prediction of the phase behavior and production in confined nanoscale tight or shale reservoirs from that of conventional oilfields are needed to be calibrated [9–11].In macroflow, fluid can be assumed as continuous medium due to its characteristic size of the liquid molecules is bigger than the mean free path. In microscale flow [12, 13], when the characteristic size and the mean free path are in the same magnitude, some macroscopic and laws based on the continuous medium would not work anymore [14]. The viscosity coefficient is also needed to be discussed again. Newtonian fluid presents the properties of non-Newtonian fluid in nanochannels [15], such as the boundary layer effect. The fluid in porous media can be classified into body fluid and boundary layer fluid. Body fluid which is close to the axial of porous media channel is not affected by the boundary effect. The fluid which is affected by boundary effect is easy to form a boundary layer on the pore wall [16].Many researches on the formation of boundary layer and its influence on porous media seepage have been carried out. Engelhardt [17] explains the reason the oil seepage deviating from Darcy’s law is mainly due to the active component adsorption on the surface of rock wall. Karimov [18] confirmed that the active surface substances in fluid could be adsorbed on rock particle surface by experiments. Secondly, tight rock such as shale and mudstone can adsorb the salt components in water. Bhasin [19] found that oil viscosity of boundary layer is 10–15 times higher than the body oil. The distribution of crude oil composition is orderly in the boundary layer, and the viscosity increased. The water presents non-Newton characteristics in the tiny pore under the effect of liquid-solid interface [20–22]. Lilatov pointed out that the boundary layer thickness of polymer solution is 160 nanons [23]. Zhang and Xiangan described the rheological properties of the oil film by transforming the constitutive Bingham equation, in view of the increasing viscosity and yield stress of oil film along the radial direction [24]. Hara and Tsuchiya examined water-(pyroclastic) rock interactions using flow-through experiments to deduce the effect of mass transport phenomena on the reaction process [25]. Bea et al. presented a reactive transport formulation coupled to thermohydrodynamic and simplified mechanical processes, taking the heat and mass transport and water-rock interactions into account [26].For the Newtonian fluid, exact solutions and transformations of equations as well as different problems on the hydrodynamic boundary layer have been addressed by many researchers. Pavlovskii began to obtain invariant solutions to the steady-state boundary layer equations in 1961 [27]. And invariant-irreducible, partially invariant solutions of steady-state boundary layer equations have been worked out [28]. Vereshchagina fibered the spatial unsteady boundary layer equations [29]. Then, many similarity solutions of the two-dimensional unsteady and steady-state boundary-layer equations came out [30, 31]. The reduction method of partial differential equations has been proposed to describe a laminar stationary flat boundary layer [32, 33]. Explicit solutions of the boundary-layer equations have been found [34–39]. Prandtl’s boundary layer equation for radial flow is derived [40]. General solution and quasiexact solution of equation is found, and norm of the residual function is presented. Equations of unsteady axisymmetric boundary layer are studied [41]. The solutions are obtained with a new method (direct method of functional separation of variables) based on using particular solutions to an auxiliary ODE [42]. Su et al. put forward a method based on a mixed control volume-discontinuous finite element formulation to accurately simulate multiphase flow in fractured shale reservoir. It has solved the problem of discontinuous or near discontinuous behavior of saturation in real oilfield [43].Previous researches have studied flow behavior in big channel considering the liquid and solid interaction, and many numerical reservoir simulation models in terms of shale or tight oil and gas have been investigated. Monteiro et al. completed mathematical modeling of the flow in nanoporous rocks based on the hypothesis that the permeability of the inclusions substantially depends on the pressure gradient; boundary layers are considered in this model [44]. Li applied engineering density functional theory (DFT) combined with the Peng-Robinson equation of state (EOS) to investigate the adsorption and phase behavior of pure substances and mixtures in nanopores [45]. Wang incorporated the capillary pressure effect and pore space compaction in a compositional reservoir simulator to evaluate their effects on production [46]. Mi et al. has established the DFN model with fracture network to study the effect of the stress-dependent fracture conductivity on the shale gas well performance [47, 48]. Wang et al. developed a new semianalytical model for multiple-fractured horizontal wells (MFHWs) with stimulated reservoir volume (SRV) in tight oil reservoirs by combining source function theory with boundary element idea [49]. However, there are no researches about the seepage description in confined nanopores considering the liquid and solid interaction. This paper establishes a two-dimensional numerical model of seepage flow in confined nanopores. Vortex equation has been applied to discretize the complex model. The numerical results are compared with the approximate mathematical results. The numerical model has also been validated by two types of nanopore flow: from pore to throat and inversely from throat to pore.


## 2. Liquid/Solid Interaction in Confined Pores

Molecular interactions include three parts, orientating effect, inducing effect, and dispersing effect. The force of water molecules in boundary layer which is affected by inner interaction action of water molecules and solid pore wall would influence the viscosity. Since viscosity can represent molecular attraction, when the molecular attraction is greater, the viscosity will be greater. According to this, we can assume viscosity is proportional to intermolecular attraction. The viscosity value in boundary layer is composed of two parts: bulk fluid viscosity and additional viscosity imposed by solid surface. In formula, the viscosity coefficient of fluid [50] in nanoflow boundary layer can be arranged as
(1)μ=μ0+ϕ′yn,(2)ϕ′=23μ12μ22kT+α1μ′22+α2μ′12+32α1α2I1I2I1+I2.From (1), on the solid surface where y→0, the viscosity coefficient of water molecule is infinity, and the water molecules would stand still, which is consistent with the classical boundary layer theory without slip conditions. While on the solid surface where y→∞, the water viscosity is just normal.


## 3. Mathematical Modeling Description

### 3.1. Governing Equation

Rheological equation of nanoscale porous media fluid in boundary layer has been revised based on the Navier-Stokes equation considering the viscosity coefficient formula of boundary layer fluid. Then, the governing equation of boundary layer fluid in the nanochannel has established.Velocity componentw∗ is equal to zero in two-dimensional flow; therefore all velocity variables has nothing to do with z∗, that is, its partial derivative is equal to zero. The continuity equation is as follows:
(3)∂ρ∂t+∂ρu∗∂x∗+∂ρv∗∂y∗=0.Momentum equations:(4a)∂u∗∂t+u∗∂u∗∂x∗+v∗∂u∗∂y∗=1ρ∂∂x∗−p∗+2μ∂u∗∂x∗−23μ∂u∗∂x∗+∂v∗∂y∗+∂∂y∗μ∂v∗∂x∗+∂u∗∂y∗,(4b)∂v∗∂t+u∗∂v∗∂x∗+v∗∂v∗∂y∗=1ρ∂∂x∗μ∂v∗∂x∗+∂u∗∂y∗+∂∂y∗−p∗+2μ∂v∗∂y∗−23μ∂u∗∂x∗+∂v∗∂y∗.For steady incompressible fluid, (3), (4a), and (4b) are as follows:
(5a)∂u∗∂x∗+∂v∗∂y∗=0,(5b)u∗∂u∗∂x∗+v∗∂u∗∂y∗=−1ρ∂p∗∂x∗+μρ∂2u∗∂x∗2+∂2u∗∂y∗2,(5c)u∗∂v∗∂x∗+v∗∂v∗∂y∗=−1ρ∂p∗∂v∗+μρ∂2v∗∂x∗2+∂2v∗∂y∗2,where μ=μ0+ϕ′/y∗n, (5a), (5b), and (5c) are the fluid flow control equations considering solid-liquid interaction.


### 3.2. Approximate Solution

If the effects of viscosity and heat diffusion in boundary layer simply react in thin layer close to the surface, the change of velocity and temperature mainly exists in the thin layer of this object plane. These two kinds of thin layer thickness can be expressed as velocity boundary layer thicknessδV and temperature boundary layer thickness δT. Compared with x∗, they are so small that the order quantity is δ, so quantity order y∗ within the boundary layer is δ. Velocity component u changes from the speed 0 in the boundary layer to U the mainstream speed beyond boundary layer. Taking x∗ and U as characteristic values, the order quantity is one. At the same time, the continuity equation can be approximatively written as
(6)∂ρ∂t+∂ρu∗∂x∗+∂ρv∗∂y∗=0.Of the viscous force and inertial force item in (4a), truncating the order items with small quantity, respectively, it can be obtained:
(7)∂u∗∂t+u∗∂u∗∂x∗+v∗∂u∗∂y∗=1ρ−∂p∗∂x∗+∂∂y∗μ∂u∗∂y∗.According to (7), if the viscous force and inertial force items are in the same amount of order, the order quantity of μ/ρ should be δ2, so it is still the highest item in (4b):
(8)∂p∗∂y∗=0.For two-dimension nanoflow which is steady and incompressible in boundary layer, when substituting (1) (μ=μ0+ϕ′/y∗n) into (7), the simplified continuity equation and momentum equation can be obtained as follows:
(9)∂u∗∂x∗+∂v∗∂y∗=0.(10)u∗∂u∗∂x∗+v∗+nϕ′y∗n+1∂u∗∂y∗=−1ρ∂p∗∂x∗+1ρμ0+ϕ′y∗n∂2u∗∂y∗2.Boundary conditions arey∗=0,u∗=v∗=0;y∗=∞,u∗=Ux∗, (9) and (10) are approximation equations of two-dimension nanoflow which is steady and incompressible.Taking the flow of Newtonian fluid in circular cross capillary tube as an example, the velocity distribution considering liquid-solid effect and not have been obtained by solving the approximation equations.As shown in Figure1, take the tube axes as x∗ axis, radial coordinates as y∗ axis, both the circumferential and radial velocity is zero, horizontal velocity component is u∗ which is the only function of y∗, and the pressure on each cross section is a constant value. When nanoflow boundary layer equations are expressed in cylindrical coordinates, (9) and (10) can be simplified to:
(11)dp∗dx∗=1y∗ddy∗μ0+ϕ′R∗−y∗ny∗du∗dy∗.Figure 1
Flow in circular section of capillary.Boundary conditions arey∗=R∗,u∗=0.For capillary, the diameter is so small that the fluid flow is severely affected by tube surface and the viscosity of fluid itself, so (11) can be simplified to
(12)dp∗dx∗=1y∗ddy∗ϕ′R∗−y∗ny∗du∗dy∗.Integrating (11), combined with the boundary conditions, the formula of fluid velocity distribution within the capillary can be calculated as follows:
(13)u∗=−12ϕ′1n+1n+2R∗−y∗n+1n+1y∗+R∗dp∗dx∗.Wheny∗ is equal to zero, the maximum speed in capillary is as follows:
(14)u∗max=−12ϕ′1n+1n+2R∗n+2dp∗dx∗.The flow rate in capillary is:(15)Q=−1ϕ′πn+1n+21n+2n+3+2n+1n+2n+3n+4R∗n+4dp∗dx∗,where ϕ′=μ0 and n=0, and (15) is just the Poiseuille formula.Figures2 and 3 are velocity distribution curves when liquid-solid interaction coefficient is ϕ′=μ0, pressure gradient is dp∗/dx∗=1Pa/m, and capillary radius R∗=50μm, considering the solid wall gravitation to water molecules. The figure of velocity distribution curve in Figure 3 close to solid wall is concave, explains that the fluid viscosity increased under the interaction of solid and liquid, and when the additional viscosity close to the wall is greater, the velocity gradient is smaller. While the figure of velocity distribution curves is convex without considering the solid wall gravitation. The effect of solid molecules near the surface to liquid is not obvious; the fluid viscosity is relatively small, and the velocity gradient is relatively big. In addition, fluid flow rate decreases significantly with the increase of liquid-solid interaction index n as shown in Figure 4. That is to say, the effect of solid molecules to liquid cannot be ignored especially in the nanochannel.Figure 2
Velocity distribution curves in the channel considering the solid-liquid effect (n≥1). (R∗=50μm, dp∗/dx∗=1Pa/m, μ0=0.0006Pa.s, ϕ′=μ0).Figure 3
Velocity distribution curves in the channel considering the solid-liquid effect (n<1). (R∗=50μm, dp∗/dx∗=1Pa/m, μ0=0.0006Pa.s, ϕ′=μ0).Figure 4
Velocity distribution curves in the channel considering the solid-liquid effect. (R∗=0.05μm, dp∗/dx∗=2500Pa/m, μ0=0.0006Pa.s, ϕ′=μ0).Figure5 shows the relationship between the average velocity ratio of fluid Vn (which considers solid-liquid interaction) and V0 (not considering the solid-liquid interaction) in different capillary radiuses and solid-liquid index n. It can be seen from the curves, with the increase of solid-liquid index n, the fluid velocity exponentially decreases, and when the capillary radius is smaller, the reduction gets bigger. On the condition of same solid-liquid index n, small capillary radius would result in obvious effect of solid-liquid reaction. In addition, when the liquid and solid interaction coefficient n is bigger than 0.1, the liquid and solid interaction obviously cannot be neglected in the flow calculation if the capillary radius is small than 0.05 μm.Figure 5
Relationship of average velocity ratio (Vn/V0) and the solid-liquid index n in the capillary.


## 3.1. Governing Equation

Rheological equation of nanoscale porous media fluid in boundary layer has been revised based on the Navier-Stokes equation considering the viscosity coefficient formula of boundary layer fluid. Then, the governing equation of boundary layer fluid in the nanochannel has established.Velocity componentw∗ is equal to zero in two-dimensional flow; therefore all velocity variables has nothing to do with z∗, that is, its partial derivative is equal to zero. The continuity equation is as follows:
(3)∂ρ∂t+∂ρu∗∂x∗+∂ρv∗∂y∗=0.Momentum equations:(4a)∂u∗∂t+u∗∂u∗∂x∗+v∗∂u∗∂y∗=1ρ∂∂x∗−p∗+2μ∂u∗∂x∗−23μ∂u∗∂x∗+∂v∗∂y∗+∂∂y∗μ∂v∗∂x∗+∂u∗∂y∗,(4b)∂v∗∂t+u∗∂v∗∂x∗+v∗∂v∗∂y∗=1ρ∂∂x∗μ∂v∗∂x∗+∂u∗∂y∗+∂∂y∗−p∗+2μ∂v∗∂y∗−23μ∂u∗∂x∗+∂v∗∂y∗.For steady incompressible fluid, (3), (4a), and (4b) are as follows:
(5a)∂u∗∂x∗+∂v∗∂y∗=0,(5b)u∗∂u∗∂x∗+v∗∂u∗∂y∗=−1ρ∂p∗∂x∗+μρ∂2u∗∂x∗2+∂2u∗∂y∗2,(5c)u∗∂v∗∂x∗+v∗∂v∗∂y∗=−1ρ∂p∗∂v∗+μρ∂2v∗∂x∗2+∂2v∗∂y∗2,where μ=μ0+ϕ′/y∗n, (5a), (5b), and (5c) are the fluid flow control equations considering solid-liquid interaction.


## 3.2. Approximate Solution

If the effects of viscosity and heat diffusion in boundary layer simply react in thin layer close to the surface, the change of velocity and temperature mainly exists in the thin layer of this object plane. These two kinds of thin layer thickness can be expressed as velocity boundary layer thicknessδV and temperature boundary layer thickness δT. Compared with x∗, they are so small that the order quantity is δ, so quantity order y∗ within the boundary layer is δ. Velocity component u changes from the speed 0 in the boundary layer to U the mainstream speed beyond boundary layer. Taking x∗ and U as characteristic values, the order quantity is one. At the same time, the continuity equation can be approximatively written as
(6)∂ρ∂t+∂ρu∗∂x∗+∂ρv∗∂y∗=0.Of the viscous force and inertial force item in (4a), truncating the order items with small quantity, respectively, it can be obtained:
(7)∂u∗∂t+u∗∂u∗∂x∗+v∗∂u∗∂y∗=1ρ−∂p∗∂x∗+∂∂y∗μ∂u∗∂y∗.According to (7), if the viscous force and inertial force items are in the same amount of order, the order quantity of μ/ρ should be δ2, so it is still the highest item in (4b):
(8)∂p∗∂y∗=0.For two-dimension nanoflow which is steady and incompressible in boundary layer, when substituting (1) (μ=μ0+ϕ′/y∗n) into (7), the simplified continuity equation and momentum equation can be obtained as follows:
(9)∂u∗∂x∗+∂v∗∂y∗=0.(10)u∗∂u∗∂x∗+v∗+nϕ′y∗n+1∂u∗∂y∗=−1ρ∂p∗∂x∗+1ρμ0+ϕ′y∗n∂2u∗∂y∗2.Boundary conditions arey∗=0,u∗=v∗=0;y∗=∞,u∗=Ux∗, (9) and (10) are approximation equations of two-dimension nanoflow which is steady and incompressible.Taking the flow of Newtonian fluid in circular cross capillary tube as an example, the velocity distribution considering liquid-solid effect and not have been obtained by solving the approximation equations.As shown in Figure1, take the tube axes as x∗ axis, radial coordinates as y∗ axis, both the circumferential and radial velocity is zero, horizontal velocity component is u∗ which is the only function of y∗, and the pressure on each cross section is a constant value. When nanoflow boundary layer equations are expressed in cylindrical coordinates, (9) and (10) can be simplified to:
(11)dp∗dx∗=1y∗ddy∗μ0+ϕ′R∗−y∗ny∗du∗dy∗.Figure 1
Flow in circular section of capillary.Boundary conditions arey∗=R∗,u∗=0.For capillary, the diameter is so small that the fluid flow is severely affected by tube surface and the viscosity of fluid itself, so (11) can be simplified to
(12)dp∗dx∗=1y∗ddy∗ϕ′R∗−y∗ny∗du∗dy∗.Integrating (11), combined with the boundary conditions, the formula of fluid velocity distribution within the capillary can be calculated as follows:
(13)u∗=−12ϕ′1n+1n+2R∗−y∗n+1n+1y∗+R∗dp∗dx∗.Wheny∗ is equal to zero, the maximum speed in capillary is as follows:
(14)u∗max=−12ϕ′1n+1n+2R∗n+2dp∗dx∗.The flow rate in capillary is:(15)Q=−1ϕ′πn+1n+21n+2n+3+2n+1n+2n+3n+4R∗n+4dp∗dx∗,where ϕ′=μ0 and n=0, and (15) is just the Poiseuille formula.Figures2 and 3 are velocity distribution curves when liquid-solid interaction coefficient is ϕ′=μ0, pressure gradient is dp∗/dx∗=1Pa/m, and capillary radius R∗=50μm, considering the solid wall gravitation to water molecules. The figure of velocity distribution curve in Figure 3 close to solid wall is concave, explains that the fluid viscosity increased under the interaction of solid and liquid, and when the additional viscosity close to the wall is greater, the velocity gradient is smaller. While the figure of velocity distribution curves is convex without considering the solid wall gravitation. The effect of solid molecules near the surface to liquid is not obvious; the fluid viscosity is relatively small, and the velocity gradient is relatively big. In addition, fluid flow rate decreases significantly with the increase of liquid-solid interaction index n as shown in Figure 4. That is to say, the effect of solid molecules to liquid cannot be ignored especially in the nanochannel.Figure 2
Velocity distribution curves in the channel considering the solid-liquid effect (n≥1). (R∗=50μm, dp∗/dx∗=1Pa/m, μ0=0.0006Pa.s, ϕ′=μ0).Figure 3
Velocity distribution curves in the channel considering the solid-liquid effect (n<1). (R∗=50μm, dp∗/dx∗=1Pa/m, μ0=0.0006Pa.s, ϕ′=μ0).Figure 4
Velocity distribution curves in the channel considering the solid-liquid effect. (R∗=0.05μm, dp∗/dx∗=2500Pa/m, μ0=0.0006Pa.s, ϕ′=μ0).Figure5 shows the relationship between the average velocity ratio of fluid Vn (which considers solid-liquid interaction) and V0 (not considering the solid-liquid interaction) in different capillary radiuses and solid-liquid index n. It can be seen from the curves, with the increase of solid-liquid index n, the fluid velocity exponentially decreases, and when the capillary radius is smaller, the reduction gets bigger. On the condition of same solid-liquid index n, small capillary radius would result in obvious effect of solid-liquid reaction. In addition, when the liquid and solid interaction coefficient n is bigger than 0.1, the liquid and solid interaction obviously cannot be neglected in the flow calculation if the capillary radius is small than 0.05 μm.Figure 5
Relationship of average velocity ratio (Vn/V0) and the solid-liquid index n in the capillary.


## 4. Numerical Modeling and Solutions

### 4.1. Numerical Equation

#### 4.1.1. Constitutive Equation

Constitutive equation of Newtonian fluid in the form of a tensor is as follows:(16)τij∗=2μεij∗+−p∗+λ∇⋅V∗δij=μ∂ui∗∂xj∗+∂uj∗∂xi∗+−p∗+λ∂uk∗∂xk∗δij,(17)δij=1i=j0i≠j.Under Cartesian coordinate system, the component forms ofτij and εij are as follows:
(18)τ∗=τxx∗τxy∗τxz∗τyx∗τyy∗τyz∗τzx∗τzy∗τzz∗,ε∗=εxx∗·εxy∗εxz∗εyx∗εyy∗·εyz∗εzx∗εzy∗εzz∗.The component form of incompressible Newtonian fluid constitutive equation (16) under the Cartesian coordinate system can be written as
(19)τxx∗=2μ∂u∗∂x∗,τyy∗=2μ∂v∗∂y∗,τzz∗=2μ∂w∗∂z∗,(20)τxy∗=τyx∗=μ∂v∗∂x∗+∂u∗∂y∗,τyz∗=τzy∗=μ∂w∗∂y∗+∂v∗∂z∗,τzx∗=τxz∗=μ∂u∗∂z∗+∂w∗∂x∗.In the two-dimension flow,τyy∗=τyz∗=τzy∗=τzx∗=τxz∗=0.Fluid constitutive equations considering liquid-solid interaction: the apparent performance of solid wall gravitation to fluid is the increase of fluid viscosity in the reservoir nanoporous media flow. Resulting in changes of fluid rheological properties, the linear relationship between fluid viscosity and strain rate is no longer existed, that is, the fluid viscosity is not a constant value but is the function of the distance to the solid wall. Substitute (4a) and (4b) into (20); 2D boundary layer fluid constitutive equations can be obtained as follows:
(21a)τxx∗=2μ+ϕ′y∗n∂u∗∂x∗,(21b)τyy∗=2μ+ϕ′y∗n∂v∗∂y∗,(21c)τxy∗=τyx∗=μ+ϕ′y∗n∂v∗∂x∗+∂u∗∂y∗.


#### 4.1.2. Dimensionless of Governing Equation

Due to the orders of magnitude of the characteristic value of flow parameters in the governing equation is generally within the ranges of nanoscale, the governing equation can be non-dimensionalized in order to proceed calculations more easily and accurately. Therefore, introducing the characteristic length and the characteristic velocity, the dimensionless parameters of equation can be written as follows:(22)x=x∗L,y=y∗L,u=u∗U,v=v∗U,τ=Lη0Uτ∗,ϕ′=1η0Lnϕ′∗,p=LμUp∗.It can be revised as "The characteristic length and characterisic velocity can be defined as 1μm and 1 μm/s respectively in the calculation. After dimensionless, the governing equation is turned into continuity equation (mass conservation equation):
(23)∂u∂x+∂v∂y=0.Dynamic equation (momentum conservation equation):(24a)Reu∂u∂x+v∂u∂y=−∂p∂x+∂τxx∂x+∂τxy∂y,(24b)Reu∂v∂x+v∂v∂y=−∂p∂y+∂τxy∂x+∂τyy∂y.Constitutive equation (rheological equation):(25a)τxx=21+ϕ′yn∂u∂x,(25b)τyy=21+φ′yn∂v∂y,(25c)τxy=τyx=1+ϕ′yn∂v∂x+∂u∂y.


### 4.2. Vortex Equation

In order to reduce the number of unknowns in the solving procedure, derivation ofy and x in momentum (24a) and (24b), respectively, can be taken to eliminate the pressure term in the momentum equation; the results can be transformed as follows:
(26a)Reu∂2u∂x∂y+∂u∂x∂u∂y+∂v∂y∂u∂y+v∂2u∂y2=−∂2p∂x∂y+∂2τxx∂x∂y+∂2τxy∂y2,(26b)Reu∂2v∂x2+∂u∂x∂v∂x+∂v∂x∂v∂y+v∂2v∂x∂y=−∂2p∂x∂y+∂2τxy∂x2+∂2τyy∂x∂y.Subtracting (26b) from (26a), the results can be further simplified to
(27)Reu∂2u∂x∂y−∂2v∂x2+v∂2u∂y2−∂2v∂x∂y=∂2τxx∂x∂y+∂2τxy∂y2−∂2τxy∂x2−∂2τyy∂x∂y.The introduction of stream function and vorticity:(28)u=∂ψ∂y,v=−∂ψ∂x,Ω=∂v∂x−∂u∂y.There is one curve AB between the two streamlines as shown in Figure6, and then the volume flow of unit thickness is
(29)qV=∫y1y2udy+∫x1x2v−dx=∫y1y2∂ψ∂ydy+∫x1x2∂ψ∂xdx=∫x1,y1x2,y2dψ=ψ2−ψ1=∫x1,y1x2,y2dψ=ψ2−ψ1.Figure 6
Sketch map of the planar flow field and flow line.The stream function equation:(30)∂2ψ∂x2+∂2ψ∂y2=−Ω.Momentum equation can be written into vorticity equation in the following form:(31)Re∂ψ∂x∂Ω∂y−∂ψ∂y∂Ω∂x=∂2τxx∂x∂y+∂2τxy∂y2−∂2τxy∂x2−∂2τyy∂x∂y.After finishing, oval vorticity equation can be obtained as follows:(32)∂2Ω∂x2+∂2Ω∂y2+Re∂ψ∂x∂Ω∂y−∂ψ∂y∂Ω∂x=∂2τxx∂x∂y+∂2τxy∂y2−∂2τxy∂x2−∂2τyy∂x∂y+∂3v∂x3−∂3u∂x2∂y+∂3v∂y2∂x−∂3u∂y3.Since the above equation has third-order derivative, it is difficult to determine the boundary conditions when applying difference methods to get the solution. Additional stress tensor caused by non-Newtonian fluid can be introduced to solve this problem, and the additional stress terms of Newtonian fluid is zero.(33)S=τ−L+LT,where
(34)L=∂u∂x∂v∂x0∂u∂y∂v∂y0000,LT=∂u∂x∂u∂y0∂v∂x∂v∂y0000,τ=τxxτxy0τxyτyy0000.Therefore, (32) can be written as
(35)F=∂2∂x∂yτxx−2∂u∂x−τyy−2∂v∂y+∂2∂y2−∂2∂x2τxy−∂u∂y+∂v∂x=∂2∂x∂ySxx−Syy+∂2∂y2−∂2∂x2Sxy.Vorticity equation (32) can be written as
(36)∂2Ω∂x2+∂2Ω∂y2+Re∂ψ∂x∂Ω∂y−∂ψ∂y∂Ω∂x=∂2∂x∂ySxx−Syy+∂2∂y2−∂2∂x2Sxy.


### 4.3. Discretization of Governing Equation

#### 4.3.1. Discretization of Stream Function

Stream function equation is the elliptic equations; it can be dispersed by arrays of central difference with second-order accuracy. The main diagonal of the discrete equation coefficient is dominant, which guarantee the convergence of algorithm, the discrete equation is:(37)ψi+1,j−2ψi,j+ψi−1,jΔx2+ψi,j+1−2ψi,j+ψi,j−1Δy2=−Ωi,j,(38)ψn+1i,j=121+βψni+1,j+ψni−1,j+βψni,j+1+βψni,j−1+Δx2Ωi,j,where Δx,Δy is the grid step length of x and y direction, β=Δx2/Δy2; n is the number of iterations.


#### 4.3.2. Discretization of Vorticity Function

The left item of vorticity equation (38) applies the first kind of windward array of difference. That is, when coefficient of the first derivative of vorticity is greater than zero, then discrete the equation by array of backward difference; when the coefficient is less than or equal to zero, then discrete the equation by array of forward difference. Although difference equation has only first-order accuracy, this array of difference is equipped with main diagonally dominant unconditionally, thus the stability of numerical calculation can be ensured. The right item of vorticity equation (38) applies central difference methods. The discrete form is as follows:
(39)k0Ωi,j+k1Ωi+1,j+k2Ωi−1,j+k3Ωi,j+1+k4Ωi,j−1=14ΔxΔySi+1,j+1xx−Si−1,j+1xx−Si+1,j−1xx+Si−1,j−1xx−Si+1,j+1yy−Si−1,j+1yy−Si+1,j−1yy+Si−1,j−1yy+1Δy2Si,j+1xy−2Si,jxy+Si,j−1xy.The coefficients in the equation can be expressed as(40)b1=Reψi+1,j−ψi−1,j2Δx,b2=−Reψi,j+1−ψi,j−12Δy,k0=−2Δx2−2Δy2−b1Δy−b2Δx,k1=1Δx2+b22Δx+b22Δx,k2=1Δx2+b22Δx−b22Δx,k3=1Δy2+b12Δy+b12Δy,k4=1Δy2+b12Δy−b12Δy.


#### 4.3.3. Discretization of Constitutive Function

The first derivative of constitutive equation in (21a), (21b), and (21c) applies first-order array of difference, and the discrete equations are as follows:
(41a)τxxi,j=2ui,j−ui−1,jΔx,(41b)τyyi,j=2vi,j−vi,j−1Δy,(41c)τxyi,j=τyxi,j=vi,j−vi−1,jΔx+ui,j−ui,j−1Δy.


### 4.4. Boundary Conditions

#### 4.4.1. Velocity

In this paper, the flow velocity distribution at the nanochannel can be determined by (13) considering effect of solid and liquid. After dimensionless, it can be written as follows:
(42)u=−12φ′1n+1n+2R−yn+1n+1y+Rdpdx.


#### 4.4.2. Velocity Derivative

According to the flow characteristic of Poiseuille, the constitutive equation and the flow velocity distribution (4), the first-order derivative in the channel can be written as:
(43)∂u∂y=1ϕ′∂p∂xR−yn2y.


#### 4.4.3. Stress Value

(44)τxx=0,τxy=1+ϕyn1ϕ′∂p∂xR−yn2y,τyy=0.Substitute (44) into constitutive equation and partial stress equation, the boundary value of additional stress S are as follows:
(45)Sxx=Syy=0Sxy=∂p∂xR−yn2yn−1.


#### 4.4.4.ψ and Ω

Based on the definition of vorticity and stream function, the distribution of stream function and vorticity at the nanochannel can be determined by velocity distribution as follows:(46)ψ=∂p∂x12ϕ′n+1n+2n+3n+1y+2RR−yn+2−2Rn+3,Ω=−1ϕ∂p∂xR−yn2y.


### 4.5. Sample Calculation and Results Comparison

Numerical calculation of the capillary flow has been proceeded in the capillary with radius 5μm, the solid-liquid effect index n is 0, 0.1, and 0.2, respectively, based on the established 2D numerical boundary layer equations. Radial velocity distribution curves in the capillary radius of 5 μm are showed in Figures 7 and 8.Figure 7
Velocity distribution at different solid-liquid index in capillary profiles (n=0,0.1).Figure 8
Velocity distribution at different solid-liquid index in capillary profiles (n=0.1,0.2).The interaction effect between solid and liquid phases in nanocapillary can be obviously seen by comparing two curves at different solid-liquid index in Figure7. When n is zero, the velocity distribution curve presents the convex type, the flow characteristics of the Newtonian fluid. When n is 0.1 or 0.2, the curve near the wall presents the concave type, the flow characteristics of the non-Newtonian fluid, which illustrates that the interaction effect is stronger close to wall surface. And the velocity is getting smaller and the increase amplitude is smaller from pore wall to the capillary center, that is to say, fluid viscosity coefficient has big changes close to the wall surface. Figures 7 and 8 are obtained by numerical calculation, the shape of curves and indicated results are in consistent with the approximate solutions obtained by the governing equation. Results proved that the numerical seepage model established in this paper can fulfill the precision and requirements of flow in confined nanopores.


## 4.1. Numerical Equation

### 4.1.1. Constitutive Equation

Constitutive equation of Newtonian fluid in the form of a tensor is as follows:(16)τij∗=2μεij∗+−p∗+λ∇⋅V∗δij=μ∂ui∗∂xj∗+∂uj∗∂xi∗+−p∗+λ∂uk∗∂xk∗δij,(17)δij=1i=j0i≠j.Under Cartesian coordinate system, the component forms ofτij and εij are as follows:
(18)τ∗=τxx∗τxy∗τxz∗τyx∗τyy∗τyz∗τzx∗τzy∗τzz∗,ε∗=εxx∗·εxy∗εxz∗εyx∗εyy∗·εyz∗εzx∗εzy∗εzz∗.The component form of incompressible Newtonian fluid constitutive equation (16) under the Cartesian coordinate system can be written as
(19)τxx∗=2μ∂u∗∂x∗,τyy∗=2μ∂v∗∂y∗,τzz∗=2μ∂w∗∂z∗,(20)τxy∗=τyx∗=μ∂v∗∂x∗+∂u∗∂y∗,τyz∗=τzy∗=μ∂w∗∂y∗+∂v∗∂z∗,τzx∗=τxz∗=μ∂u∗∂z∗+∂w∗∂x∗.In the two-dimension flow,τyy∗=τyz∗=τzy∗=τzx∗=τxz∗=0.Fluid constitutive equations considering liquid-solid interaction: the apparent performance of solid wall gravitation to fluid is the increase of fluid viscosity in the reservoir nanoporous media flow. Resulting in changes of fluid rheological properties, the linear relationship between fluid viscosity and strain rate is no longer existed, that is, the fluid viscosity is not a constant value but is the function of the distance to the solid wall. Substitute (4a) and (4b) into (20); 2D boundary layer fluid constitutive equations can be obtained as follows:
(21a)τxx∗=2μ+ϕ′y∗n∂u∗∂x∗,(21b)τyy∗=2μ+ϕ′y∗n∂v∗∂y∗,(21c)τxy∗=τyx∗=μ+ϕ′y∗n∂v∗∂x∗+∂u∗∂y∗.


### 4.1.2. Dimensionless of Governing Equation

Due to the orders of magnitude of the characteristic value of flow parameters in the governing equation is generally within the ranges of nanoscale, the governing equation can be non-dimensionalized in order to proceed calculations more easily and accurately. Therefore, introducing the characteristic length and the characteristic velocity, the dimensionless parameters of equation can be written as follows:(22)x=x∗L,y=y∗L,u=u∗U,v=v∗U,τ=Lη0Uτ∗,ϕ′=1η0Lnϕ′∗,p=LμUp∗.It can be revised as "The characteristic length and characterisic velocity can be defined as 1μm and 1 μm/s respectively in the calculation. After dimensionless, the governing equation is turned into continuity equation (mass conservation equation):
(23)∂u∂x+∂v∂y=0.Dynamic equation (momentum conservation equation):(24a)Reu∂u∂x+v∂u∂y=−∂p∂x+∂τxx∂x+∂τxy∂y,(24b)Reu∂v∂x+v∂v∂y=−∂p∂y+∂τxy∂x+∂τyy∂y.Constitutive equation (rheological equation):(25a)τxx=21+ϕ′yn∂u∂x,(25b)τyy=21+φ′yn∂v∂y,(25c)τxy=τyx=1+ϕ′yn∂v∂x+∂u∂y.


## 4.1.1. Constitutive Equation

Constitutive equation of Newtonian fluid in the form of a tensor is as follows:(16)τij∗=2μεij∗+−p∗+λ∇⋅V∗δij=μ∂ui∗∂xj∗+∂uj∗∂xi∗+−p∗+λ∂uk∗∂xk∗δij,(17)δij=1i=j0i≠j.Under Cartesian coordinate system, the component forms ofτij and εij are as follows:
(18)τ∗=τxx∗τxy∗τxz∗τyx∗τyy∗τyz∗τzx∗τzy∗τzz∗,ε∗=εxx∗·εxy∗εxz∗εyx∗εyy∗·εyz∗εzx∗εzy∗εzz∗.The component form of incompressible Newtonian fluid constitutive equation (16) under the Cartesian coordinate system can be written as
(19)τxx∗=2μ∂u∗∂x∗,τyy∗=2μ∂v∗∂y∗,τzz∗=2μ∂w∗∂z∗,(20)τxy∗=τyx∗=μ∂v∗∂x∗+∂u∗∂y∗,τyz∗=τzy∗=μ∂w∗∂y∗+∂v∗∂z∗,τzx∗=τxz∗=μ∂u∗∂z∗+∂w∗∂x∗.In the two-dimension flow,τyy∗=τyz∗=τzy∗=τzx∗=τxz∗=0.Fluid constitutive equations considering liquid-solid interaction: the apparent performance of solid wall gravitation to fluid is the increase of fluid viscosity in the reservoir nanoporous media flow. Resulting in changes of fluid rheological properties, the linear relationship between fluid viscosity and strain rate is no longer existed, that is, the fluid viscosity is not a constant value but is the function of the distance to the solid wall. Substitute (4a) and (4b) into (20); 2D boundary layer fluid constitutive equations can be obtained as follows:
(21a)τxx∗=2μ+ϕ′y∗n∂u∗∂x∗,(21b)τyy∗=2μ+ϕ′y∗n∂v∗∂y∗,(21c)τxy∗=τyx∗=μ+ϕ′y∗n∂v∗∂x∗+∂u∗∂y∗.


## 4.1.2. Dimensionless of Governing Equation

Due to the orders of magnitude of the characteristic value of flow parameters in the governing equation is generally within the ranges of nanoscale, the governing equation can be non-dimensionalized in order to proceed calculations more easily and accurately. Therefore, introducing the characteristic length and the characteristic velocity, the dimensionless parameters of equation can be written as follows:(22)x=x∗L,y=y∗L,u=u∗U,v=v∗U,τ=Lη0Uτ∗,ϕ′=1η0Lnϕ′∗,p=LμUp∗.It can be revised as "The characteristic length and characterisic velocity can be defined as 1μm and 1 μm/s respectively in the calculation. After dimensionless, the governing equation is turned into continuity equation (mass conservation equation):
(23)∂u∂x+∂v∂y=0.Dynamic equation (momentum conservation equation):(24a)Reu∂u∂x+v∂u∂y=−∂p∂x+∂τxx∂x+∂τxy∂y,(24b)Reu∂v∂x+v∂v∂y=−∂p∂y+∂τxy∂x+∂τyy∂y.Constitutive equation (rheological equation):(25a)τxx=21+ϕ′yn∂u∂x,(25b)τyy=21+φ′yn∂v∂y,(25c)τxy=τyx=1+ϕ′yn∂v∂x+∂u∂y.


## 4.2. Vortex Equation

In order to reduce the number of unknowns in the solving procedure, derivation ofy and x in momentum (24a) and (24b), respectively, can be taken to eliminate the pressure term in the momentum equation; the results can be transformed as follows:
(26a)Reu∂2u∂x∂y+∂u∂x∂u∂y+∂v∂y∂u∂y+v∂2u∂y2=−∂2p∂x∂y+∂2τxx∂x∂y+∂2τxy∂y2,(26b)Reu∂2v∂x2+∂u∂x∂v∂x+∂v∂x∂v∂y+v∂2v∂x∂y=−∂2p∂x∂y+∂2τxy∂x2+∂2τyy∂x∂y.Subtracting (26b) from (26a), the results can be further simplified to
(27)Reu∂2u∂x∂y−∂2v∂x2+v∂2u∂y2−∂2v∂x∂y=∂2τxx∂x∂y+∂2τxy∂y2−∂2τxy∂x2−∂2τyy∂x∂y.The introduction of stream function and vorticity:(28)u=∂ψ∂y,v=−∂ψ∂x,Ω=∂v∂x−∂u∂y.There is one curve AB between the two streamlines as shown in Figure6, and then the volume flow of unit thickness is
(29)qV=∫y1y2udy+∫x1x2v−dx=∫y1y2∂ψ∂ydy+∫x1x2∂ψ∂xdx=∫x1,y1x2,y2dψ=ψ2−ψ1=∫x1,y1x2,y2dψ=ψ2−ψ1.Figure 6
Sketch map of the planar flow field and flow line.The stream function equation:(30)∂2ψ∂x2+∂2ψ∂y2=−Ω.Momentum equation can be written into vorticity equation in the following form:(31)Re∂ψ∂x∂Ω∂y−∂ψ∂y∂Ω∂x=∂2τxx∂x∂y+∂2τxy∂y2−∂2τxy∂x2−∂2τyy∂x∂y.After finishing, oval vorticity equation can be obtained as follows:(32)∂2Ω∂x2+∂2Ω∂y2+Re∂ψ∂x∂Ω∂y−∂ψ∂y∂Ω∂x=∂2τxx∂x∂y+∂2τxy∂y2−∂2τxy∂x2−∂2τyy∂x∂y+∂3v∂x3−∂3u∂x2∂y+∂3v∂y2∂x−∂3u∂y3.Since the above equation has third-order derivative, it is difficult to determine the boundary conditions when applying difference methods to get the solution. Additional stress tensor caused by non-Newtonian fluid can be introduced to solve this problem, and the additional stress terms of Newtonian fluid is zero.(33)S=τ−L+LT,where
(34)L=∂u∂x∂v∂x0∂u∂y∂v∂y0000,LT=∂u∂x∂u∂y0∂v∂x∂v∂y0000,τ=τxxτxy0τxyτyy0000.Therefore, (32) can be written as
(35)F=∂2∂x∂yτxx−2∂u∂x−τyy−2∂v∂y+∂2∂y2−∂2∂x2τxy−∂u∂y+∂v∂x=∂2∂x∂ySxx−Syy+∂2∂y2−∂2∂x2Sxy.Vorticity equation (32) can be written as
(36)∂2Ω∂x2+∂2Ω∂y2+Re∂ψ∂x∂Ω∂y−∂ψ∂y∂Ω∂x=∂2∂x∂ySxx−Syy+∂2∂y2−∂2∂x2Sxy.


## 4.3. Discretization of Governing Equation

### 4.3.1. Discretization of Stream Function

Stream function equation is the elliptic equations; it can be dispersed by arrays of central difference with second-order accuracy. The main diagonal of the discrete equation coefficient is dominant, which guarantee the convergence of algorithm, the discrete equation is:(37)ψi+1,j−2ψi,j+ψi−1,jΔx2+ψi,j+1−2ψi,j+ψi,j−1Δy2=−Ωi,j,(38)ψn+1i,j=121+βψni+1,j+ψni−1,j+βψni,j+1+βψni,j−1+Δx2Ωi,j,where Δx,Δy is the grid step length of x and y direction, β=Δx2/Δy2; n is the number of iterations.


### 4.3.2. Discretization of Vorticity Function

The left item of vorticity equation (38) applies the first kind of windward array of difference. That is, when coefficient of the first derivative of vorticity is greater than zero, then discrete the equation by array of backward difference; when the coefficient is less than or equal to zero, then discrete the equation by array of forward difference. Although difference equation has only first-order accuracy, this array of difference is equipped with main diagonally dominant unconditionally, thus the stability of numerical calculation can be ensured. The right item of vorticity equation (38) applies central difference methods. The discrete form is as follows:
(39)k0Ωi,j+k1Ωi+1,j+k2Ωi−1,j+k3Ωi,j+1+k4Ωi,j−1=14ΔxΔySi+1,j+1xx−Si−1,j+1xx−Si+1,j−1xx+Si−1,j−1xx−Si+1,j+1yy−Si−1,j+1yy−Si+1,j−1yy+Si−1,j−1yy+1Δy2Si,j+1xy−2Si,jxy+Si,j−1xy.The coefficients in the equation can be expressed as(40)b1=Reψi+1,j−ψi−1,j2Δx,b2=−Reψi,j+1−ψi,j−12Δy,k0=−2Δx2−2Δy2−b1Δy−b2Δx,k1=1Δx2+b22Δx+b22Δx,k2=1Δx2+b22Δx−b22Δx,k3=1Δy2+b12Δy+b12Δy,k4=1Δy2+b12Δy−b12Δy.


### 4.3.3. Discretization of Constitutive Function

The first derivative of constitutive equation in (21a), (21b), and (21c) applies first-order array of difference, and the discrete equations are as follows:
(41a)τxxi,j=2ui,j−ui−1,jΔx,(41b)τyyi,j=2vi,j−vi,j−1Δy,(41c)τxyi,j=τyxi,j=vi,j−vi−1,jΔx+ui,j−ui,j−1Δy.


## 4.3.1. Discretization of Stream Function

Stream function equation is the elliptic equations; it can be dispersed by arrays of central difference with second-order accuracy. The main diagonal of the discrete equation coefficient is dominant, which guarantee the convergence of algorithm, the discrete equation is:(37)ψi+1,j−2ψi,j+ψi−1,jΔx2+ψi,j+1−2ψi,j+ψi,j−1Δy2=−Ωi,j,(38)ψn+1i,j=121+βψni+1,j+ψni−1,j+βψni,j+1+βψni,j−1+Δx2Ωi,j,where Δx,Δy is the grid step length of x and y direction, β=Δx2/Δy2; n is the number of iterations.


## 4.3.2. Discretization of Vorticity Function

The left item of vorticity equation (38) applies the first kind of windward array of difference. That is, when coefficient of the first derivative of vorticity is greater than zero, then discrete the equation by array of backward difference; when the coefficient is less than or equal to zero, then discrete the equation by array of forward difference. Although difference equation has only first-order accuracy, this array of difference is equipped with main diagonally dominant unconditionally, thus the stability of numerical calculation can be ensured. The right item of vorticity equation (38) applies central difference methods. The discrete form is as follows:
(39)k0Ωi,j+k1Ωi+1,j+k2Ωi−1,j+k3Ωi,j+1+k4Ωi,j−1=14ΔxΔySi+1,j+1xx−Si−1,j+1xx−Si+1,j−1xx+Si−1,j−1xx−Si+1,j+1yy−Si−1,j+1yy−Si+1,j−1yy+Si−1,j−1yy+1Δy2Si,j+1xy−2Si,jxy+Si,j−1xy.The coefficients in the equation can be expressed as(40)b1=Reψi+1,j−ψi−1,j2Δx,b2=−Reψi,j+1−ψi,j−12Δy,k0=−2Δx2−2Δy2−b1Δy−b2Δx,k1=1Δx2+b22Δx+b22Δx,k2=1Δx2+b22Δx−b22Δx,k3=1Δy2+b12Δy+b12Δy,k4=1Δy2+b12Δy−b12Δy.


## 4.3.3. Discretization of Constitutive Function

The first derivative of constitutive equation in (21a), (21b), and (21c) applies first-order array of difference, and the discrete equations are as follows:
(41a)τxxi,j=2ui,j−ui−1,jΔx,(41b)τyyi,j=2vi,j−vi,j−1Δy,(41c)τxyi,j=τyxi,j=vi,j−vi−1,jΔx+ui,j−ui,j−1Δy.


## 4.4. Boundary Conditions

### 4.4.1. Velocity

In this paper, the flow velocity distribution at the nanochannel can be determined by (13) considering effect of solid and liquid. After dimensionless, it can be written as follows:
(42)u=−12φ′1n+1n+2R−yn+1n+1y+Rdpdx.


### 4.4.2. Velocity Derivative

According to the flow characteristic of Poiseuille, the constitutive equation and the flow velocity distribution (4), the first-order derivative in the channel can be written as:
(43)∂u∂y=1ϕ′∂p∂xR−yn2y.


### 4.4.3. Stress Value

(44)τxx=0,τxy=1+ϕyn1ϕ′∂p∂xR−yn2y,τyy=0.Substitute (44) into constitutive equation and partial stress equation, the boundary value of additional stress S are as follows:
(45)Sxx=Syy=0Sxy=∂p∂xR−yn2yn−1.


### 4.4.4.ψ and Ω

Based on the definition of vorticity and stream function, the distribution of stream function and vorticity at the nanochannel can be determined by velocity distribution as follows:(46)ψ=∂p∂x12ϕ′n+1n+2n+3n+1y+2RR−yn+2−2Rn+3,Ω=−1ϕ∂p∂xR−yn2y.


## 4.4.1. Velocity

In this paper, the flow velocity distribution at the nanochannel can be determined by (13) considering effect of solid and liquid. After dimensionless, it can be written as follows:
(42)u=−12φ′1n+1n+2R−yn+1n+1y+Rdpdx.


## 4.4.2. Velocity Derivative

According to the flow characteristic of Poiseuille, the constitutive equation and the flow velocity distribution (4), the first-order derivative in the channel can be written as:
(43)∂u∂y=1ϕ′∂p∂xR−yn2y.


## 4.4.3. Stress Value

(44)τxx=0,τxy=1+ϕyn1ϕ′∂p∂xR−yn2y,τyy=0.Substitute (44) into constitutive equation and partial stress equation, the boundary value of additional stress S are as follows:
(45)Sxx=Syy=0Sxy=∂p∂xR−yn2yn−1.


## 4.4.4.ψ and Ω

Based on the definition of vorticity and stream function, the distribution of stream function and vorticity at the nanochannel can be determined by velocity distribution as follows:(46)ψ=∂p∂x12ϕ′n+1n+2n+3n+1y+2RR−yn+2−2Rn+3,Ω=−1ϕ∂p∂xR−yn2y.


## 4.5. Sample Calculation and Results Comparison

Numerical calculation of the capillary flow has been proceeded in the capillary with radius 5μm, the solid-liquid effect index n is 0, 0.1, and 0.2, respectively, based on the established 2D numerical boundary layer equations. Radial velocity distribution curves in the capillary radius of 5 μm are showed in Figures 7 and 8.Figure 7
Velocity distribution at different solid-liquid index in capillary profiles (n=0,0.1).Figure 8
Velocity distribution at different solid-liquid index in capillary profiles (n=0.1,0.2).The interaction effect between solid and liquid phases in nanocapillary can be obviously seen by comparing two curves at different solid-liquid index in Figure7. When n is zero, the velocity distribution curve presents the convex type, the flow characteristics of the Newtonian fluid. When n is 0.1 or 0.2, the curve near the wall presents the concave type, the flow characteristics of the non-Newtonian fluid, which illustrates that the interaction effect is stronger close to wall surface. And the velocity is getting smaller and the increase amplitude is smaller from pore wall to the capillary center, that is to say, fluid viscosity coefficient has big changes close to the wall surface. Figures 7 and 8 are obtained by numerical calculation, the shape of curves and indicated results are in consistent with the approximate solutions obtained by the governing equation. Results proved that the numerical seepage model established in this paper can fulfill the precision and requirements of flow in confined nanopores.


## 5. Application in Two Typical Seepage Types: from Pore to Throat and Inversely from Throat to Pore

### 5.1. Pore to Throat

Parameters of pore to throat model are shown in Table1. Pressure gradient in the inlet is 0.0075 MPa/m; numerical equation is used to calculate the parameter field distribution of different sizes in pore to throat flow channel is obtained. Pore radius is 0.05 μm; the throat radius can be automatically calculated according to pore-throat ratio.Table 1
Four kinds of parameters in pore to throat model. Considering L/S interaction (μm).
Order number
Pore-throat ratio
Throat radius
Pore radius
Throat length
Pore length
ϕ′
n
1
5
0.01
0.05
20.0
50.0
3.0
0.1
2
10
0.005
0.05
20.0
50.0
3.0
0.1
3
20
0.0025
0.05
20.0
50.0
3.0
0.1
4
50
0.001
0.05
20.0
50.0
3.0
0.1


#### 5.1.1. The Characteristics of Velocity Distribution

Figure9 shows velocity distribution contours in pore to throat channel nonconsidering L/S interaction and considering L/S interaction. It shows that the flow velocity distribution of the fluid altered obviously when considering L/S interaction under the same inlet pressure gradient. First, the fluid velocity in the throat or inside the pore is affected by the throat and the pore wall, so the velocity dramatically reduces near the wall, and it is mainly influenced by the throat size. Secondly, when the pore-throat ratio is the same, the ratio of the velocity considering L/S interaction to the velocity without considering L/S interaction increases. When the pore-throat ratio is equal to five, the dimensionless velocity of the former near the throat of the lower pore is near 1, and the dimensionless velocity of the latter is 0.15, the ratio of the two is about 7, and the ratio is about 9 when the pore-throat ratio is 10. Thirdly, the velocity distribution of fluid changes obviously from the throat to the pore, that is, the fluid velocity considering the L/S interaction in the pores is very small. When the pore-throat ratio increases to a certain degree, almost no fluid can flow, as shown in Figure 9h, the dimensionless velocity in the pore reaches the order of magnitude of 10−3, the characteristic speed is converted into the actual flow rate, and its order of magnitude is 10−12 m/s (10−3 nm/s). In the real displacement process, the flow velocity cannot mobile or displace the oil.Figure 9
The velocity distribution of different pore-throat ratios in pore to throat channels nonconsidering L/S interaction (left) and considering L/S interaction (right).


#### 5.1.2. Streamline distribution

Figure10 shows flow diagrams of different pore-throat ratios, as pore-throat ratio increases, which means that when the throat radius reduces, the value of the stream function inside the pores successively decreases; when the pore throat ratio ranges from 5 to 50, the values of the stream function at the center of the inner pores are 300, 35, 3.5, and 0, respectively. Fluid influx through the pores is very small, the pore volume of corresponding fluid in the pores is very low, and there is almost no effective injection pore volume.Figure 10
Streamline graph with different pore-throat ratios in pore to throat channels considering L/S interaction.


### 5.2. Throat to Pore

#### 5.2.1. The Characteristics of Velocity Distribution

The velocity distribution contours in throat to pore channel nonconsidering L/S interaction and considering L/S interaction are shown in Figure11.Figure 11
The velocity distribution of different pore-throat ratios in throat to pore channels nonconsidering L/S interaction (left) and considering L/S interaction (right).


#### 5.2.2. Streamline Distribution

According to the velocity distribution, the distribution of velocity (Figure11) and streamlines (Figure 12) in throat to pore channels considering L/S interaction is the same with channels not considering L/S interaction, and the value is slightly lower, but the magnitude is not very clear. Because the effect of solid boundary on fluid is not obvious since the fluid inlet are pores with relatively larger radius.Figure 12
Streamline graph with different pore-throat ratio in contraction-expansion channels considering L/S interaction.


## 5.1. Pore to Throat

Parameters of pore to throat model are shown in Table1. Pressure gradient in the inlet is 0.0075 MPa/m; numerical equation is used to calculate the parameter field distribution of different sizes in pore to throat flow channel is obtained. Pore radius is 0.05 μm; the throat radius can be automatically calculated according to pore-throat ratio.Table 1
Four kinds of parameters in pore to throat model. Considering L/S interaction (μm).
Order number
Pore-throat ratio
Throat radius
Pore radius
Throat length
Pore length
ϕ′
n
1
5
0.01
0.05
20.0
50.0
3.0
0.1
2
10
0.005
0.05
20.0
50.0
3.0
0.1
3
20
0.0025
0.05
20.0
50.0
3.0
0.1
4
50
0.001
0.05
20.0
50.0
3.0
0.1


### 5.1.1. The Characteristics of Velocity Distribution

Figure9 shows velocity distribution contours in pore to throat channel nonconsidering L/S interaction and considering L/S interaction. It shows that the flow velocity distribution of the fluid altered obviously when considering L/S interaction under the same inlet pressure gradient. First, the fluid velocity in the throat or inside the pore is affected by the throat and the pore wall, so the velocity dramatically reduces near the wall, and it is mainly influenced by the throat size. Secondly, when the pore-throat ratio is the same, the ratio of the velocity considering L/S interaction to the velocity without considering L/S interaction increases. When the pore-throat ratio is equal to five, the dimensionless velocity of the former near the throat of the lower pore is near 1, and the dimensionless velocity of the latter is 0.15, the ratio of the two is about 7, and the ratio is about 9 when the pore-throat ratio is 10. Thirdly, the velocity distribution of fluid changes obviously from the throat to the pore, that is, the fluid velocity considering the L/S interaction in the pores is very small. When the pore-throat ratio increases to a certain degree, almost no fluid can flow, as shown in Figure 9h, the dimensionless velocity in the pore reaches the order of magnitude of 10−3, the characteristic speed is converted into the actual flow rate, and its order of magnitude is 10−12 m/s (10−3 nm/s). In the real displacement process, the flow velocity cannot mobile or displace the oil.Figure 9
The velocity distribution of different pore-throat ratios in pore to throat channels nonconsidering L/S interaction (left) and considering L/S interaction (right).


### 5.1.2. Streamline distribution

Figure10 shows flow diagrams of different pore-throat ratios, as pore-throat ratio increases, which means that when the throat radius reduces, the value of the stream function inside the pores successively decreases; when the pore throat ratio ranges from 5 to 50, the values of the stream function at the center of the inner pores are 300, 35, 3.5, and 0, respectively. Fluid influx through the pores is very small, the pore volume of corresponding fluid in the pores is very low, and there is almost no effective injection pore volume.Figure 10
Streamline graph with different pore-throat ratios in pore to throat channels considering L/S interaction.


## 5.1.1. The Characteristics of Velocity Distribution

Figure9 shows velocity distribution contours in pore to throat channel nonconsidering L/S interaction and considering L/S interaction. It shows that the flow velocity distribution of the fluid altered obviously when considering L/S interaction under the same inlet pressure gradient. First, the fluid velocity in the throat or inside the pore is affected by the throat and the pore wall, so the velocity dramatically reduces near the wall, and it is mainly influenced by the throat size. Secondly, when the pore-throat ratio is the same, the ratio of the velocity considering L/S interaction to the velocity without considering L/S interaction increases. When the pore-throat ratio is equal to five, the dimensionless velocity of the former near the throat of the lower pore is near 1, and the dimensionless velocity of the latter is 0.15, the ratio of the two is about 7, and the ratio is about 9 when the pore-throat ratio is 10. Thirdly, the velocity distribution of fluid changes obviously from the throat to the pore, that is, the fluid velocity considering the L/S interaction in the pores is very small. When the pore-throat ratio increases to a certain degree, almost no fluid can flow, as shown in Figure 9h, the dimensionless velocity in the pore reaches the order of magnitude of 10−3, the characteristic speed is converted into the actual flow rate, and its order of magnitude is 10−12 m/s (10−3 nm/s). In the real displacement process, the flow velocity cannot mobile or displace the oil.Figure 9
The velocity distribution of different pore-throat ratios in pore to throat channels nonconsidering L/S interaction (left) and considering L/S interaction (right).


## 5.1.2. Streamline distribution

Figure10 shows flow diagrams of different pore-throat ratios, as pore-throat ratio increases, which means that when the throat radius reduces, the value of the stream function inside the pores successively decreases; when the pore throat ratio ranges from 5 to 50, the values of the stream function at the center of the inner pores are 300, 35, 3.5, and 0, respectively. Fluid influx through the pores is very small, the pore volume of corresponding fluid in the pores is very low, and there is almost no effective injection pore volume.Figure 10
Streamline graph with different pore-throat ratios in pore to throat channels considering L/S interaction.


## 5.2. Throat to Pore

### 5.2.1. The Characteristics of Velocity Distribution

The velocity distribution contours in throat to pore channel nonconsidering L/S interaction and considering L/S interaction are shown in Figure11.Figure 11
The velocity distribution of different pore-throat ratios in throat to pore channels nonconsidering L/S interaction (left) and considering L/S interaction (right).


### 5.2.2. Streamline Distribution

According to the velocity distribution, the distribution of velocity (Figure11) and streamlines (Figure 12) in throat to pore channels considering L/S interaction is the same with channels not considering L/S interaction, and the value is slightly lower, but the magnitude is not very clear. Because the effect of solid boundary on fluid is not obvious since the fluid inlet are pores with relatively larger radius.Figure 12
Streamline graph with different pore-throat ratio in contraction-expansion channels considering L/S interaction.


## 5.2.1. The Characteristics of Velocity Distribution

The velocity distribution contours in throat to pore channel nonconsidering L/S interaction and considering L/S interaction are shown in Figure11.Figure 11
The velocity distribution of different pore-throat ratios in throat to pore channels nonconsidering L/S interaction (left) and considering L/S interaction (right).


## 5.2.2. Streamline Distribution

According to the velocity distribution, the distribution of velocity (Figure11) and streamlines (Figure 12) in throat to pore channels considering L/S interaction is the same with channels not considering L/S interaction, and the value is slightly lower, but the magnitude is not very clear. Because the effect of solid boundary on fluid is not obvious since the fluid inlet are pores with relatively larger radius.Figure 12
Streamline graph with different pore-throat ratio in contraction-expansion channels considering L/S interaction.


## 6. Conclusions

In summary, we have established a two-dimensional numerical model of seepage flow considering the liquid and solid interaction in confined nanopores. This complex model has been discretized by the introduction of vortex equation. And it has been applied in two types of nanopore flow: from pore to throat and inversely from throat to pore. The results and recognitions can be concluded as follows:(1)
The seepage behavior in confined nanopores is different from bulk space, since the liquid and solid interaction would increase the viscosity of hydrocarbons near the pore wall.
(2)
Vortex equation can efficiently reduce the number of unknowns in solving complex numerical equation.
(3)
The liquid and solid interaction mainly affected by capillary radius. It completely cannot be neglected in the seepage model if the capillary radius is smaller than 50 nm when liquid and solid interaction index is bigger than 0.1.
(4)
When the hydrocarbons flow from nanoscale pore to throat, the liquid and solid interaction cannot be neglected. Nevertheless, the effect of liquid and solid interaction is not critical when it flows from throat to pore.In addition, the complexity of wettability has not been included in this model, since there are many situations, such as water wet, oil wet, and mixed wet. This work has taken the most common situation “water wet” in tight or shale reservoir into consideration. The other situations of wettability combined with liquid and solid interaction effect will be further investigated in the next work.



---

*Source: 8302782-2018-08-19.xml*