Surface permeability and surface flow tortuosity of particulate porous media

TL;DR

This paper develops a numerical approach to determine surface permeability and flow tortuosity in particulate porous media, enhancing the understanding of dispersion processes at low saturation levels.

Contribution

It extends the Laplace-Beltrami boundary-value problem methodology to interconnected particle ensembles using surface finite element techniques.

Findings

Accurate estimation of surface permeability in particle ensembles.

Quantitative analysis of flow tortuosity effects.

Validation of the numerical approach for complex geometries.

Abstract

The dispersion process in particulate porous media at low saturation levels takes place over the surface elements of constituent particles and, as we have found previously by comparison with experiments, can be accurately described by super-fast non-linear diffusion partial differential equations. To enhance the predictive power of the mathematical model in practical applications, one requires the knowledge of the effective surface permeability of the particle-in-contact ensemble, which can be directly related with the macroscopic permeability of the particulate media. We have shown previously that permeability of a single particulate element can be accurately determined through the solution of the Laplace-Beltrami Dirichlet boundary-value problem. Here, we demonstrate how that methodology can be applied to study permeability of a randomly packed ensemble of interconnected particles. Using surface finite element techniques we examine numerical solutions to the Laplace-Beltrami problem set in the multiply-connected domains of interconnected particles. We are able to rigorously estimate tortuosity effects of the surface flows in a particle ensemble setting.

Tables

Table 1. Table 1: Tabulated values of the system parameters in the random configurations used in the study. Here, N 1 , N 2 subscript 𝑁 1 subscript 𝑁 2 N_{1},N_{2} and N 3 subscript 𝑁 3 N_{3} are the number of particles in the configuration with radii R 1 = 1.3 ​ R 0 , R 2 = R 0 formulae-sequence subscript 𝑅 1 1.3 subscript 𝑅 0 subscript 𝑅 2 subscript 𝑅 0 R_{1}=1.3\,R_{0},R_{2}=R_{0} and R 3 = 0.7 ​ R 0 subscript 𝑅 3 0.7 subscript 𝑅 0 R_{3}=0.7\,R_{0} respectively, ϕ S subscript italic-ϕ 𝑆 \phi_{S} is porosity of the sample and R ¯ = N k ​ R k ∑ j = 1 3 N j ¯ 𝑅 subscript 𝑁 𝑘 subscript 𝑅 𝑘 superscript subscript 𝑗 1 3 subscript 𝑁 𝑗 \bar{R}=\frac{N_{k}R_{k}}{\sum_{j=1}^{3}N_{j}} .

	Parameters of the configurations
Configuration	$N_1$	$N_2$	$N_3$	$\overline{R}/R_0$	${\varphi }_S$ (%)	$L_x^B/R_0$	$L_y^B/R_0$	$L_z^B/R_0$
1	5	6	4	1	52	5.2	4.6	5.5
2	5	7	5	1	47	5.7	4.8	4.9
3	5	4	4	1	47	5.2	3.7	5.3

Equations

$$\frac{\partial s}{\partial t}=\nabla \cdot \left\{D(s)\nabla s\right\},t>0,$$

$$D(s)=\frac{D_0(s)}{(s-s_0{)}^{3/2}},s>s_0,$$

$$\frac{\partial (\varphi s)}{\partial t}+\nabla \cdot \mathbf{Q}=0,$$

$$p=-p_0\frac{A_c}{(s-s_0{)}^{1/2}}$$

$$-\frac{{\kappa }_m}{\mu }\nabla u=\mathbf{q}.$$

$$k_m\propto {\delta }_R^2.$$

$$D_0(s)=A_c\frac{K(s)}{\mu }\frac{p_0}{2\varphi },$$

$$\mathbf{Q}=-\frac{K}{\mu }\nabla p.$$

$$K(s)\propto \frac{1}{|\ln (s-s_0)|}.$$

$$\nabla \cdot \mathbf{q}=0.$$

$${\Delta }_{\Gamma }u=0.$$

$${\delta }_R\ll R,$$

$${\Delta }_{{\Omega }_0}u=0.$$

$${u|}_{\partial {\Gamma }_1}=U_1=const,{u|}_{\partial {\Gamma }_2}=U_2=const,$$

$$Q_T={\delta }_R\frac{{\kappa }_m}{\mu }{\int }_{\partial {\Gamma }_1}\frac{\partial u}{\partial n_s}dl=-{\delta }_R\frac{{\kappa }_m}{\mu }{\int }_{\partial {\Gamma }_2}\frac{\partial u}{\partial n_s}dl,$$

$$Q=\frac{Q_T}{D^2}=-\frac{K_1}{\mu }\frac{U_2-U_1}{D},$$

$$\frac{1}{\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial u}{\partial \theta }\right)=0,{\theta }_0<\theta <\pi -{\theta }_1,$$

$${u|}_{\theta ={\theta }_0}=U_1,{u|}_{\theta =\pi -{\theta }_1}=U_2.$$

$$u={\Psi }_0(U_2-U_1)\ln \left\{\frac{\sin \theta }{\sin {\theta }_0}\frac{1+\cos {\theta }_0}{1+\cos \theta }\right\}+U_1,$$

$${\Psi }_0=\frac{1}{\ln \left\{\frac{\sin {\theta }_1}{\sin {\theta }_0}\frac{1+\cos {\theta }_0}{1-\cos {\theta }_1}\right\}}.$$

$$Q_T=-\frac{K_1}{\mu }D(U_2-U_1)=-2\pi \sin {\theta }_0{\delta }_R\frac{k_m}{\mu }{\frac{\partial u}{\partial \theta }|}_{\theta ={\theta }_0}$$

$$=-(U_2-U_1)2\pi {\delta }_R{\Psi }_0\frac{k_m}{\mu }.$$

$$K_1=\pi {\Psi }_0\frac{{\delta }_R}{R}{k}_m.$$

$$K_1=\frac{{\delta }_R}{2R}\frac{\pi k_m}{|\ln {\theta }_0|}+o\left(\frac{1}{|\ln {\theta }_0|}\right).$$

$${\sin }^2{\theta }_0\approx {\theta }_0^2=\sqrt{s-s_0}.$$

$$K(s)\approx 2\frac{{\delta }_R}{R}\frac{\pi k_m}{|\ln (s-s_0)|}.$$

$$D(s)\propto \frac{1}{|\ln (s-s_0)|(s-s_0{)}^{3/2}}.$$

$$u_1{|}_{\partial {\Gamma }_1^{(2)}}=u_2{|}_{\partial {\Gamma }_2^{(1)}}=const$$

$${\oint }_{\partial {\Gamma }_1^{(2)}}\nabla u_1\cdot {\mathbf{n}}_{s_1}{|}_{\partial {\Gamma }_1^{(2)}}dl=-{\oint }_{\partial {\Gamma }_2^{(1)}}\nabla u_2\cdot {\mathbf{n}}_{s_2}{|}_{\partial {\Gamma }_2^{(1)}}dl$$

$$\frac{1}{\sin \theta }\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial u_1}{\partial \theta }\right)=0,{\theta }_0\le \theta \le \pi -{\theta }_0$$

Full text

I Introduction

Liquid transport in particulate porous media, such as sand, is customarily classified into fully saturated, funicular and pendular regimes of spreading Bear-Book ; Herminghaus-2005 ; Herminghaus-2008 ; Herminghaus-2008-2 . The first two regimes of the liquid dispersion occur at relatively high saturation levels $s>s_{c}\approx 10\%$, where saturation $s$ is defined as the ratio of the liquid volume $V_{L}$ to the volume of available voids $V_{E}$ in a sample volume element $V$, $s=\frac{V_{L}}{V_{E}}$. At high saturation levels, above the critical value $s_{c}$, liquid transport takes place in the pore space either fully or partially filled by the liquid.

Our prime concern here is the special case of liquid dispersion at low saturation levels. As the saturation level drops below the critical value, $s\leq s_{c}$, that is to the value relevant to the pendular regime of spreading, the liquid volumes in the porous matrix become isolated Herminghaus-2005 ; Herminghaus-2008 ; Herminghaus-2008-2 . As a result, at low saturation levels, the liquid is only contained in the pendular rings formed at the locations of the particle contacts and on the particle rough surfaces, and the liquid transport can only occur over the matrix surface elements, as is illustrated in Fig. 1.

Our main concern here is the wetting cycle, when the liquid spreads over a dry porous matrix or over a matrix with a very low background saturation level up to $s_{r}\approx 2\%$. These conditions are similar to those in the case studied previously experimentally and theoretically in Lukyanov2019 . The main driving force of the dispersion process, as is often the case during the wetting cycle, is capillary pressure developed at the moving front in the process of wetting of dry porous matrix, while the liquid bridges play a role of variable liquid reservoirs of uniform surface curvature.

The analysis of this regime of wetting, which is crucial for studies of biological processes and spreading of non-volatile liquids in arid natural environments and industrial installations, has shown that the liquid dispersion has many distinctive features and can be accurately described by the so-called superfast non-linear diffusion equation Lukyanov2012 ; Lukyanov2019 .

Theoretically, the superfast non-linear diffusion equation belongs to a special class of mathematical models. Unlike in the standard porous medium equation Vazquez-Book , in this special case, the non-linear coefficient of diffusion $D(s)$ demonstrates divergent behaviour as a function of saturation $s$, $D(s)\propto(s-s_{0})^{-3/2}$, where $s_{0}$ is some minimal saturation level ($s_{0}\approx 0.5\%$), which could be only achieved in a state when the liquid bridges cease to exist completely Lukyanov2012 ; Herminghaus-2008 ; Herminghaus-2008-2 ; Lukyanov2019 . Note, in that respect, that in the domain of spreading liquid bridges are supposed to never vanish, so that the condition $s>s_{0}$ is always fulfilled in the model, and there is no actual singularity of the mathematical description Lukyanov2012 ; Lukyanov2019 .

Specifically, in the macroscopic approximation, that is after averaging over some volume element containing many particles of the porous medium, the diffusion process in the slow creeping flow conditions can be described by the following non-linear diffusion equation

[TABLE]

where

[TABLE]

for $D_{0}>0$.

The details of derivation of (1) can be found in Lukyanov2012 ; Lukyanov2019 , here we note that, the resultant governing non-linear equation (1) directly follows from the conservation of mass principle

[TABLE]

where $\phi$ is porosity defined as $\displaystyle\phi=\frac{V_{E}}{V}$, which is further assumed to be constant, and $\bf Q$ is the macroscopic flux density. The macroscopic flux density $\bf Q$ is defined in such a way that the total flux through the surface of a macroscopic sample volume element is given by the surface integral $\int{\bf Q\cdot n}\,dS$, where $\bf n$ is the normal vector to the surface of the sample volume element.

To obtain (1) from (2), one needs to apply the capillary pressure-saturation relationship Halsey1998 ; Lukyanov2012 ; Lukyanov2019 dictated by the liquid bridges behaviour

[TABLE]

and the local Darcy’s law Yost-1998 ; Tuller-2000 describing the surface flow in the rough layer of the particle elements

[TABLE]

Here, $A_{c}=\sqrt{\frac{3}{4}\,\frac{1-\phi}{\phi}\frac{N_{c}}{\pi}}$, $N_{c}$ is the coordination number, that is the average number of bridges per a particle, $p_{0}=\frac{2\gamma}{R}\cos\phi_{c}$, $\gamma$ is the coefficient of the surface tension of the liquid, $\phi_{c}$ is the contact angle made by the free surface of the liquid bridge with the rough solid surface of the constituent particles, $R$ is an average radius of the porous medium particles, ${\bf q}$ and $u$ are the averaged local flux density and pressure in the rough surface layer, $\mu$ is liquid viscosity and $k_{m}$ is the local coefficient of permeability of the rough surface, which is proportional to the average amplitude of the surface roughness $\delta_{R}$, that is the width of the surface layer conducting the liquid flux

[TABLE]

One needs to emphasise here that two levels of averaging are involved in obtaining the final governing equation (1). While equations (1), (2) and (3) are ’truly’ macroscopic, that is obtained by averaging using a volume element $V$ containing many grain particles, equation (4) is only an average over some rough area of a single particle containing many surface irregularities, so that quantities ${\bf q}$ and $u$ are also only local averages over that sample surface area.

Therefore, to transit from (4) to the macroscopic description, the spatial averaging theorem formulated in Whitaker-1969 should be applied. That is, using intrinsic liquid averaging $\langle...\rangle^{l}=V_{l}^{-1}\int_{V_{l}}d^{3}x$, where $V_{l}$ is liquid volume within the sample volume $V$, one has $\langle u\rangle^{l}=p$ and $\langle{\bf q}\rangle^{l}\frac{S_{e}}{S}={\bf Q}$. Here, $S$ is the surface area of the sample volume $V$ with the effective area of entrances and exits $S_{e}$. Note, the ratio $S_{e}/S$ is not just a geometric property, but also takes into account the connectivity of the porous elements. For example, the effective area of entrances and exits $S_{e}$ is only defined by the pathways open to the flow.

As a result of the two-level averaging

[TABLE]

where $K(s)=\kappa_{m}\frac{S_{e}}{S}$ is the coefficient of permeability defined by

[TABLE]

The global surface permeability of the particles $K$ as a function of saturation is one of the main elements of the model to accurately represent liquid dispersion at low saturation levels. It is fully defined by the particle geometry and the geometry of the liquid bridge contact areas, Fig. 1 and Fig. 2.

In particular, the disposition and the size of the liquid bridges on the particle surface, that is the size of the domains $\Omega_{1,2}$ and the angle $\alpha$, should play a leading role in defining the resistance to the surface flow. It is not difficult to discern that any variations of the contact area covered by the liquid bridges (pendular rings), that is areas $\Omega_{1,2}$ shown in Fig. 2, or the value of the bridge volume, should affect the global permeability.

Previously, we have shown that permeability of a single particle element can be determined by means of a solution to the equivalent Laplace-Beltrami boundary value problem formulated in the flow domain $\Omega_{0}$ with the boundaries $\partial\Gamma_{1,2}$ in Fig. 2 Penpark2018 . We briefly formulate that problem and summarise the previous results in the next part. Here we note that, based on the analysis of the problem, we have been able to show that in a special azimuthally symmetric case of spherical particles, when the two areas covered by the liquid bridges, domains $\Omega_{1}$ and $\Omega_{2}$ in Fig. 2, are oriented symmetrically to each other, that is at $\alpha=\pi$, the permeability $K$ is supposed to follow the scaling

[TABLE]

We have studied several generalisations of the symmetric problem, such as arbitrary oriented domains, $\alpha\neq\pi$, on the surface of the spherical particle, and particles of arbitrary shapes emulating the shape of a real sand grain. While variations of the particle shape was found to produce a relatively modest effect on the particle surface permeability, the orientation of the boundaries, emulating tortuosity effects, was found to produce a stronger impact due to the substantial variation of the distance, on average, between the boundary contours $\partial\Gamma_{1,2}$. It became clear that while the previously obtained scaling was a good first step to estimate the surface permeability of particulate porous media, a more general case of an ensemble of interconnected particles should be analysed to enhance the model predictive power and at the same time to estimate rigorously the effects of tortuosity of the surface flow in the particle assembly. In this study, we will simulate a general case of an ensemble of many particles linked by liquid bridges. We will concentrate on the bunch of spherical particles, but of different radii and randomly arranged in configurations. We compare the random pack configuration results with some symmetric case to estimate the effects of tortuosity and formulate practical recipes to apply the super-fast diffusion model.

II Microscopic model of the surface permeability of the elements

Microscopically, the liquid creeping flow through the surface roughness of each particle can be described by a local Darcy-like relationship (4) between the surface flux density ${\bf q}$ and averaged (over some area containing many surface irregularities) pressure in the grooves $u$ Yost-1998 ; Tuller-2000 . Assuming incompressibility of the liquid and that the liquid layer thickness is constant $\delta_{R}=const$, one has

[TABLE]

Equation (4) taking into account (6) can then be transformed into the Laplace-Beltrami equation defined on the surface $\Gamma$ of the particle

[TABLE]

Here, $\Delta_{\Gamma}$ designates the Laplace-Beltrami operator, which is defined on the surface element $\Gamma$ through the surface gradient $\nabla_{\Gamma}$ tangential to the surface. Formally, let ${\bf{n}}_{\Gamma}$ denote the unit normal to the surface $\Gamma$, Fig. 2. Then, one can define the surface gradient of a smooth function $u$ as $\nabla_{\Gamma}u:=\nabla u-(\nabla u\cdot{\bf{n}}_{\Gamma}){\bf{n}}_{\Gamma}$ and then the Laplace-Beltrami operator is defined as $\Delta_{\Gamma}u=\nabla_{\Gamma}\cdot\nabla_{\Gamma}u$.

The second assumption $\delta_{R}=const$ implies that the surface layer is fully saturated, that is its content is not changing on the particle surface. The approximation of the fully saturated rough surface layer is well fulfilled, if the characteristic pressure amplitude $|u|$ is less than the capillary pressure amplitude defined on the length scale of the surface roughness $\delta_{R}$, which is of the order of $\delta_{R}\sim 1\,\mu\mbox{m}$ in typical sands Alshibli2004 , as is demonstrated in Yost-1998 . That is, $|u|<u_{c}=\frac{\gamma}{\delta_{R}}$, and, for example for water ($\gamma=72\,\mbox{mN}/\mbox{m}$) at $\delta_{R}=1\,\mu\mbox{m}$, this results in $|u|<7.2\times 10^{4}\,\mbox{Pa}$.

Alternatively, if the surface layer somehow is not fully saturated, parameter $\delta_{R}$ should be interpreted as the characteristic width of the liquid layer within the rough surface layer and one needs to presume that variations of the pressure $|\delta u|$ are negligible $|\delta u|\ll u_{c}$. This is usually the case in slow, creeping flow conditions in porous media, and in fact, it is a criterion for the use of macroscopic approximation to such flows Bear-Book . As is shown in Lukyanov2019 , strong negative capillary pressure on the level of $u_{c}$ are only expected at the moving front, so that the approximation is well fulfilled in the macroscopic flow domain. Note also that, it is always assumed throughout this study that

[TABLE]

that is the amplitude of the surface roughness (or the width of the liquid layer) is always much smaller than the particle size.

II.1 Permeability of a single particle element

Consider, as the simplest example, a spherical particle of radius $R$ with a closed surface $\Gamma$, which is split into three sub-domains $\Omega_{0}$, $\Omega_{1}$ and $\Omega_{2}$ with the surface boundaries between them $\partial\Gamma_{1}$ and $\partial\Gamma_{2}$, as is shown in Fig. 2. The location of the sub-domains $\Omega_{1}$ and $\Omega_{2}$ to each other on the surface is fixed by the tilt angle $\alpha$. The sub-domains $\Omega_{1}$ and $\Omega_{2}$ correspond to the contact area covered by the liquid in the bridges, while the surface flow, described by (4), takes place in $\Omega_{0}$. Our prime concern is permeability of the surface elements, so that we only consider steady state problems.

The distribution of liquid pressure $u$, as it follows from (7), should satisfy the Laplace-Beltrami equation now defined on the surface of the sub-domain $\Omega_{0}$

[TABLE]

Note that, in fact, the condition of the fully saturated surface layer is not essential in calculation of the flows over one particle element of the porous media. It is sufficient to presume that the variation of the capillary pressure on the length scale of the particle $|\delta u|$ is negligible, that is $|\delta u|\ll u_{c}$. In the case when the surface layer is not fully saturated, parameter $\delta_{R}$ should be interpreted as the effective thickness of the layer filled by the liquid.

At the same time, liquid pressure variation in the bridges is negligible in slow, creeping flows in comparison to that in $\Omega_{0}$. So that, one can assume that

[TABLE]

which are the boundary conditions to the Laplace-Beltrami Dirichlet boundary value problem. The Dirichlet boundary value problem (8)-(9) has at least a unique weak solution, if the domain $\Omega_{0}$ and the boundaries $\partial\Gamma_{1,2}$ are smooth Dziuk1988 ; AMS2005 ; Dziuk2013 , which, if it is found, allows to calculate the total flux through the particle element

[TABLE]

where $\bf n_{s}$ is the normal vector to the domain boundaries $\partial\Gamma_{1,2}$ on the surface, $\delta_{R}$ is the average amplitude of the surface roughness, that is the width of the surface layer conducting the liquid flux and the line integral is taken along a closed curve in $\Omega_{0}$, for example the boundary $\partial\Gamma_{1}$.

If the total flux $Q_{T}$ is determined, one can define the global permeability coefficient of a single particle $K_{1}$. This can be done, if we assume that the particle has a characteristic size $D$ and so that it can be enclosed in a volume element $V=D^{3}$ with the characteristic side surface area $D^{2}$. Then, the effective flux density $Q$ can be represented in terms of $K_{1}$ (and the total flux $Q_{T}$)

[TABLE]

if the flow is driven by the constant pressure difference $U_{2}-U_{1}$ applied to the sides of the volume element.

II.2 Surface permeability of a sphere in the case of azimuthally symmetric domain boundaries

Consider now a spherical particle in an azimuthally symmetric case, when the domain boundaries $\partial\Gamma_{1}$ and $\partial\Gamma_{2}$ are oriented at the reflex angle $\alpha=\pi$ and have a circular shape. We use a spherical coordinate system with its origin at the particle centre and the polar angle $\theta$ counted from the axis of symmetry passing through the centre of the circular contour $\partial\Gamma_{1}$. In this case, the Dirichlet boundary value problem (8)-(9) admits an analytical solution, so that particle permeability can be determined explicitly. Indeed, problem (8)-(9), if we assume that the liquid pressure distribution $u$ is a function of $\theta$ only and independent of the azimuthal angle, is equivalent to

[TABLE]

with the boundary conditions

[TABLE]

The analytic solution to problem (12)-(13) after applying the boundary conditions can be represented in the following form

[TABLE]

where

[TABLE]

One can now calculate the total flux and the permeability, using its definition (11),

[TABLE]

[TABLE]

So that, taking $D=2R$,

[TABLE]

Parametrically, the coefficient of permeability (16) is inversely proportional to the particle radius $R$, so that larger particles create stronger resistance to the flow. Noticeably, the coefficient demonstrates strong dependence on the surface layer thickness $\delta_{R}$, that is $K_{1}\propto\delta_{R}^{3}$ since it is anticipated that $k_{m}\propto\delta_{R}^{2}$, so that evaluation of this parameter in applications is crucial for the accurate estimates of the liquid dispersion rates.

One can see, if we take $\theta_{1}=\theta_{0}$, in fact assuming small variations of the bridge size and the pressure over one particle diameter, and $\theta_{0}\ll 1$, in fact considering small values of saturation, $s\ll 1$, that the permeability coefficient $K_{1}$ tends to zero, that is

[TABLE]

How does the result affect the super-fast diffusion model (1), and basically how can it be incorporated into the main diffusion equation? If we approximate the permeability coefficient $K$ by $K_{1}$ obtained in the azimuthally symmetric case at $\theta_{1}=\theta_{0}$, (17), and, using an approximate relationship between the radius of curvature $R\sin\theta_{0}$ of the boundary contour $\partial\Gamma_{1}$ and the pendular ring volume Herminghaus-2005 , one can show that

[TABLE]

Therefore, finally

[TABLE]

As it follows from (18), the distinctive particle shape results in logarithmic correction to the main non-linear superfast-diffusion coefficient $D(s)=\frac{D_{0}(s)}{(s-s_{0})^{3/2}}$, such that

[TABLE]

Apparently, the correction will mitigate to some extent the divergent nature of the dispersion at the very small saturation levels $s\approx s_{0}$, smoothing out the characteristic dispersion curves.

II.3 Surface permeability of a chain of spheres in the case of azimuthally symmetric domain boundaries

Consider now how the problem can be formulated in the case of several particles arranged in a single chain, as is illustrated in Fig. 3 in the case of two coupled by the bridge particles. To create the flow in the system of two coupled particles, one can set pressure difference between $\partial\Gamma_{1}^{(1)}$ and $\partial\Gamma_{2}^{(2)}$. Mathematically, this is equivalent of setting Dirichlet boundary conditions on $\partial\Gamma_{1}^{(1)}$ and $\partial\Gamma_{2}^{(2)}$ as in the previous case of a single particle. The boundaries $\partial\Gamma_{1}^{(2)}$ and $\partial\Gamma_{2}^{(1)}$ are ’internal’, that is common to the bridge linking the flow between the two particles. Apparently, the pressure is supposed to be the same on the two contours

[TABLE]

and due to conservation of mass in steady state conditions in the absence of sinks and sources of the liquid one has

[TABLE]

where ${\bf n}_{s_{1}}$ and ${\bf n}_{s_{2}}$ are the outward tangential normal vectors to the boundary contours $\partial\Gamma_{1,2}^{(2,1)}$, and $u_{1}$ and $u_{2}$ designate distribution of pressure on each particle respectively.

As a result, the problem to define the flow and the permeability of the system corresponds to a system of two Laplace-Beltrami equations

[TABLE]

and

[TABLE]

but with a slightly different set of the boundary conditions

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

where $\theta_{0}$, as before, defines the size of the bridge footprint on the particle surface in the spherical coordinate system with the axis of symmetry passing through the centre of the bridge area, Fig. 3. Since we assumed, due to relatively small variations of pressure over a few grain particles, that all bridges are roughly identical, we have only one parameter $\theta_{0}$ to describe the bridge size.

Apparently, equations (21) and (22) can be integrated twice, similar to the previous problem of a single particle (12), to obtain

[TABLE]

[TABLE]

where $C_{1,2}$ and $B_{1,2}$ are free constant parameters to be found from the boundary conditions.

It is not difficult to see from (26), that one has $C_{0}=B_{0}$ implying continuity of the contact flux. Applying the remaining boundary conditions (23)-(25), from (27) and (28)

[TABLE]

[TABLE]

where

[TABLE]

One can now calculate total flux and define permeability of the coupled spherical particles $K_{2}$

[TABLE]

[TABLE]

where $D$ is the characteristic length scale of the cross-section in the problem, $\kappa_{m}$ is local permeability of the surface layer, $\delta_{R}$ is the layer width and $\mu$ is liquid viscosity.

So that, taking simply $D=2R$,

[TABLE]

One can see that the permeability of a system of two coupled particles $K_{2}$ is identical to that of a single particle (16), basically from (16) and (32)

[TABLE]

It is not difficult to discern by deduction that in a general case of $N$ coupled particles in a chain

[TABLE]

where

[TABLE]

Note, experimentally, the setup of many beads coupled by liquid bridges is often used in microfluidics to create flexible water channels Nanoscale2016 . If the radius of curvature of the particle chain is much larger than the particle size, the transport through such a microfluidic system should be defined by the permeability of a single particle, relationship (16), if the particle shape can be approximated by a sphere.

One can conclude in this part, that if the porous media configuration is made of parallel chains of particles oriented symmetrically to each other, and the flow is generated along the chains, the surface permeability given by (16) is the exact result.

III Surface permeability of a randomly packed particle ensemble

In real systems, the particles are interconnected randomly, so that the effects of tortuosity should substantially affect the permeability of the system Tortuosity1937 ; Tortuosity1961 ; Bear-Book ; Tortuosity-Review2013 . To analyse those effects, we consider an ensemble of spherical particles randomly packed, as is shown in Fig. 4. The randomly packed configuration of approximately $3000-7000$ particles has been generated by means of a molecular dynamics technique by applying a constant force to every particle placed in a box with reflecting boundaries (in the perpendicular direction to the box side), and interacting via the Lennard-Jones potential with different characteristic length scales $R$ distributed normally, that is with the probability of the particle radius $W(R)\propto\exp\left(-\frac{(R-R_{0})^{2}}{\Delta R^{2}}\right)$ at $\Delta R/R_{0}=0.3$. In this study, there were particles with three different characteristic dimensions $R_{1}=1.3\,R_{0}$, $R_{2}=R_{0}$ and $R_{3}=0.7\,R_{0}$. The resultant porosity in the configurations was about $48\%$.

To obtain the configuration, the particle temperature controlled by the thermostat has been gradually reduced to bring the system to a minimum energy, frozen state. A representative sample volume element with dimensions $L_{x}^{B},L_{y}^{B},L_{z}^{B}$ then was cut off the system, as is illustrated in Fig. 4, containing $N_{S}=13-17$ particles, see Table 1 for details. We have generated several statistically independent sample configurations, and, as in the previous examples, set constant pressure difference $U_{2}-U_{1}$ at the boundaries of the sample elements, Figs. 4 and 6.

The Laplace-Beltrami method then has been applied after establishing the position of the liquid bridges coupling the particles in the sample. Two particles (of radii $R_{1}$ and $R_{2}$) are assumed to be coupled by a liquid bridge if the distance between their centres $r$ was only slightly larger than the sum of their radii

[TABLE]

The size of a single liquid bridge footprint $H_{B}$ on the particle surface can be characterised, as before, by the polar angle $\theta_{0}$ in the polar coordinate system with the symmetry axis passing through the centre of the circular contour, the boundary of the area covered by the bridge, as is shown in the symmetric case in Fig. 3. That is, $H_{B}=2R_{k}\sin\theta_{0}^{(k)}$. Due to the specific geometric properties of the pendular rings (constant mean curvature surface), we assume that even in the case of a distribution of particles with different radii $R_{k}$, the size of the bridge area in the sample is approximately the same in the low saturation limit $s\ll 1$ ($\theta_{0}^{(k)}\ll 1$) Herminghaus-2005 ; Halsey1998 ; Willett2000 .

Indeed, when $s\ll 1$, the pressure in the pendular ring $p$ is defined by the smallest radius of curvature $r_{1}$, Fig. 5, $p\approx-\gamma\cos\phi_{c}/r_{1}$, which is related with the second radius

[TABLE]

so that when $s\ll 1$, one has $r_{2}\gg r_{1}$. Obviously, $r_{2}$ defines the size of the area covered by the bridge, $H_{B}=2r_{2}=2R_{k}\sin\theta_{0}^{(k)}$.

If we have two particles of different radii, say $R_{1}$ and $R_{2}$, in contact, the size of the bridge area will be approximately the same $r_{2}^{(1)}\approx r_{2}^{(2)}$ at low saturation levels, $s\ll 1$, with the difference being proportional to $r_{1}$, that is

[TABLE]

Apparently, in a general case, no analytic solution is expected to the Laplace-Beltrami problem and a well established surface finite element technique Dziuk1988 ; AMS2005 ; Dziuk2013 is applied after the tessellation of the domains, as is shown in Fig. 7 for one particulate element with two boundary contours. The numerical method has been validated against analytical solutions previously demonstrating prescribed order of accuracy and numerical convergence, see details in Penpark2018 .

The number of particles in the sample volume element was negotiated between computational efficiency of the surface finite element method (so that, practically, any mesh resolution can be used to deal with any details on the boundary contours $\partial\Gamma_{k}^{(l)}$ and on the particle surfaces) and fluctuations of the averaged quantities obtained using the sample element, which are proportional to $N_{S}^{-1/2}\approx 25\%$. Moderate increase of the number of particles in the sample may significantly increase computational time to obtain highly resolved numerical solutions, while at the same time would not substantially reduce the effect of particle number fluctuations.

As one can see, problem (21) - (26) and hence total flux $Q_{T}$ through a particle or a chain of particles, (15) or (31), are invariant under the transformation of the particle dimension $R$ provided that the angular size of the bridge $\theta_{0}$ is fixed. In what follows, we change to non-dimensional description by normalising length scales by the average radius $R_{0}$ of the particles in the sample and pressure by the characteristic capillary pressure $p_{0}=2\gamma\cos\phi_{c}/R_{0}$. The flux $Q_{T}$ will be normalised by the characteristic value

[TABLE]

inspired by the analytical result (15) and by the non-dimensional sample box surface area $S_{0}=\bar{L}_{x}^{B}\bar{L}_{y}^{B}$, where non-dimensional quantities $\bar{L}_{x,y,z}^{B}=L_{x,y,z}^{B}/R_{0}$ and $\bar{U}_{1,2}=U_{1,2}/p_{0}$. The latter normalisation allows to bring simulation results in slightly different geometric settings, as is detailed in Table 1, into equivalent conditions suitable for comparison, that is basically providing the non-dimensional permeability $\bar{K}=\frac{K}{\kappa_{m}}\frac{R_{0}}{\delta_{R}}$.

Schematically, the simulation domains for the Laplace-Beltrami problem are shown in Fig. 6. As in the previous case of particles coupled in a chain, there are internal, common boundaries, contours $\partial\Gamma_{k}^{(l)}$, $k\neq l$, where the continuity boundary conditions are applied and external boundaries, contours $\partial\Gamma_{k}^{(k)}$, where the Dirichlet boundary conditions are set to generate a flow through the system. The pressure value on the contours facing the bottom of the simulation box (for example, $\partial\Gamma_{2}^{(2)},\partial\Gamma_{3}^{(3)}$ and $\partial\Gamma_{4}^{(4)}$ in Fig. 6) is set to $U_{2}$ and on the contours facing the top side of the pack (for example, $\partial\Gamma_{0}^{(0)}$ and $\partial\Gamma_{1}^{(1)}$ in Fig. 6) is set to $U_{1}$. The values of the boundary pressure $U_{1,2}$ were identical in the simulations involving different configurations.

Geometrically, the external boundary contours are oriented in the flow direction, as is illustrated in Fig. 6. While this particular orientation seems to be arbitrary or may even look artificial, within the statistical approach, the choice of the boundary contour orientation should not render any excessive (in excess of the statistical errors due to the particle number fluctuations) influence upon the results, that is the value of the total flux and the ’macroscopic’ permeability. A posteriori, one can see that this seemed to be the case, Fig. 8, as in different configurations, Table 1, the resultant curves are close and parallel to each other.

There are two main questions, we would like to answer in this part of the study. First, how does permeability of the particle sample depend on the composition? Basically, how strong are there fluctuations? Secondly, what is the contribution of the tortuosity effects? To obtain statistically meaningful results, we consider several randomly generated configurations, as is summarised in Table 1. We would like to stress here, that all configurations have been cut off from statistically independent particle distributions generated with the help of random initial distributions of larger number of particles, as we have described.

As before, we are going to find a weak solution to a system of the Laplace-Beltrami equations

[TABLE]

defined on each particle domain $\Omega^{(k)}_{0}$, as in in Fig. 6. On the internal boundaries of the domains we set up continuity conditions, for example on $\partial\Gamma_{3}^{(2)}$ and $\partial\Gamma_{2}^{(3)}$

[TABLE]

[TABLE]

While on a few external boundaries, Dirichlet boundary conditions are set.

The numerical solution allows to calculate the total flux through the system by summing up the fluxes passing through the external contours, where the Dirichlet boundary conditions are set, either at the top of the pack or at the bottom using (10). The results are summarised in Fig. 8.

Remarkably, the reduced flux $Q_{T}/S_{0}Q_{0}$ as a function of

[TABLE]

where $H_{B}$ is the bridge size $H_{B}=2R_{0}\sin\theta_{0}$, behaves linearly in all configurations. This behaviour mirrors the flux dependence observed in azimuthally symmetric analytical solutions, see (15) or (31). The variations in the dependencies between different configurations are observed to be well within the statistical error expected in this case, error bar in Fig 8. At the same time, a comparison with a similar, but a regular arrangement, as in Fig. 3 at $R=R_{0}$ demonstrates that there is a clear cut contribution from the effects of tortuosity, solid line in Fig. 8.

Indeed, given identical porosity ($\phi\approx 50\%$) and mean particle size ($R/R_{0}=1$) in the regular, symmetric and randomly generated configurations, the normalised flux values differ by a factor of two, which is consistent with the tortuosity values obtained in porous media in different conditions and configurations Tortuosity-Review2013 . For example, both hydraulic $\tau_{h}$ and diffusive $\tau_{d}$ tortuosity estimated in unsaturated porous media using different permeability models (often used in applications, for example Mualem1976 ; Mualem1978 ) was found in between $1.5\leq\tau_{h,d}\leq 2$ at $\phi=50\%$ Tortuosity-Review2013 .

It is important that the result, that is the ratio of the total flux in the random and regular configurations does not practically depend on the size of the contour $H_{B}$, basically the size of the liquid bridge, and hence the value of saturation in the porous media. This implies, that the observed effect is purely down to the distribution of contacts between the particles, but not the particular pathway on each single particle surface. That is, fundamentally, tortuosity in the surface diffusion processes is a geometric factor independent of the particular surface flow regime. At the same time, the pathways, on average, of course, does depend on the bridge size value $H_{B}$ leading to smaller permeability as the size of the contact area diminishes. This trend is expected, but essentially, the correction to the effective coefficient of diffusion

[TABLE]

is only down to a single universal factor of two representing the tortuosity effects in surface diffusion in particular porous media at low values of saturation. Note, this value is also in agreement with experimental observations and a comparison of the super-fast diffusion model with the data, where the tortuosity effects were estimated to reduce the effective permeability twofold Lukyanov2019 .

This is the main result of this study, which can be used in practical applications to calculate permeability in particular porous media. Basically, as the first step, one can calculate permeability of a single, representative element of the media or several elements to obtain some mean value and its dispersion. This way permeability $K(s)$, via (18), and the diffusion coefficient $D(s)$ for the macroscopic model can be established in the first approximation. Macroscopic permeability $K(s)$ or the diffusion coefficient $D(s)$ then should be corrected by the universal factor of two in the macroscopic diffusion model.

IV Conclusions

We have demonstrated that the Laplace-Beltrami method can be used to obtain permeability of particulate porous media at low saturation levels and to estimate contribution from the effects of tortuosity. Essentially, analytical results obtained using azimuthally oriented coupled particles can be used with a universal correcting prefactor to estimate permeability of particle ensembles. That is, from the practical point of view, results obtained by analysing single representative element of particulate porous media can be translated into permeability of a particle composition.

PS was supported through the Royal Thai Government scholarship. TP was partially supported through the EPSRC grant EP/P000835/1.