Fractal Analytical Solutions for Nonlinear Two-Phase Flow in Discontinuous Shale Gas Reservoir

Abstract

The paths of a two-phase flow are usually non-linear and discontinuous in the production of shale gas development. To research the influence mechanism between shale gas and water, several integer two-phase flow models have been studied but few analytical solutions have been obtained on shale gas and water pressure. This study first developed a local fractional mathematical model for gas and water two-phase flow in shale gas production. The model thus created considers the effects of capillary pressure, the fractal dimension of the flow pipe, and the discontinuity of the flow path. Second, the local fractional traveling wave method and variational iteration method were applied to this model for the development of iterative analytical solutions. Both shale gas and water pressure were analytically derived. Third, the depressurization process of the shale gas and water was analyzed, and a parametric study was conducted to explore the impacts of fractional dimension, entry capillary pressure, and travel wave velocity on shale gas pressure. Finally, our conclusions are drawn, based on the results of these studies.

1. Introduction

The formation of a complex fracture system through fracturing modification is key to the efficient development of shale gas. To create volumetric fracture networks, a number of emerging technologies are currently being used and studied in the field of shale reservoir fracturing, such as the in situ combustion-explosion fracturing of methane, heat treatment stimulation on the basis of hydraulic fracturing, and so on. As shown in Figure 1, a complex fracture network will form after fracturing, and the permeability of the shale reservoir is greatly enhanced. Meanwhile, after diffusing into the pore surface, the shale gas will be desorbed into the intergranular macropores and fractures and will then flow into the well-bore [1]. During the process of flowing into the well-bore, the shale gas in fractures must overcome the capillary pressure at the gas-water interface in the pores and fractures for normal production. Therefore, the two-phase flow cannot be ignored throughout the whole process of shale gas production.

Two-phase flows have been extensively utilized in rock and power engineering, in scenarios such as geothermal extraction, tight sandstone gas, and so on. As shown in Figure 2, the fractures of a shale reservoir are filled with two-phase flow [1]. Both shale gas and water flow to the well-bore in the two-phase flow stage. The mechanism and transport behaviors of gas-water two-phase flows are worth careful study for the prediction of the shale gas field production rate.

In order to closely simulate and predict the field production of unconventional gas reservoirs, scholars have conducted many studies on the process of two-phase flow and have proposed many two-phase flow mathematical models. For instance, Clarkson and McGovern [3] developed a model to simulate the two-phase flow of gas and water below desorption pressure. Wang et al. developed a numerical two-phase model with multiscale diffusion mechanisms at the flow-back stage after the hydraulic fracturing of a shale gas reservoir [4]. However, as the production of shale gas is actually a depressurization process, an analytical solution for gas pressure is necessary to forecast the field production of shale gas faster and more conveniently. Wang et al. presented a semi-analytical method for accurately modeling two-phase flow behavior and quickly predicting production performance [5]. An analytical solution for the oil production rate has been obtained by linearizing the mathematical model. To estimate the residual saturation of displaced fluid, M. Adibifard proposed a two-phase flow model in a capillary tube and obtained an analytical solution based on the decomposition of total pressure drop into Poiseuille and Young–Laplace terms [6]. The above analytical solutions of two-phase flow models ignore the nonlinear terms of governing equations. However, the nonlinear term is an important factor to achieve analytical solutions for two-phase flows.

In recent shale gas production models, more and more fractal theories are being connected to establish permeability models in porous media from different points of view. In 2002, Boming Yu et al. established a fractal permeability model of dual porous media, according to the fractal characteristics of pore structure, and calculated the fractal dimension of pore area and tortuousness by using the box-counting method, respectively [7]. The model was in good agreement with the experimental results. In 2010, Shou et al. established the permeability models of fibrous porous media under different fractal conditions and discussed the coupling relationship among fractal dimensions, porosity, and permeability [8]. In 2014, Yang et al. conducted a fractal study on the shale in the Sichuan Basin, China, and concluded that the pore surface of shale is not uniform and meets the fractal characteristics, with a fluctuation range of fractal dimensions between 2.86 and 2.83 [9]. Bu et al. used the surface fractal dimension, D₁, and the pore structure fractal dimension, D₂, to identify all the fractal characteristics, and discussed the relationship among the main controlling factors of the two fractal dimensions—the pore size and the adsorption capacity [10]. In 2019, Hu et al. studied the impacts of zone fractal properties on the shale gas productivity of a multiple fractured horizontal well through a fractal permeability model related to aperture distribution and fracture length [11]. Unlike previous studies, in this paper, the diameter and tortuosity fractal dimensions are introduced to the permeability model to describe the complex flow paths after the fracturing of shale reservoirs.

As two-phase flow paths in fractures are usually non-linear and discontinuous [12], the existing integer two-phase flow models [13] cannot accurately describe the heterogeneous pores and tortuous fractures of the real shale structure. Therefore, an anomalous two-phase flow may dominate the process of the two-phase flow [14,15]. The effects of these anomalous flows and the properties of a two-phase flow should be studied. A fractional derivative model is a good choice for the simulation of a two-phase flow with a tortuous path. This kind of model is intended to replace the derivatives of integer order with fractional order in the governing equations, on the basis of classical differential equations [14]. The theory of local fractional derivatives has been successfully applied to many problems in fluid mechanics [16,17,18].

This study considers the influences of capillary pressure, the fractional dimension of cataclastic rock, and seepage discontinuity on two-phase flow. A fractional coupling mathematical model for two-phase flow pressures was first developed, based on the two-phase flow interactions of multi-physical processes during shale gas production. Secondly, without ignoring the nonlinear phase, the water and gas pressures are analytically solved by the local fractional traveling wave method and Laplace variational iteration method, respectively. The fractional analytical solutions of two-phase flow pressures through time and space are analytically expressed. Then, the impacts of fractional dimensions, entry pressure, and travel wave velocity on gas pressure are investigated. Finally, conclusions are drawn, based on the results of these studies.

2. Governing Equations for Shale Gas and Water Two-Phase Flows

The shale gas-water two-phase flow occurs in the fractures of a shale gas reservoir. The governing equations of shale gas-water two-phase flows are proposed, as follows.

2.1. Capillary Pressure

$p_c$ stands for the difference between shale gas and water pressure [19]:

$$p_c=p_g-p_w$$

(1)

where $p_c$ is the capillary pressure, $p_g$ is the pressure of shale gas, and $p_w$ is the pressure of water.

The saturation of shale gas and water meets:

$$s_w+s_g=1$$

(2)

with a relationship of:

$$s_w=s_w^{*}(1-s_{rg}-s_{rw})+s_{rw}$$

(3)

where $s_w$ and $s_g$ are the saturations for water and gas phases, respectively. $s_w^{*}$ is the normalized saturation of water. $s_{rw}$ and $s_{rg}$ are the irreducible saturations of water and shale gas, respectively. The subscript g refers to the shale gas phase, and w is used to denote the water phase.

The normalized saturation of water is a function of capillary pressure [20]:

$$s_w^{*}={\left(\frac{p_e}{p_c}\right)}^{\lambda }$$

(4)

where $p_e$ is the entry capillary pressure and $\lambda$ is the pore size distribution index.

Integrating Equation (4) into Equation (3) yields

$$s_w=(1-s_{rg}-s_{rw}){\left(\frac{p_e}{p_c}\right)}^{\lambda }+s_{rw}.$$

(5)

Therefore, the fractional derivative of water saturation with respect to time is:

$$\frac{{\partial }^{\alpha }{s}_w}{\partial t^{\alpha }}=\left(1-s_{rg}-s_{rw}\right)p_e^{\lambda }·\frac{{\partial }^{\alpha }\left(p_c^{-\lambda }\right)}{\partial t^{\alpha }}.$$

(6)

Similarly, for the shale gas phase:

$$s_g=1-s_{rw}-(1-s_{rg}-s_{rw}){\left(\frac{p_e}{p_c}\right)}^{\lambda },$$

(7)

$$\frac{{\partial }^{\alpha }{s}_g}{\partial t^{\alpha }}=-\left(1-s_{rg}-s_{rw}\right)p_e^{\lambda }\frac{{\partial }^{\alpha }\left(p_c^{-\lambda }\right)}{\partial t^{\alpha }}.$$

(8)

2.2. Governing Equations of the Two-Phase Flow

For a two-phase flow, the fractional continuity equation can be written out for each phase. According to the law of mass conservation, the fractional continuity equation for gas is:

$$\frac{{\partial }^{\alpha }\left(\varphi s_g{\rho }_g(x,t)\right)}{\partial t^{\alpha }}+{\nabla }^{\alpha }({\rho }_g(x,t)·v_g(x,t))={\rho }_g{q}_{\mathrm{g}}.$$

(9)

The fractional continuity equation for water is:

$$\frac{{\partial }^{\alpha }\left(\varphi s_w{\rho }_w\right)}{\partial t^{\alpha }}+{\nabla }^{\alpha }\left({\rho }_w{v}_w(x,t)\right)={\rho }_w{q}_w$$

(10)

where $\varphi$ is the porosity of the shale reservoir. ${\rho }_{\mathrm{g}}$, $v_g$, $q_{\mathrm{g}}$ are the density, velocity, and source of gas, respectively. ${\rho }_w$, $v_w$, $q_w$ are the density, velocity, and source of water, respectively.

The water density ${\rho }_w$ is regarded as a constant in this study, and the gas density ${\rho }_g$ follows the equation of state:

$${\rho }_g=\frac{M}{ZRT}{p}_g$$

(11)

where M and Z are the molecular weight and compressibility factor, respectively. R denotes the universal gas constant, and T is the temperature of the shale gas.

Integrating Equation (11) into Equation (9) yields

$$\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}(\varphi s_{\mathrm{g}}{p}_{\mathrm{g}})+{\nabla }^{\alpha }·(p_{\mathrm{g}}{\mathit{V}}_g)=p_g{q}_{\mathrm{g}}.$$

(12)

Hence, Equation (10) is simplified into

$$\frac{{\partial }^{\alpha }}{\partial t^{\alpha }}(\varphi s_w)+{\nabla }^{\alpha }·{\mathit{V}}_w=q_w.$$

(13)

Both ${\mathit{V}}_g$ and ${\mathit{V}}_w$ in porous media can be described by the fractional Darcy law, as [21]:

$${\mathit{V}}_g=-\frac{k{k}_{rg}}{{\mu }_g}{\nabla }^{\alpha }{p}_g$$

(14)

$${\mathit{V}}_w=-\frac{k{k}_{rw}}{{\mu }_w}{\nabla }^{\alpha }{p}_w$$

(15)

where $k$ is the permeability of the shale reservoir, and $k_{rg}$ and $k_{rw}$ are the relative permeabilities of shale gas and water, respectively. The viscosity of the in situ conditions is represented by ${\mu }_g$ for shale gas and ${\mu }_w$ for water. $p_g$ and $p_w$ are the gas and water pressure, respectively.

The permeability of shale is defined as [22]:

$$k=\frac{D_T+D_f}{3+D_T-D_f}{k}_0$$

(16)

where $k_0$ is the initial permeability of the shale reservoir. $D_T$ and $D_f$ are the tortuosity fractal dimension and the fractal dimension of diameter size, respectively.

If $q_w$ and $q_{\mathrm{g}}$ are all 0, the governing equations can be derived from Equations (12)–(15).

For the shale gas phase:

$$\frac{\varphi {\partial }^{\alpha }(s_g{p}_g)}{\partial t^{\alpha }}-{\nabla }^{\alpha }\left(\frac{k{k}_{rg}{p}_g}{{\mu }_g}{\nabla }^{\alpha }{p}_g\right)=0.$$

(17)

For the water phase:

$$\frac{\varphi {\partial }^{\alpha }{s}_w}{\partial t^{\alpha }}-\frac{k}{{\mu }_w}{\nabla }^{\alpha }\left(k_{rw}{\nabla }^{\alpha }{p}_w\right)=0.$$

(18)

The relative permeabilities of $k_{rw}$ and $k_{rg}$ follow the Brooks and Corey functions as [23,24]:

$$k_{rw}=k_{rw}^{\max }{\left(s_w^{*}\right)}^{N_w}$$

(19)

$$k_{rg}=k_{rg}^{\max }{\left(s_g^{*}\right)}^{N_g}$$

(20)

where $s_w^{*}+s_g^{*}=1$ and $s_g^{*}=1-{(\frac{p_e}{p_c})}^{\lambda }$.

Integrating Equation (4) into Equation (19) yields

$$k_{rw}=k_{rw}^{\max }{p}_e^{\lambda N_w}·p_c^{-\lambda N_w}$$

(21)

and:

$$k_{rg}=k_{rg}^{\max }{\left(1-{(\frac{p_e}{p_c})}^{\lambda }\right)}^{N_g}.$$

(22)

Integrating Equations (1), (6) and (21) into Equation (18) get the governing equation for water flows:

$$\varphi \left(1-s_{rg}-s_{rw}\right)p_e^{\lambda }\frac{{\partial }^{\alpha }\left[{\left(p_g-p_w\right)}^{-\lambda }\right]}{\partial t^{\alpha }}-\frac{k{k}_{rw}^{\max }{p}_e^{\lambda N_w}}{{\mu }_w}{\nabla }^{\alpha }\left[{\left(p_g-p_w\right)}^{-\lambda N_w}{\nabla }^{\alpha }{p}_w\right]=0$$

(23)

Similarly, the governing equation for gas flow can be obtained from Equations (1), (8), (22) and (17):

$$\begin{array}{l}\varphi \left[1-s_{rw}-(1-s_{rg}-s_{rw})p_e^{\lambda }·{\left(p_g-p_w\right)}^{-\lambda }\right]\frac{{\partial }^{\alpha }{p}_g}{\partial t^{\alpha }}-\varphi \left(1-s_{rg}-s_{rw}\right)p_e^{\lambda }·p_g\frac{{\partial }^{\alpha }\left({\left(p_g-p_w\right)}^{-\lambda }\right)}{\partial t^{\alpha }}\\ -\frac{k{k}_{rg}^{\max }}{{\mu }_g}{\nabla }^{\alpha }\left({\left(1-{(\frac{p_e}{p_g-p_w})}^{\lambda }\right)}^{N_g}·p_g·{\nabla }^{\alpha }{p}_g\right)=0\end{array}.$$

(24)

3. Analytical Solutions

3.1. Conversion to the Traveling Wave Solution

For a non-linear equation, the traveling wave solution is:

$$\phi (x,t,u,u_x,u_t,u_{xx},u_{xt},u_{tt},\dots )=0$$

(25)

where x and t represent the space and time coordinates, respectively. u is the function of x and t, and $\phi$ is a suitable function for u and its derivative.

If the solution of function $\phi (\xi )$ of Equation (25) merely depends on x and t in the form of $\xi =x+ct$, and if c is a constant to represent the wave velocity, $\phi (\xi )$ will be the traveling wave solution.

Secondly, the solutions of Equations Equations (23) and (24) are converted into the following definite problems. If:

$$\xi =x^{\alpha }+y^{\alpha }+z^{\alpha }+c{t}^{\alpha }$$

(26)

then:

$$p_g(\xi )=p_g(x^{\alpha }+y^{\alpha }+z^{\alpha }+c{t}^{\alpha })=p_g(x,y,z,t)$$

(27)

$$p_w(\xi )=p_w(x^{\alpha }+y^{\alpha }+z^{\alpha }+c{t}^{\alpha })=p_w(x,y,z,t)$$

(28)

$$\frac{{\partial }^{\alpha }{p}_g}{\partial t^{\alpha }}=c\frac{{\partial }^{\alpha }{p}_g}{\partial {\xi }^{\alpha }}, \frac{{\partial }^{\alpha }{p}_w}{\partial t^{\alpha }}=c\frac{{\partial }^{\alpha }{p}_w}{\partial {\xi }^{\alpha }}$$

(29)

$${\nabla }^{\alpha }=\frac{{\partial }^{\alpha }}{\partial x^{\alpha }}+\frac{{\partial }^{\alpha }}{\partial y^{\alpha }}+\frac{{\partial }^{\alpha }}{\partial z^{\alpha }}=\frac{{\partial }^{\alpha }}{\partial {\xi }^{\alpha }}, {\nabla }^{2\alpha }=\frac{{\partial }^{2\alpha }}{\partial x^{2\alpha }}+\frac{{\partial }^{2\alpha }}{\partial y^{2\alpha }}+\frac{{\partial }^{2\alpha }}{\partial z^{2\alpha }}=\frac{{\partial }^{2\alpha }}{\partial {\xi }^{2\alpha }}.$$

(30)

Equation (23) is translated into:

$$c\varphi \left(1-s_{rg}-s_{rw}\right)p_e^{\lambda }\frac{{\partial }^{\alpha }\left[{\left(p_g-p_w\right)}^{-\lambda }\right]}{\partial {\xi }^{\alpha }}-\frac{k{k}_{rw}^{\max }{p}_e^{\lambda N_w}}{{\mu }_w}\frac{{\partial }^{\alpha }}{\partial {\xi }^{\alpha }}\left[{\left(p_g-p_w\right)}^{-\lambda N_w}·\frac{{\partial }^{\alpha }{p}_w}{\partial {\xi }^{\alpha }}\right]=0.$$

(31)

Equation (31) can be further simplified into:

$$p_g=p_w+{\left(A\frac{{\partial }^{\alpha }{p}_w}{\partial {\xi }^{\alpha }}\right)}^{\frac{1}{\lambda \left(N_w-1\right)}}$$

(32)

where:

$$A=\frac{k{k}_{rw}^{\max }{p}_e^{\lambda \left(N_w-1\right)}}{c{\mu }_w\varphi \left(1-s_{rg}-s_{rw}\right)}$$

(33)

Equation (24) is translated into:

$$\begin{array}{l}\varphi \left[1-s_{rw}-(1-s_{rg}-s_{rw})p_e^{\lambda }·{\left(p_g-p_w\right)}^{-\lambda }\right]\frac{{\partial }^{\alpha }{p}_g}{\partial {\xi }^{\alpha }}-\varphi \left(1-s_{rg}-s_{rw}\right)p_e^{\lambda }·p_g\frac{{\partial }^{\alpha }\left({\left(p_g-p_w\right)}^{-\lambda }\right)}{\partial {\xi }^{\alpha }}\\ -\frac{k{k}_{rg}^{\max }}{{\mu }_g}·\frac{{\partial }^{\alpha }}{\partial {\xi }^{\alpha }}\left({\left(1-{(\frac{p_e}{p_g-p_w})}^{\lambda }\right)}^{N_g}·p_g·\frac{{\partial }^{\alpha }{p}_g}{\partial {\xi }^{\alpha }}\right)=0\end{array}.$$

(34)

For the convenience of calculation, $u(\xi )$ is used as a substitute for $p_w(\xi )$, as:

$$p_w(\xi )=u(\xi ).$$

(35)

In this study, $N_w=2$, $N_g=1$ and $\lambda =1$; integrating Equations (32) and (35) into Equation (34) yields:

$$\frac{{\partial }^{\alpha }u\left(\xi \right)}{\partial {\xi }^{\alpha }}+A\frac{{\partial }^{2\alpha }u\left(\xi \right)}{\partial {\xi }^{2\alpha }}+Bu\left(\xi \right)+D+f\left(u\left(\xi \right),\frac{{\partial }^{\alpha }u\left(\xi \right)}{\partial {\xi }^{\alpha }},\frac{{\partial }^{2\alpha }u\left(\xi \right)}{\partial {\xi }^{2\alpha }}\right)=0$$

(36)

where:

$$B=-\frac{k{k}_{rg}^{\max }{p}_e}{A\left[{\mu }_g\varphi \left(1-s_{rw}\right)-k{k}_{rg}^{\max }{p}_e\right]}$$

(37)

$$D=-\frac{{\mu }_g\varphi \left(1-s_{rg}-s_{rw}\right)p_e}{A\left[{\mu }_g\varphi \left(1-s_{rw}\right)-k{k}_{rg}^{\max }{p}_e\right]}$$

(38)

$$\begin{array}{l}f\left(u\left(\xi \right),\frac{{\partial }^{\alpha }u\left(\xi \right)}{\partial {\xi }^{\alpha }},\frac{{\partial }^{2\alpha }u\left(\xi \right)}{\partial {\xi }^{2\alpha }}\right)=\frac{{\mu }_g\left(1-s_{rg}-s_{rw}\right)p_e}{A\left(1-s_{rw}\right)\left[{\mu }_g\varphi \left(1-s_{rg}\right)-k{k}_{rg}^{\max }{p}_e\right]}u\left(\xi \right)\frac{{\partial }^{2\alpha }u\left(\xi \right)}{\partial {\xi }^{2\alpha }}\\ /{\left(\frac{{\partial }^{\alpha }u\left(\xi \right)}{\partial {\xi }^{\alpha }}\right)}^2+\frac{k{k}_{rg}^{\max }}{{\mu }_g\varphi \left(1-s_{rg}\right)-k{k}_{rg}^{\max }{p}_e}[u\left(\xi \right)\frac{{\partial }^{\alpha }u\left(\xi \right)}{\partial {\xi }^{\alpha }}+A^2{\left(\frac{{\partial }^{\alpha }u\left(\xi \right)}{\partial {\xi }^{\alpha }}\right)}^2-Au\left(\xi \right)\frac{{\partial }^{2\alpha }u\left(\xi \right)}{\partial {\xi }^{2\alpha }}\\ +A^2\frac{{\partial }^{\alpha }u\left(\xi \right)}{\partial {\xi }^{\alpha }}\frac{{\partial }^{2\alpha }u\left(\xi \right)}{\partial {\xi }^{2\alpha }}-p_{e}u\left(\xi \right)\frac{{\partial }^{2\alpha }u\left(\xi \right)}{\partial {\xi }^{2\alpha }}/\frac{{\partial }^{\alpha }u\left(\xi \right)}{\partial {\xi }^{\alpha }}]\end{array}.$$

(39)

3.2. Local Variational Iteration Method

The local variational iteration method is used to solve the nonlinear problem (36). For a general partial differential equation, below:

$$L_{\alpha }u+N_{\alpha }u=g(\xi )$$

(40)

where $L_{\alpha }$ is the linear operator, and $N_{\alpha }$ is a nonlinear operator of u.

Therefore, the correction functional for Equation (40) is:

$$u_{n+1}(\xi )=u_n(\xi )+\frac{1}{\Gamma (1+\alpha )}{\displaystyle \underset{0}{\overset{\xi }{\int }}\gamma }\left\{L_{\alpha }{u}_n(\omega )+N_{\alpha }{\tilde{u}}_n(\omega )\right\}{\left(d\omega \right)}^{\alpha }$$

(41)

where $\gamma$ is a Lagrange multiplier, identified by the fractional variational theory. ${\tilde{u}}_n$ is a restricted variational under $\delta {\tilde{u}}_n=0$ [16].

The fractional multiplier of Equation (41) can be determined as follows [16]:

$$\gamma ={\left(-1\right)}^n\frac{{\left(\omega -\xi \right)}^{\left(n-1\right)\alpha }}{\Gamma \left(1+\left(n-1\right)\alpha \right)}$$

(42)

Integrating Equation (42) into Equation (41) yields:

$$u_{n+1}(\xi )=u_n(\xi )+{}_{0}I{{}_{\omega }}^{\left(\alpha \right)}\left\{{\left(-1\right)}^n\frac{{\left(\omega -\xi \right)}^{\left(n-1\right)\alpha }}{\Gamma \left(1+\left(n-1\right)\alpha \right)}\left[L_{\alpha }{u}_n(\omega )+N_{\alpha }{\tilde{u}}_n(\omega )\right]\right\}$$

(43)

where:

$${}_{0}I{{}_{\omega }}^{\left(\alpha \right)}=\frac{1}{\Gamma (1+\alpha )}{\displaystyle \underset{0}{\overset{\xi }{\int }}\{·\}}{\left(d\omega \right)}^{\alpha }.$$

(44)

In this study, the linear operator in Equation (36) is:

$$L_{\alpha }u(\omega )=\frac{{\partial }^{\alpha }u\left(\xi \right)}{\partial {\xi }^{\alpha }}+A\frac{{\partial }^{2\alpha }u\left(\xi \right)}{\partial {\xi }^{2\alpha }}+Bu\left(\xi \right)+D.$$

(45)

The nonlinear operator $N_{\alpha }{\tilde{u}}_n(\omega )$ can be derived using Equation (39) as:

$$N_{\alpha }{\tilde{u}}_n(\omega )=f\left(\tilde{u}\left(\xi \right),\frac{{\partial }^{\alpha }\tilde{u}\left(\xi \right)}{\partial {\xi }^{\alpha }},\frac{{\partial }^{2\alpha }\tilde{u}\left(\xi \right)}{\partial {\xi }^{2\alpha }}\right).$$

(46)

The local fractional partial differential equation can be rewritten as:

$$\frac{{\partial }^{\alpha }u\left(\xi \right)}{\partial {\xi }^{\alpha }}+A\frac{{\partial }^{2\alpha }u\left(\xi \right)}{\partial {\xi }^{2\alpha }}+Bu\left(\xi \right)+D+N_{\alpha }{\tilde{u}}_n(\omega )=0.$$

(47)

In the one-dimensional case, integrating Equations (29), (30), and (35) into Equation (47) yields:

$$\frac{1}{c}\frac{{\partial }^{\alpha }{u}_n\left(x,t\right)}{\partial t^{\alpha }}+\frac{A}{c}\frac{{\partial }^{2\alpha }{u}_n\left(x,t\right)}{\partial t^{2\alpha }}+Bu\left(x,t\right)+D+N_{\alpha }{\tilde{u}}_n\left(x,t\right)=0.$$

(48)

The local fractional Lagrange multiplier can be derived by Equation (42) as:

$${\gamma }_1=\frac{{\left(s-x\right)}^{\alpha }}{\Gamma \left(1+\alpha \right)}.$$

(49)

The corresponding local fractional iteration algorithms of Equation (43) become:

$$\begin{array}{l}{u}_{n+1}(x,t)=u_n(x,t)+{}_{0}I{{}_x}^{\left(\alpha \right)}\left\{\frac{{\left(s-x\right)}^{\alpha }}{\Gamma \left(1+\alpha \right)}\left[\left(\frac{{\partial }^{\alpha }{u}_n\left(x,t\right)}{\partial x^{\alpha }}+\frac{A}{c}\frac{{\partial }^{2\alpha }{u}_n\left(x,t\right)}{\partial t^{2\alpha }}+N_{\alpha }{\tilde{u}}_n\left(x,t\right)\right)\right]\right\}\\ =u_n(x,t)+\frac{1}{\Gamma (1+\alpha )}{\displaystyle \underset{0}{\overset{x}{\int }}\left\{\frac{{\left(s-x\right)}^{\alpha }}{\Gamma \left(1+\alpha \right)}\left[\left(\frac{{\partial }^{\alpha }{u}_n\left(x,t\right)}{\partial x^{\alpha }}+\frac{A}{c}\frac{{\partial }^{2\alpha }{u}_n\left(x,t\right)}{\partial t^{2\alpha }}+N_{\alpha }{\tilde{u}}_n\left(x,t\right)\right)\right]\right\}}{\left(ds\right)}^{\alpha }\end{array}$$

(50)

with the initial boundary conditions:

$$u\left(0,t\right)=p_0{E}_{\alpha }\left(-t^{\alpha }\right)$$

(51)

$$\frac{{\partial }^{\alpha }u\left(0,t\right)}{\partial t^{\alpha }}=0.$$

(52)

Considering the iterative formula:

$$u_{n+1}(x,t)=L_x^{\left(-2\alpha \right)}\left(-\frac{1}{c}\frac{{\partial }^{\alpha }{u}_n\left(x,t\right)}{\partial t^{\alpha }}+\frac{A}{c}\frac{{\partial }^{2\alpha }{u}_n\left(x,t\right)}{\partial t^{2\alpha }}+Bu\left(x,t\right)\right)$$

(53)

together with the initial value, yields:

$$u_0(x,t)=p_0\left(E_{\alpha }\left(-t^{\alpha }\right)+L_x^{\left(-2\alpha \right)}D\right)=p_0{E}_{\alpha }\left(-t^{\alpha }\right)+D\frac{x^{2\alpha }}{\Gamma \left(1+2\alpha \right)}.$$

(54)

Integrating Equation (54) into Equation (53) yields:

$$u_1(x,t)=L_x^{\left(-2\alpha \right)}\left(\frac{1}{c}\frac{{\partial }^{\alpha }{u}_0\left(x,t\right)}{\partial t^{\alpha }}+\frac{A}{c}\frac{{\partial }^{2\alpha }{u}_0\left(x,t\right)}{\partial t^{2\alpha }}+Bu\left(x,t\right)\right)=\left(\frac{-1+A}{c}+B\right)\frac{x^{2\alpha }}{\Gamma \left(1+2\alpha \right)}{E}_{\alpha }\left(-t^{\alpha }\right)p_0$$

(55)

$$u_2(x,t)=L_x^{\left(-2\alpha \right)}\left(\frac{1}{c}\frac{{\partial }^{\alpha }{u}_0\left(x,t\right)}{\partial t^{\alpha }}+\frac{A}{c}\frac{{\partial }^{2\alpha }{u}_0\left(x,t\right)}{\partial t^{2\alpha }}+Bu\left(x,t\right)\right)={\left(\frac{-1+A}{c}+B\right)}^2\frac{x^{4\alpha }}{\Gamma \left(1+4\alpha \right)}{E}_{\alpha }\left(-t^{\alpha }\right)p_0$$

(56)

$$u_2(x,t)=L_x^{\left(-2\alpha \right)}\left(\frac{1}{c}\frac{{\partial }^{\alpha }{u}_0\left(x,t\right)}{\partial t^{\alpha }}+\frac{A}{c}\frac{{\partial }^{2\alpha }{u}_0\left(x,t\right)}{\partial t^{2\alpha }}+Bu\left(x,t\right)\right)={\left(\frac{-1+A}{c}+B\right)}^2\frac{x^{4\alpha }}{\Gamma \left(1+4\alpha \right)}{E}_{\alpha }\left(-t^{\alpha }\right)p_0$$

(57)

$$u_4(x,t)=L_x^{\left(-2\alpha \right)}\left(\frac{1}{c}\frac{{\partial }^{\alpha }{u}_0\left(x,t\right)}{\partial t^{\alpha }}+\frac{A}{c}\frac{{\partial }^{2\alpha }{u}_0\left(x,t\right)}{\partial t^{2\alpha }}+Bu\left(x,t\right)\right)={\left(\frac{-1+A}{c}+B\right)}^4\frac{x^{8\alpha }}{\Gamma \left(1+8\alpha \right)}{E}_{\alpha }\left(-t^{\alpha }\right)p_0$$

(58)

$$u_n(x,t)=p_0{\left(\frac{-1+A}{c}+B\right)}^n\frac{x^{2n\alpha }}{\Gamma \left(1+2n\alpha \right)}{E}_{\alpha }\left(-t^{\alpha }\right)$$

(59)

$$\begin{array}{l}u(x,t)=\underset{n\to \infty }{\lim }{\displaystyle \sum _{n=0}^n{u}_n(x,t)}=p_0\left(E_{\alpha }\left(-t^{\alpha }\right){\displaystyle \sum _{n=0}^n{\left(-1\right)}^n{\left(\frac{1-A}{c}-B\right)}^n\frac{x^{2n\alpha }}{\Gamma \left(1+2n\alpha \right)}}\right)+\frac{D{x}^{2\alpha }}{\Gamma \left(1+2\alpha \right)}\\ =p_0\left(E_{\alpha }\left(-t^{\alpha }\right){\cos }_{\alpha }\left(\sqrt{\frac{1-A}{c}-B}x\right)\right)+\frac{D{x}^{2\alpha }}{\Gamma \left(1+2\alpha \right)}\end{array}.$$

(60)

3.3. Final Solution of Water-Phase Pressure

The final fractional analytical solution is obtained, as follows.

For water phase pressure:

$$p_w(x,t)=p_0\left(E_{\alpha }\left(-t^{\alpha }\right){\cos }_{\alpha }\left(\sqrt{\frac{1-A}{c}-B}x\right)\right)+\frac{D{x}^{2\alpha }}{\Gamma \left(1+2\alpha \right)}.$$

(61)

For gas phase pressure:

$$p_g(x,t)=p_0\left(1-A\right)\left(E_{\alpha }\left(-t^{\alpha }\right){\cos }_{\alpha }\left(\sqrt{\frac{1-A}{c}-B}x\right)\right).$$

(62)

4. Parametric Study and Discussions

4.1. Pressure Variation of the Two-Phase Flow

Firstly, fractional analytical solutions are used to predict the two-phase flow pressure variation. The parameters in the calculation are taken from Yang et al. [25,26] and are listed in

Table 1. Computation parameters from a shale gas reservoir in China.

Parameter	Unit	Value	Physical Meanings
$S_{rg}$		0.15	Shale gas residual saturation
$S_{rw}$		0.2	Water residual saturation
${\mu }_w$	Pa*s	$3.6\times {10}^{-4}$	Water viscosity
${\mu }_{\mathrm{g}}$	Pa*s	$2.0\times {10}^{-5}$	Shale gas viscosity
c	m/s	$5.3\times {10}^{-2}$	Travel wave velocity
$\lambda$		1	Pore size distribution index
$p_0$	MPa	27.4	Initial reservoir pressure
$p_1$	MPa	19.67	Well pressure
$k_0$	mD	0.1	Initial shale permeability
${\varphi }_0$		0.18	Initial shale porosity
$p_e$	MPa	2	Entry capillary pressure
$k_{rw}^{\max }$		0.004	End-point relative permeability for water
$k_{rnw}^{\max }$		1	End-point relative permeability for shale gas
$N_w$		2	Water reference parameter
$N_g$		1	Shale gas reference parameter
$\alpha$		$\ln 2/\ln 3$	Fractional dimension
$D_T$		1	Tortuosity fractal dimension
$D_f$		2	Diameter size fractal dimension

.

As shown in Figure 3, the observations of water and gas pressure are presented with red and black edge colors, respectively. It can be seen that the pressure of shale gas and water both decrease over time, respectively. The analytical fractional gas pressure is always higher than water pressure with the same displacement compared with the zero-displacement location in the production time.

In addition, the pressures of shale gas and water both decrease over the period of displacement, respectively. The analytical fractional gas pressure is also always higher than water pressure over the displacement at the same moment in the production time. Taking the position of zero displacement as an example, both the water and the shale gas experience the depressurization process. At the initial moment, the shale gas pressure is 30.94 MPa and the water pressure is 28.93 MPa. By the 100th day, the shale gas pressure is 21.05 MPa and the water pressure is 19.53 MPa, with a reduction magnitude of 31.97% and 32.49%, respectively. Furthermore, the gas and water pressures reach 13.57 MPa and 12.87 MPa by the 300th day, with a reduction magnitude of 56.14% and 55.51%, respectively. The overall downward trend of shale gas and water decreases slowly throughout the production process, but the downward rate of water decreases more slowly. This shows that the effect of the detonation wave is gradually weakened and can be transferred to other forms of wave, with the increase in explosion distance.

4.2. Impact of Fractional Dimension on Gas Pressure

The impact of the fractional dimension on the shale gas depressurization process is investigated in this section. The fractional dimension describes the discontinuity of the two-phase flow path, among the various fracture properties. It is important to study favorable shale gas production conditions and the fracture networks of a reasonable degree after treatment of the shale gas reservoir. This parametric study is intended to explore the effects of the fractional dimension of the flow path on the enhancement of shale gas production.

The fractional dimension is a parameter by which to characterize the non-linear, tortuous, and fractal flow path of the gas-water two-phase flow in the shale reservoir. A sensitivity analysis of the fractional dimension is conducted in this section. The fractional dimension is taken as 0.65, 0.75, and 0.85, respectively. Figure 4 presents the decline rate of shale gas pressure over the first 300 days, in all three cases. Taking the zero-displacement location as an example, the gas pressure on the 100th day is 24.35 MPa, 26.86 MPa, and 28.94 MPa, respectively. On the 300th day, the shale gas pressure is 10.11 MPa, 12.18 MPa, and 12.97 MPa, respectively. These data show that the gas pressure corresponding to a larger fractional dimension drops faster in the gas production stage. The closer the fractional dimension is to 1, the closer the flow path is to a linear and continuous flow, and the faster the movement of the two-phase flow.

4.3. The Impact of Entry Capillary Pressure on Gas Pressure

The occurrence of a gas–water two-phase flow is decided by the entry capillary pressure. The entry pressure varies with rock type and directly depends on the diameter of the flow pipe [4]. It is clear that the value of the entry capillary pressure is within a typical range of 0.1–48.3 MPa [4]. In this study, the impact of entry capillary pressure on shale gas pressure is investigated when the entry capillary pressure is taken as 2, 3, and 4 MPa, respectively.

As shown in Figure 5, the shale gas pressure decreases in all three cases with the increase in entry capillary pressure. After production for 100 days, the gas pressure declines to about 18.73 MPa, 19.67 MPa, and 21.58 MPa, respectively. For example, when the entry pressure Pe = 2 MPa, the gas pressure is lower by 13.2% than that at Pe = 4 MPa. The higher the entry capillary pressure, the faster the gas pressure declines. This is because more gas can flow in a flow pipe of the same diameter. These results show that the increase in entry capillary pressure makes a substantially increasing impact on the shale gas depressurization process.

4.4. Impact of Traveling Wave Velocity on Gas Pressure

At the shale gas and water two-phase flow stage, when a given time and particular location are concerned, the impact of the flow border is too late to reach, and so the wave moves forward and forms a traveling wave.

Three different kinds of traveling wave velocities are considered, to investigate the impact on shale gas pressure. The traveling wave velocity is taken as 0.05, 0.06, and 0.07 m/s, respectively. Figure 6 presents the shale gas variation, with time and displacement. The shale gas pressure shows a decreasing trend in the whole production time, in all three cases. On the 100th day, the shale gas pressure became 24.54 MPa, 23.51 MPa, and 22.47 MPa, respectively. Furthermore, on the 300th day, the rotation angle is 16.57 MPa, 13.48 MPa, and 11.2 MPa, respectively. With the same displacement, the smaller the traveling wave velocity, the faster the gas pressure drops during the same production period. This is reasonable because the traveling wave velocity characterizes the relationship between the time and displacement in the direction of the wave vibration. When the traveling wave velocity is smaller, the wave of the shale gas flow has a weaker displacement in the direction of vibration; all the flow energy is used to propagate the wave along the direction of flow. These results indicate that the traveling wave velocity has a certain degree of influence on the production of shale gas.

5. Conclusions

This study analytically investigated the pressure variations in a two-phase flow under the effects of a nonlinear and discontinuous flow path in a shale gas reservoir. A fully coupled fractional mathematical model of the two-phase flow in shale gas production was developed by considering capillary pressure, the fractal dimension of tortuosity, and diameter size. This coupled mathematical model was then analytically solved using the local fractional traveling wave method and the variational iteration method. The analytical solution of shale gas pressure was employed to predict the shale gas depressurization process from a Chinese shale gas reservoir. The effects of fractional dimension, entry capillary pressure and traveling wave velocity on shale gas pressure were analyzed. Based on these studies, the following conclusions can be drawn.

First, the fractional dimension has a significant impact on the depressurization process of shale gas. A higher fractional dimension leads to a faster decline in gas pressure and a lower shale gas pressure when at the same gas production stage. The closer the fractional dimension is to 1, the closer the gas flow path is to a linear and continuous path, and the faster the shale gas flows.

Second, entry capillary pressure impacts the production of shale gas. The extraction of shale gas is a depressurization process. In our computational case, a higher entry capillary pressure leads to a faster decline in shale gas pressure and a lower shale gas pressure in the same time period.

Third, the traveling wave velocity has a certain degree of influence on the shale gas depressurization process. At the same location, if the traveling wave velocity is smaller, the wave of shale gas flow will have a weaker displacement in the vibration direction of gas flow and all the energy is used to propagate the flow toward the well-head, leading to a faster gas pressure drop within the same production period.