Entropic stress of grafted polymer chains in shear flow

Polymer brushes are systems of surface-grafted polymer chains. They have been used in various applications, such as stabilizing colloidal suspensions,^1–3 anti-fouling coatings,^4,5 and adjusting surface interactions and wetting behavior.⁶ They naturally occur on the interface of complex biological or bio-inspired systems, such as cartilage⁷ or hydrogels,^8–17 and are known for their excellent lubricating properties.^18–20 

In this paper, we apply molecular simulations to study the relation between interfacial shear stress and polymer conformation in brushes exposed to solvent flow. We model systems of grafted chains with varying chain length and grafting density [see Figs. 1(a) and 1(b)] and analyze their microscopic response to shear flow. We explore the polymer’s role by computing the contribution of the conformational entropy to the shear stress in the nonequilibrium chain-conformations imposed by the external flow.^21,22

Prior simulation studies on grafted polymer chains discussed polymer conformations during flow,^23–29 but those studies did not discuss mechanical stresses in the system and focused on high grafting density. Here, we explicitly analyze stresses with the goal to connect polymer conformation to external drivers, as in frictional contacts. We also include brushes with low grafting density, such as those encountered on the surfaces of hydrogels,^15–17 in the hope that our work contributes to the ongoing debate on the roles of solvent hydrodynamics and nonequilibrium polymer conformations in the tribology of those systems.^10–17 Our work is hence most representative of a sliding polymeric system in a hydrodynamic regime.

Equilibrium chain conformations depend on grafting density. At low grafting densities, in what is called the “mushroom” regime, the chains do not interact with each other. With increasing grafting density Γ, the chains start to overlap and adopt increasingly straight, brush-like conformations due to excluded volume effects.¹⁹ The transition between both regimes occurs at Γ ≈ Γ^∗, with $Γ∗=1/(πRg2)$⁠, where R_g is the radius of gyration of the individual chain inside the brush.

When exposed to shear flow, the solvent’s velocity gradient penetrates a certain distance into the brush before it is screened out.³⁰ This momentum transport into the brush reorients and elongates chains in the direction of shear.^23–29 As a result, a grafted polymer chain in shear flow is subject to three forces: an osmotic force due to the gradient in polymer concentration across the interface, a shear force due to the solvent’s shear flow, and a restoring force due to conformational entropy.³¹ Studies have shown that the competition between these forces in nonlinear flows can drive chains to split into two coexisting populations, with one population elongating in flow, while the other remains unperturbed within the brush.^29,32–34

Experiments commonly do not have access to λ and, therefore, τ_r. The Weissenberg number $Wi∗=γ̇τr∗$ is then defined with respect to the equilibrium relaxation time of the full free chain, $τr∗$⁠.^{11,12,39–41} The key question that we are trying to answer is what the consequences of this approximation are for the observation of macroscopic properties, in particular the shear stress. Our results confirm that only the part of the brush with height λ actually exposed to solvent shear determines the brush’s shear response. The value of λ varies with shear rate and grafting density, information that is not experimentally accessible. This implies that the rate-dependence of macroscopic properties, such as the external stress, will depend on how Wi is defined. In the regime where chains elongate, the external stress appears to depend sublinearly on Wi only if a shear-rate-independent relaxation time, such as an equilibrium relaxation time, is used to compute the Weissenberg number, as is typical for experiments measuring friction in self-mated hydrogels and polymer brushes.^11,12,39,41

We can now insert into Eq. (3) the equation for the non-linear, entropic force of a freely jointed polymer chain $|f⃗|=kBTL−1(h)/lk$⁠,²² in which l_k is the Kuhn length, $L−1$ is the inverse Langevin function,³⁸ and $h=|r⃗|/nb$ is the extension ratio with a bond length b. The extension ratio is proportional to n^−1/2 for an ideal, Gaussian chain and unity for a fully extended chain. The Kuhn length is defined as l_k = C_∞b, in which C_∞ is Flory’s characteristic ratio. We note that this is an approximation for the systems that we are studying, as it assumes that each chain is under uniform tension which is not the case in most real environments. However, past simulation work has shown that this approximation works well (see, for example, Ref. 22).

This derivation ignores changes to the chain’s internal energy that might arise at large extensions and, hence, only captures the effects of changes in configurational entropy, which is usually the dominant source of stress in shearing flows. For the purpose of clarity and brevity, we henceforth call the polymer’s contribution to the stress tensor the entropic stress $σαβentr$⁠. Two different molecular effects go into the entropic stress. First, the hardening factor $L−1hh$ captures the mechanical effect of chain elongation. Second, the orientation factor $r̂αr̂β$ captures the alignment of polymer chains.

Our molecular system, depicted in Figs. 1(a) and 1(b), consists of polymer chains bonded to a face-centered cubic wall with a lattice spacing of 1.2 σ and the [111] crystallographic direction facing the interface. The chains are bonded to the exposed wall beads following a hexagonal pattern with grafting density Γ. The polymers are described by a coarse-grained bead-spring model⁴² in which the excluded volume is modeled by a Lennard-Jones (LJ) potential with interaction energy ɛ = 1.0, length σ = 1.0, and cutoff radius r_cut = 1.6 σ. The masses are set to unity so that the time scale of the simulation is given by τ = σ(m/ɛ)^1/2. Consecutive monomers (polymer beads) are connected by a nonlinear FENE bond with spring constant k = 30.0 ɛσ⁻² and maximum extension R₀ = 1.5 σ.^42,43 The cutoff radius of the LJ interaction of bonded monomers is set to r_cut = 2^1/6 σ, making the LJ interaction purely repulsive. The monomers interact with wall beads via an LJ interaction with ɛ = 0.6, σ = 1.3, and r_cut = 1.6 σ. The monomers are bonded to the bottom wall via FENE bonds with spring constant k = 30.0 ɛσ⁻² and maximum extension R₀ = 2.1 σ. The LJ interaction energy between the first bead near the wall and the wall is ɛ = 0.8 with length σ = 1.5 and cutoff radius r_cut = 2^1/6 σ. The top wall does not interact with the polymer chains and only serves to impose a shear flow on the solvent.

We explicitly model a dimer solvent made of chains of length N = 2.^44,45 This type of explicit solvent inhibits layering effects at the walls and permits a larger time step than a monomer solvent.³⁹ The dimers are connected by FENE bonds with k = 30.0 ɛσ⁻² and maximum extension R₀ = 1.3 σ. The LJ parameters are ɛ = 0.8, σ = 0.8, and r_cut = 2^1/6 σ. The non-bonded dimers interact via an LJ potential with parameters ɛ = 0.5, σ = 1.0, and r_cut = 2.5 σ. The parameters for dimer–polymer interactions are set to ɛ = 1.2, σ = 1.0, and r_cut = 2.5 σ. Finally, the dimers interact with the wall with ɛ = 0.6, σ = 1.0, and r_cut = 2.5 σ. The aforementioned parameters are identical to the parameters used in a study by de Beer and Müser⁴⁶ to induce good-solvent conditions. The combined polymer and solvent system is simulated at a fixed surface concentration of 0.8 beads σ⁻³. A downward pressure of 5 × 10⁻⁴ mσ⁻¹ τ⁻² acts on the top wall (which has no grafted polymer chains), slightly compressing the system. A representative system has the dimensions L_x = 100 σ, L_y = 87 σ, and L_z = 64 σ and contains 465 632 beads.

We use a DPD thermostat with friction coefficient γ = 5 mτ⁻¹ and cutoff r_c = 2^1/6 σ to thermalize the solvent and the polymers to a temperature of k_BT = 0.6 ɛ,^47,48 where k_B is the Boltzmann constant. The top wall is moved at constant sliding speed v along the x axis in order to induce a shear flow in the solvent. We use a time step of Δt = 5 × 10⁻³ τ in our simulations.

We start by analyzing the penetration of solvent flow into the polymer brush. Between the brush and outer wall [see Fig. 1(c)], we find an ideal Couette flow with a velocity profile of the solvent that is linear in channel position. The solvent is sheared at a rate $γ̇=v/Lh$⁠, where L_h is the height of the solvent’s Couette region. Note that for the definition of L_h, we extrapolate the linear Couette region into the brush toward the point where the velocity becomes zero. We use two different approaches to compute $γ̇$ and L_h: while we determine $γ̇$ from our simulations directly at high velocities, the shear rates obtained at low sliding velocities are dominated by thermal fluctuations. For sliding speeds v < v′, we compute $γ̇$ by assuming that L_h does not change from its value at v′ = 5 × 10⁻³ στ⁻¹.

Within the brush, hydrodynamic interactions are screened by the polymer chains—in other words, the flow velocity of the solvent vanishes inside the bulk of the brush. Given the overall simulation cell of height $Lh+Lh∗$ [see Fig. 1(c)], we obtain $Lh∗$ as the height of the quiescent part of the brush. Comparing the quiescent height $Lh∗$ with the actual height H of the brush yields the hydrodynamic-penetration length $λ=H−Lh∗$⁠. We compute the brush height H from the position at which the density has decayed to 1% of its maximum value [see Fig. 1(d) and Ref. 49].

The influence of the grafting density on hydrodynamic penetration can be qualitatively observed in the logarithmic depiction of the shear rate [Figs. 1(e) and 1(f)]. For systems with N = 60 beads/chain, the grafting density at which the chains begin to overlap is $Γ*=1/(πRg2)≈9⋅10−3σ−2$⁠. For Γ = 10⁻² σ⁻² ≈ Γ^∗, the shear rate gradually increases to its maximum value in the Couette region, while the transition is much sharper for the denser system with Γ = 10⁻¹ σ⁻² > Γ^∗.

For a more quantitative analysis, we plot in Fig. 2(a) the hydrodynamic-penetration length λ as a function of shear rate $γ̇$ for three grafting densities Γ. While for our densest system, λ is independent of shear rate over the range of rates, λ decreases at high shear rates for semi-dilute and dilute systems. However, for all systems and shear rates, the hydrodynamic-penetration length λ is proportional to the brush height H. Figure 2(b) shows the normalized hydrodynamic-penetration length, λ/H. The dependence on shear rate disappears for λ/H, showing that the change in penetration length with shear rate of Fig. 2(a) is simply because the chains tilt and, hence, the brush shrinks at high shear rates. In other words, the number (or mass) of monomers in the region where shear flow penetrates is constant. For the lowest grafting densities with Γ ≪ Γ^∗, the solvent flow extends almost throughout the whole brush, while at Γ ≫ Γ^∗, the flow is unable to penetrate the brush significantly.

Solvent flow in our system is driven by a non-zero shear stress at the wall boundaries. Our simulation geometry consists of two parallel plates with an entrained molecular system and is mechanically extremely simple: Mechanical equilibrium means that the shear stress σ_xz is constant along the simulation cell’s height (because ∂σ_xz/∂z = 0). It follows that the external shear stress measured from the constraint forces of our rigid walls $σxzwall$ is equal to the Newtonian viscous shear stress in the Couette-like region, $σxzvisc=ηsγ̇$⁠.

In the brush region, the shear stress will additionally contain a contribution from the polymer chains. In the simplest approximation (assuming independent chains), this will be the entropic stress in addition to the hydrodynamic, viscous stress of the solvent.³⁷ The entropic stress can be approximated from the individual chain conformations by Eq. (4). Note that Eq. (4) contains a volume average, which means it yields the mean entropic stress within the brush.

The key parameter mediating the reaction of the brush to the shear flow is the local shear rate $γ̇$ of the solvent rather than the sliding velocity v. The dynamic response of driven polymeric systems is, therefore, commonly nondimensionalized through the Weissenberg number $Wi∗=γ̇τr∗$⁠, where $τr∗$ is the equilibrium relaxation time of a single, dilute chain of length N. Examples include experiments and simulations in polymer brushes^39,40 or self-mated hydrogel friction.^11,12,41 The Weissenberg number compares an equilibrium measure (the relaxation time) with a non-equilibrium measure (the shear rate) and is proportional to $γ̇$ as $τr∗$ does not depend on shear rate.

We compute $τr∗$ from independent molecular dynamics simulations of free chains via the time-autocorrelation function of the trace of the gyration tensor $⟨Rg(0)Rg(t)⟩∝exp(−t/τr∗)$⁠.²¹ We have confirmed that the relaxation time computed from Zimm dynamics³⁸ [see Eq. (1)] is close (within a factor of 2) to the relaxation time obtained from single-chain simulations. Figure 4 shows the entropic stress as a function of Wi^∗. At low Wi ^∗, the shear stress increases linearly with Wi^∗. This regime crosses over to a sublinear regime at high Wi^∗.

We have pointed out already in the Introduction that the dynamics of the brush is controlled by an outmost blob of volume ξ³ rather than the dynamics of the full free chain; see, for example, Ref. 33. The relaxation time τ_r used in the definition of the Weissenberg number should, therefore, be the polymer relaxation time of this outermost blob, not the relaxation time $τr∗$ of the free chain. While experiments typically have access to $τr∗$ but not τ_r, we can compute τ_r from microscopic observations on our simulations. We assume that the individual blobs follow Zimm dynamics³⁸ and compute τ_r from the analytic expression Eq. (1). However, we do not know a priori the size ξ of the outermost blob. A crucial assumption of our model is, therefore, that ξ is proportional to the hydrodynamic-penetration length λ. In particular, we use $ξ=λ/6$⁠, since λ measures the penetration from the outermost part of the brush and corresponds to an end-to-end distance rather than a radius of gyration. Note that similar assumptions were made in early treatments of dense brushes.^33,37 Since λ changes with shear rate $γ̇$⁠, this implies that $τr(γ̇)$ also depends on shear rate through Eq. (1). This nonequilibrium Weissenberg number Wi $=γ̇τr(γ̇)$ hence depends nonlinearly on $γ̇$⁠.

We now analyze the entropic shear stress as a function of the nonequilibrium Weissenberg number Wi. Figure 5(a) shows a linear dependence of the shear stress on Weissenberg number above Wi ≈ 0.01. The linearity of the shear response exists independently of chain length, as evident from Fig. 5(b) that shows the entropic stress as a function of Wi for chains of different lengths at the same grafting density. Crucially, the sublinear scaling regime disappears when using the nonequilibrium Wi rather than the equilibrium Wi^∗ as the control parameter.

We note that the shear stress increases with grafting density, but the functional dependence of stress on Wi is independent of grafting density. Estimating the density for transition between mushroom and brushy regimes yields Γ^∗ ≈ 3.3 × 10⁻² σ⁻² for N = 20 and Γ^∗ ≈ 4 × 10⁻³ σ⁻² for N = 100. This means that the short chains are in the mushroom regime, while the longer chains are brushy, yet both show a qualitatively identical rate-dependence of the stress response.

A more detailed view is offered by computing the orientation and hardening factors that comprise the entropic stress. In Figs. 6(c) and 6(e), we observe that the orientation factors fall onto the same curve, increasing linearly with Wi above approximately Wi = 10⁻¹. For the densest brush analyzed, the hardening factor is independent of Wi over the range analyzed. In contrast, our less dense systems show an increase in the hardening factor starting at approximately Wi ≈ 1, consistent with the onset of chain elongation. Furthermore, their hardening factor is lower in the Wi-independent regime. When varying the chain length, we observe that the hardening factor increases with decreasing chain length in the Wi-independent regime. Furthermore, systems with longer chains show a stronger increase in the hardening factor at high Wi.

These trends are consistent with the expected conformations for flexible chains. At low grafting densities, the chains adopt random coil conformations with an equilibrium extension ratio h ∼ N^−1/2 that decreases with increasing N. At higher grafting densities, the overlapping chains elongate, increasing both h and the hardening factor in equilibrium.

Finally, we investigate how chain conformation behaves as a function of distance to the wall. We do so by computing the internal distance of chain subsegments with n_s bonds, normalized by the contour length L_c = n_sb of the said segment. The segments are chosen such that the first bead is tethered to the wall. Figure 7 shows that for the lowest grafting density, the normalized end-to-end distance is reasonably close to results for a single, dilute chain at low Wi. (The result for a dilute chain is shown by the dashed black lines in Fig. 7.) By comparison, as the grafting density increases, the normalized end-to-end distance in equilibrium deviates increasingly from the behavior of the single chain. At low Wi, the solvent flow does not affect chain conformations in all grafting densities analyzed. As Wi increases, the chains elongate significantly in both dilute and semi-dilute systems. While the flow’s effect can be observed close to the wall for the dilute system, the point at which differences can be first noticed shifts further away from the wall (higher n_s) for the semi-dilute system. The dense system’s morphology appears independent of Wi in the range analyzed.

The central observations from our simulations are as follows: The solvent flow penetrates into the brush up to a distance λ that increases with decreasing grafting density but is proportional to the rate-dependent brush-height, $λ∝H(γ̇)$⁠; the stress on the wall equals the entropic stress of the individual polymers, even at high grafting densities and high shear rates; the shear stress increases linearly with shear rate at low rates and sublinearly at high shear rates; and the shear stress increases linearly with the Weissenberg number Wi for all shear rates but only if the relaxation time that enters Wi is computed for the portion of the polymer chain penetrated by the shear flow.

We note that hydrodynamic penetration into brushes has widely been described in the literature (e.g., Ref. 30), and our results agree with those reports. Our observation $λ∝H(γ̇)$ suggests that solvent screening is mediated by a roughly constant amount of polymer mass both in and out of equilibrium.

The implication of solvent flow, namely the transfer of momentum to the boundaries of the system, has, to the best of our knowledge, not been correlated with the microscopic morphology of the chains. It is noteworthy that the chain conformational entropy fully accounts for the external wall stress in all cases considered. This implies that even at low grafting densities, viscous solvent forces and intermolecular interactions contribute negligibly relative to the tension transmitted through chain backbones. This observation is consistent with the recent results on extensional stress in dense polymer melts, which can be related to conformation using the same freely jointed chain model.²² The model even holds in highly nonequilibrium situations as encountered during extensional flow or at high shear rates.

At low shear rates, the stress increases linearly with shear rate as expected from linear response theory. However, we note that shear-rate independent regimes have been reported in simulations of brush⁵⁰ and hydrogel⁴¹ friction. We do not observe such a rate-independent friction (also called Coulomb or adhesive friction) regime, which typically requires some form of elastic instability for its emergence.⁵¹ While we cannot rule out that it exists, it is likely that Coulomb friction only emerges if the grafted chains can attach and detach from a counterbody,^52–56 which they cannot do in our case.

To quantify whether a shear-rate $γ̇$ is “high” or “low,” we need to compare it to a characteristic relaxation time of the polymers. The relevant time scale is the polymer relaxation time. Typically, the Weissenberg number is defined as $Wi∗=γ̇τr∗$ using the equilibrium relaxation time of a free chain. Our simulations show a transition to a nonlinear increase in the shear stress with shear rate at Wi^∗ ≈ 0.01. We point out that a similar sublinear behavior at high Wi^∗ has been found in the friction coefficients of monodisperse brushes^39,40 and self-mated hydrogels,^10,11 with a transition also occurring near Wi^∗ ≈ 0.01.^12,41 For self-mated hydrogels, this regime has been attributed to both non-equilibrium polymer effects¹⁰ and hydrodynamic lubrication.⁹ Our results indicate that the shear response of interfacial polymer chains might lie behind the observed frictional behavior of self-mated hydrogels.

The decomposition of the entropic stress into the orientation and hardening factors [Figs. 6(c)–6(f)] allows us to connect the shear response to morphological changes in the brush. At low shear rates, the brushes remain close to equilibrium and there is no clear signature in orientation and hardening factors. The chains reorient roughly at shear rates high enough that the stress starts to depend sublinearly on Wi^∗, which occurs at Wi^∗ ≪ 1. Our simulations indicate that the process occurring at Wi^∗ = 1 appears to be chain stretching, which can be directly seen in the hardening factors. The reason that we see chain stretching at this Weissenberg number is that the standard relaxation time, computed from the time-autocorrelation function of radius of gyration, corresponds to the slowest stretch relaxation time of the chain. Orientational relaxation typically occurs at longer time scales.⁵⁷ 

In a polymer brush, the shear response should depend on the time scale of the relaxation τ_r of the brush’s outermost boundary layer that is penetrated by the shear flow³³ and not the equilibrium relaxation time of the full chain $τr∗$⁠. This length scale is the hydrodynamic-penetration length λ, which is independent of the shear rate for our densest system only [Fig. 2(a)]. As a result, it is valid to approximate the system’s behavior with a constant relaxation time for a dense brush. However, for the semi-dilute and dilute systems, the hydrodynamic-penetration length and thereby the polymer relaxation time decrease significantly at high shear rates, which can be attributed to a decrease in the brush’s height [Fig. 2(b)].

The relaxation time τ_r of the outmost part of the brush can be estimated from Zimm dynamics, Eq. (1), with a blob size ξ ∼ λ. Using this (shear rate-dependent) relaxation time to compute the Weissenberg number, $Wi=γ̇τr(γ̇)$⁠, we observe that the sublinear regime at high Wi^∗ becomes linear. This behavior emerges at all grafting densities [Fig. 5(a)] and chain lengths [Fig. 5(b)]. These results are in agreement with the assumption made in previous theoretical models that polymer brushes in shear flow can be described as strings of Pincus blobs.^{31,32,35,37,58}

Detailed morphological insights can be obtained from the conformations within the chain. The normalized internal distance (Fig. 7) becomes susceptible to shear flow closer to the wall in the dilute system than in the semi-dilute system. This can be attributed to the larger solvent flow within the chains in the dilute system, resulting in a larger hydrodynamic-penetration length. Our results show that chain conformations become increasingly linear due to the influence of the shear flow. The lack of Wi-dependency in the dense system can be attributed to the lack of shear flow penetration into the brush.

The normalized internal distance also sheds some light onto the relation between the orientation and hardening factors in dilute and semi-dilute systems. Even though in shear flow, the strongest shear force acting on the chains is located at the outermost blob, the reorientation of the chain involves the relaxation of the full chain. However, in Fig. 7, we observe that chain elongation is not distributed equally among the chains’ backbone but becomes stronger further away from the grafting site, i.e., it becomes stronger toward the location of the outermost blob.

In summary, our results show that the brush’s outermost “blob,” the portion of the polymer chain that is exposed to shear of the solvent, determines the system’s shear response. The flow penetrates less into dense than into dilute brushes. For dilute brushes, the flow penetrates less at high shear rates because the chains tilt. This leads to a shrinking, and effectively also a densification, of the brush. As a result, the relaxation time relevant to the system’s shear response decreases at high shear rates. If this mechanism is not considered in the definition of the Weissenberg number—the quantity typically used to characterize the effect of solvent shear on polymeric systems—then the shear stress will appear to depend sublinearly on Wi. This phenomenology has been observed in polymer brush and self-mated hydrogel simulations and experiments at high Wi.^{12,14,17,39–41}

We used LAMMPS⁵⁹ for all simulations and OVITO Pro⁶⁰ for visualization and post-processing. Computations were carried out on NEMO (University of Freiburg, DFG Grant No. INST 39/963-1 FUGG) and HOREKA (project HydroFriction). T. O. acknowledges startup funding from Carnegie Mellon University. We acknowledge the Deutsche Forschungsgemeinschaft for providing funding within Grant Nos. PA 2023/2 and EXC-2193/1-390951807.

The authors have no conflicts to disclose.

Jan Mees: Conceptualization (equal); Data curation (equal); Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal). Thomas C. O’Connor: Conceptualization (equal); Methodology (equal); Validation (equal); Writing – review & editing (equal). Lars Pastewka: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

We compute Flory’s characteristic ratio from polymer melt simulations with chains of length N. The systems were first equilibrated at a temperature of T = 0.6 ɛ and a pressure of P = 5 × 10⁻⁴ mσ⁻¹τ⁻². This procedure produced systems with a monomer density of approximately ρ = 0.87 σ⁻³. A second equilibration at constant volume followed at an increased temperature of T = 2.0 ɛ to increase diffusivity. Finally, at T = 0.6 ɛ, the end-to-end distance of each chain was sampled every 10³ time steps for 10⁵ time steps. Flory’s characteristic ratio was computed for each system using the relation C_n = ⟨R²⟩/nb², where n = N − 1 is the number of bonds in a chain.³⁸ 

We see in Fig. 8 that Flory’s ratio increases rapidly with N for N < 30 beads and saturates for higher N at ∼1.7. It follows that C_∞ = 1.7.