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

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

(1)$\rho ={\rho }_i{e}^{-C_f\left({\overline{p}}_i-\overline{p}\right)},\varphi ={\varphi }_i{e}^{-C_{\varphi }\left({\overline{p}}_i-\overline{p}\right)},$

where $\rho$ is the density of the fluid, ${\rho }_i$ is the initial density of the fluid, $\overline{p}$ is the formation pressure, ${\overline{p}}_i$ is the original formation pressure, $\varphi$ is the porosity, ${\varphi }_i$ is the initial porosity, $C_f$ is the fluid compression coefficient, and $C_{\varphi }$ 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

(2)$\overline{v}=\{\begin{array}{l}\begin{array}{cccc}0,& & & \begin{array}{ccc}& & \begin{array}{cc}& 0\le \left|\frac{\partial \overline{p}}{\partial \overline{x}}\right|\le \overline{\lambda }\end{array}\end{array}\end{array}\\ -\frac{k}{\mu }\frac{\partial \overline{p}}{\partial \overline{x}}\left(1-\frac{\overline{\lambda }}{\left|\frac{\partial \overline{p}}{\partial \overline{x}}\right|}\right)\begin{array}{cc},& \left|\frac{\partial \overline{p}}{\partial \overline{x}}\right|>\overline{\lambda }\end{array}\end{array},$

where $\overline{v}$ is the seepage velocity, $k$ is the permeability, $\mu$ is the fluid viscosity, $\overline{x}$ is the distance, and $\overline{\lambda }$ is the threshold pressure gradient.

There is a continuity equation for seepage flow in one-dimensional porous media, which is

(3)$-\frac{\partial \left(\rho \overline{v}\right)}{\partial \overline{x}}=\frac{\partial \left(\rho \varphi \right)}{\partial \overline{t}}\begin{array}{cc},& 0\le \overline{x}\le \overline{s}\left(\overline{t}\right)\end{array},$

where $\overline{s}$ is the moving boundary and $\overline{t}$ is time.

If we substitute Equations (1) and (2) into Equation (3), the governing equation is deduced to be

(4)$\frac{{\partial }^2\overline{p}}{\partial {\overline{x}}^2}=\frac{\mu {\varphi }_i{C}_t}{k}\frac{\partial \overline{p}}{\partial \overline{t}}\begin{array}{cc},& 0\le \overline{x}\le \overline{s}\left(\overline{t}\right)\end{array},$

where the initial conditions, the internal boundary conditions under variable levels of pressure, and the moving boundary conditions are, respectively,

(5)${\overline{s}|}_{\overline{t}=0}=0{\begin{array}{cc},& \overline{p}|\end{array}}_{\overline{t}=0}={\overline{p}}_i,$

(6)${\overline{p}|}_{\overline{x}=0}=f\left(\overline{t}\right),$

(7)${\overline{p}|}_{\overline{x}=\overline{s}\left(\overline{t}\right)}={\overline{p}}_i{\begin{array}{cc},& \frac{\partial \overline{p}}{\partial \overline{x}}|\end{array}}_{\overline{x}=\overline{s}(\overline{t})}=\overline{\lambda },$

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

(8)$x=\frac{\overline{x}}{{\overline{x}}_w},t=\frac{k}{\mu {\varphi }_i{C}_t{\overline{x}}_w^2}\overline{t},\delta =\frac{\overline{s}}{{\overline{x}}_w},p=\frac{k}{{\overline{v}}_w{\overline{x}}_w\mu }\left({\overline{p}}_i-\overline{p}\right),\lambda =\frac{k\overline{\lambda }}{{\overline{v}}_w\mu },$

where ${\overline{x}}_w$ is the distance constant, $x$ is the non-dimensional distance, $t$ is the non-dimensional time, $\delta$ is the non-dimensional moving boundary, $p$ is the non-dimensional formation pressure, and $\lambda$ is the non-dimensional threshold pressure gradient.

Equations (4)–(7) can be changed into

(9)$\frac{{\partial }^{2}p}{\partial x^2}=\frac{\partial p}{\partial t}\begin{array}{cc},& 0\le x\le \delta \left(t\right)\end{array},$

(10)$\delta \left(0\right)=0,{p|}_{t=0}=0,{p|}_{x=0}=f\left(t\right),{p|}_{x=\delta \left(t\right)}=0,{\frac{\partial p}{\partial x}|}_{x=\delta \left(t\right)}=-\lambda ,$

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:

(11)$f\left(t\right)=2U\sqrt{t},$

(12)$g=\frac{p}{2\sqrt{t}},\xi =\frac{x}{2\sqrt{t}},\epsilon =\frac{\delta }{2\sqrt{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

(13)$\frac{1}{2}\frac{{\partial }^{2}g}{\partial {\xi }^2}+\xi \frac{\partial g}{\partial \xi }-g=0\begin{array}{cc},& 0\le \xi \le \epsilon \end{array},$

(14)${g|}_{\xi =0}=U,{g|}_{\xi =\epsilon }=0,{\frac{\partial g}{\partial \xi }|}_{\xi =\epsilon }=-\lambda ,$

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:

(15)$p=U\left[2\sqrt{t}{e}^{-\frac{x^2}{4t}}+\sqrt{\pi }x\mathrm{erf}\left(\frac{x}{2\sqrt{t}}\right)-\sqrt{\pi }x\mathrm{erf}\left(\epsilon \right)-\frac{x}{e^{{\epsilon }^2}\epsilon }\right]\begin{array}{cc},& x\end{array}\in [0,\delta ]\begin{array}{cc},& \lambda >0\end{array},$

(16)$p=U\left[2\sqrt{t}{e}^{-\frac{x^2}{4t}}+\sqrt{\pi }x\mathrm{erf}\left(\frac{x}{2\sqrt{t}}\right)-\sqrt{\pi }x\right]\begin{array}{cc},& x\in [0,+\infty )\end{array}\begin{array}{cc},& \lambda =0\end{array},$

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\left(\xi \right)$ can be expressed by

(17)$\left\{{\epsilon }^k{\xi }^i\exp \left(-j\beta \xi \right)|k\ge 0,i\ge 0,j\ge 0\right\},$

(18)$g(\xi ,\epsilon )={\displaystyle \sum _{k=0}^{+\infty }{\displaystyle \sum _{i=0}^{+\infty }{\displaystyle \sum _{j=0}^{+\infty }{a}_{i,j}^k{\epsilon }^k{\xi }^i\exp \left(-j\beta \xi \right)}}},$

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

(19)$g_0(\xi ,\epsilon )=U\left(1-\frac{\xi }{\epsilon }\right)e^{-\beta \xi },$

and the auxiliary linear operators

(20)$L_g=\frac{{\partial }^{2}G}{\partial {\xi }^2}-{\beta }^{2}G,$

which have the following properties

(21)$L_g\left[C_0\exp \left(-\beta \xi \right)+C_1\exp \left(\beta \xi \right)\right]=0,$

where $C_0$, $C_1$ are the integral coefficients.

The HAM deformation equation is constructed as follows

(22)$\begin{array}{l}\left(1-q\right)L_g\left[G\left(\xi ,\epsilon ;q\right)-g_0\left(\xi ,\epsilon \right)\right]=q{h}_g{N}_g\left[G\left(\xi ,\epsilon ;q\right)\right]\\ \begin{array}{cccc}& & & \end{array}\begin{array}{cccc}& & & \end{array}\begin{array}{ccc}& & =\end{array}qh\left[\frac{1}{2}\frac{{\partial }^{2}G\left(\xi ,\epsilon ;q\right)}{\partial {\xi }^2}+\xi \frac{\partial G\left(\xi ,\epsilon ;q\right)}{\partial \xi }-G\left(\xi ,\epsilon ;q\right)\right]\end{array}$

with the boundary conditions

(23)${G\left(\xi ,\epsilon ;q\right)|}_{\xi =0}={U\begin{array}{cc},& {G\left(\xi ,\epsilon ;q\right)|}_{\xi =\epsilon }=\end{array}\begin{array}{cc}0,& \frac{\partial G\left(\xi ,\epsilon ;q\right)}{\partial \xi }\end{array}|}_{\xi =\epsilon }=-\lambda ,$

where the embedding parameter is $q\in [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:

(24)${G\left(\xi ,\epsilon ;q\right)|}_{q=0}=g{}_0\left(\xi ,\epsilon \right)\begin{array}{cc},& {G\left(\xi ,\epsilon ;q\right)|}_{q=1}=g\left(\xi ,\epsilon \right)\end{array}.$

The embedded parameter $q$ contributes to mapping. Namely, as $q$ changes from 0 to 1, mapping ensures that $G\left(\xi ,\epsilon ;q\right)$ from the initial guess $g_0\left(\xi ,\epsilon \right)$ continues to converge to the exact solution $g\left(\xi ,\epsilon \right)$. According to Taylor’s theorem, $G\left(\xi ,\epsilon ;q\right)$ for power series expansion is

(25)$G\left(\xi ,\epsilon ;q\right)=G\left(\xi ,\epsilon ;0\right)+{\displaystyle \sum _{n=1}^{+\infty }{g}_n\left(\xi ,\epsilon \right)}{q}^n,$

where

(26)$g_n\left(\xi ,\epsilon \right)={\frac{1}{n!}\frac{{\partial }^{n}G\left(\xi ,\epsilon ;q\right)}{\partial 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\left(\xi ,\epsilon \right)$, 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):

(27)$g\left(\xi ,\epsilon \right)=g_0\left(\xi ,\epsilon \right)+{\displaystyle \sum _{n=1}^{+\infty }{g}_n\left(\xi ,\epsilon \right)}.$

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$

(28)$L_g\left[g\left(\xi ,\epsilon \right)-{\chi }_m{g}_{m-1}\left(\xi ,\epsilon \right)\right]=h{R}_m\left(\xi ,\epsilon \right)\begin{array}{cc},& m\ge 1\end{array},$

for which the boundary conditions are

(29)${g_m\left(\xi ,\epsilon \right)|}_{\xi =0}={0\begin{array}{cc},& {g_m\left(\xi ,\epsilon \right)|}_{\xi =\epsilon }=\end{array}\begin{array}{cc}0,& \frac{\partial g_m\left(\xi ,\epsilon \right)}{\partial \xi }\end{array}|}_{\xi =\epsilon }=0,$

where

(30)$R_m\left(\xi ,\epsilon \right)=\frac{1}{2}\frac{{\partial }^2{g}_{m-1}}{\partial {\xi }^2}+{\displaystyle \sum _{n=1}^{m-1}{\xi }_{m-1-n}}\frac{\partial g_{m-1}}{\partial \xi }-g_{m-1},$

(31)${\chi }_m=\{\begin{array}{l}\begin{array}{cc}1& m>1\end{array}\\ \begin{array}{cc}0& m=0\end{array}\end{array}.$

The right-hand side of Equation (22) is taken from Equation (13), and the left-hand side is independent of $\epsilon$. Therefore, Equation (22) is actually an ordinary differential equation for $\xi$, so it is easy to solve. The particular solution of Equation (22) is as follows:

(32)$g_m^{*}\left(\xi ,\epsilon \right)={\chi }_m{g}_{m-1}\left(\xi ,\epsilon \right)+\hslash {\displaystyle {\int }_0^{\xi }\left\{{\displaystyle {\int }_0^{\xi }\left[e^{-\beta \xi }{R}_m\left(\xi ,\epsilon \right)e^{\beta s}ds\right]d\xi }\right\}}d\xi .$

According to Equation (21), its general solution is

(33)$g_m\left(\xi ,\epsilon \right)=g_m^{*}\left(\xi ,\epsilon \right)+C_{0,m}\exp \left(-\beta \xi \right)+C_{1,m}\exp \left(\beta \xi \right).$

Two integral coefficients are determined by the boundary condition in Equation (29):

(34)$C_{0,m}=\frac{1}{\beta }\frac{\partial f_m^{*}}{\partial \xi }\begin{array}{cc},& C_{1,m}=0\end{array}.$

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 ${\partial }^{2}g/\partial {\xi }^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 $\beta$. 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_{\xi }\left(\xi ,\epsilon \right)$ should also be convergent. For simplicity, consider the convergence of $g_{\xi }\left(\xi ,\epsilon \right)$ series solutions. First, let $h=-1$ and consider $\beta$ to be an unknown variable. As shown in Figure 2a, it was found that for large values of $\beta$, for example, when $\beta \ge 2$, the series solutions of $g_{\xi }\left(\xi ,\epsilon \right)$ converged to the same value when $\xi =0,\epsilon =1$. Now, $g_{\xi }\left(\xi ,\epsilon \right)$ is only dependent on $\beta$ and parameter $h$ if $\xi ,\epsilon$ are fixed. Given $\beta =2$, $g_{\xi }\left(\xi ,\epsilon \right)$ is just a power series of $h$, and its convergence depends on $h$ when $\xi =0,\epsilon =1$, as shown in Figure 2b.

Therefore, $\beta =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\left(x,t\right)$ of 10th-order HAM with exact analytical solutions were made for $\lambda =0,U=0.5,t=5000$. As shown in Figure 2, a comparison was made between the results for dimensionless formation pressure $p\left(x,t\right)$ of 10th-order HAM and the exact analytical solutions of the non-dimensional formation pressure distributions depending on the non-dimensional threshold pressure gradients $\lambda =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 $\lambda =0$, the pressure drop of the Darcy seepage model can be instantaneously transmitted to infinity. When $\lambda >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 $\delta \left(t\right)$ compared with the transient distance with moving boundaries depending on different values of non-dimensional threshold pressure gradients ($\lambda =0.2,0.5,0.8,1$) when $U=0.5$. As $\lambda$ rises from 0.2 to 1, the dimensionless transient distance $\delta \left(t\right)$ increases. The value of the dimensionless transient distance $\delta \left(t\right)$ 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 $\delta \left(t\right)$ is not obvious; however, when the value of the dimensionless time $t$ is large, the value of the dimensionless transient distance $\delta \left(t\right)$ will increase significantly. As shown in Figure 4b, the distribution of the dimensionless transient distance $\delta \left(t\right)$ 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 $\lambda =0.5$. As $U$ decreased from 1 to 0.2, the dimensionless transient distance $\delta \left(t\right)$ increased. As the dimensionless time $t$ increased, the value of the dimensionless transient distance $\delta \left(t\right)$ increased. When the value of the dimensionless time $t$ was small, the increased in the value of the dimensionless transient distance $\delta \left(t\right)$ was not obvious; however, when the value of the dimensionless time $t$ was large, the value of the dimensionless transient distance $\delta \left(t\right)$ increased significantly.

As shown in Figure 5a, distributions of non-dimensional formation pressure $p\left(x,t\right)$ 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 $\lambda =0,t=5000$. As $U$ increased from 0.2 to 1, the dimensionless formation pressure $p\left(x,t\right)$ increased. As the dimensionless distance $x$ increased, the value of the dimensionless formation pressure $p\left(x,t\right)$ decreased. When the value of the dimensionless distance $x$ was small, the increase in the dimensionless formation pressure $p\left(x,t\right)$ 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\left(x,t\right)$ was not obvious as $U$ increased from 0.2 to 1. When $\lambda =0$, the dimensionless formation pressure $p\left(x,t\right)$ 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\left(x,t\right)$ 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 $\lambda =0.5,t=5000$. As $U$ increased from 0.2 to 1, the dimensionless formation pressure $p\left(x,t\right)$ increased. As the dimensionless distance $x$ increased, the value of the dimensionless formation pressure $p\left(x,t\right)$ decreased. When the non-dimensional distance $x$ was small, the increase in the non-dimensional formation pressure $p\left(x,t\right)$ 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\left(x,t\right)$ was not obvious as $U$ increased from 0.2 to 1. When the non-dimensional threshold pressure gradient was $\lambda >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 $\lambda =0$, the pressure drop of Darcy seepage models can be instantaneously transmitted to infinity. When $\lambda >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 $\lambda$.

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.

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

$\overline{p}$	Formation pressure
${\overline{p}}_i$	Original formation pressure
$C_f$	Fluid compression coefficient
$C_{\varphi }$	Pore compression coefficient
$\overline{v}$	Seepage velocity
$k$	Permeability
$\overline{x}$	Distance
$\overline{t}$	Time
$\overline{s}$	Moving boundary
$C_t$	Comprehensive compression coefficient
$f$	Bottom-hole pressure function that changes with time
${\overline{x}}_w$	Distance constant
$x$	Non-dimensional distance
$t$	Non-dimensional time
$\delta$	Non-dimensional moving boundary
$p$	Non-dimensional formation pressure
$U$	Non-dimensional positive fitting parameter of production pressure data
Greek Symbols
$\rho$	Fluid density
${\rho }_i$	Initial density of fluid
$\varphi$	Porosity
${\varphi }_i$	Initial porosity
$\mu$	Fluid viscosity
$\overline{\lambda }$	Threshold pressure gradient
$\lambda$	Non-dimensional threshold pressure gradient

Footnotes

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.

Figures and Table

Figure 1. Schematic description of an oil reservoir.

Figure 2. Curves of [Forumla omitted. See PDF.] against [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] when [Forumla omitted. See PDF.]: (a) 10th-order and 15th-order HAM approximations when [Forumla omitted. See PDF.]; (b) 10th-order and 15th-order HAM approximations when [Forumla omitted. See PDF.].

Figure 3. Comparisons between the results of dimensionless formation pressure [Forumla omitted. See PDF.] of 10th-order HAM and the exact analytical solutions of the non-dimensional formation pressure distribution depending on [Forumla omitted. See PDF.] when [Forumla omitted. See PDF.].

Figure 4. The distributions of the non-dimensional transient distance [Forumla omitted. See PDF.] of moving boundaries: (a) comparisons of the transient distance of moving boundaries depending on various values of [Forumla omitted. See PDF.] when [Forumla omitted. See PDF.]; (b) comparisons of the transient distance of moving boundaries depending on various values of [Forumla omitted. See PDF.] when [Forumla omitted. See PDF.].

Figure 5. The distribution of non-dimensional formation pressure [Forumla omitted. See PDF.]: (a) Comparison of the distributions of non-dimensional formation pressure depending on different values of [Forumla omitted. See PDF.] when [Forumla omitted. See PDF.]; (b) Comparison of the distributions of non-dimensional formation pressure depending on different values of [Forumla omitted. See PDF.] when [Forumla omitted. See PDF.].

Table 1

Comparison of the results of dimensionless formation pressure $p\left(x,t\right)$ of 10th-orderHAM and the exact analytical solutions when $\lambda =0,U=0.5,t=5000$.

$x$	$p\left(x,5000\right)\left(\mathrm{Analytical}\right)$	$p\left(x,5000\right)$ (10th-Order HAM)	Relative Error (%)
0	70.7170	70.7107	0.0089
5	66.3679	66.3589	0.0135
10	62.2017	62.1848	0.0271
15	58.2113	58.1879	0.0402
20	54.3957	54.3672	0.0524
30	47.2822	47.2482	0.0719
50	35.0585	35.0304	0.0802
100	14.7673	14.7895	0.1501
150	5.1944	5.2046	0.1948
200	1.5049	1.4704	2.3490
250	0.3552	0.3463	2.5577
300	0.0677	0.0660	2.6023

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. [DOI: https://dx.doi.org/10.1088/1674-1056/23/4/044701]

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. [DOI: https://dx.doi.org/10.1142/S0218348X18500779]

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, pp. 427-431. [DOI: https://dx.doi.org/10.1016/j.petrol.2019.04.039]

4. Hayat, T.; Yasmin, H.; Ahmad, B.; Chen, B. Simultaneous effects of convective conditions and nanoparticles on peristaltic motion. J. Mol. Liq.; 2014; 193, pp. 74-82. [DOI: https://dx.doi.org/10.1016/j.molliq.2013.12.036]

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. [DOI: https://dx.doi.org/10.1155/2021/2638647]

6. Prada, A.; Civan, F. Modification of Darcy’s law for the threshold pressure gradient. J. Petrol. Sci. Eng.; 1999; 22, pp. 237-240. [DOI: https://dx.doi.org/10.1016/S0920-4105(98)00083-7]

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, pp. 221-229. [DOI: https://dx.doi.org/10.1016/j.fuel.2018.03.192]

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, pp. 232-236.

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, pp. 1111-1117. [DOI: https://dx.doi.org/10.1021/ef1014633]

10. Yao, T.; Huang, Y.; Li, J. Nonliner flow equations for heavy oil in porous media. Chin. J. Theor. Appl. Mech.; 2012; 44, pp. 106-110.

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. [DOI: https://dx.doi.org/10.1016/j.petrol.2019.106518]

12. Wei, Q.; Zhou, H.; Yang, S. Non-Darcy flow models in porous media via Atangana-Baleanu derivative. Chaos Solitons Fractals; 2020; 141, 110335. [DOI: https://dx.doi.org/10.1016/j.chaos.2020.110335]

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. [DOI: https://dx.doi.org/10.1016/j.amc.2020.125674]

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. [DOI: https://dx.doi.org/10.1016/j.jngse.2021.104269]

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. [DOI: https://dx.doi.org/10.1016/j.seta.2022.102048]

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, pp. 2004-2013. [DOI: https://dx.doi.org/10.1016/j.petsci.2022.05.015]

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, pp. 15-17.

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, pp. 605-612.

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, pp. 6017-6022. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2012.06.012]

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, pp. 17-30. [DOI: https://dx.doi.org/10.1016/j.ijnonlinmec.2019.03.011]

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, pp. 931-953. [DOI: https://dx.doi.org/10.1016/j.apm.2020.01.028]

22. Guo, Y.; Zeng, Q.; Wang, Z. Numerical simulation of low permeability flow with moving-boundary using meshless methods. Eng. Mech.; 2006; 23, pp. 188-192.

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, pp. 879-883.

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, pp. 597-601.

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. [DOI: https://dx.doi.org/10.1155/2013/384246]

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, pp. 1112-1127. [DOI: https://dx.doi.org/10.1016/j.petrol.2019.04.033]

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. [DOI: https://dx.doi.org/10.3390/en15228505]

28. Liao, S. The Proposed Homotopy Analysis Technique for the Solution of Nonlinear Problems. Ph.D. Thesis; Shanghai Jiao Tong University: Shanghai, China, 1992.

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.

30. Liao, S. Homotopy Analysis Method in Nonlinear Differential Equations; 1st ed. Springer: Berlin/Heidelberg, Germany, 2012; pp. 40-86.

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, pp. 2501-2512. [DOI: https://dx.doi.org/10.1016/j.cnsns.2010.09.015]

32. Vajravelu, K.; Van Gorder, R.A. Nonlinear Flow Phenomena and Homotopy Analysis: Fluid Flow and Heat Transfer; Springer: Berlin/Heidelberg, Germany, 2012.

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, pp. 484-495. [DOI: https://dx.doi.org/10.1016/j.amc.2014.12.005]

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. [DOI: https://dx.doi.org/10.1063/1.4986822]

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, pp. 21-28. [DOI: https://dx.doi.org/10.15415/mjis.2018.71004]

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. [DOI: https://dx.doi.org/10.1103/PhysRevE.103.023206] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33736054]

37. Yang, X.Y.; Li, Y. On bi-chromatic steady-state gravity waves with an arbitrary included angle. Phys. Fluids; 2022; 34, 032107. [DOI: https://dx.doi.org/10.1063/5.0083705]

38. Liao, S. A new non-perturbative approach in quantum mechanics for time-independent Schrödinger equations. Sci. China Phys. Mech.; 2020; 63, 234612. [DOI: https://dx.doi.org/10.1007/s11433-019-9430-4]

39. Liao, S. Avoiding small denominator problems by means of the Homotopy Analysis Method. Adv. Appl. Math. Mech.; 2023; 15, pp. 267-299. [DOI: https://dx.doi.org/10.4208/aamm.OA-2022-0260]

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, pp. 5537-5545. [DOI: https://dx.doi.org/10.1016/j.aej.2021.04.061]

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, pp. 505-522. [DOI: https://dx.doi.org/10.1016/j.matcom.2021.11.027]

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, pp. 486-498. [DOI: https://dx.doi.org/10.1016/j.apm.2021.12.001]

43. Mehmet, P.; Muhammet, Y. Similarity transformations for partial differential equations. SIAM Rev.; 1998; 40, pp. 96-101.

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, pp. 5387-5390. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2004.07.007]

45. Liao, S. An explicit, totally analytic approximate solution for Blasius’ viscous flow problems. Int. J. Non Linear Mech.; 1999; 34, pp. 759-778. [DOI: https://dx.doi.org/10.1016/S0020-7462(98)00056-0]

46. Liao, S. On the analytical solution of magnetohydrodynamic flows of non-Newtonian fluids over a stretching sheet. J. Fluid Mech.; 2003; 488, pp. 189-212. [DOI: https://dx.doi.org/10.1017/S0022112003004865]

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, pp. 507-514. [DOI: https://dx.doi.org/10.1002/zamm.200700132]

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. [DOI: https://dx.doi.org/10.1115/1.3173764]

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, pp. 102-108. [DOI: https://dx.doi.org/10.1016/j.icheatmasstransfer.2014.04.016]

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, pp. 69-75. [DOI: https://dx.doi.org/10.1016/j.compfluid.2015.01.005]

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, pp. 70-76. [DOI: https://dx.doi.org/10.5560/zna.2011-0057]

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, pp. 224-238. [DOI: https://dx.doi.org/10.4236/ajcm.2015.52019]

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. [DOI: https://dx.doi.org/10.3390/nano12193289]

54. Yang, Y.; Liao, S. Comparison between homotopy analysis method and homotopy renormalization method in fluid mechanics. Eur. J. Mech. B Fluids; 2023; 97, pp. 187-198. [DOI: https://dx.doi.org/10.1016/j.euromechflu.2022.10.005]