Fiber-based materials are prevalent around us. While microscopically these systems resemble a discrete assembly of randomly interconnected fibers, the network architecture varies from one system to another. To identify the role of the network architecture, we study here cellular and fibrous random networks in tension and compression, and in the context of large strain elasticity. We observe that, compared to cellular networks of same global parameter set, fibrous networks exhibit in tension reduced strain stiffening, reduced fiber alignment, and reduced Poisson's contraction in uniaxial tension. These effects are due to the larger number of kinematic constraints in the form of cross-links per fiber in the fibrous case. The dependence of the small strain modulus on network density is cubic in the fibrous case and quadratic in the cellular case. This difference persists when the number of cross-links per fiber in the fibrous case is rendered equal to that of the cellular case, which indicates that the different scaling is due to the higher structural disorder of the fibrous networks. The behavior of the two network types in compression is similar, although softening induced by fiber buckling and strain localization is less pronounced in the fibrous case. The contribution of transient interfiber contacts is weak in tension and important in compression.
Fibrous materials have, generally, random fiber network-like microstructure. Examples are widespread, ranging from the cytoskeleton of eukaryotic cells , the extracellular matrix and various connective tissues , as well as man-made materials such as open cell polymeric foams , rubber, and biomaterials such as mycelium . Network materials can be classified in terms of their architecture into two generic groups: cellular and fibrous. In cellular networks, fibers are arranged at the edges of polygonal cells and are connected only through their extreme end points. Such networks are identical to open cell foams. Fibrous networks are formed by randomly placing fibers in space and connecting them at places of minimum approach. Randomness in cellular networks is largely associated with the size and shape of the polygonal cells, whereas fibrous networks inherit randomness from the fiber distribution and the potentially spatially nonuniform degree of cross-linking. Fibrous networks are structurally similar to various biopolymer networks . As such, a systematic comparison of these two types of networks provides useful insights into the rich mechanics of the corresponding material systems.
Cellular networks were studied extensively over the past decades. Most of the earlier analytical and numerical studies considered simple unit cell models, such as the Kelvin cell or lattice-based cells, and investigated the small strain response [6–8]. In the small strain regime, the deformation is dominated by fiber bending and the behavior can be adequately defined by a single structural parameter—the network volume fraction (). It has been shown that both modulus () and the ultimate strength () scale with the volume fraction with exponents 2 and 1.5 respectively for open cell foams . The Poisson's ratio is observed to be independent of the volume fraction. Several numerical models have been developed to investigate the effect of realistic microstructural randomness associated with cell irregularity , nonuniform fiber cross section , and cell wall microstructure  on the overall network behavior.
The nonlinear mechanics of fibrous networks is of interest in particular in the case of biopolymer networks. Most biological filaments are semi-flexible, i.e., have large persistence length, comparable to their contour length, but are flexible enough to undergo bending deformation . Experiments show that networks of semi-flexible polymers are soft under small deformations and stiffen significantly with increasing strain. The stress–strain curve in uniaxial tension and volume-preserving shear has a characteristic “J-shape” . Strain stiffening is commonly attributed to either or both the nonlinearity of individual fiber constitutive behavior  and/or to the cooperative alignment of the network during straining [14,15]. In addition, biopolymer networks exhibit significant Poisson's contraction coupled with strong preferential alignment under stretch .
Cellular networks subjected to compression exhibit softening associated with strain localization in bands that run predominantly perpendicular to the compression axis. As bands form throughout the sample volume, the stress–strain curve develops a plateau. Further compression leads to densification and rapid strain stiffening. In the biopolymer literature, compression was studied primarily in the context of fibrin networks [17,18]. The general behavior is similar to that of cellular networks: gradual softening at small strains due to fiber bending and buckling followed by stiffening at large strains, due to densification.
The mechanical behavior of random fiber networks has been studied by modeling using two-dimensional (2D) and three-dimensional (3D) models. Mikado networks, generated by random deposition of fibers, are the most widely used models in 2D, whereas 3D networks can be generated by Voronoi tessellation , dilution of lattice structures , and self-organization . Such models can capture the key characteristics observed experimentally and provide useful insights into the effect of various geometric parameters (i.e., network density, coordination number, fiber bending and axial rigidities) on the overall mechanical behavior [20,22].
Despite this progress, the role of network architecture on its mechanical response remains poorly understood. This is due, in part, to the difficulty of describing quantitatively “the structure” of random networks. Heussinger and Frey  studied the effect of network architecture on the small strain elasticity comparing cellular and fibrous structures in 2D. They considered both thermal and athermal systems and observed that thermal networks, irrespective of architecture, are more sensitive to heterogeneity. They showed that elasticity of cellular networks can be defined by a serial arrangement of bending and stretching deformation modes, where the softer mode controls the modulus. These observations allowed establishing scaling laws for the small strain shear modulus in terms of network parameters. Licup et al.  also investigated this issue comparing 2D Mikado and lattice based networks (2D and 3D) having equivalent mean coordination number. Most importantly, they showed that nonlinear stiffening (onset strain and functional form) is independent of the geometry for 2D Mikado and 2D lattice networks, provided the comparison is performed for networks of same mean connectivity, z (z < 4). A similar conclusion is reached when comparing 2D and 3D lattice-based random networks produced by the specific generation procedure used in the study.
In this work, we compare the mechanical behavior of 3D networks with cellular and fibrous structures having equivalent network parameters. These types of structures have been compared in 2D in Ref. . We compare the small and large strain response in tension and compression and observe that, although the general behavior is qualitatively similar, the fibrous architecture is more sensitive to the variation of the density, and strain stiffens at a smaller rate than the cellular architecture. Further, we study the contribution of cross-link density and of fiber–fiber interactions to the overall mechanical behavior. The results indicate that the functional form of the relation between the network modulus and density is different in cellular and fibrous networks, and the way the structure stiffens under tension depends on the degree of cross-linking.
Modeling and Methods
We consider 3D cross-linked networks of straight fibers having uniform cross section and identical material properties. Fibers are distributed in space with random orientation and positions of their centers of mass. Fibers transmit axial forces and bending and torsional moments, and are cross-linked rigidly, i.e., the cross-links transmit both forces and moments. The mean coordination number (number of fibers emerging from a node) is in all cases . Central force networks with a similar value are subisostatic. Our models are stable due to the presence of the bending deformation mode.
The important network parameters are:
The network density (), defined as the total length of fiber per unit volume of the network. For networks made from fibers of identical length , the density is defined as where is number of fibers per unit volume. The network density is also related to the volume fraction (), as where is the fiber cross-sectional area.
The mean number of cross-links per fiber (). For a cellular network, the number of cross-links per fiber is always two (). Thereby, the number of fibers () and the number of fiber segments () are equal for cellular networks. In fibrous networks, the number of cross-links per fiber () can be varied independently and . Further, the total number of cross-links in the network () is . If is the average length of the segment defined by two successive cross-links, is also given by
The fiber material is considered linear elastic. Fibers have axial and bending rigidities, and , where is Young's modulus of the fiber and I is the axial moment of inertia of the fiber cross section. is used as the unit of stress here. It has been discussed in the literature [24–26] that the mechanical behavior of networks made from the same type of fiber does not depend on the two rigidities separately, rather it depends on their ratio, which is referred to as the bending length . The torsion deformation mode of fibers makes little contribution to the overall network mechanics.
In order to evidence the role of the architecture, we consider networks with similar parameter values. Table 1 shows geometric parameters of cellular and fibrous networks with volume fraction which are considered in the study. These cellular and fibrous (fully cross-linked) structures differ mainly in terms of the mean number of cross-links per fiber, . To investigate the effect of variable , we dilute the cross-link density from a “fully cross-linked” (cross-linked at all possible locations) fibrous network to a “sparsely cross-linked” network, by retaining the cross-links closest to the fiber ends and deleting the middle ones. We introduce a parameter p which represents the percentage of the total number of cross-links retained in given structure (relative to fully cross-linked state). p = 1 represents the fully cross-linked case and represents the degree of cross-link reduction in sparsely cross-linked cases. When p = 0.69, of such fibrous network is equal to that of the cellular network. This case is considered in order to allow the direct comparison of cellular and fibrous architectures since these structures become equivalent from the connectivity point of view to disordered cellular networks. Additionally, for each type of network, we investigate the dependence of the network response on the network volume fraction, , and bending length, , by varying these parameters in the range shown in Table 2.
It is important to observe that, when fibers have circular cross section of diameter D, cellular networks have only two characteristic lengths: the mean segment length, , which is inversely proportional to the density, and the fiber diameter, D. These form a single nondimensional parameter: . Hence, the network behavior is expected to depend exclusively on this parameter. Fibrous networks have three characteristic lengths: , , and the fiber length, . These form two nondimensional groups, for example and .
Cellular networks are generated using a Voronoi tessellation of the problem domain. Specifically, random seeds are used to generate a Voronoi tessellation and the network results by placing fibers along all edges of the resulting polyhedra. A typical cellular network is shown in Fig. 1(a). Each fiber has a cross-link at each end and coordination number (z = 4) by construction. However, in the actual networks since the Voronoi structure obtained by tessellation was modified slightly by eliminating a small number of very short segments which are always introduced by the tessellation procedure, but are inconsequential for network mechanics. The network density is controlled by the number of seed points and the fiber lengths are Poisson's distributed. The model size is selected to be large enough (approximately 15 times the mean segment length of the network) to minimize the effect of model size on model predictions. If size effects would be present, the boundary conditions used in this work would lead to stress estimates which are below the values predicted by an infinitely large model.
Fibrous networks are generated using a fiber packing algorithm described in Ref. . We use the random sequential adsorption technique to generate sparse fiber assemblies and then utilize dynamic finite element simulations to bring the fibers within crosslink-able distance. During the finite element simulation we discretize each fiber into multiple Timoshenko beam elements and enforce surface-based contact constraints between fiber-to-fiber surfaces in order to avoid interpenetration. These packing simulations are performed using the abaqus/explicit solver (version 6.13-1) . The compacted fiber assembly is then transformed into a network structure by introducing cross-links at all sites where the interfiber distance is smaller than , where D is the fiber diameter. The mean coordination number is also 3.8 for this system. After the cross-linking process, isolated fibers may remain unconnected to the main network. These, as well as the dangling ends of fibers, are removed. If this reduces the density below the desired value, new fibers are added and the procedure is repeated until the target density is reached. The remaining structure is a fully connected graph (Fig. 1(b)). The simulation cell size is approximately . We considered models of smaller and larger dimensions and observe that this size is sufficient to eliminate model size effects associated with .
The generated cellular and fibrous networks are discretized using multiple Timoshenko beam elements (B32 quadratic elements in abaqus version 6.13-1) per fiber so as to maintain an element aspect ratio of 5. In the fibrous network case, cross-links are modeled using rigid connector elements which transmit both moments and forces. Both cellular and fibrous networks become subisostatic and their small strain stiffness vanishes if the cross-links are rendered pin joints (transmit only forces from fiber to fiber). Fibrous networks remain rigid if the cross-links are represented as rotating links, which transmit only forces between fibers, but allow forces and moments to be transmitted along fibers. In 2D models of Mikado networks, it was seen that the difference between rotating and welded cross-links is small and the two models lead to the same dependence of the network modulus on network parameters . In this work, we tested that 3D fibrous networks behave similarly. To this end, we compared the stress–strain curves of models with welded and rotating cross-links and observe that the difference is small.
The models are subjected to uniaxial tension and compression. In tension, displacements are applied on two opposing faces of the cubic model in the direction of loading. The lateral faces of the model, which are not subjected to displacement boundary condition, are constrained to remain planar, but allowed to move in the direction of their normal under traction free condition. This arrangement facilitates the computation of the transverse strains. All other degrees-of-freedom on the model surface are free. Networks are compressed uniaxially by placing two rigid surfaces at the top and bottom boundaries of the models and imposing relative displacements on these surfaces. Fibers cannot penetrate through these loading planes during model deformation. The solution is obtained using the finite element solver abaqus/explicit (version 6.13-1). Additionally, surface-based contact between fibers and between fibers and loading surfaces is used in compression simulations in order to prevent unphysical fiber crossing. The general contact algorithm implemented in abaqus is used for this purpose.
The stress and strain measures used in this work are the nominal stress () and the stretch ratio (). The nominal stress is computed as the total reaction force in the loading direction, at nodes located on surfaces with imposed displacements, divided by the initial model cross-sectional area. The Poisson's ratio of the fiber material has no effect on the mechanics of the network and is considered to be 0.3 in this work.
Results and Discussion
Generic Stress–Strain Behavior in Tension and Compression.
Figure 2(a) shows the tensile stress–strain response of cellular and fibrous networks for three network volume fractions, %. The stress is normalized by the fiber modulus (). Models correspond to the parameter set in Table 1 and hence differ with respect to their architecture and the mean number of cross-links per fiber ( for cellular and for fibrous). The generic behavior of the two types of networks is similar to that of other fibrous systems reported in literature [30,31]. Three distinct regimes can be observed in uniaxial tension. In regime I the behavior is linear elastic, with effective stiffness . This regime occurs because geometric nonlinearity is vanishingly small at small strains and the behavior of individual fibers is linear elastic. Regime I ends at for both networks. Regime II is characterized by fiber re-orientation in the loading direction. A subset of fibers which are initially oriented roughly perpendicular to the loading direction deform significantly in bending to allow the alignment of the other fibers. This strong geometric nonlinearity causes global strain stiffening. A strong Poisson's contraction associated with this process is also observed. This regime ends at about . In regime III, the stress–strain curve is again linear, this time due to the formation of stress paths, i.e., chains of fiber segments which are almost straight and connect the opposite loaded faces of the model. These chains carry most of the applied load. Once they form, further network structural evolution is inhibited.
Compressive stress–strain curves for three volume fractions of are shown in Fig. 2(b) for both types of network. Three regimes are observed in compression as well. Regime I is linear elastic with an effective stiffness equal to that measured in tension, . This regime ends at . Fiber buckling and rearrangement takes place in regime II, while the network loses global carrying capacity and softens. Cellular materials also exhibit a similar plateau regime, which is attributed to global localization in the form of a collapse band that advances across the gauge of the sample . If the collapse band is well defined and of small width, the stress remains essentially constant during band propagation. This type of behavior is observed in cellular materials with periodic structure and no (or small) imperfections. In both types of networks, the buckling collapse is not spatially localized due to the randomness of the structure. This causes the stress to increase continuously during regime II. No advancing localization band occurs, but small scale collapse regions which do not percolate to form a global band are observed. This softening effect is more pronounced in cellular networks due to their less disordered, cellular architecture. Regime III corresponds to rapid stiffening and is observed beyond . This is associated with densification and formation of a large number of fiber–fiber contacts. This issue is discussed further in Sec. 3.5.
This generic comparison of the behavior of the two types of networks shows that for the same network parameters () the fibrous network is stiffer than the cellular network, but the transition between regimes takes place in the same strain range.
Small Strain Modulus.
While the general comparison discussed in Sec. 3.1 seems to indicate that the behavior of the two types of networks is similar, a detailed analysis of the parameters describing the stress–strain curves and their dependence on network parameters indicates substantial differences. We focus first on the small strain modulus, , and its dependence on network density () and fiber bending length (). These parameters were varied in a broad range ( was varied by a factor of 5 and by 3 orders of magnitude, Table 2) and the resulting network stiffness values, , are shown in Fig. 3 in a nondimensional form. The vertical axis shows normalized by . Note that is the expected modulus if the deformation would be affine. The nondimensional group is used on the horizontal axis for cellular networks. A combination of the two nondimensional groups and is used on the horizontal axis for the fibrous case. With this normalization, the data collapse onto a “master curve.” Similar curves defining the relationship between network parameters and the small strain modulus have been presented in the literature for other network architectures [29,33].
As previously seen for other networks [29,33], two distinct regimes are observed (marked by the vertical dashed line in Fig. 3). At low value of and , the data points define a straight line of slope 1. For cellular networks, this leads to the scaling relationship- which indicates that modulus varies quadratically with the density (or the volume fraction) and the small strain behavior is dominated by fiber bending (hence, ). A similar scaling is reported for open cell foams .
For fibrous networks in the nonaffine regime (low value of and ), the scaling relationship results , which demonstrates that fibrous network deformation is also bending dominated in the small strain regime I, but it is more sensitive to the variation of density ( as opposed to for cellular networks). can be also written as , which evidences the contribution of the volume fraction and of the two nondimensional groups of network parameters. The only other possible dependence of is . We verified that using this group on the horizontal axis in Fig. 3(b) does not lead to proper data collapse. A similar scaling of with the density was observed experimentally for collagen gels  and aegagropilae fiber network  which have similar fibrous network like microstructure .
At large values of and , the network modulus () converges to a plateau defined by . In this limit, the network deformation is dominated by the axial deformation mode of fibers and is predominantly affine. As expected based on analytic models of affine deformation , the scaling function of modulus () is independent of architecture in this regime. The theoretical volume fraction at which the network transitions from the bending-dominated, nonaffine regime to the axially dominated affine regime is approximately = 30% . The present simulation results also indicate that the transition occurs at and for cellular and fibrous networks, respectively.
Fibrous networks have multiple cross-links per fiber. Therefore, it is of interest to determine to what extent this geometric feature affects the master plot. As discussed in Sec. 2, we analyze this effect by reducing the cross-link density from the fully cross-linked state (p = 1). Cellular networks and fibrous networks with p = 0.69 are equivalent in terms of , Table 1. The only difference between the cellular networks and the fibrous network with p = 0.69 is the increased disorder of the structure. The variation of the small strain modulus with network parameters for the fibrous network with p = 0.69 is shown in Fig. 3(b). The network modulus decreases when p decreases, but the functional form of the dependence on density and , remains unchanged. This indicates that the stronger dependence of on observed in the case of fibrous networks is related to the structure and not with the increased number of cross-links per fiber. This observation and the scaling relations and for fibrous and cellular networks, respectively, are the first important results of this work.
Nonlinear Behavior and Structural Evolution.
To investigate the nonlinear stiffening of the two types of networks under tension, we evaluate the tangent modulus () as a function of nominal stress () for various network volume fractions (Fig. 4). The vertical axis is normalized with the small strain modulus () and the horizontal axis is normalized by the stress value at onset of stiffening () in order to render the curves directly comparable. The curves without this normalization are shown in the inset. The three regimes of the stress–strain curve are recovered in this representation; the transitions between regimes are indicated by red arrows. The strain stiffening regime II of Fig. 2 corresponds to the branch bounded by the two red arrows. Cellular networks exhibit a slope of 1 in this log-log plot (), which indicates exponential stiffening () in regime II. However, fibrous networks exhibit a different scaling, , which indicates quadratic stiffening, i.e., . This is the second important result of this work.
Experimental observations on collagen gels and soft connective tissues [36,37] indicate that the tangent stiffness is proportional to the stress in this regime which is consistent with the cellular network behavior. Stiffening in this regime is primarily associated with the rapid re-orientation of load carrying fibers. We interpret this result based on the observation that in the fibrous network case additional constraints are imposed through random cross-linking along each fiber. As such, fiber re-orientation is more difficult and the stiffening rate is consequently smaller.
To test this interpretation, we compare the stiffening rate of fibrous networks of same volume fraction, but different values of p. The results are shown in Fig. 4(c). The strain stiffening rate is seen to depend markedly on p. The slope varies smoothly from 1/2 to 1 as p increases from 0.69 to 1. The analysis provides support to the assumption that reducing the degree of internal constraints in the network favors fiber alignment which is associated with increased strain stiffening. From a design perspective, the results indicate that networks with strain stiffening anywhere between exponential and quadratic can be produced by modifying the number of cross-links per fiber, while keeping all other parameters (including and <z>) constant. The observation that the degree of strain stiffening depends on the number of cross-links per fiber and is independent of the structural disorder (compare fibrous networks with p = 0.69 with cellular networks) is the second important result of this work.
Poisson's Effect and Fiber Alignment.
In the linear elastic regime I, is smaller than 0.5 for both networks and does not vary much with strain or network volume fraction. In regime II, it increases nonlinearly with strain. The curves reach a peak at strains close to the onset of regime III and the Poisson's effect decreases in magnitude gradually with further straining, in regime III. The occurrence of the peak and the underlying mechanics have been discussed in Ref. , where a comparison with experimental data from tests on various connective tissue is presented. Regime III corresponds to the formation of stress paths and hence to strong network alignment. Once this strain-induced structural re-organization has been achieved, further lateral contraction is limited by the fibers which are not part of the stress paths and run, more or less, perpendicular to these. Continued Poisson's contraction requires further compressive/bending deformation of these filaments, which, however, have been already strongly distorted during the process of stress path formation. It is seen that the Poisson's effect is weaker in the fibrous network, but the overall trends are similar to those observed in cellular networks. We also observe that in compression, the incremental Poisson's ratio is essentially constant (and smaller than 0.5) during deformation.
where is the angle between the stretch direction and the axis of a fiber segment, while indicates ensemble averaging over all fibers in the network. corresponds to randomly oriented fibers, whereas indicates complete alignment in the stretch direction. In the undeformed state, values for both networks are close to zero, indicating that fibers are randomly oriented in space and that these networks are isotropic in the unloaded state. The variation of with the stretch during uniaxial tensile testing is shown in Fig. 5(b). Cellular networks exhibit more pronounced alignment than the fibrous networks. This is, again, due to the fact that fibers are less likely to align in the loading direction in presence of multiple cross-links per fiber (fibrous case). However, the orientation index follows the same trend in both network types.
Effect of Interfiber Contacts.
It is of interest to determine to what extent interfiber contacts that form dynamically during loading influence the overall mechanical behavior of the network. This is easily achievable in our models by comparing responses obtained with and without imposing the nonpenetration restriction at interfiber contacts. The interfiber contacts have no effect on the stress–strain curve in tension. This is actually expected since these rather open structures with low are composed from straight, unentangled fibers. In tension, contacts form only at large strains, in regime III, and have a weak effect on the lateral contraction, but have no effect on the stress–strain curve.
Interfiber contacts play a central role in compression in both network types. Figures 6(a) and 6(b) show compression stress–strain curves for the two network types, evaluated while accounting for interfiber contacts and without this provision. The difference between the two cases is negligible for in cellular networks. In fibrous networks, this limit is smaller, . Beyond this limit, the case with contacts presents continuous stiffening due to the rapid formation of new interfiber contacts (Fig. 6(b)). The case without contacts exhibits softening beyond this threshold strain. The more pronounced increase of the number of contacts in the fibrous networks is expected, given the more pronounced structural disorder in these networks.
Figure 6(b) shows the variation of the density of interfiber contacts, , with stretch during compression for both types of networks and for . increases nonlinearly with strain for both networks, but the rate of increase is higher in fibrous networks.
In Eq. (4), , the initial volume fraction, and . Expression (4) is plotted in Fig. 6(b) along with the numerical data. Toll's equation predicts the evolution of the number of contacts accurately, especially in the small strain regime.
Figure 7 shows the variation of the mean of the distribution of contact forces during uniaxial compression of cellular and fibrous networks of same volume fraction (). For , the mean contact force increases at the same (small) rate in the two networks. The total number of contacts in this regime is small (Fig. 6(b)). Further compression leads to a rapid rise of the mean contact force in the fibrous network, while the increase is much more gradual in the cellular network. At the same time, the total number of contacts is approximately identical in the two networks at any stage of the deformation, Fig. 6(b). This indicates the more important role of contacts in the fibrous case compared with the cellular case, which is also demonstrated by the data in Fig. 6(a).
The mechanical behavior of cellular and fibrous networks is compared in this study. To this end, networks of same model parameters are constructed and tested in uniaxial tension and compression. Although the general behavior is similar, we observe important differences. Specifically, we observe that the strain stiffening of cellular networks is more pronounced than that of fibrous networks. This is due to the higher number of cross-links per fiber in the fibrous case which impose kinematic restrictions that limit the ability of the network to re-orient under strain. Furthermore, the small strain stiffness of fibrous and cellular networks increases with the density as and , respectively. This difference persists when the number of cross-links per fiber in the fibrous case is reduced to 2, i.e., it is rendered equal to that of the cellular case. Both networks exhibit a strong Poisson's effect associated with fiber preferential alignment. The less kinematically constrained cellular networks align more and exhibit a more pronounced Poisson's effect at all strains. Deformation-induced interfiber contacts are important only in compression when induce significant stiffening.
U.S. National Science Foundation (NSF), Division of Civil, Mechanical and Manufacturing Innovation (Grant Nos. CMMI 1362234 and CMMI-1634328).