An Investigation of Parallel Post-Laminar Flow through Coarse Granular Porous Media with the Wilkins Equation

Abstract

Behaviour of flow resistance with velocity is still undefined for post-laminar flow through coarse granular media. This can cause considerable errors during flow measurements in situations like rock fill dams, water filters, pumping wells, oil and gas exploration, and so on. Keeping the non-deviating nature of Wilkins coefficients with the hydraulic radius of media in mind, the present study further explores their behaviour to independently varying media size and porosity, subjected to parallel post-laminar flow through granular media. Furthermore, an attempt is made to simulate the post-laminar flow conditions with the help of a Computational Fluid Dynamic (CFD) Model in ANSYS FLUENT, since conducting large-scale experiments are often costly and time-consuming. The model output and the experimental results are found to be in good agreement. Percentage deviations between the experimental and numerical results are found to be in the considerable range. Furthermore, the simulation results are statistically validated with the experimental results using the standard ‘Z-test’. The output from the model advocates the importance and applicability of CFD modelling in understanding post-laminar flow through granular media.

1. Introduction

The porous media flow is commonly characterized by the linear relation between the superficial velocity and hydraulic gradient, proposed by Henry Darcy as:

$$V=ki$$

(1)

where V is the superficial velocity (m/s); i is the hydraulic gradient (head loss per unit length in the direction of flow), and k is the coefficient of permeability (m/s).

However, it is observed that Equation (1) predicts the flow satisfactorily only when the Reynolds number is less than 10 [1,2]. At higher Reynolds numbers, the velocity and hydraulic gradient do not exhibit a linear relationship. Therefore, use of Darcy’s law can cause considerable error in the post-laminar conditions such as the calculation of well productivity [3] and the discharge measurement and design of pumping wells, hydraulic structures such as rock fill dams [4,5], water treatment filters and fissured rocks. Therefore, the characterization and modeling of flow through porous media is pertinent since the issue aptly addresses the abovementioned challenges. Furthermore, the studies and knowledge about the porous media flow can also largely influence efficiency and understanding of the oil and gas exploration which are directly related to the energy sector.

The flow, in the post-laminar regime is mostly described using the Forchheimer equation, as:

i = aV + bV²

(2)

where a (s/m) and b (s²/m²) are the Darcy and Non-Darcy coefficients, respectively.

From a thorough review of the literature, one can find enough evidence confirming the complex variation pattern of the Darcy and non-Darcy coefficients with field and media conditions such as porosity, media size, convergent angle, permeability, and so on [6,7,8,9,10,11]. Efforts are made to relate the Darcy and non-Darcy coefficients with the parameters influencing the flow, but they are observed to be empirical therefore of limited use [12,13,14,15,16,17]. The complex variation pattern of the coefficients makes it difficult to predict their values for a given set of media and field conditions. Therefore, modelling of the post-laminar flow, using the Forchheimer type equations often prove to be very difficult.

Since there is no proper equation for representing flow in the non-laminar regime, the search for a single constitutive relation continues even now [18,19].

The Wilkins equation was developed to represent flow in both laminar and non-laminar flow regimes. The equation can be presented as [20]:

$$V=C_1{\mu }^{\alpha }{r}^{\beta }{i}^{\gamma }$$

(3)

where r is characteristic length; α and β are constants; C₁ is a coefficient; γ is an exponent having a value between 0.5 and 1.0, depending on the nature of the regime; and μ is the dynamic viscosity of the fluid. For a constant viscosity, Equation (3) can be written as:

$$V=W{r}^{\beta }{i}^{\gamma }\mathrm{with},W=C_1{\mu }^{\alpha }$$

(4)

where W is the modified Wilkins coefficient for a constant viscosity. The coefficients of the Wilkins equation are reported to be relatively non-deviating for media with different hydraulic radius (discussed in Section 2) [21,22,23], unlike the Forchheimer equation [6,20]. The non-deviating nature of Wilkins coefficients encouraged the authors to further investigate the behaviour of these coefficients when subjected to independently varying media sizes and porosities. Further, an effort is made in the present study to simulate the post-laminar flow through porous media with a comprehensive CFD model using ANSYS FLUENT. The model outputs are compared with the experimentally obtained results and further validated statistically using the standard ‘Z-test’. The model outputs and its correlation with the experimental results may advocate the use of such CFD models to analyse and understand the nature and characteristics of post-laminar flow through granular porous media.

2. Experiments and Methodology

A specially conceived parallel flow permeameter with a diameter of 250 mm and a length of 1100 mm is used after packing it with irregularly shaped media having a volume diameter (diameter of a sphere having the same volume as the irregular media) of 29.8 mm, 34.78 mm, and 41.59 mm (Figure 1). To investigate the effect of porosity variation on the flow resistance, velocity behaviour; the same media is repacked three times to achieve separate porosities. The media is packed between two perforated plates and water is allowed to flow at the maximum possible rate for 4 to 5 h to make sure no reorientation of media occurs during the experimentation.

Water is supplied through a header tank (1300 mm × 300 mm × 300 mm) to the permeameter at a constant volumetric flow rate (m³/s). The flow rate is measured as the average of three volumetric flow measurements (m³/s) with an accuracy of ±2.35%. The measured volumetric flow rate m³/s) is divided by the cross sectional area (m²) of the permeameter to obtain the superficial velocity (m/s). For every flow rate, eight separate piezometric head differences are noted from tapings placed vertically over the permeameter at a regular interval of 50 mm. The average hydraulic gradient calculated from the recorded head differences is used for further analysis with an accuracy of ±1.89%. Such an arrangement eliminates any error due to the non-uniformity of packing.

To ensure a post laminar flow condition, experimentations are performed at a higher range of Reynolds number (1736–7194) defined by Kovacs [1,2] as:

$$\mathrm{Re}=\frac{V{d}_k}{\nu }$$

(5)

where ν = kinematic viscosity (m²/s); d_k = characteristic length, defined as d(1 − f)α_s/4f; d and α_s are the volume diameter (m) and shape factor of the media, respectively.

The ‘hydraulic radius’ (r) [24] (ratio of the void ratio (e) to specific surface (S₀)) is used as the characteristic length while calculating the Wilkins coefficients since it includes both parameters (porosity and media size) that influence the flow. The average surface areas and volumes which are prerequisites to calculate the specific surface (surface area per unit volume) are measured as described by Banerjee et al. [23] and presented in

Table 1. Properties of the porous media used.

Passing and Retaining Sieve Sizes (cm)	Volume (cm3)	Volume Diameter (cm)	Avg. Surface Area (cm2)	Specific Surface (/cm)
2.50–3.15	14.46	2.98	44.68	2.88
3.15–3.75	23.10	3.48	58.44	2.53
3.75–5.00	38.55	4.16	82.11	2.13

. After measuring the surface area and volume, the volume diameter (diameter of the sphere having the same volume of the media) and the specific surface is calculated.

3. Correction Factors

The experimental setups used in the laboratory are confined in nature, which creates uneven packing in the vicinity of the wall and in the interior of the bed. However, in reality, most of the flows are unconfined in nature. Similarly, during the experimentation, the flow path is assumed to be linear whereas, in reality, it is very tortuous. As a consequence of such constraints, an exact simulation of the field condition in the laboratory is very difficult. In view of the fact that the results obtained in the laboratory cannot be directly applicable to the field, the wall, tortuosity and porosity corrections are applied to the obtained results.

3.1. Wall Correction

The confined boundary wall of the permeameters used in the laboratory causes uneven packing between the media near the walls, which is loosely packed, compared to the media inside. Therefore, the velocity of the flow inside the media and near the wall becomes different. Given that, the media used in the present study are larger in size and the ratio of the permeameter to media diameter is low. Therefore, a wall correction factor developed by Thiruvengadam and Kumar [6] is used to correct the data from any wall effect as follows:

$$V_w=\frac{V}{C_w}\mathrm{with}{C}_w={\left[\frac{(D+4.83\times \frac{d}{2})(D-0.83\times \frac{d}{2})}{D^2}\right]}^{-1}$$

(6)

where V_w is the corrected velocity after wall correction, D is the diameter of the permeameter section (mm), and d is the diameter of the porous media used.

3.2. Tortuosity Correction

Generally, the hydraulic gradient is calculated as the ratio of the total head loss to the distance between the tapings, assuming that the flow path is linear in the direction of the pressure loss. However in reality, the fluid flows through the interlinked void spaces within the media. Owing to that, the flow path is mostly non-linear which can cause inaccuracies in the calculation of hydraulic gradient. Therefore, tortuosity corrections are helpful to reduce that error. By definition, tortuosity (τ) is presented as the ratio between the actual lengths of flow between two tapings to the linear distance between them.

The complexity of the porous structure poses a very challenging task while deriving an expression of tortuosity correction factor. Therefore, most tortuosity corrections are empirical in nature obtained from experiments. The Yu and Li [25] model is used here because it theoretically represents a well-defined relationship between tortuosity and porosity. The model presents the tortuosity correction as:

$$\tau =1+\sqrt{1-f}(\frac{1}{2}+\frac{\sqrt{{(\frac{1}{\sqrt{1-f}}-1)}^2+\frac{1}{4}}}{1-\sqrt{1-f}})$$

(7)

where f is the porosity of the media.

3.3. Porosity Correction

The velocities obtained from the test are calculated by dividing discharge with the cross-section area of the permeameter (superficial velocity), which is different from the velocity inside the media (pore velocity). In order to estimate the actual pore velocity, porosity corrections [26] are introduced as follows:

$$V_v=\frac{V}{f}$$

(8)

where V_ν is the pore velocity and f is the porosity.

The results are analysed after incorporating all the corrections to their respective parameters as:

$$V_c=\frac{V.\tau }{f{C}_w}\mathrm{and}{i}_c=\frac{i}{\tau }$$

(9)

where V_c and i_c are the corrected velocity and hydraulic gradients, respectively.

4. Data Analysis

Before the analysis of the data, the following assumptions were made: the flow is single phase and unidirectional, the medium is homogeneous, and the porosity of the medium is uniform. After applying all the corrections in to the experimentally obtained superficial velocities (V) and hydraulic gradients (i), the corrected velocities (V_c) and corrected hydraulic gradients (i_c) are plotted in Figure 2 and a relationship in the form of a power law is obtained. This type of equations are also referred as Izbash equation [27] or Missbach equation [20] in the literature:

$$i_c=P{V}_c{}^j$$

(10)

Values of the coefficient P and the exponent j from Equation (10) are calculated from Figure 2. The equation can be modified as follows:

$$V_c=(\frac{i_c{}^{(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$j$}\right.)}}{\sqrt[j]{P}})$$

(11)

The Wilkins equation for a constant viscosity can be written as:

$$V_c=W{r}^{\beta }{i}_c{}^{\gamma }$$

(12)

Thus, judging by the similarity of Equations (11) and (12):

$$\sqrt[j]{P}=\frac{1}{W{r}^{\beta }}\mathrm{with}\gamma =\frac{1}{j}$$

(13)

$$\frac{1}{j}\log P=-(\log W+\beta \log r)$$

(14)

When Equation (14) is plotted in linear form (Figure 3) with −(1/j)log P on the y axis and log r on the x axis, the slope represents the value of the coefficient β and the intercept represents the value of log W. The obtained values of W and β are compared for different media sizes and porosities to understand their effect on the Wilkins coefficients.

5. The Behavior of the Wilkins Coefficients with Independently Varying Media Size and Porosity

Experimentally obtained values for the Wilkins coefficients (W and β) are presented in

Table 2. The values of the Wilkins coefficient with variation in media size.

Media Size (mm)	Porosity (%)	W (m-s)	β	γ
29.80	40.59	6.15	0.39	0.54
	43.34
	45.69
34.78	41.72	5.52	0.38	0.56
	44.70
	46.34
41.59	41.03	5.55	0.38	0.55
	43.62
	46.15

. The values of the coefficients are found to be rather constant for all the media experimented with. However, values of γ depend on the flow regime, as mentioned earlier. When compared with the results obtained from earlier reported studies (

Table 3. The values of the Wilkins coefficients as reported by earlier researchers.

Proposed by	Media	Volume Diameter (mm)	Porosity (%)	W (m/s)	β	γ
Wilkins (1956) [21]	Crushed stone	51.00	40.00	5.24	0.50	0.54
Garga et al. (1990) [22]	Crushed stone	24.60	47.00	5.39	0.50	0.53
Pradeep Kumar (1994) [27]	Crushed stone	13.10	47.00	4.94	0.51	0.52
		20.10	45.88
		28.90	48.73
		39.50	48.26

), almost similar values of the coefficients are observed. However, it is worth mentioning that the coefficients in the present study vary a little from the ones presented by Banerjee et al. [23]. This may be due to the application of the correction factors. To compare the present experimental results with some of the earlier reported data [28], correction factors are directly applied to the measured velocity and hydraulic gradient instead of modifying the Wilkins equation itself. Since some of the studies have used correction factors where others did not, and the minor differences in the presented results can be attributed to application of correction factors and their different values used to correct the velocity and hydraulic gradient along with the uncertainty and complexity of the porous structure and the human errors during experimentation and calculation. However, one can agree that, for all the reported data, the coefficients of the Wilkins equation present similar values. Furthermore, an attempt is made in the present study to understand the nature of the Wilkins coefficient (W) when the size of the media (volume diameter) and porosity vary independently, after having the coefficients calculated for different hydraulic radius. The values of W are calculated from Equation (12) and plotted for different media sizes whilst keeping the porosity of the packing constant (Figure 4a) and similarly, for different porosities whilst keeping the media size constant (Figure 4b). Figure 4a,b indicate a constant nature of W for media size variation, whereas it shows a slightly increasing trend with porosity when media size is constant.

6. Numerical Modelling

The CFD approach is a very powerful tool for modelling and visualising flows under different conditions. The CFD approach numerically solves three basic equations: the continuity (mass conservation) equation, the motion equation and the energy equation at each point of the computational domain. The ANSYS FLUENT CFD solver, which is based on the Finite Volume method, is used to create and simulate laboratory conditions and to solve the CFD problem.

6.1. Turbulence Modelling

The continuity equation and the momentum conservation (Navier–Stokes) equation for incompressible flow are presented as Equations (15) and (16) respectively:

$$\frac{\partial u_i}{\partial x_j}=0$$

(15)

where u_i is the velocity component in the x_i direction:

$$\rho u_j\frac{\partial u_i}{\partial x_j}=\frac{\partial p}{\partial x_j}+\frac{\partial }{\partial x_j}(\mu \frac{\partial u_i}{\partial x_j})+\rho g_i$$

(16)

where ρ is the density of the fluid, p the pressure, μ is the molecular viscosity and g is the gravitational acceleration.

The turbulent flow is usually simulated after modification of the continuity equation and the momentum conservation (Navier–Stokes) equation for laminar flow with Reynolds decomposition.

The generated quantities known as “Reynolds Stresses” ($-\rho \overline{u_i^{\prime }{u}_j^{\prime }}$) can be defined by the Boussinesq hypothesis using turbulent viscosity (μ_t) and turbulent kinetic energy (k) [29].

Finally, two additional transport equations need to be incorporated in order to close the equation system representing the turbulent kinetic energy and the turbulence dissipation rate (the k-ε model). The “realizable” k-ε model [30] with a standard wall function is used after a bibliographical study concerning the validation of the CFD model including the realizable k-ε model for a wide range of flows, as well as flows that include the boundary layers and flows over obstacles [30,31].

The two equations for turbulent modelling kinetic energy (k) and turbulence dissipation rate ($\epsilon$) can be written as [29]:

$$\rho [\frac{\partial }{\partial x_j}(k{u}_j)]=\frac{\partial }{\partial x_j}[(\mu +\frac{{\mu }_t}{{\sigma }_k})\frac{\partial k}{\partial x_j}]+G_k-\rho \epsilon$$

(17)

$$\rho [\frac{\partial }{\partial x_j}(\epsilon u_j)]=\frac{\partial }{\partial x_j}[(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }})\frac{\partial \epsilon }{\partial x_j}]+\rho C_{1}S\epsilon -\rho C_2\frac{{\epsilon }^2}{k+\sqrt{\nu \epsilon }}$$

(18)

where $C_1=\max [0.43;\frac{\eta }{\eta +5}]$; $C_2$ = 1.9; ${\sigma }_k$ = 1; ${\sigma }_{\epsilon }$ = 1.2; $\eta$ is the time scale ratio defined as $=S\frac{k}{\epsilon }$ with S is the modulus of the mean rate of strain tensor defined by $S=\sqrt{2S_{ij}{S}_{ij}}$; $G_k$ represents the production of the turbulence kinetic energy.

It is worth mentioning that the authors have opted for the standard values of the model constants in the present study since the required model outputs (pressure loss and superficial velocity) are found to be in judicious agreement with the experimental ones. However, the standard values of the model constants may not always provide satisfactory output especially while modelling with no or little turbulence. In such cases, some low Reynolds number closure [32,33,34,35] may be used.

6.2. Model Description

The 3D geometric model with an inlet section, a tank, a parallel flow permeameter, and an outlet section of similar dimensions as the experimental set-up is created using the geometry feature in the ANSYS FLUENT workbench. Owing to the high degree of randomness and the fact that no two porous media shows similar structures at the pore level, it is unnecessary to construct the whole porous media at the fine level as it may not affect the overall physical property of the media. Therefore, the permeameter section is identified as the porous zone with uniform porosity similar to that achieved during experimentation. The porous media is numerically accounted as a momentum sink or a resistance in the momentum equation composed of viscous and inertial terms as given in Equation (19) [29]:

$$S_i=-({\displaystyle \sum _{j=1}^3{D}_{ij}\mu v_j+{\displaystyle \sum _{j=1}^3{C}_{ij}\frac{1}{2}\rho \left|v\right|v_j)}}$$

(19)

where $S_i$ is the source term for i^th (x, y or z) momentum equation; $\left|v\right|$ is the magnitude of the velocity; D and C are the prescribed matrices. For a homogeneous media:

$$S_i=-(\frac{\mu }{\alpha }{v}_i+C_2\frac{1}{2}\rho \left|v\right|v_i)$$

(20)

where α is the permeability; C₂ is the inertial resistance factor simply specifying D and C as diagonal matrices; and μ and ρ are the dynamic viscosity and the density of the fluid, respectively. The viscous and inertial resistance are measured based on the pressure loss observed from the experimental setup [29,36]. Furthermore, the model is solved by the finite volume method with the velocity pressure coupling done by the SIMPLE algorithm. The governing equations are discretized using the second order upwind differencing scheme. Convergence is attained when the scaled residuals are less than 10⁻⁴ times of their initial values.

6.3. Meshing

By creating the mesh, a continuous domain describing the flow with partial differential equations is replaced with a finite number of volumes (meshes). The discretization is performed using structured tetrahedral meshes. Nine different grid sizes are selected arbitrarily having element number ranging from 8867 to 6,463,169. Iterative convergence is obtained for each of these grid sizes whilst keeping all other input conditions identical. The output parameters, hydraulic gradient, and velocity are compared for all these meshes which have different numbers of elements. Figure 5 shows that the output parameters become constants as the grid size is refined. Furthermore, a grid sensitivity analysis is performed using the Grid Convergence Index (GCI) method to authenticate the precision of the numerical model [37].

The GCI is the representation of the discretization error between the numerical solutions obtained from the finer grid and the coarser one. For three different selected grids (h₁, h₂, h₃) simulations are carried out to estimate the key variables ($\varphi$); in the present study these variables are Hydraulic gradient, Velocity (m/s) at 0.285 m and Velocity (m/s) at 1.11 m. The apparent order (p) of the method is then calculated using the following equations:

$$p=\frac{1}{\ln (r_{21})}\left|\ln \left|{\epsilon }_{32}/{\epsilon }_{21}\right|+q(p)\right|$$

(21)

$$q(p)=\ln (\frac{r_{21}^p-[1.\mathrm{sgn}({\epsilon }_{32}/{\epsilon }_{21})]}{32-[1.\mathrm{sgn}({\epsilon }_{32}/{\epsilon }_{21})]})$$

(22)

where ${\epsilon }_{32}={\varphi }_3-{\varphi }_2$, ${\epsilon }_{21}={\varphi }_2-{\varphi }_1$, and $r_{21}=h_2/h_1$. The equations are solved using fixed-point iteration using the first term as an initial guess [37]. The relative error ($e_{\alpha }^{21}$, $e_{\alpha }^{32}$) and the extrapolated relative error ($e_{ext}^{21}$, $e_{ext}^{32}$) are calculated using Equations (23) and (24):

$$e_a^{21}=\frac{{\varphi }_1-{\varphi }_2}{{\varphi }_1}$$

(23)

$$e_{ext}^{21}=\frac{{\varphi }_{ext}^{21}-{\varphi }_2}{{\varphi }_{ext}^{21}}\mathrm{where}{\varphi }_{ext}^{21}=(r_{21}^p{\varphi }_1-{\varphi }_2)/(r_{21}^p-1)$$

(24)

Finally, the fine Grid Convergence Index (GCI) is calculated using the following equation:

$$GC{I}_{fine}^{21}=\frac{1.25e_a^{21}}{r_{21}^p-1}$$

(25)

Table 4. Sample calculations of numerical uncertainty using the GCI method [37].

Parameter	Hydraulic Gradient	Velocity (m/s) at 0.285 m	Velocity (m/s) at 1.11 m
h1	0.005216	0.005216	0.005216
h2	0.006937	0.006937	0.006937
h3	0.009226	0.009226	0.009226
r21	1.33	1.33	1.33
r32	1.33	1.33	1.33
${\varphi }_1$	0.01992	0.01083	0.01156
${\varphi }_2$	0.01984	0.01080	0.01159
${\varphi }_3$	0.01973	0.01077	0.01096
P	1.32772	1.30788	10.5924
${\varphi }^{21}{}_{ext}$	0.02010	0.01088	0.01156
${\varphi }^{32}{}_{ext}$	0.02010	0.01088	0.01163
$e_{\alpha }^{21}$	0.40%	0.20%	0.27%
$e_{\alpha }^{32}$	0.59%	0.29%	5.45%
$e_{ext}^{21}$	0.86%	0.45%	0.01%
$e_{ext}^{32}$	1.26%	0.65%	0.28%
$GC{I}_{fine}^{21}$	1.09%	0.56%	0.02%

represents values of these parameters for three selected grids (227,414; 72,709; and 32,301). The output parameters are mentioned as ${\varphi }_1$, ${\varphi }_2$, ${\varphi }_3$ or grid sizes h₁, h₂ and h₃ in Table 4. A similar analysis is performed for other selected grid sizes. Numerical uncertainty in the fine grid solution is observed to be 1.09%, 0.56%, and 0.02% respectively for the hydraulic gradient, the velocity at a vertical distance of 0.285 m from the entrance of the permeameter, and the velocity at a vertical distance of 1.11 m from the entrance of the permeameter (Table 4). The GCI values indicate a negligible discretization error and therefore a grid independent solution for the selected grid. Considering the simulation time and precision, a grid system of 1,288,576 elements is adopted for the model.

6.4. Boundary Conditions

The boundary conditions are very important for the numerical solution of the problem. The type and the numerical values of the boundary conditions are carefully chosen. The inlet velocity is calculated by dividing the experimental discharge with the area of the inlet pipe. The turbulent quantities (k and ε) at the inlet are calculated from the Equations (26) and (28) using the turbulent intensity (I); turbulent length scale (l),which depends on the hydraulic diameter (D_h); and the inlet velocity(V) [38]. Equation (26) is used to calculate the turbulent intensity at the core of a fully developed duct [29]:

$$I=0.16{(\mathrm{Re})}^{-(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$8$}\right.)}$$

(26)

The value of k is then calculated using the turbulent intensity from the following equation:

$$k=\frac{3}{2}{(V\cdot I)}^2$$

(27)

where Re is the Reynolds number for a pipe flow defined as $\frac{\rho Vd}{\mu }$.

$$\epsilon =C_{\mu }^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}\frac{k^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{l}\mathrm{with}l=0.07D_h$$

(28)

The outlet is kept as an outflow boundary, assuming that the flow is completely developed; thereby causing the diffusion flux for all flow variables in the exit direction to be zero [29].

7. Results and Discussion

The values of the pressure and velocity along the direction of the flow (length of the permeameter) for the media sizes of 29.8, 34.78 and 41.59 mm are presented in Figure 6a,b, Figure 7a,b and Figure 8a,b, respectively. It can be observed that the total pressure decreases as the flow passes through the length of the permeameter, whereas the velocity plots illustrate a constant superficial velocity throughout the length of permeameter. The pressures and velocities at the tank section (−0.3 m to 0 m) are not used for further analysis since only the permeameter section (0 m to +1.1 m) is identified as a porous zone and subjected to the study. It is worth mentioning that, during the simulation, the operating pressure is set to zero in order to cut down the rounding errors, resulting in negative values in the pressure distribution diagram. However, it does not affect the simulation because of two main reasons; (a) the flow is assumed to be incompressible and (b) the pressure difference is taken into account by the Navier–Stokes equation which drives the flow. The validation of the results and the importance of the model are discussed further in the following sections.

7.1. Comparison between the Experimental and Simulation Data and Statistical Validation of the Simulation

Experimentally and numerically calculated hydraulic gradients are plotted as exponential function (Equation (10)) of the respective velocities in Figure 9a–c, Figure 10a–c and Figure 11a–c. The plots suggest an acceptable correlation between the experimental results and the simulation results for the range of the data tested.

Figure 12a–c represent the percentage deviations between experimental and simulation results applied to the different velocities for all media sizes and porosities. Simulation results are found to be in good agreement with the experimental results for the range of velocities presented. However, the deviation between the simulation and experimental results is not constant. The deviation reduces to a minimum value at a certain velocity and then further increases (Figure 12). This is due to the change in the flow regime with the variation of the velocity. With the change in velocity, the magnitudes of the resistive forces affecting the flow differ, causing a variation in the total resistance thereafter in the velocity-flow resistance relationship. This finally results in a great deviation between the simulation and experimental results. Furthermore, the simulation results have been validated statistically with the experimental results using the standard “Z-test” [39,40], which is explained further in the next section.

7.2. Validation of the Simulated Data Using the Z-Test

In order to test the validity of the simulation result, first the hypothesis is introduced [39,40,41] as follows:

$$\begin{array}{l}{H}_0:{\mu }_1={\mu }_2\\ H_1:{\mu }_1\ne {\mu }_2\end{array}$$

(29)

where μ₁ and μ₂ are the population mean of the experimental and simulation result. In order to validate the simulated result, the null hypothesis, given by H₀, must be accepted.

The level of significance (α) is taken to be 5%, which implies that a confidence level of 95% is considered [41]. As the present study deals with two-tailed hypothesis, α = 0.05/2 = 0.025. The corresponding critical value of $Z$ for α = 0.025 is ±1.96 (from the Z table). Therefore, for the hypothesis given as Equation (29) to be accepted, the calculated Z value should be within the interval of [−1.96, 1.96]. The calculated Z values are obtained by employing the Z-test, as given by Equation (30):

$$Z_{calculated}=\frac{(\overline{x_1}-\overline{x_2})-({\mu }_1-{\mu }_2)}{\sqrt{\frac{{\sigma }_1^2}{N_D}+\frac{{\sigma }_2^2}{N_D}}}\mathrm{for}{N}_D\ge 30$$

(30)

where $\overline{x_1}$, $\overline{x_2}$ and ${\mu }_1$, ${\mu }_2$ are the sample and population means for the experimental and simulation results. N_D is the total number of data points (sample size) and σ₁, σ₂ are the standard deviations for the experimental and simulation results.

The results of the Z-test for the three different media sizes are presented in

Table 5. Values from the Z-test for different media sizes and porosities.

Media Size (mm)	Porosity (%)	Range of Velocity (m/s)	$\overline{{\mathit{x}}_1}$	$\overline{{\mathit{x}}_2}$	${\mathit{\sigma }}_1$	${\mathit{\sigma }}_2$	${\mathit{N}}_{\mathit{D}}$	Z Value
29.80	40.59	0.01–0.757	0.643	0.599	0.540	0.359	250	−0.940
	43.34	0.01–0.757	0.586	0.606	0.497	0.517	250	0.429
	45.69	0.01–0.757	0.554	0.535	0.464	0.446	250	−0.459
34.78	41.72	0.01–0.757	0.565	0.554	0.463	0.452	250	−0.276
	44.70	0.01–0.757	0.527	0.510	0.435	0.420	250	−0.425
	46.34	0.01–0.757	0.496	0.471	0.407	0.383	250	−0.694
41.59	41.03	0.01–0.757	0.511	0.489	0.422	0.406	250	−0.574
	43.62	0.01–0.757	0.464	0.440	0.391	0.368	250	−0.684
	46.15	0.01–0.757	0.439	0.418	0.369	0.349	250	−0.661

. The calculated values of Z for all three media sizes are found to be within the acceptable range [−1.96, 1.96] (Figure 13a–c). Therefore, the null hypothesis H₀ is accepted, validating the simulation results with the experimental results. Result from the statistical analysis along with the percentage deviation between the simulation and experimental data signify the accuracy of the CFD model used.

This type of CFD model can be very useful to predict the nature of the flow for any hydraulic structure. For a given media size and porosity, this type of simulation can provide a complete velocity and pressure loss profile throughout any given hydraulic structure which can immensely aid the designers and engineers. Furthermore, the obtained velocity and pressure loss data from such models can also be used to analyse the behaviour of different post-laminar flow equations used in porous media flow including the Wilkins equation.

8. Conclusions

A number of non-linear equations were proposed to predict post-laminar flow through porous media at a macroscopic level. However, the complex variation pattern of their respective coefficients limits their applicability in field conditions. The Wilkins equation is studied both theoretically and experimentally with respect to parallel flow. The obtained results are presented as follows:

• The Wilkins equation can be satisfactorily used to represent post-laminar flow through porous media.
• The Wilkins coefficients are found to have a non-deviating nature with varying hydraulic radius. The obtained results from the present study are similar to the results reported in the literature.
• When subjected to variation in media size, the coefficients of the Wilkins equation are constant, given that the porosity is constant. However, variations in the porosity result in small variations of the coefficient W.
• The flow condition inside the experimental set up is simulated with a CFD model in the ANSYS FLUENT software. Trends similar to the experimental ones are obtained from the simulation results. The percentage deviation between the simulation and experimental results are within the acceptable range.
• For further validation, the experimental results are statistically compared with the simulation results using the standard Z-test. The values of Z calculated are found to be within the acceptable region for all the experimental results.

Finally, the experimental and simulation results in the present study conclude that the Wilkins equation can represent non-linear flow through porous media with a satisfactory accuracy, compared to Darcy’s equation. The coefficients of the Wilkins equations are found to be constant for variations of porosity and media size, unlike the Forchheimer type equations. Hence, this equation can be regarded as a convenient tool for designing and measuring the discharge for parallel flow through porous media in the non-laminar regimes such as rock fill dams, water filters, aquifers, water, oil and gas wells, and so on.