Estimating the interfacial permeability for flow into a poroelastic medium

Abstract

Boundary conditions between a porous solid and a fluid has been a long-standing problem in modeling porous media. For deformable poroelastic materials such as hydrogels, the question is further complicated by the elastic stress from the solid network. Recently, an interfacial permeability condition has been developed from the principle of positive energy dissipation on the hydrogel–fluid interface. Although this boundary condition has been used in flow computations and yielded reasonable predictions, it contains an interfacial permeability η as a phenomenological parameter. In this work, we use pore-scale models of flow into a periodic array of solid cylinders or parallel holes to determine η as a function of the pore size and porosity. This provides a means to evaluate the interfacial permeability for a wide range of poroelastic materials, including hydrogels, foams and biological tissues, to enable realistic flow simulations.

---

Introduction

Flow through a deformable porous medium can be described by a poroelasticity model. This is essentially a mixture model that views the solid and the fluid as interpenetrating continua, each phase being described by its own volume fraction, velocity and stress that obey continuity and momentum equations for each phase.^1–3 Since the model erases all pore-scale geometric information about the fluid-solid boundaries, and instead poses momentum equations for both phases over the entire porous medium, additional boundary conditions (BCs) must be supplied to make the problem solvable.^4,5 This issue arises in soft porous media such as hydrogels, foams, flocculated fiber suspensions and biological tissues. For brevity, we will refer to our media only as “hydrogels” or “gels” hereafter. Even for rigid porous media, where the solid stress becomes indeterminate, a similar question arises, albeit in a different form.^6,7

The mathematical structure of the problem can be appreciated from a simplified one-dimensional (1D) flow (Fig. 1). The general vectorial forms of the governing equations and BCs can be found in Xu et al.⁸ A uniform flow of a solvent of viscosity μ encounters a layer of hydrogel whose downstream surface is fixed in space. On the upstream interface, the fluid experiences a contraction flow into the pores. Inside the gel layer, the fluid exerts a Darcy drag on the solid, and compresses the gel layer into a steady state. The poroelastic model for this steady state can be reduced to the following 1D form:⁸

(ϕ_fv_f)′ = 0, (1)

p′ = (ϕ_fσ_f + ϕ_sσ_s)′, (3)

where the prime indicates the spatial derivative d/dx, ϕ_f and ϕ_s = 1 − ϕ_f are the fluid and solid volume fractions, p is the pressure inside the gel, shared by the fluid and the solid phases, v_f is the fluid pore velocity, and σ_f and σ_s are respectively the fluid and solid normal stresses in the x direction. Eqn (1) is the fluid continuity equation. Eqn (2) is the Brinkman extension of Darcy's law for the fluid flow. The permeability k is a function of the local porosity, a commonly used form being the Kozeny–Carman model.^3,9 Finally, eqn (3) is the total momentum balance on the two phases. For the viscous solvent, . For the solid, σ_s is given by a linearly elastic or hyperelastic constitutive equation.⁸

Fig. 1 A one-dimensional (1D) flow into a hydrogel for illustrating the interfacial permeability boundary condition of eqn (6).

On the upstream interface, volume conservation and traction balance give us two BCs:

ϕ_fv_f = V₀, (4)

p − ϕ_fσ_f − ϕ_sσ_s = P, (5)

where V₀ and P are the uniform velocity and pressure upstream of the interface. For a hydrogel, an additional BC is required, and Young et al.^8,10,11 have derived the following based on the principle of positive energy dissipation on the fluid–gel interface:

V₀ − v_f = η(P − p + σ_s), (6)

with η > 0 being the interfacial permeability. At the downstream interface, we impose zero solid displacement (u_s = 0) and zero fluid stress. The permeability BC of eqn (6) is intuitive with the normal velocity jump proportional to the normal stress jump across the interface. It has also produced reasonable results in flow computations through hydrogels.^5,8,12 Nevertheless, η remains a phenomenological coefficient that has never been measured experimentally. For realistic flow simulations, it will be desirable to know its value for specific poroelastic media.

Pore-scale models

To estimate η, we adopt an idea from pore-scale models for rigid porous media. Many authors have calculated the bulk permeability of rigid porous media by the tube-bundle model^9,13,14 or by flow around periodic arrays of solid cylinders or spheres.^9,15–18 The solid-array representation has also been used to explore the interfacial conditions.^{6,7,16,17,19–21} We can use such a rigid pore model to estimate η for our poroelastic medium insofar as we seek a constitutive relation for η as a function of the steady-state local porosity ϕ_f and pore size r. The compression that has produced such a steady state and the fact that ϕ_f and r may vary in the bulk of the medium are of no immediate import to our purpose.

Another issue is that the BC of eqn (6) involves the solid stress σ_s, which is undefined for rigid solids. We bypass this obstacle by eliminating σ_s using the total traction balance of eqn (5). Furthermore, we eliminate v_f in favor of V₀ by the volume conservation of eqn (4) to arrive at the following interfacial permeability BC:

where we have simplified the notation by introducing the penetration length [small eta, Greek, tilde] = ημϕ_f²/ϕ_s². This removes the viscosity dependence and leaves [small eta, Greek, tilde] a geometric parameter, akin to the Navier slip length or the “transpiration length”.⁷ We have further dropped the viscous normal stress σ_f as it is typically much smaller than the friction that gives rise to the Darcy drag and the entry pressure drop P − p.^22,23 Now we can use a pore model to compute the entry pressure drop and back out the interfacial penetration length [small eta, Greek, tilde] and the interfacial permeability η.

It is interesting to note the similarity between the permeability BC of eqn (7) and the pressure jump BCs for rigid porous media.^6,7,20 For example, Lācis et al.⁷ proposed the following form (in our notation):

where the “resistance coefficient” f represents the friction force against the solvent passing through the interface. Although eqn (7) and (8) have come from different physical arguments, they are algebraically similar, and both recognize the need for an interfacial pressure drop to overcome the resistance to the fluid entering the pores. This idea deviates from the long-standing assumption of traction continuity or pressure continuity at the interface.^{21,22,24–26} Lācis et al.⁷ further used numerical computation in a periodic pore model to evaluate f, and we will be able to compare our result quantitatively with theirs.

We use four geometric representations for the porous material (Fig. 2a): semi-infinite arrays of circular and square cylinders aligned in the flow direction, and parallel circular and square holes in a semi-infinite solid. In either case, the cylinders or holes are arranged in a periodic square lattice on the cross-section of the medium. Thus, the solvent coming from the left of Fig. 1 experiences a contraction flow to enter the “pores”.

Fig. 2 (a) The interface for four pore geometries, with periodic arrays of circular and square cylinders and circular and square holes. The grey and white areas represent the solid and pore space, the latter occupying an area fraction that equals ϕ_f. The dashed squares, of edge length r, represent the cross-sections of our computational domain. (b) The computational domain for Stokes flow entering a pore between square cylinders, with its cross-section at the interface corresponding to the dashed square above. The flow is depicted by velocity vectors.

These geometries are inspired by prior models for porous medium, e.g. ref. 15–19, but with one important difference. Prior models specify a periodicity along the depth of the porous medium (i.e., the x-axis of Fig. 1). Then ambiguities arise as to where to place the interface, and how many periods to use for spatial averaging. Depending on whether the nominal interface cuts through the solid objects,^15,16 stays tangential to their apex in the top foremost row,^19,27 or is placed further in the bulk fluid,^7,17 the continuum-level result may differ greatly. Furthermore, James and Davis¹⁹ demonstrated that a shear flow penetrates only the first row of cylinders. This challenges the scale-separation assumption underlying the continuum Brinkman or Darcy model.⁷ Our setups avoid both difficulties; the semi-infinite cylinders or holes present a clear-cut interface, and they have “pre-averaged” the spatial variations in the depth direction. Experimentally, Tachie et al.^28,29 tested brush models that resemble our setup, although they only considered shear flows past the top of the transverse posts. Our results will turn out to be largely insensitive to the geometric details of the model (Fig. 3 below), and thus we will not consider other solid shapes or other spatial arrangements.^15,19

Fig. 3 The function ψ(ϕ_f) obtained by computing the Stokes flow into the model porous medium of Fig. 2 with different geometries. The solid curve represents the ψ₁(ϕ_f) function of eqn (12), while the dashed curve ψ₂(ϕ_f) of eqn (13).

Numerical and theoretical results

Taking advantage of spatial periodicity and symmetry, we solve the 3D Stokes flow in a computational domain that corresponds to a quarter of a single pore; Fig. 2(b) shows an example for the square cylinder. The 3D pore-scale flow will be processed to yield the interfacial penetration of the 1D setup on the continuum level (Fig. 1). We impose no-slip conditions on the solid surfaces and symmetry conditions on all fluid boundaries. The upstream boundary has a uniform incoming flow at velocity V₀ and the downstream boundary has stress-free conditions. The finite-element code is based on the open-source finite-element library deal.II. Mesh refinement has confirmed adequate numerical resolution. The up- and downstream boundaries are each 4r away from the opening of the pore so they have negligible effect on the contraction flow and the final result.

To compute the interfacial pressure drop P − p of eqn (7), we identify P with the upstream pressure P₀ in our pore model, and p with the averaged pore pressure at the opening of the pore: , where A_c is the cross-section of the interface excluding the solid portion (the shaded area of Fig. 2b). Dimensional analysis dictates

where ψ(ϕ_f) is a dimensionless function that we have determined numerically and plotted in Fig. 3. As expected, ψ decreases toward zero as ϕ_f → 0 and increases with ϕ_f. Furthermore, ψ differs little among the four geometries tested, suggesting that they are all robust and realistic representations of the interfacial property of a porous medium. On the continuum level, therefore, ϕ_f is the main geometric parameter that matters. The numerical data in Fig. 3 embody the idea that an interfacial pressure drop is needed to propel the fluid into the pores; it is essentially the excess pressure in entry flow.³⁰

To complement the numerical results, we can estimate theoretically the viscous dissipation of the entry flow, under the following two assumptions:

• The dissipation mostly occurs in an entry region Ω extending a distance of r upstream from the interface, based on the perturbed flow field of Fig. 2(b).

• The dissipation is mostly attributable to v_x, the x-component of the velocity.

The longitudinal velocity accelerates from V₀ at the upstream surface of Ω to V₀/ϕ_f at the interface, suggesting a strain rate ∂v_x/∂x ∼ (V₀/ϕ_f − V₀)/r = V₀ϕ_s/(rϕ_f). To account for the 3D nature of the flow, we introduce a fitting parameter C such that ∂v_x/∂x = CV₀ϕ_s/(rϕ_f). Similarly, we estimate the shear rates using the averaged velocity V₀/ϕ_f divided by the pore size r: ∂v_x/∂y = ∂v_x/∂z = CV₀/(ϕ_fr). Therefore, the total viscous dissipation inside Ω is

Meanwhile, the pressure work on the entrance of Ω is PV₀r², and that on the exit of Ω is −p_f(V₀/ϕ_f)(r²ϕ_f). Equating the total work (P − p_f)V₀r² to the viscous dissipation, we haveand in turn the following approximation to the ψ(ϕ_f) function of eqn (9): gives a good fit to the numerical data of Fig. 3. Considering the simplicity of the estimation, the semi-analytical ψ₁(ϕ_f) captures the numerical data well. The greatest discrepancy appears at the dilute limit ϕ_f → 1; some numerical data turn upward whereas the formula remains concave. This discrepancy can be understood as an effect of the pore geometry in the limit of ϕ_f → 1, in an exception to our earlier statement on the insensitivity to pore geometry. For flow along a cylinder, the cross-sectionally averaged shear stress vanishes as the cylinder becomes thinner. This explains the upturn in the numerical data for cylinders in Fig. 3. Conversely, for flow in a hole, the average shear stress stays finite when the wall become thinner. This difference is noticeable among the numerical data of Fig. 3. Our estimation of the shear rates ∂v_x/∂y = ∂v_x/∂z = CV₀/(ϕ_fr) remains finite as ϕ_f → 1. Thus, it reflects internal flow in holes, not external flow along cylinders.

As an alternative, we can estimate the average shear rate in the external axial flow along an infinite circular cylinder, with its analytical logarithmic velocity profile. Recalculating the viscous dissipation in Ω, we arrive at the following formula:

This ψ₂(ϕ_f) function is plotted as the dashed line in Fig. 3 by taking the same fitting parameter . It accurately captures the numerical results for the cylinder setup, including the sharp increase near the limit ϕ_f → 1. But for lower porosities, it underpredicts the data for the hole setup. Either ψ₁(ϕ_f) or ψ₂(ϕ_f) serves as an adequate representation of the interfacial permeability. They are counterparts of the two formulas for the bulk permeability k, the Kozeny–Carman formula based on flow in conduits, and the Brinkman formula on external flow around solid objects.⁹ As a further test, we have numerically computed the viscous dissipation incurred by the contraction flow, and the ψ function thus obtained agrees with the numerical data of Fig. 3 to within the difference among the pore geometries.

From ψ we can estimate the interfacial penetration length of eqn (7),

[small eta, Greek, tilde] = ψ(ϕ_f)r, (14)

and in turn the interfacial permeability η. Thus, ψ = [small eta, Greek, tilde]/r is the dimensionless interfacial penetration length scaled by the pore size. A slightly different scaling, , can be interpreted as the ratio between the interfacial permeability and the bulk permeability. We can now use eqn (14) and Fig. 3 to estimate [small eta, Greek, tilde] for real materials. Most hydrogels have small pores (nm to micron scale)^31–33 and high porosity ϕ_f > 0.9.^33–36 Thus, [small eta, Greek, tilde] can range from r up to 4r. For low-swelling-ratio gels,³⁷ biological tissues,³⁸ larger-scale porous media such as foam³⁹ and packed beds of soft particles,⁴⁰ ϕ_f can be as low as 0.5, and [small eta, Greek, tilde] will be on the order of 0.1r.

Lācis et al.⁷ computed the f factor in eqn (8) numerically in a lid-driven flow over a model porous medium. A direct comparison can thus be made with our ψ of Fig. 3. For ϕ_f = 0.75, they reported f = −5.215μ/r for a square array of circular cylinders aligned in the vortex direction of the flow outside. This translates to a penetration length [small eta, Greek, tilde]/r = 0.192, smaller than our ψ(ϕ_f = 0.75) = 0.367, but on the same order of magnitude. The quantitative discrepancy may have been due to differences in the geometric and flow setup. For example, theirs was a planar 2D flow around the cross-sections of cylinders, while ours is a 3D axial flow. Their interface was placed in the bulk of the clear fluid, 0.2r outside the plane tangential to the first row of cylinders, whereas our interface is flush with the ends of the semi-infinite cylinders and holes. Their pressure difference was volume-averaged, going 4 periodic cells into the depth of the porous medium, while our p_f is surface-averaged on the interface. Finally, their lid-driven flow is 2D on the continuum level, with variations along the interface, whereas our setup mimics the simpler 1D compression of Fig. 1. Such differences notwithstanding, both studies support the argument for an interfacial pressure drop. The numerical and semi-analytical results of Fig. 3 can be used to determine the interfacial penetration [small eta, Greek, tilde] and interfacial permeability η for a wide range of porous materials, deformable and rigid alike.

Data availability

The numerical data of this work have been generated by finite-element computations using an open-source library deal.II, which can be found at https://github.com/dealii/dealii. The version of the code employed for this study is version 9.5.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

J. J. F. acknowledges financial support by Natural Sciences and Engineering Research Council of Canada (Discovery Grants No. 2019-04162, 2024-03982, Alliance International Grant No. 586462-23). P. Y. acknowledges the financial support by National Science Foundation (Grants DMS-2012480, DMS-2309732). We thank Neil Balmforth, Uǧis Lācis, Arun Ramchandran and Boris Stoeber for helpful discussions.

References

1. R. Burridge and J. B. Keller, Poroelasticity equations derived from microstructures, J. Acoust. Soc. Am., 1981, 70, 1140–1146.
2. J. M. Carcione, Biot theory for porous media, in Wave Fields in Real Media, Elsevier, Oxford, 3rd edn, 2015, ch. 7, pp. 299–420.
3. C. W. MacMinn, E. R. Dufresne and J. S. Wettlaufer, Large deformations of a soft porous material, Phys. Rev. Appl., 2016, 5, 044020.
4. M. Minale, Momentum transfer within a porous medium. Stress boundary condition, Phys. Fluids, 2014, 26, 123102.
5. L. Li, J. Zhang, Z. Xu, Y.-N. Young, J. J. Feng and P. Yue, An arbitrary Lagrangian-Eulerian method for simulating interfacial dynamics between a hydrogel and a fluid, J. Comput. Phys., 2022, 451, 110851.
6. T. Carraro, C. Goll, A. Marciniak-Czochra and A. Mikelic, Pressure jump interface law for the Stokes-Darcy coupling: confirmation by direct numerical simulations, J. Fluid Mech., 2013, 732, 510–536.
7. U. Lācis, Y. Sudhakar, S. Pasche and S. Bagheri, Transfer of mass and momentum at rough and porous surfaces, J. Fluid Mech., 2020, 884, A21.
8. Z. Xu, J. Zhang, Y.-N. Young, P. Yue and J. J. Feng, A comparison of four boundary conditions for the fluid-hydrogel interface, Phys. Rev. Fluids, 2022, 7, 093301.
9. F. A. L. Dullien, Single-phase transport phenomena in porous media, Porous Media: Fluid Transport and Pore Structure, Academic Press, San Diego, 2nd edn, 1992, ch. 3, pp. 237–317.
10. Y.-N. Young, Y. Mori and M. J. Miksis, Slightly deformable Darcy drop in linear flows, Phys. Rev. Fluids, 2019, 4, 063601.
11. J. J. Feng and Y.-N. Young, Boundary conditions at a gel-fluid interface, Phys. Rev. Fluids, 2020, 5, 124304.
12. Z. Xu, P. Yue and J. J. Feng, Poroelastic modeling reveals the cooperation between two mechanisms for albuminuria, J. R. Soc., Interface, 2023, 20, 20220634.
13. J. T. Bartley and D. W. Ruth, Relative permeability analysis of tube bundle models, Transp. Porous Media, 1999, 36, 161–188.
14. H. K. Dahle, M. A. Celia and S. M. Hassanizadeh, Bundle-of-tubes model for calculating dynamic effects in the capillary-pressure-saturation relationship, Transp. Porous Media, 2005, 58, 5–22.
15. R. E. Larson and J. J. L. Higdon, Microscopic flow near the surface of two-dimensional porous media. Part 1. Axial flow, J. Fluid Mech., 1986, 166, 449–472.
16. R. E. Larson and J. J. L. Higdon, Microscopic flow near the surface of two-dimensional porous media. Part 2. Transverse flow, J. Fluid Mech., 1987, 178, 119–136.
17. A. S. Sangani and S. Behl, The planar singular solutions of Stokes and Laplace equations and their application to transport processes near porous surfaces, Phys. Fluids A, 1989, 1, 21–37.
18. N. Martys, D. P. Bentz and E. J. Garboczi, Computer simulation study of the effective viscosity in Brinkman's equation, Phys. Fluids, 1994, 6, 1434–1439.
19. D. F. James and A. M. Davis, Flow at the interface of a model fibrous porous medium, J. Fluid Mech., 2001, 426, 47–72.
20. T. Carraro, C. Goll, A. Marciniak-Czochra and A. Mikeli, Effective interface conditions for the forced infiltration of a viscous fluid into a porous medium using homogenization, Comput. Methods Appl. Mech. Eng, 2015, 292, 195–220.
21. U. Lācis and S. Bagheri, A framework for computing effective boundary conditions at the interface between fluid and a porous medium, J. Fluid Mech., 2017, 812, 866–889.
22. S. Haber and R. Mauri, Boundary conditions for Darcy's flow through porous media, Int. J. Multiphase Flow, 1983, 9, 561–574.
23. O. Coussy, Poromechanics, Wiley, 2004.
24. J. Koplik, H. Levine and A. Zee, Viscosity renormalization in the brinkman equation, Phys. Fluids, 1983, 26, 2864–2870.
25. J. H. Masliyah, G. Neale, K. Malysa and T. G. M. Van De Ven, Creeping flow over a composite sphere: Solid core with porous shell, Chem. Eng. Sci., 1987, 42, 245–253.
26. U. Lācis, G. A. Zampogna and S. Bagheri, A computational continuum model of poroelastic beds, Proc. R. Soc. A, 2017, 20160932.
27. M. F. Tachie, D. F. James and I. G. Currie, Velocity measurements of a shear flow penetrating a porous medium, J. Fluid Mech., 2003, 493, 319–343.
28. M. Agelinchaab, M. F. Tachie and D. W. Ruth, Velocity measurement of flow through a model three-dimensional porous medium, Phys. Fluids, 2006, 18, 017105.
29. J. K. Arthur, D. W. Ruth and M. F. Tachie, PIV measurements of flow through a model porous medium with varying boundary conditions, J. Fluid Mech., 2009, 629, 343–374.
30. M. A. Rogerson and Y. L. Yeow, Representation of Stokes flow through a planar contraction by Papkovich-Fadle eigenfunctions, Trans. ASME, 1999, 66, 940–944.
31. M. M. Chui, R. J. Phillips and M. J. McCarthy, Measurement of the porous microstructure of hydrogels by nuclear magnetic resonance, J. Colloid Interface Sci., 1995, 174, 336–344.
32. Y.-C. Chiu, M.-H. Cheng, H. Engel, S.-W. Kao, J. C. Larson, S. Gupta and E. M. Brey, The role of pore size on vascularization and tissue remodeling in PEG hydrogels, Biomaterials, 2011, 32, 6045–6051.
33. A. Salerno, R. Borzacchiello and P. A. Netti, Pore structure and swelling behavior of porous hydrogels prepared via a thermal reverse-casting technique, J. Appl. Polym. Sci., 2011, 122, 3651–3660.
34. N. Annabi, J. W. Nichol, X. Zhong, C. Ji, S. Koshy, A. Khademhosseini and F. Dehghani, Controlling the porosity and microarchitecture of hydrogels for tissue engineering, Tissue Eng., Part B, 2010, 16, 371–383.
35. Y. Li, O. S. Sariyer, A. Ramachandran, S. Panyukov, M. Rubinstein and E. Kumacheva, Universal behavior of hydrogels confined to narrow capillaries, Sci. Rep., 2015, 5, 17017.
36. S. Kovačič and M. S. Silverstein, Superabsorbent, high porosity, PAMPS-based hydrogels through emulsion templating, Macromol. Rapid Commun., 2016, 37, 1814–1819.
37. J. Yoon, S. Cai, Z. Suo and R. C. Hayward, Poroelastic swelling kinetics of thin hydrogel layers: comparison of theory and experiment, Soft Matter, 2010, 6(23), 6004.
38. R. Penta, L. Miller, A. Grillo, A. Ramírez-Torres, P. Mascheroni and R. Rodríguez-Ramos, Porosity and diffusion in biological tissues. Recent advances and further perspectives, in Constitutive Modelling of Solid Continua, ed. J. Merodio and R. Ogden, Springer International Publishing, Cham, 2020, pp. 311–356.
39. T. Lutz, L. Wilen and J. Wettlaufer, A method for measuring fluid pressure and solid deformation profiles in uniaxial porous media flows, Rev. Sci. Instrum., 2021, 92, 025101.
40. D. R. Hewitt, J. S. Nijjer, M. G. Worster and J. A. Neufeld, Flow-induced compaction of a deformable porous medium, Phys. Rev. E, 2016, 93, 023116.

---