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 superfast nonlinear 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 directly estimate tortuosity effects of the surface flows in the particle ensemble setting.
Similar content being viewed by others
1 Introduction
Liquid transport in particulate porous media, such as sand, is customarily classified into fully saturated, funicular and pendular regimes of spreading (Bear 1972; Herminghaus 2005; Scheel et al. 2008a, b). 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\le 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; Scheel et al. 2008a, b). 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.
The analysis of this regime of wetting, which is crucial for studies of biological processes and spreading of nonvolatile 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 nonlinear diffusion equation (Lukyanov et al. 2012, 2019).
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 Lukyanov et al. (2019). 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. Note that, unlike at high saturation values \(s>s_c\), at low saturation \(s<s_c\) the capillary pressure is developed on the length scale of the surface roughness.
Theoretically, the superfast nonlinear diffusion equation belongs to a special class of mathematical models. Unlike in the standard porous medium equation (Vazquez 2006), where the coefficient of diffusion D(s) is either constant (leading to normal, Gaussian dispersive behaviour) or a monotonically vanishing function of saturation (\(D(s)\propto s^m, \, m>0\), leading to hypodispersive behaviour), in this special case, the nonlinear coefficient of diffusion demonstrates divergent behaviour as a function of saturation, \(D(s)\propto (s-s_0)^{-3/2}\), leading to hyperdispersive character of the spreading process. Here, \(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 (Lukyanov et al. 2012; Scheel et al. 2008a, b; Lukyanov et al. 2019). 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 (Lukyanov et al. 2012, 2019). The character of the dispersion process may often change from hypodispersive to hyperdispersive depending on the spreading scenario, that is on the driving forces of the spreading process (de Gennes 1985; Bacri et al. 1985; Novy et al. 1989; Toledo et al. 1993; Meleá 2003), especially in multi-phase systems during forced imbibition cycles used in petroleum engineering (Meleá 2003). If the forcing, external pressure at the system inlet is strong enough, the low saturation zone may not have a chance to be developed, and hypodispersive character of the spreading front is observed.
Returning to our mathematical model at low saturation values, 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 nonlinear diffusion equation
where
for \(D_0 > 0\).
The details of derivation of (1) can be found in Lukyanov et al. (2012) and Lukyanov et al. (2019), here we note that, the resultant governing nonlinear Eq. (1) directly follows from the conservation of mass principle
where \(\phi \) is porosity defined as \(\displaystyle \phi =\frac{V_E}{V}\), which is further assumed to be constant, and \(\mathbf Q\) is the macroscopic flux density. The macroscopic flux density \(\mathbf 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 \mathbf{\!\! Q\cdot n}\, \hbox {d}S\), where \(\mathbf 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 (Halsey and Levine 1998; Lukyanov et al. 2012, 2019) dictated by the liquid bridges behaviour
and the local Darcy’s law (Rye et al. 1998; Or and Tuller 2000) describing the surface flow in the rough layer of the particle elements
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, \(\mathbf{q}\) and u are the averaged local flux density and pressure in the rough surface layer, \(\mu \) is liquid viscosity, and \(\kappa _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
One needs to emphasise here that two levels of averaging are involved in obtaining the final governing Eq. (1). While Eqs. (1)–(3) are ‘truly’ macroscopic, that is obtained by averaging using a volume element V containing many grain particles, Eq. (4) is only an average over some rough area of a single particle containing many surface irregularities, so that quantities \(\mathbf{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 \cdots \rangle ^l=V_l^{-1}\!\int _{V_l} \ldots \hbox {d}^3 x\), where \(V_l\) is liquid volume within the sample volume V, one has \(\langle u \rangle ^l=p\) and \(\langle \mathbf{q}\rangle ^l \frac{S_e}{S} = \mathbf{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
where \(K(s)=\kappa _m \frac{S_e}{S}\) is the coefficient of permeability defined by
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, Figs. 1 and 2.
In particular, the disposition and the size of the liquid bridges on the particle surface, that is the size of the domains \(\varOmega _{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 \(\varOmega _{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 \(\varOmega _0\) with the boundaries \(\partial \varGamma _{1,2}\) in Fig. 2 (Sirimark et al. 2018). 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 \(\varOmega _1\) and \(\varOmega _2\) in Fig. 2, are oriented symmetrically to each other, that is at \(\alpha =\pi \), the permeability K is supposed to follow the scaling
We have studied several generalisations of the symmetric problem, such as arbitrary oriented domains, \(\alpha \ne \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 were 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 \varGamma _{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 superfast diffusion model.
2 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 \(\mathbf{q}\) and averaged (over some area containing many surface irregularities) pressure in the grooves u (Rye et al. 1998; Or and Tuller 2000). Assuming incompressibility of the liquid and that the liquid layer thickness is constant \(\delta _R=\mathrm{const}\), one has
Equation (4) taking into account (6) can then be transformed into the Laplace–Beltrami equation defined on the surface \(\varGamma \) of the particle
here \(\Delta _{\varGamma }\) designates the Laplace–Beltrami operator, which is defined on the surface element \(\varGamma \) through the surface gradient \(\nabla _{\varGamma }\) tangential to the surface. Formally, let \(\mathbf{{n}}_{\varGamma }\) denote the unit normal to the surface \(\varGamma \), Fig. 2. Then, one can define the surface gradient of a smooth function u as \(\nabla _{\varGamma } u:=\nabla u - (\nabla u \cdot \mathbf{{n}}_{\varGamma }) \mathbf{{n}}_{\varGamma }\), and then, the Laplace–Beltrami operator is defined as \(\Delta _{\varGamma } u = \nabla _{\varGamma }\cdot \nabla _{\varGamma } u\).
The second assumption \(\delta _R=\hbox {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\,\upmu \text{ m }\) in typical sands (Alshibli and Alsaleh 2004), as is demonstrated in Rye et al. (1998). That is, \(|u| < u_c=\frac{\gamma }{\delta _R}\), and, for example for water (\(\gamma =72\, \text{ mN }/\text{ m }\)) at \(\delta _R =1\,\mu \text{ m }\), this results in \(|u| < 7.2\times 10^4\, \text{ 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 1972). As is shown in Lukyanov et al. (2019), strong negative capillary pressure on the level of \(u_c\) is 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
that is the amplitude of the surface roughness (or the width of the liquid layer) is always much smaller than the particle size.
2.1 Permeability of a Single Particle Element
Consider, as the simplest example, a spherical particle of radius R with a closed surface \(\varGamma \), which is split into three sub-domains \(\varOmega _0\), \(\varOmega _1\) and \(\varOmega _2\) with the surface boundaries between them \(\partial \varGamma _1\) and \(\partial \varGamma _2\), as is shown in Fig. 2. The location of the sub-domains \(\varOmega _1\) and \(\varOmega _2\) to each other on the surface is fixed by the tilt angle \(\alpha \). The sub-domains \(\varOmega _1\) and \(\varOmega _2\) correspond to the contact area covered by the liquid in the bridges, while the surface flow, described by (4), takes place in \(\varOmega _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 \(\varOmega _0\)
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 with that in \(\varOmega _0\). So that, one can assume that
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 \(\varOmega _0\) and the boundaries \(\partial \varGamma _{1,2}\) are smooth (Dziuk 1988; Pigola et al. 2005; Dziuk and Elliott 2013), which, if it is found, allows to calculate the total flux through the particle element
where \(\mathbf n_s\) is the normal vector to the domain boundaries \(\partial \varGamma _{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 \(\varOmega _0\), for example the boundary \(\partial \varGamma _1\).
If the total flux \(Q_\mathrm{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_\mathrm{T}\))
if the flow is driven by the constant pressure difference \(U_2-U_1\) applied to the sides of the volume element.
2.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 \varGamma _1\) and \(\partial \varGamma _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 \varGamma _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
with the boundary conditions
The analytic solution to problem (12)–(13) after applying the boundary conditions can be represented in the following form
where
One can now calculate the total flux and the permeability, using its definition (11),
So that, taking \(D=2R\),
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 \(\kappa _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
How does the result affect the superfast 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 \varGamma _1\) and the pendular ring volume (Herminghaus 2005), one can show that
Therefore, finally
As it follows from (18), the distinctive particle shape results in logarithmic correction to the main nonlinear superfast diffusion coefficient \(D(s)=\frac{D_0(s)}{(s-s_0)^{3/2}}\), such that
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.
2.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 \varGamma _1^{(1)}\) and \(\partial \varGamma _2^{(2)}\). Mathematically, this is equivalent of setting Dirichlet boundary conditions on \(\partial \varGamma _1^{(1)}\) and \(\partial \varGamma _2^{(2)}\) as in the previous case of a single particle. The boundaries \(\partial \varGamma _1^{(2)}\) and \(\partial \varGamma _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
and due to conservation of mass in steady-state conditions in the absence of sinks and sources of the liquid one has
where \(\mathbf{n}_{s_1}\) and \(\mathbf{n}_{s_2}\) are the outward tangential normal vectors to the boundary contours \(\partial \varGamma _{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
and
but with a slightly different set of the boundary conditions
and
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, Eqs. (21) and (22) can be integrated twice, similar to the previous problem of a single particle (12), to obtain
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)
where
One can now calculate total flux and define permeability of the coupled spherical particles \(K_2\)
where D is the characteristic length scale of the cross section in the problem, \(\kappa _m\) is local permeability of the surface layer, and \(\delta _R\) is the layer width and \(\mu \) is liquid viscosity.
So that, taking simply \(D=2R\),
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 Eqs. (16) and (32)
It is not difficult to discern by deduction that in a general case of N coupled particles in a chain
where
Note, a setup of many beads coupled by liquid bridges can be potentially used in microfluidics to create flexible water channels (Chen et al. 2016). 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 Eq. (16) is the exact result.
3 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 (Carman 1937; Lorenz 1961; Bear 1972; Ghanbarian et al. 2013). 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 5.
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
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; Halsey and Levine 1998; Willett et al. 2000).
Indeed, when \(s\ll 1\), the pressure in the pendular ring p is defined by the smallest radius of curvature \(r_1\), Fig. 6, \(p\approx - \gamma \cos \phi _c/r_1\), which is related with the second radius
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
Apparently, in a general case, no analytic solution is expected to the Laplace–Beltrami problem and a well established surface finite element technique (Dziuk 1988; Pigola et al. 2005; Dziuk and Elliott 2013) 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 Sirimark et al. (2018).
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 \varGamma _{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, problems (21)–(26) and hence total flux \(Q_\mathrm{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_\mathrm{T}\) will be normalised by the characteristic value
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. 5. As in the previous case of particles coupled in a chain, there are internal, common boundaries, contours \(\partial \varGamma _{k}^{(l)}\), \(k\ne l\), where the continuity boundary conditions are applied and external boundaries, contours \(\partial \varGamma _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 \varGamma _2^{(2)}, \partial \varGamma _3^{(3)}\) and \(\partial \varGamma _4^{(4)}\) in Fig. 5) is set to \(U_2\) and on the contours facing the top side of the pack (for example, \(\partial \varGamma _0^{(0)}\) and \(\partial \varGamma _1^{(1)}\) in Fig. 5) 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. 5. 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
defined on each particle domain \(\varOmega ^{(k)}_0\), as in Fig. 5. On the internal boundaries of the domains, we set up continuity conditions, for example on \(\partial \varGamma _3^{(2)}\) and \(\partial \varGamma _2^{(3)}\)
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_\mathrm{T}/ S_0 Q_0\) as a function of
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 \(\tau _s\approx 2\), which can be regarded as the tortuosity of the surface flow. The value of the surface flow tortuosity \(\tau _s\approx 2\) is consistent with the tortuosity values obtained in porous media in different conditions and configurations (Ghanbarian et al. 2013). For example, both hydraulic \(\tau _h\) and diffusive \(\tau _d\) tortuosities estimated in unsaturated porous media using different permeability models [often used in applications, for example Mualem (1976, 1978)] were found in between \(1.5 \le \tau _{h,d}\le 2\) at \(\phi =50\%\) (Ghanbarian et al. 2013).
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
is only down to a single universal factor of two, \(\tau _s\), representing the tortuosity effects in surface diffusion in particular porous media at low values of saturation, where \(K_1(s)\propto |\ln (s-s_0)|^{-1}\) defined by (18) is the surface permeability of a single porous matrix element. Note, the value of \(\tau _s\approx 2\) is also in agreement with experimental observations and a comparison of the superfast diffusion model with the data, where the tortuosity effects were estimated to reduce the effective permeability twofold (Lukyanov et al. 2019).
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, \(\tau _s\), in the macroscopic diffusion model.
4 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 \(\tau _s\), as in (39), 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.
References
Alshibli, K.A., Alsaleh, M.I.: Characterizing surface roughness and shape of sands using digital microscopy. J. Comput. Civ. Eng. 18, 36–45 (2004)
Bacri, J.C., Leygnac, C., Salin, D.: Evidence of capillary hyperdiffusion in two-phase fluid flows. J. Phys. Lett. 46, 467–473 (1985)
Bear, J.: Dynamics of Fluids in Porous Media. Dover, New York (1972)
Carman, P.C.: Fluid flow through granular beds. Trans. Inst. Chem. Eng. 15, 150–166 (1937)
Chen, J., Wang, C., Wei, N., Wan, R., Gao, Y.: 3D flexible water channel: stretchability of nanoscale water bridge. Nanoscale 8, 5676–5681 (2016)
de Gennes, P.G.: Partial filling of a fractal structure by a wetting Fluid. In: Adler, D., Fritzsche, E., Ovshirisky, S.R. (eds.) Physics of Disordered Material, pp. 227–241. Plenum Press, New York (1985)
Dziuk, G.: Finite elements for the Beltrami operator on arbitrary surfaces. Partial Differ. Equ. Calc. Var. 1357, 142–155 (1988)
Dziuk, G., Elliott, C.M.: Finite element methods for surface PDEs. Acta Numer. 22, 289–396 (2013)
Ghanbarian, B., Hunt, A.G., Ewing, R.P., Sahimi, M.: Tortuosity in porous media: a critical review. Soil Sci. Soc. Am. J. 77, 1461–1477 (2013)
Halsey, T.C., Levine, A.J.: How sandcastles fall. Phys. Rev. Lett. 80, 3141–3144 (1998)
Herminghaus, S.: Dynamics of wet granular matter. Adv. Phys. 54, 221 (2005)
Lorenz, P.B.: Tortuosity in porous media. Nature 189, 386–387 (1961)
Lukyanov, A.V., Sushchikh, M.M., Baines, M.J., Theofanous, T.G.: Superfast nonlinear diffusion: capillary transport in particulate porous media. Phys. Rev. Lett. 109, 214501 (2012)
Lukyanov, A.V., Mitkin, V.V., Theofanous, T.G., Baines, M.J.: Capillary transport in particulate porous media at low levels of saturation. J. Appl. Phys. 125, 185301 (2019)
Meleán, Y., Broseta, D., Blossey, R.: Imbibition fronts in porous media: effects of initial wetting fluid saturation and flow rate. J. Pet. Sci. Eng. 39, 327–336 (2003)
Mualem, Y.: A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12, 513–522 (1976)
Mualem, Y.: Hydraulic conductivity of unsaturated porous media: generalized macroscopic approach. Water Resour. Res. 14, 325–334 (1978)
Novy, R.A., Toledo, P.G., Davis, H.T., Scriven, L.E.: Capillary dispersion in porous media at low wetting phase saturations. Chem. Eng. Sci. 44, 1785–1797 (1989)
Or, D., Tuller, M.: Flow in unsaturated fractured porous media: hydraulic conductivity of rough surfaces. Water Resour. Res. 36, 1165–1177 (2000)
Pigola, S., Rigoli, M., Setti, A.G.: Maximum principles on Riemannian manifolds and applications. Mem. Am. Math. Soc. 174, 1–36 (2005)
Rye, R.R., Yost, F.G., O’Toole, E.J.: Capillary flow in irregular surface grooves. Langmuir 14, 3937 (1998)
Scheel, M., Seemann, R., Brinkmann, M., Michiel, M.D.I., Sheppard, A., Breidenbach, B., Herminghaus, S.: Morphological clues to wet granular pile stability. Nat. Mater. 7, 189 (2008a)
Scheel, M., Seemann, R., Brinkmann, M., Michiel, M.D.I., Sheppard, A., Herminghaus, S.: Liquid distribution and cohesion in wet granular assemblies beyond the capillary bridge regime. J. Phys. Condens. Matter 20, 494236 (2008b)
Sirimark, P., Lukyanov, A.V., Pryer, T.: Surface permeability of porous media particles and capillary transport. Eur. Phys. J. E 41, 106 (2018)
Toledo, P.G., Davis, H.T., Scriven, L.E.: Capillary hyperdispersion of wetting liquids in fractal porous media. Transp. Porous Med. 10, 81–94 (1993)
Vazquez, J.L.: The Porous Medium Equation: Mathematical Theory. Oxford University Press, Oxford (2006)
Whitaker, S.: Advances in theory of fluid motion in porous media. Ind. Eng. Chem. 61, 14–28 (1969)
Willett, C.D., Adams, M.J., Johnson, S.A., Seville, J.P.K.: Capillary bridges between two spherical bodies. Langmuir 16, 9396–9405 (2000)
Acknowledgements
PS was supported through the Royal Thai Government scholarship. TP was partially supported through the EPSRC Grant EP/P000835/1.
Author information
Authors and Affiliations
Contributions
All the authors were involved in the preparation of the manuscript. All the authors have read and approved the final manuscript.
Corresponding author
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
About this article
Cite this article
Sirimark, P., Lukyanov, A.V. & Pryer, T. Surface Permeability of Particulate Porous Media. Transp Porous Med 130, 637–654 (2019). https://doi.org/10.1007/s11242-019-01332-9
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s11242-019-01332-9