3.1. Macroscale fluid velocity

To understand the properties of these hairy channels, I first adopt a macroscale view and model the hair array as an isotropic porous medium of permeability $k$ partially obstructing the channel, as shown in figure 3(a). Several approaches exist to model such flows and the conditions at the porous–free interface (Le Bars & Worster Reference Le Bars and Worster2006). For simplicity, I use the Beavers & Joseph (Reference Beavers and Joseph1967) formulation that assumes a uniform Darcy flow $u_D=-(k/\eta )\boldsymbol{\nabla }P$ in the porous domain and a Poiseuille flow $u_P$ with a slip length of size $\sqrt {k}/\alpha$ in the free domain of height $h$ , with $\alpha$ an empirical coefficient of order $1$ characterising the porous structure

(3.1) \begin{equation} u_P(y)=-\frac {\boldsymbol{\nabla }P}{2\eta } \left [h^2-(y-(H-h))^2-\frac {\alpha \left (h^2-2k\right )}{\sqrt {k}+\alpha h}\left (H-y\right )\right ]. \end{equation}

Here, $y$ is the vertical direction whose origin is at the base of the hairs (see figure 1 d). The channel hydraulic resistance (per unit width and length) $R=-\boldsymbol{\nabla }P/q$ can then be computed by integrating the flow profile $q=\int _{0}^{H-h} u_D \:{\textrm{d}}y+ \int _{H-h}^H u_P \:{\textrm{d}}y$

(3.2) \begin{equation} R(h)=\frac {12\eta \left (\sqrt {k}+\alpha h\right )}{\alpha h^4+4\sqrt {k}h^3+6\alpha k h (2H-h)+12(H-h)k^{3/2}}. \end{equation}

Figure 3.

(a) Schematic of flow modelling at the macroscale. (b) Schematic of hair mechanics modelling at the microscale. (c) Numerical solution of the fluid–structure interaction model, for the experiment in figure 1(d). (df) Normalised horizontal and vertical tip deflections $({\textrm{d}}x/L,{\textrm{d}}y/L)$ (d,e), and normalised hydraulic resistance $\hat {R}/\hat {R}_0$ (f) as a function of the dimensionless drag force $\hat {f}_0$ for all the experiments (see Table S1). The parameters are coded through the marker inside and outside colours, shape and size; see legend figure 2. The black curves are the result of the model for the same parameters.

Equation (3.2) can be readily used to predict the pressure drop for rigid hair beds for which $h=H{-}L$ , provided the permeability $k={\textrm {f}}(d,\delta )$ is known. For regular arrays of cylinders, the permeability can be computed through asymptotic models (Happel Reference Happel1959; Drummond & Tahir Reference Drummond and Tahir1984), scaling laws (Sobera & Kleijn Reference Sobera and Kleijn2006) or homogenisation methods (Zampogna & Bottaro Reference Zampogna and Bottaro2016). For its simplicity and validity over a wide range of packing fractions, I use the scaling $k=c\delta ^2 (1{-}d/\delta )^3$ with $c=0.1475$ (Sobera & Kleijn Reference Sobera and Kleijn2006). Figure 2(c,f,i,l) shows the pressure drop predicted by (3.2) with $\alpha =1$ as dashed black lines. It captures reasonably well the pressure drop measured for ‘rigid’ hairs (see figure 2 c), as well as the initial slope of the pressure–flow rate relationship of soft hairs. Note that figure 2(c) shows three curves because the channel height $H$ varies slightly due to the roof installation procedure, resulting in a deviation of $\pm 0.2$ mm from the target (see measured values in Table S1). I also compare the predicted flow profiles, $u_D$ and $u_P$ , with simulations of the flow through rigid hair beds in the Supplementary material and find a reasonable agreement. I therefore set $\alpha =1$ in the following, which is equivalent to setting the slip length at the porous interface to $\sqrt {k}$ .

3.2. Microscale fibre deformation

To compute the hydraulic resistance of soft hair beds, I need to predict the evolution of the gap $h$ which results from the deformation of individual hairs at the microscale (see figure 3 b). Linear beam theory does not allow vertical deformations, and I thus turn to nonlinear beam theory. Denoting the coordinates of the hair centreline as $(x_h(s),y_h(s))$ , with $s$ the arclength, and $\theta (s)$ the tangent angle with respect to the vertical, the force and moment balances on a hair read

(3.3) \begin{equation} \begin{gathered} EI \theta ''(s)+\left (F(h)+f(h)(L-s)\right )\cos \theta (s)=0,\\ y^{\prime}_{h}(s)=\cos \theta (s),\qquad x^{\prime}_h(s)=\sin \theta (s),\\ x_h(0)=0,\; y_h(0)=0,\; \theta (0)=0,\; \theta ^{\prime}(L)=0. \end{gathered} \end{equation}

Here, $(\boldsymbol{\cdot })'={\textrm{d}}/{\textrm{d}}s$ , $I=\pi d^4/64$ , $f(h)$ is a distributed load coming from the fluid drag inside the hair bed, and $F(h)$ is a concentrated load at the tip from the fluid shear stress in the free region. Both forces a priori depend on the gap $h$ and thus the beam deflection itself since $h=H-y_h(L)=H{-}\int _0^L \cos (\theta (s))\:{\textrm{d}}s$ .

To compute these microscale forces, the averaged quantities computed in § 3.1 are not adequate and one needs to consider the details of the flow at the hair scale. Inside the hair bed, studies on two-dimensional cylinder arrays (Happel Reference Happel1959; Drummond & Tahir Reference Drummond and Tahir1984; Sobera & Kleijn Reference Sobera and Kleijn2006) have shown that the fluid drag coming from the microscale flow is related to the permeability through momentum conservation

(3.4) \begin{equation} f=\frac {\delta ^2}{k}\eta u_D=\delta ^2 q R(h). \end{equation}

At the hair tip, the wall shear stress $\tau _s=\eta \left . ({{\textrm{d}}u}/{{\textrm{d}}y}) \right \vert _{H-h}$ integrated over the tip surface generates the force $F=\int \tau _s \:{\textrm{d}}S$ . I am not aware of previous studies calculating the local fluid velocity in this configuration. However, it is clear that using the averaged macroscale velocity $u_P$ which is non-zero at the hair surface is inadequate. Guided by the numerical simulations of the flow through rigid hair beds shown in the Supplementary material, I propose to approximate the local velocity profile as parabolic, $u(y)= ({6\bar {u}}/{h^2}) (H{-}y ) (y{-}(H{-}h))$ , with $\bar {u}$ the local average flow velocity. I then use the macromodel to estimate this average local velocity $\bar {u}= {1}/{h}\int _{H{-}h}^H u_P \:{\textrm{d}}y$ . Using this approximate local velocity profile to compute the shear stress and integrating it over the hair tip yields the tip force

(3.5) \begin{equation} F=\frac {\pi d^2}{4}\eta \left . \frac {{\textrm{d}}u}{{\textrm{d}}y} \right \vert _{H{-}h}=\frac {\pi }{8} d^2 q R(h) \frac {\alpha h^2+4\sqrt {k}h+6\alpha k}{\alpha h +\sqrt {k}}. \end{equation}

3.3. Non-dimensionalisation

Combining (3.2)–(3.5) yields a reduced-order, fully coupled, fluid–structure interaction model that has the form of a one-dimensional integro-differential boundary value problem. This model is made dimensionless by rescaling all lengths by $L$ in the beam equation, (3.3)

(3.6) \begin{equation} \begin{gathered} \theta ''(\hat {s})+\left (\hat {F}+\hat {f}(1-\hat {s})\right )\cos \theta (\hat {s})=0,\\ \hat {y_h}'(\hat {s})=\cos \theta (\hat {s}),\qquad \hat {x_h}'(\hat {s})=\sin \theta (\hat {s}),\\ \hat {x_h}(0)=0,\; \hat {y_h}(0)=0,\; \theta (0)=0,\; \theta '(1)=0. \end{gathered} \end{equation}

Here, $\hat {s}=s/L$ , $\hat {x_h}=x_h/L$ , $\hat {y_h}=y_h/L$ and the dimensionless fluid forces are $\hat {F}=FL^2/EI$ , $\hat {f}=fL^3/EI$ . Using the relevant aspect ratios, $\hat {L}=L/H$ , $\hat {d}=d/\delta$ , $\hat {\delta }=\delta /H$ , the dimensionless permeability reads $\hat {k}=k/H^2=c\hat {\delta }^2(1-\hat {d})^{3/4}$ and the dimensionless gap $\hat {h}=h/H=1-\hat {L}\int _0^1\cos \theta \:{\textrm{d}}\hat {s}=1-\hat {L}\hat {y_h}(1)$ . Using these along with the elasto-viscous parameter $\hat {q}=\eta q/E d^2$ yields the dimensionless forces

(3.7) \begin{equation} \hat {f}=\frac {768}{\pi }\frac {\hat {L}^3\hat {q}}{\hat {d}^2}\hat {R}, \quad \hat {F}=96\hat {L}^2\hat {q}\hat {R}\frac {\alpha (1-\hat {L}\hat {y_h}(1))^2+4\sqrt {\hat {k}}(1-\hat {L}\hat {y_h}(1))+6\alpha \hat {k}}{\alpha (1-\hat {L}\hat {y_h}(1))+\sqrt {\hat {k}}},\\ \end{equation}

with $\hat {R}=RH^3/12\eta$ the hydraulic resistance rescaled by that of an empty channel

(3.8) \begin{align} \hat {R}=\frac {\sqrt {\hat {k}}+\alpha (1-\hat {L}\hat {y_h}(1))}{\alpha (1-\hat {L}\hat {y_h}(1))^4+4\sqrt {\hat {k}}(1-\hat {L}\hat {y_h}(1))^3+6\alpha \hat {k} (1-(\hat {L}\hat {y_h}(1))^2)+12\hat {L}\hat {y_h}(1)\hat {k}^{3/2}}.\\[-18pt]\nonumber \end{align}

I solve the dimensionless equations, (3.6), (3.7), (3.8) numerically with the Matlab solver bvp5c and show in figure 3(c) the hair profile as the flow rate is increased for the experimental parameters of figure 1(d). In figure 2, I compare the tip deflection predicted by the model with the one measured experimentally. I find a quantitative agreement for most experimental conditions without any adjustable parameters. In panels (g,h), the permeability spans an order of magnitude, and the quantitative discrepancies in the model are likely due to inaccuracies in the scaling used to calculate $k$ . Despite these limits, the model captures the surprisingly weak influence of the hair spacing $\delta$ on the hair deflection seen in these experiments. Using the computed tip deflection for $h$ in (3.2) yields a prediction for the pressure drop across the channel which is compared with the experimental pressure drop in figure 2(c,f,i,l). Here also, the agreement is excellent.

3.4. Asymptotic expansion for small flow rates

To get analytical insights on the role of the various length scales in the problem, I first notice that the tip force $\hat {F}$ is usually negligible compared with the drag inside the hair bed $\hat {f}$

(3.9) \begin{equation} \frac {\hat {f}}{\hat {F}}=\frac {8}{\pi }\frac {\hat {L}}{\hat {d}^2}\frac {\alpha (1-\hat {L}\hat {y_h}(1))+\sqrt {\hat {k}}}{\alpha (1-\hat {L}\hat {y_h}(1))^2+4\sqrt {\hat {k}}(1-\hat {L}\hat {y_h}(1))+6\alpha \hat {k}}\gg 1. \end{equation}

Since $\hat {q}\lt 2.3\times10^{-3}$ in experiments, I then introduce the asymptotic expansion $\theta (\hat {s})=\sum _n \hat {q}^n \theta _n(\hat {s})$ , with $\hat {q}$ treated as a small parameter (Alvarado et al. Reference Alvarado, Comtet, de Langre and Hosoi2017) in (3.6) with $\hat {F}=0$ . Restricting myself to the case of straight hairs for simplicity, the solution should be symmetric with respect to the flow direction and even terms are thus zero. In particular, at order zero $\theta _0(\hat {s})=0$ . At the first order in $\hat {q}$ , $\theta (\hat {s})= \hat {q} \theta _1(\hat {s})$ and thus $\cos \theta (\hat {s})=\cos (\hat {q} \theta _1(\hat {s}))= 1$ , and $\hat {y_h}(1)=\int _0^1\cos \theta (\hat {s})\:{\textrm{d}}\hat {s}= 1$ . Inserting this in the beam equation without tip force yields at first order

(3.10) \begin{equation} \hat {q} \theta ^{\prime \prime}_{1}(\hat {s})+\hat {f}_0(1-\hat {s})=0, \end{equation}

with

(3.11) \begin{equation} \hat {f}_0=\frac {768}{\pi }\frac {\hat {L}^3\hat {q}}{\hat {d}^2}\hat {R}_0,\quad \hat {R}_{0}=\frac {\sqrt {\hat {k}}+\alpha (1-\hat {L})}{\alpha (1-\hat {L})^4+4\sqrt {\hat {k}}(1-\hat {L})^3+6\alpha \hat {k} (1-\hat {L}^2)+12\hat {L}\hat {k}^\frac {3}{2}}. \end{equation}

Here, $\hat {R}_0$ is the dimensionless hydraulic resistance in the absence of deformations, i.e. for $\hat {y_h}(1)=1$ . Equation (3.10) is mathematically equivalent to a linear beam model under a uniform distributed load whose solution is

(3.12) \begin{equation} \theta _1(\hat {s})=\frac {\hat {f}_0}{6 \hat {q}}\left (\hat {s}^3-3\hat {s}^2+\hat {s}\right ). \end{equation}

Since $\theta _2(\hat {s})=0$ by symmetry, $\theta (\hat {s})= \hat {q} \theta _1(\hat {s})$ with $\theta _1(\hat {s})$ from (3.12) is valid up to second order in $\hat {q}$ .

This analysis suggests that the four dimensionless parameters $\hat {q}$ , $\hat {L}$ , $\hat {\delta }$ and $\hat {d}$ can be combined in a single one: $\hat {f}_0$ (see (3.11)). This small $\hat {q}$ /small slope solution can then be used to compute the tip deflection

(3.13) \begin{equation} \begin{aligned} \hat {x_h}(1) &=\int _0^1 \sin \theta \:{\textrm{d}}\hat {s}\approx \int _0^1 \hat {q} \theta _1 \:{\textrm{d}}\hat {s}\approx \frac {\hat {f}_0}{8}, \\ \hat {y_h}(1) &=\int _0^1 \cos \theta \:{\textrm{d}}\hat {s}\approx \int _0^1 \left (1-\frac {(\hat {q} \theta _1)^2}{2}\right ) \:{\textrm{d}}\hat {s}\approx 1- \frac {\hat {f}^2_0}{112}. \end{aligned} \end{equation}

Figure 3(d,e) shows the tip deflection for all the experiments and the model plotted as a function of $\hat {f}_0$ . Surprisingly, the data not only collapse for small deformations, as expected, but also for larger deformations where the asymptotic expansion is expected to fail. The hydraulic resistance $\hat {R}$ , computed from the measured pressure loss $\Delta P$ for experiments and calculated from the tip deflection using (3.8) for the model, is rescaled by its value for rigid hairs $\hat {R}_0$ and shown as a function of $\hat {f}_0$ in figure 3( f). Here also the data collapse, even for large deflection. The hair flexibility starts to impact the hydraulic resistance when $\hat {f}_0\sim 1$ , and I observe empirically $\hat {R}/\hat {R}_{0}\sim \hat {f}^{-2/3}_0$ for large deformations. The rigid hair bed drag force $\hat {f}_0$ is thus an efficient control parameter encompassing both the material and geometrical properties of the system and can be used to predict the hydraulic resistance of channels obstructed by dense hair-like structures.