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

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


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][3][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][6][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][10][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][21][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 boundarylayer 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][35][36][37][38][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.

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 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.Velocity component w * is equal to zero in twodimensional 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:

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 Of the viscous force and inertial force item in (4a), truncating the order items with small quantity, respectively, it can be obtained: 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): 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) 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 Figure 1, 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: 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 Integrating (11), combined with the boundary conditions, the formula of fluid velocity distribution within the capillary can be calculated as follows: When y * is equal to zero, the maximum speed in capillary is as follows: The flow rate in capillary is: where ϕ′ = μ 0 and n = 0, and ( 15) is just the Poiseuille formula.
Figures 2 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-   4 Geofluids 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 5 shows the relationship between the average velocity ratio of fluid V n (which considers solid-liquid interaction) and V 0 (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 solidliquid 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.

Numerical Modeling and Solutions
Under Cartesian coordinate system, the component forms of τ ij and ε ij are as follows: The component form of incompressible Newtonian fluid constitutive equation (16) under the Cartesian coordinate system can be written as 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:   .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: 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): The introduction of stream function and vorticity: There is one curve AB between the two streamlines as shown in Figure 6, and then the volume flow of unit thickness is The stream function equation: 6 Geofluids Momentum equation can be written into vorticity equation in the following form: After finishing, oval vorticity equation can be obtained as follows: 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.Therefore, (32) can be written as Vorticity equation ( 32) can be written as where Δx, Δy is the grid step length of x and y direction, β = Δx 2 / Δy 2 ; n is the number of iterations.7 Geofluids applies central difference methods.The discrete form is as follows: The coefficients in the equation can be expressed as 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: 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: Table 1: Four kinds of parameters in pore to throat model.Considering L/S interaction (μm).

Geofluids
The interaction effect between solid and liquid phases in nanocapillary can be obviously seen by comparing two curves at different solid-liquid index in Figure 7.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.

Application in Two Typical Seepage
Types: from Pore to Throat and Inversely from Throat to Pore the throat radius can be automatically calculated according to pore-throat ratio.
5.1.1.The Characteristics of Velocity Distribution.Figure 9 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 porethroat 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 porethroat 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.
5.1.2.Streamline distribution.Figure 10 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  10 Geofluids small, the pore volume of corresponding fluid in the pores is very low, and there is almost no effective injection pore volume.

Streamline Distribution.
According to the velocity distribution, the distribution of velocity (Figure 11) 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.

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 11 Geofluids 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.The first and second fluid viscosity, Pa.s δ ij : Ronnie Kerr mark.

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 2 Geofluids governing equation of boundary layer fluid in the nanochannel has established.

Figure 1 :
Figure 1: Flow in circular section of capillary.

4. 1 .
Numerical Equation 4.1.1.Constitutive Equation.Constitutive equation of Newtonian fluid in the form of a tensor is as follows:

36 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:

1 𝜓 2 𝜐Figure 6 :
Figure 6: Sketch map of the planar flow field and flow line.

R − y n 2 y, τ yy = 0 44 Substitute ( 44 )
into constitutive equation and partial stress equation, the boundary value of additional stress S are as follows: S xx = S yy = 0 S xy = ∂p ∂x R − y n 2y n−1 45 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 4. 5 .
Sample Calculation and ResultsComparison.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 Figures7 and 8.

Pore-throat ratio is 5 Pore-throat ratio is 10 Pore-throat ratio is 20 PoreFigure 9 :
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).
ratio is 10Pore-throat ratio is 20Pore-throat ratio is 50

Figure 10 :
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 Figure 11.

Pore-throat ratio is 5 Pore-throat ratio is 10 Pore-throat ratio is 20 PoreFigure 11 :
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).

Nomenclature μ 0 :
Bulk fluid viscosity, Pa.s ϕ′/y n : Additional viscosity imposed by solid surface, Pa.s n: Liquid/solid interaction index y: Distance to a solid surface, meter ϕ ′ : Interaction coefficient of liquid/solid molecule μ 1 ′ : Dipole moment of water molecule μ 2 ′ : Dipole moment of surface molecule α 1 : Polarizability energy of water molecule α 2 : Polarizability energy of surface molecule I 1 : Ionization energy of water molecule I 2 : Ionization energy of surface molecule μ: Fluid viscosity coefficient of boundary layer, Pa.s μ * : Velocity component of fluid particle, m/s v * : Velocity component of fluid particle, m/s p * : Positive pressure of fluid particle, Pa R * : Capillary radius, m R : Dimensionless radius of channel τ * ij : Stress tensor, Pa ε * ij : Strain rate tensor, m/t λ: 5.1.Pore to Throat.Parameters of pore to throat model are shown in Table 1.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;