Length-Controlled Elasticity in 3D Fiber Networks
LLength-Controlled Elasticity in 3D Fiber Networks
C. P. Broedersz,
1, 2
M. Sheinman,
1, 2 and F. C. MacKintosh
1, 2 Department of Physics and Astronomy, Vrije Universiteit, Amsterdam, The Netherlands Kavli Institute for Theoretical Physics, University of California, Santa Barbara, California 93106, USA (Dated: October 5, 2018)We present a model for disordered 3D fiber networks to study their linear and nonlinear elasticityover a wide range of network densities and fiber lengths. In contrast to previous 2D models, these3D networks with binary cross-links are under-constrained with respect to fiber stretching elasticity,suggesting that bending may dominate their response. We find that such networks exhibit a fiberlength-controlled bending regime and a crossover to a stretch-dominated regime for lengths beyonda characteristic scale that depends on the fiber’s elastic properties. Finally, by extending the modelto the nonlinear regime, we show that these networks become intrinsically nonlinear with a vanishinglinear response regime in the limit of floppy or long filaments.
PACS numbers: 83.10.Tv, 62.20.de, 87.16.Ka, 83.60.Df
Materials ranging from paper and textiles to the struc-tural components of living cells and tissues[1] all exhibitnetworks of fibers or stiff polymers. Such networks haveextraordinary mechanical properties[2–4]. Their elastic-ity depends in part on their connectivity[5, 6], in analogywith jammed matter[7, 8] and random network glasses[9].The mechanics of the constituent fibers, and specifi-cally their bending rigidity can also strongly impact net-work elasticity. However, the relative importance of fiberstretching versus bending is not understood, especially in3D. Prior work has mostly focused on 2D networks[10–15]since simulations in 3D have been proven to be challengingand have usually been limited to small system size[16–18]. Significant qualitative differences are expected be-tween 2D and 3D networks: for the typical case of binaryfiber interactions, the high-molecular weight limit in 2Dactually corresponds to the Maxwell central-force (CF)isostatic threshold, where stretching interactions begin tocompletely constrain network deformations. In contrast,3D networks with binary interactions remain well belowCF isostaticity. Thus, owing to their marginal stability,real 3D fiber networks are expected to be fundamentallymore bending-dominated and more prone to collectivenonaffine deformations[16, 18].Here we develop a numerical model for the elasticity ofrandom 3D fiber networks with binary cross-links. Thismodel provides access to a wide range of network densi-ties below the CF isostatic threshold, as well as the previ-ously inaccessible high-molecular weight limit. These net-works exhibit various qualitatively distinct linear elasticregimes: a critical regime governed by the rigidity perco-lation point, a length-controlled bending regime and anaffine stretching regime, as illustrated in Fig. 1a. We pro-vide a scaling analysis for insight into the origins of theseregimes. Network mechanics is determined by the ratiobetween a nonaffinity length scale and molecular weight.The high-molecular weight limit exhibits behavior rem-iniscent of zero-temperature critical behavior, includingdivergent strain fluctuations; these ultimately govern a nonaffinelength-controlledbendingcriticalbendingcriticalstretching affinestretching -2 linear nonaffinelength-controlledbending linear affinestretchingnonlinearbend-stretchstiffening a) b) Small strain elasticity Strain dependence
FIG. 1: (Color online) Phase diagrams for the linear (a) andnonlinear elasticity (b) of 3D fiber networks on the PhantomFCC lattice, where L is the average filament length, z is net-work connectivity, γ is strain and κ is the fiber bending rigid-ity. All lengths are measured in units of the lattice spacing (cid:96) and κ in units of µ(cid:96) . crossover from the length-controlled bending regime to anaffine, stretching dominated network response. This is re-markable, since the network remains well below Maxwell’sCF isostatic connectivity threshold in this limit. Thus,paradoxically, although such networks can only be rigidat non-zero fiber bending stiffness, we find that no mat-ter how weak this bending rigidity is, the network elas-ticity becomes insensitive to fiber bending in the limitof high—yet finite—molecular weight. Moreover, in thelimit of floppy filaments with weak bending rigidity orhigh molecular weight, these networks become intrinsi-cally nonlinear with a vanishing linear response regime(Fig. 1b).Much has been learned about stiff polymer gels fromminimal models, such as 2D Mikado networks of randomlyplaced straight filaments with binary cross-links[10, 11].The elasticity of such Mikado networks is governedby nonaffine fiber bending deformations at low densi-ties, while higher density networks exhibit predominantly affine stretching elasticity of single fiber segments. Thisnonaffine-affine (NA-A) transition can be understood asbeing the result of increasing fiber-length. However, for a r X i v : . [ c ond - m a t . s o f t ] A ug −2 −1 −6 −5 −4 −3 −2 −1 a)b) FIG. 2: (Color online) a) The Shear modulus as a functionof L in units of (cid:96) for various κ in units of µ(cid:96) . G is the(affine) shear modulus of the undiluted Ph-FCC lattice. Theinset illustrates the phantom principle: At each lattice vertex3 independent binary cross-links are formed between randomlychosen fiber pairs labeled by color. b) Non-affinity parameterΓ as a function of L . Dashed black lines indicate a slope of 2. such 2D networks, this high molecular weight limit actu-ally coincides with Maxwell’s CF isostatic connectivity, z CF = 2 d in d dimensions[5], which can also give riseto a NA-A transition[6]; it is thus unclear whether theobserved transition in 2D is controlled by CF stretchingconstraints, or by filament length, as previously suggestedby scaling arguments and floppy mode theory[10–12].However, 3D networks with binary cross-links—like mostbiopolymer systems—are qualitatively different; in thiscase the high-molecular weight limit corresponds to net-work connectivities well below z CF . In the absence of fiberbending resistance, such networks exhibit zero-energy de-formation modes and hence, they do not resist shearstresses. Thus, there are reasons to question the existenceof an affine limit in realistic 3D networks with fibers thatare softer to bending than to stretching[12, 18, 19]. Thisis still subject of debate since studies in 3D have so farbeen limited to small systems[18] or to networks with highconnectivities[6, 19].To provide insight in network mechanics at densitiesranging from the rigidity percolation point to the high-molecular weight limit, we develop a 3D lattice-basedfiber network model with binary cross-links. Our net-work’s consists of straight fibers organized geometricallyon a face centered cubic (FCC) lattice. However, we limitthe maximum coordination number to four by randomlychoosing three independent pairs of cross-linked fibers outof the six fibers crossing at every vertex. Although the dif-ferent binary cross-links may overlap geometrically, theydo not constrain each other[21] (inset Fig. 2a). Therefore, we term this the Phantom FCC (Ph-FCC) lattice. Thismodel is similar to a generalized Kagome lattice in 3D[20],although it has a lower symmetry than the Ph-FCC lat-tice. By randomly cutting bonds with a probability, 1 − p ,we tune the average molecular weight, L = (cid:96) / (1 − p ),where (cid:96) is the distance between vertices[6, 21].The elastic energy of the 3D Ph-FCC network involvesstretching and bending contributions of the constituentfibers, characterized by their stretching modulus µ andbending rigidity κ . Each lattice vertex consists of 3 in-dependent freely-hinging binary cross-links ranked by h .For small displacements, denoted by u hi , the stretchingenergy of the network is expressed as E S = 12 µ(cid:96) (cid:88) h =1 (cid:88) (cid:104) ij (cid:105) g hij (cid:0) u hij · ˆr ij (cid:1) , (1)where the second sum extends over neighboring pairs ofvertices, u hij = u hj − u hi and ˆr ij is the bond direction inthe undeformed lattice. Bond-dilution is implemented bysetting g hij = 1 for present bonds and g hij = 0 for removedbonds. Fibers form straight chains that resist angulardeflections, leading to a total bending energy[6, 14], E B = 12 κ(cid:96)
30 3 (cid:88) h =1 (cid:88) (cid:104) ijk (cid:105) g hij g hjk (cid:2)(cid:0) u hij − u hjk (cid:1) × ˆr ij (cid:3) . (2)Since the cross-links themselves do not provide a torsionalstiffness, the second sum only extends over coaxial nearestneighbor triplets along the same fiber. The shear mod-ulus, G , is determined numerically by applying a shearstrain along the 111-plane with Lees-Edwards periodicboundary conditions and energy minimizations are per-formed using a conjugate gradient algorithm.Crucially, to avoid effects due to filaments spanning thenetwork, thereby making unphysical stretch contributionsto the elasticity of the sample, at least one bond is re-moved along every fiber. Our network sizes range from W = 20 to 150 unit cells, with up to three times thatmany cross-links. Due to the finite system size, this modelcan only approach z = 4 asymptotically from below. Linear regime
We find that these networks have a fi-nite shear rigidity only if κ >
0, even though the perfect( z = 4) Ph-FCC lattice deforms affinely and has a fi-nite shear modulus for κ = 0. Nonetheless, over a broadrange of L , the system can be bending dominated G ∼ κ (low κ ), or stretching dominated, G ∼ µ (high κ ), asshown in Fig. 2a. Interestingly, there appear to be twodistinct regimes well above the rigidity percolation point:a bending-dominated regime where G depends stronglyon L (low κ and L ) and an L - and κ -independent affinestretching regime (high κ and L ).These results can be understood as follows. In the high- κ limit, the system deforms affinely, with a shear modulus G A ∼ µ(cid:96) z . However, in the critical regime—controlledby the bending rigidity percolation point z b — G vanishescontinuously with ∆ z = z − z b [6, 9, 11, 20] as G cs ∼ µ(cid:96) | ∆ z | f , G cb ∼ κ(cid:96) | ∆ z | f , (3)for high and low κ , respectively. We find z b ≈ . f ≈ .
65 for a system size W = 30 , as demonstratedin the inset in Fig. 3 by showing that G | ∆ z | − f /κ reachesa plateau for low values of ∆ z . The rigidity threshold issimilar to observations in prior 3D models[18], although f is considerably lower here, which is more consistentwith findings on the generalized 3D Kagome lattice[20].The rigidity threshold can be estimated using a Maxwellcounting argument[5, 6, 18]; this connectivity-thresholdoccurs when per cross-link the number of stretching con-straints, n b z/
4, and bending constraints, n b ( d − z / d . Here,the number of bonds per cross-link n b = 2 in the undi-luted network ( z = 4). This yields z b ≈ .
6, in reasonableagreement with the numerical results.Since the CF isostatic point lies far beyond the physicalconnectivity range of this model, the naive expectationwould be that this percolation regime extends over thewhole range z <
4. This would imply bending dominatednetwork elasticity for low κ , such that G (cid:28) G A as z → z = 4 toaffine stretching dominated behavior due to fibers thatspan the whole network. However, this argument ignorespossible effects due to filament length. In networks ofstraight fibers with binary interactions, the average fiberlength diverges as z → L may lead to nonaffinedisplacements over greater length scales[12]. The effectsof high L on the deformation field have been discussedin the context of 2D Mikado networks using both scalingarguments[10, 11] and floppy mode theory[12].Here, we investigate the effects of molecular weighton the nonaffine deformation field and their implicationsfor the mechanics of 3D fiber networks. Network nodesalong a fiber can undergo independent nonaffine defor-mations scaling as γL to avoid stretching of the otherfibers to which they are connected. This direct scalingof nonaffine displacements with L constitutes one of thecentral assumptions of floppy mode theory that was ap-plied to Mikado networks[12]. To test this assumptionhere, we investigate the strain fluctuations using the non-affinity measure[6, 10, 22], Γ = (cid:96) γ (cid:10) ( δ u NA ) (cid:11) , where δ u NA = u − u A denotes the nonaffine displacement ofa cross-link and the brackets represent a network aver-age. This nonaffinity measure exhibits a cusp at the bend-ing rigidity percolation point, reflecting the criticality ofthe network’s mechanics in this regime[6, 7], as shown inFig. 2b. Furthermore, there appears to be a regime forsufficiently low κ where Γ ∼ L independent κ , lendingcredence to the basic assumption that δu NA ∼ Lγ [12]. −3 −2 −1 −6 −5 −4 −3 −2 −1 −3 −2 −1 −3 −2 −1 FIG. 3: (Color online) The shear modulus scaled with itsaffine prediction versus L scaled with λ NA = (cid:96) /(cid:96) b for vari-ous values of κ in units of µ(cid:96) . The open symbols indicatedata ranges in the rigidity percolation regime where we ob-serve different scaling. The inset shows G scaled with | ∆ z | f as a function of | ∆ z | . Such nonaffine deformations imply a bending regime, G LC ∼ κ(cid:96) (cid:18) δu NA (cid:96) (cid:19) ∼ κ(cid:96) L . (4)This prediction for the L -dependence of G is in fact bornout by the numerical data, as shown in Fig. 2a. This hasthe important implication that as L becomes large, G LC will exceed the affine shear modulus, and hence, an affinenetwork deformation becomes more favorable. Thus, aNA-A transition[10, 14, 15] occurs when G LC (cid:39) G A ; thisis satisfied when L becomes comparable to the length scale λ NA = (cid:96) /(cid:96) b , (5)where (cid:96) b = (cid:112) κ/µ . Indeed, by plotting G/G A as a func-tion of L/λ NA we find a good collapse of the data to auniversal curve, for which G/G A ≈ L/λ NA > ∼ λ NA ∼ (cid:96) − αb , with α ≈ . − . L limit for the Mikado model. nonlinear regime The length-controlled bending me-chanics also has important implications for the nonlin-ear elasticity of 3D fiber networks. Even in a bendingdominated regime, stretching modes are excited at finitenetwork deformations[13], but to a higher order in theapplied strain[7, 12, 16]. Specifically, assuming length-controlled nonaffine deformations, a transverse bend withan amplitude ∼ γL results in a stretch energy in the as-sociated bond, δE S ∼ µ(cid:15) (cid:96) , where (cid:15) ∼ ( γL/(cid:96) ) + O ( γ ).The onset of nonlinear network elasticity due to this ef-fect occurs at a strain γ , at which δE S becomes compa-rable to the leading order bending contribution, δE B ∼ −1 −4 −3 −2 −1 −1 −4 −3 −2 −1 a)b) linearbending nonlinear linearstretching FIG. 4: (Color online) a) The differential modulus scaled withits affine prediction together with γ and γ A for networks with L = 18 .
3. b) Scaling predictions (solid lines) and numericaldata for γ (open symbols) and γ A (closed symbols). L ismeasured in units of (cid:96) and κ in units of µ(cid:96) . The red andgreen points in the lower panel correspond to κ = 10 − and κ = 10 − , respectively. κL γ /(cid:96) c . This stiffening saturates at a strain γ A , set bythe condition δE B + δE S ∼ µ(cid:96) γ , at which the network’sresponse becomes affine. Thus, the onset and completionof the stiffening regime are expected to scale as γ ∼ (cid:96) b L , and γ A ∼ (cid:96) L (cid:113) − L (cid:96) b /(cid:96) . (6)These predictions are in excellent agreement with the nu-merical results over a broad range of L and κ , as shownin Fig. 4. Interestingly, the nonlinear threshold γ corre-sponds to a stress σ ∼ κ / L , which is distinct from theEuler buckling threshold [13] and the stiffening thresholdin sub-isostatic spring networks [7, 23].The Phantom FCC model developed here, provides apowerful numerical model to probe the mechanics of 3Dfiber networks with large system sizes. Using this modeltogether with a scaling analysis, we have shown that eventhough the mechanical stability of 3D networks relies onthe bending resistance of the constituent fibers, surpris-ingly for any κ >
0, network mechanics becomes affineand independent of κ when L > λ NA [25]. This NA-A transition is induced because it becomes energeticallyunfavorable to accommodate length-controlled non-affinedeformations[12] at large molecular weights. Below thisNA-A transition, in the length-controlled bending regime,we expect G ∼ ρ / , for semiflexible polymers for whichthe cross-linking lengthscale is expected to scale with ρ − / [24], or G ∼ ρ for stiff polymers, where ρ is thepolymer length-density. These predictions may accountfor a recent report of G ∼ ρ . in collagen networks[17].We conjecture that the results for the length-controlledregime also apply to models with additional interactions,other than fiber bending, that stabilize the network below the CF-threshold—including next-nearest neighbor inter-actions or bond-bending interactions for cross-links thatfix a preferred bond-angle. Specifically, we expect thatsuch networks exhibit an affine high-molecular weightlimit even for arbitrarily weak additional interactions.This research was supported by the National ScienceFoundation under Grant No. NSF PHY05-51164 and byFOM/NWO. It’s a pleasure to acknowledge discussionswith E. Frey, M. Das, C. Heussinger, O. Stenull and T.C. Lubensky. We also thank OS and TCL for sharingtheir manuscript (Ref. [20]) prior to publication. [1] D. A. Fletcher and R. D. Mullins, Nature , 485 (2010).[2] P. A. Janmey, S. Hvidt, J. Lamb, and T. P. Stossel, Nature , 89 (1990).[3] A. R. Bausch and K. Kroy, Nature Phys. , 231 (2006).[4] P. A. Janmey et al., Nature Materials, : 48 (2007).[5] J. C. Maxwell, Philos. Mag. , 294 (1864).[6] C. P. Broedersz, X. Mao, T. C. Lubensky and F.C. MacK-intosh. arXiv:1011.6535 (2011).[7] M. Wyart, H. Liang, A. Kabla and L. Mahadevan, Phys.Rev. Lett. , 215501 (2008).[8] A. J. Liu and S. R. Nagel, Nature , 21 (1998).[9] M. F. Thorpe, J. Non-Cryst. Solids , 355 (1983).[10] D. A. Head, A. J. Levine and F. C. MacKintosh, Phys.Rev. Lett. , 108102 (2003); Phys. Rev. E , 061907(2003).[11] J. Wilhelm and E. Frey, Phys. Rev. Lett. , 108103(2003).[12] C. Heussinger and E. Frey, Phys. Rev. Lett. , 105501(2006); C. Heussinger, B. Shaefer and E. Frey, Phys. Rev.E , 031906 (2007).[13] P. R. Onck, T. Koeman, T. van Dillen and E. van derGiessen, Phys. Rev. Lett. , 178102 (2005).[14] M. Das, F. C. MacKintosh and A. J. Levine, Phys. Rev.Lett. , 038101 (2007).[15] A. R. Missel M. Bai, W. S. Klug and A. J. Levine, Phys.Rev. E , 041907 (2010).[16] O. Lieleg et al., Phys. Rev. Lett. , 088102 (2007).[17] S. B. Lindstr¨om, D. A. Vader, A. Kulachenko, and D. A.Weitz. Phys. Rev. E, , 051905 (2010).[18] E. M. Huisman, T. van Dillen, P. R. Onck, and E. Van derGiessen, Phys. Rev. Lett. , 208103 (2007); E. M. Huis-man and T. C. Lubensky, Phys. Rev. Lett. , 088301(2011).[19] G. A. Buxton and N. Clarke, Phys. Rev. Lett. , 238103(2007).[20] O. Stenull and T. C. Lubensky, (to appear).[21] C. P. Broedersz and F. C. MacKintosh, Soft Mater , 3186(2011).[22] B. A. Didonna and T. C. Lubensky, Phys. Rev. E ,066619 (2005).[23] M. Sheinman, C. P. Broedersz and F. C. MacKintosh (toappear).[24] F. C. MacKintosh, J. K¨as and P. Janmey Phys. Rev. Lett. , 4425 (1995).[25] For thermal semiflexible polymers we expect λ NA ≈ (cid:112) (cid:96) (cid:96) pp