Open Access Article

This Open Access Article is licensed under a

Creative Commons Attribution 3.0 Unported Licence

DOI: 10.1039/C9SM01678C
(Paper)
Soft Matter, 2020, Advance Article

Marco Edoardo Rosti*^{a},
Satyajit Pramanik^{b},
Luca Brandt^{a} and
Dhrubaditya Mitra^{b}
^{a}Linné Flow Centre and SeRC (Swedish e-Science Research Centre), KTH Mechanics, SE 100 44 Stockholm, Sweden. E-mail: merosti@mech.kth.se
^{b}Nordita, Royal Institute of Technology and Stockholm University, SE 106 91 Stockholm, Sweden

Received
19th August 2019
, Accepted 9th December 2019

First published on 17th December 2019

We perform direct numerical simulations of the flow through a model of 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 (Δ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 material with a linear constitutive relation but can undergo large deformation. Our theory further predicts that the flow permeability is an universal function of ΔP/G, which is confirmed by the present simulations.

(1) |

Poroelasticity play an important role in understanding the transport through a wide range of materials ranging from individual cells,^{2,3} to biological tissues, e.g., soft-tissues,^{4–6} bones,^{7} even to hydraulic fracture.^{8,9} 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.^{10}

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^{11} but has never been observed before in simulations or experiments.

(2) |

σ_{ij} = ϕσ^{f}_{ij} + (1 − ϕ)σ^{s}_{ij},
| (3) |

σ^{f}_{ij} = −pδ_{ij} + 2μ_{ij},
| (4a) |

σ^{s}_{ij} = −pδ_{ij} + 2μ_{ij} + τ^{e}_{ij}.
| (4b) |

(5) |

(6) |

We use a rectangular domain of size MD × 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 = ρqD/μ ∈ [10^{−5}:10^{−4}] (in particular, we used ρ = 1, D = 0.22 and μ = 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.

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.

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 (), 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 () between the elastic cylinders, such that ≪ ≪ 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 be the length of this channel, Δp the pressure-difference across the channel and h(ξ) the width of the channel as a function of the stream-wise coordinate, ξ. 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.^{1} In particular, we assume that within a channel the pressure is a function of the stream-wise direction alone, i.e., p = p(ξ), 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

q = −[h^{3}/(12μ)](dp/dξ).
| (7) |

= h_{0}^{3}f(β),
| (8) |

(9) |

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 × 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, , Φ_{}, is different in different sub-domains. In the inset of Fig. 3 we show a representative plot of Φ_{}, the porosity averaged over a domain of size, , extracted from the DNS with the largest domain. In particular, we perform a volume average of the local fluid fraction ϕ on a domain of size resulting in the porosity Φ_{}. We incorporate this randomness into our model by choosing different values of h_{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 Φ_{}. 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, _{}, where the flow rate in each channel is given as function of the pressure drop by q = [(Δp)/μ](Δp/), where the permeability of each channel (Δp) is the nonlinear function, f(β) in eqn (9). Given a nonlinear function f(β) there is no general method of attack known to us. We proceed therefore by assuming that is independent of Δ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, Δp/ to the voltage drops across the bonds of the network, and /μ to the conductivity of each of the bond, in bus-bar geometry^{32} – 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 _{} = γ(/W) where W is the width of the sub-domain and γ is a constant that depends on m and n. In particular, we use m = 8,n = 9 and obtain γ = 3.047.

Thus, we obtain _{} = CD^{2}Φ_{}^{3}f(β), with a constant C = γ(3/8)^{3}(/W). Here we have used = 8D, and /W = 1.125. The form of the function in eqn (9) suggests that the effective permeability, _{}, is solely a function of the dimensionless parameter, β = ΔP/G, where ΔP is the pressure-difference across the sub-domain. In Fig. 3 we show a scatter plot of _{} versus Φ_{} for different sub-domains and observe that _{} grows as Φ_{}^{3}, although with some scatter of the data. Notwithstanding this, when we plot the values of _{} normalized by CD^{2}Φ_{}^{3} as a function of β 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 with probability or infinite conductance (zero resistance) with probability 1 − This probability, , is different in each sub-domain. Using the expression for effective resistance of infinite but random lattices,^{32} we obtain the effective permeability to be , where α 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

∇·[_{}(Φ_{})∇p] = 0,
| (10) |

(11) |

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.^{11} – these models start from the intermediate scale denoted as here – show the possibility of super-linear response because of the nonlinear elastic behavior of the solid skeleton. These two models further assume that _{} is independent of Φ_{}.

In experiments,^{34,35} 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^{36} – of the porous material. Two of the models in Macminn et al.,^{11} which include the Kozeny–Carman formula, show nonlinear but sub-linear behavior. Recent experimental measurements by Song et al.^{35} also exhibit such behavior. In our case, both theory and the DNS approximate the behavior of the Kozeny–Carman function for small Φ_{}, _{} ∼ Φ_{}^{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,^{19–24} 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.

- G. Batchelor, An Introduction to Fluid Dynamics, Cambridge University Press, Cambridge, UK, 1967 Search PubMed.
- G. T. Charras, J. C. Yarrow, M. A. Horton, L. Mahadevan and T. Mitchison, Nature, 2005, 435, 365 CrossRef CAS PubMed.
- E. Moeendarbary, L. Valon, M. Fritzsche, A. R. Harris, D. A. Moulding, A. J. Thrasher, E. Stride, L. Mahadevan and G. T. Charras, Nat. Mater., 2013, 12, 253 CrossRef CAS PubMed.
- W. M. Lai, J. Hou and V. C. Mow, J. Biomech. Eng., 1991, 113, 245–258 CrossRef CAS.
- M. Yang and L. A. Taber, J. Biomech., 1991, 24, 587–597 CrossRef CAS PubMed.
- L. C. Auton and C. W. MacMinn, Proc. R. Soc. A, 2017, 473, 20160753 CrossRef CAS PubMed.
- S. C. Cowin, J. Biomech., 1999, 32, 217–238 CrossRef CAS PubMed.
- E. Detournay and A.-D. Cheng, International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1988, pp. 171–182.
- V. M. Yarushina, D. Bercovici and M. L. Oristaglio, Geophys. J. Int., 2013, 194, 1514–1526 CrossRef.
- C. W. MacMinn, E. R. Dufresne and J. S. Wettlaufer, Phys. Rev. X, 2015, 5, 011020 Search PubMed.
- C. W. MacMinn, E. R. Dufresne and J. S. Wettlaufer, Phys. Rev. Appl., 2016, 5, 044020 CrossRef.
- J. Bonet and R. D. Wood, Nonlinear continuum mechanics for finite element analysis, Cambridge University Press, 1997 Search PubMed.
- G. Astarita and G. Marrucci, Principles of non-Newtonian fluid mechanics, McGraw-Hill, 1974 Search PubMed.
- R. G. Larson, Constitutive equations for polymer melts and solutions, Elsevier, 1988 Search PubMed.
- K. Sugiyama, S. Ii, S. Takeuchi, S. Takagi and Y. Matsumoto, J. Comput. Phys., 2011, 230, 596–627 CrossRef CAS.
- M. E. Rosti and L. Brandt, J. Fluid Mech., 2017, 830, 708–735 CrossRef CAS.
- M. E. Rosti, L. Brandt and D. Mitra, Phys. Rev. Fluids, 2018, 3, 012301(R) CrossRef.
- M. E. Rosti and L. Brandt, J. Non-Newtonian Fluid Mech., 2018, 262, 3–11 CrossRef CAS.
- S. Whitaker, Transp. Porous Media, 1986, 1, 3–25 CrossRef.
- S. Whitaker, Transp. Porous Media, 1986, 1, 127–154 CrossRef CAS.
- C. C. Mei and J. L. Auriault, Proc. R. Soc. A, 1989, 426, 391–423 CrossRef.
- C. C. Mei and J.-L. Auriault, J. Fluid Mech., 1991, 222, 647–663 CrossRef.
- J. L. Auriault and C. Boutin, Transp. Porous Media, 1992, 7, 63–82 CrossRef CAS.
- J. Collis, D. L. Brown, M. E. Hubbard and R. D. O'Dea, Proc. R. Soc. A, 2017, 473, 20160755 CrossRef CAS PubMed.
- Y. Wang, L. Tham and Y. Cheung, Progr. Struct. Eng. Mater., 2005, 7, 174–182 CrossRef.
- D. A. Dillard, B. Mukherjee, P. Karnal, R. C. Batra and J. Frechette, Soft Matter, 2018, 14, 3669–3683 RSC.
- S. Rubinow and J. B. Keller, J. Theor. Biol., 1972, 35, 299–313 CrossRef CAS PubMed.
- R. H. Davis, J.-M. Serayssol and E. J. Hinch, J. Fluid Mech., 1986, 163, 479–497 CrossRef.
- J. B. Grotberg and O. E. Jensen, Annu. Rev. Fluid Mech., 2004, 36, 121–147 CrossRef.
- M. Gomez, D. E. Moulton and D. Vella, Phys. Rev. Lett., 2017, 119, 144502 CrossRef PubMed.
- I. C. Christov, V. Cognet, T. C. Shidhore and H. A. Stone, J. Fluid Mech., 2018, 841, 267–286 CrossRef CAS.
- S. Redner, Encyclopedia of complexity and systems science, Springer, 2009, pp. 3737–3754 Search PubMed.
- U. Frisch, Turbulence the legacy of A.N. Kolmogorov, Cambridge University Press, Cambridge, 1996 Search PubMed.
- D. R. Hewitt, J. S. Nijjer, M. G. Worster and J. A. Neufeld, Phys. Rev. E, 2016, 93, 023116 CrossRef PubMed.
- R. Song, H. A. Stone, K. H. Jensen and J. Lee, J. Fluid Mech., 2019, 871, 742–754 CrossRef CAS.
- P. Carman, Chem. Eng. Res. Des., 1937, 15, S32–S48 Search PubMed.

## Footnotes |

† Actually 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. |

‡ Note 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. |

This journal is © The Royal Society of Chemistry 2020 |