## Abstract

Surfaces covered by bristles, hairs, polymers and other filamentous structures arise in a variety of natural settings in science such as the active lining of many biological organs, e.g. lungs, reproductive tracts, etc., and have increasingly begun to be used in technological applications. We derive an effective field theory for the elastohydrodynamics of ordered brushes and disordered carpets that are made of a large number of elastic filaments grafted on to a substrate and interspersed in a fluid. Our formulation for the elastohydrodynamic response of these materials leads naturally to a set of constitutive equations coupling bed deformation to fluid flow, accounts for the anisotropic properties of the medium, and generalizes the theory of poroelasticity to these systems. We use the effective medium equations to study three canonical problems—the normal settling of a rigid sphere onto a carpet, the squeeze flow in a carpet and the tangential shearing motion of a rigid sphere over the carpet, all problems of relevance in mechanosensation in biology with implications for biomimetic devices.

## 1. Introduction

Soft, porous beds made of arrayed, slender, deformable, filamentous structures, anchored to a substrate and permeated by a viscous liquid are ubiquitous in nature, and increasingly seen in technology. These composites span a spectrum of scales in terms of structure and are versatile in terms of function. Examples in nature range from ciliary beds enabling respiratory, circulatory, urogenital and ambulatory functions (Rensch *et al.* 1983; Schurch *et al.* 1990; Damiano & Stace 2005; Smith *et al.* 2007), motion sensing otoliths and stereocilia in ears (Howard & Ashmore 1986), lubricating filamentous aggrecan brushes (Han *et al.* 2007) and lubricin macromolecules (Zappone *et al.* 2007) in human joints. Figure 1 depicts some examples of poroelastic filamentous beds in nature. In the technological world, analogous fibrous interfaces have been employed both as sensors as well as actuators in machines (Tung & Kim 2006), fluidic pumps driven by beating elastic filaments (Kim & Netz 2006) and in nanorod arrays used in DNA analysis and separation (Evans *et al.* 2007). These biphasic material interfaces are constituted primarily of an elastic skeleton permeated by an incompressible fluid. Imposition of a pressure gradient drives fluid flow through such a composite and causes it to deform elastically; alternately deforming the composite actively results in fluid flow.

We illustrate the coupling between deformation and fluid flow using an example motivated by biological mechanosensing structures which resemble the schematic in figure 1*a*. A sphere of radius *R* moves through a liquid in the vicinity of a filamentous bed comprising thin elastic bristles. The bed is completely infiltrated by the same fluid. Each filament is of length ℓ, radius *a* and spaced apart such that the inter-filament distance is Δ so that the pore size ℓ_{p}∼*a*(Δ/*a*−2) and the solid area fraction *ϕ*_{s}∼*O*(*a*/Δ)^{2}. The macroscopic or system scale is ℓ_{s}. The bed height ℓ_{b}∼ℓ in the undeformed state.^{1} The filaments are anchored at the base to an impermeable surface. The permeability of the bed depends on the pore size and the details of the arrangement of fibres and can be estimated as *K*∼ℓ^{2}_{p}∼*a*^{2}*f*_{1}(*ϕ*_{s}), the function *f*_{1} depending on the specifics of the array geometry. The bending modulus of each filament *B* depends on ℓ, *a* and the Young’s modulus, *E*. The effective elastic modulus of the bed *C*_{eff} however, depends on *E*, ℓ, *a*, as well as on *ϕ*_{s} and on the details of the anchoring at the base shoring up the fibres.^{2} For simplicity, we neglect the intrinsically anisotropic nature of the bed.

Consider first the sphere translating parallel to the bed with speed *V* _{0} at a height ℓ_{b}+*h*_{0}∼ℓ+*h*_{0}. As the sphere moves, it sets into motion the fluid in the interstitial region. This results in an externally applied stress acting on the bed. The time for which the bed is forcibly deformed is roughly the convective time scale *τ*_{c}=ℓ_{s}/*V* _{0}. The deformation of the bed initiates fluid flow through the bed away from the deformed region. If the sphere moves very fast, the fluid trapped in the deforming region is unable to escape and thus the composite bed resists the deformation as if it were incompressible with a large elastic modulus. If however, the sphere moves very slowly, the fluid has time to escape through the porous structure and equilibrate with fluid far away from the interaction region. The elastic solid skeleton then supports the bulk of the imposed stress and the apparent composite modulus is lower. This behaviour may be quantified using a *poroelastic time scale* *τ*_{p}=(*μ*/*C*_{eff})(ℓ^{2}/*K*)∼(*μ*/*C*_{eff})(ℓ/ℓ_{p})^{2}*f*_{2}(*ϕ*_{s},*a*/ℓ). The ratio *τ*_{c}/*τ*_{p} controls the apparent stiffness of the bed and the magnitude of the resulting deformation and consequently the lift force acting on the translating sphere.

For a sphere translating towards the surface with a normal speed *V* _{0}, fluid is driven into the poroelastic bed and then eventually streams out radially due to the no-flux condition at the impermeable base. When the sphere moves very slowly, fluid has time to move through the pores and the bed deformation, fluid flux and the position of the sphere are all quasi-statically coupled.

For both of the above-mentioned cases, one may ask the following questions. How does the bed deform in response to the forcing? What are the resultant forces on the sphere resisting its motion? To answer these questions, it is essential to accurately take into account the deformation–flow coupling as well the effects of anisotropy. A natural foundation underlying models of such soft, porous, fluid infiltrated composites is the theory of poroelasticity first proposed by Biot (1941, 1955) more than half a century ago and put on a firm footing using both homogenization techniques and the intuitive mixture-theory formulation for bulk and low-dimensional objects (Bowen 1973; Rice & Cleary 1976; Burridge & Keller 1981; Mei & Auriault 1989; Wang 2000; Skotheim & Mahadevan 2004). Building on these studies, we describe and characterize carpet-like interfaces on length (spatial) scales that are large compared with the microscopic length scales characterizing the pore structure using an effective theory that takes into account the two-way coupling between fluid flow and elasticity and intrinsically anisotropic characteristics of the material. Experimental support for this approach comes from studies of naturally occurring macromolecular brushes (Han *et al.* 2007; Zappone *et al.* 2007) demonstrating that these brushes exhibit both viscoelastic and poroelastic responses under loading.

We depart from earlier attempts to study these materials in three important ways. (i) We account for the elastic nature of the bed and its ability to store and release energy in response to fluid flow, improving earlier efforts (Keller *et al.* 1975; Feng & Weinbaum 2000; Du *et al.* 2004; Damiano & Stace 2005; Han *et al.* 2005; Smith *et al.* 2007) that neglect the effects of normal deformations, their concomitant elasticity and the complete tensorial nature of the coupling between local flow and local deformation. (ii) We illustrate how the anisotropic nature of both bed permeability and elasticity due to the alignment of the fibres may be properly accounted in the theory. Our work builds on previous studies of anisotropic poroelasticity (Biot 1941, 1955; Barry & Holmes 2001), but includes and clarifies the role of boundary conditions in controlling the elastic properties of the bed. (iii) Finally, we show that matching of shear and normal stress at the composite-free liquid interface and complete flow–deformation coupling are both accomplished naturally using a Darcy-type model (Darcy 1856) describing relative flow between the solid and fluid phases in the soft medium. This naturally resolves the problem of boundary conditions at soft, porous interfaces and alleviates the need for ad hoc models based on the Brinkman equation (Fredrickson & Pincus 1991; Feng & Weinbaum 2000; Du *et al.* 2004) that introduce higher order corrections that are asymptotically inconsistent.

It is useful to classify the wet poroelastic beds into two general categories illustrated in figure 2*a*—*carpets* and *brushes*. Carpets have stiff ordered fibres with a reference state that consists of straight fibres, while brushes are made of soft, dense, disordered beds typically made up of permanently anchored polymeric chains with a reference state given by the initial bed height and spatial polymer concentration field. We will assume that individual elastic filaments in the brush and carpet typically do not touch each other even when the bed deforms since the lubrication forces associated with the permeating fluid prevents such contact for sufficiently smooth surfaces. Then fibre–fibre contact forces are negligible.

In §2, we formulate the equations of motion and describe the physical characteristics embodied in the formulation, while in §§3 and 4 we evaluate the elastic and permeability tensors that characterize the bed in term of microscopic parameters for disordered brushes and ordered bristles and carpets. These enable us to close the set of equations formulated in §2. Section 5 is dedicated to a description of the boundary conditions at the impermeable and permeable interfaces, correcting earlier ad hoc approximations. In §6, we analyse model problems in the context of these poroelastic carpets: the settling of a sphere, compression of a brush and deformation due to a sphere translating past a bed. Finally, in §7, we conclude with a discussion of our results in the context of a broader class of phenomena where our framework might be relevant.

## 2. An effective poroelastic field theory

Formalizing these ideas in the context of an asymptotic theory valid for both anisotropic carpets and brushes, we note that soft porous beds are typically characterized at the microscale level by a pore size and at the macroscale level by a externally determined system size. If **∇** and **∇**_{p} are the gradients on the system scale and the pore scale and *P* and *P*_{p} be characteristic pressure scales, fluid stress balance at the pore scale implies that the sum of the macroscopic pressure gradient driving the flow, **∇***P*, and the microscopic pressure gradient, **∇**_{p}*P*_{p}, balances the viscous resistance of the pore fluid . Momentum balance in the fluid at the pore scale then yields
2.1
so that when ℓ_{p}/ℓ_{s}≪1, and the dominant contribution to the fluid stress tensor *σ*^{F} comes from the pressure, so that *σ*^{F}=−*P*. Consequently on scales large compared with the distance between the fibres, but small compared with the system size, we may model the composite bed as an effective continuum, although it might be anisotropic and inhomogeneous.

In order to make progress and formulate the governing equations for the elastohydrodynamic response, we limit ourselves to a consideration of linearly elastic behaviour corresponding to small strains, wherein the linear strain ** ϵ** is related to the displacement field

**u**, measured relative to the reference state via 2.2 We now define two elastic limits of the test volume element—the jacketed limit where the composite is enclosed by an impermeable surface so that fluid cannot flow out of the test element and the freely drained (or draining) limit which corresponds to the fluid being able to freely move in and out of the test volume. The simplest form for the averaged stress tensor characterizing the composite then arises by considering the linear superposition of the fluid stress which, at leading order, is due to the pore pressure and the elastic stress due to the deformation of the solid skeleton. Then we may write 2.3 Here

**C**is the fourth-order tensor characterizing the effective elastic moduli of the poroelastic skeleton and is in fact the drained modulus corresponding to the deformation response—when fluid can move freely or equivalently when there is no interstitial fluid. Thus,

**C**is determined by solving canonical microscale problems that involve analysis of independent deformation modes of a unit cell, and must clearly vanish as the solid volume fraction,

*ϕ*

_{s}vanishes. The tensorial quantity multiplying the pore pressure in equation (2.3)

*viz*.,

**, is a measure of fluid fraction changes due to the compressibility of the solid making up the skeleton and is thus dependent on both material properties as well as microstructural properties. When inertial effects are negligible, the composite effective stress**

*α***and external volumetric body forces**

*σ***F**

_{B}are in balance. If the constituent solid (subscript s) comprising the skeleton satisfies the constitutive law

*σ*_{s}=

**C**

_{s}:

*ϵ*_{s}, then following steps in Rice & Cleary (1976), we may show that 2.4 with

**being the identity tensor.**

*δ*^{3}Note that for an incompressible skeleton and interstitial fluid, we get

**=**

*α*

*α*_{IN}=

**. For an isotropic, compressible porous medium, equation (2.4) can be simplified to yield,**

*δ***=**

*α*

*α*_{C}=

*α*

_{C}

**=(1−**

*δ**K*

_{d}/

*K*

_{s})

**with**

*δ**K*

_{d}being the bulk modulus of the drained composite.

The local fluid content within material elements of the composite changes due to stress that results in deformation of the skeleton. For small inertia-less deformations, Darcy’s law describes the relative motion of the fluid and the solid phases at the pore level. In the reference rest state, the porosity and solid volume fraction are homogeneous on pore-averaged system length scales. When the porous solid is elastic, changes in liquid content per unit volume result due to change in the volume of the solid content as the solid skeleton moves. To leading order, the pore-scale-averaged discharge velocity of the fluid *relative* to the deforming solid skeleton, **Q** is then related to the system scale pressure gradient by
2.5
where **K** is cell-averaged anisotropic tensorial version of the classical Darcy permeability. To understand variations in flow, we first define the porosity, *ϕ*_{s}=*V* _{F}/*V* and fluid mass content, *m*_{F}=*M*_{F}/*V* =*ρ*_{F}*ϕ*_{F}, where *V* _{F} is the volume of fluid, and *M*_{F} is the mass of fluid contained in unit cell of porous material which would occupy volume *V* in an unstressed reference state. With *K*_{F} being the bulk modulus of the pure pore liquid, let variations in density and pressure be small relative to the initial state so that the fluid exhibits adiabatic compressibility, conservation of fluid mass in an elemental volume of the composite then yields
2.6
where
2.7
Fluid–solid coupling due to compressibility of either constituent enters via
2.8
To uncover the meaning of the different terms in equations (2.7) and (2.8) we first consider the limit of an isotropic, compressible composite. In this limit and for a jacketed composite, *D*(*m*_{F})=0. In the general unjacketed case, equation (2.7) reduces to *D*(*m*_{F})=*ρ*_{F}(*α*_{C}**∇**⋅**u**+*β*_{C}*P*). The parameter *β*_{C}≡(*ϕ*_{F}*K*_{s}+(*α*_{C}−*ϕ*_{F})*K*_{F})/(*K*_{F}*K*_{s}) can be rewritten as and is then seen to be a measure of the jacketed modulus, *K*_{u}. In other words, for incompressible constituents (, ) we find *β*=*β*_{IN}=0 with ∂(**∇**⋅**u**)/∂*t*=(*K*/*μ*)∇^{2}*P* matching the analysis of Barry & Holmes (2001). We slao note that in the limit , equation (2.7) matches expressions for isotropic media derived by Rice & Cleary (1976), while the expression for *β*, equation (2.8), can be rewritten in the form derived earlier by Loret *et al.* (2001). We emphasize that equation (2.7) embodies the role played by skeleton compressibility (due to porosity), fluid compressibility and the compressibility of the bulk solid comprising the skeleton.

Finally, combining equations (2.5)–(2.8), we obtain the generalization of Darcy’s law to anisotropic, compressible media
2.9
In deriving this equation we have ignored nonlinear terms arising from the product *ρ*_{F}**Q** in equation (2.6). Equations (2.2)–(2.4), (2.8) and (2.9) involve the deformation field, **u** and the pore pressure *P* as unknowns.^{4} To complete the formulation of the problem, we need to specify the effective elastic modulus tensor **C** and the effective fluid permeability tensor **K** as well as the boundary conditions that accompany the field equations.

## 3. Disordered brushes and carpets

We start with a specialization of our theory to a polymer brush obtained by grafting polymers with a characteristic monomer size *a*_{m} onto surfaces with density *Γ*—(the number of monomers per unit area) in the presence of a good solvent of viscosity *μ*. We assume that the equilibrium brush thickness is much larger than the mesh size, the monomer size and the grafting length scale (distance between chains at the base), so that the polymers form a spongy, compressible fluid-infiltrated bed. Furthermore, in good solvent conditions the polymer coil swells and the various filaments tend to intersperse so that the brush is transversely isotropic.

For a Newtonian, incompressible solvent and incompressible polymer units, the total composite stress in the soft medium reduces to a simple form. Thus, while, ** α**=

**, the elasticity tensor**

*δ***C**satisfies 3.1 Substituting equation (3.1) in equations (2.3) and (2.4) and setting external body forces to zero yields

^{5}, 3.2 Here,

*P*

_{T}is the net fluid pore pressure that includes concentration-dependent osmotic effects. The pore pressure at steady state corresponds to a local equilibrium between the osmotic pressure, ambient equilibrium fluid pressure and the elastic stress. In the absence of any external force or flow, the elastic stress due to network deformation balances the osmotic pressure. The fluid flux through the soft brush couples to deformation via a permeability that now depends on polymer concentration, 3.3 where

*ζ*

_{H}(

*c*) is the pore size, or equivalently a hydrodynamic screening length that depends on

*c*(Fredrickson & Pincus 1991), the local polymer concentration. The shear and bulk moduli for the brush follow the simple scaling laws 3.4 where

*χ*

_{1}and

*χ*

_{2}are constant with

*χ*

_{1}≪

*χ*

_{2}.

The hydrodynamic screening length *ζ*_{H}(*c*) is a function of local mean-field polymer concentration, *c*(*z*) and the monomer size *a*_{m}. We now restrict our descriptions to macroscopic length scales *L*_{c} much larger than the mesh size of the brush and to small and moderate volume fractions. For grafted polymers in a good solvent with volume fractions in the semi-dilute regime a reasonable estimate for *ζ*_{H}(*c*) is the equilibrium concentration correlation length (deGennes 1979)—thus . For grafted polymer brushes, as a first approximation (Alexander 1977) we can choose *c* to be an appropriate spatially averaged value. For instance, when considering the compression of a polymer brush of thickness *h*, we could replace *c*(*z*) with the spatial average *c*∼*Γ*/*h*. Adsorbed polymer brushes are complicated to treat (deGennes 1979) as *ζ*_{e}∼*z*, where *z* is the axial position perpendicular to the grafting surface. In general, the hydrodynamic screening length as well as elastic constants in equations (3.2) and (3.3) are functions of the coordinate normal to the plane of transverse isotropy.

## 4. Ordered bristles and carpets

We now turn to a bed of bristles or a carpet of cylindrical fibres made of an isotropic elastic material, each of radius *a* and length ℓ≫*a* arranged in a regular hexagonal lattice as shown in figure 2*a*. The fibres are anchored to a rigid, impermeable fixed substrate (corresponding to the *x*–*y* plane). We assume that the fibres are not perpendicular to the bed since that is a very non-generic situation, and will instead assume that they are tilted on average at a small angle, *θ*_{0} with the *z*-axis (). The inter-fibre distance Δ≡2*a*(1+*ϵ*_{c}), and the system length in the lateral direction, . For *z*≫ℓ, provided ℓ≫Δ, the fibre–fluid composite is effectively a permeable, poroelastic bed with 2*ϵ*_{c}=ℓ_{p} being the characteristic pore size. The coordination number for this geometry *Z*_{c}=6 and thus the maximum packing fraction is given by . The solid area fraction, or equivalently, the volume fraction, is given by in terms of which both the elastic modulus tensor **C** and the permeability tensor **K** may be expressed.

Symmetry dictates that *C*_{ijkl}=*C*_{jikl}, *C*_{ijkl}=*C*_{jilk} and *C*_{ijkl}=*C*_{klij}. Defining the column vector (*ϵ*)^{T}=(*ϵ*_{1},*ϵ*_{2},*ϵ*_{3},*ϵ*_{4},*ϵ*_{5},*ϵ*_{6}) = (*ϵ*_{11},*ϵ*_{22},*ϵ*_{33},2*ϵ*_{23},2*ϵ*_{13},2*ϵ*_{12}) and analogously, we write the stress–strain relationship characterizing the drained response of the skeleton concisely as and where we have specifically used *τ* to denote the skeleton stress, as distinguished from the composite stress ** σ**. Note that and are

*not*the fourth-order elasticity and compliance tensors but matrices that

*do not*transform like tensors. The stiffness and compliance matrix can be written in terms of just five constants—the transverse (in-plane—with 1 and 2 denoting the

*x*and

*y*axes) Young’s moduli,

*E*

_{p}=

*E*

_{11}=

*E*

_{22}, the longitudinal Young’s modulus,

*E*

_{⊥}=

*E*

_{3}, the two transverse in-plane Poisson ratios,

*ν*

_{p}=

*ν*

_{12}=

*ν*

_{21}, the Poisson ratio,

*ν*

_{⊥}=

*ν*

_{zp}that denotes the effect of in-plane deformations due to axial loading and the in-plane shear modulus,

*G*

_{⊥}. The compliance matrix takes the form 4.1

To determine the non-zero elastic constants in terms of the solid volume fraction *ϕ*_{s} and evaluate the dependence on boundary conditions of the composite material,^{6} it suffices to consider the mechanics of a single fibre. The unit vector perpendicular to the substrate plane, **e**_{z} is also approximately the axis about which the fibres are arrayed symmetrically as shown in figure 3. A typical situation involves fibres with bending stiffness *B* clamped at the base to the underlying substrate by torsional springs of stiffness *K*_{s}, as seen often in natural and artificial systems. For example, hairs and cilia are strongly anchored at the basal body but can bend, while actin bundles and bristles often are fairly stiff but can easily pivot about their base. Since the fibres are unlikely to be precisely normal to the base, we assume that the fibres are inclined at a small angle *θ*_{0} with the *z*-axis (*θ*_{0}≪1). The small offset angle also introduces out of plane bending forces, but these are an order of magnitude smaller than the dominant bending moment.

Consider first a test fibre subject to shearing forces as in figure 3 (cases I and II). Denoting derivatives with respect to *z* by primes, the equation for the shape of a fibre *h*(*z*) relative to the vertical can be written *Bh*′′′′=*q*_{ℓ}, *h*(0)=0, *K*_{s}(*h*′(0)−*θ*_{0})+*Bh*′′(0)=0,*h*′′(ℓ)=0 and *h*′′′(ℓ)=*P*_{ℓ}, where *P*_{ℓ} is the load at the tip, and *q*_{ℓ} is any distributed load along its entire length, ℓ. At a scaling level, the tip deflects by an amount that scales as *P*_{ℓ}ℓ^{3}/*B*, when *q*=0 a result that follows from a simple energy balance. In contrast, a rigid fibre that is supported by a torsional spring of stiffness *K*_{s} at the base when subject to a point load will deflect by amount that scales as *P*_{ℓ}ℓ^{2}/*K*_{s}. Similar scaling relations may be derived for the case when *q*_{ℓ}≠0, and also when the test fibre is compressed by a perpendicular force *P*^{⊥} (figure 3, case III). Expressions for the transverse displacement *h*^{∥}=*x*_{Tip} and normal displacement are listed in Table 1.

To convert the response of a single fibre to that for the entire bed, we consider the elastic response of the composite in the fully drained limit—a useful way to visualize this is to consider a slab of the carpet encapsulated by a membrane that is completely permeable to the fluid but retains the fibres inside. We first consider the shear response of the bed associated with a stress per unit area——acting parallel to the (*x*−*y*) plane at the upper edge of the brush. The typical shear strain associated with this load . The number of fibres per unit area that support this load is *ϕ*_{s}(*πa*^{2})^{−1} (here we consider the fully drained limit to calculate the response of the skeleton). The effective bed shear modulus is then given by . For instance, when , one obtains, *G*_{⊥}=*ϕ*_{s}*K*_{s}(2*π*ℓ*a*^{2})^{−1} while when , we obtain *G*_{⊥}=3*Bϕ*_{s}(2*π*ℓ^{2}*a*^{2})^{−1}. Next, to calculate *E*_{3}=*E*_{⊥}, we consider a compressive stress, acting on the same elemental test area that leads to a compressive normal strain so that the elastic modulus to leading order is given by . When , , while when the torsional spring is very stiff and bending is the dominant mode of deflection, . The composite in-plane Poisson ratio, *ν*_{p}=0 is zero. When the bed of fibres is subjected to a compressive stress at the top boundary, each fibre individually buckles in response. However, since there is no fibre–fibre contact and hence no correlation between the fibres, so that *ν*_{⊥}=*ν*_{zp}=*ν*_{pz}=0. Finally, we consider the in-plane modulus *E*_{p}, which characterizes the response of the carpet as it is squeezed by normal forces laterally.^{7} Consider a normal stress acting on an area of unit width ℓ_{y} and height equal to the bed height ℓ_{z}=ℓ. A calculation similar to that carried out earlier then yields the estimates for and for .

We now turn to estimates for the pore-averaged macroscale permeability tensor **K**. As the bed is *not* unbounded, one can speak of only apparent permeabilities defined to establish a relationship between the mean flow and the macroscopic pressure gradient. For small fibre inclinations 0<*θ*_{0}≪1 and small deformations, the permeability is close to that for perfectly aligned straight fibres. In the interior of the poroelastic bed outside of the boundary layers (see figure 2), transverse isotropy at scales larger than the pore size implies that
4.2
where is the in-plane permeability and is the vertical permeability. In the limit of very small Reynolds number, the Stokes equations characterize the flow both at pore scales as well as at scales of the order of the fibre radius. The linearity of the Stokes equations and time reversibility implies that the mean flow is linear in the imposed pressure gradient and thus the permeability tensor can be obtained from a single canonical calculation. Since the pressure drop across a unit cell is in effect balanced by the net force per unit length per cylinder in the cell multiplied by the number of cylinders per unit area, one can write, with *F* being the axial (∥) or transverse (⊥) mean drag force per unit length per cylinder, *U*, the appropriate superficial (mean) velocity in the bed and *ϕ*_{s} the areal/volume solid fraction. The exact forms of the reduced force, *f* and the permeabilities, *K*_{∥} and *K*_{⊥} depend on the details of the packing and on the solid volume fraction, *ϕ*_{s}. In the following, we focus on the case of hexagonally packed fibres and summarize the single fibre mean-field problems that have to be solved to evaluate the permeabilities.

In the dilute (D) limit as 0<*ϕ*_{s}≪1, the fluid velocity field governed by the Stokes equations varies logarithmically with distance from the cylindrical fibre with the divergence cut-off by finite pore size effects via the dependence of velocity on the parameter (Δ/*a*). Matching the Stokes solution valid for distances that are smaller than *O*(*a*) with a (renormalized) mean field far field velocity field that takes into account the presence of other cylinders yields (Sangani & Acrivos 1982)
4.3
A similar analysis for the axial permeability yields (Sparrow & Loeffler 1959)
4.4

In the dense limit (C) as , with *ϕ*^{m}_{s} being the critical (maximum) packing at the solid fraction at which the fluid cannot flow through the bed, most of the pressure drop required to maintain the flow is associated with the resistance due to the small gap between adjacent cylinders. Lubrication theory implies that the fluid velocity through this gap is and the pressure drop across the gap, obtained by matching with contributions from viscous shear, is . Using this result and accounting for the hexagonal packing as in Sangani & Acrivos (1982), we obtain
4.5
The axial permeabilities for a hexagonal array are obtained from an adaptation of the Hagen–Poiseuille relationship and yeild
4.6
While one may use an interpolation for moderate volume fractions , equations (4.3) and (4.4) are observed to provide good estimates of the permeability for volume fractions as large as *ϕ*_{s}∼0.2. For comparison, we note that the Carman–Kozeny relationship (Happel & Brenner 1973) yields for an isotropic porous medium.

## 5. Boundary conditions

The correct boundary conditions at a porous interface have been a source of long-standing debate owing to the problem of matching of both normal and tangential stresses at the interface (see Beavers & Joseph 1967; Saffman 1971; Taylor 1971; Hou *et al.* 1989; Barry & Holmes 2001; Jager & Mikeli 2009; Nield 2009). Previous models for the flow through porous media (Fredrickson & Pincus 1991; Feng & Weinbaum 2000; Du *et al.* 2004) use the Brinkman equation and its variants that take the form
where the exact value of the micro-structure dependent effective viscosity *μ*_{eff} is the subject of significant disagreement. While this equation is applicable in many scenarios, using the Brinkman equation to model macroscale flow through a poroelastic bed is not mathematically correct, since it violates the asymptotic separation between the viscous contribution and the pressure in the fluid stress tensor, as indicated in §2. In fact, it is also not required since shear stresses can be borne by the deformation of the solid skeleton.

Since the filamentous bed—be it a carpet or a brush—is bounded from below by a rigid, impenetrable wall, the mean-field deformation and the normal component of the macroscopic pore-averaged flux vanish at the wall, so that 5.1 Note that this implies that the pore-averaged normal fluid velocity component is zero as well and is thus consistent with the point-wise fluid velocity field vanishing at the impermeable rigid wall.

At the upper boundary that separates the brush from the pure fluid, we must impose continuity of the averaged traction as well as a condition on the fluid flux. For a fluid flow above the carpet with velocity field with the concomitant fluid stress field , continuity of traction at the interface I implies that
5.2
This stress balance yields a vector equation with the two independent components for normal stress balance *and* the tangential stress balance. Tangential stresses exerted by interfacial fluid flowing past the bed are borne by the deformation of the solid skeleton—the pore fluid itself contributes only to the isotropic stress via the pore pressure. Details of the velocity field in the porous medium are averaged into the effective fluid flux flowing in accordance with the system scale pressure gradient and viscous effects at the pore scale *are already* taken into account via the effective pore-averaged permeability in the Darcy term. Thus the matching condition (5.2) involves tangential and normal stresses of the composite and the pore pressure and *not individual velocity components*.

The general boundary condition on the *fluid flux* at the fluid–bed interface may be formulated by assuming the upper fluid–bed interface to be bounded by a semipermeable thin skin of thickness *d*_{f} and isotropic permeability *K*_{f}. Imposing continuity of the fluid flux at the interface, we obtain
5.3
where *P**|_{I} is the pressure difference across the skin. This expression thus includes both the case of an impermeable skin for which *K*_{f}=0, but is more general, being a consequence of flux continuity across the skin.

Equations (5.1)–(5.3) differ significantly from the boundary conditions based on the Brinkman model for the porous bed that involve point-wise velocity fields (described in appendix B). Equations (2.2)–(2.4), (2.8) and (2.9) along with the boundary conditions (5.1)–(5.3) complete the formulation of the problem for the effective field-theoretic description of the poroelastic bed of ordered or disordered bristles or brushes, with the relevant elasticity and permeability tensors given by equations (3.1) and (3.4) for disordered brushes and by equations (4.1)–(4.5) for ordered carpets.

## 6. Model problems: normal and tangential deformations of a bed

We now discuss canonical problems that characterize the response of a carpet or a brush to normal and tangential forcing. Schematic sketches of these scenarios are depicted in figure 4. In all cases, bed deformation drives an anisotropic fluid flow through the bed.

### (a) Compression of a poroelastic bed by a sedimenting sphere

We consider first the gravity-driven sedimentation of a small rigid sphere of density *ρ*_{s} and radius *R* onto a poroelastic fibrous bed of thickness *H*_{ℓ}, infiltrated by a Newtonian fluid with viscosity *μ* and density *ρ*_{f}, as shown in figure 4*a*. The density difference is Δ*ρ*=*ρ*_{s}−*ρ*_{f} with Δ*ρ*/*ρ*_{f}≪1. We choose our coordinate system so that the centre of mass of the sphere lies a height *R*+*h*(*t*)+*H*_{ℓ} away from the rigid base of the carpet. Measuring radial position from the axis of symmetry, let *H*_{p}(*r*,*t*) be the position of the carpet interface as measured from the rigid base, the deformation of the interface relative to the initial thickness be *δ*_{b}(*r*,*t*)=*H*_{ℓ}−*H*_{p}(*r*,*t*), and the velocity of the sphere be d*h*/d*t*. The width of the gap between the bed and the sphere is thus *h*+*δ*_{b}(*r*,*t*) with *h*>0 initially, but eventually we have *h*<0.

The settling of the sphere and resulting bed deformation are coupled quasi-statically, a situation of some relevance in mechanosensation, as statoliths, otoliths and amyloplasts settle on sensory cilia in a variety of animal organs and plants. For negligible sphere and fluid inertia, gravitational sedimentation is resisted by viscous and elastic forces. As the particle moves close to the wall with *H*_{ℓ}≪*h*≪*R*, the fluid pressure in the gap and flow through the bed resists motion. Power balance then reads:
6.1
a form applicable to both thin (*R*≫*H*_{ℓ}) and thick (*R*≪*H*_{ℓ}) beds.

The first term on the right-hand side of equation (6.1) corresponds to viscous dissipation due to the radial flow of liquid at a speed *O*((d*h*/d*t*)(*R*/(*h*+*δ*_{b}))^{1/2}) through a gap of size *O*(*h*+*δ*_{b}). The second term corresponds to the viscous dissipation due to radial Darcy flow as the fluid is squeezed between the fibres whose number scales as *O*(*ϕ*_{s}(*Rδ*_{b})(*πa*^{2})^{−1}) over a volume . The third term corresponds to the power associated with the deformation of the elastic bed with a compressive modulus *O*(*E*_{⊥}) which is strained by an amount , over a volume . We assume here that the deformation itself is small compared with the thickness of the bed *δ*_{b}≪*H*_{ℓ}.

#### (i) Viscous regime

During the initial stages of approach, *δ*_{b}∼0 so that the second and third terms on the right-hand side of equation (6.1) are negligible compared to the first. Then the gap *h*≥*R*, the first term reduces to *μ*(d*h*/d*t*)^{2}*R*, and we find that the particle falls at the Stokes settling speed |d*h*/d*t*|=*u*_{s}∼Δ*ρR*^{2}*g*/*μ*. When *H*_{ℓ}≤*h*≪(*Rh*)^{1/2}, the bed is deformed weakly, and the gap is *O*(*h*). Provided the flux of fluid through the interface into the bed is small compared with flow through the gap i.e. , the first term on the right-hand side of equation (6.1) reads *μ*(d*h*/d*t*)^{2}*R*^{2}/*h* so that the settling velocity now becomes^{8} |d*h*/d*t*|∼Δ*ρRgh*/*μ*.

#### (ii) Poroelastic regime

As the sphere sinks into the bed *h*=−*δ*_{b}<0, and fluid flow occurs mainly through the bed. Now the settling velocity is determined by balancing the gravitational power with dissipation rate due to the Darcy flow and the energy stored in the bed. Initially, for , i.e. , viscous dissipation dominates over the elastic storage term. Then (for ) and (for ). Since the rate at which the bed is being compressed decreases with increasing *δ*_{b}, the sphere slows down so that——then a larger fraction of the gravitational work is stored as elastic energy, and a progressively lower volume of fluid is squeezed out.

#### (iii) Elastic regime

When i.e. , elastic power balances gravitational power, and at equilibrium, the first two terms on the right-hand side of equation (6.1) vanish, and we may integrate the equation once. For a thin bed, when *H*_{ℓ}/*δ*_{b}≪1, the characteristic strain so that (6.1) now yields . In contrast, for a thick carpet , the characteristic strain stores energy in a volume so that consistent with the classical Hertzian scaling for the indentation of an elastic half space by a rigid sphere of radius *R* under its own weight. These different scaling regimes are shown schematically in figure 4*b*.

### (b) Compressive deformation of a soft brush

Analysis of the final stages of approach becomes easier to analyse when the soft bed is isotropic; indeed a typical assay to measure the mechanical properties of poroelastic biological tissues uses the compression of a thin disc of the material as shown in figure 4*c*. Here we consider a circular block of the soft brush of length *H*_{ℓ} and radius *R*≫*H*_{ℓ} attached to a rigid lower plate of the same radius. The brush is compressed locally at the top by a rigid impermeable circular punch also of radius *R* positioned as shown. At time *t*=0, we apply a force to the centre of the upper plate. For a short time interval thereafter the bed behaves as a solid impermeable block with the response corresponding to the undrained elastic response. This is followed by equilibration of the pressure across the thin gap (since *R*≫*H*_{ℓ}) after which the reduction in the thickness of the gel is accompanied by a compensatory flow directed radially outwards. Scalings for the steady and transient deformations may be obtained by ignoring the (small) boundary layer at *r*=*R*. Denoting by *u*_{z} and *u*_{r} the *z* and *r* components of the deformation field and the pore pressure by *P*, one can write for this latter stage in the transient response, (∂*P*/∂*r*)∼*G*(∂^{2}*u*_{r}/∂*z*^{2}), (∂^{2}*u*_{z}/∂*z*^{2})∼0, and . We thus obtain *u*_{z}∼*A*(*t*)*z* where *A* evolves according to
6.2
When the moduli and the permeability are constant, *A*(*t*) exhibits exponential decay with a characteristic relaxation time reflecting the rate at which the interstitial fluid can diffuse out through the pores. In general for a polymer brush, *K*=*K*(*A*), *G*=*G*(*A*) and with the displacement of the top of the bed being *u*_{z}(*H*_{ℓ})=*AH*_{ℓ}. At steady state, the bed height attains a constant value given by *H*_{ℓ} with the elastic constants and the final steady state value of determined self-consistently.

### (c) Tangential motion of sphere parallel to a squishy bed

We now consider steady tangential motion of a rigid sphere over squishy beds as shown in figure 4*d*. A sphere of radius *R* moves with velocity *V* _{0}**e**_{x} such that its centre lies at a distance *R*+*h*_{0} from the surface of the undeformed bed of thickness *H*_{ℓ}. We restrict ourselves to the small deformation limit so that *h*_{0}/*R*≪1, (*R*/*h*_{0})^{1/2}≫1 and *H*_{ℓ}≪(*Rh*_{0})^{1/2}. Of interest here is the coupling between tangential velocity and normal lift force on the sphere.

Note that anisotropy in either the permeability or elasticity results in multiple pressure/deformation relaxation scales. Specifically, if the in-plane and transverse permeabilities are different, the time for elastic stresses to diffuse through the deformable layer is while the time for diffusion across the contact region is . Here, we have assumed that the shear modulus and compressive modulus are of the same order. One can then ask if the coupling of tangential motion to normal lift as elucidated by Skotheim & Mahadevan (2005) is still qualitatively valid.

To enable the analysis of both polymer, we study a two-dimensional problem of a circular cylinder translating past a model thin, dense anisotropic brush. To enable reduction to both polymer brushes (*α*_{c}=1 and *β*_{c}=0) as well as more general compressible disordered media such as filamentous gels (*α*_{c}<1 and *β*_{c}>0), we choose a simpler version of equations (2.2) and (2.9):
6.3
and
6.4
with the elastic and permeability constants *G*, *λ*, *α*_{c}, *β*_{c}, and now functions of position *z*. We differentiate between the axial and in-plane permeability to include the case of adsorbed brushes and structured gels. The elastic and permeability constants are to leading order equal to the values evaluated for the static, undeformed base state.

Recasting the equations (6.3) and (6.4) in a frame translating with the cylinder and scale quantities as in Skotheim & Mahadevan (2005), we find that the solution depends on three dimensionless parameters—the slenderness ratio , flow anisotropy , and the deformation . In the physically interesting limit *χ*≪1 and *Ψ*≪1 valid when the flux driven through the bed is small (compared with the total flux driven by the cylinder) and the in-plane and out of plane permeabilities are roughly the same, we find that the two parameters and determine the transient response. The lift force per unit length, *F*_{L} is given by . When *Ψ*∼*O*(1) anisotropy in the permeability begins to strongly influence the deformation and this expression for the lift force is no longer valid.

Extending these results to three dimensions, we find that the lift force on a sphere translating past the bed (for *ν*_{1}≫1), scales as with the prefactor being a function of *ν*_{1} and *ν*_{2}. Thus we find that the equation for the normal lift derived by Skotheim & Mahadevan (2005) is still qualitatively valid provided that the brush comprising the bed is dense.

## 7. Discussion

We have shown that a simple, asymptotically correct, anisotropic generalization of linear poroelasticity theory suffices to describe soft squishy beds of both ordered as well as disordered elastic microstructures. Our theory allows for a coupling between deformation of the solid skeleton and flow through the pores constituting the network in a self-consistent manner. Except in boundary layers close to interfaces that are of the order the pore size, momentum exchange between the fluid and solid constituents is adequately modelled using a modified Darcy’s law. Since the solid skeleton comprising the bed can deform, and matching of stresses at interfaces between the bed and pure fluid is easily accomplished without having to invoke ad hoc viscous terms as in previous approaches. Furthermore, we have shown how anisotropy in both elastic properties of the effective soft composite and the permeability due to the aligned fibres or bristles might be incorporated to explicitly calculate the material constants needed to close the governing equations in terms of the micro-geometry of the bed and the intrinsic material properties of the solid and liquid phases.

While linear poroelasticity does not correctly take into account boundary layers (of spatial extent comparable to the pore size) near interfaces, it does provide a consistent set of boundary conditions that enables matching of tangential and normal tractions and normal flow/pressure at the interface albeit with a jump in the liquid velocity gradient. In this spirit, our model resembles attempts to treat interfaces between a rigid porous medium and a free-fluid via a Darcy model supplemented with modified boundary conditions rather than invoking the Brinkman model. While the latter model is indeed valid in boundary layers close to porous(and rigid) and fluid interfaces as used in its original form, arbitrary use of the Brinkman term in the bulk is tantamount to double counting the viscous contributions. Finally, we emphasize that previous models that invoke the Brinkman approximation to study filamentous porous beds do not correctly account for the coupling between flow of permeating fluid and deformation of the solid skeleton. So we need a poroelastic model—as here. Note that for very dilute systems, the poroelastic model itself breaks down. In this case, the position and motion of each deforming filament in the bed has to be tracked in time and the resulting forces on the permeating fluid calculated explicitly—as for example in the study by Favier *et al.* (2009) of flow at high Reynolds number around a cylinder coated with rigid, slender filaments.

We conclude by commenting on two possible extensions of our theory. In the first, we consider the frictional interaction between contacting polymer brushes permeated with a good solvent. Studies on compressed and sheared polymer brushes (Rubinstein & Obukhov 1993; Galuschko *et al.* 2010) may be used to systematically account for the effect of additional non-hydrodynamic interactions on the collective dynamics of sliding between the brushes. Second, while we considered only passive elastic media, active poroelastic media that couple solid deformation and fluid flow to the concentration *c* of a species such as motor concentration may also be modelled using similar principles. The composite stress tensor for such a composite will include the fluid pore pressure, elastic response of the skeleton as well as additional terms incorporating the coupling of mechanics to chemistry—for example volumetric forces generated by active motors and moving fibres such as natural or artificial cilia via appropriate volumetric body force terms. Possible applications of such descriptions would be in modelling the mechanics of biological tissue (Teff *et al*. 2007) as well as the transport and dispersion properties of ciliary beds. We envision combining microscale experimental measurements of the geometry and material properties with macroscale poroelastic descriptions to yield a multi-scale characterization of such active beds.

## Acknowledgements

The authors thank the Harvard NSF MRSEC, the Kavli Center for NanoBio Science and Technology, the US National Institutes of Health R21 (HL091331-02 and the MacArthur Foundation for partial support.

## Footnotes

↵1 For microvilli in the endothelial glycocalyx, typical dimensional values are

*a*≈45 nm, ℓ≈1.5−3 μm,*ϕ*_{s}≈0.23–0.32.↵2 For instance, a fibre that is not clamped at the base will tend to rotate rather than deflect when acted on by a force at its tip.

↵3 Here

is to be interpreted as a poroelastic constant with the property that if pore pressure is increased by Δ*α**P*, and all normal stresses acting on the element are decreased (the compression is increased) by*α*Δ*P*, there is no change in strain.↵4 The correspondence between equation (2.9) and the analogous equation obtained using mixture theory is given in the electronic supplementary material.

↵5 It is interesting to compare these equations with the two-fluid equations characterizing the coupling between network stress

*σ*^{(n)}and composition in polymer blends (Brochard & deGennes 1983; Doi & Onuki 1992). The stress balance in these systems is usually written as↵6 Evaluating the terms in the stiffness tensor, we find , , , , and

*ν*_{13}=(*E*_{1}/*E*_{3})*ν*_{31}with constants*A*=1+*ν*_{12}and .↵7 Given that the bed is anchored at the base, one might argue that it makes no sense to talk of an in-plane modulus. Nonetheless, it is not zero as the anchoring shores up the bed. Visualize a strip of the carpet being squeezed or compressed by normal forces acting in the

*x*direction—it is this response we want to characterize. This mode of deformation exists because the fibres are fixed at the base.↵8 For a rigid non-porous bed, classical lubrication theory would predict that the speed of settling vanishes with the gap; in this degenerate limit, the sphere never makes contact. When the layer is poroelastic, the functional dependence of the force on height

*h*changes significantly once*h*∼*O*(*H*_{ℓ}).

- Received April 30, 2010.
- Accepted December 2, 2010.

- This journal is © 2011 The Royal Society