Fractal Permeability Model of Newtonian Fluids in Rough Fractured Dual Porous Media

Based on the statistical self-similar fractal characteristics of the microstructure of porous media, the total flow rate and permeability of Newtonian fluids in the rough fracture network and rough matrix pores are derived, respectively. According to the connection structure between fractures and pores, the permeability analysis model of fluids in a matrix-embedded fracture network is established. The comparison between the predicted values of the model and the experimental data shows that the predicted values of the permeability of the rough fracture network and the rough matrix pores decrease with the increase in the relative roughness of the fractures and matrix pores, and are lower than the experimental data. Meanwhile, the predicted total flow rate of a rough fractured dual porous media is lower than that of a smooth fractal model and experimental data. In addition, it is also found that the larger the average inclination angle and the relative roughness of the fracture network, the smaller the permeability of the fractured dual porous media, and the relative roughness of the fracture network has a far greater influence on fluid permeability in the fractured dual porous media than the relative roughness of the matrix pores.


Introduction
Natural porous media, such as oil and gas reservoirs, coal reservoirs, mineral resources, groundwater resources, and rock and soil, not only contain micropores, but also have a large number of large-scale channels and fractures. The study on the seepage characteristics of the fluid in this kind of rough fracture dual porous media has important significance and wide practical application fields, such as the spontaneous self-priming of fluids in rough porous media [1], effective thermal conductivity and heat and mass transfer [2,3], development of tight oil and gas reservoirs, transportation of natural gas in complex pores of tight sandstone [4,5], control of working fluid leakage and reservoir protection in fractured oil and gas reservoirs [6], regular packing in industrial separation processes [7], fractal study on electrolyte diffusion in charged porous media [8], and the improvement of fuel cell performance [9].
Previous studies have shown that the pore size distribution of porous media satisfies the fractal relationship [10,11]. Yu et al. [12] gave the statistical characteristics of porous media based on fractal theory, proposed a unified model suitable for fractal porous media, and gave the criteria for judging whether porous media can be studied by fractal theory. Since the seepage characteristics of fluids in porous media are affected by many factors. Considering the shape, the spatial distribution of pores, and the characteristics of the fluids in rock formations, scholars have proposed various models and methods to simulate fluid flow in the porous medium. Based on fractal features, Yu et al. [13] considered factors such as tortuosity, particle size, pore area, and porosity, and used the box counting method to calculate the fractal dimension, established a bidisperse porous media model and verified its effectiveness. Based on fractal theory, Wang et al. [14] equated fractured rocks with rectangular and elliptical pipes with cross-sections, respectively, and established a prediction model for the total flow rate and effective permeability of power-law fluid in porous media. Zhu et al. [15] used the tree-like fractal bifurcation network to characterize the fracturing fracture morphology of low permeability reservoirs, established the permeability model of the fracturing fracture network with elliptical and rectangular cross-sections, and discussed the influence of fracture bifurcation series, fracture aspect ratio, cross-section shape and other factors on the permeability of the fracture network. The analytical expression of the relative permeability of wet and non-wet phases was derived by Cai et al. [16], based on the fractal characteristics of pore size distribution and tortuosity in porous media. These models have been established based on fractal theory. Zhang et al. [17] combined the multi-scale finite element method with deep learning to simulate the flow of fluids in porous media based on a set of coarse grids and a set of fine grids. Zeng et al. [18] based on capillary bundles and tree-like network models, established apparent permeability models of gas in single capillaries and fractures by considering such factors as dynamic variation of fracture width, water saturation and actual gas effect, and provided the influence of matrix structural parameters on permeability. Although the above studies have made some progress, only one smooth feature is considered.
The fracture network is embedded into the matrix porous medium to form the socalled dual porous media, and dual porous media is widely present in nature. When exploiting natural gas resources or solving surface pollution problems, dual porosity and dual porous media are involved. Therefore, the transport characteristics of dual-porosity, dual porous media has aroused extensive and lasting research interest. Miao et al. [19] extended the fractal porous media theory to fracture networks and deduced the permeability research model of fracture-pore porous media and found that the maximum pore diameter and maximum trace length have a great influence on fluid seepage characteristics. Miao et al. [20] also considered the transfer flow behavior of fluids from matrix pores to fractured media, and proposed a dimensionless permeability model for porous media, finding that the ratio of fracture pressure difference to matrix pore pressure difference has an important impact on fluid flow from matrix pores to fractured media. Based on the topology model and complex network theory, Zhu et al. [21] proposed the power-law distribution of edges with node degree and established the permeability analytical formula of the dual porous medium, and found that the self-similarity index and the clustering coefficient of pore nodes play a decisive role in the permeability characteristics of the dual porous medium. By considering the channeling between the capillary bundle and the fracture, Wu et al. [22] established the effective permeability model of a power-law fluid in a fracture-pore type medium and found that the effective permeability of power-law fluid in a dual porous medium would be affected by the characteristics of the power-law fluid, structural parameters of fractures and pores, and pressure difference. Xiong et al. [23] represented rock fractures and pores as elliptical cross-section microtubules, calculated the permeability of a three-dimensional fracture network under the action of periodic pressure, and analyzed the variation rule of permeability. These are all studies of smooth porous media.
Since smooth surfaces are almost non-existent in nature, any surface (such as skin, desktop, road, glass, and tiny devices), regardless of its smoothness, is microscopically displayed as rough because of a large number of structural defects or geometric defects. Rough topography has an important influence on other transport mechanisms such as the flow resistance of fluids along rough surfaces, so it is necessary to characterize and describe the characteristics of rough surfaces accurately and simply. Qu [24] summarized the description methods of rough fracture structures, including raised height characterization method, joint roughness coefficient method, and fractal dimension method, and got the rough fracture surface characteristic description parameters such as tortuosity, roughness, and inclination. In practical research, fractal dimension, tortuosity, roughness, and inclination are generally used to characterize the rough surface of the rock stratum. Yang et al. [25] Materials 2022, 15, 4662 3 of 26 simplified the rough element on the surface of the microchannel into a cone and proposed the permeability fractal model of the rough capillary bundle, the effect of relative roughness on microchannel permeability was studied. Xiao et al. [26] derived the fractal model of fluid flow through a single rough curved capillary in fibrous porous media, elucidating the physical properties of fluid flow in a rough fiber porous medium. Liang et al. [27] calculated the immersion depth of the fluid in the curved capillary and used fractal theory to discuss the effect of relative roughness with immersion depth. He et al. [28] introduced tortuosity based on a fractal binary tree model, established a fractal network model of tortuosity fractal binary tree, and obtained the fluidity change law of gas in the shale reservoir fracture network. All the above models studied the permeability of the fluid in a single rough fracture but did not study in-depth the effect of roughness on fluid seepage in dual porous media.
Zhang et al. [29] wrote a calculation program based on the Boltzmann method to simulate the fluid flow in a single fracture with different JRC values and openings, and compared the simulated values with the results of the tortudic modified cubic law without considering roughness, and found that the tortudic modified cubic law predicted the flow with a certain deviation. Wang et al. [30] made 30 fractured rock samples with 10-level roughness (JRC value) and 3 different gap widths made by 3D printing technology into interpenetrating and filled fractured rock samples, and tested their permeability. Experiment data show that when the confining pressure is small, the greater the roughness, the greater the difference in fracture permeability of rock samples with different gap widths. Based on ten Barton standard lines, Li et al. [31] used COMSOL numerical simulation software to establish a rough single-fracture seepage model to detect the change of seepage velocity under different roughness. The results showed that the larger the roughness coefficient (JRC), the smaller the maximum and average flow velocity of fluid seepage. They described the effect of roughness on fluid seepage qualitatively, and considered that roughness is one of the important factors to be considered in the study of fluid seepage characteristics in fractures. According to the actual physical significance, roughness can reduce the actual fluid flow space, increase the flow resistance, and then reduce the permeability. Therefore, roughness is a factor that cannot be ignored.
Considering the roughness factor can result in a more comprehensive analysis of the permeability variation trend of porous media, and can also improve the accuracy of the model. At the same time, the roughness model can also be extended to other research fields, such as the development of shale gas, the study of thermal conductivity of porous media, and the development of tight oil and gas reservoirs. The literature [29][30][31] did not quantitatively describe the roughness. In order to quantitatively study the influence of roughness on permeability variation trend, based on the fractal theory and the basic formula of fluid mechanics, the permeability model of rough dual porous media is established. In Section 2, the fractal theory of fracture and pore size distribution in fractured dual porous media is presented. In Section 3, the rough surface of the fractured dual-porous media is characterized, and the specific expression of relative roughness is given. In Section 4, the permeability model of coarse fractured dual porous media is derived. In Section 5, the effectiveness of the model is verified by comparing the predicted values with experimental data, and the influence of structural parameters on the permeability of fractured porous media is studied. In Section 6, the conclusions of this study are given.

Fractal Theory of Fractured Dual Porous Medium
Many studies have shown that the pore size distribution of porous media meets the fractal scale law [10][11][12]32,33]. Yu et al. [33] gave the fractal scale relation of pore size distribution. Miao et al. [19] extend the fractal scale rate into the fracture network and give the fractal scale relationship of the distribution of fracture trace length. In this section, the fractal scale law is further extended to the fractured dual porous media, and the fractal relations satisfied by pore diameter and fracture trace length are given, respectively.
A fractured double porous medium is composed of a fracture network and porous matrix, because the representative elementary volume of the fracture network is much larger than that of the matrix pore [15,19], so the following research is based on the representative elementary volume of the fracture network. Figure 1 is the cross-section structure diagram of the representative elementary volume of the fracture network. A fractured double porous medium is composed of a fracture network and porous matrix, because the representative elementary volume of the fracture network is much larger than that of the matrix pore [15,19], so the following research is based on the representative elementary volume of the fracture network. Figure 1 is the cross-section structure diagram of the representative elementary volume of the fracture network. Due to the difference in geological structure, the diagenetic environment, and other factors, the fractures and matrix pores have different morphological characteristics, and their cross-sections may be: rectangular, polygonal, circular, ellipse, etc. For fractured dual porous media, to simplify the model, it is assumed that the cross-section of the fracture is rectangular and the shape of the fracture is represented by a cuboid. The crosssection of the matrix pores is circular, and its structure is characterized by a capillary bundle model. Figure 2 is the characterization diagram of the rough fracture, the fracture opening is a , the fracture trace length is l , and the fracture inclination angle is θ .  Due to the difference in geological structure, the diagenetic environment, and other factors, the fractures and matrix pores have different morphological characteristics, and their cross-sections may be: rectangular, polygonal, circular, ellipse, etc. For fractured dual porous media, to simplify the model, it is assumed that the cross-section of the fracture is rectangular and the shape of the fracture is represented by a cuboid. The cross-section of the matrix pores is circular, and its structure is characterized by a capillary bundle model. Figure 2 is the characterization diagram of the rough fracture, the fracture opening is a, the fracture trace length is l, and the fracture inclination angle is θ. Figure 3 shows the characterization diagram of the rough capillary, the diameter of the capillary is λ, the characteristic length is L 0 , and the actual length is L t . Additionally, the fractal theory of fractures and pores is established based on Figures 2 and 3, respectively. A fractured double porous medium is composed of a fracture network and porous matrix, because the representative elementary volume of the fracture network is much larger than that of the matrix pore [15,19], so the following research is based on the representative elementary volume of the fracture network. Figure 1 is the cross-section structure diagram of the representative elementary volume of the fracture network. Due to the difference in geological structure, the diagenetic environment, and other factors, the fractures and matrix pores have different morphological characteristics, and their cross-sections may be: rectangular, polygonal, circular, ellipse, etc. For fractured dual porous media, to simplify the model, it is assumed that the cross-section of the fracture is rectangular and the shape of the fracture is represented by a cuboid. The crosssection of the matrix pores is circular, and its structure is characterized by a capillary bundle model. Figure 2 is the characterization diagram of the rough fracture, the fracture opening is a , the fracture trace length is l , and the fracture inclination angle is θ .

Fractal Theory of Fracture Networks
Since the trace length of fractures satisfies the fractal scaling law [12,19], its power law expression can be written as [33]: where the negative sign indicates that the number of fractures decreases with the increase in fracture trace length, which is the same as the actual situation. The total number of fracture trace lengths between min l and max l is denoted as t N . According to Equation (3), there holds: where ( ) f l is the trace length of the fracture, that is: normalized by the probability density function, there holds:

Fractal Theory of Fracture Networks
Since the trace length of fractures satisfies the fractal scaling law [12,19], its power law expression can be written as [33]: (1) where N(l) represents the total number of fractures whose trace lengths are greater than or equal to l, and l min and l max represent the minimum and maximum lengths of fracture traces in the fracture network, respectively. D f represents the fractal dimension of fracture trace length in two dimensions, 0 < D f < 2, in three dimensions, 0 < D f < 3 . According to the fractal scaling relationship [13,19], Equation (1) can be further rewritten as: where k is the scale factor and an unknown quantity, l min ≤ l ≤ l max . Due to the large number of fractures in the fracture network, Equation (2) can be approximated as a continuously differentiable equation, then the number of fractures in [l, l + dl] is: where the negative sign indicates that the number of fractures decreases with the increase in fracture trace length, which is the same as the actual situation. The total number of fracture trace lengths between l min and l max is denoted as N t . According to Equation (3), there holds: where f (l) is the trace length of the fracture, that is: normalized by the probability density function, there holds: and in porous media there is l min /l max < 10 −2 , i.e., l min << l max , there holds: the expression of k is: substituting Equation (8) into Equation (5), the exact expression for the probability density function is: Yu et al. [13,34] gave the total number of pores when the size distribution of pore diameters satisfies the fractal scaling rate, Majumdar et al. [35] gave the total number of rough spots on the engineering surface, When studying the pores of porous media, it can be compared with the contact points on engineering surfaces, so the total number of fractures is: substituting Equation (10) into Equation (8), we obtain the exact value of the proportionality coefficient k: then substitute Equation (11) into Equation (3) to obtain the exact fractal scaling rate satisfied by the fracture trace length: Yu et al. [12,34] presented the relationship between porosity and fractal dimension when studying fractal scaling rate of porous media. Then, the fractal dimension of the fracture trace length can be defined as: represents the porosity of the fracture network, and refers to the ratio of the total area A p f of fracture pores to the total area A f of cross-section in the cross-section of the fracture network representative elementary volume. According to Equation (13), the expression of porosity φ f of the fracture network can be written as: at the same time, according to the physical meaning of fracture network porosity φ f , φ f can also be defined as [36]: where A f is the cross-sectional area of the representative elementary volume, A p f represents the total area of fractured pores in that area. According to Equation (12), the total area A p f of fracture pores on the cross-section of the representative elementary volume can be calculated as [36]: then substitute Equation (16) into Equation (15), and the cross-sectional area A f of the fracture network representative elementary volume is obtained as [36]: where β = a/l, l is the length of the fracture, a is the opening of the fracture. According to Equation (12), the total trace length of the fractures in the feature element body is: according to Equation (14), Equation (18) can be further simplified as: In a two-dimensional fracture network, the fracture areal density D is defined as [37]: substituting Equations (17) and (19) into Equation (20), the expression of fracture areal density D is [36]: Equation (21) shows that the areal density of the fracture network is a function of the maximum length l max of the fracture, and the proportional coefficient β, the fracture porosity φ f , and the fractal dimension D f of the fracture length.

Fractal Theory of Matrix Pores
Considering the matrix porous medium as a curved capillary whose pore diameter satisfies the fractal scaling law [13,34], the total number N(λ) of pores with diameters greater than or equal to λ meets the fractal scale rate: where λ min and λ max are the minimum and maximum values of pore diameter, respectively. D p is the fractal dimension of pore diameter, and its value range is the same as that of fracture length. The total number of pores is: due to the large number of matrix pores, Equation (23) can also be regarded as continuously differentiable functions. The number of pores in [λ, λ + dλ] is: where the negative sign indicates that the total number of pores decreases with the increase in pore diameter, which is the same as the actual situation. According to Equations (23) and (24), the probability density function of pore diameter is: According to the relationship between porosity and fractal dimension [12,34], the fractal dimension D p of the pore diameter can be defined as: where d E is Euclidean dimension, φ m is the porosity of the matrix pores, which refers to the ratio of the total pores area A p to the cross-sectional area A m in the cross-section of matrix representative elementary volume, according to Equation (26). The effective porosity of matrix porous media can be expressed as [12,34]: at the same time, according to the physical meaning of matrix porosity φ m , φ m can be expressed as [36]: According to Equation (24), the total pore area on the cross-section of matrix representative elementary volume can be calculated as: according to Equation (28), the cross-sectional area of matrix representative elementary volume is expressed as: The average diameter of pores according to Equations (24) and (25) is: The actual length L t (λ) and characteristic length L 0 of the capillaries satisfy the relation [12,34]: The tortuosity is defined as [38]: The average tortuosity is defined as [39]: The tortuosity fractal dimension D T can be defined as [34,40]: where τ is the average tortuosity of the capillary, and L 0 is the characteristic length of the capillary. λ is the average diameter of the pores, 1 < D T < 2 in two-dimensional space, 1 < D T < 3 in three-dimensional space. When D T = 1, the capillary is straight, the greater the D T , the greater the bending degree of the capillary.

Characterization of Rough Surfaces of Fractured Double Porous Media
Qu et al. [24] presented the description methods of rough surface: fractal dimension method, convex height characterization method, tortuosity, roughness and inclination angle. Yang et al. [25] used the convex height method to represent rough elements with cones. He et al. [28] use tortuosity to characterize rough surfaces, They both quantitatively studied the effect of rough surfaces on permeability. In this section, the cone-shaped roughness element model is introduced into the fractured dual porous media, and the roughness representation of fracture and pore surface is given, respectively, according to the distribution of roughness element on the fractured rock surface meets the fractal scale rate [41,42].

Characterization of Rough Surfaces of Fracture Networks
In the following research, Rough elements of rough fracture surfaces are approximated as cones [25,43,44], analogous to Equations (1) and (22), the distribution of the bottom diameter d of the cone follows the fractal scaling law [41,42]: Equation (36) is the fractal scaling law satisfied by the cumulative number N(d) of a rough element whose bottom diameter is greater than or equal to d. Due to a large number of rough elements, Equation (36) where D s is the fractal dimension of the bottom diameter of the cone, and the larger D s , the larger the bottom diameter of the rough element, in two dimensions, 1 < D s < 2, in three dimensions, 1 < D s < 3. By analogy with Equations (13) and (26), D s can be defined as [12,34]: where φ s denotes the percentage of the total bottom area S 1 of all cones to the total surface area S 0 of fractures in a fracture network representative elementary volume, and d E is the Euclidean dimension. According to Equation (38), φ s is [12,34]: At the same time, according to the physical meaning of φ s , φ s can also be defined as [25]: where S 1 is the total undersurface area of all cones in the representative elementary volume of the fracture network, and S 0 is the total area of the fracture surface. The ratio of the height of the cone to the diameter of the base is denoted as σ, that is: where d is the bottom diameter of the cone, and h is the height of the cone. The bottom area S i and volume V i of a small cone can be expressed as: According to Equation (37), the total volume V t and the total bottom area S 1 of all cones in the fracture network representative elementary volume are: According to Equations (40) and (45), the total area of the fracture surface is: According to the total volume expression by Equation (44) for all cones and the total area expression by Equation (36) for the fracture surface, the effective average height of the cone can be expressed as [25,[43][44][45]: Then, the relative roughness of the fracture surface is: In Equation (48) are functions of fracture opening a, fractal dimension D s of rough element bottom diameter, maximum value d max and minimum value d min of rough element bottom diameter, and proportional coefficient σ. When the effective average height of the small cone h = 0, i.e., σ = 0, the fracture surface is smooth, i.e., ε r = 0. This corresponds to the actual physical meaning.

Characterization of Rough Surfaces of Matrix Pores
Similar to the characterization of the rough elements of the fracture network, the surface rough masses of the curved capillary channel can also be approximated as cones [25,43,44], according to Equations (44) and (46). The effective average height of rough elements on rough capillaries can be expressed as: The relative roughness of the rough capillary bundle microchannel is [25,43,44]: Equation (50) shows that the greater the effective average height of the rough element, the greater the relative roughness of the capillary channel. When the effective average height of the rough element is zero, the capillary channel surface is smooth.

Fractal Model of Seepage in Rough Fractured Dual Porous Media
Miao et al. [19] established a fracture network permeability model considering the spatial orientation of fractures. The distribution of fractures in space is generally characterized by azimuth, dip angle, and strike [19,46]. The spatial distribution of fractures is shown in Figure 4 below [47], where AB is the fracture strike, CD is the maximum inclined direction, θ is the fracture inclination, and α is the fracture azimuth. Due to a large number of fractures in complex fracture networks and random distribution patterns in space, to simplify the study of fluid permeability in fractures, in this study, the average inclination angle and average azimuth angle of fractures are generally taken to represent the spatial orientation of the fracture network. It is assumed that the average inclination angle of fractures is θ and the average azimuth angle is α. Yang et al. [25] proposed a permeability model in a single horizontal rough pore by assuming that the fluid moves along the same direction. In this section, relative roughness, spatial orientation and tortuosity are introduced into the fractured dual porous media to obtain the seepage theory of fractured dual porous media.

Seepage Characteristics of Rough Fracture Networks
Firstly, the horizontal rough single fracture is studied. It is assumed that the fluid flows in the same direction, namely the X direction, and only the velocity component in the X direction is not zero, and the pressure gradient along the X direction remains unchanged. The spatial Cartesian coordinate system shown in Figure 5 is established for research [25]. The fracture is approximately regarded as a cuboid, and the width of the cuboid is the trace length l , the fracture opening is a .  Yang et al. [25] proposed a permeability model in a single horizontal rough pore by assuming that the fluid moves along the same direction. In this section, relative roughness, spatial orientation and tortuosity are introduced into the fractured dual porous media to obtain the seepage theory of fractured dual porous media.

Seepage Characteristics of Rough Fracture Networks
Firstly, the horizontal rough single fracture is studied. It is assumed that the fluid flows in the same direction, namely the X direction, and only the velocity component in the X direction is not zero, and the pressure gradient along the X direction remains unchanged. The spatial Cartesian coordinate system shown in Figure 5 is established for research [25]. The fracture is approximately regarded as a cuboid, and the width of the cuboid is the trace length l, the fracture opening is a. Yang et al. [25] proposed a permeability model in a single horizontal rough pore by assuming that the fluid moves along the same direction. In this section, relative roughness, spatial orientation and tortuosity are introduced into the fractured dual porous media to obtain the seepage theory of fractured dual porous media.

Seepage Characteristics of Rough Fracture Networks
Firstly, the horizontal rough single fracture is studied. It is assumed that the fluid flows in the same direction, namely the X direction, and only the velocity component in the X direction is not zero, and the pressure gradient along the X direction remains unchanged. The spatial Cartesian coordinate system shown in Figure 5 is established for research [25]. The fracture is approximately regarded as a cuboid, and the width of the cuboid is the trace length l , the fracture opening is a .  According to the Hagen-Poiseuille law, the fluid motion equation in fracture microchannel is [48]: where dp/dx is the pressure gradient along the flow direction X, u is the flow velocity along the X direction, and µ is the viscosity coefficient. Since the fracture opening is much smaller than the fracture trace length, namely a l, Equation (51) is simplified to: Assuming that the fracture microchannel is symmetric, the non-slip boundary condition is [25,48]: where h s is the effective average height of the roughness element on the fracture surface, u R is the flow velocity of the fluid in the fracture microchannel, and the solution is: where dp/dx is the absolute value of the pressure gradient along the flow direction. Equation (54) shows that the flow velocity in the fracture network is a function of the absolute value of the pressure gradient dp/dx, the effective average height of the rough surface h s , and the distance |z| of the fluid from the center of the pipe wall. At the same time, it can be found that the closer the distance |z| of the fluid to the center of the pipe, the greater the flow velocity of the fluid, and the greater the effective average height of the roughness element, the smaller the flow velocity of the fluid. This is consistent with objective facts. The average velocity of the fluid through the rough fracture microchannel is: If h s = 0, then Equation (55) is simplified to: This is the average velocity of the fluid in the smooth fracture network. According to Equation (56), through the horizontal rough single fracture, the volume of fracture trace length l is: Secondly, considering the spatial orientation of the fracture [19,47], the flow rate in the rough single fracture in the representative elementary volume of the fracture network is: Equation (58) shows that the volume flow rate in a single rough fracture is the function of effective opening h s of the fracture surface rough element, pressure difference p, average inclination angle θ, azimuth angle α, and fracture opening a, in the representative elementary volume of the fracture network. At the same time, it can be found that when other quantities are constant, the greater the effective opening of the rough element, the smaller the volume flow in the microchannel of a single rough fracture. In addition, when the volume flow in a single rough fracture is constant, the greater the effective opening of the rough element, the greater the pressure difference between the two ends of the fracture, and the greater the flow resistance of the fluid. According to Equation (12), the total flow of fluids through the rough fractures network representative elementary volume as: where ε r is determined by Equation (47). Equation (59) shows that the total flow in the fracture network is the function of average inclination angle θ and azimuth angle α, pressure difference p, the relative roughness ε r of fracture, the fracture length fractal dimension D f , and the maximum length l max of the fracture. When other quantities are constant, the larger the relative roughness of the fracture, the smaller the flow of the fracture network, which is in line with objective physical facts. According to Darcy's Law: The permeability of the fracture network is: Substituting the A f value in Equation (17) into Equation (61), the permeability of the fracture network can be simplified as: By substituting Equation (21) into Equation (61), the permeability of the fracture network can be expressed as: Equation (63) shows that the permeability of the rough fracture network is a function of the structural parameters of the fracture medium (fracture length fractal dimension D f , fracture area density D, porosity φ f , relative roughness ε r , maximum fracture trace length l max , average fracture inclination angle θ and azimuth angle α, proportional coefficient β), and is negatively correlated with the relative roughness ε r of the fracture network. When ε r = 0, Equation (63) can be simplified as: This is the permeability when the microchannel wall of the fracture network is smooth.

Seepage Characteristics of Pores in Rough Matrix
Firstly, the seepage characteristics of the fluid in the rough straight capillary are studied. Assuming that the fluid flows in the X direction, only the velocity component in the X direction is not zero, and the pressure gradient along the X direction is constant, and the spatial rectangular coordinate system shown in Figure 6 is established [25]. The diameter of the capillary is λ and the radius is r 0 .
This is the permeability when the microchannel wall of the fracture network is smooth.

Seepage Characteristics of Pores in Rough Matrix
Firstly, the seepage characteristics of the fluid in the rough straight capillary are studied. Assuming that the fluid flows in the X direction, only the velocity component in the X direction is not zero, and the pressure gradient along the X direction is constant, and the spatial rectangular coordinate system shown in Figure 6 is established [25]. The diameter of the capillary is λ and the radius is 0 r .
The no-slip boundary condition for flow is [25]: Analogous to the calculation of fluid flow rate in rough fracture, the equation of motion of fluids in rough straight capillaries in fracture microchannels is [25]: The no-slip boundary condition for flow is [25]: where h m is the effective average height of roughness element on matrix pore surface, and u m is the flow rate of the fluid in pore microchannel. The flow rate of fluids in the rough capillary microchannel is: The average flow rate is: The flow rate through a single rough straight capillary in a matrix representative elementary volume is: Second, considering the tortuosity characteristics of rough capillary, according to Equation (23) and Equation (32), the total flow of fluids through the matrix representative elementary volume is: Equation (70) is the function of fractal dimension of tortuosity, the fractal dimension of pore diameter, pressure difference, the maximum diameter of the pore, and the relative roughness of matrix pores. According to the definition of Darcy's law and Equation (70), the permeability of rough matrices pores can be expressed as: Equation (71) shows that the permeability of rough matrix pores is the function of structure parameters of matrix pore (the maximum diameter λ max of the pore, fractal dimension D p of the pore diameter, fractal dimension D T of the tortuosity, the relative roughness ε c of the pore). When other quantities are constant, the permeability of the rough matrix decreases with the increase in the relative roughness of the pores. When ε c = 0, Equation (71) can be simplified to: This is the permeability when the pore wall of the matrix is smooth.

Permeability Model of Rough Fractured Double Porous Media
Because the characterization unit of the fracture network is larger than that of matrix pore and contains multiple matrix units [15,19], the unit body of the dual porous media is selected as the unit body of the fracture network, there holds: That is: Since the flow into the porous media is equal to the flow out of the porous media [49], the total flow of fluids in the fractured dual porous media is equal to the sum of the flow in the matrix and the fracture network: The permeability of Newtonian fluids in rough fractured double porous media is defined by Darcy's law, as: Equation (76) is a function of fracture network structure parameters (cross-sectional area A f of the fracture network characteristic unit body, fracture average inclination θ and azimuth α, fractal dimension D f of fracture length, relative roughness ε r of fracture, and proportional coefficient β) and matrix pore structure parameters (fractal dimension D p of pore diameter, tortuosity fractal dimension D T , maximum diameter λ max of pore, cross-sectional area A m of matrix pore characteristic unit body, and relative roughness ε c of pore). Meanwhile, the increase in the relative roughness of the fracture network and the relative roughness of the matrix pores lead to the decrease in permeability of fractured dual porous media.
To study the contribution of porous matrix and the fracture network to fluid permeability, dimensionless permeability K + is defined, namely [19,20]: where K m and K f represent the permeability of matrix pores and fracture networks, respectively.

Results and Discussion
Firstly, the reliability of the rough fracture network model is verified. Figure 7 compares the predicted value of the permeability model of the rough fracture network with the numerical simulation value in the literature [50] under different relative roughness. The authors of [50] selected 22 fracture types with different scales in southwestern Turkey, and the digital fracture model was imported into commercial modeling software, and the equivalent fracture network permeability was calculated by constructing a 3 D model with a network size of 100 m × 100 m × 10 m. Among them, the largest fracture length is 1 m, and the fracture inclination angle θ = 0. In the calculation, take β = 0.002, and calculate the predicted value of the model according to Equation (61).
with the numerical simulation value in the literature [50] under different relative roughness. The authors of [50]  this is because the larger the roughness, the smaller the volume of the fracture pores, the space for fluid to flow is reduced, and the flow resistance is increased. Through the comparison of Figure 7, it is found that the established rough fracture network model is effective, and the influence of roughness on permeability is also very large.  According to Figure 7, the permeability of the rough fracture network increases with the increase in fractal dimension D f of fracture trace length, because the increase in fractal dimension of fracture trace length leads to the increase in fracture cross-sectional area. In addition, when the relative roughness ε r of the fracture network is 0 and 0.1, the predicted permeability of the fracture network calculated by Equation (61) is 78.23% and 29.93% higher than that of the experimental simulation value. When the relative roughness ε r of the fracture network is 0.3 and 0.5, respectively, the predicted permeability of the fracture network calculated by Equation (61) is 38.87% and 77.72% lower than the experimental simulation value. This shows that under the same fracture fractal dimension D f , the larger the relative roughness ε r , the smaller the fracture permeability, this is because the larger the roughness, the smaller the volume of the fracture pores, the space for fluid to flow is reduced, and the flow resistance is increased. Through the comparison of Figure 7, it is found that the established rough fracture network model is effective, and the influence of roughness on permeability is also very large.
The influence of fracture geometry parameters on the permeability of rough fracture model is studied. Figure 8 shows the relationship between fracture network permeability and porosity under different roughness and fracture network inclination. In the calculation, β = 0.01 , fracture maximum length l max = 10 mm, and fracture spatial azimuthal angle α = 0, according to Equations (13) and (62), the corresponding rough fracture network permeability is calculated. According to Figure 8, the permeability of the rough fracture network increases with the increase in fracture porosity, and the greater the inclination and relative roughness of the fracture, the smaller the permeability of the rough fracture. This is due to the increase in fracture porosity, which will increase the volume of fracture pores and increase the flowable volume of fluid, so that the fluid flow resistance decreases. The increase in inclination angle will increase the flow resistance and the increase in fracture relative roughness will reduce the effective opening of the fracture. The influence of fracture geometry parameters on the permeability of rough fracture model is studied. Figure 8 shows the relationship between fracture network permeability and porosity under different roughness and fracture network inclination. In the calculation, 0.01 β = , fracture maximum length max 10 mm l = , and fracture spatial azimuthal angle =0 α , according to Equations (13) and (62), the corresponding rough fracture network permeability is calculated. According to Figure 8, the permeability of the rough fracture network increases with the increase in fracture porosity, and the greater the inclination and relative roughness of the fracture, the smaller the permeability of the rough fracture. This is due to the increase in fracture porosity, which will increase the volume of fracture pores and increase the flowable volume of fluid, so that the fluid flow resistance decreases. The increase in inclination angle will increase the flow resistance and the increase in fracture relative roughness will reduce the effective opening of the fracture.  β . According to Equation (63), the corresponding rough fracture network permeability is calculated. According to Figure 9, the permeability of the rough fracture network increases with the increase in the fracture surface density. The greater the relative roughness of the rough fracture network, the smaller the permeability, and the effect of relative roughness on the fracture permeability is very obvious. This is due to the increase in fracture surface density will increase the length of fracture trace and reduce the fluid flow resistance. The increase in the relative roughness of the fracture will increase the flow resistance of the fluid.  Figure 9 shows the relationship between fracture network permeability and fracture surface density under different relative roughness, in the calculation, the maximum fracture length l max = 10 mm , the fracture dip angle θ = 45 • , the average porosity of the fracture network is 0.018, the average fractal dimension of the fracture network is 1.3, and β = 0.01. According to Equation (63), the corresponding rough fracture network permeability is calculated. According to Figure 9, the permeability of the rough fracture network increases with the increase in the fracture surface density. The greater the relative roughness of the rough fracture network, the smaller the permeability, and the effect of relative roughness on the fracture permeability is very obvious. This is due to the increase in fracture surface density will increase the length of fracture trace and reduce the fluid flow resistance. The increase in the relative roughness of the fracture will increase the flow resistance of the fluid. . Relationship between the permeability of rough fracture network and fracture surface density and relative roughness. Figure 10 is a comparison of the predicted pore permeability of the rough matrix with the numerical simulation values of smooth matrix permeability in [51]. The authors of [51] measured the permeability and other data of 30 natural mudstones collected from oil wells at depths of two to three km, among them, the permeability was measured on disks with a thickness of 0.005-0.008 m and a diameter of 0.0254 m. The porosity is between 0.06 and 0.27, the average radius of the pore is between 2.8 and 1403.4 nm, and no fractures are found. So the maximum pore diameter is 2806.8 nm and the minimum is 5.6 nm. Equation (26) can be used to calculate p D , Equations (31), (34) and (35) can be used to calculate T D , and Equation (71) can be used to calculate the corresponding matrix permeability prediction value. Figure 10 shows that when the relative roughness of matrix pores is 0.1, 0.2, and 0.4, the predicted values of matrix pore permeability calculated by Equation (71) are 66.4%, 79.05%, and 93.37% lower than the experimental simulation data, respectively. Through comparison, it is found that the predicted value of the rough pore permeability model is smaller than the experimental data, and the greater the relative roughness of the pore, the smaller the pore permeability, which further illustrates the effectiveness of this model. Figure 9. Relationship between the permeability of rough fracture network and fracture surface density and relative roughness. Figure 10 is a comparison of the predicted pore permeability of the rough matrix with the numerical simulation values of smooth matrix permeability in [51]. The authors of [51] measured the permeability and other data of 30 natural mudstones collected from oil wells at depths of two to three km, among them, the permeability was measured on disks with a thickness of 0.005-0.008 m and a diameter of 0.0254 m. The porosity is between 0.06 and 0.27, the average radius of the pore is between 2.8 and 1403.4 nm, and no fractures are found. So the maximum pore diameter is 2806.8 nm and the minimum is 5.6 nm. Equation (26) can be used to calculate D p , Equations (31), (34) and (35) can be used to calculate D T , and Equation (71) can be used to calculate the corresponding matrix permeability prediction value. Figure 10 shows that when the relative roughness of matrix pores is 0.1, 0.2, and 0.4, the predicted values of matrix pore permeability calculated by Equation (71) are 66.4%, 79.05%, and 93.37% lower than the experimental simulation data, respectively. Through comparison, it is found that the predicted value of the rough pore permeability model is smaller than the experimental data, and the greater the relative roughness of the pore, the smaller the pore permeability, which further illustrates the effectiveness of this model. Figure 11 shows the relationship between the total flow and the pressure difference for a Newtonian fluid, when the temperature is 25 • C and the confining pressure is 500 Kpa, the total flow prediction of a Newtonian fluid through the fractured double porous media, and the flow prediction through the fracture network and the matrix pore are compared with the experimental data [52], which simulated the smooth dual porous media. In the calculation, the average inclination angle of the fracture is θ = 0 • , the minimum diameter of the matrix capillary is 2 nm, the viscosity coefficient of the water is µ = 1.12 × 10 −3 Pas, and the maximum length of the fracture is 10 mm. The relationship between the maximum and minimum pore diameter and the fracture length is λ min /λ max = 0.001, and l min /l max = 0.001. According to the literature [52], the matrix porosity is φ m = 0.25, and the fracture porosity is 1/100 of the matrix, that is φ f = 0.0025. Then, the corresponding flow rates of the fracture network, matrix pores, and fractured dual porous media are calculated according to Equations (59), (70), and (75).   (70) is 7.94% and 32.98% lower than the experimental simulation value, respectively, and it can be found that the flow of the fracture network is almost the same as that of the fractured dual porous media, and the flow of the matrix pore is much smaller than the flow of the fracture network, and the relative roughness of the matrix pore has little effect on the matrix flow. The influence of structural parameters on permeability of fractured dual porous media is studied below. When the structural parameters ( f L ) of the fractured dual porous media are known, the permeability prediction can be calculated using Equation (76). The structural parameters in the calculation process take Figure 11. Comparison of predicted and experimental data [52] of Newtonian fluid flow at 25 • C temperature and 500 kPa confining pressure.
According to Figure 11, when the relative roughness of the matrix pore and the fracture network are ε r = ε c = 0 and ε r = ε c = 0.1, the total flow prediction value of the fractured dual porous media calculated by Equation (75) is 7.53% and 32.62% lower than the experimental simulation values, respectively, This indicates that the flow rate of rough fractal model is lower than the total flow of smooth fractal model and the experimental data, and further shows that the roughness of the porous medium has a greater influence on fluid flow. This is because when the relative roughness increases, the flow resistance of the fluid increases, resulting in a decrease in the total flow through it, this is the same as the actual physical situation and further indicates the validity of the established model. At the same time, when the relative roughness of the fracture network is 0, 0.1, the predicted value of the total fracture network flow calculated by Equation (70) is 7.94% and 32.98% lower than the experimental simulation value, respectively, and it can be found that the flow of the fracture network is almost the same as that of the fractured dual porous media, and the flow of the matrix pore is much smaller than the flow of the fracture network, and the relative roughness of the matrix pore has little effect on the matrix flow.
The influence of structural parameters on permeability of fractured dual porous media is studied below. When the structural parameters (φ f , φ m , D p , D T , β, α, θ, L 0 ) of the fractured dual porous media are known, the permeability prediction can be calculated using Equation (76). The structural parameters in the calculation process take the values in the literature [19], where the maximum diameter of the matrix pores is λ max = 2 µm , the average inclination angle of the fracture network is θ = π/6, the maximum value of the fracture length is l max = 0.010 m, the average azimuth angle of the fracture network is α = 0, the linear length of the matrix pore is L 0 = 0.1 m, and the ratio of the minimum and the maximum pore diameters λ min /λ max and the ratio of the minimum and maximum fracture length l min / l max are all 0.001.
Firstly, when the structural parameters of matrix porous medium are determined, the influence of fracture network structural parameters on the permeability of fractured dual porous media is studied. Figure 12 shows the change of permeability K of fractured dual porous media when the relative roughness ε r and ε c of the fracture network and matrix pores and the fractal dimension D f of fracture length take different values. In the calculation, take β = 0.01, φ m = 0.25. According to Equations (14), (26) and (35), the fractal dimension D p of the pore diameter, the porosity φ f of the fracture network, and the fractal dimension D T of the tortuosity can be determined, and the cross-sectional area A f , A m of the fracture and matrix characterization element can be determined by Equations (17) and (30). Then, the corresponding permeability of fractured dual porous media can be determined according to Equation (76). According to Figure 12, the permeability of fractured dual porous media increases with the increase in fractal dimension of fracture trace length. When the fractal dimension of fracture trace length is constant, the permeability of dual porous media decreases with the increase in the relative roughness of fracture and matrix. At the same time, it can also be found that when the values of ε r , ε c are 0.3, 0.1 and 0.3, 0.3, respectively, the permeability of the fractured dual porous media is almost equal, indicating that the relative roughness of the matrix pores has little effect on the seepage of the fluid. When ε r , ε c are 0.1 and 0.3, respectively, the permeability of fractured dual porous media is far greater than that of ε r , ε c is 0.3 and 0.3, respectively, and is almost equal to that of ε r , ε c is 0.1 and 0.1, respectively. This indicates that the relative roughness of the fracture network has a very large influence on the seepage of the fluid, which further indicates that the main function of matrix pores is storage, and the fracture network can be the main channel of fluid seepage. Thus, when studying the permeability of rough fractured dual porous media, the influence of matrix roughness can be ignored, considering that the matrix pores are smooth. and is almost equal to that of , r c ε ε is 0.1 and 0.1, respectively. This indicates that the relative roughness of the fracture network has a very large influence on the seepage of the fluid, which further indicates that the main function of matrix pores is storage, and the fracture network can be the main channel of fluid seepage. Thus, when studying the permeability of rough fractured dual porous media, the influence of matrix roughness can be ignored, considering that the matrix pores are smooth. Figure 12. Relationship between permeability K of fractured dual porous media and fracture length fractal dimension D f with different relative roughness. Figure 13 shows the variation law of the permeability of the fractured dual porous media when the ratio β of fracture opening to trace length, the relative roughness ε r of the fracture network, and the average inclination θ change. In the calculation, take φ m = 0.25, φ f = 0.0025. According to Figure 13, the permeability K of fractured dual porous media increases with the increase in the ratio β of fracture opening to the fracture length, and the greater the β, the faster the increase in K . At the same time, the greater the average inclination angle θ and relative roughness ε r of the fracture network, the smaller the permeability K of the fractured dual porous media, and the slower its increase rate with the increase in β. This indicates that the greater the average inclination angle and the relative roughness of the fracture network, the greater the flow resistance of the fluid.
Secondly, when the structural parameters of the fracture parameters are determined, the influence of matrix pore structural parameters on fluid permeability is studied. According to Equations (76) and (77), when the structural parameters of the fracture network are determined, the dimensionless permeability K + is only affected by matrix pore structural parameters. Figure 14 shows the changing trend of dimensionless permeability K + with matrix porosity φ m , fracture network and matrix relative roughness ε r , ε c . In the calculation, the structural parameters of the fracture network are taken from the values in [19]. As can be seen from Figure 14, when the relative roughness ε r , ε c of the fracture network and the matrix are constant, K + increases with the increase in the matrix porosity φ m , and when relative roughness ε r , ε c increases at the same time, the dimensionless permeability K + decreases, indicating that the decreases rate of matrix permeability is less than that of fracture network permeability. This further shows that relative roughness has a greater impact on fracture networks. In addition, the permeability of matrix pores is much smaller than that of the fracture network, no matter how the relative roughness changes, in other words, the main function of the matrix pores is storage, and fracture networks serve as the main channels for fluid flow. Additionally, because the fluid has a dissolution effect, the pores of the porous medium will gradually become large caves, and the fractures will also become large fractures. ture length, and the greater the , the faster the increase in . At the same time, the greater the average inclination angle θ and relative roughness r ε of the fracture network, the smaller the permeability K of the fractured dual porous media, and the slower its increase rate with the increase in β . This indicates that the greater the average inclination angle and the relative roughness of the fracture network, the greater the flow resistance of the fluid. Figure 13. Relationship between the permeability K of the fractured dual porous media and the ratio β of fracture opening to trace length, the relative roughness r ε of the fracture network, and the average inclination angle θ .
Secondly, when the structural parameters of the fracture parameters are determined, the influence of matrix pore structural parameters on fluid permeability is studied. According to Equations (76) and (77), when the structural parameters of the fracture network are determined, the dimensionless permeability K + is only affected by matrix pore structural parameters. Figure 13. Relationship between the permeability K of the fractured dual porous media and the ratio β of fracture opening to trace length, the relative roughness ε r of the fracture network, and the average inclination angle θ.
Materials 2022, 15,4662 25 of 29 Figure 14 shows the changing trend of dimensionless permeability K + with matrix porosity m φ , fracture network and matrix relative roughness r ε , c ε . In the calculation, the structural parameters of the fracture network are taken from the values in [19]. As can be seen from Figure 14, when the relative roughness r ε , c ε of the fracture network and the matrix are constant, K + increases with the increase in the matrix porosity m φ , and when relative roughness r ε , c ε increases at the same time, the dimensionless permeability K + decreases, indicating that the decreases rate of matrix permeability is less than that of fracture network permeability. This further shows that relative roughness has a greater impact on fracture networks. In addition, the permeability of matrix pores is much smaller than that of the fracture network, no matter how the relative roughness changes, in other words, the main function of the matrix pores is storage, and fracture networks serve as the main channels for fluid flow. Additionally, because the fluid has a dissolution effect, the pores of the porous medium will gradually become large caves, and the fractures will also become large fractures.

Conclusions
Based on the hypothesis that the fracture network is composed of rectangular pipes Figure 14. Relationship between dimensionless permeability K + and matrix porosity φ m , fracture network and matrix relative roughness ε r , ε c .

Conclusions
Based on the hypothesis that the fracture network is composed of rectangular pipes of different sizes and satisfies the fractal distribution, the matrix pores are composed of curved capillaries with circular cross-sections, and the roughness is characterized by small cones. The fractal theory is used to deduce the total flow rate and permeability model of a Newtonian fluid in rough fractured dual porous media, and the relationship between the permeability of dual porous media and the fractal dimension D f , D p , D T , relative roughness (ε r , ε c ), fracture network inclination θ, fracture porosity φ f , matrix porosity φ m and other parameters of the fracture network and matrix pores is obtained. By comparing with the experimental data, it is found that the total flow rate of rough fractured dual porous media decreases with the increase in relative roughness and is lower than the experimental data, indicating the reliability of the established model. Through comparison, it is found that the permeability of the fractured dual porous media increases with the increase in the fracture porosity, the ratio of fracture opening to trace length, and the porosity of the matrix; and the greater the relative roughness of the fracture and matrix and fracture inclination angle, the smaller the permeability of the fractured dual porous media. In addition, the relative roughness of the fracture network has a much greater influence on the permeability of fractured dual porous media than the relative roughness of matrix pores; therefore, the roughness of the matrix pores can be ignored in this study, assuming that it is smooth.  Acknowledgments: The authors would like to thank the associate editors and anonymous reviewers for their numerous constructive comments that will improve the content of this paper.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature a
The opening of the fracture l The trace length of the fracture The total flow of fluids through the rough fractures network representative elementary volume Q m The total flow of fluids through the matrix representative elementary volume K f The permeability of the fracture network K m The permeability of rough matrices pores L 0 The characteristic length of the capillary L t The actual length of the capillary d The bottom diameter of the cone h The height of the cone σ The ratio of the height of the cone to the diameter of the base D T The tortuosity fractal dimension D s The fractal dimension of the bottom diameter of the cone S 1 The total bottom area of all cones in representative elementary volume S 0 The total area of fracture surface in representative elementary volume A p f The total area of fracture pores A p The total matrix pores area A f The cross-section area of the representative elementary volume of the fracture network A m The cross-sectional area of matrix representative elementary volume S i The bottom area of a small cone V t The total volume of all cones in the representative elementary volume V i The volume of a small cone D The fracture areal density Q The total flow of fluids in the fractured dual porous media K The permeability of a Newtonian fluid in rough fractured double porous media K + The dimensionless permeability