Fractional Derivative Model on Physical Fractal Space: Improving Rock Permeability Analysis

Abstract

As challenges in gas extraction from coal mines increase, precise measurement of permeability becomes crucial. This study proposes a novel pulse transient method based on a fractional derivative model derived on physical fractal space, incorporating operator algebra and the mechanics–electricity analogy to derive a new control equation that more accurately delineates the permeability evolution in coal. To validate the approach, permeability experiments were conducted on coal samples under mining stress conditions. The results showed that the adoption of a physically meaningful fractional-order relaxation equation provides a more accurate description of non-Darcy flow behaviour in rocks than traditional integer-order control equations. Additionally, the method proved effective across different rock types, verifying its broad applicability. By establishing a new theoretical foundation, this approach illustrates how the microscale fractal structure of rocks is fundamentally linked to their macroscale fractional responses, thereby enhancing the understanding of fractional modelling methods in rock mechanics and related domains.

1. Introduction

For coal mines, gas is both the culprit of coal mine disasters and an irreplaceable clean energy source [1,2,3]. With the increasing depth of coal mining, gas pressure and content gradually increase while the permeability decreases, posing great challenges to gas extraction [4,5,6]. Permeability is a key parameter in the theory of gas extraction in coal mines, and it is a prerequisite for establishing permeability models and studying the law of gas flow [7,8,9]. Therefore, accurately measuring the permeability of coal samples in the laboratory is of great significance for guiding gas extraction technology and preventing coal mine gas disasters.

Numerous scholars have conducted indoor experiments to investigate the evolution of coal permeability [10]. Permeability testing methods are typically classified into two categories: steady-state and transient methods. However, for rocks with low permeability, such as coal and granite, achieving steady-state fluid flow in experiments often requires a considerable amount of time, making the use of steady-state methods for permeability determination highly time-consuming. In response to this challenge, Brace et al. [11] proposed a novel testing approach to assess the permeability of tight rocks, known as the transient pulse method. Leveraging Darcy’s law, Brace et al. [11] derived a simplified solution in the form of a negative exponential by neglecting rock porosity. Additionally, Hsieh et al. [12] presented an analytical solution in the form of a series; however, its complexity often leads to the preference for simplified solutions in practical research.

Nevertheless, the loaded coal specimens form numerous pores and fractures after damage, making them difficult to ignore [13,14]. Additionally, the roughness of the pore and fracture structure (PFS) and the tortuosity of fluid flow increase solid–liquid interactions, resulting in non-Darcy flow characteristics [15,16,17]. This leads to pressure experimental data exhibiting a thick-tailed feature, with a relaxation rate slower than that of negative exponential relaxation. This phenomenon has been observed in many experiments, resulting in lower fitting effectiveness of negative exponential relaxation [18,19,20,21,22,23]. Given the limitations of fitting negative exponential relationships, Zhao et al. [20] proposed a polynomial fitting method, but it involves too many parameters and lacks physical significance. Yang et al. [21] used the Mittag–Leffler function to simulate the pressure relaxation characteristics in transient pulse experiments in order to determine the permeability of rock samples. This method is primarily based on mathematical theory. Inspired by this, An et al. [22] employed the Mittag–Leffler function to simulate the pressure relaxation characteristics in transient pulse experiments and then established a relationship between fractional order and fracture connectivity, reflecting the original rock sample’s density and the degree of connectivity of internal fracture networks under triaxial compression. Zhou et al. [23] further elaborated on the relaxation of coal pressure in the context of mining stress by developing a novel fractional order relaxation model. They provided a preliminary discussion on the physical significance of the fractional order and demonstrated that this model can more accurately depict the impact of mining stress on non-Darcy flow of coal.

In the domain of rock mechanics, fractional calculus plays a critical role. It efficaciously addresses a wide array of complex nonlinear issues, including anomalous diffusion phenomena [24], accelerated phases of creep [25], and pore compression models [14], thereby demonstrating its unparalleled advantage in tackling both low- and high-speed non-Darcy flow challenges [26]. The uniqueness of this mathematical tool lies in its ability to simulate memory effects that traditional calculus struggles to capture [27,28]. Despite achieving significant success in describing and predicting experimental data, many fractional order models still face a significant challenge: the lack of a sufficient physical basis. Researchers often rely on phenomenological modeling based on observed experimental phenomena, yet struggle to delve into the physical essence behind these phenomena. This limitation not only affects the interpretability of the models, but also restricts the application of fractional calculus methods in addressing a broader range of scientific issues.

Excitingly, recent research in the field of biomechanics has proposed an innovative perspective, starting from the physical fractal space to directly derive fractional order responses. For instance, Guo et al. [29] revealed that the multilevel self-similar structure of ligaments and tendons could induce a 1/2 order viscoelastic response through analyzing their physical fractal space. Similarly, in the study of neural spike signal transmission, the physical fractal layout of multilevel dendritic structures was observed to cause a 1/2 order neural electrical signal conduction characteristic [30]. Furthermore, Jian et al. [31,32] presented an abstract mechanical description of bones and derived fractal operators within the physical fractal space, thereby establishing a fractional order relaxation model. These studies indicate that the time fractional order fractal operators, originating from the motion in physical fractal space structures or on them, are a common phenomenon in biological entities, whether it be the fractal characteristics of myofibrils, nerve fibers, or bones. Based on this, we have reasons to believe that there exists a close connection between the microscopic fractal structure of coal and its macroscopic fractional order response, offering a new theoretical basis for developing fractional transient methods with actual physical significance. More importantly, Peng et al. [33] introduced the physical fractal space structure in the study of arterial hemodynamics and revealed through the force–electricity analogy method that the flow of blood in small arteries is regulated by fractional calculus, further proving the potential of this approach in explaining complex seepage processes.

Therefore, to reveal the evolution law of permeability more accurately in coals, this study proposes an improvement to the control equation of the pulse transient method by introducing a physically meaningful fractional order relaxation equation to describe the physical fractal characteristics of fluid flow within coal bodies. To verify the accuracy of the proposed model, we conducted permeability tests on coal samples under mining stress conditions and applied the fractional order relaxation equation to process these data. Moreover, by fitting and analyzing the pressure difference curves of other rock samples in the existing literature, this study further explores the broad applicability of the proposed model. Finally, our discussion emphasizes the general interrelation between fractional order dynamic response characteristics and physical fractal space structures in rock mechanics problems.

2. Permeability Measurements by Transient Pulse Testing

2.1. Traditional Perspective: Relaxation of Integer-Order Control Equations with Negative Exponentials

Brace et al. [11] developed a transient method founded on Darcy’s Law, which involves a rock specimen interfaced with two reservoirs at its extremities, as shown in Figure 1. Initially, these reservoirs are in a state of pressure equilibrium. At the beginning of the experiment, an instantaneous pressure pulse is applied to the upstream reservoir, disrupting this equilibrium and inducing a pressure differential across the rock sample’s extremities, thereby commencing fluid permeation. The permeability can be derived by analyzing the relaxation process of the pressure differential between the upstream and downstream reservoirs. Brace et al. [11] propose that the reduction in pressure differential over time, neglecting rock porosity, can be described by a negative exponential function, which is given as:

$$P\left(t\right)=P_u\left(t\right)-P_d\left(t\right)=P_0\exp (-\gamma t)$$

(1)

$$\gamma ={\displaystyle \frac{kA}{\mu {\beta }_{w}L}}\left({\displaystyle \frac{1}{V_d}}+{\displaystyle \frac{1}{V_u}}\right)$$

(2)

where P(t) represents the pressure difference (Pa), P₀ is the instantaneous pressure pulse or referred to as the initial pressure difference (Pa), k denotes the permeability of the rock specimen (m²), A is the cross-sectional area of the rock specimen (m²), L is the height of the rock specimen (m), μ is the dynamic viscosity of the permeating fluid (Pa·s), β_w is the fluid compressibility coefficient (Pa⁻¹), and V_u and V_d, respectively, represent the volumes of the upper and lower reservoirs (m³).

Therefore, the widely used formula for calculating permeability using the transient method can be expressed as:

$$k={\displaystyle \frac{\mu {\beta }_{w}L}{At\left(\frac{1}{V_d}+\frac{1}{V_u}\right)}}\ln \left({\displaystyle \frac{P_0}{P\left(t\right)}}\right)$$

(3)

The most crucial aspect of permeability calculation is the determination of the parameter γ in Equation (3) through optimal fitting analysis of experimental pressure difference data. However, numerous experimental studies have shown that the pressure relaxation characteristics in transient pulse tests sometimes deviate from the negative exponential relaxation process due to the non-Darcy flow characteristics of fluid percolation. For instance, the fitting accuracy for experimental data from coal samples’ permeability tests using the transient method, as described in Zhou et al.’s study [34], is relatively low (R² = 0.7583), as demonstrated in Figure 2.

2.2. Fractal Perspective: Relaxation of Fractional-Order Control Equations with Mittag–Leffler Functions

In fact, during mechanical loading experiments conducted concurrently with permeability testing, the damage caused by loading often makes the porosity of the rock difficult to ignore. Additionally, rocks like coal typically exhibit a wide range of pore size distributions, complex PFS, and distinct fractal characteristics [13,14]. These features lead to more complex flow geometries and rough contacts, as well as stronger solid–liquid interactions, thus exhibiting non-Darcy flow characteristics. Therefore, this section utilizes fractal theory to establish a new transient method permeability control equation.

2.2.1. Mechanical–Electrical Analogy of Rock Permeability Parameters

To establish a physical fractal model simulating fluid flow behaviour within rocks, a mechanical–electrical analogy of rock permeability parameters is initially conducted. Instantaneous pressure pulses, or pressure differentials P(t), are equated to voltage u, while flow rates Q are equated to current i. The rock and the reservoirs at either end are treated as components within an electrical circuit. The resistance encountered by fluid flowing through the rock, stemming from pores, particles, and other obstacles, is dependent on fluid viscosity, flow path, and PFS and can be described using Darcy’s law.

$$\Delta P=-{\displaystyle \frac{\mu L}{kA}}Q$$

(4)

This behaviour can be equated to electrical resistance, with its resistance value given by:

$$R={\displaystyle \frac{\mu L}{kA}}$$

(5)

The fluid can be stored in the reservoirs at both ends. When the fluid pressure increases, the fluid volume in the reservoir increases, and when the fluid pressure decreases, the fluid volume in the reservoir decreases. This behaviour is analogous to a capacitor, where the equivalent capacitance depends on the reservoir volume and the compressibility of the fluid:

$$C_u={\beta }_w{V}_u$$

(6)

$$C_d={\beta }_w{V}_d$$

(7)

The two reservoirs can be viewed as two capacitors connected in series, with their equivalent capacitance given by:

$${\displaystyle \frac{1}{C_r}}={\displaystyle \frac{1}{C_u}}+{\displaystyle \frac{1}{C_d}}={\displaystyle \frac{1}{{\beta }_w}}\left({\displaystyle \frac{1}{V_d}}+{\displaystyle \frac{1}{V_u}}\right)$$

(8)

If the water storage capacity of the pores in the rock is ignored, then transient pulse testing can be equivalent to the discharge process of a simple RC circuit. The discharge process of an RC circuit is a typical negative exponential relaxation process:

$$u=u_0\exp (-{\displaystyle \frac{1}{R{C}_{\mathrm{r}}}}t)$$

(9)

Note that in Equation (9), γ = 1/RC_r. This equivalence arises because, by substituting the expressions for R (Equation (5)) and C_r (Equation (8)) into the RC circuit model, we obtain the same decay constant γ as in Equation (2). This shows that the mechanical–electrical analogy holds, where voltage u corresponds to the pressure differential P(t). Consequently, the result derived using the electrical analogy (Equation (9)) matches the traditional method’s control equation for permeability testing (Equation (1)), confirming the consistency and applicability of the analogy.

In experiments, the observed negative exponential relaxation does not hold true. This is primarily due to the significant water storage capacity of rock pores in samples with larger porosity. Consequently, the rock cannot be accurately represented merely as a resistor. Moreover, as loading advances, numerous fractures emerge, intensifying the internal connectivity of the rock [13,14]. This increased connectivity further underscores the significance of the water storage capacity within the pores, rendering it impossible to disregard.

2.2.2. Abstraction of a Physical Fractal Circuit for Pore-Fracture Structures

In the study of fluid flow within rocks, the “capillary bundle model” is commonly employed to simplify the complex PFS. This model assumes that the structure within the rock comprises a set of cylindrical capillaries with uniform diameters, simulating fluid flow channels (Figure 3a). Each capillary’s flow is assumed to occur independently, without interaction between the capillaries. However, numerous experiments, including those by Liu et al. [35] based on CT scans (Figure 3b), indicate that the capillaries are not truly independent, suggesting the need for model modification.

Modifications to the capillary bundle model include transforming a single straight capillary into multiple interconnected capillaries (Figure 3c). From the perspective of electro-hydraulic analogy, fluid flow in capillaries follows Darcy’s law, akin to electrical resistance in circuits. Branching of flow paths results in fluid leakage into other pores, analogous to the capacitive effect in circuits.

Considering that the PFS of rocks typically exhibits clear fractal characteristics [13,14]; to simplify the model, this study makes the following assumptions:

• Branches leading to other capillaries are uniformly distributed along the current capillary’s axis.
• Only the axial resistance effect of the current capillary is considered, ignoring its water storage capacity.
• The branch capillaries of the flow path follow the same pattern.

Based on these assumptions, in each unit segment, a resistor and a capacitor are connected in parallel (Figure 3d), and each capacitor connects to a new resistor–capacitor network. When these unit segments are connected in series, they form a tree-like topology with infinite branches, creating a physical fractal circuit (Figure 4). The connection rules between elements when analysing fluid transport in the capillary network include:

• Parallel connection of a resistor element R and a capacitor element C to form a basic unit.
• Replication of the basic unit to create a secondary structure.
• Series connection of the secondary structure with the resistor R and capacitor C elements in the basic structure.
• Replication of this structure again to form a tertiary structure and series connection with the elements in the secondary structure.

Repeating this construction according to these rules can form a multi-level branched tree-like physical fractal structure (Figure 4). Notably, this physical fractal circuit does not depend on geometric scale or fractal dimension, but manifests as an infinitely self-similar structure with inherent physical functions. This structure can simulate complex physical processes, surpassing the simplicity of geometric repetition. However, it is important to note that the assumptions made in this model, such as uniform branch distribution, may limit its accuracy in representing natural systems.

2.2.3. Fractal Admittance Operator Based on Operator Algebra Method

To utilize circuit simulation for the rock seepage process, it is necessary to solve the infinite-level fractal circuit first. However, conventional calculation methods would entail enormous computational costs and might even fail to yield a theoretical solution. Nonetheless, the infinite-level fractal circuit of capillaries possesses fractal characteristics, providing us an opportunity to compute the circuit response.

The control equation of the circuit is usually a complex system of differential equations. To simplify the computation, we adopt an operational method to handle the differential equations that describe the admittance characteristics of the physical fractal circuit by introducing a fractal admittance operator to describe these equations. This method draws from Heaviside’s technique [36,37] to streamline calculus operations, converting the complex system of differential equations into straightforward algebraic equations involving operators, allowing us to use operators to directly describe the modulation mechanism of the circuit.

To represent the admittance characteristics of the basic physical components in fractal circuits, we utilize the operator T(p):

$$T(p)={\displaystyle \frac{i(t)}{u(t)}}$$

(10)

where i(t) and u(t) denote the time-varying current and voltage of the electrical component, respectively, and p is the differential operator with respect to time. The function T(p) encapsulates the admittance properties as a function of the operator p.

In the framework of operator algebra, the differential operator pp is defined for functions f(t) with continuous derivatives for t ≥ 0. The action of p on such a function f(t) is given by [36,37]:

$$pf\left(t\right)={\displaystyle \frac{\mathrm{d}f(t)}{\mathrm{d}t}}$$

(11)

For a capacitive component C, its admittance operator T_C is:

$$T_C=Cp$$

(12)

The equivalent capacitance is determined by the pore compressibility under effective stress and the sample volume:

$$Cs={\beta }_{s}AL$$

(13)

where β_s represents the pore compression coefficient under the action of effective stress (Pa⁻¹).

For a resistive component R, its admittance operator T_R is:

$$T_R={\displaystyle \frac{1}{R}}$$

(14)

The second step involves deducing the overall admittance operator TF for fractal circuits. To note, the components (resistors and capacitors) in each segment of the physical fractal structure studied in this paper are linear, thus allowing us to use operator algebra to handle these equations. This approach is valid if the governing differential equations are linear, which is an inherent limitation to be considered [36,37].

The fractal tree structure composed of resistor elements R and capacitor elements C can be viewed as a composite fractal component F (Figure 5a). The premise of this model is that the overall fractal structure possesses infinite levels of self-similarity; that is, only when the structure reaches an infinite level does it exhibit equivalence to the fractal element F [38]. When examining the first resistor or capacitor element (Figure 4), one can immediately observe another fractal tree structure (or another fractal element F) in the second-level structure. Thus, the entire fractal step structure is equivalent to a simplified structure: a resistor element R in series with a fractal element F, and similarly, a capacitor element C in series with another fractal element F. These two series structures are then connected in parallel. This equivalent simple structure is referred to as a fractal cell [29] (Figure 5b).

Based on the equivalence between fractal cell (Figure 5b) and fractal component (Figure 5a), an operator method can be used to solve the admittance operator of the fractal component. For the fractal circuit, through the series and parallel connections of the elements, the relationship between the admittance operator of the fractal cell and the admittance operator of the fractal component can be derived as follows:

$$T_F={\displaystyle \frac{T_F{T}_R}{T_F+T_R}}+{\displaystyle \frac{T_F{T}_c}{T_F+T_c}}$$

(15)

It is important to emphasize that Equation (15) is identified from the concept of infinite-level fractal structures rather than being derived through rigorous mathematical proof. The concept of infinite structural hierarchy is not a natural physical manifestation, but emerges from a logical and abstract framework. Fortunately, Yu et al. [39] successfully demonstrated within the scope of operator theory the monotonic boundedness of the fractal structure operator sequence, verifying the existence of a unique limit for such sequences. They discovered that the overall mechanical behavior of structures with more than three finite levels can be approximated by an infinite-level fractal structure. Solving the algebraic Equation (15) yields the precise formula for the fractal admittance operator:

$$T_F=\sqrt{T_C{T}_R}$$

(16)

Equation (16) takes the form of a quadratic root solution, indicating that the fractal admittance operator T_F exhibits characteristics of a fractional order derivative of order 1/2.

2.2.4. Fractional Order Control Equations Derived from Fractal Admittance Operators

Considering the rock’s porosity and its ability to store water, the rock can be equivalently modelled as a fractal tree-like structural component T_F whose behaviour follows the principles outlined in Equation (16). Additionally, the water troughs at both ends of the rock can be regarded as two capacitors in series, represented here as T_C2 (Figure 6).

When an instantaneous pressure pulse is applied to the upstream reservoir, this is equivalent to the discharge process of the capacitor T_C2 with an initial voltage of u₀. During this process, for the capacitor T_C2, it is given by:

$$i=-C_{r}p{u}_{\mathrm{c}}$$

(17)

where u_C represents the voltage across the capacitor (V). For the fractal tree component T_F, it is given by:

$$u_T={\displaystyle \frac{1}{T_F}}i=-C_{r}p{u}_{\mathrm{c}}{\displaystyle \frac{1}{T_F}}=-C_r{u}_{\mathrm{c}}\sqrt{{\displaystyle \frac{Rp}{C_S}}}$$

(18)

where u_T represents the voltage across T_F (V). According to the Kirchhoff Voltage Law, it can be derived that the sum of voltage drops around any closed loop in the circuit must be zero:

$$C_r{u}_{\mathrm{c}}\sqrt{{\displaystyle \frac{Rp}{C_S}}}+u_{\mathrm{c}}=0$$

(19)

This law is applicable here because the pressure differences and flow rates in the rock samples can be analogously treated as voltage and current in an electrical circuit. From this, we derive:

$$\sqrt{p}{u}_c=-\alpha u_c\iff {\displaystyle \frac{{\mathrm{d}}^{1/2}{u}_c}{d{t}^{1/2}}}=-\alpha u_c$$

(20)

$$\alpha ={\displaystyle \frac{1}{C_r}}\sqrt{{\displaystyle \frac{C_S}{R}}}={\displaystyle \frac{1}{{\beta }_w}}\left({\displaystyle \frac{1}{V_d}}+{\displaystyle \frac{1}{V_u}}\right)\sqrt{{\beta }_{\mathrm{s}}{\displaystyle \frac{kA}{\mu L}}}$$

(21)

We employ the definition of the Caputo fractional derivative [40], given by:

$$\begin{array}{l}{}_0{}^{C}D_x^{\alpha }f(x)={\displaystyle \frac{1}{\Gamma (n-\alpha )}}{\displaystyle \underset{0}{\overset{x}{\int }}\frac{f^{(n)}(t)}{{(x-t)}^{\alpha -n+1}}dt},\hspace{1em}\alpha \notin \mathrm{N}\\ {}_0{}^{C}D_x^{\alpha }f(x)=f^{(\alpha )}(x),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \alpha \in \mathrm{N}\end{array}$$

(22)

where N is a set of positive integers and Γ is the Gamma function. We then apply the Laplace transform and its inverse to obtain the solution as:

$$u_c=u_0{E}_{0.5}\left(-\alpha t^{0.5}\right)$$

(23)

where E_λ(z) represents the Mittag–Leffler function [41]. Therefore, the permeability calculation formula incorporating physical fractal characteristics can be expressed as:

$$k={\displaystyle \frac{{\beta }_w^2}{t{\left(\frac{1}{V_d}+\frac{1}{V_u}\right)}^2}}{\left[-E_{0.5}^{-1}\left({\displaystyle \frac{P\left(t\right)}{P_0}}\right)\right]}^2{\displaystyle \frac{\mu L}{{\beta }_{S}A}}$$

(24)

Consequently, rocks with a fractal tree-like structure follow the Mittag–Leffler relaxation process for the variation of osmotic pressure, rather than the traditional negative exponential relaxation process.

3. Experimental Validation of the Fractional Order Pulse Transient Method

3.1. Specimens and Experimental Equipment

To validate the fractional order pulse transient method model delineated in Section 2, this section presents the execution of permeability tests on coal specimens procured from the No. 8 coal mining face at Dongqu Mine, Shanxi, China. Located approximately 400 m underground with an initial geostress of around 10 MPa, the chosen site provided samples for this study. To ensure the integrity of the coal samples during transit, the samples were encased in multiple layers of cling film and anti-vibration foam before being dispatched to the facility for preparation. Subsequent processes involved drilling, coring, cutting, and polishing the coal to yield six cylindrical specimens, each standing 100 mm tall and 50 mm in diameter. Three of the most intact samples were selected and transported to the laboratory for advanced experimentation.

The experimental phase employed the American MTS815 Flex Test GT rock mechanics testing system, acclaimed for its incorporation of triaxial loading, permeability pressure control, and temperature management, among other capabilities, positioning it as a leading solution in global rock mechanics examination. The entire experimental procedure was governed by a pre-defined computing program, eliminating the need for manual adjustment across different loading modes. Utilization of high-precision sensors and computing technology facilitated the real-time monitoring and automatic documentation of critical physical parameters such as load, displacement, pressure, and time throughout the experimental timeline.

3.2. Experimental Procedure

During the coal extraction process, the coal body at the mining face experiences a comprehensive mechanical evolution process. This starts from the original in situ stresses, progressing through an increase in vertical stress and a slow reduction in horizontal stress, and culminating in instability and rupture. This sequence of events modifies the coal’s PFS, along with its permeability [42]. In the lab, vertical stress is simulated by axial pressure and horizontal stress by confining pressure, with the designed stress path for coal sample specimens under simulated mining conditions depicted in Figure 7.

The experimental procedure is detailed as follows:

(1)

Axial and confining pressures are simultaneously increased at a rate of 10 kN/min to a hydrostatic pressure of 10 MPa.

(2)

The axial pressure is raised to its peak value, while the confining pressure is gradually decreased in 1 MPa steps. Permeability is measured three times before the peak stress is reached.

(3)

After achieving the specimen’s peak stress, a permeability measurement is conducted, and the method of axial loading is switched to displacement control at a rate of 0.06 mm/min.

(4)

Axial pressure is progressively decreased in 20% stress decrements, with confining pressure also being reduced in 1 MPa steps until the specimen fails completely, with permeability measured three more times after the peak stress.

3.3. Test Results and Analysis

In this section, the fractional order relaxation control equation was employed for fitting and analyzing the transient permeability test data under different experimental conditions (Equation (24)) to verify its applicability. Figure 8, taking coal sample D2 as an example, illustrates the non-Darcy flow pressure relaxation curve derived from transient method measurements of permeability evolution. A total of six permeability measurements were conducted on coal sample D2 during the triaxial compression failure process. It is evident from Figure 8 that the exponential relaxation deviates from the non-Darcy flow pressure test data. In contrast, the fractional order relaxation control equation based on physical fractal space can describe the pressure relaxation in transient pulse tests more accurately than the traditional integer order relaxation control equation used in transient methods, thereby ensuring higher precision in the calculation of coal sample permeability.

The same approach was applied to fit and compare the pressure test data from other coal permeability tests, and the results also indicate that the fractional order relaxation control equation based on physical fractal space can accurately determine the permeability of coal. More specifically, Figure 9 presents the curves of axial strain ε₁ versus deviatoric stress σ₁–σ₃ for different coal samples, along with the permeability–strain relationship curves. It can be observed from Figure 9 that, compared to the reference permeability k calculated based on Equation (3), the permeability k_F considering the physical fractal space of coal samples and calculated based on Equation (24) is more precise. Moreover, the permeability evolution patterns of coals calculated by Equations (3) and (24) are essentially consistent, suggesting that the fractional order transient method proposed in this study can fully replace the traditional transient method to more accurately determine the evolution characteristics of rock permeability.

The permeability measurements obtained from coal samples show some variability. However, the coal samples generally exhibit an initial permeability around the magnitude between 10⁻¹⁷ and 10⁻¹⁶ m². Upon reaching the peak stress, the permeability rises by one order of magnitude, reaching between 10⁻¹⁷ and 10⁻¹⁶ m². It is also noted that fewer permeability readings are available after the peak, a challenge attributed to the difficulty in controlling the samples post-peak. The rapid decline in post-peak stress results in overlapping measurement points, causing the control program of the MTS815 experimental system to crash. Despite these challenges, the data reveal a distinct pattern of permeability evolution.

4. Discussion

4.1. Generalized Model of Permeability Evolution

To better understand the evolution of coal permeability under mining-induced stress, Figure 10 presents a generalized model of permeability evolution. This model is based on the experimental results shown in Figure 9. Throughout the entire loading process, the permeability of the coal sample decreases by 1–2 times during the compaction and elastic deformation stages and then increases rapidly. Before reaching the peak, the permeability increases slowly; after reaching the peak, the permeability rises sharply until the coal sample completely fractures, with permeability increasing by 2–8 times compared to the initial value. In the post-peak stage, the permeability k_F calculated based on the physical fractal space of the coal sample, and that calculated by Equation (24) is more accurate than the reference permeability k calculated by Equation (3), with the difference being particularly significant.

The reason for this difference lies in the fact that the control equation of the transient pulse method (Equation (3)) proposed by Brace et al. is based on dense rock granite. Moreover, Brace and his contemporaries primarily focused on the process of permeability reduction under the influence of confining pressure. During this stage, the rock is in a state of compression or elastic deformation, and the porosity, as well as permeability, decrease with increasing effective stress. Therefore, the assumption that the pore volume is negligible in Equation (3) is essentially valid.

However, with the establishment of low-carbon energy strategic goals, an increasing number of scholars are focusing on the permeability evolution process of unconventional reservoirs with high porosity, such as coal, especially under the influence of mining activities or coal body fracturing and damage. These studies not only cover the compression and elastic deformation stages, but also involve the plastic stage and the post-peak failure stage. In these two stages, many new PFS form within the coal sample, making the porosity non-negligible, and the phenomenon of non-Darcy flow becomes more apparent. Consequently, the integer-order relaxation control equation used in traditional transient methods becomes severely invalid.

In this context, the fractional-order control equation from the perspective of physical fractals more accurately describes the permeability evolution process (Figure 10).

4.2. Application Scope of the Fractional Order Pulse Transient Method

Section 3.3 validates the fractional order relaxation model. To further explore its applicability, we fitted the model to experimental data from Kato et al. [43] and Zhao et al. [20], as shown in Figure 11. The material used by Kato et al. [43] was clay, while Zhao et al. [20] conducted their experiments with limestone. As demonstrated in Figure 11, the fractional order relaxation control equation based on physical fractal space accurately describes the experimental data from both Kato et al. [43] and Zhao et al. [20]. To date, three sets of experimental data have validated the fractional order relaxation model, indicating its applicability for describing non-Darcy flow in various types of rocks. However, to strengthen the model’s generalizability, it is essential to conduct further experiments on a wider variety of rock types from different geographical regions. This will help to establish a more robust understanding of the model’s performance across diverse geological settings.

4.3. Fractional Order Responses in Rock Physical Fractal Space

Fractional calculus, as an effective mathematical tool for addressing nonlinear problems, has shown exceptional efficacy in fields such as rock mechanics. This approach has achieved remarkable results in areas including creep behaviour [44], anomalous diffusion [24], non-Darcy flow [26], permeability testing methods [21], porosity evolution models [14], and permeability evolution models [45]. These models not only excel in data fitting, but also succinctly capture the essence of the problems, sparking interest in the underlying principles of fractional derivatives. However, despite the excellent fitting results of these models, the order of the fractional derivatives often lacks a clear physical interpretation and is typically determined through data fitting, highlighting a significant gap in the research. Thus, a theoretical framework is needed to provide theoretical justification and physical explanations for fractional derivative modelling.

In this regard, the studies by Guo et al. [29,30], Jian et al. [31,32], and Peng et al. [33] are noteworthy. They conducted abstract analyses of the physical fractal space of biological entities, deriving fractal operators and successfully establishing fractional derivative models. These findings suggest that time-fractional fractal operators can originate from the dynamics of physical fractal spaces or their motion, a phenomenon widely observed in nature. For example, the fractal characteristics of myofibrils, nerve fibres, bones, and blood flow all validate this theory. This study further explores the relationship between the permeability behaviour of coal in physical fractal spaces and its non-Darcy fractional response, providing a new perspective for establishing fractional models with practical physical significance. We believe that the microscopic physical fractal space of rocks fundamentally determines their macroscopic fractional order response, offering important theoretical support for the understanding and application of fractional derivatives in rock mechanics and related fields.

5. Conclusions

This study introduces an innovative fractional order derivative-based pulse transient method for measuring the permeability of rocks, particularly for those with complex PFS like coal. By integrating fractional order derivatives with physical fractal space, our research not only improves the accuracy of permeability measurements, but also offers new insights and methodologies within the field of rock mechanics. The main conclusions drawn from this study are as follows:

(1)

Compared to traditional integer order control equations, the fractional order control equations more accurately describe the non-Darcy flow behaviour in coals, thereby improving the precision of permeability measurements.

(2)

Through fitting experimental data from other rock types, such as clay and limestone, this study validates the wide applicability of the fractional order method. This suggests that fractional order derivatives hold potential advantages in describing complex seepage processes, serving as a powerful analytical tool for the fields of rock mechanics and geological engineering.

(3)

By investigating the relationship between the microscale fractal structure of rocks and their macroscale fractional order response, this paper lays a new theoretical foundation for the physical significance of fractional order derivatives. The extension of this theory not only deepens our understanding of the complex internal flow mechanisms in rocks, but may also facilitate the application of fractional models in rock mechanics and related domains.

The concept of physical fractal space structures is not only applicable to the flow channels in rocks, but also to other phenomena with a physical fractal background, which warrant further exploration for their fractional order characteristics. For instance, similar physical processes can be observed in micro- and nanoscale heat transfer in thermodynamics, which can also be described using a bifurcated physical functional fractal structure model. This model may exhibit a fractional order time response of 1/2, providing a broader perspective for the application of fractional order theories in various scientific fields.