The breakdown of Darcy's law in a soft porous material

TL;DR

This study reveals that flow through deformable porous media deviates from Darcy's law, showing nonlinear behavior influenced by material softness, and introduces a universal relation involving pressure difference and shear modulus.

Contribution

The paper develops a theoretical model explaining the super-linear flow behavior in soft porous materials and confirms a universal relation involving pressure and shear modulus.

Findings

Flow flux is nonlinear and steeper than linear with pressure difference.

Flow permeability depends on the ratio of pressure difference to shear modulus.

Theoretical predictions match numerical simulation results.

Abstract

We perform direct numerical simulations of the flow through a model of a deformable porous medium. Our model is a two-dimensional hexagonal lattice, with defects, of soft elastic cylindrical pillars, with elastic shear modulus $G$, immersed in a liquid. We use a two-phase approach: the liquid phase is a viscous fluid and the solid phase is modeled as an incompressible viscoelastic material, whose complete nonlinear structural response is considered. We observe that the Darcy flux ($q$) is a nonlinear function -- steeper than linear -- of the pressure-difference ($\Delta P$) across the medium. Furthermore, the flux is larger for a softer medium (smaller $G$). We construct a theory of this super-linear behavior by modelling the channels between the solid cylinders as elastic channels whose walls are made of a material with a linear constitutive relation but can undergo large deformation. Our theory further predicts that the flow permeability is a universal function of $\Delta P/G$, which is confirmed by the present simulations.

Figures

Equations

$$q=-\frac{k}{\mu }\frac{\Delta P}{L}\text{/},$$

$$\frac{\partial u_i}{\partial t}+\frac{\partial u_i{u}_j}{\partial x_j}=\frac{1}{\rho }\frac{\partial {\sigma }_{ij}}{\partial x_j}\text{and}\frac{\partial u_i}{\partial x_i}=0,$$

$${\sigma }_{ij}=\varphi {\sigma }_{ij}^{\mathrm{f}}+\left(1-\varphi \right){\sigma }_{ij}^{\mathrm{s}}\text{/},$$

$${\sigma }_{ij}^{\mathrm{f}}$$

$${\sigma }_{ij}^{\mathrm{s}}$$

$$\frac{\partial {\mathcal{B}}_{ij}}{\partial t}+\frac{\partial u_k{\mathcal{B}}_{ij}}{\partial x_k}$$

$$\frac{\partial \varphi }{\partial t}+\frac{\partial u_k\varphi }{\partial x_k}$$

$$q=-[h^3/(12\mu )](dp/d\xi )\text{/}.$$

$$\mathcal{K}=h_0^{3}f(\beta ),\text{mbox}where$$

$$\nabla \cdot \left[{\mathcal{K}}_{\mathcal{L}}({\Phi }_{\mathcal{L}})\nabla p\right]=0,$$

$$\frac{1}{k}=\int d{\Phi }_{\mathcal{L}}\frac{P({\Phi }_{\mathcal{L}})}{{\mathcal{K}}_{\mathcal{L}}({\Phi }_{\mathcal{L}})}\text{/},$$

Full text

I Introduction

Percolation of water through soil is one of the oldest problems in hydrodynamics. The fluid passes through a network of irregularly arranged interstices between solid objects. Typically, each individual thread of water passes through a narrow channel in which the equations of viscous flow can be applied. The difficulty arises from the fact that the detailed knowledge of the channels is neither available nor useful due to their complexity. We therefore typically take coarse-grained approaches, averaging over a length-scale much larger than the individual channels but still small compared to the scale of the medium. We thus define a relation between the flux, $q$, and the pressure-difference, $\Delta P$. In the simplest case of a rigid isotropic medium, this gives rise to Darcy’s law Batchelor (1967)

[TABLE]

where $\mu$ is the dynamic viscosity of the fluid, $k$ is the permeability of the porous medium, and $L$ its length in the flow direction. Henceforth we shall call $q$ the Darcy flux. The permeability has the same status as all transport coefficients in hydrodynamics – for a real system it is very difficult to calculate from first principles, but can be calculated in a model system by first solving a problem at the pore-scale and then by either analytical or numerical coarse-graining. This problem develops an additional degree of complexity if we consider that under the fluid stress the solid obstacles can move, i.e., the flow itself can form channels. We treat the complication wherein the solid skeleton is deformable, i.e., poroelasticity, the simplest example of which is the kitchen sponge.

Poroelasticity play an important role in understanding the transport through a wide range of materials ranging from individual cells Charras et al. (2005); Moeendarbary et al. (2013), to biological tissues, e.g., soft-tissues Lai et al. (1991); Yang and Taber (1991); Auton and MacMinn (2017), bones Cowin (1999), even to hydraulic fracture Detournay and Cheng (1988); Yarushina et al. (2013). The simplest poroelastic problem is that of linear poroelasticity where we assume that the flow of the liquid is governed by the Darcy’s law and the solid skeleton not only has linear constitutive relation but also undergoes small deformation. In reality, often the deformation of the solid matrix is large consequently nonlinear elastic effects have to be taken into account even if the constitutive relation is linear. Such systems are notoriously difficult to study both experimentally and numerically MacMinn et al. (2015).

The central question in this paper is how a coarse-grained description of the Darcy type emerges from a pore-scale model. As our model we choose a bed, a two-dimensional hexagonal lattice with defects, of soft elastic cylinders immersed in a liquid. Using both direct numerical simulations – a set of fully coupled equations for a viscoelastic solid in contact with a Newtonian fluid – and theory, we show that at scales that are large compared to the diameter of a cylinder the flux versus pressure-difference relationship in the system is a Darcy-like equation. When the deformability is small, as measured by the shear modulus of the solid, we obtain the Darcy equation exactly: the permeability $k$ is a constant, independent of the pressure-difference. However, as the solid skeleton becomes more deformable, the permeability becomes a nonlinear function of the pressure difference, namely for the same pressure drop we get a larger flux. Our theoretical calculations suggest that this result is largely model independent. This behavior has been already predicted from theoretical modelling at a coarse-grained level MacMinn et al. (2016) but has never been observed before in simulations or experiments.

II Numerical method

We first describe our Direct Numerical Simulations (DNS). The deformable cylinders in the Newtonian fluid are modeled with a two-phase approach, defined by a variable $\phi=0$ inside the viscoelastic solid phase and $\phi=1$ in the fluid phase, with an evolving interface. The cylinders are organized on a hexagonal lattice, see Fig. 1. If all the lattice sites are filled we reach the maximum solid volume fraction 111Actually we reach a packing fraction slightly lower than the maximum possible one by reducing the diameter of the cylinders such that they are initially placed in a manner that the surfaces of neighboring cylinders do not touch – i.e. there is a small gap between their surfaces.. In the rest of this paper we use a porosity $\Phi$, which is the fraction of the total volume occupied by the fluid, equal to $0.42$ by removing a certain number of randomly selected cylinders. The cylinders are made of a hyper-elastic Mooney-Rivlin Bonet and Wood (1997) material characterized by a shear elastic modulus, $G$. We emphasize that the full non-linear structural response of the elastic solid is included in the simulations. The theoretical model considered later in this paper, however, is simpler. The motion of the fluid and of the viscoelastic material are governed by the conservation of momentum and the incompressibility constraint :

[TABLE]

where the Cauchy stress tensor $\sigma_{ij}$ has contributions from both solid and fluid stresses with a weight set by the phase variable $\phi$, i.e.,

[TABLE]

with suffixes f,s used to distinguish the two phases, fluid f and elastic solid s. The fluid is assumed to be Newtonian and the solid is an incompressible viscous hyper-elastic material with constitutive equations:

[TABLE]

Here $p$ is the pressure, $\rho$ and $\mu$ are respectively the density and the dynamic viscosity both of which are assumed to be the same in the two phases222Note that, deformation of the solid becomes stationary and the value of its viscosity does not matter anymore since the solid viscosity only contributes when the solid is undergoing deformation prior to reaching the steady state., $\mathcal{D}_{ij}$ the rate-of-strain tensor and $\delta_{ij}$ is the Kronecker delta. The last term in $\sigma_{ij}^{\rm s}$, Eq. (4b), is the hyper-elastic contribution $\tau^{\rm e}_{ij}$, here modeled as a neo-Hookean Mooney-Rivlin material with the constitutive relation $\tau_{ij}^{\rm e}=G\mathcal{B}_{ij}$, where $\mathcal{B}_{ij}$ is the left Cauchy-Green tensor sometimes also called the Finger tensor. The full set of equations can be closed in a purely Eulerian manner by updating $\mathcal{B}_{ij}$ and $\phi$ with the following transport equations (Astarita and Marrucci, 1974; Larson, 1988; Bonet and Wood, 1997; Sugiyama et al., 2011)

[TABLE]

The algorithms have been described, used and validated against standard test cases in several earlier publications Rosti and Brandt (2017); Rosti et al. (2018); Rosti and Brandt (2018), and more details can be found in these references.

We use a rectangular domain of size $MD\times ND$, where $D$ is the diameter of an undeformed cylinder, for three different sets of values for $M$ and $N$. We apply periodic boundary conditions in the stream-wise $x$-direction, and no-slip/no-penetration boundary conditions on the two rigid walls bounding the domain in the $y$-direction. We consider $6$ different values for the shear elastic modulus $G$, and for each one of them we impose $6$ different pressure differences from $0.5$ to $50$ to drive the flow and measure the resulting flux, resulting in a Reynolds number varying in the range $Re=\rho qD/\mu\in\left[10^{-5}:10^{-4}\right]$ (in particular, we used $\rho=1$, $D=0.22$ and $\mu=1$ in our simulations). In all the cases, the numerical domain is discretised with $68$ grid points per diameter $D$. After a short transient time, the flow and the deformation of the solid skeleton reach a stationary state (see Fig. 1). As $G$ decreases – the solid skeleton is more deformable – the flow changes the width and nature of the channels through which the liquid flows. We observe a general tendency of the flow to exploit the defects of the underlying lattice, by generating preferential channels which transport most of the fluid.

III Numerical Results

In Fig. 2 we show how the Darcy flux, $q$, depends on the pressure-difference, $\Delta P$, across the domain: the most rigid case (black line) shows a linear increase of the flow rate with the pressure difference – the standard Darcy’s law for rigid porous materials. By contrast, as the material becomes more elastic, we observe a non-linear growth of the Darcy flux steeper than for its rigid counterpart, i.e., a super-linear dependence of the Darcy flux on the pressure-difference. A different way to interpret the same result is to say that the permeability of the porous medium, $k$, is itself a nonlinear function of the pressure difference (or the flow rate).

We check how robust this result is in the following ways. We run simulations in three domain sizes, the smallest one being, $M=8$ and $N=9$, in the next one we double the size in each direction, $M=16$ and $N=18$ and then obtain the largest one by again doubling the size, $M=32$ and $N=36$. We obtain the same result in all the three cases. Next, in the largest domain, we select sub-domains of the same size as the smallest domain and the relationship shown in Fig. 2 remains the same for these sub-domains too. Note that, we can reach the same porosity in different ways, depending on the position of the defects in the material; we have checked that our results remain unchanged in two different realisations of the random defects.

IV Theoretical Model

Next we construct a theory for this behavior. Our specific numerical simulations act as a motivation for the theory but the theory is not necessarily limited to our numerical model. Typically, the theory of porous media involves multiple scales, Whitaker (1986a, b); Mei and Auriault (1989, 1991); Auriault and Boutin (1992); MacMinn et al. (2016); Collis et al. (2017), ours is no exception.

In Fig. 1 we show a two-dimensional cross-section of our domain on three different spatial scales, from left to right we go from the full domain (size $L$) down to the scale of sub-domains ($\mathcal{L}$), smaller than the full domain but still much larger than the size of a single deformable circle, down to the scale of a single channel ($\ell$) between the elastic cylinders, such that $\ell\ll\mathcal{L}\ll L$. Our first step is to derive a relation between the flux and pressure difference across a deformable two-dimensional channel, e.g., we solve the microscale problem sketched in the rightmost panel of Fig. 1.

Let $\ell$ be the length of this channel, ${\Delta p}$ the pressure-difference across the channel and $h(\xi)$ the width of the channel as a function of the stream-wise coordinate, $\xi$. Within the range of parameters used in our simulations, we can safely assume that within this channel the Reynolds number is so low that the flow can be described by the Stokes equations. In fact, we shall go one step further and assume that it is safe to use the lubrication approximation Batchelor (1967). In particular, we assume that within a channel the pressure is a function of the stream-wise direction alone, i.e., $p=p(\xi)$, the wall-normal component of the velocity is zero, and the velocity gradient along the stream-wise direction is much smaller than that in the wall-normal direction. In addition, the flow velocity must go to zero at the boundaries of the channel, hence the flow-rate through the channel is given by

[TABLE]

The flow-rate must be a constant, independent of $\xi$. The width of the channel, $h(\xi)$, is determined by the mutual interaction between the flow and the elastic property of the walls of the channel. To make further progress, we assume a Hookean response of the boundary of the channel, often called Winkler foundation Wang et al. (2005); Dillard et al. (2018) in other context. In this framework, the undeformed width of the channel is $h_{\rm 0}$, which together with the deformation $w(\xi)$ sets the total width of the channel, $h(\xi)=h_{\rm 0}+w(\xi)$; the elastic property of the channel walls is parameterized by a Hookean spring with a spring constant $\kappa$, such that the force-per-unit-area necessary to generate a deformation $w$ is given by $\kappa w$. Hence the pressure and the deformation are related by $w(\xi)=p(\xi)/\kappa$. This allows us to write a differential equation for $p(\xi)$ where the flow-rate, $q$, appears as a parameter – the equation has the general form of $q=-\sigma(p)(dp/d\xi)$ Rubinow and Keller (1972). We integrate it and enforce the result to conform to the form of Darcy’s law, $q=(\mathcal{K}/\mu)({\Delta p}/\ell)$ with

[TABLE]

Here, $\beta\equiv{\Delta p}/(\kappa h_{\rm 0})$ and we have used the boundary conditions $p(0)={\Delta p}$ and $p(\ell)=0$. To build a connection to our simulations it is appropriate to choose $\kappa$ such that $\kappa h_{\rm 0}=G$. The use of the lubrication approximation coupled with the elastic properties of solid is quite commonly used to analyze flows in deformable channels, see e.g., Davis et al. (1986); Grotberg and Jensen (2004); Gomez et al. (2017); Christov et al. (2018). Christov et al. (2018) contains derivation of a similar relationship using a systematic application of asymptotics for a three dimensional channel. The only difference is that in Christov et al. (2018) a different model for the elastic wall – isotropic quasi-static bending of a plate under a transverse load due to the fluid pressure – is used.

In the next step we consider the mesoscopic scale, larger than the size of single cylinders but still smaller than the scale of the whole bed, see the middle panel in Fig. 1. The large domain contains many such mesoscopic domains of the same size. In the event of no defects, each of these subdomain contains $m\times n$ cylinders organized on a regular hexagonal lattice and the channels between the cylinders form a regular honeycomb lattice. In this case, each subdomain has exactly the same porosity and the same permeability.

Recall that we have randomly removed few cylinders from a regular hexagonal lattice to create our porous medium. Thus, the porosity at the scale $\mathcal{L}$, $\Phi_{\mathcal{L}}$, is different in different sub-domains. In the inset of Fig. 3 we show a representative plot of $\Phi_{\mathcal{L}}$, the porosity averaged over a domain of size $\mathcal{L}$, extracted from the DNS with the largest domain. In particular, we perform a volume average of the local fluid fraction $\phi$ on a domain of size $\mathcal{L}$ resulting in the porosity $\Phi_{\mathcal{L}}$. We incorporate this randomness into our model by choosing different values of $h_{\rm 0}$ – undeformed width of the channel – in different sub-domains. Thus at this mesoscale our model is a honeycomb network of channels. The length of the channels is same in all subdomains and within each subdomain all the channels have the same width set by the subdomain porosity $\Phi_{\mathcal{L}}$. Consequently, the width of the channels is different in different subdomains.

For a single subdomain, our task is to calculate the effective permeability of a network of channels, $\mathcal{K}_{\mathcal{L}}$, where the flow rate in each channel is given as function of the pressure drop by $q=[\mathcal{K}({\Delta p})/\mu]({\Delta p}/\ell)$, where the permeability of each channel $\mathcal{K}({\Delta p})$ is the nonlinear function, $f(\beta)$ in Eq. (9). Given a nonlinear function $f(\beta)$ there is no general method of attack known to us. We proceed therefore by assuming that $\mathcal{K}$ is independent of ${\Delta p}$. In this case, the problem corresponds to that of the effective conductivity of a honeycomb network of resistors by mapping the $q$ to current, ${\Delta p}/\ell$ to the voltage drops across the bonds of the network, and $\mathcal{K}/\mu$ to the conductivity of each of the bond, in bus-bar geometry Redner (2009) – the network is connected to two parallel lines and the battery is connected across the two lines. We solve this linear problem by matrix inversion to obtain $\mathcal{K}_{\mathcal{L}}=\gamma\mathcal{K}(\mathcal{L}/W)$ where $W$ is the width of the sub-domain and $\gamma$ is a constant that depends on $m$ and $n$. In particular, we use $m=8,\;n=9$ and obtain $\gamma=3.047$.

Thus, we obtain $\mathcal{K}_{\mathcal{L}}=CD^{2}\Phi_{\mathcal{L}}^{3}f(\beta)$, with a constant $C=\gamma(3/8)^{3}(\mathcal{L}/W)$. Here we have used $\mathcal{L}=8D$, and $\mathcal{L}/W=1.125$. The form of the function in Eq. (9) suggests that the effective permeability, $\mathcal{K}_{\mathcal{L}}$, is solely a function of the dimensionless parameter, $\beta=\Delta P/G$, where $\Delta P$ is the pressure-difference across the sub-domain. In Fig. 3 we show a scatter plot of $\mathcal{K}_{\mathcal{L}}$ versus $\Phi_{\mathcal{L}}$ for different sub-domains and observe that $\mathcal{K}_{\mathcal{L}}$ grows as $\Phi_{\mathcal{L}}^{3}$, although with some scatter of the data. Notwithstanding this, when we plot the values of $\mathcal{K}_{\mathcal{L}}$ normalized by $CD^{2}\Phi_{\mathcal{L}}^{3}$ as a function of $\beta$ for the different sub-domains, we obtain a reasonable data-collapse in agreement with our theory.

There is another, equivalent, way to calculate the effective permeability of a mesoscale subdomain, which also requires the assumption of linearity. In particular, we incorporate the random removal of cylinders by mapping to the problem of a honeycomb network of resistors such that each bond in the network has a conductance of $\mathcal{K}$ with probability $\mathcal{P}$ or infinite conductance (zero resistance) with probability $1-\mathcal{P}$. This probability, $\mathcal{P}$, is different in each sub-domain. Using the expression for effective resistance of infinite but random lattices Redner (2009), we obtain the effective permeability to be $\mathcal{K}_{\mathcal{L}}=\frac{\alpha}{\left(\alpha+1\right)\mathcal{P}-1}\mathcal{K}$, where $\alpha$ is a geometric parameters that depends on the lattice.

In the last and final step, we average over different sub-domains to obtain an effective permeability for the whole domain. Taking the divergence of the Darcy’s flux in a sub-domain, we obtain

[TABLE]

which is a steady-state heat equation with variable diffusivity $\mathcal{K}_{\mathcal{L}}(\Phi_{\mathcal{L}})$, function of a fast variable $\Phi_{\mathcal{L}}$. Straightforward application of the method of multiple scales (see, e.g., Frisch, 1996, section 9.6.2) shows that the effective permeability $k$ is the harmonic mean of the effective permeability of each sub-domain:

[TABLE]

where $P(\Phi_{\mathcal{L}})$ is the probability density function of the porosity of a sub-domain, $\Phi_{\mathcal{L}}$. The data in Fig. 3 justify treating $P(\Phi_{\mathcal{L}})$ as a Gaussian with mean value equal to the mean porosity calculated over the full domain so that Eq. (11) is integrable. As the Gaussian is sharply peaked, the integral is well-approximated by its leading order contribution, i.e., $k\approx\mathcal{K}_{\mathcal{L}}(\Phi)$ where $\Phi=\left\langle\Phi_{\mathcal{L}}\right\rangle$ is the mean porosity of the whole domain. This implies that the collapse we have observed for each sub-domain should also work if we plot the permeability $k$ of the full domain as a function of $\beta=\Delta P/G$ where the pressure-difference, $\Delta P$, now is across the whole domain. This is confirmed by the results in the inset of Fig. 2.

V Discussion and Conclusions

Several comments are now in order. For the analytical calculations we have used a simple model, the Winkler foundation, whereas we have used the hyper-elastic Mooney-Rivlin model for the cylinders in our DNS. The qualitative agreement between the two shows the robustness of our results. Different elastic models will result in different expression for the function $f(\beta)$ in Eq. (9). A crucial result of our work is to show that such a function exists, which implies that data on permeability collapse to a single function when plotted as a function of $\beta$.

The enhanced flux, which is the most striking result of our work, has not been observed in experiments but is found analytically for certain classes of models. In particular, two among the five models discussed in MacMinn et al. (2016) – these models start from the intermediate scale denoted as $\mathcal{L}$ here– show the possibility of super-linear response because of the nonlinear elastic behavior of the solid skeleton. These two models further assume that $\mathcal{K}_{\mathcal{L}}$ is independent of $\Phi_{\mathcal{L}}$.

In experiments Hewitt et al. (2016); Song et al. (2019), when a fluid is forced through a deformable porous medium, the boundary between the porous material and the fluid on the inlet is normally left unconstrained. Hence under fluid pressure the boundary moves and squeezes the porous material. This decreases the permeability – often modeled by the empirical Kozeny-Carman formula Carman (1937) – of the porous material. Two of the models in MacMinn et al. (2016), which include the Kozeny-Carman formula, show nonlinear but sub-linear behavior. Recent experimental measurements by Song et al. (2019) also exhibit such behavior. In our case, both theory and the DNS approximate the behavior of the Kozeny-Carman function for small $\Phi_{\mathcal{L}}$, $\mathcal{K}_{\mathcal{L}}\sim\Phi_{\mathcal{L}}^{3}$. Indeed, in our simulations the boundaries are held fixed, hence by construction fluid-driven compaction is missing from our simulations. Hence, we also expect that it is possible to observe the enhanced Darcy flux in experiments, but not for very large pressure-differences where fluid-driven compaction dominates. We hope our work will encourage further experimental and numerical explorations.

Most studies in this field using homogenization to understand the fluid flow through rigid/deformable/active porous media Whitaker (1986a, b); Mei and Auriault (1989, 1991); Auriault and Boutin (1992); Collis et al. (2017), adopt a continuum description for both the solid and fluid phase and couple them through the kinematic interface conditions. The constitutive equations for the pore-scale description of the problem are coarse-grained to obtain a description of the equivalent fluid-solid interaction. The transport coefficients of the resultant equations depend on the solvability conditions (closure problem) of the homogenization techniques. However, the closure problem remained unsolved in all these models, thus no explicit Darcy-like relation was obtained. Our theory stands apart from such models. We write down a Darcy-like relation between flow-rate and pressure-difference with an explicit expression for the permeability that depends on the shear modulus of the solid skeleton. Furthermore, we show that the nonlinear flow-rate versus pressure-difference relation is an intrinsic property of the medium rooted in the pore scale rearrangement induced by fluid flow. The weakest link in our theory is the assumption of linearity to calculate the effective permeability of a network of channels.

Acknowledgements

Our work is inspired by experiments being done by John Wettlaufer and his group. We thank John for helping us at every stage of this work. We thank Dominic Vella for introducing us to the Winkler foundations. MR and LB are supported by the European Research Council Grant no. ERC-2013-CoG-616186, TRITOS and by the Swedish Research Council Grant no. VR 2014-5001. SP acknowledges the support of the Swedish Research Council Grant no. 638-2013-9243. DM acknowledges the support of the Swedish Research Council Grant no. 638-2013-9243 as well as 2016-05225. We gratefully acknowledge computer time provided by SNIC (Swedish National Infrastructure for Computing).