© A. Oukhlef et al., published by EDP Sciences 2020

1 Introduction

Porous media are found almost everywhere around us [1–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₂ 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–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 pressure difference ΔP_lg must be imposed so that mercury penetrates the pores whose radii r_p are greater than $r_p^{*}=2{\sigma }_{lg}\cos \theta /\Delta 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 $\left(\nabla P\right)$ is greater than a critical value depending on the fluid yield stress $\left({\tau }_0\right)$ and the pore radius $\left(r_p\right)$. In other words, for a given pressure gradient magnitude ∇P, only the pores whose radius is greater than the critical radius r₀ = 2τ₀∕∇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.

2 Models and procedure

2.1 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 $\varphi =\pi \sqrt{3}/6\approx 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 non-inertial 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): $$Q\left(\nabla P\right)={\displaystyle \underset{r_0=\frac{2{\tau }_0}{\nabla P}}{\overset{\infty }{\int }}q}\left(\nabla P,r,n\right)p\left(r\right)dr$$(1)

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 τ₀, 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 τ₀, the consistency of the fluid k and power-law index n. The rheological behavior law describing this model is given by: $$\text{}\text{}\{\text{}\text{}\begin{array}{l}\overline{\overline{\tau }}=2\left(\text{}k\text{}{\left(\text{}\text{}\sqrt{\text{}2\overline{\overline{D}}:\overline{\overline{D}}}\right)}^{n-1}\text{}\text{}+\text{}\frac{{\tau }_0}{\sqrt{2\overline{\overline{D}}:\overline{\overline{D}}}}\text{}\right)\text{}\overline{\overline{D}}\text{for}\sqrt{\frac{\overline{\overline{\tau }}:\overline{\overline{\tau }}}{2}}>{\tau }_0\hfill \\ \overline{\overline{D}}=0\text{for}\sqrt{\frac{\overline{\overline{\tau }}:\overline{\overline{\tau }}}{2}}\le {\tau }_0\hfill \end{array}$$(2)

with $\overline{\overline{\tau }}$ the shear stress tensor and $\overline{\overline{D}}$ the rate of deformation tensor. For the flow in a tube with circular cross-section, these equations take the simpler form: $$\{\begin{array}{l}{\tau }_{r^{\prime }z}={\tau }_0+k{\left|\frac{\partial u_z}{\partial r^{\prime }}\right|}^n\text{for}{\tau }_{r^{\prime }z}>{\tau }_0\hfill \\ \frac{\partial u_z}{\partial r^{\prime }}=0\text{for}{\tau }_{r^{\prime }z}\le {\tau }_0\hfill \end{array}$$(3)

where ${\tau }_{r^{\prime }z}$ is the shear stress, $\left(\partial u_z/\partial r^{\prime }\right)$ 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]: $$\begin{array}{l}q\left(\nabla P,r,n\right)\\ =\text{}\{\begin{array}{l}n{\pi }_{}{r}^3{\left(\frac{r\nabla P}{2k}\right)}^{\frac{1}{n}}{\left(1-\frac{r_0}{r}\right)}^{\frac{n+1}{n}}\hfill \\ \left[\frac{{\left(1-\frac{r_0}{r}\right)}^2}{\left(3n+1\right)}+\frac{\frac{2r_0}{r}\left(1-\frac{r_0}{r}\right)}{\left(2n+1\right)}+\frac{{\left(\frac{r_0}{r}\right)}^2}{\left(n+1\right)}\right]\text{for}r>r_0\hfill \\ \hfill \\ 0\text{for}r\le r_0\hfill \end{array}\end{array}$$(4)

For r ≤ r₀ the fluid does not flow. The critical radius r₀ = 2τ₀∕∇P is also the radius of the core zone of the plug flow. The pore size distribution $p\left(r\right)$ can be obtained through differential operators applied to $Q\left(\nabla P\right)$ [13]: $$\begin{array}{ll}p\left(r\right)=\hfill & \frac{2^{\left(1+3n\right)/n}{k}^{1/n}{\left(\nabla P\right)}^2}{16\pi {\tau }_0{r}^{\left(1+3n\right)/n}\left(1/n\right)!}[\left(\frac{1+4n}{n}\right)\frac{{\partial }^{\left(1+n\right)/n}.}{\partial {\left(\nabla P\right)}^{\left(1+n\right)/n}}\hfill \\ \hfill & +\nabla P\frac{{\partial }^{\left(1+2n\right)/n}.}{\partial {\left(\nabla P\right)}^{\left(1+2n\right)/n}}]Q {|}_{\nabla P=\frac{2{\tau }_0}{r}}\hfill \end{array}$$(5)

It should be noted (i) that equation (5) is valid only for $n\in \mathbb{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: $$\left(\frac{1}{n}\right)!=\Gamma \left(1+\frac{1}{n}\right)={\displaystyle \underset{0}{\overset{\infty }{\int }}{e}^{-x}}{x}^{1/n}dx$$

2.2 Analytical inverse method to determine the PSD

In the case where the power-law index is n = 1∕q with $q\in {\mathbb{N}}^{*}$, a general relationship is obtained between the probability density p(r) and the integer partial derivatives of the total flow ${\partial }^{i}Q/\partial {\left(\nabla P\right)}^i$. In the case where q is non-integer 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 τ₀ ∕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: $$\begin{array}{l}{q}^{+}\left(\nabla P^{+},r^{+},n\right)\\ =\text{}\{\begin{array}{l}H{e}^{H-B}\frac{1}{2^{\frac{1}{n}}}\left(\frac{1}{2n+1}\right)\frac{1}{{\left(\nabla P^{+}\right)}^3}{\left(r^{+}\nabla P^{+}-2\right)}^{1+\frac{1}{n}}\hfill \\ [8n^2+4n\left(n+1\right)r^{+}\nabla P^{+}\hfill \\ +\left(2n^2+3n+1\right){\left(r^{+}\nabla P^{+}\right)}^2]\text{for}{r}^{+}>r_0^{+}\hfill \\ \hfill \\ 0\text{for}{r}^{+}\le r_0^{+}\hfill \end{array}\end{array}$$(6)

with $H{e}^{H-B}=\rho {\tau }_0^{\left(2-n\right)/n}{L}^2/k^{\left(2/n\right)}$ the Hedström number in the case of a Herschel Bulkley fluid [19], ρ the density of fluid and $r_0^{+}=2/\nabla P^{+}$ the non-dimensional critical radius. In equation (6), the flow rate is normalized by the characteristic flow rate: $$q_c=\frac{n\pi }{\left(3n+1\right)\left(n+1\right)}\left(\frac{L{k}^{1/n}{\tau }_0^{\left(n-1\right)/n}}{\rho }\right)$$(7)

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: $$\begin{array}{ll}{p}^{+}\left(r^{+}\right)\hfill & =\frac{n{\left(\nabla P^{+}\right)}^{\left(1+5n\right)/n}}{16H{e}^{H-B}\left(3n+1\right)\left(1+n\right)\left(1/n\right)!}\hfill \\ \hfill & [\left(\frac{1+4n}{n}\right)\frac{{\partial }^{\left(1+n\right)/n}.}{\partial {\left(\nabla P^{+}\right)}^{\left(1+n\right)/n}}\hfill \\ \hfill & +\nabla P^{+}\frac{{\partial }^{\left(1+2n\right)/n}.}{\partial {\left(\nabla P^{+}\right)}^{\left(1+2n\right)/n}}]Q^{+} {|}_{\nabla P^{+}=\frac{2}{r^{+}}}\hfill \end{array}$$(8)

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]: $$\frac{d^{\nu }\left(x^m\right)}{d{x}^{\nu }}=\frac{\Gamma \left(m+1\right)}{\Gamma \left(m-\nu +1\right)}{x}^{m-\nu }$$(9)

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\left(\nabla P^{+}\right)$ is given (in the present study, all the characteristic curves $Q^{+}=f\left(\nabla P^{+}\right)$ 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\in {\mathbb{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.

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

2.2.1 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^{+}\left(r^{+}\right)=p\left(r\right)L$), this distribution is: $$p^{+}\left(r^{+}\right)=\frac{1}{{\sigma }^{+}\sqrt{2\pi }}\exp \left[-\frac{{\left(r^{+}-{\mu }^{+}\right)}^2}{2{\sigma }^{+}{}^2}\right]$$(10)

and the normalized Volterra equation (Eq. (1)) is: $$Q^{+}\left(\nabla P^{+}\right)={\displaystyle \underset{r_0^{+}}{\overset{\infty }{\int }}{q}^{+}}\left(\nabla P^{+},r^{+},n=1\right)p^{+}\left(r^{+}\right)d{r}^{+}$$(11)

Using equations (6), (10) and (11), it is possible to calculate the characteristic curve $Q^{+}\left(\nabla P^{+}\right)$. In Figure 2 above, it is shown for μ⁺ = 10⁻⁴, σ⁺ = 2 × 10⁻⁵ 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^{+}\left(\nabla P^{+}\right)$ and let us try to calculate $p^{+}\left(r^{+}\right)$. For n = 1, equation (8) becomes: $$\begin{array}{ll}{p}^{+}\left(r^{+}\right)=\hfill & \frac{{\left(\nabla P^{+}\right)}^6}{128H{e}^{H-B}}\left[5\frac{{\partial }^2.}{\partial {\left(\nabla P^{+}\right)}^2}+\nabla P^{+}\frac{{\partial }^3.}{\partial {\left(\nabla P^{+}\right)}^3}\right]\hfill \\ \hfill & \times Q^{+} {|}_{\nabla P^{+}=\frac{2}{r^{+}}}\hfill \end{array}$$(12)

Because equation (12) has only integer derivatives, it can be directly applied to the non-dimensional characteristic curve (Fig. 2) to calculate $p^{+}\left(r^{+}\right)$. 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^{+}\left(\nabla P^{+}\right)$ 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.

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

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

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

2.2.2 Inversion in the case of a Herschel-Bulkley fluid for n = 1/q with $q\in {\mathbb{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\in {\mathbb{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: $$\begin{array}{ll}{p}^{+}\left(r^{+}\right)=\hfill & \frac{{\left(\nabla P^{+}\right)}^7}{240H{e}^{H-B}}\left[6\frac{{\partial }^3.}{\partial {\left(\nabla P^{+}\right)}^3}+\nabla P^{+}\frac{{\partial }^4.}{\partial {\left(\nabla P^{+}\right)}^4}\right]\hfill \\ \hfill & \times Q^{+} {|}_{\nabla P^{+}=\frac{2}{r^{+}}}\hfill \end{array}$$(13)

Following the same procedure as in the previous section and for the same initial Gaussian distribution, we have calculated the characteristic curve $Q^{+}\left(\nabla P^{+}\right)$ 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.

Fig. 5 Total flow rate vs. pressure gradient magnitude for a Gaussian distribution and a Herschel-Bulkley fluid for n = 1∕2 and HeH−B = 0.02.

Fig. 6 Comparison between the initial and the analytically calculated PSD for a Herschel-Bulkley fluid with n = 0.5.

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

2.2.3 Inversion in the case of a Herschel-Bulkley fluid for n = p∕q with $p\in {\mathbb{N}}^{*}$\{1} and with $q\in {\mathbb{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 inFigures 8 and 9 for a normal distribution whose parameters are again μ⁺ = 10⁻⁴, σ⁺ = 2 × 10⁻⁵ 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.

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

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

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

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

2.3 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\left(\nabla P\right)$, the kernel $q\left(\nabla P,r,n\right)$ 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: $$Q_j={\sum }_{i=1}^N{q}_{ij}{p}_i\Delta rforj=1,\mathrm{...},N$$(14)

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: $$p_i=\frac{{\left(q_{ij}\right)}^{-1}{Q}_j}{\Delta r}fori=1,\mathrm{...},N$$(15)

where ${\left(q_{ij}\right)}^{-1}$ is the inverse matrix of the kernel.

2.3.1 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 $\Delta {\left(\nabla P\right)}^{+}$ such that: $$\Delta r^{+}=(r_{\text{max}}^{+}-r_{\text{min}}^{+})/(N-1)$$ $$\Delta {\left(\nabla P\right)}^{+}=(\nabla P_{\text{max}}^{+}-\nabla P_{\text{min}}^{+})/(N-1)$$

with $\nabla P_{\text{min}}^{+}=2/r_{\text{max}}^{+}$ and $\nabla P_{\text{max}}^{+}=2/r_{\text{min}}^{+}$ for given values of $r_{\text{min}}^{+}$ and $r_{\text{max}}^{+}$. The discretized non-dimensional radius and pressure gradient magnitude now become: $$r_i^{+}=r_{\text{min}}^{+}+\left(i-1\right)\Delta r^{+}fori=1,\mathrm{...},N$$ $$\nabla P_j^{+}=\nabla P_{\text{min}}^{+}+\left(j-1\right)\Delta {\left(\nabla P\right)}^{+}forj=1,\mathrm{...},N$$

The kernel matrix is then built with the previously discretized elements: $$q_{ij}^{+}=q^{+}\left(\nabla P_j^{+},r_i^{+},n\right)$$

with $$\begin{array}{ll}{q}^{+}\left(\nabla P_j^{+},r_i^{+},n\right)=\hfill & H{e}^{H-B}\frac{1}{2^{\frac{1}{n}}}\left(\frac{1}{2n+1}\right)\frac{1}{{\left(\nabla P_j^{+}\right)}^3}\hfill \\ \hfill & \times {\left(r_i^{+}\nabla P_j^{+}-2\right)}^{1+\frac{1}{n}}\hfill \\ \hfill & \times [\left(2n^2+3n+1\right){\left(r_i^{+}\nabla P_j^{+}\right)}^{2}8n^2+4n\left(n+1\right)r_i^{+}\nabla P_j^{+}\hfill \\ \hfill & +\left(2n^2+3n+1\right){\left(r_i^{+}\nabla P_j^{+}\right)}^2]\hfill \end{array}$$

At this point, after having experimentally obtained the N components of the source vector Q_j and calculated the inverse of the kernel matrix ${\left(q_{ij}\right)}^{-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 pressure 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.

2.3.2 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⁻⁴, σ⁺ = 2 × 10⁻⁵) 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 ${\mu }_1^{+}={\mu }_1/L={10}^{-4}$ and ${\mu }_2^{+}=2{\mu }_1^{+}$ and two different standard deviations ${\sigma }_1^{+}={\sigma }_1/L=2\times {10}^{-5}$ and ${\sigma }_2^{+}=2{\sigma }_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).

• Herschel-Bulkley fluid with n = 0.9 (pseudo-plastic) and n = 1.2 (dilatant)

Finally, we can test this method on pseudo-plastic or dilatant yield stress fluids. Let us recall that these two cases generally require the use of the fractional derivatives method leading to some noise when the pore radius becomes very small. The two fluids tested here are characterized by the following parameters: He^H−B = 0.02, n =0.9 (pseudo-plastic) and n =1.2 (dilatant). A tri-modal Gaussian distribution is tested with ${\mu }_1^{+}={10}^{-4}$, ${\mu }_2^{+}=2{\mu }_1^{+}$, ${\mu }_3^{+}=3{\mu }_1^{+}$, ${\sigma }_1^{+}=2\times {10}^{-5}$, ${\sigma }_2^{+}=3{\sigma }_1^{+}/2$ and ${\sigma }_3^{+}=2{\sigma }_1^{+}$.

Figures 16 and 17 validate the pseudo-plastic case (n = 0.9) and Figures 18 and 19 validate the dilatant case (n = 1.2). Therefore, this numerical method seems very effective and robust and can be applied to all kind of yield stress fluids such as the Herschel-Bulkley for any power-law index n < 1 (pseudo-plastic behavior) orn > 1 (dilatant behavior) and for all kinds of distributions. All these tests show that the numerical inversion method is relevant for determining the pore size distribution p(r) provided that a sufficient number of points and well-adapted radius and pressure gradient magnitude intervals are taken.

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

Fig. 13 Comparison between the initial and the numerically calculated PSD for a Bingham model with N = 25.

Fig. 14 Total flow rate vs. pressure gradient magnitude for a bimodal distribution (${\mu }_1^{+}={10}^{-4}$, ${\mu }_2^{+}=2{\mu }_1^{+}$, ${\sigma }_1^{+}=2\times {10}^{-5}$, ${\sigma }_2^{+}=2{\sigma }_1^{+}$) and a Bingham fluid (n = 1) for N = 1000 points.

Fig. 15 Comparison between the initial and the numerically calculated PSD for a bimodal distribution (${\mu }_1^{+}={10}^{-4}$, ${\mu }_2^{+}=2{\mu }_1^{+}$, ${\sigma }_1^{+}=2\times {10}^{-5}$, ${\sigma }_2^{+}=2{\sigma }_1^{+}$) and a Bingham fluid (n = 1) for N = 1000 points.

Fig. 16 Total flow rate vs. pressure gradient magnitude for a tri-modal distribution with (${\mu }_1^{+}={10}^{-4}$, ${\mu }_2^{+}=2{\mu }_1^{+}$, ${\mu }_3^{+}=3{\mu }_1^{+}$, ${\sigma }_1^{+}=2\times {10}^{-5}$, ${\sigma }_2^{+}=3{\sigma }_1^{+}/2$ and ${\sigma }_3^{+}=2{\sigma }_1^{+}$) and a Herschel-Bulkley pseudo-plastic fluid (n = 0.9).

Fig. 17 Comparison between the initial and the numerically calculated PSD for a tri-modal distribution (${\mu }_1^{+}={10}^{-4}$, ${\mu }_2^{+}=2{\mu }_1^{+}$, ${\mu }_3^{+}=3{\mu }_1^{+}$, ${\sigma }_1^{+}=2\times {10}^{-5}$, ${\sigma }_2^{+}=3{\sigma }_1^{+}/2$ and ${\sigma }_3^{+}=2{\sigma }_1^{+}$) and a Herschel-Bulkley pseudo-plastic fluid (n = 0.9).

Fig. 18 Total flow rate vs. pressure gradient magnitude for a tri-modal distribution (${\mu }_1^{+}={10}^{-4}$, ${\mu }_2^{+}=2{\mu }_1^{+}$, ${\mu }_3^{+}=3{\mu }_1^{+}$, ${\sigma }_1^{+}=2\times {10}^{-5}$, ${\sigma }_2^{+}=3{\sigma }_1^{+}/2$ and ${\sigma }_3^{+}=2{\sigma }_1^{+}$) and a Herschel-Bulkley dilatant fluid (n = 1.2).

Fig. 19 Comparison between the initial and the numerically calculated PSD for a tri-modal distribution (${\mu }_1^{+}={10}^{-4}$, ${\mu }_2^{+}=2{\mu }_1^{+}$, ${\mu }_3^{+}=3{\mu }_1^{+}$, ${\sigma }_1^{+}=2\times {10}^{-5}$, ${\sigma }_2^{+}=3{\sigma }_1^{+}/2$ and ${\sigma }_3^{+}=2{\sigma }_1^{+}$) and a Herschel-Bulkley dilatant fluid (n = 1.2).

3 Conclusion

This work presents a new method for determining and identifying the pore size distribution of porous materials. It is based on the capillary bundle model like most of the other experimental techniques. This method uses non-Newtonian yield stress fluids. The existence of a yield stress gives the possibility to scan the pore distribution. Analytically, the total flow rate of such fluids into a porous medium appears as a Volterra integral equation of the first kind whose kernel is analytically known in the case of the flow into a straight capillary. The mathematical determination of the probability density p(r) function is possible using the partial derivatives (of integer or non-integer order) of the total flow rate of fluid through the porous medium as a function of the pressure gradient magnitude or by a numerical method based on the inverse matrix of the kernel. These techniques are successfully tested for Herschel-Bulkley fluids in the case of mono, bi or tri-modal Gaussian distributions but any other distribution or yield stress fluid (Casson, Robertson-Stiff, etc.) could be used potentially. However, in some cases the resolution of the Volterra equation should be done numerically and not analytically. We hope that this method could become in the future an alternative to the toxic and expensive existing methods for the characterization of porous materials.

References

All Figures

	Fig. 1 Elementary flow rate vs. pressure gradient magnitude for different power-law index and HeH−B = 0.02.
In the text

	Fig. 2 Total flow rate vs. pressure gradient magnitude for a Gaussian distribution and a Bingham fluid for HeH−B = 0.02.
In the text

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

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

	Fig. 5 Total flow rate vs. pressure gradient magnitude for a Gaussian distribution and a Herschel-Bulkley fluid for n = 1∕2 and HeH−B = 0.02.
In the text

	Fig. 6 Comparison between the initial and the analytically calculated PSD for a Herschel-Bulkley fluid with n = 0.5.
In the text

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

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

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

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

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

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

	Fig. 13 Comparison between the initial and the numerically calculated PSD for a Bingham model with N = 25.
In the text

	Fig. 14 Total flow rate vs. pressure gradient magnitude for a bimodal distribution (${\mu }_1^{+}={10}^{-4}$, ${\mu }_2^{+}=2{\mu }_1^{+}$, ${\sigma }_1^{+}=2\times {10}^{-5}$, ${\sigma }_2^{+}=2{\sigma }_1^{+}$) and a Bingham fluid (n = 1) for N = 1000 points.
In the text

	Fig. 15 Comparison between the initial and the numerically calculated PSD for a bimodal distribution (${\mu }_1^{+}={10}^{-4}$, ${\mu }_2^{+}=2{\mu }_1^{+}$, ${\sigma }_1^{+}=2\times {10}^{-5}$, ${\sigma }_2^{+}=2{\sigma }_1^{+}$) and a Bingham fluid (n = 1) for N = 1000 points.
In the text

	Fig. 16 Total flow rate vs. pressure gradient magnitude for a tri-modal distribution with (${\mu }_1^{+}={10}^{-4}$, ${\mu }_2^{+}=2{\mu }_1^{+}$, ${\mu }_3^{+}=3{\mu }_1^{+}$, ${\sigma }_1^{+}=2\times {10}^{-5}$, ${\sigma }_2^{+}=3{\sigma }_1^{+}/2$ and ${\sigma }_3^{+}=2{\sigma }_1^{+}$) and a Herschel-Bulkley pseudo-plastic fluid (n = 0.9).
In the text

	Fig. 17 Comparison between the initial and the numerically calculated PSD for a tri-modal distribution (${\mu }_1^{+}={10}^{-4}$, ${\mu }_2^{+}=2{\mu }_1^{+}$, ${\mu }_3^{+}=3{\mu }_1^{+}$, ${\sigma }_1^{+}=2\times {10}^{-5}$, ${\sigma }_2^{+}=3{\sigma }_1^{+}/2$ and ${\sigma }_3^{+}=2{\sigma }_1^{+}$) and a Herschel-Bulkley pseudo-plastic fluid (n = 0.9).
In the text

	Fig. 18 Total flow rate vs. pressure gradient magnitude for a tri-modal distribution (${\mu }_1^{+}={10}^{-4}$, ${\mu }_2^{+}=2{\mu }_1^{+}$, ${\mu }_3^{+}=3{\mu }_1^{+}$, ${\sigma }_1^{+}=2\times {10}^{-5}$, ${\sigma }_2^{+}=3{\sigma }_1^{+}/2$ and ${\sigma }_3^{+}=2{\sigma }_1^{+}$) and a Herschel-Bulkley dilatant fluid (n = 1.2).
In the text

	Fig. 19 Comparison between the initial and the numerically calculated PSD for a tri-modal distribution (${\mu }_1^{+}={10}^{-4}$, ${\mu }_2^{+}=2{\mu }_1^{+}$, ${\mu }_3^{+}=3{\mu }_1^{+}$, ${\sigma }_1^{+}=2\times {10}^{-5}$, ${\sigma }_2^{+}=3{\sigma }_1^{+}/2$ and ${\sigma }_3^{+}=2{\sigma }_1^{+}$) and a Herschel-Bulkley dilatant fluid (n = 1.2).
In the text