A Unified View of Nonlinear Resistance Formulas for Seepage Flow in Coarse Granular Media

Abstract

There are many studies on the nonlinear relationship between seepage velocity and hydraulic gradient in coarse granular materials, using different approaches and variables to define the resistance formula applicable to that type of granular media. On the basis of an analysis of the existing formulations developed in different studies, we propose an approach for comparing the results obtained by some of the most important studies on state-of-the-art seepage flow in coarse granular media.

1. Introduction and Objectives

It is essential to know the relationship between filtration speed and hydraulic gradient to understand the interactions that occurs between infiltrated water and structures composed of gravel or rockfill, such as dams, levees, drainage structures, and coastal dikes. These porous media, made up of coarse particles, possess special characteristics due to their large pores which, under certain conditions, may give rise to non-laminar seepage flow, invalidating Darcy’s Law. Calculating the seepage flow through this type of porous media requires the use of so-called nonlinear resistance formulas. The flow regime in which this is applied is known as non-Darcy flow.

Various nonlinear relationships have been proposed to describe the flow in coarse porous media. They can be grouped into two types of equations as follows:

$$i=a·V^b$$

(1)

$$i=r·V+s·V^2$$

(2)

where V is the seepage velocity, defined as the average fluid velocity in the whole transversal section; $i$ is the hydraulic gradient; $a$ and $r$ are parameters depending on the characteristics of the porous medium and of the fluid; $s$ is only a function of the characteristics of the porous medium; and $b$ is a function parameter of the conditions of the flow.

Equation (1) is the so-called exponential equation and Equation (2) is the quadratic equation. To obtain relationships among the parameters $a$, $b$, $r$, and $s$, different physical parameters have been considered; experimental data with different intervals of size, shape, and particle angularity have been used, as well as a wide range of gradient intervals. However, there is currently no formula that can completely create Equations (1) and (2) and combine the determination criteria of their coefficients based on physical parameters.

In order to have a unified vision of the relationships between gradients and seepage velocity developed in different studies, we conduct an analysis of the existing formulations and identify similarities by mathematically comparing the physical parameters considered in these studies. These relationships make it easier to compare the results and formulations obtained by each of the studies considered in this article.

2. Review of Resistance Formulas in Nonlinear Porous Media

2.1. Conceptual Approach

Various studies have been completed with the aim of developing expressions for parameters $r$ (linear coefficient) and $s$ (quadratic coefficient) of Forchheimer’s Law (1901) [1] Equation (2), which defines a macroscopic hydraulic behaviour. The first term of Equation (2) (i.e., the term with coefficient $r$) represents the loss of energy due to the viscous forces and depends on the properties of both the porous medium and the fluid. The second term of Equation (2) (i.e., the term with coefficient s) considers the loss of energy due to the forces of inertia, and depends only on the properties of the porous medium.

In most of these studies, the development of formulas was based on the analogy of the flow in pipes, through the application of two dimensionless groups that, in this paper, are referred to as the generalised friction factor $f$ of the Darcy–Weisbach Equation (3) and generalised Reynolds number $R_e$ Equation (4) represented as:

$$f=L_c·2g·\frac{i}{V_p^2}$$

(3)

$$R_e=\frac{L_c·V_p}{v}$$

(4)

where $L_c$ is the characteristic length adopted in each case, $g$ is the gravitational acceleration, $i$ is the hydraulic gradient, $v$ is the kinematic viscosity, and $V_p$ is the pore velocity, determined by Equation (5):

$$V_p=\frac{V}{n}$$

(5)

being $n$ the porosity of the porous medium.

The characteristic lengths, on which most studies have been based, can be grouped into three types:

(a)

the representative size of the particle ($D$)

(b)

the square root of the intrinsic permeability ($K_0$), as a macroscopic property of the porous medium. For the laminar regime, it is determined by Equation (6):

$$K_0=\frac{v}{g}·\frac{1}{i}·V$$

(6)

(c)

the hydraulic mean radius $R_h$ that was first defined by Taylor (1948) [2] and determined by Equation (7):

$$R_h=\frac{n}{S_e·\left(1-n\right)}$$

(7)

where $S_e$ is the average specific surface area of the solid particles that make up the porous medium and depends on the shape, angularity, and surface roughness of the particles (Crawford, C.W. et al., 1986 [3]; Sabin G.C.W. and Hansen D., 1994 [4]) and is defined by the Equation (8):

$$S_e=\frac{SP}{VP}$$

(8)

where SP is the average particle surface area and VP is the average particle volume.

Among the first studies are those developed by Blake (1922) [5] and S.P. Burke and W.B. Plumer (1928) [6], who used dimensionless groups defined by Equations (3) and (4), adopting $L_c$ as the characteristic length, and $S$ as the specific surface area of the packing of porous medium. The relationship with the average particle specific surface area $S_e$ determined by Equation (9) as follows:

$$S=S_e·\left(1-n\right)$$

(9)

In their tests, they used commercial porous media (glass beads, glass rings, solid glass cylinders, and Raschig rings), and the experimental data corresponding to a porous medium with the same geometry and size. The results were represented a generalised diagram [$R_e,f$, Equation (10) and the values fitted to a smooth curve:

$$f=\frac{A_1}{R_e}+A_2$$

(10)

where $A_1$ is referred to as a linear generalised dimensionless coefficient, and $A_2$ as a quadratic generalised dimensionless coefficient.

Equation (10) is referred to as a generalised equation $\left[R_e,f\right]$.

Figure 1 shows a [$R_e,f$], schematic diagram similar to that used by Bear (1988) [7] with asymptotic values for increased Reynolds numbers $R_e$.

As indicated by Sabri Ergun and A. A. Orning (1949) [8] “This transition from the dominance of viscous to kinetic effects, for most packed systems, is smooth, indicating that there should be a continuous function relating pressure drop to flow rate.”

Ward (1964) [9] stated the same thing when he asserted, “The smooth transition from laminar to turbulent flow in porous media is expected. In the laminar flow region, the flow is laminar in all parts of the porous media. In the laminar transition region, the flow is laminar in most parts of the porous medium, but there are parts where the flow is turbulent. In the turbulent transition region, the flow is turbulent in most parts of the porous medium, but there are still parts where laminar flow conditions persist. Finally, turbulent flow exists in all parts of the porous medium at high values of $R_k$. Simultaneous existence of laminar and turbulent flow in different parts of a porous medium is possible because of the irregularities and variation in pore size.”

Dudgeon (1966) [10] indicated that “the only likely solution to the problem of a generalised $\left[R_e,f\right]$ plot is in terms of a set of graphs for each family of geometrically similar porous media.”

Other authors (Blake (1922) [5] and S.P. Burke and W.B. Plumer (1928) [6]; Morcon A.R. (1946) [11]; Ergun (1952) [12], Kadlec, H.R., and Knight, L.R. (1966) [13]; Ahmed and Sunada (1969) [14]; Kovacs (1969) [15], Arbhabhirama and Dinoy (1973) [16]; Stephenson (1979) [17]; Li, B et al. (1998) [18];and more recently Sidiropuolou et al. (2007) [19]; Moutsopoulos et al. (2009) [20]; Sedghi-Asl and Rahimi (2013) [21]; and Salahi et al. (2015) [22]), obtained continuous curves such as those given by Equation (10) through the corresponding adjustment of the experimental data used, always within the range of the Reynolds number $R_e$ on which the tests were developed.

The existence of these continuous curves in porous media contrasts with the turbulent flow in pipes where there are sudden jumps in the Reynolds number interval $R_e$ between 2000 and 4000 (White, F.W. (2003) [23]).

If we substitute the generalised values of $f$ and $R_e$, given by Equations (3) and (4) in Equation (10) we get Equation (11):

$$L_c·2·g·\frac{i}{V_p^2}=\frac{A_1·v}{L_c·V_p}+A_2$$

(11)

Resolving $i$ in Equation (11) we get what we refer to as a quadratic generalised equation, i.e., Equation (12):

$$i=\frac{v}{2·g}·A_1·\frac{1}{L_c^2}·V_p+\frac{1}{2·g}·A_2·\frac{1}{L_c}·V_p^2$$

(12)

In accordance with the Equation (12) the parameters $r$ and $s$ of the Forchheimer equation, i.e., Equation (2) determine Equations (13) and (14):

$$r=\frac{v}{2·g}·A_1·\frac{1}{L_c^2}·\frac{1}{n}$$

(13)

$$s=\frac{1}{2·g}·A_2·\frac{1}{L_c}·\frac{1}{n^2}$$

(14)

In the case of the fully developed turbulent regime, it may be possible to disregard the linear expression $r$, in such a way that we could obtain the exponential equation, Equation (2), considering the coefficient a as equal to the quadratic expression $s$ and the exponent $b$ as equal to 2.

Various researchers have worked on the exponential law, Equation (1), in regimes of transition, and therefore with the exponent values b below 2: Wilkins (1956) [24], b = 1.85; Dudgeon (1966) [10], 1.2 < b < 1.91; Parkin (1991) [25], b = 1.85; Moutsopoulos K. N. et al. (2009) [20], 1.280 < b < 1.687; Sedhi-Asl et al. (2013) [21] 1.479 < b < 1.804.

In this respect, Stephenson (1979) [17] pointed out the fact that not obtaining an adjustment of $b$ equal to 2, with the experimental data, is due to “the tests being carried out on a small scale and as a result of low Reynolds numbers.” Ferdos, F. et al. (2015) [26] obtained an exponent $b$ equal to 2 in their tests as a result of using elevated Reynolds numbers $R_e$. These studies both used the particle Reynolds number $R_d$ determined by Equation (15) as follows:

$$R_d=\frac{D·V_p}{v}$$

(15)

They achieved values of $R_d$ = 220,000 for the size interval of 100–160 mm and $R_d$ = 320,000 for the range of 160–240 mm. The values were both very much higher as compared with the value proposed by Stephenson (i.e., R_d = 10,000) in order to reach the fully developed turbulent flow.

It is important to consider that the flow through for highly permeable porous media may contain various flow regimes, as the generalised Reynolds number $R_e$ increases. In general, most of the studies have shown that there are four flow regimes: laminar, nonlinear laminar, turbulent transition, and fully developed turbulent flow (Ward, 1964 [9]; Wright (1968) [27]; Kovacks (1969) [15], Dybbs and Edwards (1975) [28]; R.M. Fand et al. (1987) [29]; H. Huang and S. Ayoub (2007) [30]) and, consequently, there is a debate regarding whether the dimensionless coefficients $A_1$ and $A_2$ considered in Equations (10) and (12) for a determined porous medium, are constants throughout a wide range of seepage velocity spanning the four flow regimes, in other words maintaining the relationship determined by Equation (10) without producing sudden jumps as occurred in the flow in tubes.

For that purpose, it is interesting to linearise Equation (10) multiplying both parts of the equation by $R_e$ to obtain Equation (16):

$$\lambda =A_1+A_2·R_e$$

(16)

where $\lambda$, the linearised generalised friction faction, is determined by Equation (17):

$$\lambda =f·R_e$$

(17)

Equation (16) is referred to as a linearised generalised equation $\left[R_e,\lambda \right]$.

Figure 2 shows a $\left[R_e,\lambda \right]$ diagram similar to that used by R.M Fand et al. (1987) [29]. In accordance with Equation (16), each porous material is represented by a line whose curve corresponds with the quadratic generalised dimensionless coefficient $A_2$; the cut on the ordinate axis corresponds to the linear generalised dimensionless coefficient $A_1$. The laminar regime is represented by a horizontal line. This type of graph is used to check if, for the same porous material, each non-Darcy flow regime remains determined by their corresponding value pairs $\left[A_1,A_2\right]$ or, on the contrary, such values are constant for the three non-Darcy flow regimes.

In this debate, the following studies stand out: McCorquodale et al. (1978) [31] and Fand et al. (1987) [29]. The first group used granular materials, of various sizes, shapes, and angularity, whereas the second worked with spheres with a range of sizes from 2.00 mm < D < 4.00 mm. The studies both proposed different coefficient values, $A_1$ and $A_2$, for the transition zones that they detected in their tests: nonlinear laminar and turbulent transition.

2.2. Resistance Formulas

Next, we describe the main nonlinear resistance formulas based on the adopted characteristic length $L_c$.

2.2.1. Resistance Formulas Based on the Representative Diameter of the Particles

In 1952, Sabri Ergun [12] carried out a study to develop a formula for general application based on the dimensionless groups:

$$f_E=D_p·g·\frac{n^3}{\left(1-n\right)}·\frac{i}{V^2}$$

(18)

$$R_E=\frac{D_p·V}{v}$$

(19)

where $f_E$ is Ergun’s friction factor and $R_E$ is Ergun’s Reynolds number. The parameter $D_p$ was adopted as the representative diameter of the particles, that is, the average effective diameter of the granular material and determined by Equation (20):

$$D_p=\frac{6}{S_e}$$

(20)

where $D_p$ corresponds with the diameter of a sphere that has the same specific surface area, $S_e$, as the particle.

The continuous curve, as determined by Equation (10) and defined by the author, we refer to it as the Ergun equation $\left[R_E,f_E\right]$ as follows:

$$f_E=72·\alpha ·\frac{\left(1-n\right)}{R_E}+\frac{6}{8}·\beta$$

(21)

where $\alpha$ is the linear dimensionless coefficient and $\beta$ is the quadratic dimensionless coefficient (by Sabri Ergun and A. A. Orning).

Finally, substituting the values of $f_E$ and $R_E$ of Equations (18) and (19) into Equation (21) we get the quadratic equation from Ergun (1952) [12] as Equation (22):

$$i=\frac{v}{g}·72·\alpha ·\frac{1}{D_p^2}·\frac{{\left(1-n\right)}^2}{n^3}·V+\frac{1}{g}·\frac{6}{8}·\beta ·\frac{\left(1-n\right)}{n^3}·\frac{1}{D_p}·V^2$$

(22)

This equation proposed by Ergun (1952) [12] had, in fact, been previously developed by Sabri Ergun and A. A. Orning (1949) [8] who adopted as characteristic length $L_c$, the specific surface area of the particles $S_e$. The equation is obtained by substituting in Equation (22) the value of the effective diameter $D_p$ defined by Equation (20). These authors developed this equation based on Kozeny’s (1927) hypothesis [32], which is based on a capillary model that, as the authors show “the granular bed is equivalent to a group of parallel and equal-sized channels, such that the total internal surface and the free internal volume, are equal to the total packing surface area and the void volume, respectively, of the randomly Packed bed.” However, they added that “For a packed bed, the flow path is sinuous, and the stream lines frequently converge and diverge. The kinetic losses, which occur only once for the capillary, occur with a frequency that is statistically related to the number of particles per unit length. For these reasons, a correction factor must be applied to each term. These factors may be designated as α and β.” For these reasons, in agreement with the authors, it is essential to consider some correction factors in each linear and quadratic expression of Equation (22) to consider these characteristics of the porous medium. The authors designated them as $\alpha$ and $\beta$, respectively.

Equation (22) considers the existence of a function of porosity for each expression: $f_L$ for the linear expression ($r$) and $f_T$ for the quadratic expression ($s$) of Equation (2) where:

$$f_L=\frac{{\left(1-n\right)}^2}{n^3}$$

(23)

$$f_T=\frac{\left(1-n\right)}{n^3}$$

(24)

As per Ergun, Leva M., and Grimmer M. (1947) [33] they confirmed these functions of porosity. Ergun himself carried out various tests on a wide range of porosity variation $n$ to check the validity of these functions.

With respect to the shape, angularity, and surface roughness of the particles, these physical parameters seem implicit in the effective diameter, $D_p$, which relate to the specific surface area S_e through Equation (20).

Finally, with respect to the packed porous material he indicates that “the orientation of the randomly packed beds is not susceptible to exact mathematical formulation”, and, consequently, “the effect of the orientation was not included.”

Ergun, worked on 640 experiments including his own, made up of crushed porous solids, and those obtained through other authors such as Burke and Plummer (1928) [6], who worked with lead shot (spherical shape) of reduced sizes 1.48, 3.08, and 6.34 mm and Morcon (1949) [11], who worked with capsules, cylinders, nodules, and spheres. Ergun represented the three series of data in the diagram $\left[R_E,f_E\right]$ checking that they correctly matched Equation (21) The fluids used in this case were the gases of CO₂, N₂, CH₄, and H₂. As a result of this, he obtained the universal values of $\alpha$ = 2.08 and $\beta$ = 2.33 applicable to all porous media.

According to Ergun, most dispersions happen with porous materials that include a mixture of sizes (non-uniform materials), and with those in which the relationship between the diameter of the permeameter ($D_x$) and the representative size of the particle $D$ is less than 10 (influence of the wall effect). In such cases, the corresponding tests were not considered.

Although the author did not analyse the theoretical significance of the parameters $\alpha$ and $\beta$, he did confirm that over a wide range of porosities he found no relationships of $\alpha$ and $\beta$ with the porosity.

Later, Frank Engelund (1953) [34] carried out a study to analyse the influence of turbulence in subterranean waters on uniform limestone sands, working on his own tests and those of other authors (Lindquist, 1933 [35] and Chardabellas, 1940 [36,37]). For his tests, he worked only on three samples, two of them were with $D$ = 2.6 mm and one with $D$ = 1.4 mm. The proposed quadratic equation Equation (25) was:

$$i=\frac{v}{g}·{\alpha }_0·\frac{1}{D_e^2}·\frac{{\left(1-n\right)}^3}{n^2}·V+\frac{1}{g}·{\beta }_0\frac{\left(1-n\right)}{n^3}·\frac{1}{D_e}·V^2$$

(25)

He used as the representative diameter of particle $D$ the parameter $D_e$, which is the equivalent diameter that corresponds to the diameter of a sphere with the same volume of the particle; α₀ is the linear dimensionless coefficient of Engelund and ${\beta }_0$ is the quadratic dimensionless coefficient of Engelund. According to the author, “${\alpha }_0$ and ${\beta }_0$ are dimensionless numerical constants depending, for uniform soil, on the structure and the grain shape.”

The structure of Equation (25) is remarkably similar to Equation (22) given by Ergun (1952) [12]. The difference in the representative diameter $D$ of the adopted particle, $D_p$ by Ergun (1952) [12] and $D_e$ by Engelund (1953) [34], and the function of linear porosity $f_L^{’}$ obtained as:

$$f_L^{’}=\frac{{\left(1-n\right)}^3}{n^2}$$

(26)

According to Frank Engelund (1953) [34], this porosity function $f_L^{’}$ of the linear expression fits better than the function $f_L$ from Ergun Equation (23) to the measurements obtained by Rose (1953) [37] and Franzini (1951) [38], especially the latter, who worked with a wide interval of porosities 0.270 < n < 0.476.

Ergun (1952) [12] proposed universal values for the dimensionless coefficients $\alpha$ and $\beta$. However, Engelund (1953) [34] did consider the shape and angularity of the particles through the dimensionless coefficients ${\alpha }_0$ and ${\beta }_0$ for which he proposed the following values:

${\alpha }_0$ = 780 for uniform spherical particles (using data from Lindquist (1933)); ${\alpha }_0$ = 1000 for uniform rounded sands (using data from Chardabellas (1940)), and ${\alpha }_0$ = 1500 or above for angular sands, (using own experimental data).

${\beta }_0$ = 1.8 for uniform spherical particles (using data from Lindquist (1933)); ${\beta }_0$ = 2.8 for uniform rounded sands (using data from Chardabellas (1940)), and ${\beta }_0$ = 3.6 or above for angular and uniform sands, (using own experimental data).

Finally, we must point out that it is simpler to measure as representative diameter of the particle $D$ the value of $D_e$ than the value of $D_p$, which requires the measurement of the particle specific surface area $S_e,$ in accordance with Equation (20).

In 1979, Stephenson carried out a study whose main purpose was to research the fully developed turbulent flow regime. On the basis of an analogy with the flow in pipes (Darcy–Weisbach equation), he considered that the gradient $i$ was proportional to the expression Equation (27):

$$i\propto \frac{1}{R_h}·\frac{V_p^2}{2·g}$$

(27)

As the hydraulic mean radius $R_h$ is proportional to the size of the particle $D$ (Leps (1973)) [39] and, furthermore, this variable is more easily measured, he suggested using the exponential equation, i.e., Equation (1) with coefficient b equal to 2 for Equation (28):

$$i=\frac{1}{g}·f_d·\frac{1}{D}·\frac{1}{n^2}·V^2$$

(28)

where $f_d$ is the particle friction factor and, as the author noted, “is actually the function of the Reynolds number $R_d$”. The author used as the representative diameter of the particle $D$, the parameter D₅₀ for lack of other data in the bibliography. However, for conceptual reasons, we will continue working with the representative diameter of the particles $D$.

He considered the dimensionless groups [$R_d,f_d$] where:

$$f_d=D·g·\frac{i}{V_p^2}$$

(29)

and $R_d$ is defined by the Equation (15).

He represented the data of three types of porous materials: smooth spheres, river gravel, and crushed aggregates on the diagram from Stephenson (1979) [17] [$R_d,f_d$]. This data came from tests carried out by the same author and by various researchers: Dudgeon (1966) [10], Volker (1969) [40], Leps (1973) [39], and Cedergren (1977) [41]. With this experimental data and their corresponding ranges of particle sizes and hydraulic gradients, he matched three smooth curves whose lines coincided in the laminar regime. From this line, they gradually separate in the transition regime towards turbulence with three asymptotic values: $f$_d = 1 for smooth spheres, $f$_d = 2 for river gravel, and $f$_d = 4 for crushed aggregate (see Figure 3).

For laminar flow, considering the analogy with the flow in pipes and based on data from Dudgeon (1966) [10] and Cedergren (1977) [41], he proposed the relationship:

$$f_d=\frac{K_l}{R_d}$$

(30)

where $K_l$ is the linear dimensionless coefficient from Stephenson (1979) [17]. The value of $K_l$ obtained from the experimental data was 800. The limit of the Reynolds number $R_d$ for the laminar regime proposed by Stephenson was $R_d$ = 10⁻⁴.

In the transition zone, he proposed Equation (31) of the type given by Equation (10):

$$f_d=\frac{K_l}{R_d}+K_t$$

(31)

where $K_t$ is the quadratic dimensionless coefficient from Stephenson (1979) [17]. According to the author, Reynolds numbers of $R_d$ > 10,000 produce the fully developed turbulent flow regime and, in this case, the author indicated that it was only a function of the shape and angularity of the particles.

Ultimately, Equation (31) (similar to the generalised equation Equation (10)), the Ergun equation, i.e., Equation (21), and by implication Equation (25) by Frank Engelund, represent continuous curves for each porous medium in a generalised diagram [R_e,$f$] that tend to an asymptotic value for fully developed turbulent regime (see Figure 1 and Figure 3).

Finally, substituting Equations (15), (29) and (31), we get the quadratic equation Equation (32):

$$i=\frac{v}{g}·K_l·\frac{1}{D^2}·\frac{1}{n}·V+\frac{1}{g}·K_t·\frac{1}{D}·\frac{1}{n^2}·V^2$$

(32)

Equation (32), according to Li B. et al. (1998), was not initially proposed by Stephenson (1979) [17] and these authors termed it as the modified Stephenson equation.

2.2.2. Resistance Formulas Based on Intrinsic Permeability

Ward (1964) [9] obtained a quadratic formula for general application taking $\sqrt{K_0}$ as the characteristic length. For this, he worked with the dimensionless groups expressed in Equations (33) and (34):

$$f_k=\sqrt{K_0}·g·\frac{i}{V^2}$$

(33)

$$R_k=\frac{\sqrt{K_0}·V}{v}$$

(34)

where $f_k$ is Ward’s friction factor and $R_k$ is Ward’s Reynolds number.

For the transition zone, in a $\left[R_k,f_k\right]$ diagram of the type produced by Equation (10), Ward (1964) [9] proposed Equation (35):

$$f_k=\frac{1}{R_k}+C$$

(35)

where $C$ is the quadratic dimensionless coefficient of Ward which, as indicated by Ahmed and Sunada (1969) [14], “Previous investigators have interpreted the term $C$ as a constant reflecting geometric properties of the medium.”

Considering the three previous equations, the author obtained the quadratic equation, i.e., Equation (36):

$$i=\frac{v}{g}·\frac{1}{K_0}·V+\frac{1}{g}·C·\frac{1}{\sqrt{K_0}}·V^2$$

(36)

Ward (1964) [9] had the same aim as Ergun (1952) [12], to adjust a single value of $C$ so that Equation (35) might be generally applicable. Engelund (1953) [34], and Stephenson (1979) [17], in fact, proposed three equations also for general application, but depending on the structure of the porous material and the shape of the particles.

Ward (1964) [9] used 20 different porous media: glass bead, ion exchange resin, sands, gravel, granular activated carbon, and anthracite. The size of each porous medium was defined through the geometric mean of the particle sizes $M_g$ in each porous medium. The size interval was 0.27 mm < $M_g$ < 16.10 mm. In total, 53 tests were carried out and the fluid used was water. The Reynolds number interval was 0.122 < $R_k$ < 18.10, with the size interval of $M_g$ studied, he obtained an adjusted value for $C$ of 0.55 with a standard deviation ${\sigma }_s$ of 0.024.

The author noted that “D.K. Todd (1959) [42] showed a plot similar to Figure 1, except that an average grain diameter is used in place of the square root of the permeability in Equations (10) and (16). Because the average grain diameter is not sufficient to characterize a porous medium, there is considerable scatter in the plotted points.”

Regardless of the existence of the smooth transition curve (see Figure 1), Ward (1964) [9] proposed, for engineering applications, the division into four types of flow regimes: laminar, nonlinear laminar, turbulent transition, and fully developed turbulent; with the following limits of the Reynolds number $R_k$ among them: laminar $R_k$ < 0.0182, nonlinear laminar $R_k$ < 1.82, and turbulent transition $R_k$ < 182.

To conclude, we must point out that Ward attributes, in some cases, the deviations seen in respect of Equation (35) to an inadequate determination of the intrinsic permeability $K_0$, which must be obtained with very low gradients to be in a laminar regime in accordance with Equation (6).

Continuing with Ward’s study (1964) [9], Ahmed and Sunada (1969) [14] developed the same quadratic equation, i.e., Equation (36) from the Navier–Stokes’ equations, based on the hypothesis that the porous medium is homogeneous and isotropic on a macroscopic scale, and that the chemical and thermodynamic effects are small. His study focused specifically on the nonlinear laminar regime.

Ahmed and Sunada (1969) [14] stated that “for flow through porous media, convective accelerations are always present, whereas turbulence is a random phenomenon dependent upon of flow velocity and space geometry.” According to the authors, “This fact was demonstrated experimentally by Schneebeli (1955) [43] who used dye to identify the flow. He injected dye into the flow at various velocities (steady-state conditions) and found that, even though measurements of gradients and velocities indicated nonlinear flow, the dye assumed laminar characteristics, i.e., stream-lined flow. Increasing the flow velocities approximately four times caused the dye from one channel to mix with the dye of another, indicating that departure from Darcy’s Law should be the result of convective acceleration of the fluid within the pores space.”

More recently, H. Huang and J. Ayoub (2007) [30], concluded, “Derivation of the Forchheimer equation from the Navier–Stokes’ equation reveals that the nature of the Forchheimer flow regime is laminar with inertial effect. The inertia resistance factor β can be used to characterize this flow regime and is therefore an intrinsic property of the porous media.”

Along the same lines, Balhoff, M.T. et al. (2009) [44] indicated, “The constant, β, is referred to us as the non-Darcy coefficient and, like permeability, is an empirical value specific to the porous medium. It represents the additional inertial resistance caused by the converging/diverging and tortuous medium geometry.”

The phenomenon of the turbulence in porous media is very complex. Recently, Sidiripoulou, M.G. et al. (2007) [19] referred to studies by Skjetne and Auriault (1998) [45], Panfilov et al. (2003) [46], and Fourar et al. (2004) [47] on the mechanisms of turbulence, which were related to the separation of the layer limit and the recircularisation of the vortices formed.

2.2.3. Resistance Formulas Based on the Hydraulic Mean Radius

In 1998, Li B et al. (1998) [18], based on an analogy with the flow in pipes and using as characteristic length L_c and the hydraulic mean radius $R_h$ defined by Equation (7), proposed the dimensionless groups Equations (37) and (38):

$$f_p=4·R_h·2g·\frac{i}{V_p^2}$$

(37)

$$R_p^{’}=\frac{R_h·V_p}{v}$$

(38)

where $f_p$ is the pore friction factor and $R_p^{’}$ is the Reynolds number based on $R_h$.

The transition curve proposed in the pore diagram $\left[R_p^{’},f_p\right]$ was determined by Equation (39):

$$f_p=\frac{\alpha ’}{R_p}+\beta ’$$

(39)

where $\alpha$’ and $\beta$’ were dimensionless coefficients of the pores from the linear and quadratic expressions, respectively.

Substituting Equations (37)–(39), the authors obtained the quadratic equation:

$$i=\frac{v}{8·g}·{\alpha }^{\prime }·\frac{1}{R_h^2}·V_p+\frac{1}{8·g}·{\beta }^{\prime }·\frac{1}{R_h}·V_p^2$$

(40)

With the experimental data provided by the University of Ottawa (Hansen (1992)) [48], which included materials for rockfill dams in intervals of 16.0 mm < $D$ < 40.0 mm, they obtained values of 98 for $\alpha$’ and 3 for $\beta$’. The permeameter used had a diameter $D_x$ of 300 mm. The pore Reynolds number $R_p^{’}$, from which the total turbulent regime was developed, was 200.

They subsequently extended this data, with that supplied by Stephenson (1979) [17] and Li and Hu (1988) [18] obtaining a value of 1279 for the nonlinear dimensionless coefficient $K_l$ and 3.84 for the quadratic dimensionless coefficient $K_t$ which is shown in the modified Stephenson equation, i.e., Equation (32). These values are in the same order of magnitude as those obtained by Stepheson for crushed aggregate, i.e., $K_l=800$ and $K_t=4.00$.

More recently, Mohammad-Bagher Salahi et al. (2015) [22], working on rounded granular materials (2.10 mm < $D$ < 17.78 mm) and aggregate (1.77 mm < $D$ < 16.62 mm), obtained values $K_l=488 \mathrm{and} K_t=4725$ for rounded aggregate $K_l=588 \mathrm{and} K_t=5550$ for crushed aggregate. Both values of $k_t$ are higher than those proposed by Stephenson ($k_t$ = 2.00 for rounded aggregate and $K_t$ = 4.00 for crushed aggregate). The Reynolds number interval, $R_d$, with which they developed the tests, was 10 < $R_d$ < 1882.

Additionally, Li B et al. (1998) [18] represented the experimental data in the diagram from Stephenson [$R_d,f_d$]. According to an analysis of this diagram, they proposed that the fully developed turbulent regime (asymptotic curve) should be obtained for values of $R_d$ > 2000. This value is consistent with the data previously obtained for the pore Reynolds number $R_p^{’}$ > 200 if we consider that the hydraulic mean radius $R_h$ is approximately 10% of the representative size of the particle $D$ (Parkin, 1991) [25]. However, this limit value of $R_d$ is lower to that proposed by Stephenson of $R_d$ = 10,000.

2.2.4. On the Physical Parameters $r$ and $s$ of the Forchheimer Equation

Sidiropoulou et al. (2007) [19] obtained empirical relationships for a general application to any porous medium for the coefficients r and s in the function of physical parameters such as representative particle size $D$ and porosity n on the basis of an analysis of multiple regression using the data from various authors: Ward (1964) [9], Ahmed and Sunada (1969) [14], Arbhabhirama and Dinoy (1973) [16], Ranganadha Rao and Suresh (1979) [49], Tyagi and Todd (1970) [50] who used the data from Dudgeon (1966) [10], Venkataraman and Rao (1988) [51], and Bordier and Zimer (2000) [52]. The number from the available experimental data (N) was 115. However, in many cases, the complete data needed for the adjustment was not available, that is, $r,s,D$ and $n$.

They obtained three different empirical relationships for the expressions $r$ and s to study how the porosity influenced parameters r and s:

$$r = 0.00333·D^{-1.500403}·n^{0.060350} (R^2 = 0.9108; N = 55)$$

(41a)

$$s = 0.194325·D^{-1.265775}·n^{-1.141417} (R^2 = 0.8715; N = 49)$$

(41b)

$$r = 0.0002789·D^{-1.502361}·{(1-n)}^{-0.216014} (R^2 = 0.9142; N = 55)$$

(42a)

$$s = 1.228873·D^{-1.263314}·{(1-n)}^{-1.532475} (R^2 = 0.8762; N = 49)$$

(42b)

$$r = 6.527953·{10}^{-15}·D^{-1.54745}·n^{-16.068711}·{(1-n)}^{-23.157232} (R^2 = 0.9188; N = 55)$$

(43a)

$$s = 1.107768·{10}^{-10}·D^{-1.30182}·n^{-13.836369}·{(1-n)}^{-18.365290} (R^2 = 0.8806; N = 49)$$

(43b)

where $D$ is given in metres, $r$ in seconds per metre, and $s$ in seconds squared per metre squared.

The best adjustment obtained was for the Equation (43a,b).

The issue with the previous empirical relationships is that they were based on porous materials made up of particles with different geometries: glass beads, granular activated carbon, ion exchange resin, sand, gravel, anthracite coal, angular gravel, round river gravel, blue metal, river gravel, marbles, and glass spheres.

In fact, these equations, were an attempt to provide a general application equation as per the studies made by Ergun (1952) [12] and Ward (1964) [9]. Frank Engelund (1953) [34] went further by proposing three general equations in the function of the shape and angularity of particles through the dimensionless coefficients ${\alpha }_0$ and ${\beta }_0$. Stephenson (1979) [17] also proposed three different equations: smooth spheres, river gravel, and crushed aggregate (see Figure 3).

3. Analysis of the Relationships among Parameters of the Different Formulas of Resistance

To adequately define the filtration phenomenon through a porous medium, we use three types of equations: the generalised equation $\left[R_e,f\right]$ determined by Equation (10), the quadratic generalised equation, i.e., Equation (12), and the linear generalised equation $\left[R_e,\lambda \right]$ determined by Equation (16), all of them using the general characteristic length $L_c$.

The purpose of this section is to standardise the formulas described in the previous section. For this, we have chosen the quadratic generalised equation as a base in accordance with Equation (12), it includes the generalised characteristic length $L_c$ and the generalised dimensionless coefficient $A_1$ (linear term) and $A_2$ (quadratic term).

To arrive at this standardisation, the following mathematical relationships were studied:

(a)

Among the characteristic lengths, R_h, $\sqrt{K_0}$, and D;

(b)

Among the Reynolds numbers $R_p,R_k \mathrm{and} R_d$;

(c)

Among the different laminar dimensionless coefficients $\alpha$, ${\alpha }_0$, $K_l$, and $\alpha$′, and quadratic dimensionless coefficients $\beta$, ${\beta }_0$, $K_t,\beta$′, and $C$.

These relationships should allow us to compare the results proposed by different authors who use different characteristic lengths $L_c$ and, as a result, different values for the laminar dimensionless coefficients, (α, $\alpha$′, ${\alpha }_0$, and $K_l$), quadratic dimensionless coefficients ($\beta ,$ ${\beta }_0$, $\beta$′, $C$, and $K_t$), and also different values for the limit of the Reynolds number $R_e$ that define the zones of the flow regime: laminar, nonlinear laminar, turbulent transition, and fully developed turbulent.

The relationships that we have obtained in this section are based on uniform granular materials where the uniformity coefficient is not included $C_u$.

3.1. Equations with Characteristic Length Based on the Hydraulic Diameter

In accordance with the analogy of the flow in pipes, for noncircular sections, we can take the characteristic length to be $L_c$ and the hydraulic mean diameter as $D_h$:

$$D_h=4·R_h$$

(44)

If we substitute Equation (44) in Equation (12) we get:

$$i=\frac{v}{2·g}·A_1^{’}·\frac{1}{D_h^2}·V_p+\frac{1}{2·g}·A_2^{’}·\frac{1}{D_h}·V_p^2$$

(45)

where $A_1^{’}$ is the linear dimensionless coefficient corresponding to $L_c=D_h$ and $A_2^{’}$ is the quadratic dimensionless coefficient corresponding to $L_c=D_h$.

Substituting Equation (44) in Equation (45) we obtain Equation (46):

$$i=\frac{v}{2·g}·\frac{A_1^{’}}{16}·\frac{1}{R_h^2}·V_p+\frac{1}{2·g}·\frac{A_2^{’}}{4}·\frac{1}{R_h}·V_p^2$$

(46)

If we compare Equation (46) with Equation (40) proposed by Li, B. et al. (1998) [18], we get the dimensionless coefficients $A_1^{’}\left(\mathrm{linear}\right) \mathrm{and} A_2^{’}$(quadratic):

$$A_1^{’}=4·\alpha ’$$

(47)

$$A_2^{’}={\beta }^{\prime }$$

(48)

Considering Equation (7), which defines the hydraulic mean radius $R_h$ we obtain the expression for the specific surface area $S_e$:

$$S_e=\frac{n}{\left(1-n\right)}·\frac{1}{R_h}$$

(49)

Considering Equations (5) and (49), which defines the pore velocity $V_p$ and substituting Equation (46) we get:

$$i=\frac{v}{32·g}·A_1^{’}·\frac{v}{g}·\frac{{\left(1-n\right)}^2}{n^3}·S_e^2·V+\frac{1}{8·g}·A_2^{’}·\frac{1}{g}·\frac{\left(1-n\right)}{n^3}·S_e·V^2$$

(50)

Equation (50) is the same as that developed by Sabri Ergun and A. A. Orning (1949) [8] where the dimensionless coefficients $A_1^{’}$ and $A_2^{’}$ have the values:

$$A_1^{’}=64\alpha$$

(51)

$$A_2^{’}=\beta$$

(52)

Considering the Equations (47), (48), (51), and (52):

$$\alpha =\frac{\alpha ’}{16}$$

(53)

$$\beta =\beta ’$$

(54)

In addition, the specific surface area of the packing $S$ is related to the average specific surface area of the particles $S_e$ through Equation (9) with which we get:

$$S_e^2=\frac{S^2}{{\left(1-n\right)}^2}$$

(55)

Substituting this last value in Equation (50), we obtain the quadratic equation with characteristic length $L_c$ the surface of the packing $S$:

$$i=2\alpha ·\frac{v}{g}·\frac{1}{n^3}·S^2·V+\frac{1}{8}·\beta ·\frac{1}{g}·\frac{1}{n^3}·S·V^2$$

(56)

Alternatively, S.P. Burke and W.B. Plumer (1928) [6] developed the exponential equation with exponent (2−$b$):

$$i=K_b·\frac{1}{g}·\frac{V^2}{n^3}·{\left(\frac{\mu ·S}{\rho ·V}·\frac{n}{\left(1-n\right)}\right)}^{2-b}$$

(57)

where $K_b$ is a dimensionless coefficient, which according to the authors includes the shape of the porous material and the symmetry of the packing.

The authors did not consider taking a quadratic equation from this. However, this can be achieved by doing no more than taking $b$ = 1 for the linear component and $b$ = 2 for the quadratic component. Through this approach we get the equation known as the quadratic equation from Burke and Plumer (1928) [6]:

$$i={\left[K_b\right]}_L·\frac{v}{g}·\frac{1}{n^3}·\frac{n}{\left(1-n\right)}·S^2·V+{\left[K_b\right]}_T·\frac{1}{g}·\frac{1}{n^3}·S·V^2$$

(58)

where ${\left[K_b\right]}_L$ is the linear dimensionless coefficient by Burke and Plumer ${\left[K_b\right]}_T$ is the quadratic dimensionless coefficient by Burke and Plumer.

Now, seeing Equation (58) and comparing it with Equation (56) we see that the dimensionless coefficients ${\left[K_b\right]}_L$ and ${\left[K_b\right]}_T$ are determined by the expression:

$${\left[K_b\right]}_L=2\alpha ·\frac{n}{\left(1-n\right)}$$

(59)

$${\left[K_b\right]}_T=\frac{\beta }{8}$$

(60)

Accordingly, the previous equations show that a mathematical relationship exists among the formulas of the authors cited.

3.2. Equations with Characteristic Length Based on the Intrinsic Permeability

Now, if we consider, as characteristic length in Equation (12), the quadratic root of the intrinsic permeability $\sqrt{K_0}$ adopted by Ward (1964) [9], Ahmed and Sunada (1969) [14], and Arbhabhirama and Dinoy (1973) [16], among others, we obtain:

$$i=\frac{v}{2·g}·A_1^{’’}·\frac{1}{K_0}·V_p+\frac{1}{2·g}·A_2^{’’}·\frac{1}{\sqrt{K_0}}·V_p^2$$

(61)

where $A_1^{’’}$ is the laminar dimensionless coefficient corresponding to $L_c=\sqrt{K_0}$ and $A_2^{’’}$ is the quadratic dimensionless coefficient corresponding to $L_c=\sqrt{K_0}.$

If we compare Equation (61) with Equation (30) by Ward (1964), we get the dimensionless coefficients $A_1^{’’}$ and $A_2^{’’}$:

$$A_1^{’’}=2·n$$

(62)

$$A_2^{’’}=2·C·n^2$$

(63)

3.3. Equations with Characteristic Length Based on the Representative Size of Particles

We are going to consider in general the representative size of the particle $D$, as obviously, the formulas developed in this case can be applied for any representative size such as $D_e,D_a,D_p,D_{50},M_g$, etc.

For this general case of parameter $D$ and applying Equation (12) we obtain the equation:

$$i=\frac{v}{2·g}·A_1^{’’’}·\frac{1}{D^2}·V_p+\frac{1}{2·g}·A_2^{’’’}·\frac{1}{D}·V_p^2$$

(64)

where $A_1^{’’’}$ is the linear dimensionless coefficient corresponding to $L_c=D$ and $A_2^{’’’}$ is the quadratic dimensionless coefficient corresponding to $L_c=D$.

If we observe Equation (64) and compare it with modified Stephenson equation, i.e., Equation (36), we get the dimensionless coefficients $A_1^{’’’}$ and $A_2^{’’’}$:

$$A_1^{’’’}=2K_l$$

(65)

$$A_2^{’’’}=2K_t$$

(66)

Once the quadratic equations have been obtained for each characteristic length $L_c$ based on Equation (12), we can determine the mathematical relationships between the main parameters.

3.4. Relationships among Characteristic Lengths

Next, we determine the relationships between $R_h$ and $\sqrt{K_0}$ and between $R_h$ and $D$. With regard to the first relationship, being the linear components of Equations (45) and (61):

$$\frac{1}{2}·A_1^{’}·\frac{v}{g}·\frac{1}{D_h^2}=\frac{1}{2}·A_1^{’’}·\frac{v}{g}·\frac{1}{K_0}$$

(67)

Substituting the values of $A_1^{’}$ defined in Equation (51) and $A_1^{’’}$ defined in Equation (62) into Equation (67) we obtain:

$$32·\alpha ·\frac{1}{D_h^2}=n·\frac{1}{K_0}$$

(68)

$$D_h^2=\frac{32·\alpha ·K_0}{n}$$

(69)

Substituting the value of the hydraulic mean diameter $D_h$ defined in Equation (44) into Equation (69), we obtain the relationship between the hydraulic mean radius R_h and the intrinsic permeability $K_0$:

$$R_h=\frac{\sqrt{2·\alpha }}{n^{0.5}}·\sqrt{K_0}$$

(70)

We determine the relationship between the hydraulic mean radius $R_h$ and the representative size of the particle $D$, by equaling the linear components of the Equation (45) and the Equation (64):

$$\frac{1}{2}·A_1^{’}·\frac{v}{g}·\frac{1}{D_h^2}=\frac{1}{2}{A}_1^{’’’}·\frac{v}{g}·\frac{1}{D^2}$$

(71)

Substituting the values of $A_1^{’}$ defined in Equation (51) and $A_1^{’’’}$ defined in Equation (65) in the previous expression we obtain:

$$32·\alpha ·\frac{1}{D_h^2}=K_l·\frac{1}{D^2}$$

(72)

$$D_h^2=\frac{32·\alpha }{K_l}·D^2$$

(73)

With which we obtain the relationship:

$$R_h=\frac{\sqrt{2·\alpha }}{\sqrt{K_l}}·D$$

(74)

We determine the relationship between the hydraulic mean radius $R_h$ and the representative size of the particle $D$ by applying Equation (7). We know that, for a sphere, the specific surface area $S_e$ has a value of 6/$D$, where $D$ is the diameter of the sphere. To consider the shape, angularity, and surface roughness of the particles that may affect its specific surface area $S_e$ (Crawford, C.W. et al. (1986) [3]; Sabin G. C. W. and Hansen D. (1994) [4], we can define a coefficient $F$ that is determined by the expression:

$$S_e=\frac{6·F}{D}$$

(75)

In accordance with Equation (7) that defines hydraulic mean radius $R_h$ and considering Equation (72) we get the relationship:

$$R_h=\frac{n}{6·F·\left(1-n\right)}·D_{}$$

(76)

Loudon (1953) [53] used a coefficient similar to $F$ to determine the permeability of sands. The author used $D$ as the representative size of the particle, ($D_g$) defined as the average geometry between two consecutive sieves:

$$S_e=F·\frac{6}{D_g}$$

(77)

where:

$$D_g=\sqrt{D_n·D_{n-1}}$$

(78)

where $D_n$ and $D_{n-1}$ are the apertures in two consecutive sieves.

Loudon proposed the following values for the coefficient $F$:

(a)

Round sand, F = 1.10;

(b)

Semi angular sand, F = 1.25;

(c)

Angular sand, F = 1.40.

In accordance with Martins (1990) [54] and considering the coefficient of the shape $c$’ we can obtain the relationship:

$$c’ = 6·F$$

(79)

In accordance with the values of $c’$ proposed by Linford, A. and Saunders, D. (1967) [55] and Martins, R. and Escarameria, M. (1989) [56] (according to Martins), the equivalent values of the coefficient $F$ are:

(a)

Angular Particles $c$’ = 8.5 (F = 1.47).

(b)

Round particles $c$’ = 6.3 (F = 1.05).

These values are similar to those proposed by Loudon (1953) [53].

With these values and applying Equation (76), we can estimate the value of the hydraulic mean radius $R_h$.

Finally, we determine the relationships between $\sqrt{K_0}$ and $D$ considering the linear components of Equations (61) and (64):

$$\frac{1}{2}·A_1^{’’}·\frac{v}{g}·\frac{1}{K_0}=\frac{1}{2}·A_1^{’’’}·\frac{v}{g}·\frac{1}{D^2}$$

(80)

Substituting the values of $A_1^{’’}$ defined in Equation (62) and $A_1^{’’’}$ defined in Equation (65) in the previous expression we get:

$$n·\frac{v}{g}·\frac{1}{K_0}=K_l·\frac{v}{g}·\frac{1}{D^2}$$

(81)

Developing:

$$n·D^2=K_l·K_0$$

(82)

With which we finally get the relationship:

$$\sqrt{K_0}=\frac{n^{0.5}}{\sqrt{K_l}}·D$$

(83)

Summarizing above,

Table 1. Relationships among characteristic lengths L_c.

${\mathit{L}}_{\mathit{c}}$	${\mathit{A}}_{\mathbf{1}}$	${\mathit{A}}_{\mathbf{2}}$	${\mathit{R}}_{\mathit{e}}$	$\mathit{f}$	${\mathit{R}}_{\mathit{h}}$	$\sqrt{{\mathit{K}}_{\mathbf{0}}}$	$\mathit{D}$
Dh	64α	β	$\frac{4·R_h·V_p}{v}$	$8R_h·g·\frac{i}{V_p^2}$	4·Rh	$4·\frac{\sqrt{2\alpha }}{n}·\sqrt{K_0}$	$4·\frac{\sqrt{2\alpha }}{\sqrt{K_l}}·D$
$\sqrt{{\mathit{K}}_{\mathbf{0}}}$	2n	2·C·n2	$\frac{\sqrt{K_0}·V_p}{v}$	$\sqrt{K_0}·2g·\frac{i}{V_p^2}$	$\frac{n^{0.5}}{\sqrt{2\alpha }}·R_h$	1	$\frac{n^{0.5}}{\sqrt{K_l}}·D$
D	2Kl	2·Kt	$\frac{D·V_p}{v}$	$D·2g·\frac{i}{V_p^2}$	$\frac{\sqrt{K_l}}{\sqrt{2\alpha }}·R_h$	$\frac{\sqrt{K_l}}{n^{0.5}}·\sqrt{K_0}$	1

shows the relationships among the characteristic lengths $R_h,\sqrt{K_0}$ and D.

3.5. Relationships among Reynolds Numbers

First, we determine the relationships between the pore Reynolds number $R_p$ and the Ward Reynolds number $R_k$. Considering the pore Reynolds number $D_h$:

$$R_p=\frac{4·R_h·V_p}{v}$$

(84)

And substituting Equation (70) in Equation (84) we obtain:

$$R_p=4·\left(\frac{\sqrt{2·\alpha }}{n^{0.5}}·\sqrt{K_0}\right)·\frac{V}{n·v}$$

(85)

If we compare the last equation with Equation (28) that defines the Reynolds number from Ward ($R_k$) we obtain the relationship:

$$R_p=\frac{4·\sqrt{2·\alpha }}{n^{1,5}}·R_k$$

(86)

Now, we determine the relationships between the pore Reynolds number $R_p$ based on $D_h$ and the Reynolds number of the particles $R_d$.

Considering Equation (74) that relates to the hydraulic mean radius $R_h$ with the representative size of the particles $D$, and substituting in Equation (84) that defines pore Reynolds number $R_p$ based on $D_h$ we obtain:

$$R_p=4·\left(\frac{\sqrt{2·\alpha }}{\sqrt{K_l}}·D\right)·\frac{V_p}{v}=4·\frac{\sqrt{2·\alpha }}{\sqrt{K_l}}·\frac{D}{v}·\frac{V}{n}$$

(87)

Alternatively, in accordance with Equation (15) that defines the Reynolds number of the particles, and substituting this last expression into the previous expression, we obtain the relationship between $R_p$ and $R_d$:

$$R_p=\frac{4·\sqrt{2·\alpha }}{\sqrt{K_l}}·R_d$$

(88)

In addition, in accordance with Equation (76) that relates to the hydraulic mean radius $R_h$ and the representative size of the particles $D$, and substituting in Equation (84) that defines $R_p$ as:

$$R_p=4·\left[n/\left(6·F·\left(1-n\right)\right)·D\right]·\frac{V_p}{v}$$

(89)

In accordance with Equation (15) we finally obtain another relationship between $R_p$ and $R_d$:

$$R_p=\frac{4·n}{6·F·\left(1-n\right)}·R_d$$

(90)

Finally, we determine the relationships between the Reynolds number of Ward $R_k$ and the Reynolds number of the particles $R_d$.

In accordance with Equation (28) that defines the Reynolds number of Ward and considering Equation (83) that relates the intrinsic permeability $K_0$ with the representative size of the particles D, and substituting this last equation into the definition of the Reynolds number of Ward $R_k,$ we obtain:

$$R_k=\frac{D}{\sqrt{K_l}}·n^{0.5}·\frac{V}{v}=\frac{D}{\sqrt{K_l}}·n^{0.5}·\frac{n}{n}·\frac{V}{v}$$

(91)

Finally, in accordance with the definition of the Reynolds number of the particles $R_d$ defined by the Equation (15) we obtain the relationship:

$$R_k=\frac{n^{1.5}}{\sqrt{K_l}}·R_d$$

(92)

In accordance with the above,

Table 2. Relationships among the Reynolds numbers R_e.

${\mathit{L}}_{\mathit{c}}$	${\mathit{R}}_{\mathit{e}}$	${\mathit{R}}_{\mathit{p}}$	${\mathit{R}}_{\mathit{k}}$	${\mathit{R}}_{\mathit{d}}$
4Rh	Rp	1	$\frac{4\sqrt{2\alpha }}{n^{1.5}}·R_k$	$\frac{4\sqrt{2\alpha }}{\sqrt{K_l}}·R_d$
$\sqrt{{\mathit{K}}_{\mathbf{0}}}$	Rk	$\frac{n^{1,5}}{4\sqrt{2\alpha }}·R_p$	1	$\frac{n^{1,5}}{\sqrt{K_l}}·R_d$
D	Rd	$\frac{\sqrt{K_l}}{4\sqrt{2\alpha }}·R_p$	$\frac{\sqrt{K_l}}{n^{1.5}}·R_k$	1

shows the relationships among the Reynolds numbers $R_p,R_k \mathrm{and} R_d$.

3.6. Relationships among the Laminar Dimensionless Coefficients $\alpha ,{\alpha }_0,K_l and \alpha ’$

Considering Equations (74) and (76) and equaling both equations, we obtain:

$$\frac{\sqrt{2\alpha }}{\sqrt{K_l}}=\frac{n}{6·F·\left(1-n\right)}$$

(93)

Developing Equation (93), we obtain:

$$2\alpha ·36·F^2{\left(1-n\right)}^2=K_l·n^2$$

(94)

Solving the linear dimensionless coefficient α, we finally determine the relationship between the linear dimensionless coefficients $\alpha$ and $K_l$ Equation (95).

$$\alpha =\frac{n^2}{72·F^2·{\left(1-n\right)}^2}·K_l$$

(95)

Next, we relate the coefficients α and α₀. If we equal the linear components of Engelund Equation (25) considering the representative diameter of the particle $D$ instead of $D_e$ and Equation (64), we obtain:

$$\frac{v}{g}·{\alpha }_0·\frac{1}{D^2}·\frac{{\left(1-n\right)}^2}{n^2}=\frac{v}{2·g}·\left(2·K_l\right)·\frac{1}{D^2}·\frac{1}{n}$$

(96)

We get:

$$K_l={\alpha }_0·\frac{{\left(1-n\right)}^3}{n}$$

(97)

Substituting the values of $K_l$ in Equation (95) we obtain Equation (98):

$$\alpha =\frac{{\alpha }_0·\left(1-n\right)·n}{72·F^2}$$

(98)

Coefficient α′ relates to coefficient α through Equation (53)

Summarizing,

Table 3. Relationships among the laminar dimensionless coefficients.

${\mathit{L}}_{\mathit{c}}$	Laminar	$\mathit{\alpha }$	${\mathit{K}}_{\mathit{l}}$	${\mathit{\alpha }}_{\mathbf{0}}$
4Rh	α	1	$\frac{n^2}{72·F^2·{\left(1-n\right)}^2}·K_l$	$\frac{\left(1-n\right)·n}{72·F^2}·{\alpha }_0$
$\sqrt{{\mathit{K}}_{\mathbf{0}}}$	1
D	Kl	$\frac{72·F^2·{\left(1-n\right)}^2}{n^2}·\alpha$	1	$\frac{{\left(1-n\right)}^2}{n}·{\alpha }_0$
D	${\alpha }_0$	$\frac{72·F^2}{\left(1-n\right)·n}·\alpha$	$\frac{n}{{\left(1-n\right)}^3}·K_l$	1

shows the relationships among the laminar dimensionless coefficients $\alpha ,K_l \mathrm{and} {\alpha }_0$.

3.7. Relationships among Turbulent Dimensionless Coefficients $\beta ,{\beta }_0,K_t,{\beta }^{\prime } and C$

If we equal the quadratic components of Equations (45) and (61):

$$\frac{1}{2}·A_2^{’}·\frac{1}{g}·\frac{1}{4R_h}·\frac{1}{n^2}=\frac{1}{2}·A_2^{″}·\frac{1}{g}·\frac{1}{\sqrt{K_0}}·\frac{1}{n^2}$$

(99)

Considering the values of $A_2^{’}$ and $A_2^{″}$ defined by Equations (52) and (63) and substituting them in the first expression we get:

$$\beta ·\frac{1}{4R_h}=2·C·n^2·\frac{1}{\sqrt{K_0}}$$

(100)

In accordance with Equation (70) that defines the relationship between the hydraulic mean radius $R_h$ and the intrinsic permeability $K_0$ and substituting this equation in Equation (100):

$$\frac{\beta }{4}·\frac{n^{0.5}}{\sqrt{2\alpha }}·\frac{1}{\sqrt{K_0}}=2·C·n^2·\frac{1}{\sqrt{K_0}}$$

(101)

We finally get:

$$\beta =8·\sqrt{2·\alpha }·C·n^{1.5}$$

(102)

This last derived equation shows the dependence between the linear dimensionless coefficient $\alpha$ and the quadratic dimensionless coefficient $\beta$. As indicated by Huang H. and Ayoub J. (2007), “both expressions are intrinsically related, and the division is not arbitrary.”

If we equal the quadratic components of Equations (45) and (64) we obtain the relationship turbulent coefficients $\beta$ and $K_t$:

$$\frac{1}{2}·A_2^{’}·\frac{1}{g}·\frac{1}{4·R_h}·\frac{1}{n^2}=\frac{1}{2}·A_2^{’’’}·\frac{1}{g}·\frac{1}{D}·\frac{1}{n^2}$$

(103)

If we consider Equation (76) that relates the characteristic lengths $R_h$ and $D$ the values $A_2^{’}$ and $A_2^{’’’}$ given by Equations (52) and (66), substituting in the Equation (103) we obtain:

$$\beta ·\frac{6·F·\left(1-n\right)}{4·n}·\frac{1}{D}=2·K_t·\frac{1}{D}$$

(104)

Therefore, we finally obtain the relationship:

$$\beta = 86·nF·\left(1-n\right)·K_t$$

(105)

Finally, equaling Equations (104) and (105) that define the coefficient $\beta$:

$$8\sqrt{2·\alpha }·C·n^{1,5}=\frac{8}{6}·\frac{n}{F·\left(1-n\right)}·K_t$$

(106)

By clearing C in Equation (106) we obtain the relationship between the turbulent dimensionless coefficients $C$ and $K_t$:

$$C=\frac{1}{6}·\frac{1}{\sqrt{2\alpha }}·\frac{1}{F·n^{0.5}\left(1-n\right)}·K_t$$

(107)

Finally, if we equal the quadratic components of Equation (25) of Engelund (1953) and Equation (64):

$$\frac{1}{g}·{\beta }_0·\frac{\left(1-n\right)}{n^3}·\frac{1}{D}=\frac{1}{2·g}·\left(2·K_t\right)·\frac{1}{D}·\frac{1}{n^2}$$

(108)

We obtain the relationship:

$$K_t={\beta }_0·\frac{\left(1-n\right)}{n}$$

(109)

And substituting in Equation (108) we obtain:

$$\beta =\frac{8}{6}·\frac{1}{F}·{\beta }_0$$

(110)

Coefficient ${\beta }^{\prime }$ is equal to coefficient $\beta$ in accordance with Equation (54).

Accordingly,

Table 4. Relationships among the turbulent coefficients.

${\mathit{L}}_{\mathit{c}}$	Turbulent	$\mathit{\beta }$	$\mathit{C}$	${\mathit{K}}_{\mathit{t}}$	${\mathit{\beta }}_{\mathbf{0}}$
4Rh	β	1	$8·\sqrt{2\alpha }·n^{1.5}·C$	$\frac{8}{6}·\frac{n}{F\left(1-n\right)}·K_t$	$\frac{8}{6}·\frac{1}{F}·{\beta }_0$
$\sqrt{{\mathit{K}}_{\mathbf{0}}}$	C	$\frac{1}{8·\sqrt{2\alpha }·n^{1.5}}·\beta$	1	$\frac{1}{6·\sqrt{2\alpha }·F{n}^{0.5}\left(1-n\right)}·\beta$	$\frac{1}{6·\sqrt{2\alpha }·F·n^{1.5}}·{\beta }_0$
D	Kt	$\frac{6·F\left(1-n\right)}{8n}·\beta$	$6·\sqrt{2\alpha }·F·n^{0.5}·\left(1-n\right)·C$	1	$\frac{\left(1-n\right)}{n}·{\beta }_0$
D	${\beta }_0$	$\frac{6}{8}·F·\beta$	$6·\sqrt{2\alpha }·F·n^{1.5}·C$	$\frac{n}{\left(1-n\right)}·K_t$	1

shows the relationships among the dimensionless coefficients $\beta ,C,K_t \mathrm{and} {\beta }_0$.

4. Conclusions

The seepage process in a coarse porous medium can be represented by three types of equation: a generalized equation $\left[R_e,f\right]$, i.e., Equation (10), a generalized quadratic equation, i.e., Equation (12), and a generalized linear equation $\left[R_e,\lambda \right]$, i.e., Equation (16); all these equations consider a characteristic length $L_c$.

The linearized equation, i.e., Equation (16), allows one to check whether the coefficients $A_1$ and $A_2$ remain constant or, on the contrary, vary throughout the three non-Darcy flow regimes: nonlinear laminar, turbulent transition, and fully developed turbulent.

On the basis of the analysis of the different relationships between the gradient and the seepage velocity developed by different authors, all of them based on Equation (1) proposed by Forchheimer (1901) [1], it can be concluded that all the physical parameters considered in the different formulations are related to each other. Such parameters are the characteristic length ($L_c$, Table 1), the Reynolds number ($R_e$, Table 2), the dimensionless coefficient of the linear term r ($A_1$, Table 3), and the dimensionless coefficient of the quadratic term s ($A_2$, Table 4).

In this paper, we establish the equations that relate these parameters, thus, facilitating comparisons among the main studies carried out to date by different authors.