Next Article in Journal
A Survey of Renewable Energy, Climate Change, and Policy Awareness in Israel: The Long Path for Citizen Participation in the National Renewable Energy Transition
Next Article in Special Issue
Microscopic Study of Shale Anisotropy with SEM In Situ Compression and Three-Point Bending Experiments
Previous Article in Journal
A Review of the Optimization Strategies and Methods Used to Locate Hydrogen Fuel Refueling Stations
Previous Article in Special Issue
Impact of Terrigenous Organic Matter Input on Organic Matter Enrichment of Paleocene Source Rocks, Lishui Sag, East China Sea
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Study of the Non-Linear Seepage Problem in Porous Media via the Homotopy Analysis Method

State Key Laboratory of Petroleum Resources and Prospecting, College of Petroleum Engineering, China University of Petroleum-Beijing, Beijing 102249, China
*
Author to whom correspondence should be addressed.
Energies 2023, 16(5), 2175; https://doi.org/10.3390/en16052175
Submission received: 27 January 2023 / Revised: 16 February 2023 / Accepted: 21 February 2023 / Published: 23 February 2023
(This article belongs to the Special Issue Advances in Petroleum Geology and Unconventional Oil and Gas)

Abstract

:
A non-Darcy flow with moving boundary conditions in a low-permeability reservoir was solved using the homotopy analysis method (HAM), which was converted into a fixed-boundary mathematical model via similarity transformation. Approximate analytical solutions based on the HAM are guaranteed to be more accurate than exact analytical solutions, with relative errors between 0.0089% and 2.64%. When λ = 0, the pressure drop of the Darcy seepage model could be instantaneously transmitted to infinity. When λ > 0, the pressure drop curve of the non-Darcy seepage model exhibited the characteristics of tight support, which was clearly different from the Darcy seepage model’s formation pressure distribution curve. According to the results of the HAM, a moving boundary is more influenced by threshold pressure gradients with a longer time. When the threshold pressure gradients were smaller, the moving boundaries move more quickly and are more sensitive to external influences. One-dimensional, low-permeability porous media with a non-Darcy flow with moving boundary conditions can be reduced to a Darcy seepage model if the threshold pressure gradient values tend to zero.

1. Introduction

With the reduction in traditional oil resources, global oil companies have increasingly focused on developing heavy oil reservoirs and tight reservoirs (such as shale reservoirs). As a result of the boundary layer effects and fluid dynamics, theories [1,2,3,4,5] and experiments [6,7,8] have shown that there are threshold pressure gradients [8] when a fluid flows through low-permeability porous media.
Studies of non-Darcy seepage models in porous media considering threshold pressure gradients [9,10,11,12,13,14,15,16] have many practical applications in the development of low-permeability oil and gas reservoirs, hydraulic fracturing proppant suspensions, and polymer-enhanced oil recovery.The existence conditions of threshold pressure gradients of gas flow in tight gas reservoirs with water was investigated experimentally by Zhu et al. [9]. Using laboratory tests, the threshold pressure gradient, water saturation, and absolute permeability were determined, and the threshold pressure gradient was proven. With the law of conservation of mass and momentum, Yao et al. [10] developed a comprehensive mathematical model of gas flow and adsorption in adsorption beds. An analysis of the typical cycle process of axial flow rapid pressure swing adsorption was performed numerically. A new physics-based non-Darcy equation for low-velocity flow was presented by Zhao et al. [11]. A modified Buckingham-Rehner equation was used in order to derive the new equation strictly. Darcy equations are modified by using boundary layer parameters and non-Darcy coefficients. Wei et al. [12] proposed A-B Swartzendruber model connecting non-Darcy flows in porous media at different scales. Laplace transform was used to obtain the analytical solution, and provided a unified description of non-Darcy flow in low-permeability and high- permeability porous media. A one-dimensional non-Darcy flow originating from linear sources was studied by Zhou et al. [13] by taking into account the bilinear relationship between fluid velocity and pressure gradient. A comparison of non-Darcy flow with threshold pressure gradient and non-Darcy flow with bilinear relationship was conducted, and the error caused by neglecting permeability under low pressure gradients was investigated. The Izbash equation was used to construct a radial two-zone composite non-Darcy flow model by Nie et al. [14]. For simulating Darcy flow, the outer region was considered as a homogeneous medium, whereas the inner region was considered as a dual porous medium. By combining Laplace transform, linearization, and Stehfest numerical inversion, this model could be used to study pressure transient behavior of non-Darcy flows in homogeneous systems. Moghimi et al. [15] found the porosity of the foam structure and the number of pores per inch influence non-Darcy flow initiation. A study of the Forchheimer coefficient and permeability value of the system was also conducted. Despite low Reynolds numbers, non-Darcy flow was observed due to the complexity of pore flow mode. The impact of interlayer interfaces on Darcy and non-Darcy flow characteristics in layered porous media was investigated by Zhang et al. [16] using pore scale flow data. Darcy permeability of layers of porous media could be estimated more accurately using the effective permeability factoring in interface effect. The moving boundary models of non-Darcy flows in low-permeability reservoirs have gained a considerable amount of attention in recent years [17,18,19,20,21,22,23,24,25,26]. At present, the research methods for such moving boundary problems mainly include analytical methods [17,18,19,20,21] and numerical methods [22,23,24,25,26]. In a comprehensive study of non-Darcy single-phase flows, Chen et al. [17] examined the effects of threshold pressure gradients on the distribution of pressure at the free boundaries. By using the similarity transformation method and the finite difference method, non-Darcy seepage models in porous media with low permeability were numerically and analytically solved by Liu et al. [18,19,20,21,24,25,26]. To solve the problem of low permeability with a moving boundary, Guo et al. [22] produced a non-dimensional seepage control equation, the initial conditions, and the boundary conditions using a non-grid numerical simulation. Cheng et al. [23] developed numerical models to test low-velocity non-Darcy seepage in micro-compressors with threshold pressure gradients and moving boundaries, discretized the time and space variables, and obtained the numerical solutions. Permeability was numerically simulated using non-Darcy models by Xu et al. [27]. In the new calculation of the relative permeability of oil and water, the value was slightly increased compared with that in Darcy seepage models.
There are many non-linear problems in the mechanics of petroleum engineering. It is of great scientific value to obtain approximate analytical solutions to non-linear equations. Traditional analytic approximate methods cannot provide effective solutions for all the physical parameters, which means that, in essence, they are only applicable to weak non-linear problems. The homotopy analysis method (HAM) is widely used to solve strong non-linear problems in many different fields [28,29,30,31,32,33,34,35,36,37,38,39,40,41,42]. Unlike the perturbation method, the results of HAM do not depend on small or large physical parameters. A further advantage is that the homotopy analysis method (HAM) allows for a choice in the form of the expression for higher-order approximate series solutions. Simple control can also be obtained over the convergence of a series solution using the homotopy analysis method (HAM).
The consideration of threshold pressure gradients in non-Darcy models of porous media is a moving boundary problem. The exact analytical solutions of moving boundary models are difficult to obtain unless the mathematical model is transformed equivalently by similarity variables. For flows in porous media with threshold pressure gradients, there are few studies on the exact analytical solutions of correlative moving boundary models. It is also more complicated to solve this type of problem using numerical methods, fractal modeling, and pore scale network modeling. In this study, by applying the similarity transformation [43] widely used in heat transfer systems to solve classical Stefan moving boundary problems [44], we present approximate analytical solutions based on the homotopy analysis method for one-dimensional non-Darcy flow in porous media with low permeability. In addition to the introduction in Section 1, Section 2 contains a mathematical description of the problem and a derivative using the homotopy analysis method (HAM). Section 3 includes the results and discussion, covering graphic illustrations and tables. Finally, Section 4 contains the conclusions, highlighting the main findings of this work.

2. Mathematical Models and Methods

The research object of this study was to investigate non-Darcy models in one-dimensional, low-permeability porous media with threshold pressure gradients [16]. It was assumed that: (1) the flow had one phase, and was slightly compressible and not influenced by gravity, considering a low-pressure gradient for their formation; and (2) the porous media were homogeneous, isotropic, isothermal, and slightly compressible.
The state equation for the density of a fluid and the porosity of rock [16] is
ρ = ρ i e C f ( p ¯ i p ¯ ) ,   ϕ = ϕ i e C ϕ ( p ¯ i p ¯ ) ,
where ρ is the density of the fluid, ρ i is the initial density of the fluid, p ¯ is the formation pressure, p ¯ i is the original formation pressure, ϕ is the porosity, ϕ i is the initial porosity, C f is the fluid compression coefficient, and C ϕ is the pore compression coefficient.
When the threshold pressure gradient is taken into account, the non-Darcy flow equation for porous media with low permeability is
v ¯ = { 0 , 0 | p ¯ x ¯ | λ ¯ k μ p ¯ x ¯ ( 1 λ ¯ | p ¯ x ¯ | ) , | p ¯ x ¯ | > λ ¯ ,
where v ¯ is the seepage velocity, k is the permeability, μ is the fluid viscosity, x ¯ is the distance, and λ ¯ is the threshold pressure gradient.
There is a continuity equation for seepage flow in one-dimensional porous media, which is
( ρ v ¯ ) x ¯ = ( ρ ϕ ) t ¯ , 0 x ¯ s ¯ ( t ¯ ) ,
where s ¯ is the moving boundary and t ¯ is time.
If we substitute Equations (1) and (2) into Equation (3), the governing equation is deduced to be
2 p ¯ x ¯ 2 = μ ϕ i C t k p ¯ t ¯ , 0 x ¯ s ¯ ( t ¯ ) ,
where the initial conditions, the internal boundary conditions under variable levels of pressure, and the moving boundary conditions are, respectively,
s ¯ | t ¯ = 0 = 0 , p ¯ | t ¯ = 0 = p ¯ i ,
p ¯ | x ¯ = 0 = f ( t ¯ ) ,
p ¯ | x ¯ = s ¯ ( t ¯ ) = p ¯ i , p ¯ x ¯ | x ¯ = s ¯ ( t ¯ ) = λ ¯ ,
where C t is the comprehensive compression coefficient and f represents the bottomhole pressure function, which changes with time.
The moving boundary condition is the major difference between the non-Darcy seepage model and the classical Darcy seepage model, as shown in Figure 1. Therefore, seepage only occurs in the range where the formation pressure gradients near the well exceed the threshold pressure gradients. There is a smaller pressure gradient outside the moving boundary than there is inside the moving boundary; consequently, seepage cannot occur, and the formation pressure remains the same.
To transform them into a general form, the dimensionless variables are defined as
x = x ¯ x ¯ w , t = k μ ϕ i C t x ¯ w 2 t ¯ , δ = s ¯ x ¯ w , p = k v ¯ w x ¯ w μ ( p ¯ i p ¯ ) , λ = k λ ¯ v ¯ w μ ,
where x ¯ w is the distance constant, x is the non-dimensional distance, t is the non-dimensional time, δ is the non-dimensional moving boundary, p is the non-dimensional formation pressure, and λ is the non-dimensional threshold pressure gradient.
Equations (4)–(7) can be changed into
2 p x 2 = p t , 0 x δ ( t ) ,
δ ( 0 ) = 0 , p | t = 0 = 0 , p | x = 0 = f ( t ) , p | x = δ ( t ) = 0 , p x | x = δ ( t ) = λ ,
Equations (9) and (10) jointly constitute a dimensionless moving boundary for non-Darcy flows in one-dimensional, low-permeability porous media considering a threshold pressure gradient under variable levels of internal boundary pressure. As a result of the moving boundary conditions, seepage will only occur when the formation pressure gradient near the well exceeds the threshold pressure gradient. There will be no seepage beyond the moving boundary, as the formation pressure gradient will be smaller than the threshold pressure gradient. This is the difference between the non-Darcy seepage model considering the influence of the threshold pressure gradient and the classical Darcy seepage model [16].
The function of the variable internal boundary pressure and the similarity variables are as follows:
f ( t ) = 2 U t ,
g = p 2 t , ξ = x 2 t , ε = δ 2 t .
In fact, the production pressure is more difficult to control than the internal boundary during the production process. Interpretation of the well test shows that the slope of the double logarithmic curve of dimensionless pressure with respect to dimensionless time is 1/2. In addition, Equation (11) can also make the mathematical model maintain adequately similar characteristics [16].
For a mathematical model with a non-dimensional moving boundary, Equations (9) and (10) can be transformed into
1 2 2 g ξ 2 + ξ g ξ g = 0 , 0 ξ ε ,
g | ξ = 0 = U , g | ξ = ε = 0 , g ξ | ξ = ε = λ ,
where U is a dimensionless positive fitting parameter of the production pressure data. Equations (13) and (14) set up a closed form of an ordinary differential equation that has fixed boundary conditions, which can be analytically approximated more easily.
Under the conditions of variable internal boundary pressure, 1D non-Darcy flows in low-permeability porous media considering the threshold pressure gradients have the following exact analytical solutions:
p = U [ 2 t e x 2 4 t + π x erf ( x 2 t ) π x erf ( ε ) x e ε 2 ε ] , x [ 0 , δ ] , λ > 0 ,
p = U [ 2 t e x 2 4 t + π x erf ( x 2 t ) π x ] , x [ 0 , + ) , λ = 0 ,
These exact analytical solutions can verify the correction of approximate analytical solutions based on the HAM.
As a result of the HAM approach, a non-linear problem can be reduced to an infinite number of linear problems. Most boundary layer flows decay exponentially after infinity from a physical perspective [45,46,47,48,49,50,51,52,53,54]. According to Equations (13) and (14), g ( ξ ) can be expressed by
{ ε k ξ i exp ( j β ξ ) | k 0 , i 0 , j 0 } ,
g ( ξ , ε ) = k = 0 + i = 0 + j = 0 + a i , j k ε k ξ i exp ( j β ξ ) ,
where a i , j is the constant coefficient to be determined using the HAM.
Using Equations (13) and (14), Equation (17) provides a convenient way to select the initial guess
g 0 ( ξ , ε ) = U ( 1 ξ ε ) e β ξ ,
and the auxiliary linear operators
L g = 2 G ξ 2 β 2 G ,
which have the following properties
L g [ C 0 exp ( β ξ ) + C 1 exp ( β ξ ) ] = 0 ,
where C 0 , C 1 are the integral coefficients.
The HAM deformation equation is constructed as follows
( 1 q ) L g [ G ( ξ , ε ; q ) g 0 ( ξ , ε ) ] = q h g N g [ G ( ξ , ε ; q ) ] = q h [ 1 2 2 G ( ξ , ε ; q ) ξ 2 + ξ G ( ξ , ε ; q ) ξ G ( ξ , ε ; q ) ]
with the boundary conditions
G ( ξ , ε ; q ) | ξ = 0 = U , G ( ξ , ε ; q ) | ξ = ε = 0 , G ( ξ , ε ; q ) ξ | ξ = ε = λ ,
where the embedding parameter is q [ 0 , 1 ] and the non-linear operator based on the governing Equation (13) is N g . Through use of the HAM, approximate analytical solutions can be fully determined:
G ( ξ , ε ; q ) | q = 0 = g ( ξ , ε ) 0 , G ( ξ , ε ; q ) | q = 1 = g ( ξ , ε ) .
The embedded parameter q contributes to mapping. Namely, as q changes from 0 to 1, mapping ensures that G ( ξ , ε ; q ) from the initial guess g 0 ( ξ , ε ) continues to converge to the exact solution g ( ξ , ε ) . According to Taylor’s theorem, G ( ξ , ε ; q ) for power series expansion is
G ( ξ , ε ; q ) = G ( ξ , ε ; 0 ) + n = 1 + g n ( ξ , ε ) q n ,
where
g n ( ξ , ε ) = 1 n ! n G ( ξ , ε ; q ) q n | q = 0 .
As mentioned above, there is considerable freedom of choice regarding the auxiliary linear operator L g , the initial solution g 0 ( ξ , ε ) , and the convergent control parameters. Assuming that all of them are correctly chosen so that the series of Equation (25) converges at q = 1 , the following series solutions can be obtained from Equation (23):
g ( ξ , ε ) = g 0 ( ξ , ε ) + n = 1 + g n ( ξ , ε ) .
Taking the series of Equation (25) and substituting it into the zero-order deformation equation (Equation (22)) and the boundary conditions of Equation (23), the m order deformation equation is obtained by using the equal power coefficients of q
L g [ g ( ξ , ε ) χ m g m 1 ( ξ , ε ) ] = h R m ( ξ , ε ) , m 1 ,
for which the boundary conditions are
g m ( ξ , ε ) | ξ = 0 = 0 , g m ( ξ , ε ) | ξ = ε = 0 , g m ( ξ , ε ) ξ | ξ = ε = 0 ,
where
R m ( ξ , ε ) = 1 2 2 g m 1 ξ 2 + n = 1 m 1 ξ m 1 n g m 1 ξ g m 1 ,
χ m = { 1 m > 1 0 m = 0 .
The right-hand side of Equation (22) is taken from Equation (13), and the left-hand side is independent of ε . Therefore, Equation (22) is actually an ordinary differential equation for ξ , so it is easy to solve. The particular solution of Equation (22) is as follows:
g m * ( ξ , ε ) = χ m g m 1 ( ξ , ε ) + 0 ξ { 0 ξ [ e β ξ R m ( ξ , ε ) e β s d s ] d ξ } d ξ .
According to Equation (21), its general solution is
g m ( ξ , ε ) = g m * ( ξ , ε ) + C 0 , m exp ( β ξ ) + C 1 , m exp ( β ξ ) .
Two integral coefficients are determined by the boundary condition in Equation (29):
C 0 , m = 1 β f m * ξ , C 1 , m = 0 .
Thus, there is an infinite number of linear ordinary differential equations based on Equation (28) with coefficients that are constants derived from the original non-linear partial differential equation (Equation (13)) with variable coefficients.
Obviously, it is much easier to solve a linear ordinary differential equation with constant coefficients than it is to solve a non-linear partial differential equation with variable coefficients. It is worth noting that, unlike perturbation methods, this kind of transformation does not require any small or large physical parameters. The auxiliary linear operator of Equation (20) chosen in this problem had no obvious relationship with the linear term 2 g / ξ 2 in the original Equation (13). This is mainly because, unlike other analytic methods, the homotopy analysis method (HAM) provides a great deal of freedom to choose the auxiliary linear operators appropriately. Without this freedom, there can be no transformation of a non-linear partial differential equation such as Equation (13) with variable coefficients into a linear ordinary differential equation (Equation (28)) with constant coefficients.

3. Results and Discussion

On the basis of the results of the HAM, non-dimensional moving boundary mathematical models of non-Darcy flows in one-dimensional, low-permeability porous media, considering threshold pressure gradients under varying production pressures and internal boundary values, were solved approximately and analytically. It is noteworthy that the series solutions of Equation (27) contained the convergence control parameter h and the auxiliary parameter β . Although there is no physical significance, this can be used to ensure the convergence of the series solution. Obviously, if the series solution of Equation (27) is convergent, then the corresponding first derivative g ξ ( ξ , ε ) should also be convergent. For simplicity, consider the convergence of g ξ ( ξ , ε ) series solutions. First, let h = 1 and consider β to be an unknown variable. As shown in Figure 2a, it was found that for large values of β , for example, when β 2 , the series solutions of g ξ ( ξ , ε ) converged to the same value when ξ = 0 , ε = 1 . Now, g ξ ( ξ , ε ) is only dependent on β and parameter h if ξ , ε are fixed. Given β = 2 , g ξ ( ξ , ε ) is just a power series of h , and its convergence depends on h when ξ = 0 , ε = 1 , as shown in Figure 2b.
Therefore, β = 2 and h = 1 can be simply selected for HAM calculations, as shown in Table 1. In order to verify results of the HAM, as shown in Table 1, comparisons of the results for dimensionless formation pressure p ( x , t ) of 10th-order HAM with exact analytical solutions were made for λ = 0 , U = 0.5 , t = 5000 . As shown in Figure 2, a comparison was made between the results for dimensionless formation pressure p ( x , t ) of 10th-order HAM and the exact analytical solutions of the non-dimensional formation pressure distributions depending on the non-dimensional threshold pressure gradients λ = 0 , 0.2 , 0.5 , 0.8 , 1 when U = 0.5 , t = 5000 . There was no significant difference between the 10th-order HAM results and the exact analytical solutions, with relative errors between 0.0089% and 2.6023% when the dimensionless distance was taken into account. The results of the 10th-order HAM matched the exact analytical solutions well, as shown in Table 1 and Figure 2. Therefore, the correctness and effectiveness of solving a moving boundary model by using the homotopy analysis method (HAM) can be verified.
As can be seen in Figure 3, the threshold pressure gradients have a profound influence on one-dimensional porous media seepage models. When the non-dimensional threshold pressure gradient was λ = 0 , the pressure drop of the Darcy seepage model can be instantaneously transmitted to infinity. When λ > 0 , the formation pressure corresponding to the non-zero non-dimensional threshold pressure gradient value was greater than zero within the moving boundary, while the formation pressure outside the moving boundary remained at zero. The pressure drop curve of the non-Darcy seepage model showed the characteristics of tight support, which clearly showed a different distribution of the formation pressures in contrast to the Darcy seepage model. Thus, the greater the non-dimensional threshold pressure gradients, the greater the difference between the solutions of Darcy’s seepage model and the moving boundary model.
Figure 4a illustrates distributions of the dimensionless transient distance δ ( t ) compared with the transient distance with moving boundaries depending on different values of non-dimensional threshold pressure gradients ( λ = 0.2 , 0.5 , 0.8 , 1 ) when U = 0.5 . As λ rises from 0.2 to 1, the dimensionless transient distance δ ( t ) increases. The value of the dimensionless transient distance δ ( t ) increases with an increase in the non-dimensional time t . When the value of the dimensionless time t is small, the increase in the dimensionless transient distance δ ( t ) is not obvious; however, when the value of the dimensionless time t is large, the value of the dimensionless transient distance δ ( t ) will increase significantly. As shown in Figure 4b, the distribution of the dimensionless transient distance δ ( t ) was compared with the transient distance of a moving boundary, depending on different values of non-dimensional positive fitting parameters within the production pressure data ( U = 0.2 , 0.5 , 0.8 , 1 ) when λ = 0.5 . As U decreased from 1 to 0.2, the dimensionless transient distance δ ( t ) increased. As the dimensionless time t increased, the value of the dimensionless transient distance δ ( t ) increased. When the value of the dimensionless time t was small, the increased in the value of the dimensionless transient distance δ ( t ) was not obvious; however, when the value of the dimensionless time t was large, the value of the dimensionless transient distance δ ( t ) increased significantly.
As shown in Figure 5a, distributions of non-dimensional formation pressure p ( x , t ) were compared with those of non-dimensional formation pressure depending on different values of the non-dimensional positive fitting parameter within the production pressure data ( U = 0.2 , 0.5 , 0.8 , 1 ) when λ = 0 , t = 5000 . As U increased from 0.2 to 1, the dimensionless formation pressure p ( x , t ) increased. As the dimensionless distance x increased, the value of the dimensionless formation pressure p ( x , t ) decreased. When the value of the dimensionless distance x was small, the increase in the dimensionless formation pressure p ( x , t ) was significant as U increased from 0.2 to 1; however, when the value of the dimensionless distance x was large, the increase in the dimensionless formation pressure p ( x , t ) was not obvious as U increased from 0.2 to 1. When λ = 0 , the dimensionless formation pressure p ( x , t ) tended to zero when the dimensionless distance x was limited to infinity. Namely, the pressure drop in the Darcy seepage model could be instantaneously transmitted to infinity. As shown in Figure 5b, the distributions of dimensionless formation pressure p ( x , t ) were compared with the distributions of non-dimensional formation pressure depending on different values of the non-dimensional positive fitting parameter within the production pressure data ( U = 0.2 , 0.5 , 0.8 , 1 ) when λ = 0.5 , t = 5000 . As U increased from 0.2 to 1, the dimensionless formation pressure p ( x , t ) increased. As the dimensionless distance x increased, the value of the dimensionless formation pressure p ( x , t ) decreased. When the non-dimensional distance x was small, the increase in the non-dimensional formation pressure p ( x , t ) was significant as U increased from 0.2 to 1; however, when the non-dimensional distance x was large, the increase in the non-dimensional formation pressure p ( x , t ) was not obvious as U increased from 0.2 to 1. When the non-dimensional threshold pressure gradient was λ > 0 , the pressure drop curves of non-Darcy models showed the characteristics of tight support, which are different from the formation pressure distribution curves corresponding to the Darcy seepage models. Whether the moving boundary conditions are considered makes a great deal of difference to the model’s calculation of the results. Therefore, the influence of moving boundaries should be considered when studying unsteady flows in low-permeability porous media, especially when threshold pressure gradients in tight porous media (such as tight oil, shale oil, etc.) are large.

4. Conclusions

The homotopy analysis method (HAM) was used to solve the approximate analytical solutions of non-Darcy flows with moving boundary conditions of a low-permeability reservoir, which were converted into a fixed boundary mathematical model via similarity transformation. Compared with the exact analytical solutions, the accuracy of the analytical approximate solutions based on the HAM is guaranteed. The numerical stability was high, and the relative errors of the calculated case were between 0.0089% and 2.6023%. When λ = 0 , the pressure drop of Darcy seepage models can be instantaneously transmitted to infinity. When λ > 0 , the pressure drop curves of the non-Darcy seepage model showed the characteristics of tight support, which made them obviously different from the distribution curves of the formation pressure corresponding to the Darcy seepage models. Therefore, the influence of a moving boundary should be considered when studying the unsteady flows in low-permeability porous media, especially when the threshold pressure gradients in tight porous media (such as tight oil, shale oil, etc.) are large. The results of HAM showed that the greater the non-dimensional time, the greater the influence of the threshold pressure gradients on the moving boundary. The smaller the threshold pressure gradients, the faster the moving boundaries will move and the stronger the sensitivity affecting the moving boundaries will be. When the threshold pressure gradient data tend to zero, one-dimensional, low-permeability, porous media non-Darcy flows with moving boundary conditions can be reducible to Darcy seepage models. In future research, it is suggested that the present research should be compared with real cases to determine the non-dimensional positive fitting parameter of the production pressure data U and the non-dimensional threshold pressure gradient λ .

Author Contributions

Writing—original draft preparation, X.Y.; project administration, S.L., L.K. and L.C. All authors have read and agreed to the published version of the manuscript.

Funding

National Natural Science Foundation of China (12002390, 52174011).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The article includes all relevant data.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

p ¯ Formation pressure
p ¯ i Original formation pressure
C f Fluid compression coefficient
C ϕ Pore compression coefficient
v ¯ Seepage velocity
k Permeability
x ¯ Distance
t ¯ Time
s ¯ Moving boundary
C t Comprehensive compression coefficient
f Bottom-hole pressure function that changes with time
x ¯ w Distance constant
x Non-dimensional distance
t Non-dimensional time
δ Non-dimensional moving boundary
p Non-dimensional formation pressure
U Non-dimensional positive fitting parameter of production pressure data
Greek Symbols
ρ Fluid density
ρ i Initial density of fluid
ϕ Porosity
ϕ i Initial porosity
μ Fluid viscosity
λ ¯ Threshold pressure gradient
λ Non-dimensional threshold pressure gradient

References

  1. Cai, J.C. A fractal approach to low velocity non-Darcy flow in a low permeability porous medium. Chin. Phys. B 2014, 23, 044701. [Google Scholar] [CrossRef]
  2. Wang, F.Y.; Liu, Z.C.; Cai, J.C.; Gao, J. A fractal model for low-velocity non-Darcy flow in tight oil reservoirs considering boundary-layer effect. Fractals 2018, 26, 1850077. [Google Scholar] [CrossRef]
  3. Ye, W.; Wang, X.; Cao, C.; Yu, W. A fractal model for threshold pressure gradient of tight oil reservoirs. J. Petrol. Sci. Eng. 2019, 179, 427–431. [Google Scholar] [CrossRef]
  4. Hayat, T.; Yasmin, H.; Ahmad, B.; Chen, B. Simultaneous effects of convective conditions and nanoparticles on peristaltic motion. J. Mol. Liq. 2014, 193, 74–82. [Google Scholar] [CrossRef]
  5. Yasmin, H.; Iqbal, N. Convective Mass/Heat Analysis of an Electroosmotic Peristaltic Flow of Ionic Liquid in a Symmetric Porous Microchannel with Soret and Dufour. Math. Probl. Eng. 2021, 2021, 2638647. [Google Scholar] [CrossRef]
  6. Prada, A.; Civan, F. Modification of Darcy’s law for the threshold pressure gradient. J. Petrol. Sci. Eng. 1999, 22, 237–240. [Google Scholar] [CrossRef]
  7. Tian, W.; Li, A.; Ren, X.; Josephine, Y. The threshold pressure gradient effect in the tight sandstone gas reservoirs with high water saturation. Fuel 2018, 226, 221–229. [Google Scholar] [CrossRef]
  8. Xiong, W.; Lei, Q.; Gao, S.; Hu, Z.; Xue, H. Pseudo threshold pressure gradient to flow for low permeability reservoirs. Petrol. Explor. Dev. 2009, 36, 232–236. [Google Scholar]
  9. Zhu, W.; Song, H.; Huang, X.; Liu, X. Pressure characteristics and effective deployment in a water-bearing tight gas reservoir with low velocity non-Darcy flow. Energy Fuel 2011, 25, 1111–1117. [Google Scholar] [CrossRef]
  10. Yao, T.; Huang, Y.; Li, J. Nonliner flow equations for heavy oil in porous media. Chin. J. Theor. Appl. Mech. 2012, 44, 106–110. [Google Scholar]
  11. Zhao, L.; Jiang, H.; Wang, H.; Yang, H.; Sun, F.; Li, J. Representation of a new physics-based non-Darcy equation for low-velocity flow in tight reservoirs. J. Petrol. Sci. Eng. 2020, 184, 106518. [Google Scholar] [CrossRef]
  12. Wei, Q.; Zhou, H.; Yang, S. Non-Darcy flow models in porous media via Atangana-Baleanu derivative. Chaos Solitons Fractals 2020, 141, 110335. [Google Scholar] [CrossRef]
  13. Zhou, Y.; Zhang, L.; Wang, T. Analytical solution for one-dimensional non-Darcy flow with bilinear relation in porous medium caused by line source. Appl. Math. Comput. 2021, 392, 125674. [Google Scholar] [CrossRef]
  14. Nie, R.; Fan, X.; Li, M.; Chen, Z.; Deng, Q.; Lu, C.; Zhou, Z.; Jiang, D.; Zhan, J. Modeling transient flow behavior with the high velocity non-Darcy effect in composite naturally fractured-homogeneous gas reservoirs. J. Nat. Gas Sci. Eng. 2021, 96, 104269. [Google Scholar] [CrossRef]
  15. Moghimi, H.; Siavashi, M.; Nezhad, M.M.; Guadagnini, A. Pore-scale computational analyses of non-Darcy flow through highly porous structures with various degrees of geometrical complexity. Energy Technol. Assess. 2022, 52, 102048. [Google Scholar] [CrossRef]
  16. Zhang, X.; Dou, Z.; Wang, J.; Zhou, Z.; Zhuang, C. Non-Darcy flows in layered porous media (LPMs) with contrasting pore space structures. Pet. Sci. 2022, 19, 2004–2013. [Google Scholar] [CrossRef]
  17. Cheng, S.; Zhang, S.; Huang, Y.; Zhu, W. An integral solution of free-boundary problem of non-Darcy flow behavior. Mech. Eng. 2002, 24, 15–17. [Google Scholar]
  18. Liu, W.; Yao, J.; Chen, Z.; Liu, Y.; Sun, H. Research on analytical and numerical solution of a moving boundary model of seepage flow in low-permeable porous media. Chin. J. Theor. Appl. Mech. 2015, 47, 605–612. [Google Scholar]
  19. Liu, W.; Yao, J.; Wang, Y. Exact analytical solutions of moving boundary problems of one-dimensional flow in semi-infinite long porous media with threshold pressure gradient. Int. J. Heat Mass Transf. 2012, 55, 6017–6022. [Google Scholar] [CrossRef]
  20. Liu, W. Analytical study on a moving boundary problem of semispherical centripetal seepage flow of Bingham fluid with threshold pressure gradient. Int. J. Non-Linear Mech. 2019, 113, 17–30. [Google Scholar] [CrossRef]
  21. Liu, W. Exact analytical solution of a generalized multiple moving boundary model of one-dimensional non-Darcy flow in heterogeneous multilayered low-permeability porous media with a threshold pressure gradient. Appl. Math. Model. 2020, 81, 931–953. [Google Scholar] [CrossRef]
  22. Guo, Y.; Zeng, Q.; Wang, Z. Numerical simulation of low permeability flow with moving-boundary using meshless methods. Eng. Mech. 2006, 23, 188–192. [Google Scholar]
  23. Chen, L.; Ren, S.; Lian, P. Well test analysis on low velocity and non-Darcy flow in dual-porosity reservoir with dynamic boundary. Chin. J. Comput. Mech. 2011, 28, 879–883. [Google Scholar]
  24. Yao, J.; Liu, W. New method for solution of the model of non-Darcy seepage flow in low-permeability reservoirs with moving boundary. Chin. Q. Mech. 2012, 33, 597–601. [Google Scholar]
  25. Yao, J.; Liu, W.; Chen, Z. Numerical solution of a moving boundary problem of one-dimensional flow in semi-Infinite long porous media with threshold pressure gradient. Math. Probl. Eng. 2013, 2013, 384246. [Google Scholar] [CrossRef] [Green Version]
  26. Liu, W.; Zhang, Q.; Zhu, W. Numerical simulation of multi-stage fractured horizontal well in low-permeable oil reservoir with threshold pressure gradient with moving boundary. J. Petrol. Sci. Eng. 2019, 178, 1112–1127. [Google Scholar] [CrossRef]
  27. Xu, H.; Liu, N.; Chen, Y.; Tian, Y.; Guo, Z.; Jiang, W.; He, Y. A novel equivalent numerical simulation method for non-Darcy seepage flow in low-permeability reservoirs. Energies 2022, 15, 8505. [Google Scholar] [CrossRef]
  28. Liao, S. The Proposed Homotopy Analysis Technique for the Solution of Nonlinear Problems. Ph.D. Thesis, Shanghai Jiao Tong University, Shanghai, China, 1992. [Google Scholar]
  29. Liao, S. Beyond Perturbation: Introduction to the Homotopy Analysis Method, 1st ed.; Chapman and Hall/CRC: New York, NY, USA, 2003; pp. 154–196. [Google Scholar]
  30. Liao, S. Homotopy Analysis Method in Nonlinear Differential Equations, 1st ed.; Springer: Berlin/Heidelberg, Germany, 2012; pp. 40–86. [Google Scholar]
  31. Nassar, C.J.; Revelli, J.F.; Bowman, R.J. Application of the homotopy analysis method to the Poisson-Boltzmann equation for semiconductor devices. Commun. Nonlinear Sci. Numer. Simul. 2011, 16, 2501–2512. [Google Scholar] [CrossRef]
  32. Vajravelu, K.; Van Gorder, R.A. Nonlinear Flow Phenomena and Homotopy Analysis: Fluid Flow and Heat Transfer; Springer: Berlin/Heidelberg, Germany, 2012. [Google Scholar]
  33. Sardanyés, J.; Rodrigues, C.; Januário, C.; Martins, N.; Gil-Gómez, G.; Duarte, J. Activation of effector immune cells promotes tumor stochastic extinction: A homotopy analysis approach. Appl. Math. Comput. 2015, 252, 484–495. [Google Scholar] [CrossRef] [Green Version]
  34. Ramzan, M.; Bilal, M.; Chung, J.D.; Lu, D.C.; Farooq, U. Impact of generalized Fourier’s and Fick’s laws on MHD 3D second grade nanofluid flow with variable thermalconductivity and convective heat and mass conditions. Phys. Fluids 2017, 29, 093102. [Google Scholar] [CrossRef]
  35. Patel, M.A.; Desai, N.B. Homotopy analysis approach of Boussinesq equation for infiltration phenomenon in unsaturated porous medium. Math. J. Interdiscip. Sci. 2018, 7, 21–28. [Google Scholar] [CrossRef]
  36. Liu, L.; Rana, J.; Liao, S. Analytical solutions for the hydrogen atom in plasmas with electric, magnetic, and Aharonov-Bohm flux fields. Phys. Rev. E 2021, 103, 023206. [Google Scholar] [CrossRef] [PubMed]
  37. Yang, X.Y.; Li, Y. On bi-chromatic steady-state gravity waves with an arbitrary included angle. Phys. Fluids 2022, 34, 032107. [Google Scholar] [CrossRef]
  38. Liao, S. A new non-perturbative approach in quantum mechanics for time-independent Schrödinger equations. Sci. China Phys. Mech. 2020, 63, 234612. [Google Scholar] [CrossRef] [Green Version]
  39. Liao, S. Avoiding small denominator problems by means of the Homotopy Analysis Method. Adv. Appl. Math. Mech. 2023, 15, 267–299. [Google Scholar] [CrossRef]
  40. Hayat, T.; Alsaedi, A.; Ahmad, B. Thermo diffusion and diffusion thermo impacts on bioconvection Walter-B nanomaterial involving gyrotactic microorganisms. Alex. Eng. J. 2021, 60, 5537–5545. [Google Scholar] [CrossRef]
  41. Al-Qudaha, A.; Odibatb, Z.; Shawagfeha, N. A linearization-based computational algorithm of homotopy analysis method for nonlinear reaction-diffusion systems. Math. Comput. Simul. 2022, 194, 505–522. [Google Scholar] [CrossRef]
  42. Bottona, E.; Greenberga, J.B.; Arad, A.; Katoshevski, D.; Vaikuntanathan, V.; Ibach, M.; Weigand, B. An investigation of grouping of two falling dissimilar droplets using the homotopy analysis method. Appl. Math. Model 2022, 104, 486–498. [Google Scholar] [CrossRef]
  43. Mehmet, P.; Muhammet, Y. Similarity transformations for partial differential equations. SIAM Rev. 1998, 40, 96–101. [Google Scholar]
  44. Voller, V.R.; Swenson, J.B.; Paola, C. An analytical solution for a Stefan problem with variable latent heat. Int. J. Heat Mass Tran. 2004, 47, 5387–5390. [Google Scholar] [CrossRef]
  45. Liao, S. An explicit, totally analytic approximate solution for Blasius’ viscous flow problems. Int. J. Non Linear Mech. 1999, 34, 759–778. [Google Scholar] [CrossRef]
  46. Liao, S. On the analytical solution of magnetohydrodynamic flows of non-Newtonian fluids over a stretching sheet. J. Fluid Mech. 2003, 488, 189–212. [Google Scholar] [CrossRef] [Green Version]
  47. Xu, H.; Pop, I. Homotopy analysis of unsteady boundary-layer flow started impulsively from rest along a symmetric wedge. Z. Angew. Math. Mech. 2008, 88, 507–514. [Google Scholar] [CrossRef]
  48. You, X.; Xu, H.; Liao, S. On the non-similarity boundary-layer flows of second-order fluid over a stretching sheet. J. Appl. Mech. 2010, 77, 021002. [Google Scholar] [CrossRef]
  49. You, X.; Xu, H.; Pop, I. Free convection along a convectively heated vertical flat sheet embedded in a saturated porous medium. Int. Commun. Heat Mass 2014, 55, 102–108. [Google Scholar] [CrossRef]
  50. Farooq, U.; Zhao, Y.L.; Hayat, T.; Alsaedi, A.; Liao, S. Application of the HAM-based Mathematica package BVPh 2.0 on MHD Falkner–Skan flow of nano-fluid. Comput. Fluids 2015, 111, 69–75. [Google Scholar] [CrossRef]
  51. Mustafa, M.; Hayat, T.; Pop, I.; Hendi, A. Stagnation-point flow and heat transfer of a Casson fluid towards a stretching sheet. Z. Nat. A 2015, 67, 70–76. [Google Scholar] [CrossRef]
  52. Ali, A.; Zaman, H.; Abidin, M.Z.; Shah, S.I.A. Analytic solution for fluid flow over an exponentially stretching porous sheet with surface heat flux in porous medium by means of Homotopy Analysis Method. Am. J. Comput. Math. 2015, 5, 224–238. [Google Scholar] [CrossRef] [Green Version]
  53. Liu, L.; Li, J.; Liao, S. Explicit solutions of MHD flow and heat transfer of Casson fluid over an exponentially shrinking sheet with suction. Nanomaterials 2022, 12, 3289. [Google Scholar] [CrossRef]
  54. Yang, Y.; Liao, S. Comparison between homotopy analysis method and homotopy renormalization method in fluid mechanics. Eur. J. Mech. B Fluids 2023, 97, 187–198. [Google Scholar] [CrossRef]
Figure 1. Schematic description of an oil reservoir.
Figure 1. Schematic description of an oil reservoir.
Energies 16 02175 g001
Figure 2. Curves of g ξ ( ξ , ε ) against β and h when ξ = 0 , ε = 1 : (a) 10th-order and 15th-order HAM approximations when h = 1 ; (b) 10th-order and 15th-order HAM approximations when β = 1 .
Figure 2. Curves of g ξ ( ξ , ε ) against β and h when ξ = 0 , ε = 1 : (a) 10th-order and 15th-order HAM approximations when h = 1 ; (b) 10th-order and 15th-order HAM approximations when β = 1 .
Energies 16 02175 g002
Figure 3. Comparisons between the results of dimensionless formation pressure p ( x , t ) of 10th-order HAM and the exact analytical solutions of the non-dimensional formation pressure distribution depending on λ when U = 0.5 , t = 5000 .
Figure 3. Comparisons between the results of dimensionless formation pressure p ( x , t ) of 10th-order HAM and the exact analytical solutions of the non-dimensional formation pressure distribution depending on λ when U = 0.5 , t = 5000 .
Energies 16 02175 g003
Figure 4. The distributions of the non-dimensional transient distance δ ( t ) of moving boundaries: (a) comparisons of the transient distance of moving boundaries depending on various values of λ when U = 0.5 ; (b) comparisons of the transient distance of moving boundaries depending on various values of U when λ = 0.5 .
Figure 4. The distributions of the non-dimensional transient distance δ ( t ) of moving boundaries: (a) comparisons of the transient distance of moving boundaries depending on various values of λ when U = 0.5 ; (b) comparisons of the transient distance of moving boundaries depending on various values of U when λ = 0.5 .
Energies 16 02175 g004
Figure 5. The distribution of non-dimensional formation pressure p ( x , t ) : (a) Comparison of the distributions of non-dimensional formation pressure depending on different values of U when λ = 0 , t = 5000 ; (b) Comparison of the distributions of non-dimensional formation pressure depending on different values of U when λ = 0.5 , t = 5000 .
Figure 5. The distribution of non-dimensional formation pressure p ( x , t ) : (a) Comparison of the distributions of non-dimensional formation pressure depending on different values of U when λ = 0 , t = 5000 ; (b) Comparison of the distributions of non-dimensional formation pressure depending on different values of U when λ = 0.5 , t = 5000 .
Energies 16 02175 g005
Table 1. Comparison of the results of dimensionless formation pressure p ( x , t ) of 10th-orderHAM and the exact analytical solutions when λ = 0 , U = 0.5 , t = 5000 .
Table 1. Comparison of the results of dimensionless formation pressure p ( x , t ) of 10th-orderHAM and the exact analytical solutions when λ = 0 , U = 0.5 , t = 5000 .
x p ( x , 5000 )   ( Analytical ) p ( x , 5000 ) (10th-Order HAM)Relative Error (%)
070.717070.71070.0089
566.367966.35890.0135
1062.201762.18480.0271
1558.211358.18790.0402
2054.395754.36720.0524
3047.282247.24820.0719
5035.058535.03040.0802
10014.767314.78950.1501
1505.19445.20460.1948
2001.50491.47042.3490
2500.35520.34632.5577
3000.06770.06602.6023
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

You, X.; Li, S.; Kang, L.; Cheng, L. A Study of the Non-Linear Seepage Problem in Porous Media via the Homotopy Analysis Method. Energies 2023, 16, 2175. https://doi.org/10.3390/en16052175

AMA Style

You X, Li S, Kang L, Cheng L. A Study of the Non-Linear Seepage Problem in Porous Media via the Homotopy Analysis Method. Energies. 2023; 16(5):2175. https://doi.org/10.3390/en16052175

Chicago/Turabian Style

You, Xiangcheng, Shiyuan Li, Lei Kang, and Li Cheng. 2023. "A Study of the Non-Linear Seepage Problem in Porous Media via the Homotopy Analysis Method" Energies 16, no. 5: 2175. https://doi.org/10.3390/en16052175

APA Style

You, X., Li, S., Kang, L., & Cheng, L. (2023). A Study of the Non-Linear Seepage Problem in Porous Media via the Homotopy Analysis Method. Energies, 16(5), 2175. https://doi.org/10.3390/en16052175

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop