Effects of disorder and chain stiffening on the elasticity of flexible polymer networks
We examine how the distribution of contour lengths and the high-stretch stiffening of individual chain segments affect the macroscopic shear modulus of flexible polymer gels, using a 2D numerical model, in which polymer segments form a triangular network and disorder is introduced by varying their contour lengths. We show that in the relevant parameter range: (i) the non-affine contribution to the shear modulus is negligible, i.e. the Born approximation is satisfactory; (ii) the shear modulus is dominated by the contribution originating from equilibrium chain tensions. Moreover, mechanical equilibration at the nodes induces specific correlations between the end-to-end distances and contour lengths of chain segments, which must be properly accounted for to construct reasonable estimates of chain pressure and shear modulus.
In the last decades, a lot of attention has been dedicated to understanding how the macroscopic elasticity of polymer hydrogels is controlled by the properties of single polymer strands and by network disorder. Most of these studies Kroy (2006); Broedersz and MacKintosh (2014); Pritchard et al. (2014) have focussed on the case of filamentous networks, in which the constitutive polymers have persistence lengths larger that the network mesh size – which is relevant to many biopolymers. They have shown, in particular, that large chain stiffnesses give rise to spectacular mechanical responses, such as pre-stress-induced stiffening or softening, negative normal stress, etc. Few works Carrillo et al. (2013); Meng and Terentjev (2016), however, have addressed recently the case – which is relevant to many hydrogels, such as gelatin or synthetic polymer gels – where the persistence length of polymer strands lies in an intermediate range between the monomer size and the mesh size .
The elasticity of such “flexible polymer networks” is usually treated in the framework of the classical theory of rubber elasticity, which suggests that their shear modulus is (in 3D). This expression is often used in practice to estimate network mesh-sizes from measurements of the shear modulus, but only arises from approximate scaling arguments de Gennes (1979), or from a mean-field theory Treloar (1949) where the end-to-end distances of chain segments are supposed to follow the same Gaussian statistics as if they were independent. This assumption overlooks the fact that chains fluctuate around a state of mechanical equilibrium; as we shall see, it amounts to a very stringent ansatz about the correlations between end-to-end and contour lengths.
In this paper, we study in detail how the distribution of contour lengths and the high-stretch stiffening of individual chain segments affect the macroscopic shear modulus. We address these issues using a 2D numerical model of gel, in which disorder is introduced by varying the contour lengths of polymer strands that form a triangular network. The model takes into account the existence of excluded volume effects à la Flory-Rehner Onuki (2002), and is tested for different expressions of the single chain response to large extensions. The shear modulus of such a system can be decomposed into a sum of three terms:
which are respectively the opposite of the chain pressure , the elastic “constant” , and the non-affine term . Even though the respective importance of these three terms depends on disorder strength and swelling level as well as on the importance of stretch-stiffening, we show that, in the range of parameter values relevant for usual flexible gels, these three contributions are ordered according to the following hierarchy: . Namely, the non-affine contribution to the shear modulus is small enough for the Born approximation to provide, for all practical purposes, a very satisfactory estimate of the shear modulus. The elastic constant contributes at most a fraction of less than typically 15–20% and is all the weaker that the persistence length is small. The chain pressure term is always dominant, but its value depends on both the distribution of chain segment stiffnesses and on the accomodation of elastic disorder by mechanical equilibration. This hierarchy of the contributions is specific to flexible polymer gels, in constrast with semi-flexible or rigid polymer networks.
We choose to specialize to 2D networks with a fixed, triangular topology, as illustrated on Fig. 1. The nodes represent permanent crosslinks and the links flexible polymer strands that rotate freely at the nodes. The number of monomers on the strand connecting nodes and is denoted .
ii.1 System free-energy
We define the system free-energy as , where the two terms account respectively for:
(i) the elastic free-energy of individual strands in an ideal solvent,
where the sum runs over all strands , with the distance between nodes and , and
the subscript in accounts for its dependence on .
(ii) a Flory-Rehner-like contribution Onuki (2002) that accounts for excluded volume effects. To model it, we consider that the number of monomers lying in a given network triangle is . The average monomer concentration in a network triangle is thus
with the triangle area. In the spirit of the Flory-Rehner mean-field approximation Onuki (2002), we then set:
where the sum runs over all triangles, and is the 2D excluded volume parameter.
Now we need to specify the expression of the elastic free-energy of individual chains that enters equation (2). We consider that the persistence length is small compared with the average mesh size. So long as the end-to-end distance of any strand is much smaller than its contour length , elasticity is purely entropic, and we can use for the standard Gaussian expression Rubinstein and Colby (2003):
It is well-known, however, that this expression is insufficient at high stretch, i.e. when approaches 1.
A standard model for the large stretch response is the freely jointed chain Rubinstein and Colby (2003), for which the relation between elongation and applied force is provided, in 3D, by a Langevin function. Although the associated free-energy does not possess an explicit analytic expression, a satisfactory approximation is provided by the Cohen expression Cohen (1991), which interpolates between the small stretch () Gaussian and the nearly taut () limits. In Appendix A we derive the analogous expression for a 2D freely-jointed chain containing Kuhn segments; it reads:
As needed, reduces to the Gaussian expression in the limit .
It is known that the logarithmic growth of the freely-jointed chain free-energy underestimates the hardening of the single chain response at high stretch levels. This effect certainly becomes all the more important that the persistence length increases and, in the limit of high stiffnesses, the worm-like chain (WLC) model provides a much better description of the stretch-force relation. To tackle the problem of intermediate persistence lengths, Blundell and Terentjev Blundell and Terentjev (2009) have proposed a model of the single chain response that interpolates between the Gaussian and the WLC relations. This provides a third model for the elastic free-energy111In order to ensure that matches in the , limit, our definition of the persistence length differs by a factor of from Blundell and Terentjev’s:
In the following, we will probe the linear and non-linear elastic response of networks using these three expressions for the single chain free-energy.
ii.2 Units and parameter ranges
We will use parameter values in ranges that are reasonable for the thoroughly investigated hydrogels of gelatin, for which the monomer size , persistence length , and Bohidar and Jena (1994). Typical mesh sizes lie in the range, i.e. . A rough evaluation of the average number of monomers per strand is , on the order of 100.
When presenting numerical data, the monomer size and will be taken as units of length and energy.
ii.3 Network disorder
To construct a triangular network, the node points are initially placed on a Bravais lattice with vectors and , in a biperiodic cell of extension , with and integers. There are nodes for a cell area , i.e. an average areal density of nodes . The number of strands is and the number of triangles .
The network topology being fixed, disorder is introduced via the values of the monomer numbers , which we take to be random and uncorrelated variables. For the sake of simplicity, and to facilitate the qualitative analysis of disorder effects, we assume their distribution to be bimodal:
For any realization of the set , mechanical equilibrium is then found by minimizing the total free-energy of the system, which results in a distorted network. This is illustrated on Fig. 1, which displays two mechanically equilibrated configurations with the bimodal disorder defined by and .
Iii Elastic network response
iii.1 Affine vs non-affine contributions
To compute elastic constants, we rely on the general formalism developed in Lemaître and Maloney (2006) for the elastic response of disordered solids. Its main lines are briefly summarized as follows.
Let us consider some initial (reference) state about which we compute the elastic response. The externally imposed macroscopic deformation about this initial state is specified via the strain tensor . In the initial state, and the nodes assume equilibrium positions denoted . Under deformation, an arbitrary configuration of the system is defined by the macroscopic strain and the node positions . For any , we define its zero strain antecedent as . Clearly in the initial configuration, where .
The Born approximation of affine deformation amounts to assuming that node positions vary with as with fixed antecedents that coincide with the initial node positions. To separate the affine and non-affine contributions to the elastic response, it is convenient to write formally the total free-energy of the deformed system in terms of the antecendents of the (a priori arbitrary) node positions. This is realized by writing
which defines .
We are interested in the static elastic response, which entails that, under loading, the system remains at mechanical equilibrium, i.e. that the force on each node vanishes at all times:
We denote the node positions at mechanical equilibrium under strain . In a disordered system, as varies, the follow trajectories that, in general, are not affine. It means that their antecedents are not fixed, but vary with . Their trajectories are specified by derivating, with respect to each component of the strain tensor, the condition
which is an immediate consequence of Eq. (11). In the limit , it comes:
(note that we use the convention of implicit summation on repeated indices). Here,
is the Hessian matrix in the reference (initial) state, and the vector field is defined as:
Note that taking partial derivatives of with respect to strain components amounts to varying at constant , i.e. to performing affine deformations about state . Since moreover, in the limit , tends to the initial (reference) configuration, it turns out that can be interpreted as the force induced by an infinitesimal affine deformation Lemaître and Maloney (2006). In view of equation (13), the non-affine displacement field characterized by can be interpreted as the linear elastic response of the system to this field of virtual forces.
This framework can be used to write explicit expressions for stresses and elastic stiffnesses, which are first and second derivatives of the free-energy with respect to strain. In particular, it has been shown that elastic stiffnesses Barron and Klein (1965), defined as
with the system area, can be decomposed as Lemaître and Maloney (2006):
is the Born approximation for the stiffness tensor, which assumes that the nodes follow affine trajectories, and the non-affine contribution
results from the non-affinity of the displacement field.
iii.2 Microscopic expression of the stress tensor
Since , the macroscopic Cauchy stress tensor is a sum of a chain and a Flory contribution
The “chain stress” is given by the classical expression for systems with pair interactions, namely:
with . The “Flory stress”, derived in Appendix LABEL:app:details, which reads
is diagonal under any state of deformation, as expected from the microscopically isotropic character of the Flory interaction; it hence only contributes to the osmotic pressure
but not to deviatoric stresses.
iii.3 Response to simple shear
Since hydrogels are incompressible on time scales where poroelasticity is irrelevant, their elasticity is characterized, for all practical purposes, by their response to shear. We thus focus here on the case of simple shear deformation, for which
We first calculate explicitly the Born contribution , which is obtained under the assumption of affine node displacement. Denoting the initial node positions, any strand end-to-end vector is transformed into , i.e.:
Under such an affine displacement, which preserves areas, the Flory free-energy , defined by Eq. (4), remains invariant. The Born modulus [Eq. (18)] hence reduces to . It is computed by writing the shear stress at arbitrary :
and derivating once more. It will be useful to decompose the result as follows:
is the component of the stress carried by the chain network in the initial, undeformed, state.
Let us note that the chain free-energy , which determines , only depends on the distances between connected nodes. In such a case, the free-energy under deformation can be written as a function of the Green-Saint-Venant tensor , since , and the general elastic theory for discrete systems shows that the elastic stiffness tensor can be written as Lemaître and Maloney (2006) , where is called the tensor of elastic constants, and is the stress in the undeformed system. Expression (27) corresponds exactly to this decomposition since is precisely the elastic constant . It should be emphasized that, because the Flory free-energy is invariant under simple shear and, consequently, determined by the variations of , the stress term is only but not the total stress.
which leads to Eq. (1) when the stress tensor is isotropic, a condition which, as we will see shortly, is satisfied by our triangular networks. The non-affine contribution [Eq. (19)] involves the field of virtual forces and the Hessian matrix . The explicit expression of is provided in Appendix LABEL:app:Hessian. Concerning , we note that, in view of Eq. (15) and since the Flory free-energy is invariant under affine simple shear (), there is no Flory contribution to it.
It is worth noting that, since the Flory contributions to both and vanish, excluded volume effects impact the linear elastic response only indirectly via the role they play in defining the equilibrium structure.
Iv A simple case: the homogeneous network
We focus in this section on the case when all strands have an equal number of monomers . Then, the nodes lie on a regular lattice at any level of deformation and non-affine effects are absent.
iv.0.1 General expressions
As expected for a 2D triangular lattice, is isotropic: . Note that changing the swelling state of our gel network amounts to varying at fixed . Accordingly is the osmotic pressure.
From equation (27) we obtain for the shear modulus (which reduces to its Born estimate):
iv.0.2 The Gaussian chain model
where , the mesh size at swelling equilibrium (), reads:
It is seen to be proportional to , the standard 2D scaling behavior expected in the Flory framework. When , then and the gel is underswollen.
Turning to the value of the shear modulus, we first note that, since is quadratic, the elastic constant [Eq. (34)] vanishes so that:
We recover here the classical expression, deduced from rubber elasticity Rubinstein and Colby (2003), which states that the shear modulus of a polymer gel is proportional to the chain number density times the elastic free-energy per chain.
iv.0.3 Dependence on swelling level for the three chain models
Note that varying at fixed amounts to changing the swelling level of a gel of fixed network structure. So, the above expression means that, for the Gaussian chain model, is independent of the swelling level.
To compare the different chain models, we study a homogeneous network of fixed , and plot on Fig. (2) the four quantities , , , and which characterize the gel mechanical state, versus the mesh size , which characterizes the swelling level. The sharp drop of the osmotic pressure with swelling at small is essentially due to the decay of the Flory pressure [Eq. (32)], which is identical for all models.
Let us recall that the FJ and BT chain free-energies account for the stretch hardening of the chain segments while matching the Gaussian expression at small . Expectedly, the values of , , , and obtained with both models smoothly grow away from the Gaussian ones with increasing . The BT model data exhibit a much steeper dependence on the swelling level, in agreement with the fact that it interpolates at high stretch with the WLC behavior which is stiffer than the FJ one.
The conditions of our study should be contrasted with the numerous existing works on filamentous networksKroy (2006); Broedersz and MacKintosh (2014); Pritchard et al. (2014), which usually deal with the very stiff regime, . Here, on the contrary, we consider much more flexible gels with the persistence length much smaller than the contour length , at moderate stretch ratio . It is thus striking that both the FJ and BT models lead to a very substantial growth of the shear modulus with swelling level, so that it departs from the Gaussian prediction by amplification factors of respectively 1.5 (FJ) and 4.5 (BT) at our highest stretch ratio .
iv.0.4 Dependence on cross-link density at fixed monomer concentration
The progress of the crosslinking reaction under undrained conditions corresponds in our model to decreasing values of at fixed monomer density . Rubber elasticity theory Treloar (1949) then predicts that , with the space dimension, an expression which is commonly used to estimate mesh sizes from the tracking of the shear modulus in experiments, yet is derived under the assumptions that chains are Gaussian and deformations affine.
We expect our Gaussian chain homogeneous network model to satisfy this scaling relation since it uphelds both assumptions. Indeed, rewriting expression (38) in terms of the monomer density yields:
Since this scaling derives from purely entropic arguments, it should break down in our FJ and BT models, which take stretch hardening into account. To probe the magnitude of the expected deviations, we plot on Fig. 3 , , , and as a function of for networks of a fixed monomer density corresponding to , , . As the crosslink density increases (decreasing ), the stretching ratio grows, and so does the chain tension: is increasingly negative. Since only depends on , hence is - and model-independent, we recover the intuitive result that osmotic pressure decreases together with , i.e. that the swelling level grows with crosslinking.
Chain pressure is found to be only weakly model-dependent, i.e. is hardly affected by stretch hardening. As for the elastic constant , its FJ value weakly departs from the Gaussian limit, but grows significantly at low in the BT model, in agreement with the fact that it specifically captures the deviations from harmonicity of the chain potential [see Eq. (34)]. The resulting departure of the shear modulus from the rubber elasticity scaling is quantified by plotting vs in the insert of Fig. 3. Over the whole considered range (from to ) we find that increases from 1.2 to 1.6, which correspond to quite substantial relative deviations.
V Disordered networks
We now introduce disorder by attributing, as explained in II.3, random values to the bond monomer numbers . To isolate the effect of disorder strength, we vary the distribution width for fixed average monomer and crosslink densities, i.e. fixed average and . The results, averaged over 100 configurations, are shown, for the bimodal distribution [Eq. (9)], on Fig. 4.
v.0.1 Osmotic and chain pressures
Panel (a) displays the osmotic pressure and its two terms, the chain and Flory pressures, and . Over the considered range, is not strictly constant, but its variations are much too small to be visible on the graph. The changes of osmotic pressure hence only result from the changes of .
Visibly, decreases, i.e. the average chain tension increases with disorder. A hint to understand this trend is provided by constructing a regular network approximation (RNA) which neglects the fact that the mechanically equilibrated network is distorted with respect to the reference triangular lattice. In this approximation, all end-to-end vectors assume the same length and orientations as in the homogeneous problem, and the chain pressure reduces to an expression similar to Eq. (31): where the average is taken over the distribution. One easily checks that, for the Gaussian and FJ models, the chain tension is a concave function of at all stretch ratios (). The same holds for BT as long as , a condition satisfied by gels of flexible polymers at not too large underswelling levels. This property entails that (i) , the chain pressure of the truly homogeneous network (for which , i.e. ); (ii) at fixed , decreases (in algebraic value) with .
We test this approximation in the simplest case of the Gaussian model by plotting and the true of the inhomogeneous problem (insert of Fig. 4a): it appears that noticeably overestimates the effect of disorder. To understand the origin of this discrepancy, let us write the microscopic expression of the chain pressure:
where stands for the pair average. In the case of Gaussian chains, , while . The comparison between and points to the importance of correlations between chain end-to-end distances and monomer numbers . Indeed, neglecting these correlations would lead to: . Since, in a triangular network, always holds222 To check this inequality it suffices to notice that the energy of a network of Gaussian chains of unique is minimized, at fixed monomer density, i.e. fixed , by the homogeneous state with ., this assumption would result in , in contradiction with the data.
To evidence the correlations between chain end-to-end distances and monomer numbers, we plot on Fig. 5, for two bimodal disorder strengths, the distribution of ’s and its decomposition into the contributions of short and long chains. The two sub-distributions are clearly separated. Hence, in both cases, corresponding to large and modest elastic contrasts, chain lengths and monomer numbers are strongly correlated.
To understand why long (softer) chains are more extended than short (stiffer) ones, note that in the regular network approximation, the stretch ratios of long and short chains take the values : short chains are more taut (than long ones) and hence pull more strongly on their surroundings. Under the effect of mechanical equilibration, short chains thus tend to contract while the long chains expand, which reduces the contrast of stretch ratios (and hence of chain tensions). Nevertheless, this effect does not fully resorb the stretch ratio contrast: the data of Fig. 5 show that, in mechanical equilibrium, the stretch ratio of the short chains peaks around 0.35 for and 0.5 for , while those of the long chains peak around 0.26 for and 0.23 for . In mechanical equilibrium, short chains always remain on average more taut than long ones, and increasingly so for large disorder strength. Mechanical equilibration hence only mitigates the effect of chain stiffness disorder as captured by the regular network approximation.
This observation sharply constrasts with the assumptions of rubber elasticity, which overlooks mechanical equilibration, and postulates that chains are independent, with a Gaussian statistics for end-to-end vectors. It follows that is independent of , which leads to predicting a constant . Since [see Eq. (28)] vanishes for Gaussian chains, rubber elasticity predicts the shear modulus to be disorder-independent.
v.0.2 Contributions to shear modulus
Let us now turn to panel (b) of Fig. 4, which displays the shear modulus along with the three terms of its decomposition . In all cases, the chain pressure term provides the largest contribution. The elastic constant captures departures from Gaussianity. It vanishes for Gaussian chains. It remains weak for the FJ model, is comparable to for the BT one, and in both cases grows with disorder strength. This latter effect results primarily from the fact that in both models is a concave function of , a property which is expected to hold rather generally for stretch-hardening chain free-energies. Of course, like chain pressure, the exact value of in the FJ and BT models depends on the correlations evidenced in Fig. 5 between chain end-to-end distances and monomer numbers. Yet, as we saw above, mechanical equilibration does not destroy the constrast of stretch ratios between short (more taut) and long chains, and hence the value of grows with disorder strength.
The non-affine contribution also increases with disorder, as intuitively expected. Yet, strikingly, it only contributes a small fraction of the shear modulus: remains smaller than 1% for both the Gaussian and FJ models; it reaches at most for the highest disorder strength with the BT model. It thus turns out that the effect of disorder on the shear modulus is, for all practical purposes, captured by the Born approximation: . As both terms grow similarly with disorder, so does the shear modulus .
v.0.3 Response to swelling
In order to illustrate how disorder impacts the elastic response of the network to swelling, we compare on Fig. 6 the swelling response of a network with its homogeneous () counterpart, both using the BT model of chain free-energy. For pedagogical purposes, our swelling range extends up to a very high – of course unrealistic – level (). All the measured quantities vary much more steeply with swelling level, in the disordered network. This can be ascribed to the hardening sensitivity of the short chains. For example, at , the end-to-end distances of shorts chains are strongly peaked around , which corresponds to a very high stretch ratio . By contrast, the long chains exhibit rather modest stretch ratios , their end-to-end distances being broadly distributed about . It is thus primarily the short chains that are responsible for the disordered-induced enhancement of stiffening with swelling.
Here we have analyzed in detail the case of flexible networks, in which the persistence length remains substantially smaller than the contour length of inter-crosslinks polymer chain segments. Their modulus can be decomposed into three terms,
which present different sensitivities to chain flexibility and network disorder.
We have shown that:
(i) The non-affine contribution vanishes in the limit of either very flexible (Gaussian) chains or of homogeneous (non-disordered) networks. As a consequence, it contributes only a few percents of if we limit ourselves to a realistic range of disorder and swelling levels. Hence, the shear modulus is very well captured by the Born approximation:
(ii) The elastic constant vanishes in the limit of flexible (Gaussian) chains, but is non-zero for stretch-hardening chains even in the absence of disorder. It thus presents values that are systematically larger that . However, it contributes only a small fraction of at reasonable disorder and swelling levels.
(iii) In all cases, is dominated by the chain pressure contribution , all the more so that chains are more flexible.
This hierarchy of importance between the three contributions () is specific to flexible networks, in contrast with networks made of rigid polymers exhibiting finite average end-to-end distances on the scale of the mesh size .
Of course, the specific value of (and of ) depends on both the distribution of chain stiffnesses and on the accomodation of this elastic disorder by mechanical equilibration. A trivial effect of disorder – captured in the regular network approximation, which assumes that crosslinks lie on a regular lattice – is that shorter (stiffer) chains tend to be more taut than long (softer) ones, which entails that disorder strongly enhances the stretch-hardening effects. Mechanical equilibration mitigates this effect by introducing correlations between end-to-end distances and chain contour lengths. These joint effects bring in a noticeable contribution to even in the case of fully-flexible (Gaussian) chains, in contradiction with the prediction of rubber elasticity.
Acknowledgements.We are grateful to Tristan Baumberger for stimulating our interest in the subject and for illuminating discussions.
Appendix A Freely-jointed chain in 2D
We compute the 2D polymer entropic elasticity following the usual route. A chain is supposed to comprise monomers of size . To take into account the finite persistence length, it is decomposed into Kuhn segments of fixed size and orientation . The Kuhn segments are attached at their end points and can rotate and overlap freely. One end of the polymer is held fixed, and a force is applied at the other, which amounts to introducing a potential , with the end-to-end vector. Under equilibrium at temperature , the partition function is with
with and a modified Bessel function of the first kind. The average end-to-end vector is:
so that the norm of the end-to-end distance is , with .
In order to invert approximately this relation, we note that for small , at lowest order, , whence ; and for large arguments, , whence . The inverse relations, respectively (when ) and (when ). In the spirit of the Cohen approximation for the 3D case Cohen (1991), we interpolate between these two limits using the rational function . So we write . From this, we derive the expression for the elastic the free-energy of the freely jointed chain that appears in Eq. (6):
Appendix B Calculational details
The total free-energy is of the form: , with the chain and Flory contributions:
In these expressions, and respectively index pairs and triangles; is the norm of the difference vector between the positions of nodes and ; is the area of triangle . It is convenient, for any vector to define ; also, we assume – without loss of generality – that all triangles are oriented counterclockwise, so that the area of a triangle reads .
b.1 Microscopic forces
The force exerted on node can be decomposed into the contributions of both the chain and Flory free-energies:
The chain contribution is of the form:
where the summation counts any pair once, and where denotes the force exerted by on :
The last equation is Newton’s second law. The Flory contribution reads:
is the force caused on node , by the three body interaction between nodes , and . Of course, it is non-zero only if . Thanks to the translation and rotation invariance of , the forces it induces on the summits verify: