The use of the inverse problem methodology in analysis of fluid flow through granular beds with non-uniform grain sizes DOI: 10.31648/ts.6772

The pressure drop during water flow through two gravel beds with 2-8 and 8-16 [mm] grain size was measured across a wide range of filtration velocities, and the optimal method for calculating the coefficients for Darcy’s law and Forchheimer’s law was selected. The laws and the experimental data were used to develop a computational program based on the Finite Element Method (FEM). The results were compared, and errors were analyzed to determine which law better describes flow data. Various methods of measuring porosity and average grain diameter, representative of the sample, were analyzed. The data were used to determine the limits of applicability of both laws. The study was motivated by the observation that computational formulas in the literature produce results that differ by several orders of magnitude, which significantly compromises their applicability. The present study is a continuation of our previous research into artificial granular materials with similarly sized particles. In our previous work, the results produced by analytical and numerical models were highly consistent with the experimental data. The aim of this study was to determine whether the inverse problem methodology can deliver equally reliable results in natural materials composed of large particles. The experimental data were presented in detail to facilitate the replication, reproduction and verification of all analyses and calculations.


Introduction
Fluid flow through porous media can be described mathematically at two qualitatively different levels. The first level (historical) involves macroscopic measurements. In this approach, a porous medium is regarded as a homogeneous medium where flow resistance is generally averaged in space and time (ERGUN 1952, HELLSTRÖM, LUNDSTRÖM 2006, SIDIROPOULOU et al. 2007, SOBIESKI 2010. The structure of the material, channel shape and particle surface are not considered. Flow resistance is determined with the use of generalized parameters such as average particle diameter, porosity and tortuosity. The main disadvantage of this approach is that constant values, which usually differ for various types of porous media, have to be incorporated into mathematical models. As a result, a generally applicable mathematical model for measuring flow resistance in porous media has never been developed. Generalized Forchheimer's law appears to be the only relatively universal equation (SOBIESKI, TRYKOZKO 2011), but two coefficients have to be defined individually for every medium (SIDIROPOULOU et al. 2007). A dedicated methodology for calculating these coefficients has never been proposed, and dozens of formulas generating results that differ by many orders of magnitude have been described in the literature (SOBIESKI, TRYKOZKO 2011). The second level involves microscopic measurements where the shape of channels and local factors are taken into consideration. This approach appears to be consistent with the laws of physics, but it requires highly complex numerical models and is very difficult to implement in practice. The main disadvantage of this approach is that microscopic observations cannot be expressed on the macroscopic scale. In practice, numerical grids can be generated to define a small fragment of space and calculate pressure fields and velocity with the use of CFD techniques. However, numerical grids describing geometrically complex pores are difficult to generate and require considerable computational power, which is why this approach cannot be used to develop full-size models of real objects. Numerous simplifications have been applied to overcome these problems. Problems are analyzed in 2D space, and pore space is modeled with simple geometric shapes, including circles (HELLSTRÖM, LUNDSTRÖM 2006, PESZYŃSKA et al. 2009a, PESZYŃSKA et al. 2009b) and rectangles (KOPONNEN et al. 1996, KOPONEN et al. 1997, MATYKA et al. 2008. Pseudo-random methods are also used (NABOVATI, SOUSA 2007). Pore space is very rarely modeled in 3D space based on the actual structure of the porous medium. Such an attempt was made by CARMINATI (2007) who captured the real geometry of pore space by computed tomography and used the results to develop a classical CFD model. In other studies, the structure of porous space was described at the microscopic scale, but fluid flow was not modeled (MOTA et al. 1999, NEETHIRAJAN et al. 2006).
It should be noted that both levels can be combined by multiscale (hybrid) modeling. In this approach, flow is described with macroscale models, but the required parameters and coefficients, such as tortuosity, are determined analytically based on a simplified geometry of pore space. This approach has been rarely described in the literature, and notable examples include a study by WU et al. (2008) which presents a full mathematical model, or a study by YU and LI (2004) who proposed an analytical method for calculating tortuosity in geometrically simplified pore space. The multiscale approach can also be used in full CFD models (HELLSTRÖM, LUNDSTRÖM 2006, PESZYŃSKA et al. 2009a, PESZYŃSKA et al. 2009b) but significant progress in this area has not been made to date.
The current study relies on the macroscopic approach, but the existing empirical formulas for calculating the coefficients for Darcy's law and Forchheimer's law were not used. Empirical formulas produce highly varied results, which is why the proposed methodology was based on the inverse problem ( VAN BATTENBURG, MILTON-TAYLER 2005, HUANG, AYOUB 2008, SOBIESKI, TRYKOZKO 2011, 2014a, 2014b. The method of processing experimental data plays an important role in this approach (SOBIESKI, TRYKOZKO 2011). The proposed methodology involves a simple laboratory experiment, but it guarantees high consistency of experimental data with the analytical or numerical model across a wide range of filtration velocities. In the work of SOBIESKI and TRYKOZKO (2011), data consistency was estimated at 1-3% across the tested range of filtration velocities. Similar results were reported by SOBIESKI and TRYKOZKO (2014b), but relative error was higher for several initial measurement points (with the lowest filtration velocity), and it reached 10.25% for the lowest filtration velocity (relative error decreased rapidly for successive filtration velocities). The greatest advantage of the inverse methodology is that pressure loss can be reliably predicted when the value of the permeability coefficient and/or the Forchheimer coefficient for a given porous medium is known. In a cross-validation procedure, an experiment involving the determination of the above coefficients can be regarded as an element of the training set, and a system where flow resistance is modeled can be regarded as the validation set. The training model is one-dimensional and simple to develop, whereas predictive models can be two-or threedimensional, at least for isotropic media.
Porous media composed of glass beads with roughly identical diameters were investigated by SOBIESKI and TRYKOZKO (2011, 2014a, 2014b. Porosity and pressure loss during fluid flow were highly similar in the investigated media. These results could suggest that in less homogeneous media, experimental data will be less consistent with numerical data. To determine the magnitude of relative error in such cases, these experiments were replicated with the use of a different porous medium. The aim of this study was to determine whether the inverse problem method can produce equally reliable results in natural materials composed of large particles. The applied tools and methods are not new. However, the presented results were obtained in a new medium, which contributes to the body of knowledge on the applicability of the inverse problem method for describing the properties of porous media with variable particle size.

Experimental design
The test stand for analyzing the porous medium is presented in Figure 1. The stand is composed of a plexiglass tube filled with gravel (6). Water from a bottom tank (1) is supplied to the bottom of the plexiglass cylinder by a pump (2). The flow rate is controlled by a control valve (3) and an overflow valve (4). Water flows through gravel to the top tank (7), and it returns to the bottom tank via an overflow pipe (8). Volumetric flow rate was measured with a rotameter (5). Pressure was measured with U-tube manometers joined to pipe connectors (marked 1 to 4 in the diagram). Water temperature was measured with a thermometer (9). Measurement errors were as follows: volumetric flow rate 0.000000056 m 3 /s, piezometric head 1.0 mm H2O, temperature 0.1 K. The distance between extreme measurement points was 0.9 m, and the porous medium had a crosssectional area of 0.005 m 2 . Two porous media were tested: gravel with 8-16 mm grain size and gravel with 2-8 mm grain size.
All measurements were conducted in the test site presented in Figure 1. The correlation between volumetric flow rate and piezometric head was determined between extreme measurement points (4-1). The drop in pressure was measured at different settings of the control valve (3). In the porous medium with 2-8 mm grain size, the control valve was set to 12 volumetric flow rates from 0.007 m 3 /s to 0.188 m 3 /s. In the porous medium with 8-16 mm grain size, the control valve was set to 7 volumetric flow rates from 0.137 m 3 /s to 0.188 m 3 /s. In the porous medium with 8-16 mm grain size, the range of valve settings was small because differences in the pressure of water flowing through the medium were negligible. Flow resistance was very low in the tested porous medium. The measurements were performed in triplicate for each medium. The cylinder was then emptied, gravel was stirred and poured back into the cylinder. This operation was repeated three times to calculate the average pressure difference Δℎ. The measurements were conducted at constant water temperature of 304 [K], at which water density was = 995.34 kg/m 3 , and water viscosity was Filtration velocity is calculated with the following continuity equation:
Darcy's law can also be expressed by the following formula: where: -decrease in hydraulic head [-], calculated as the piezometric gradient between two points (in this case, the average value from several measurements).
The filtration coefficient and the permeability coefficient [m 2 ] are bound by the following correlation: The results of the measurements in the model based on Darcy's law are presented in Tables 3 and  4. The average values of the permeability coefficient for porous media with 2-8 mm and 8-16 mm grain size were determined at 2.16E-09 and 3.9938E-08 [m 2 ], respectively. The average values of the filtration coefficient for the above media were determined at 0.026892 m/s and 0.497196 m/s, respectively. The pressure drop in the compared media is presented in Figure 3.

Mathematical model based on Forchheimer's law
Flow resistance in a porous medium is also described with Forchheimer's law (1901). The Forchheimer equation contains a non-linear term for adjusting the result in high-flow rate environments. The applicability limits of both laws are discussed in successive chapters of the study. The Forchheimer equation has the following form: where: -Forchheimer coefficient (non-Darcy coefficient,  factor) [1/m]. If = 0, the flow regime is described by Darcy's law.
To determine parameters and based on the experimental data, equation (5) was transformed as follows (HUANG, AYOUB 2008): The following expressions were introduced: to produce a linear equation: Coefficients and were determined with the use of equation (7), and the coefficients in equation (8) were determined by least squares approximation. The results are shown in Figures 4 and 5. These coefficients are equivalent to the corresponding terms in equation (5). In the media with 2-8 mm and 8-16 mm grain size, the value of coefficient 1 ⁄ was determined at 0.33303E+09 and -0.758964E+07, respectively, and the value of coefficient was determined at 0.221048E+05 and 0.578304E+04, respectively. When the assumption of conservation of momentum modeled by Darcy's law is met: In equation (10), the filtration coefficient is a generalized tensor (BREUGEM et al. 2004). Equation (10) is a generalization of formula (3) for multidimensional flow.
The following boundary conditions were imposed: the Neumann condition is zero along the cylinder, volumetric flow rate was imposed as the Neumann condition at the cylinder inlet, and piezometric head was imposed as the boundary condition at the cylinder outlet. It was assumed that the reference level = 0 corresponds to the cylinder inlet (bottom); therefore, piezometric head ∆ℎ ̅̅̅̅ at the outlet was equal to the anticipated pressure, increased by the cylinder's total length.
If pressure equals zero at the outlet, piezometric head at the outlet will be equal to the cylinder's length.
The flow equation (9) was solved with the Fine Element Method (FEM) with linear basis functions (LUCQUIN, PIRONNEAU 1998). A cross-section of the computational domain and its elements is presented in Figure 6. The model was build based on a discrete mesh that had been originally developed in GAMBIT (2004) and converted to the format used in the designed computational program. Figure 6. Segment of the numerical mesh used in the designed computational program A system of linear equations generated by the numerical mesh was solved with the conjugate gradient method and incomplete Cholesky preconditioning (KAASSCHIETER 1998).
The Forchheimer equation was modified for use in the designed program based on the solution presented by FOURAR et al. (2005). An iterative algorithm is based on the effective permeability coefficient. Successive approximations of the effective permeability coefficient were determined with the use of equation (5): where the upper index denotes velocity fields in successive iterations. Experimentally measured pressure was imposed as the boundary condition, and boundary conditions remained unchanged during iteration. The flow equation (11) was solved for (0) = to approximate the velocity field ( ) . The iterative scheme was as follows: the following approximation of the permeability coefficient ( ) was derived with equation (11), and the corresponding velocity field was calculated. In each case, the calculated velocity field was compared with the field approximated in the previous iteration. Iterations were repeated until convergence was reached. This is not an optimal approach, but it offers an easy solution, including in non-linear models.

Results of numerical simulations
The results of the simulations conducted in the compared media are presented in Figure 7. The numbers in the diagram correspond to the measurements presented in Tables 1 and 2. In both cases, the pressure drops estimated with Darcy's laws differ from the experimental data. where is the numerically determined value and 0 is the exact (measured) value, are presented in Tables 5 and 6 as well as in Figure 8.

Determination of the porosity coefficient
The porosity coefficient was determined with the use of two graduated cylinders with 1 [ml] grading divisions. The first cylinder was filled with porous material with a known volume , and the second cylinder was filled with water. Total skeletal volume can be determined when the volume of porous material and the volume of water required to completely fill the porous part of the medium are known: = − .
In the porous medium with 2-8 mm grain size, the material "settled" under the influence of water, which decreased its total volume. The medium's volume decreased as flowing water transported the smallest grains which filled empty spaces between large grains. To account for the above observation, the porosity coefficient was measured twice. In the first measurement, 0.0015 m 3 of gravel was loosely poured into the cylinder and covered with water, and pore volume was measured. In the second measurement, 0.002 m 3 of gravel was lightly compacted in the cylinder, covered with water, and pore volume was measured. The "settling" of gravel particles was not observed in the medium with 8-16 mm grain size. Average porosity was determined at 8−16 = 0.411, 2−8 max = 0.358 and 2−8 min = 0.294. Results of measurements are available in Tables 7, 8 and 9.

Determination of the equivalent diameter of gravel grains
The average diameter was determined by counting the number of grains s n and measuring their total volume in a control volume . The average volume of one grain was calculated by dividing total grain volume by the number of grains. The formula for calculating the volume of a sphere was used to determine the average grain diameter [m]. The control volume was 6E-04 [m 3 ] for the medium with 8-16 mm grain size and 15E-06 [m 3 ] for the medium with 2-8 mm grain size.
The measured values of average grain diameter seem to be correct. This parameter was determined at 14.7 mm in the medium with 8-16 mm grain size, and 3.3 mm in the medium with 2-8 mm grain size. The obtained values were used in alternative calculations of the porosity coefficient. The results were highly similar to the previously obtained values. The results of the measurements and calculations performed for both media are presented in Tables 10 and 11. Grain "settlement" was negligible in the medium with 2-8 mm grain size due to small control volume, and this observation was not taken into consideration.

Reynolds number
According to the literature, Darcy's law and Forchheimer's law have certain limits of applicability, as illustrated in Figure 9. The Reynolds number is the main criterion for determining the validity of both laws. The Reynolds number is defined as (NAŁĘCZ 1991, SAWICKI et al. 2004): where is the grain diameter [m] in an ideally homogeneous medium composed of spherical particles that offer similar resistance as the real material.
In the literature, considerable variations are noted in the upper limit of applicability of Darcy's law and Forchheimer's law, for example: = 3 − 10 (BEAR 1972), = 1 − 15 (HASSANIZADEH, GRAY (1987); the suggested value is Re = 10),   Figure 10. Reynolds numbers were also calculated for extreme grain diameters in the compared media to present the full range of changes in the evaluated parameter. In the medium with 8-16 [mm] grain size, flow always occurred in a region that was separated by a considerable distance from the upper limit of applicability of Darcy's law and Forchheimer's law (regardless of the fact which limit value was adopted). In this medium, fluid flow should not be modeled with Darcy's law. This conclusion had been initially formulated in the error analysis, and it was further confirmed by the above observations. In fine-grained gravel, measurements began at much lower filtration velocity, and based on literature data, the first two or three measurements could be classified as Darcy flow. However, high percentage error does not confirm the above observation. These findings suggest that in forced-flow systems, such as the presented test stand, flow should be generally modeled with the Forchheimer equation or similar non-linear equations. The above conclusion confirms that the presented Reynolds numbers were calculated based on equivalent diameters, which does not guarantee that a given Reynolds number will describe all regions of a porous medium. Due to possible deviations, Forchheimer's law should be applied in this case because it applies to a wider range of velocities.

Conclusions
The results of the study supported the formulation of the following conclusions: • Due to the absence of strict criteria for determining the validity of the basic laws for modeling flow in porous media, fluid flow in a gravel medium with 2-16 [mm] grain size should be generally described with the mire universal Forchheimer's law which applies to a wider range of velocities. Forchheimer's law can probably be also used to model flow in gravel and sand with different particle size fractions. The above is confirmed by Figure 7 which clearly indicate that unlike Darcy's law, Forchheimer's law guarantees high consistency of modeled data with experimental data across the entire range of filtration velocity values. • The results of this study confirm that the upper limit of applicability of Darcy's law is • Forchheimer's law should be regarded as a broader version of Darcy's law, rather than its extension for flows with higher Reynolds numbers. Published sources (for example in Fig.  9) can be somewhat misleading by suggesting that Darcy's law should be applied to lower values of the Reynolds number, whereas Forchheimer's law should be applied to higher values of the Reynolds number. The results presented in Tables 5 and 6 and in Figure 8 indicate that Forchheimer's law can be used in all cases where Darcy's law could be applied. In these cases, the non-linear term of the Forchheimer equation will simply lose its significance. • Flow resistance in porous media can be reliably modeled based on the results of a simple "training" experiment which produces accurate values of the coefficients for Forchheimer's law. These coefficients should be calculated by measuring pressure drops across the widest possible range of filtration velocities. The wider the velocity range and the more accurate the measurements, the smaller the discrepancy between simulation results and experimental data.
The presented findings indicate that flow resistance in sand and gravel (and probably many other porous media) should be modeled based on Forchheimer's law and the results of a preliminary "training" experiment. The application of the inverse problem methodology significantly increases the reliability of the modeled data.