Dense brushes of stiff polymers or filaments in fluid flow

Polymer brushes have been subject to enormous interest in the last several decades due to their importance in various technological and biological systems [1–3]. However, the main focus of most experimental, theoretical, and simulation studies has been on the behavior of brushes consisting of grafted flexible polymers (i.e., with negligible flexural stiffness on the length of a polymer). Brushes with stiff chains, whose persistence length is larger or comparable to the contour length, have received limited attention to date, in spite of their frequent presence in many biophysical systems. Dense filamentous brush-like structures often serve as structural elements in the human body with the examples of lubricating aggrecan brushes in joints [4], endothelial glycocalyx layer in blood vessels [5,6], periciliary layer of lung airway [7], and hair cells in the inner ear mediating the sense of hearing and balance [8,9]. Better understanding of their non-equilibrium behavior, for instance their response to flow, will lead to insights into the functionality of the respective structures and the effects of pathological alterations. Technologically, brushes with stiff grafted elements, such as cantilevered micro- and nano-rod arrays [10] or high-density brushes [11], can be used in nanofluidic devices [12] or in the context of functionalized surfaces [13].

Equilibrium properties of stiff-polymer brushes have been investigated in a few simulation [14–16] and theoretical [16–18] studies predicting the brush height. In silico compression tests of dense semiflexible-polymer brushes [19,20] have shown that such brushes might be much softer mechanically than the brushes made out of flexible polymers making them useful in micro-flow applications. As an example, the response of a glycocalyx layer in blood vessels to fluid flow has been studied in simulations [21] predicting a flow-rate–dependent increase in the flow resistance. Recently, a mean-field approach [22] has been developed to describe the flow-induced deformation of grafted semiflexible polymers. In this theory, the steric interactions between polymers have been neglected limiting the model's applicability to the deformation of a single grafted polymer.

In this letter, we develop a theoretical model which is able to describe quantitatively the flow-induced deformation of a dense stiff-polymer brush and its effect on flow. The model is applicable to the brushes where individual grafted polymers or filaments possess a persistence length which is larger or comparable to their total contour length. To validate this model we perform corresponding mesoscopic simulations for a wide range of conditions including polymer elasticity, grafting density, and flow rate. The model allows us to identify the importance of different contributions to brush deformation such as hydrodynamic friction due to flow and steric excluded-volume interactions between polymers. This model can be employed for the quantification of experiments in biological settings and for the development of brush interfaces with specific surface properties. It can be also used to describe the behavior of a stiff-polymer brush in flow, where explicit brush modeling is difficult or unfeasible.

To illustrate the system of interest, fig. 1 presents the simulation setup, where a slit-like geometry with a height D is employed. Both the top and bottom walls are covered by stiff polymers grafted on the square lattice with a constant $a_{\mathrm{lat}}$. The polymers are built up from N = 11 bonded beads with a diameter $d=1.0$ and the potential energy given by

$$\begin{equation} U=\sum_{i=1}^{N-1}\left[\frac{k_s}{2d}\left[r_{i,i+1}-d\right]^2+\frac{EI}{d}\left[1-\cos \theta_i\right]\right], \end{equation} \tag{ 1 }$$

where k_s is the spring constant, $r_{i,i+1}$ is the distance between the beads i and i + 1, EI is the bending rigidity, $\theta_i$ is the angle between two consecutive springs, and $EI/k_s=d^2/16$. The polymers are stiffly anchored at the surface by fixing the bead i = 1 at the surface and the bead i = 0 inside the surface. To prevent overlap between beads, the purely repulsive, shifted, and truncated at $r_c=d$ Lennard-Jones (LJ) potential [23] is introduced with the parameters $\epsilon=k_BT$ and $2^{1/6}\sigma_{LJ}=d$, where k_BT is the thermal energy unit.

Zoom In Zoom Out Reset image size

To model the fluid flow, we employ the smoothed dissipative particle dynamics (SDPD) method [24] using a density of $n_s=3$ for both fluid and wall particles. Interactions of the fluid particles with the polymer beads are mediated by friction forces. Shear flow is generated by moving one of the walls with a constant velocity. The Poiseuille flow is driven by an external force f_q which acts on each fluid particle in the x-direction. This corresponds to a pressure drop $\Delta P$ along the slit length l_x such that $f_q n_s=\Delta P/l_x$.

To characterize the simulated systems we use several dimensionless parameters such as the bending rigidity of the grafted polymers,

$$\begin{equation} \frac{l_p}{L}=\frac{EI}{k_BTL}, \end{equation} \tag{ 2 }$$

grafting density $\sigma=(d/a_{\mathrm{lat}})^2$, and non-dimensional shear rate,

$$\begin{equation} \tilde{\dot{\gamma}}=\frac{\eta L^3}{k_BT}\dot{\gamma}. \end{equation} \tag{ 3 }$$

Here, $l_p/L$ is the ratio of polymer persistence and contour lengths, η is the fluid's viscosity, and $\dot{\gamma}$ is the shear rate on top of a brush. The effective brush height h is calculated using the first moment of the polymer-bead density profile similar to that in ref. [21],

$$\begin{equation} h=2\frac{\int y\rho(y)\text{d}y}{\int\rho(y)\text{d}y}. \end{equation} \tag{ 4 }$$

To verify our simulation model, we compare in fig. 2 the SDPD results with previous simulations of grafted semiflexible polymers in shear flow using the Brownian-dynamics and lattice-Boltzmann methods [22]. The SDPD results of this work are in excellent agreement with the previous simulations of similar brush systems [22]. In addition, we have also tested whether the type of shear flow (e.g., Couette or Poiseuille) over a brush may have a considerable effect. Figure 2 shows that brush deformation appears to be nearly independent of the flow applied and can be well characterized by $\tilde{\dot{\gamma}}$. To span a wide range of conditions, we have simulated systems with grafting densities σ in the range from 0.01 to 1 (which corresponds to SC close packing), polymer elasticities $l_p/L$ between 10 and 100, and shear rates $\tilde{\dot{\gamma}}$ between 10¹ and 10⁶. Note that the grafted polymers are relatively stiff with $l_p/L>1$, and, therefore, entropic conformational changes such as polymer coiling can be practically neglected. The corresponding Reynolds numbers $Re=vLn_s/\eta$, with v being the velocity at the brush top, have been kept between $10^{-6}$ and $10^{-1}$ to ensure no inertial effects. Monitoring of bond lengths has shown that no significant extension $(\lesssim1{\%})$ of polymers has been present, and thus, polymer stretching can be neglected.

Zoom In Zoom Out Reset image size

Our mean-field approach to describe the behavior of stiff polymer brushes in shear flow assumes identical deformation for all polymers. Thus, we consider a single polymer as an elastic cantilever which is subject to a hydrodynamic drag force due to fluid flow and steric interactions due to its neighboring polymers in the lattice. The polymer deformation is described by its internal curve-linear coordinate $s \in[0,L]$ and the angle $\theta(s)$ between the y-axis and a local tangent line at s. The correspondence of internal coordinates to the Cartesian coordinates is simply done as $x(s)=\int_{0}^{s}\sin\theta(s') \,\text{d}s'$ and $y(s)=\int_{0}^{s}\cos\theta(s')\,\text{d}s'$.

The torque balance for a circular rod [25], where external forces $(\textbf{F}(s)=(F_x(s);~F_y(s)))$ are counteracted by the rods's elastic resistance, leads to

$$\begin{equation} EI\frac{\text{d}\textbf{r}(s)}{\text{d}s}\times\frac{\text{d}^3\textbf{r}(s)}{\text{d}s^3}= \textbf{F}(s)\times\frac{\text{d}\textbf{r}(s)}{\text{d}s}, \end{equation} \tag{ 5 }$$

where $\textbf{r}(s)$ is the position on a rod at s such that $\text{d}\textbf{r}(s)/ \text{d}s=(\sin\theta(s);~\cos\theta(s))$. This equation formulated in terms of the local angle $\theta(s)$ reduces to

$$\begin{equation} EI\frac{\text{d}^2\theta(s)}{\text{d}s^2}=F_y(s)\sin\theta(s)-F_x(s)\cos\theta(s), \end{equation} \tag{ 6 }$$

with the boundary conditions of $\theta(s)|_{s=0}=\theta_0$ at the grafted surface and $\text{d}\theta/\text{d}s|_{s=L}=0$ at the free end. $\theta_0$ characterizes the grafting angle between the polymer beam and the surface's normal and can generally be a constant or a function of stress applied to the beam. However, in the current study we investigate perpendicular grafting of stiff filaments to a surface with $\theta_0=0$.

The force $\textbf{F}(s)$ consists of two contributions: i) drag force $(\textbf{F}^{\mathrm{d}}(s))$ from fluid flow and ii) a force due to excluded-volume interactions $(\textbf{F}^{\mathrm{v}}(s))$, and thus, $\textbf{F}(s)=\textbf{F} ^{\mathrm{d}}(s)+\textbf{F}^{\mathrm{v}}(s)$. Note that beam stretching is neglected. The force at a position s can be found as an integral over local force density $(\textbf{f} (s)=\textbf{f}\,^{\mathrm{d}}(s)+\textbf{f}^{\mathrm{v}}(s))$ such that

$$\begin{equation} \textbf{F}(s)=\int_{s}^{L}\textbf{f}(s')\text{d}s'. \end{equation} \tag{ 7 }$$

The force density of the drag force exerted on a stiff polymer is given by

$$\begin{eqnarray} \begin{array}{rcl} f^{\mathrm{d}}_x(s)&=&f^{\mathrm{d}}_{\perp}(s)\cos\theta(s)+f^{\mathrm{d}}_{\parallel} (s)\sin\theta(s)=\\ &&\hspace{-16pt}\left(\zeta_{\perp}\cos^2\theta(s)\!+\!\zeta_{\parallel}\sin^2\theta(s)\right)\eta u_x(s)\!=\! \zeta(s)\eta u_x(s), \\ f^{\mathrm{d}}_y(s)&=&f^{\mathrm{d}}_{\perp}(s)\sin\theta(s)-f^{\mathrm{d}}_{\parallel}(s) \cos\theta(s)= \\ &&\hspace{-16pt}\eta u_x(s)\sin\theta(s)\cos\theta(s)\left(\zeta_{\perp}-\zeta_{\parallel} \right), \end{array} \quad\end{eqnarray} \tag{ 8 }$$

where $f^{\mathrm{d}}_{\perp}(s)$ and $f^{\mathrm{d}}_{\parallel}(s)$ are the corresponding normal and parallel components of the force density with the convention that $f^{\mathrm{d}} _{\perp}(s)$ is directed with the x-axis and down to the wall (see fig. 1), while $f^{\mathrm{d}}_{\parallel}(s)$ has a tangential direction toward the increase of the s-coordinate. $\zeta_{\perp}$ and $\zeta_{\parallel}$ are the normal and tangential components of the friction coefficient per unit length, $\zeta(s)= \zeta_{\parallel}\sin^2\theta(s)+\zeta_{\perp}\cos^2\theta(s)$, and $u_x(s)$ is the local flow velocity and $u_y(s)=0$. The friction coefficients are approximated using the slender body theory [26] for a thin cylinder similar to that in ref. [27],

$$\begin{equation} \begin{array}{rcl} \zeta^{\perp}&=&\displaystyle\frac{8\pi}{\ln(L/d)-\frac{1}{2}+\ln(2)}, \zeta^{\parallel}&=&\displaystyle\frac{4\pi}{\ln(L/d)-\frac{3}{2}+\ln(2)}. \end{array} \end{equation} \tag{ 9 }$$

Note that these friction coefficients describe a drag on an isolated rod. If we neglect hydrodynamic interactions, the drag force on a polymer can be described using a single friction coefficient $\zeta=\zeta_{\parallel}=\zeta_{\perp}=3\pi$ as in ref. [22]. In this case, the force density $f^{\mathrm{d}}_y(s)$ vanishes. The differences in model predictions using these two limiting cases for the friction coefficients will be discussed further in text.

The local velocity $u_x(s)$ depends on hydrodynamic penetration into the brush [28] and is described by the Brinkman equation [29] for flow in porous media as

$$\begin{equation} \frac{\text{d}^2u_x(y)}{\text{d}y^2}=\zeta(y)\frac{\sigma}{d^2}\frac{u_x(y)}{\cos\theta(y)}. \end{equation} \tag{ 10 }$$

The same equation has been also used in the theory by Kim et al. [22]. Boundary conditions for eq. (10) are $u_x(y)|_{y=0}=0$ and $\text{d}u_x(y)/\text{d}y|_{y_{\textit{max}}}=\dot{\gamma}|_{s=L}$. Equation (10) in internal coordinates becomes

$$\begin{equation} \frac{\text{d}^2u_x(s)}{\text{d}s^2}=\zeta(s)u_x(s)\frac{\sigma}{d^2}\cos\theta(s)- \frac{\text{d}u_x(s)}{\text{d}s}\frac{\text{d}\theta(s)}{\text{d}s}\tan\theta(s), \end{equation} \tag{ 11 }$$

with $u_x(s)|_{s=0}=0$ and $\text{d}u_x(s)/\text{d}s|_{s=L}=\dot{\gamma}|_{s=L}\cos \theta(L)$ being boundary conditions.

To introduce the steric interactions between the filaments, we discretize the beam into $N=\left\lfloor L/d\right\rfloor$ spheres similar to the polymer representation in simulations, where the symbol $\left\lfloor*\right\rfloor$ denotes the integer floor function. Since identical deformation of all polymers is assumed, only steric interactions in the plane of beam deformation need to be considered. Therefore, the calculation of volume exclusion interactions includes only the two neighbors surrounding a beam in the flow direction. Finally, the force density $\textbf{f}^{\mathrm{v}}_n$ on a sphere $n\in [0,\ldots,N-1]$ in a discretized beam due to excluded-volume interactions is computed by a sum over all spheres j of neighboring polymers within the plane of beam deformation as

$$\begin{equation} \textbf{f}^{\mathrm{v}}_n=\sum_j\frac{\textbf{g}_{nj}}{d}, \end{equation} \tag{ 12 }$$

where $\textbf{g}_{nj}$ is the force between spheres n and j approximated by a simple repulsive force,

$$\begin{equation} \textbf{g}_{ij}=\left\{\begin{array}{l@{\quad}c} r_{ij} \leqslant d\!: &\displaystyle\epsilon_g \frac{k_BT}{d} \left[ \left( d/r_{ij} \right) ^{\alpha}-1\right] \frac{\textbf{r}_{ij}}{r_{ij}}\\ r_{ij} > d\!: &0 \end{array}\right\}, \end{equation} \tag{ 13 }$$

where $r_{ij}=|\textbf{r}_{ij}|$, and $\epsilon_g$ and α are parameters controlling the softness of inter-bead interactions. This force definition appears to be very robust, and the results are hardly affected by $\epsilon_g$ and α, if $\epsilon_g \gtrsim 50$ and $\alpha \geq 3$, since the interactions become hard enough. Thus, in order to have a larger time step for a numerical solver, we use $\epsilon_g = 50$ and $\alpha = 3$ in the theoretical model, while in simulations the employed LJ potential is harder. The volume exclusion force $\textbf{F}^{\mathrm{v}}(s)$ is calculated through the integral in eq. (7), whose discrete representation is given by

$$\begin{equation} \textbf{F}^{\mathrm{v}}(s)=\textbf{f}^{\mathrm{v}}_{\left\lfloor s/d\right\rfloor}\left( d\left\lceil s/d\right\rceil-s/d\right)+\sum_{n=\left\lceil s/d\right\rceil}^{N-1} \textbf{f}^{\mathrm{v}}_nd, \end{equation} \tag{ 14 }$$

where the symbol $\left\lceil*\right\rceil$ denotes the integer ceiling function.

Direct numerical solving of eq. (6) appeared to be difficult for high grafting densities and shear rates, i.e. in the regime where excluded-volume interactions dominate. Therefore, we substituted eq. (6) with its time-dependent version given by

$$\begin{equation} \frac{\partial\theta}{\partial t}=EI\frac{\partial^2\theta}{\partial s^2}-F_y(s)\sin\theta +F_x(s)\cos\theta, \end{equation} \tag{ 15 }$$

where θ becomes a function of s and time t. Finally, eqs. (15) and (11) are solved numerically for a given set of parameters (elasticity, grafting density, and shear rate on top of a brush) using the following iterative procedure. i) Initially, we guess a beam configuration and the polymer is discretized into 10² elements. Then, eq. (11) is solved using a boundary value problem (BVP) solver for ordinary differential equations [30] in order to obtain the velocity profile. ii) At this step, the calculation of the hydrodynamic drag force $\textbf{F}^{\mathrm{d}}(s)$ and the volume-exclusion interaction force $\textbf{F} ^{\mathrm{v}}(s)$ is performed following eqs. (8) and (14), respectively. iii) As a next step, eq. (15) is integrated using the DuFort-Frankel scheme [31], and the resulting configuration is taken as an initial guess for a new cycle. This iteration is repeated until the convergence condition $\sum_{i}\left\|\theta_i(t)-\theta_i(t-\Delta t)\right\| \leq10^{-7}$ is satisfied, where $\Delta t$ is the time step.

Figure 3(a) presents our main result, the relative brush height as a function of the shear rate on top of the brush normalized by the bending rigidity of the polymers for the case of two friction coefficients from eq. (9). The predictions of our theoretical model are in good quantitative agreement with the corresponding SDPD simulations for two different bending rigidities of $l_p/L=10$ and 100. Brush deformation for different elasticities of the grafted polymers shows a universal behavior with respect to the polymer bending rigidity, if polymers are stiff enough with $l_p/L\gtrsim10$. Clearly, the theory is expected to fail when EI/L becomes smaller or comparable to k_BT (or if $l_p/L\lesssim1$) or when the grafted polymers can be considered rather flexible. In this case entropic effects have to be necessarily included, which is out of the scope of this work.

Zoom In Zoom Out Reset image size

Some deviations between theoretical predictions and simulations are observed at high grafting densities ${(\sigma\gtrsim0.5)}$ and flow rates due to a better packing of deformed stiff polymers. Remember that the theory assumes identical deformation for all fibers such that possible deformation in the flow vorticity direction (z-axis) is not considered. Thus, the theory overpredicts the effect of excluded-volume interactions, which can be seen in fig. 3(a) for $\sigma=1.0$ where the theory estimates a larger brush height in comparison to that obtained in simulations. The curve for $\sigma=1.0$ also nicely illustrates the onset of excluded-volume effects with respect to shear rate by a significant change in the slope of brush height occurring at approximately $EI\frac{ \text{d}^2\theta}{\text{d}s^2}\sim F^{\mathrm{v}}_y\sin\theta-F^{\mathrm{v}}_x \cos\theta$. Thus, at large shear rates the brush height is mainly determined by the excluded-volume interactions.

To further support the validity of the theoretical model, we present the correspondence of fluid velocities in fig. 3(b). An excellent agreement of theoretical and simulated velocity profiles indicates that the approximation using the Brinkman equation (eq. (11)) properly captures the frictional effect of the brush on fluid flow for a wide range of grafting densities. However, we need to mention that at the lowest grafting density studied $(\sigma=0.01)$, the Brinkman equation overestimates the influence of the brush on flow. Therefore, in this single case we simply assumed a constant shear rate $\frac{\text{d}u}{\text{d}y}={\dot{\gamma}}|_{y=h}$ within the brush, which matches very well the simulated velocity profile, as shown in fig. 3(b).

The model predictions in fig. 3 correspond to the case of two friction coefficients from eq. (9), which were derived for a single isolated rod [26]. Thus, they are expected to provide a reasonable approximation for the friction on filaments at low grafting densities. As the grafting density is increased, hydrodynamic interactions are expected to get screened leading to the loss of hydrodynamic correlations between different segments of a polymer. In this case, a single friction coefficient $\zeta=\zeta_{\parallel}=\zeta_{\perp}=3\pi$ can be assumed. Figure 4 presents the comparison of model predictions using the two different choices of friction coefficients for various grafting densities and $l_p/L=10$. As expected, the friction coefficient of $3\pi$ leads to an overestimation of the applied drag for low grafting densities, and the other choice (eq. (9)) appears to be better. At high grafting densities, the difference in model predictions using different friction coefficients nearly vanishes, even though the choice of $\zeta=3\pi$ has been expected to be the best. This indicates that at high grafting densities the brush height is mainly determined by the excluded-volume interactions and moderate changes in the applied friction play a secondary role.

Zoom In Zoom Out Reset image size

The effect of stiff-polymer brush on the fluid flow can be characterized by an increase of the flow resistance due to the presence of a brush. Figure 5 shows that the relative apparent viscosity, the ratio of the apparent viscosity (obtained by fitting the Poiseuille law to the resulting flow rate) to the fluid viscosity, might strongly increase if a brush is hardly deformed, and becomes smaller as the brush gets bent by the flow. Qualitatively, this effect can be understood by a change in an effective channel diameter. The flow-induced brush deformation can be used for a flow control in microfluidics and is directly related to an increase in blood flow resistance in small vessels whose surface is covered by glycocalyx having a brush-like structure [5,6]. In vivo experiments on blood flow resistance [32] reveal a much higher resistance in small vessels $(D\lesssim35\ \mu \text{m})$ in comparison to in vitro experiments in glass tubes [33]. This effect is mainly attributed to the glycocalyx layer at vessel walls [34], and the theoretical model of this work can provide its quantitative description.

Zoom In Zoom Out Reset image size

In conclusion, the presented theoretical model is able to predict quantitatively the flow-induced deformation of a brush of stiff polymers or filaments and its effect on fluid flow. The model quantifies the relative contributions of fluid friction, polymer bending resistance, and steric excluded-volume interactions between polymers. In comparison to the analytic mean-field approach by Kim et al. [22], the presented theoretical model significantly extends the range of model applicability to high grafting densities of stiff polymers and strong deformations under flow. In particular, the current model explicitly includes the effect of excluded-volume interactions between different polymers. We expect that this model will be used for the quantification of biological and technological experiments with stiff brush-like structures, and can be extended to include direct mechanical deformations due to external forces (e.g., brush-cell interactions).

We would like to thank M. Deng and G. E. Karniadakis for useful discussions. We also gratefully acknowledge the computing time granted on the supercomputer JUROPA at Jülich Supercomputing Centre (JSC). DAF acknowledges funding by the Alexander von Humboldt Foundation.