Open Access Article
This Open Access Article is licensed under a
Creative Commons Attribution 3.0 Unported Licence

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

Marco Edoardo Rosti*a, Satyajit Pramanikb, Luca Brandta and Dhrubaditya Mitrab
aLinné Flow Centre and SeRC (Swedish e-Science Research Centre), KTH Mechanics, SE 100 44 Stockholm, Sweden. E-mail:
bNordita, 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 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, ΔP. In the simplest case of a rigid isotropic medium, this gives rise to Darcy's law1
image file: c9sm01678c-t1.tif(1)
where μ 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,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 level11 but has never been observed before in simulations or experiments.

2 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 ϕ = 0 inside the viscoelastic solid phase and ϕ = 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. In the rest of this paper we use a porosity Φ, 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–Rivlin12 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:
image file: c9sm01678c-t2.tif(2)
where the Cauchy stress tensor σij has contributions from both solid and fluid stresses with a weight set by the phase variable ϕ, i.e.,
σij = ϕσfij + (1 − ϕ)σsij, (3)
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:
σfij = −ij + 2μ[scr D, script letter D]ij, (4a)
σsij = −ij + 2μ[scr D, script letter D]ij + τeij. (4b)
Here p is the pressure, ρ and μ are respectively the density and the dynamic viscosity both of which are assumed to be the same in the two phases, [scr D, script letter D]ij the rate-of-strain tensor and δij is the Kronecker delta. The last term in σsij, eqn (4b), is the hyper-elastic contribution τeij, here modeled as a neo-Hookean Mooney–Rivlin material with the constitutive relation τeij = G[scr B, script letter B]ij, where [scr B, script letter 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 [scr B, script letter B]ij and ϕ with the following transport equations12–15
image file: c9sm01678c-t3.tif(5)
image file: c9sm01678c-t4.tif(6)
The algorithms have been described, used and validated against standard test cases in several earlier publications,16–18 and more details can be found in these references.

image file: c9sm01678c-f1.tif
Fig. 1 From left to right, a snapshot of the cross-section of one of our simulation under different levels of magnification. We magnify twice into a part of the domain to show first a sub-domain and then a channel. We apply the lubrication theory to study the flow thorough this channel, which implies that the flow inside the channel is assumed to be parabolic.

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[thin space (1/6-em)]:[thin space (1/6-em)]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.

3 Numerical results

In Fig. 2 we show how the Darcy flux, q, depends on the pressure-difference, Δ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).
image file: c9sm01678c-f2.tif
Fig. 2 The Darcy flux q as a function of the overall pressure difference ΔP/L for different values of the deformability of the medium, G ∈ [0.5,50]. In particular, the red, brown, orange, green, blue and black colors are used for G = 0.5, 1, 1.5, 3, 6 and 50, respectively. We also plot the results for different domain sizes (grey squares and cross). In the inset we plot the permeability k normalized with [scr K, script letter K]t = CD2Φ3 versus the dimensionless variable β = ΔP/G. C is a constant coming from the theory, equal to 0.1808 for the present cases.

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.

4 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,11,19–24 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 ([script 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 ([small script l]) between the elastic cylinders, such that [small script l][script L]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 [small script l] 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 = −[h3/(12μ)](dp/dξ). (7)
The flow-rate must be a constant, independent of ξ. The width of the channel, h(ξ), 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 foundation25,26 in other context. In this framework, the undeformed width of the channel is h0, which together with the deformation w(ξ) sets the total width of the channel, h(ξ) = h0 + w(ξ); the elastic property of the channel walls is parameterized by a Hookean spring with a spring constant κ, such that the force-per-unit-area necessary to generate a deformation w is given by κw. Hence the pressure and the deformation are related by w(ξ) = p(ξ)/κ. This allows us to write a differential equation for p(ξ) where the flow-rate, q, appears as a parameter – the equation has the general form of q = −σ(p)(dp/dξ).27 We integrate it and enforce the result to conform to the form of Darcy's law, q = ([scr K, script letter K]/μ)(Δp/[small script l]) with
[scr K, script letter K] = h03f(β), (8)
image file: c9sm01678c-t5.tif(9)
Here, β ≡ Δp/(κh0) and we have used the boundary conditions p(0) = Δp and p([small script l]) = 0. To build a connection to our simulations it is appropriate to choose κ such that κh0 = 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.,28 Grotberg and Jensen,29 Gomez et al.,30 Christov et al.31 Christov et al.31 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.31 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 × 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, [script L], Φ[script L], is different in different sub-domains. In the inset of Fig. 3 we show a representative plot of Φ[script L], the porosity averaged over a domain of size, [script L], 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 [script L] resulting in the porosity Φ[script L]. We incorporate this randomness into our model by choosing different values of h0 – 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 Φ[script L]. Consequently, the width of the channels is different in different subdomains.

image file: c9sm01678c-f3.tif
Fig. 3 A scatter plot of the local permeability [scr K, script letter K][script L] in each sub-domain normalized with CD2 as a function of the corresponding local porosity Φ[script L]. The inset shows a snapshot of porosity (averaged over the scale [script L]), Φ[script L]. The color scale goes from 0.3 (black) to 0.6 (white).

For a single subdomain, our task is to calculate the effective permeability of a network of channels, [scr K, script letter K][script L], where the flow rate in each channel is given as function of the pressure drop by q = [[scr K, script letter K]p)/μ](Δp/[small script l]), where the permeability of each channel [scr K, script letter K]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 [scr K, script letter K] 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/[small script l] to the voltage drops across the bonds of the network, and [scr K, script letter K]/μ to the conductivity of each of the bond, in bus-bar geometry32 – 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 [scr K, script letter K][script L] = γ[scr K, script letter K]([script L]/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,[thin space (1/6-em)]n = 9 and obtain γ = 3.047.

Thus, we obtain [scr K, script letter K][script L] = CD2Φ[script L]3f(β), with a constant C = γ(3/8)3([script L]/W). Here we have used [script L] = 8D, and [script L]/W = 1.125. The form of the function in eqn (9) suggests that the effective permeability, [scr K, script letter K][script L], 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 [scr K, script letter K][script L] versus Φ[script L] for different sub-domains and observe that [scr K, script letter K][script L] grows as Φ[script L]3, although with some scatter of the data. Notwithstanding this, when we plot the values of [scr K, script letter K][script L] normalized by CD2Φ[script L]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 [scr K, script letter K] with probability [scr P, script letter P] or infinite conductance (zero resistance) with probability 1 − [scr P, script letter P] This probability, [scr P, script letter P], 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 image file: c9sm01678c-t6.tif, 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

∇·[[scr K, script letter K][script L](Φ[script L])∇p] = 0, (10)
which is a steady-state heat equation with variable diffusivity [scr K, script letter K][script L](Φ[script L]), function of a fast variable Φ[script L]. Straightforward application of the method of multiple scales see, e.g., ref. 33, Section 9.6.2 shows that the effective permeability k is the harmonic mean of the effective permeability of each sub-domain:
image file: c9sm01678c-t7.tif(11)
where P(Φ[script L]) is the probability density function of the porosity of a sub-domain, Φ[script L]. The data in Fig. 3 justify treating P(Φ[script L]) as a Gaussian with mean value equal to the mean porosity calculated over the full domain so that eqn (11) is integrable. As the Gaussian is sharply peaked, the integral is well-approximated by its leading order contribution, i.e., k[scr K, script letter K][script L](Φ) where Φ = 〈Φ[script L]〉 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 β = ΔP/G where the pressure-difference, ΔP, now is across the whole domain. This is confirmed by the results in the inset of Fig. 2.

5 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(β) in eqn (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 β.

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 [script 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 [scr K, script letter K][script L] is independent of Φ[script L].

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 formula36 – 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 Φ[script L], [scr K, script letter K][script L]Φ[script 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,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.

Conflicts of interest

There are no conflicts to declare.


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).


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


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