Identification of the pore size distribution of a porous medium by yield stress fluids using Herschel-Bulkley model

. In this paper, we present a new method to determine the pore-size distribution (PSD) in a porous medium. This innovative technique uses the rheological properties of non-Newtonian yield stress fluids flowing through the porous sample. In a first approach, the capillary bundle model will be used. The PSD is obtained from the measurement of the total flow rate of fluid as a function of the imposed pressure gradient magnitude. The mathematical processing of the experimental data, which depends on the type of yield stress fluid, provides an overview of the pore size distribution of the porous material. The technique proposed here was successfully tested analytically and numerically for usual pore size distributions such as the Gaussian mono and multimodal distributions. The study was conducted for yield stress fluids obeying the classical Bingham model and extended to the more realistic Herschel-Bulkley model. Unlike other complex methods, expensive and sometimes toxic, this technique presents a lower cost, requires simple measurements and is easy to interpret. This new method could become in the future an alternative, non-toxic and cheap method for the characterization of porous materials.


Introduction
Porous media are found almost everywhere around us [1][2][3][4], whether in living matter (human skin, cartilage, bone, etc.), inert matter (soils, layers sedimentary rocks, etc.) or industrial materials (concretes, cements, powders, textiles, etc.).The characterization of porous media in terms of porosity, specific surface, PSD etc., is an important issue for many industrial sectors: assisted petroleum recovery, thermal insulation of buildings, CO 2 sequestration, energy storage, etc.The phenomena related to flow through porous media have occupied and continue to stimulate a strong research activity, both fundamental and applied.The main objective of this work is to propose a new non-polluting, simple and economic method to determine the pore size distribution of a porous medium without the defects of the currently existing techniques.Among them, we can cite classical methods such as: (i) the mercury intrusion porosimetry (MIP) [1][2][3]5].This technique is based on the existence of a threshold below which the pores cannot be invaded.Indeed, due to its large surface tension mercury does not wet most materials.A * e-mail: a.oukhlef.est@usms.mapressure difference ∆P lg must be imposed so that mercury penetrates the pores whose radii r p are greater than r * p = 2σ lg cos θ/∆P lg where σ lg is the liquid/gas interfacial tension and θ the contact angle.Because of the harmful effects of mercury vapor, this technique is destined to be abandoned; (ii) the method of isothermal gas adsorption [6], based on the molecular Van der Waals interactions between a condensing vapor and the internal surface of the pores is particularly adapted to the characterization of porous media with mesopores (in the range 5 nm r 100 nm).There are other alternative techniques that can characterize the PSD like thermoporometry [7], the Small Angle Neutron (SANS) or X-Ray Scattering (SAXS) [8,9], nuclear magnetic resonance (NMR) [10], not to mention the destructive methods such as stereological analysis of polished slices [11].All these techniques are either difficult to implement or expensive (microscopy, image acquisition by X-ray Computerized Micro-Tomography, etc.).The above cited classical methods are based on the existence of a threshold below which the pores cannot be invaded according to a thermodynamic control parameter.Like the aforementioned techniques, our new approach is based on the injection of non-Newtonian yield stress fluid into the porous medium (invasive method).In the present method, we propose to use the rheological properties of yield stress fluid to scan different pore scales and determine the PSD.The basic idea is the following: in order to set such fluids into motion, it is necessary to impose between both ends of a pore a pressure gradient whose magnitude (∇P ) is greater than a critical value depending on the fluid yield stress (τ 0 ) and the pore radius (r p ).In other words, for a given pressure gradient magnitude ∇P , only the pores whose radius is greater than the critical radius r 0 = 2τ 0 /∇P are invaded.Then it is possible to scan the PSD by adjusting the pressure gradient step by step and measuring the corresponding flow rate Q.From the curve Q = f (∇P ), it is possible to extract the PSD solving an inverse problem.This technique has been successfully tested on unimodal, bi or tri-modal Gaussian PSD.The problem inversion depending on the rheological behavior of the fluid, we started using ideal yield stress fluids such as Bingham fluids [12].Since it is rather difficult to find fluids obeying these rheological laws, we have generalized our method to the Herschel-Bulkley model [13] which better describes most of the non-Newtonian yield stress fluids.In this work, we will give the results of this generalized model by presenting first the analytical expression of the solution of the inverse problem coming from the resolution of the associated integral equation described below.So, we can calculate the probability density depending on the partial derivatives (integer or fractional) of the characteristic curve (Q = f (∇P )).This paper also reports an improvement to this technique: the solution of the problem for the PSD identification using a method based on a numerical approach with yield stress fluids such as Herschel-Bulkley model.

Formulation of the problem and solution
In this paper, the porous material is described by a bundle of capillary tubes [14,15], the radii of which are distributed according to an unknown probability density function p(r).This model, originally introduced by Purcell [16], is used nearly in all theoretical studies of porous media because of its amenability to analytical treatment and because variables that are difficult or impossible to measure physically can be computed directly from this model.Still, let us review the main limitations of this model: -The first one concerns the maximum porosity: for a porous medium whose pores have the same radius, the maximum 2D-packing is the hexagonal packing arrangement with a maximum porosity of φ = π √ 3/6 ≈ 0.907 in an infinite medium (this value decreases when the confinement increases and strongly depends on the shape of the boundaries).The maximum porosity for a random close packing of disks is even lesser with a maximal value φ = 0.84 for an infinite medium.These analytical values overestimate the maximal porosity because the contact of the pores would not be possible in a "real" capillary bundle porous medium in which some solid material has to remain between the pores in order to guarantee the structural robustness of the material.Some values around φ = 0.78 are considered acceptable in practice.Moreover, let us note that this model gives access only to the open porosity since it is based on the flow properties; -It neglects complex features of porous media such as tortuosity, the dead arms and it neglects interconnectivity between pores as it assumes unidirectional flow; -The contraction-expansion features of non-Newtonian yield stress fluid flows cannot be accounted for because this model assumes a single radius along each pore of the bundle with no variation in the cross-sectional area; -It cannot be representative for flows in an anisotropic medium due to its assumption of unique permeability along propagation direction; -It neglects the elongation effects associated with porous media flows.
However, some of the limitations cited above can be erased by more complex variants not considered here.
In our study, we ignore the tortuosity that only affects the characteristic curve by increasing the length of the pores [1,2] and not the non-dimensional PSD.We will use this model to derive the inversion technique which allows us to obtain the PSD from the characteristic curve for a given yield stress fluid in noninertial regimes.Under these conditions, the expression of the total flow rate Q through this capillary model, when a pressure gradient magnitude ∇P is imposed on this system, is a function of the elementary flow rate q(∇P, r, n) in each capillary and of the probability density p(r): In this study, our goal is to extract the probability density p(r) from the curve describing the evolution of total flow rate Q(∇P ) as a function of the pressure gradient magnitude of a Herschel-Bulkley fluid (H-B) through a porous medium.This measured flow rate is calculated directly from the integral equation (1) when p(r) is known.If this distribution is unknown, this relation constitutes a Volterra equation of the first kind, whose kernel q(∇P, r, n) is the elementary flow rate in each capillary and Q(∇P ) constitutes the source term obtained experimentally.The kernel q(∇P, r, n) depends on the yield stress fluid type and the Herschel-Bulkley model is used here because it describes the rheological behavior of most non-Newtonian yield stress fluids more realistically than the Bingham model.Indeed, the H-B model can take into account the pseudo-plastic or dilatant behavior of the real yield stress fluids and constitute a generalization of the Bingham model.In H-B fluids beyond the yield stress τ 0 , the linearity of the shear stress with the shear rate is replaced by a power law.It is, therefore, a model with three parameters: the yield stress τ 0 , the consistency of the fluid k and power-law index n.The rheological behavior law describing this model is given by: (2) with τ the shear stress tensor and D the rate of deformation tensor.For the flow in a tube with circular cross-section, these equations take the simpler form: where τ r z is the shear stress, (∂u z /∂r ) the rate of deformation and r is the radial coordinate.The elementary flow rate of the Herschel-Bulkley fluid through a capillary of circular cross-section with radius r, which constitutes the kernel of equation ( 1), is given by [17]: For r ≤ r 0 the fluid does not flow.The critical radius r 0 = 2τ 0 /∇P is also the radius of the core zone of the plug flow.The pore size distribution p (r) can be obtained through differential operators applied to Q (∇P ) [13]: It should be noted (i) that equation ( 5) is valid only for n ∈ Q but since the rational numbers are a dense subset of the real numbers, any real number has rational numbers arbitrarily close to it and (ii) that for n = 1 the Herschel-Bulkley model reduces to that of Bingham fluid and equation ( 5) reduces to that obtained by Ambari et al. [18] and (iii) that in equation ( 5), (1/n)! is finite even for any real value of n = 0 and is calculated by the Gamma function:

Analytical inverse method to determine the PSD
In the case where the power-law index is n = 1/q with q ∈ N * , a general relationship is obtained between the probability density p(r) and the integer partial derivatives of the total flow ∂ i Q/∂ (∇P ) i .In the case where q is noninteger it is possible to generalize this formula considering these derivatives as fractional.This will be justified later.In order to reduce the number of physical parameters involved in this problem, we will normalize the previous equations, using as characteristic scales L for the length and τ 0 /L for the pressure gradient magnitude where L is the thickness of the studied sample (directly accessible) and not the mean radius of the capillaries which is still unknown and would require a complementary measurement.Therefore, the non-dimensional kernel is then given by: with He H−B = ρτ (2/n) the Hedström number in the case of a Herschel Bulkley fluid [19], ρ the density of fluid and r + 0 = 2/∇P + the non-dimensional critical radius.In equation (6), the flow rate is normalized by the characteristic flow rate: Figure 1 below presents the evolution of the elementary flow rate q + (that is to say the kernel of Eq. ( 1)) in a pore of radius r + = 1 as a function of the imposed pressure gradient magnitude ∇P + for different power-law index and for a Hedström number He H−B = 0.02.
As for non-dimensional pore size distribution: In the case, where n is non-integer, the calculation of the fractional derivatives in equation ( 8) are very difficult unless the derived function is simple.One of the simplest solutions for non-integer derivation is when the derived function is a polynomial.it can be carried out by means of the following relation defined by Riemann [20]: where ν is the non-integer derivation order, Γ the gamma function and m the order of the derivative of a monomial, with (x = ∇P ).In our study, the above relation will be used to calculate the different fractional derivatives that occur in equation ( 8) provided a polynomial fit of the evolution Q + = f (∇P + ) is given (in the present study, all the characteristic curves Q + = f (∇P + ) were fitted by polynomial functions of maximum order 40 because this order leads to a satisfactory inversion of all the curves whatever the value of the power-law index n.However, we have noticed that the lower the value of n, the lower the order of the polynomials to use.For instance, for n = 0.9, a polynomial of order 15 starts to give satisfactory results).Remember that in the case where n = 1/q with q ∈ N * , the two successive derivatives are integers.In all cases, the distribution will be calculated by using equation ( 8).In the next section we explain the approach followed to obtain this inverse distribution.

Inversion in the case of a Bingham fluid (n = 1)
As a first attempt to test this method, let us assume that the fluid is a Bingham fluid (n = 1) and that the PSD is of Gaussian type of mean radius µ and standard deviation σ.
Written in non-dimensional form (r + = r/L, µ + = µ/L, σ + = σ/L and p + (r + ) = p (r) L), this distribution is: and the normalized Volterra equation (Eq.( 1)) is: Using equations ( 6), ( 10) and ( 11), it is possible to calculate the characteristic curve Q + (∇P + ).In Figure 2 above, it is shown for µ + = 10 −4 , σ + = 2 × 10 −5 and He H−B = 0.02.This figure shows the non-dimensional total flow rate vs. the non-dimensional pressure gradient magnitude resulting from the flow of the Bingham fluid through such a porous medium.This figure is characterized by a first region at low pressure gradients in which the flow rate is zero.This region extends until the largest pore in the material is invaded.It is followed by a second region in which the flow rate increases with the pressure gradient.Now, let us start from this characteristic curve Q + (∇P + ) and let us try to calculate p + (r + ).For n = 1, equation ( 8) becomes: Because equation ( 12) has only integer derivatives, it can be directly applied to the non-dimensional characteristic curve (Fig. 2) to calculate p + (r + ).It is plotted in Figure 3 together with the original normal distribution (Eq.( 10)) and the agreement between both curves is very good (see the article [12]).
We have also tried to use the method involving the polynomial fit of the characteristic curve Q + (∇P + ) and its fractional derivatives (Eq.( 9)).The results obtained by this method are shown in Figure 4.It shows that this technique is also relevant even though it introduces some noise for the small radii related to the polynomial fit needed here.

Inversion in the case of a Herschel-Bulkley fluid
for n = 1/q with q ∈ ∈ ∈ N * * * As a second attempt to test this method, let us make the problem more complex and consider a Herschel-Bulkley pseudo-plastic fluid whose power law index has the form  n = 1/q with q ∈ N * .In this case, the inversion procedure also involves integer partial derivatives of order 1 + q and 2 + q.For instance for n = 1/2, the normalized PSD (Eq.( 8)) writes as: Following the same procedure as in the previous section and for the same initial Gaussian distribution, we have calculated the characteristic curve Q + (∇P + ) presented in Figure 5. Now, when equation ( 13) is applied to the characteristic curve (Fig. 5), the initially injected Gaussian PSD is retrieved as we can see in Figure 6.Likewise this distribution can be obtained using the fractional derivatives method (Fig. 7) with some additional noise for the smallest radii because of the poor description of the characteristic curve by the polynomial function when r is small.

in the case of a Herschel-Bulkley fluid
for n = p/q with p ∈ ∈ ∈ N * * * \ {1} {1} {1} and with q ∈ ∈ ∈ N * * * In the case where the flow index of the fluid is a rational number whose numerator is different from unity, the identification formula (Eq.( 8)) presents also non-integer derivatives.In this case, only the fractional derivatives method is usable based on the polynomial fit of the characteristic curve.The results for a pseudo-plastic fluid (n = 0.9) are presented in Figures 8 and 9 for a normal distribution whose parameters are again µ + = 10 −4 , σ + = 2 × 10 −5 and for He H−B = 0.02.Likewise for the same set of parameters, the results concerning a dilatant fluid (n = 1.2) are given in Figures 10 and 11.We note that the inverse distribution can be obtained by the method using the fractional derivatives but like in the case of a Bingham fluid, some fluctuations appear for the  small radii due to the lack of precision of the polynomial fit used here.Therefore, we will present in the next section another numerical method effective and robust and that can be applied to all yield stress fluids such as the Herschel-Bulkley fluids and for all kinds of distributions.Note that this approach, based on the fractional derivatives and the polynomial interpolation function, fails for more complex distributions (bimodal or tri-modal) probably because of the unsuited interpolation function used here.This has motivated another numerical method based on the numerical discretization of the Volterra equation.

Numerical inverse method for determining PSD
An alternative method to solve the Volterra equation (Eq.( 1)) is to use a matrix system.For this, we adopt the following approach: the pressure gradient magnitude (j-index) and the pore radius (i-index) are discretized using N points so that the source term Q (∇P ), the kernel q (∇P, r, n) and the PSD p(r) become respectively   Q j = Q j (∇P j ), q ij = q(∇P j , r i , n) and p i = p(r i ).The Volterra equation is then rewritten as follows: with ∆r = (r max − r min )/(N − 1).Mathematically, if the kernel q ij is not singular, equation ( 14) has a unique solution which is given by the following relation: where (q ij ) −1 is the inverse matrix of the kernel.

Numerical approach
The non-dimensional form of the functions is used to solve the previous system.For this, we use the kernel in the form written in equation ( 6).This kernel is discretized using a non-dimensional radius step ∆r + and a non-dimensional pressure gradient magnitude step ∆(∇P ) + such that: with ∇P + min = 2/r + max and ∇P + max = 2/r + min for given values of r + min and r + max .The discretized non-dimensional radius and pressure gradient magnitude now become: The kernel matrix is then built with the previously discretized elements: At this point, after having experimentally obtained the N components of the source vector Q j and calculated the inverse of the kernel matrix (q ij ) −1 , we can determine the unknown distribution p i .The number of values used (here N = 1000) can naturally be reduced depending on the desired accuracy.Indeed, to numerically invert the Volterra equation, we no longer need a large number of values as previously required in the case of the analytical and fractional derivatives methods.Note that in an experimental work, this is the number of points of the gradient (or flow rate) that would command the number of points for the radius.But within the range of pressure gradients experimentally explored, if one wishes a finer discretization for the radius, it is always possible to interpolate the missing data of the pressure gradient (or the flow rate) in order to find a square matrix.

Numerical inversions of some distributions by using Herschel-Bulkley model
• Herschel-Bulkley fluid with n = 1 (Bingham model) As in most experiments it is cumbersome to obtain a large sample of data, we have tested the accuracy of the present numerical method on a sample of N = 25 points (instead of 1000 points used in the two previous methods).
Figure 12 shows the characteristic curve for a normal distribution (µ + = 10 −4 , σ + = 2 × 10 −5 ) and a Bingham fluid (n = 1) for He H−B = 0.02 and Figure 13 presents the results obtained by the numerical method compared to the initial Gaussian PSD.A good agreement is found here although the size of the sample is small which seems to be encouraging.
To verify the efficiency of this numerical technique with more complex distributions, a bimodal Gaussian distribution is considered, with two peaks at µ + 1 = µ 1 /L = 10 −4 and µ + 2 = 2µ + 1 and two different standard deviations σ + 1 = σ 1 /L = 2 × 10 −5 and σ + 2 = 2σ + 1 .For a Bingham fluid (n = 1) at He H−B = 0.02, the characteristic curve can be evaluated (Fig. 14).For this example, N = 1000 points are used to increase the accuracy.Once again, when the numerical method is applied to the data in Figure 14, the initial bimodal Gaussian distribution is perfectly found (Fig. 15).

Fig. 1 .
Fig. 1.Elementary flow rate vs. pressure gradient magnitude for different power-law index and He H−B = 0.02.

Fig. 2 .
Fig. 2. Total flow rate vs. pressure gradient magnitude for a Gaussian distribution and a Bingham fluid for He H−B = 0.02.

Fig. 3 .
Fig. 3. Comparison between the initial and the analytically calculated PSD for a Bingham fluid.

Fig. 4 .
Fig. 4. Comparison between the initial and the fractional method calculated PSD for a Bingham fluid.

Fig. 7 .
Fig. 7. Comparison between the initial and the fractional method calculated PSD for a Herschel-Bulkley fluid with n = 0.5.

Fig. 8 .
Fig. 8.Total flow rate vs. pressure gradient magnitude for a Gaussian distribution and a Herschel-Bulkley fluid for n = 0.9.

Fig. 9 .
Fig. 9.Comparison between the initial and the fractional method calculated PSD for a Herschel-Bulkley fluid with n = 0.9.

Fig. 10 .
Fig. 10.Total flow rate vs. pressure gradient magnitude for a Gaussian distribution and a Herschel-Bulkley fluid for n = 1.2.

Fig. 11 .
Fig. 11.Comparison between the initial and the fractional method calculated PSD for a Herschel-Bulkley fluid with n = 1.2.

Fig. 12 .
Fig. 12.Total flow rate vs. pressure gradient magnitude for a Gaussian distribution and a Bingham fluid (N = 25).