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

Abstract

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 law¹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,⁷ 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.¹⁰

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¹¹ 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–Rivlin¹² 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: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 = ϕσ^f_ij + (1 − ϕ)σ^s_ij, (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:

σ^f_ij = −pδ_ij + 2μ[scr D, script letter D]_ij, (4a)

σ^s_ij = −pδ_ij + 2μ[scr D, script letter D]_ij + τ^e_ij. (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 σ^s_ij, eqn (4b), is the hyper-elastic contribution τ^e_ij, here modeled as a neo-Hookean Mooney–Rivlin material with the constitutive relation τ^e_ij = 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 equations^12–15The 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.

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⁻⁵ : 10⁻⁴] (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).

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 = CD²Φ³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.¹ 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³/(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 foundation^25,26 in other context. In this framework, the undeformed width of the channel is h₀, which together with the deformation w(ξ) sets the total width of the channel, h(ξ) = h₀ + 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ξ).²⁷ 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] = h₀³f(β), (8)

whereHere, β ≡ Δp/(κh₀) 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 κh₀ = 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.,²⁸ Grotberg and Jensen,²⁹ Gomez et al.,³⁰ Christov et al.³¹ Christov et al.³¹ 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.³¹ 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 h₀ – 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.

Fig. 3 A scatter plot of the local permeability [scr K, script letter K]_{[script L]} in each sub-domain normalized with CD² 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 geometry³² – 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, n = 9 and obtain γ = 3.047.

Thus, we obtain [scr K, script letter K]_{[script L]} = CD²Φ_{[script L]}³f(β), with a constant C = γ(3/8)³([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]}³, although with some scatter of the data. Notwithstanding this, when we plot the values of [scr K, script letter K]_{[script L]} normalized by CD²Φ_{[script L]}³ 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,³² 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

∇·[[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: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.¹¹ – 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 formula³⁶ – of the porous material. Two of the models in Macminn et al.,¹¹ which include the Kozeny–Carman formula, show nonlinear but sub-linear behavior. Recent experimental measurements by Song et al.³⁵ 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]}³. 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.

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

References

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

---

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.

---