# Chemical Potential-Based Modeling of Shale Gas Transport

Authors: Jisheng Kou; Lingyun Chen; Amgad Salama; Jianchao Cai
Journal: Geofluids
 (2021)

---

## Abstract

Shale gas plays an increasingly important role in the current energy industry. Modeling of gas flow in shale media has become a crucial and useful tool to estimate shale gas production accurately. The second law of thermodynamics provides a theoretical criterion to justify any promising model, but it has been never fully considered in the existing models of shale gas. In this paper, a new mathematical model of gas flow in shale formations is proposed, which uses gas density instead of pressure as the primary variable. A distinctive feature of the model is to employ chemical potential gradient rather than pressure gradient as the primary driving force. This allows to prove that the proposed model obeys an energy dissipation law, and thus, the second law of thermodynamics is satisfied. Moreover, on the basis of energy factorization approach for the Helmholtz free energy density, an efficient, linear, energy stable semi-implicit numerical scheme is proposed for the proposed model. Numerical experiments are also performed to validate the model and numerical method.


---

## Body

## 1. Introduction

Shale gas has become a significant energy resource over the last decade. Shale gas refers to natural gas composed of primarily methane, which is trapped within the pores of fine-grained sedimentary rocks with rich micropores and relatively low permeability. The shale gas reservoirs differ from conventional natural gas reservoirs that apart from free gas in the pores and fractures, a certain amount of gas is adsorbed onto the solid surfaces, and as a result, it can not only store shale gas but also generate gas [1–4]. Experimental investigations have indicated that adsorbed gas storage capacity is primarily affected by shale reservoir conditions, such as temperature, pressure, and shale matrix pore structures [2, 3, 5, 6]. A number of adsorption models have been developed to describe the methane adsorption in shale gas reservoirs. The Langmuir model [7] and Dubinin-Radushkevich (D-R) model [8] are the most popularly employed models to describe the gas adsorption in micropore-rich materials. Methane under shale formation conditions usually stays at the supercritical state, and consequently, the classical models that require a saturation pressure cannot be applied to describe the gas adsorption under supercritical conditions [3, 9]. By the use of gas density rather than gas pressure, modified formulations of the Langmuir and D-R models have been developed for supercritical conditions [3, 9]. In this paper, it is shown that the use of gas density will also be beneficial for ensuring thermodynamical consistency of the models.Modeling of gas flow in shale media plays a crucial role in predicting shale gas production [10–19]. For gas flow in tight porous media, the most remarkable phenomenon is the so-called Klinkenberg effect [20], which results from slip flow of gas molecules through very small pores. This effect leads to the apparent permeability that is generally greater than the absolute permeability of a porous medium [14, 15]. By using the apparent permeability, the shale gas flow equation can be simply formulated as the form of Darcy’s law, which states that gas velocity is proportional to the pressure gradient [12–19, 21].As general principles, the fundamental laws of thermodynamics play a significant role in modeling of various physical problems [22–24]. Specially, the second law of thermodynamics states that any spontaneous process in an isolated system will always lead to an escalation in the entropy of this system. In terms of the second law of thermodynamics, for an isothermal system, any dynamical process should obey an energy dissipation law [22], and a promising model should preserve this property. Nevertheless, energy dissipation laws for modeling of shale gas transport have been scarcely studied so far. In this paper, a new mathematical model is proposed, which uses gas density instead of pressure as the primary variable and introduces chemical potential gradient instead of pressure gradient as the primary driving force. It is rigorously proved that the proposed model obeys an energy dissipation law, and thus, it is thermodynamically consistent (that is, it obeys the second law of thermodynamics).Numerical algorithms that preserve the energy dissipation law at the discrete level, the so-called energy stable methods, are preferred as well [25–27]. In general, it is a quite challenging issue to construct such methods since the commonly used explicit or implicit scheme could not inherit a discrete energy dissipation law. The energy factorization (EF) approach for the Helmholtz free energy density [26] is a newly developed approach to design efficient energy stable numerical methods. An appealing feature of this approach is that it leads to linear, easy-to-implement, and energy stable numerical schemes, and this advantage is more notable for numerical simulation of realistic fluids. Due to its excellent features, this approach has been successfully extended to phase-field models [28, 29]. In this work, using the EF approach, an efficient, linear, energy stable semi-implicit numerical scheme is constructed for the proposed model of shale gas transport.The rest of this paper is organized as follows. In Section2, a chemical potential-based model of gas flow in shale media is proposed, which is proved to obey an energy dissipation law. An efficient, linear, energy stable semi-implicit numerical method is proposed in Section 3. In Section 4, numerical experiments are performed to validate the proposed model and numerical scheme. Finally, some concluding remarks are provided in Section 5.


## 2. Model Equations

In this section, a chemical potential-based model of shale gas transport is proposed. The chemical potential is defined as the derivative of Helmholtz free energy density. In AppendixA, we elaborate on the Helmholtz free energy density determined by Peng-Robinson equation of state [30].Molar density (mol m-3) of methane is denoted by c. For specific temperature TK, the Helmholtz free energy density, denoted by fc, is a function of molar density c as described in Appendix A. The chemical potential is defined as the derivative of Helmholtz free energy density function fc with respect to molar density
(1)μc=f′c,where μ is the chemical potential (Pa·mol-1 m3) and fc is the Helmholtz free energy density (Pa).For specific temperature and pressure, a cubic equation is solved to obtain molar density of methane, and moreover, the solutions ofc are not unique in general. In contrast, for specific temperature and molar density, the pressure can be uniquely and explicitly calculated from the Peng-Robinson equation of state [30], which is formulated in Appendix B. Consequently, under the constant temperature, molar density is preferred to pressure as the primary variable in numerical simulation.For specific temperature, the pressure becomes a function of molar densityc only. Moreover, the pressure relates to the Helmholtz free energy density and chemical potential as [22]
(2)p=cμ−f,where p is the pressure (Pa). The chain rule gives the relationship between the pressure gradient and chemical potential gradient.
(3)∇p=∇cμ−∇f=c∇μ+μ∇c−μ∇c=c∇μ.Letϕ denote the porosity. In tight reservoirs with abundant micropores, the Klinkenberg effect cannot be neglected, and thus, the apparent permeability [14, 15] is expressed as
(4)Kapp=K01+bp,where b is the Klinkenberg slippage factor (Pa) and K0 is the intrinsic permeability (md). The slippage factor accounts for the slippage effect on permeability of gas in reservoirs. Various formulations for the slippage factor have been developed in the literature, and the following formulation [13] is used in this work:
(5)b=ησ8πRTϕMwK0,where η is the gas viscosity (Pa·s), σ is the comprehensive coefficient, R is the universal gas constant (m3·Pa·mol-1·K-1), and Mw is the molar weight of gas (g·mol-1). On the basis of the apparent permeability, the velocity can be described by the form of Darcy’s law as
(6)u=−Kappη∇p,where uis the Darcy velocity (m·s-1) and Kapp is the apparent permeability (md). Applying relation (3), the following chemical potential-based velocity formulation is obtained:
(7)u=−Kappηc∇μ.In the shale gas reservoir, there is a large amount of gas adsorbed onto micropores in addition to free gas. The adsorbed gas has no mobility unless it converts to free gas, and consequently, it is assumed to have no contribution to the free energy. In the context of shale gas,c stands for the free gas density, and meanwhile, cads represents the adsorbed gas molar density. In order to describe the adsorption of gas, there have been several models proposed in the literature, which can be classified into two classes according to their primary variables, namely, the pressure-based models and density-based models. One famous pressure-based model is the classical Langmuir isotherm adsorption model [7]. The pressure-based models are not suitable to describe the gas adsorption under supercritical conditions [3, 9]. The density-based models have been developed and extensively applied to characterize the adsorption that occurs in shale media. Two modified adsorption models using molar density as the primary variable are described in Appendix C, which are employed in this paper due to their consistency to the chemical potential-based model. The mass accumulation of both free and adsorbed gas is given by ctot=ϕc+1−ϕcads, where cads is the adsorbed gas density (mol·m-3). The mass balance equation can be expressed as
(8)∂ctot∂t+∇·uc=0.Furthermore, a density-dependent function is defined as(9)φc=∂ctot∂c=ϕ+1−ϕ∂cads∂c,and then, (8) can be rewritten as
(10)φc∂c∂t+∇·uc=0.For an isothermal dynamical system, the second law of thermodynamics leads to a certain energy dissipation law [22]. A distinctive feature of the proposed model is that it obeys an energy dissipation law, and as a consequence, it is thermodynamically consistent. Let Ω be a connected and smooth space domain with the boundary ∂Ω, and further, let n denote the normal unit outward vector to ∂Ω. Multiplying (10) by μc and then integrating it over the domain, it follows that
(11)∫Ωμcφc∂c∂tdx+∫Ωμc∇·ucdx=0.The definition of chemical potential and integration by parts lead to(12)∫Ωφc∂fc∂tdx+∫∂Ωμccu·nds−∫Ωu·c∇μcdx=0.Substituting (7) into (12) yields
(13)∫Ωφc∂fc∂tdx+∫∂Ωμccu·nds=−∫ΩKappηc2∇μ2dx.The right-hand side term of (13) represents the energy-dissipation rate. In order to demonstrate this property clearly, the adsorption effect is disregarded and the no-flow boundary condition is applied on the entire boundary of the domain, i.e., u·n=0. In this situation, total free energy within the domain is defined as
(14)Ft=∫Ωϕfcdx,and from (13), the following energy dissipation law can be derived:
(15)∂Ft∂t=−∫ΩKappηc2∇μ2dx≤0,which implies that total free energy would be dissipated over time until an equilibrium state is reached.


## 3. Numerical Method

In this section, a linear, efficient numerical scheme is proposed, and it is able to ensure the energy dissipation law at the discrete level. Denote byτ the time step size and also denote by tn=nτ the time. In the time discrete scheme, cn stands for the approximation of molar density c at time tn.


### 3.1. Energy Factorization Approach

The proposed scheme employs the EF approach for the Helmholtz free energy density [26], which leads to the linear discrete chemical potential inheriting the energy dissipation property. The Helmholtz free energy density can be expressed as a sum of three contributions
(16)fc=fidealc+frepulsionc+fattractionc,where
(17)fidealc=cRTlnc,(18)frepulsionc=−cRTln1−βc,(19)fattractionc=αTc22βln1+1−2βc1+1+2βc.Here,R is the universal gas constant and the parameters α and β are described in Appendix A. The energy factorization approach gives the discrete chemical potential as
(20)μn+1=μidealn+1+μrepulsionn+1+μattractionn,where
(21)μidealn+1=RTlncn+RTcn+1cn,(22)μrepulsionn+1=RTG′cn2Gcn+G′cncn+1−cn−λRT,(23)μattractionn=αT22βln1+1−2βcn+1−2βcn1+1−2βcn−αT22βln1+1+2βcn+1+2βcn1+1+2βcn.In (22), Gc is an intermediate energy function, which is defined as
(24)Gc=λc−cln1−βc,where λ is the dimensionless stabilization constant. Apparently, the discrete chemical potential μn+1 is a linear function of cn+1 only. From the physical point of view, the boundedness of molar density c is assumed as
(25)0<c≤ϱ,ϱβ<1,where ρ is the upper bound of molar density. Let ε0=ϱβ. If the stabilization constant λ in (24) is chosen such that
(26)λ≥ϵ01−ϵ02+ϵ021−ϵ04−2ln1−ϵ0ε01−ϵ021/2.It has been proved in [26] that the discrete chemical potential (20) satisfies the following energy inequality:
(27)fcn+1−fcn≤μn+1cn+1−cn.


### 3.2. Discrete Scheme

On the basis of the discrete chemical potential given in (20)–(23), the proposed semi-implicit time discrete scheme reads as follows:
(28)φcncn+1−cnτ+∇·un+1cn=0,(29)un+1=−Kappnηcn∇μn+1.The apparent permeability is explicitly calculated as(30)Kappn=K01+bpcn,where the pressure can be directly obtained by (2) using cn. Combining (28) and (29) yields a single linear equation of cn+1 as
(31)φcncn+1−cnτ−∇·Kappnηcn2∇μn+1=0,which is easy to be solved due to the fact that μn+1 is a linear function of cn+1.As depicted in AppendixC, it is assumed that φc>0. It is ready to prove that the proposed numerical scheme preserves a discrete energy dissipation law. Multiplying Equations (28) and (29) by μn+1 and cn∇μn+1, respectively, and then integrating them over the domain Ω, we obtain
(32)∫Ωμn+1φcncn+1−cnτdx+∫Ωμn+1∇·un+1cndx=0,(33)∫Ωun+1·cn∇μn+1dx=−∫ΩKappnηcn∇μn+12dx.It follows from (32) and (33) that
(34)∫Ωμn+1φcncn+1−cnτdx=−∫Ωμn+1∇·un+1cndx=∫Ωun+1·cn∇μn+1dx−∫∂Ωμn+1cnun+1·nds=−∫ΩKappnηcn∇μn+12dx−∫∂Ωμn+1cnun+1·nds.On the other hand, it follows from (17) that
(35)∫Ωμn+1φcncn+1−cndx≥∫Ωφcnfcn+1−fcndx.As a consequence of (34) and (35), the following discrete energy dissipation law is obtained
(36)∫Ωφcnfcn+1−fcnτdx+∫∂Ωμn+1cnun+1·nds≤−∫ΩKappnηcn∇μn+12dx.Equation (31) can be further spatially discretized using the cell-centered finite difference (CCFD) method, which is a popular discrete method in various fluid flow fields [25, 26].


## 3.1. Energy Factorization Approach

The proposed scheme employs the EF approach for the Helmholtz free energy density [26], which leads to the linear discrete chemical potential inheriting the energy dissipation property. The Helmholtz free energy density can be expressed as a sum of three contributions
(16)fc=fidealc+frepulsionc+fattractionc,where
(17)fidealc=cRTlnc,(18)frepulsionc=−cRTln1−βc,(19)fattractionc=αTc22βln1+1−2βc1+1+2βc.Here,R is the universal gas constant and the parameters α and β are described in Appendix A. The energy factorization approach gives the discrete chemical potential as
(20)μn+1=μidealn+1+μrepulsionn+1+μattractionn,where
(21)μidealn+1=RTlncn+RTcn+1cn,(22)μrepulsionn+1=RTG′cn2Gcn+G′cncn+1−cn−λRT,(23)μattractionn=αT22βln1+1−2βcn+1−2βcn1+1−2βcn−αT22βln1+1+2βcn+1+2βcn1+1+2βcn.In (22), Gc is an intermediate energy function, which is defined as
(24)Gc=λc−cln1−βc,where λ is the dimensionless stabilization constant. Apparently, the discrete chemical potential μn+1 is a linear function of cn+1 only. From the physical point of view, the boundedness of molar density c is assumed as
(25)0<c≤ϱ,ϱβ<1,where ρ is the upper bound of molar density. Let ε0=ϱβ. If the stabilization constant λ in (24) is chosen such that
(26)λ≥ϵ01−ϵ02+ϵ021−ϵ04−2ln1−ϵ0ε01−ϵ021/2.It has been proved in [26] that the discrete chemical potential (20) satisfies the following energy inequality:
(27)fcn+1−fcn≤μn+1cn+1−cn.


## 3.2. Discrete Scheme

On the basis of the discrete chemical potential given in (20)–(23), the proposed semi-implicit time discrete scheme reads as follows:
(28)φcncn+1−cnτ+∇·un+1cn=0,(29)un+1=−Kappnηcn∇μn+1.The apparent permeability is explicitly calculated as(30)Kappn=K01+bpcn,where the pressure can be directly obtained by (2) using cn. Combining (28) and (29) yields a single linear equation of cn+1 as
(31)φcncn+1−cnτ−∇·Kappnηcn2∇μn+1=0,which is easy to be solved due to the fact that μn+1 is a linear function of cn+1.As depicted in AppendixC, it is assumed that φc>0. It is ready to prove that the proposed numerical scheme preserves a discrete energy dissipation law. Multiplying Equations (28) and (29) by μn+1 and cn∇μn+1, respectively, and then integrating them over the domain Ω, we obtain
(32)∫Ωμn+1φcncn+1−cnτdx+∫Ωμn+1∇·un+1cndx=0,(33)∫Ωun+1·cn∇μn+1dx=−∫ΩKappnηcn∇μn+12dx.It follows from (32) and (33) that
(34)∫Ωμn+1φcncn+1−cnτdx=−∫Ωμn+1∇·un+1cndx=∫Ωun+1·cn∇μn+1dx−∫∂Ωμn+1cnun+1·nds=−∫ΩKappnηcn∇μn+12dx−∫∂Ωμn+1cnun+1·nds.On the other hand, it follows from (17) that
(35)∫Ωμn+1φcncn+1−cndx≥∫Ωφcnfcn+1−fcndx.As a consequence of (34) and (35), the following discrete energy dissipation law is obtained
(36)∫Ωφcnfcn+1−fcnτdx+∫∂Ωμn+1cnun+1·nds≤−∫ΩKappnηcn∇μn+12dx.Equation (31) can be further spatially discretized using the cell-centered finite difference (CCFD) method, which is a popular discrete method in various fluid flow fields [25, 26].


## 4. Numerical Examples

In this section, the performance of the proposed model and numerical method is demonstrated through some numerical tests. The spatial domain in all numerical tests is a square domain0,10×0,10, and the length unit is meter. A uniform rectangular mesh with 3600 grid cells is used to divide the domain. The gas is methane and its physical properties are listed in Table 1. The modified Dubinin’s adsorption model is used, and the initial velocity is taken to be zero for all tests.Table 1
Physical properties of methane.
Pc (bar)Tc (K)Acentric factorMw (g/mole)45.99190.560.01116.04The Hagen-Poiseuille equation is used to calculate the intrinsic permeability of gas flow in cylinders [13].
(37)K0=r28,where r refers to the averaged pore diameter (nm).


### 4.1. Example 1: Verification of the Energy Dissipation Property

In this example, the aim is to verify the energy dissipation property of the model and numerical scheme. For this purpose, the adsorption of methane is disregarded and the no-flow boundary condition is applied on the entire boundary of the domain, i.e.,u·n=0. In this situation, total free energy is expressed as
(38)Fn=∫Ωϕfcndx.From (36), the following energy dissipation law at the discrete level is obtained
(39)Fn+1−Fn≤−τ∫ΩKappnηcn∇μn+12dx≤0,which implies that total free energy would be decreasing with time steps until an equilibrium state is reached. The physical parameters used in this example are listed in Table 2.Table 2
Parameters of the model used in example 1.
ParameterValueUnitDescriptionϕ0.05—Porosityr10nmAveraged pore diameterR8.3144621m3·Pa·mol−1·K−1Ideal gas constantT300KTemperatureη10−5Pa·sViscosityσ10−3—Comprehensive coefficientThe stabilization constant in (24) is taken as λ=0.0375. The time step size is chosen as τ=20hours, and 50 time steps are simulated. The initial methane density is generated by a random way.The proposed model and numerical scheme are employed to simulate this problem. Dynamics of molar density and pressure is illustrated in Figures1 and 2. For comparison, the classical pressure-based model is also applied to simulate the same problem, and the corresponding molar density and pressure are shown in Figures 3 and 4. The classical model are described in detail in Appendix D. Although two models produce slightly different results, it is clearly observed that due to spatial inhomogeneity of molar density, the fluid system tends to reach an equilibrium state where molar density becomes uniform in space.Figure 1
Gas density distributions computed by the proposed model in example 1 at different times.
(a)
Initial(b)
100 hours(c)
1000 hoursFigure 2
Pressure distributions computed by the proposed model in example 1 at different times.
(a)
Initial(b)
100 hours(c)
1000 hoursFigure 3
Gas density distributions computed by the classical model in example 1 at different times.
(a)
Initial(b)
100 hours(c)
1000 hoursFigure 4
Pressure distributions computed by the classical model in example 1 at different times.
(a)
Initial(b)
100 hours(c)
1000 hoursFigure5 depicts total free energy profiles of the proposed model and numerical scheme. Total energy is always decreasing monotonously over time, and thus, theoretical result is validated. For comparison, total free energy computed by the classical model is also plotted in Figure 6. Despite the decrease of total energy occurring in the simulation of the classical model, total energy is obviously oscillating and not monotonously decreasing over time as shown in the zoom-in plot.Figure 5
Total energy profiles of the proposed model in example 1.
(a)(b)Figure 6
Total energy profiles of the classical model in example 1.
(a)(b)The performance of two models in preserving the property of total mole conservation is also compared. The relative error of total moles is defined as(40)Δcn=∫Ωcndx−∫Ωc0dx∫Ωc0dx.Figure7 depicts the relative errors of total moles computed by two models. It can be observed that the proposed model can guarantee this key property very well within the range of roundoff errors, whereas the classical pressure-based model can cause a nonnegligible mass loss even in this simple problem. This means that the use of molar density as the primary variable has a great advantage in guaranteeing mass conservation.Figure 7
Comparison of total moles conservation in example 1.
(a)
Proposed model(b)
Classical model


### 4.2. Example 2

A gas reservoir with multiple high-permeable layers is considered in this example. The intrinsic permeability of the reservoir is illustrated in Figure8. The layers with high permeability are denoted by
(41)Ωm=x=x,y:x∈0,Lx,y∈0.2Ly,0.3Ly∪0.45Ly,0.55Ly∪0.7Ly,0.8Ly,where Lx=Ly=10 m. The initial gas density is uniform and equal to 200 mol/m3, while the initial adsorption density is illustrated in Figure 9. The gas is flowing out from the production boundary Γ, which locates on the right end of the domain
(42)Γ=x=x,y:x=Lx,y∈0.2Ly,0.3Ly∪0.45Ly,0.55Ly∪0.7Ly,0.8Ly.Figure 8
Intrinsic permeability in example 2.Figure 9
Initial adsorption density in example 2.The boundary condition on the production end is specified as(43)cx,y,t=cb,x,y∈Γ,wherecb is the gas production density. The rest boundaries are closed, and no mass transfer across these boundaries takes place. The physical parameters used in this example are listed in Table 3. The stabilization constant in (24) is taken as λ=0.0148, and the time step size is chosen as τ=0.8 hour.Table 3
Parameters of the model used in example 2.
ParameterValueUnitDescriptionϕ0.02, 0.05—Porosityr10, 200nmAveraged pore diameterR8.3144621m3·Pa·mol−1·K−1Ideal gas constantT360KTemperatureD0.1, 0.2—Pore structure parameterη10−5Pa·sViscosityσ0.82—Comprehensive coefficientca322,302mol·m−3Adsorbed gas phase densitycM322,302mol·m−3Maximum absolute adsorption densityTests withcb=20,40,60 mol/m3 are performed. The gas density distributions at various times computed with cb=20,40,60 mol/m3 are illustrated in Figures 10–12, while the corresponding pressure contours are shown in Figures 13–15. The cumulative gas production profiles under different values of cb are also plotted in Figure 16. These results indicate that the production density has a great effect on the gas production; in fact, the increase of cb can largely reduce the cumulative gas production. The adsorption densities at different times computed with cb=20,40,60 mol/m3 are shown in Figures 17–19. The total adsorption amounts under different values of the production density are illustrated in Figure 20. It is clearly observed that the remaining total adsorption amounts are against the increase of production density. This means that a large production density may not be preferred for the sake of achieving optimal gas production.Figure 10
Gas density distributions in example 2 withcb=20 mol/m3 at different times.
(a)
20 days(b)
60 days(c)
200 daysFigure 11
Gas density distributions in example 2 withcb=40 mol/m3 at different times.
(a)
20 days(b)
60 days(c)
200 daysFigure 12
Gas density distributions in example 2 withcb=60 mol/m3 at different times.
(a)
20 days(b)
60 days(c)
200 daysFigure 13
Pressure distributions in example 2 withcb=20 mol/m3 at different times.
(a)
20 days(b)
60 days(c)
200 daysFigure 14
Pressure distributions in example 2 withcb=40 mol/m3 at different times.
(a)
20 days(b)
60 days(c)
200 daysFigure 15
Pressure distributions in example 2 withcb=60 mol/m3 at different times.
(a)
20 days(b)
60 days(c)
200 daysFigure 16
The cumulative production under different values ofcb.Figure 17
Adsorption density distributions in example 2 withcb=20 mol/m3 at different times.
(a)
20 days(b)
60 days(c)
200 daysFigure 18
Adsorption density distributions in example 2 withcb=40 mol/m3 at different times.
(a)
20 days(b)
60 days(c)
200 daysFigure 19
Adsorption density distributions in example 2 withcb=60 mol/m3 at different times.
(a)
20 days(b)
60 days(c)
200 daysFigure 20
Total adsorption curves under different values ofcb.


## 4.1. Example 1: Verification of the Energy Dissipation Property

In this example, the aim is to verify the energy dissipation property of the model and numerical scheme. For this purpose, the adsorption of methane is disregarded and the no-flow boundary condition is applied on the entire boundary of the domain, i.e.,u·n=0. In this situation, total free energy is expressed as
(38)Fn=∫Ωϕfcndx.From (36), the following energy dissipation law at the discrete level is obtained
(39)Fn+1−Fn≤−τ∫ΩKappnηcn∇μn+12dx≤0,which implies that total free energy would be decreasing with time steps until an equilibrium state is reached. The physical parameters used in this example are listed in Table 2.Table 2
Parameters of the model used in example 1.
ParameterValueUnitDescriptionϕ0.05—Porosityr10nmAveraged pore diameterR8.3144621m3·Pa·mol−1·K−1Ideal gas constantT300KTemperatureη10−5Pa·sViscosityσ10−3—Comprehensive coefficientThe stabilization constant in (24) is taken as λ=0.0375. The time step size is chosen as τ=20hours, and 50 time steps are simulated. The initial methane density is generated by a random way.The proposed model and numerical scheme are employed to simulate this problem. Dynamics of molar density and pressure is illustrated in Figures1 and 2. For comparison, the classical pressure-based model is also applied to simulate the same problem, and the corresponding molar density and pressure are shown in Figures 3 and 4. The classical model are described in detail in Appendix D. Although two models produce slightly different results, it is clearly observed that due to spatial inhomogeneity of molar density, the fluid system tends to reach an equilibrium state where molar density becomes uniform in space.Figure 1
Gas density distributions computed by the proposed model in example 1 at different times.
(a)
Initial(b)
100 hours(c)
1000 hoursFigure 2
Pressure distributions computed by the proposed model in example 1 at different times.
(a)
Initial(b)
100 hours(c)
1000 hoursFigure 3
Gas density distributions computed by the classical model in example 1 at different times.
(a)
Initial(b)
100 hours(c)
1000 hoursFigure 4
Pressure distributions computed by the classical model in example 1 at different times.
(a)
Initial(b)
100 hours(c)
1000 hoursFigure5 depicts total free energy profiles of the proposed model and numerical scheme. Total energy is always decreasing monotonously over time, and thus, theoretical result is validated. For comparison, total free energy computed by the classical model is also plotted in Figure 6. Despite the decrease of total energy occurring in the simulation of the classical model, total energy is obviously oscillating and not monotonously decreasing over time as shown in the zoom-in plot.Figure 5
Total energy profiles of the proposed model in example 1.
(a)(b)Figure 6
Total energy profiles of the classical model in example 1.
(a)(b)The performance of two models in preserving the property of total mole conservation is also compared. The relative error of total moles is defined as(40)Δcn=∫Ωcndx−∫Ωc0dx∫Ωc0dx.Figure7 depicts the relative errors of total moles computed by two models. It can be observed that the proposed model can guarantee this key property very well within the range of roundoff errors, whereas the classical pressure-based model can cause a nonnegligible mass loss even in this simple problem. This means that the use of molar density as the primary variable has a great advantage in guaranteeing mass conservation.Figure 7
Comparison of total moles conservation in example 1.
(a)
Proposed model(b)
Classical model


## 4.2. Example 2

A gas reservoir with multiple high-permeable layers is considered in this example. The intrinsic permeability of the reservoir is illustrated in Figure8. The layers with high permeability are denoted by
(41)Ωm=x=x,y:x∈0,Lx,y∈0.2Ly,0.3Ly∪0.45Ly,0.55Ly∪0.7Ly,0.8Ly,where Lx=Ly=10 m. The initial gas density is uniform and equal to 200 mol/m3, while the initial adsorption density is illustrated in Figure 9. The gas is flowing out from the production boundary Γ, which locates on the right end of the domain
(42)Γ=x=x,y:x=Lx,y∈0.2Ly,0.3Ly∪0.45Ly,0.55Ly∪0.7Ly,0.8Ly.Figure 8
Intrinsic permeability in example 2.Figure 9
Initial adsorption density in example 2.The boundary condition on the production end is specified as(43)cx,y,t=cb,x,y∈Γ,wherecb is the gas production density. The rest boundaries are closed, and no mass transfer across these boundaries takes place. The physical parameters used in this example are listed in Table 3. The stabilization constant in (24) is taken as λ=0.0148, and the time step size is chosen as τ=0.8 hour.Table 3
Parameters of the model used in example 2.
ParameterValueUnitDescriptionϕ0.02, 0.05—Porosityr10, 200nmAveraged pore diameterR8.3144621m3·Pa·mol−1·K−1Ideal gas constantT360KTemperatureD0.1, 0.2—Pore structure parameterη10−5Pa·sViscosityσ0.82—Comprehensive coefficientca322,302mol·m−3Adsorbed gas phase densitycM322,302mol·m−3Maximum absolute adsorption densityTests withcb=20,40,60 mol/m3 are performed. The gas density distributions at various times computed with cb=20,40,60 mol/m3 are illustrated in Figures 10–12, while the corresponding pressure contours are shown in Figures 13–15. The cumulative gas production profiles under different values of cb are also plotted in Figure 16. These results indicate that the production density has a great effect on the gas production; in fact, the increase of cb can largely reduce the cumulative gas production. The adsorption densities at different times computed with cb=20,40,60 mol/m3 are shown in Figures 17–19. The total adsorption amounts under different values of the production density are illustrated in Figure 20. It is clearly observed that the remaining total adsorption amounts are against the increase of production density. This means that a large production density may not be preferred for the sake of achieving optimal gas production.Figure 10
Gas density distributions in example 2 withcb=20 mol/m3 at different times.
(a)
20 days(b)
60 days(c)
200 daysFigure 11
Gas density distributions in example 2 withcb=40 mol/m3 at different times.
(a)
20 days(b)
60 days(c)
200 daysFigure 12
Gas density distributions in example 2 withcb=60 mol/m3 at different times.
(a)
20 days(b)
60 days(c)
200 daysFigure 13
Pressure distributions in example 2 withcb=20 mol/m3 at different times.
(a)
20 days(b)
60 days(c)
200 daysFigure 14
Pressure distributions in example 2 withcb=40 mol/m3 at different times.
(a)
20 days(b)
60 days(c)
200 daysFigure 15
Pressure distributions in example 2 withcb=60 mol/m3 at different times.
(a)
20 days(b)
60 days(c)
200 daysFigure 16
The cumulative production under different values ofcb.Figure 17
Adsorption density distributions in example 2 withcb=20 mol/m3 at different times.
(a)
20 days(b)
60 days(c)
200 daysFigure 18
Adsorption density distributions in example 2 withcb=40 mol/m3 at different times.
(a)
20 days(b)
60 days(c)
200 daysFigure 19
Adsorption density distributions in example 2 withcb=60 mol/m3 at different times.
(a)
20 days(b)
60 days(c)
200 daysFigure 20
Total adsorption curves under different values ofcb.


## 5. Conclusions

A new mathematical model of gas flow in shale media has been proposed. Different from the existing models, the proposed model uses gas density instead of pressure as the primary variable, and it employs chemical potential gradient rather than pressure gradient as the primary driving force. This distinctive feature brings up with thermodynamical consistency of the proposed model; that is, the model obeys an energy dissipation law, which implies the satisfaction of the second law of thermodynamics. Physically, energy stable numerical methods that preserve a discrete energy dissipation law are strongly demanded for numerical approximation of the model. For this purpose, on the basis of the energy factorization approach, an efficient, linear, energy stable semi-implicit numerical scheme is proposed for the proposed model, which inherits an energy dissipation law at the discrete level. Numerical experiments are performed, and numerical results show that the model and numerical method can not only ensure energy stability but also provide physically reasonable results.



---

*Source: 2535284-2021-10-04.xml*