An Adaptive QuasiContinuum Approach for Modeling Fracture in Networked Materials: Application to Modeling of Polymer Networks
Abstract
Materials with networklike microstructure, including polymers, are the backbone for many natural and humanmade materials such as gels, biological tissues, metamaterials, and rubbers. Fracture processes in these networked materials are intrinsically multiscale, and it is computationally prohibitive to adopt a fully discrete approach for large scale systems. To overcome such a challenge, we introduce an adaptive numerical algorithm for modeling fracture in this class of materials, with a primary application to polymer networks, using an extended version of the QuasiContinuum method that accounts for both material and geometric nonlinearities. In regions of high interest, for example near crack tips, explicit representation of the local topology is retained where each polymer chain is idealized using the worm like chain model. Away from these imperfections, the degrees of freedom are limited to a fraction of the network nodes and the network structure is computationally homogenized, using the micromacro energy consistency condition, to yield an anisotropic material tensor consistent with the underlying network structure. A nonlinear finite element framework including both material and geometric nonlinearities is used to solve the system where dynamic adaptivity allows transition between the continuum and discrete scales. The method enables accurate modelling of crack propagation without a priori constraint on the fracture energy while maintaining the influence of largescale elastic loading in the bulk. We demonstrate the accuracy and efficiency of the method by applying it to study the fracture in different examples of network structures. We further use the method to investigate the effects of network topology and disorder on its fracture characteristics. We discuss the implications of our method for multiscale analysis of fracture in networked material as they arise in different applications in biology and engineering.
Keywords: QuasiContinuum Method, Polymer Networks, Fracture
1 Introduction
Polymers are the building blocks for many natural and artificial materials. The load transfer system in polymeric materials may be abstracted as a complex network of crosslinked nonlinear chains [1]. The topological properties of polymer networks, such as local connectivity, crosslinking density, bond types, and sacrificial bonds and hidden length, may dramatically affect their mechanical response and fracture proprieties [2, 3]. Most current experimental techniques do not allow direct visualization of the deformation and damage on the network micro scale, thus, the relation between the topology and mechanical behavior of these materials remains elusive. Multiscale numerical models are thus needed to bridge this gap [2].
Modeling of soft polymeric materials has been mainly approached using continuum theories, including linear elasticity, hyperelasticity, viscoelasticity, viscoplasticity and poroelasticity [4, 5, 6, 7, 8]. Continuum representation of materials has been used in traditional solid mechanics, where the material is assumed to be a continuous medium with its response described mathematically by a set of constitutive laws. Through experiments, the parameters of the constitutive laws may be determined and then utilized to solve the boundary value problems using either analytical or computational approaches. These techniques enabled significant progress in modeling the elasticity and distributed damage behavior of polymeric materials with and without infused fluids. However, these continuum methods are usually independent of the size scale of the network structures and they may not capture localized phenomena without enrichment, which in many cases is either expensive or is based on idealized models of microstructure. In particular, when it comes to fracture, a phenomenon that depends critically on the local conditions in the vicinity of the propagating cracks, continuum methods are not suited for predicting the effect of the network local topology on its fracture resistance. Furthermore, continuum methods may also fall short in guiding the network geometric design for enhanced fracture toughness.
An exception in the recent years has been the emergence of phase field models which is increasingly used to simulate fracture in a wide class of materials including gels and rubbers [9, 10, 11]. The method may be used to simulate curved crack paths, crack kinking and branching, and crack front segmentation. Phase filed methods are based on gradient enhanced damage models where a parameter is introduced to distinguish between fully intact material and fully damaged material in a gradual manner [12]. The width of the damage zone depends on a length scale that is introduced primarily for regularization purposes but has been shown theoretically that as this length scale goes to zero, Griffith fracture criterion [13] is recovered. The damage variable changes from 1 (fully intact) to 0 (Fully damaged) over this length scale. This length scale must be carefully chosen to regularize the strainsoftening behavior during the fracture process, and to avoid mesh dependencyrelated issues during finite element simulations. The choice of this lengths scale for materials with microstructure, which possess physical length scales of their own, is not clear apriori and is left, in most cases, as a tunable parameter to fit experimental observations [14].
Recently, thanks to advances in imaging, computational power, and nano/micro fabrication technologies, the mechanics of materials takes on another challenge of modeling fundamental material behavior starting from the atomic scale [15]. For the accurate representation of the material response at this scale, the material must be modeled using discrete simulations. Adapting these discrete models to polymer networks where each polymer chain is modeled explicitly [16, 17], it is possible to accurately describe localized phenomena such as fracture and cavitation. Furthermore, these discrete models provide a pathway for exploring an explicit connection between network design and its mechanical response. However, it is computationally prohibitive to adopt a full discrete approach for large scale networks which limit the applicability of these models.
In this paper, We introduce a new adaptive numerical algorithm for modeling fracture in polymer networks using an extended version of the QuasiContinuum method [18, 19, 20] that accounts for both the nonlinear elastic nature of the polymer chains as well as the geometric nonlinearity associated with their potentially large deformation. In regions of high interest where deformation gradients are non uniform, for example near crack tips, explicit representation of the local topology is retained where each polymer chain is idealized using the worm like chain model. Away from these imperfections where the deformation gradient is sufficiently uniform, the network structure is computationally homogenized, using an equivalence of the microscale and macroscale incremental work, to yield an anisotropic material tensor consistent with the underlying network structure. Thus at any instant in time, only a fraction of the network nodes is solved. Dynamic adaptivity allows efficient transition between the two resolutions [20]. Our goal is to develop a method that enables modeling crack propagation without apriori constraint on the fracture energy or the need to assume phenomenological length scales, while maintaining the influence of largescale elastic loading in the bulk.
The remainder of the paper is organized as follows. In Section 2 we present an overview of the QC method. We introduce the problem formulation, the homogenization technique, and the numerical implementation algorithms used in this study in Sections 3 and 4. We verify the numerical algorithm for modeling both pristine and cracked samples in Section 5. In Section 6, we demonstrate the efficiency of the method by applying it to study the fracture of in model systems that are prohibitive for a fullfield discrete approach. We further use the method to study the effects of network topology on its fracture resistance. Finally, We discuss the method implications for the analysis of networked material systems in Section 7.
2 QuasiContinuum method overview:
Tadmor, Ortiz, and Phillips [18] proposed the Quasicontinuum (QC) method in 1996 to overcome the scale limitations of atomic simulations. The original QC method is a computational technique to model an atomistic system without the need to treat all atoms in the domain explicitly. The QC provides a framework where degrees of freedom are judiciously eliminated and force or energy calculations are approximated at areas of low interest, whereas exact discrete representation is used at the area of high interest such as crack tips or dislocation zones. The QC method also provides an adaptive framework for the fully resolved critical part to evolve during the simulation to balance between computational cost and error estimates.
The QC method finds the deformation field that minimizes the system potential energy while achieving the following three criteria [19]:

The number of degrees of freedom is reduced by limiting the degrees of freedoms to a small fraction of the nodes, called representative nodes (repnodes), but the full discrete description is retained in high interest regions, around crack tip for example. Hence, the computational cost of solving the system is significantly reduced.

The computation of the total energy is accurately approximated without the need to explicitly compute the site energy of all the atoms. Hence, the computational cost associated with system assembly is significantly reduced.

The fully resolved critical regions can evolve with the deformation, during the simulation. Thus, the method achieves efficiency by using adaptivity to reduce the fully resolved domain as much as possible.
The QC method applied to atomistic systems has been successfully used to investigate many localized phenomena, such as Nanoindentation, cracktip deformation, deformation and fracture of grain boundaries, and dislocation interactions [21, 22, 23].
While the QC method was originally introduced for atomic simulations, it may be extended to any discrete systems including networks or lattices. For example, Mikes and Jirasek [24] extended the quasi continuum method to irregular linear elastic lattices with only axial interactions and small deformations. Beex et al. [25] extended the method to small deformation uniform linear elastic beamlike lattices for both in plane and out of plane deformation, and Phlipot et al. [26] used the method to model linear elastic periodic beam lattices with corotational framework. The flexibility of the framework makes it a suitable candidate for applications to irregular polymer networks with both material and geometric nonlinearities, which is the focus of this study.
3 QuasiContinuum for irregular polymer networks:
In this section, we discuss the problem setup and the finite element formulation for both the discrete part and the homogenized part of the network model. For the discrete part, we use a nonlinear 1D finite element with finite deformation and rotation. For the homogenized part, we use 2D nonlinear triangular finite element with finite deformations. The material properties for the 2D elements are derived by applying homogenization rule based on the micromacro energy consistency condition that conserves the variation of local work between the macro (2D homogenized elements) and micro (underlying network) scales.
3.1 Finite Element Formulation: 2D homogenized elements
The strong form for the boundary value problem in the absence of inertia effects is given by:
(1) 
Where is the PiolaKirchhoff stress tensor, is the material point coordinate in the reference configuration, is the body force, and is the material density at the reference configuration.
The boundary conditions are given by:
(2) 
Where
Where is the applied traction on boundary , is the unit vector normal to the boundary surface in the reference configuration, is the prescribed displacement on boundary , and is the boundary surface.
The weak form in a total Lagrangian formulation may be written as:
(3) 
Where is the deformation gradient, and where the space of kinematically admissible displacements with the requirement that the displacements vanish on displacement boundaries. The finite element approximation is derived in terms of the objective PiolaKirchhoff stress tensor . After lineraization and making use of Voigt notation, the internal force vector and material tangent matrix are given by:
(4) 
Where
(5) 
Where is the straindisplacement matrix. The goal of the homogenization is to derive the macroscopic stress and material tangent tensors that map the microscale behavior with sufficient accuracy. These quantities are then provided to the nonlinear finite element framework for solving the system.
3.2 Homogenization:
3.2.1 Micromacro energy consistency condition:
We establish the micromacro scale transition relation based on the HillMandel condition [27, 28]. This condition requires the volume average of the variation of work on the micro level to be equal to the variation of local work on the macroscale:
(6) 
In terms of macro deformation gradient tensor and the first PiolaKirchhoff stress tensor, the condition reads [29]:
(7) 
Where and are the force and the length of each link in the current configuration, is the number of links. Using , where is the link length in the original configuration, Eq.7 becomes:
(8) 
Where the subscript denotes the macro level tensors, and denotes the micro level tensors. Using Voigt assumption (i.e ). The homogenized PiolaKirchhoff stress is given by:
(9) 
To calculate the tangent tensor, we start by taking the variation of :
(10) 
Substituting and rearranging:
(11) 
Where the subscript LT denotes the transpose of the left two indices. Again, using Voigt assumption, it follows that:
(12) 
Where is the homogenized order first elasticity tensor. It must be noted that the homogenized stress tensor is a twopoint tensor (i.e not symmetric), and the elasticity tensor has only major symmetry.
It is more convenient to use the objective symmetric PiolaKirchhoff Stress tensor and the corresponding elasticity tensor. By definition, the PiolaKirchhoff Stress tensor is given by:
(13) 
Using indicial notation and substituting for from Eq.9
(14) 
After using Voigt assumption:
(15) 
We note that the quantity is the (fictitious) force vector corresponding to the PK2 stress.
To derive the elasticity tensor relating the PK2 tensor and the Green strain, we start by varying the relation between PK1 and PK2:
(16) 
After substituting from Eq.12 and using Voigt assumption, the first term may then be written as: (we drop the subscript p for convenience)
(17) 
Simplifying and using , :
(18) 
The second term is written as:
(19) 
After simplification and some indices manipulation:
(20) 
(21) 
(22) 
Hence:
(23) 
Where is the homogenized order second elasticity tensor.
3.3 Assembling the tensors:
Since the links in the network have axial deformations only, the deformation gradient for each link is written as:
(24) 
Where is the stretch, is the unit vector in the current configuration, and is the unit vector in the original configuration. Similarly:
(25) 
The PK1 stress tensor is then written as
(26) 
And the PK2 stress tensor:
(27) 
Whereas the material tensors are given by
(28) 
(29) 
Here for the link, , and are the magnitudes of the link force and reference length respectively. The stress tensor is symmetric, and the elasticity tensor has both major and minor symmetries.
3.4 Polymer chain constitutive law:
In this study, we model each link in the polymer network using a nonlinear elastic force elongation relation give by the worm like chain model. The forceelongation relation is given by [30]:
(30) 
Where is the force, is the endtoend distance, is the persistence length, is Boltzmann’ constant, is the temperature, and is the contour length of the polymer chain. We assume the network is force balanced at the reference configuration. The chain stiffness is given by the slope of the tangent of the force elongation curve and it takes the form
(31) 
These expressions are used in Eq.27 and Eq.29 to compute the homogenized stress and elasticity tensor. As the endtoend distance approaches the contour length , the polymer chain response becomes highly nonlinear as both the force and the stiffness values go to infinity. This signals a limitation of this constitutive description which may be circumvented by accounting for strain energy of the chain in addition to its entropic energy [31]. Such correction will be investigated further in future work. In this initial study, such extreme limit is not probed as we adopt a maximum stretch failure criterion for the links. That is, we assume the link would fail if its stretch relative to the available contour length exceeds a threshold value that is predefined. Other failure criterion based on transition state theory for bond breakage [32] may also be used.
4 Numerical Implementation:
To couple the micro and macroscales, the network nodal displacements are related to the representative nodes (repnodes) displacement , see Fig.3c, through the interpolation functions:
(32) 
where is the interpolation function associated with repnode and is the number of repnodes. Using the compact support of the finite element shape functions, the displacement of each node is determined from the sum over the three vertices of the triangle containing this node. In this initial study, we use linear shape functions which turns out to perform satisfactorily as will be discussed shortly.
After evaluating the network nodal displacements, the stress and material tangent tensors are assembled as follows:
Since we use the total Lagrangian formulation, the reference configuration quantities, such as the reference length , and the reference unit vector , are calculated at the system initialization and stored instead of calculating them each time step. This leads to further computational savings. In the current work, we account for the contribution of all links in each continuum element in the homogenization procedure. We discuss possible alternative approaches in Section 7.
4.1 Automatic mesh adaptivity:
Automatic Mesh Adaptivity is a critical ingredient for the efficient implementation of the QC method as it enables the QC mesh to evolve dynamically based on the crack propagation path. In our proposed approach, we mark a 2D finite element for refinement if it meets 1 of the following criteria:

Link stretch criterion:
(33) Where is the stretch relative to the polymer chain contour length of any network link inside the QC element , and is a threshold value. The threshold value is chosen as a fraction of the failure stretch.

Deformation gradient error criterion [19]:
(34) Where is a scalar measure that quantifies the error introduced into the solution by the current density of representative nodes, is the volume of element , is the QC solution for the deformation gradient in element , is the projection of the QC solution for given by where is the shape function array, and is computed by averaging the deformation gradient in each element sharing a given repnode, and is a specified threshold value.
We start by constructing righttriangulated initial mesh to reduce summation errors as per [33]. Each QC mesh point is then moved to the nearest network node. All elements marked for refinement are added to a set . For each element inside . We use the standard Rivera algorithm [34], which conserves nondegeneracy, conformity, and smoothness. This algorithm tracks the longestedge propagation path (LEPP) associated with a triangular element in a backward manner. The LEPP of in a conforming triangular mesh is defined as an ordered list of triangles , ,…, such that is a neighbour by the longest edge of for each , the refinement scheme is shown in Fig.2. The algorithm, which is summarized in Algorithm 2, is then used for each and repeated till is fully resolved.
After mesh refinement, each representative node is moved to the nearest network node. All fully resolved QC elements (i.e. elements that contain only 3 network nodes) are transformed from 2D triangular elements to 1D links consistent with the underlying network structure.
4.2 Time adaptivity:
Time adaptivity is required for computational efficiency by allowing larger time steps, when possible, without compromising accuracy. In addition, time adaptivity is also critical for the accurate modeling of fracture in high interest zones where small time steps are needed. Using a constant time step may be overly restrictive if a small time step is used throughout and may lead to inaccurate results if a large time step is perpetually maintained. Time adaptivity balances computational cost and error. The simulation starts with a large time step, that accurately and efficiently compute the nonlinear elastic deformation, and then once failure initiates, the time step is reduced progressively such that further reduction in the step size won’t affect the number of elements failing at the current time step.
4.3 Nonlinear finite element framework:
The algorithm used for adaptive modeling of fracture in polymer networks is given in Algorithm 3. The homogenization is performed on the fly to provide the stress and the material tangent tensor each time step.
5 Verification
For verification of the proposed homogenization technique, we run multiple tests for both pristine and cracked samples. In each case, we verify the QC results by comparing it to the results of the fully discrete model and quantify the error in terms of the measured reaction normalized force at different levels of stretch. We define the normalized force as (Eq.30).
5.1 Uniform deformation gradient tests:
For the first test, we consider a regular network geometry shown in Fig.3a with uniform material properties. All the network links in the sample are homogenized using the procedure outlined in Section 4 . We impose an affine deformation gradient, considering both uniaxial and biaxial tension loading cases, up to a stretch value of 6. The QC results match exactly the fully discrete model in terms of the reaction force and the total strain energy values. This verifies that the homogenization scheme recovers the exact response if the microscale deformation is affine and the network structure is regular.
Next, to study the effect of disorder in the geometric and material properties on the performance of the QC method, we simulate the response of a non uniform network, as shown in Fig.3b and the zoomed in view in Fig.3c, under uniaxial tension. We introduce nonuniformity in two ways. First, we shift the position of the network nodes randomly with respect to a reference perfectly ordered lattice. The maximum shift values in and directions are and , respectively, where is the geometric nonconformity parameter. Second, we draw the polymer chain contour length randomly from a uniform distribution ranging between , where is the contour length nonconformity parameter. For each nonuniformity parameter value, ten randomly generated meshes are simulated and the average results are plotted in Fig.4. We define the relative error in the force as:
(35) 
Where is the sum of reactions in the loading direction, refers to the corresponding fully discrete simulation and refers to the QC simulations containing homogenized continuum elements only.
Fig.4a shows that as the disorder increases the error percentage from the homogenization scheme increases. However, this error decreases as the fraction of repnodes in the sample increases. This is further confirmed in Fig.4b which shows that variation of error as a function of variability in the chain lengths at two values of repnodes percentages. At low values of disorder, less than 0.2, both the homogenized and discrete models give identical results (the error is almost zero). As the disorder increases, the error increases. However, it remains less than 8 % and it decreases with mesh refinement. Fig.4c shows that the error from the homogenization depends on the stretch level. For low values of disorder, the error is negligible even at high stretch values (up to 6). As the disorder increases, the error is relatively small (few percents) at low stretch values but it rapidly increases at high stretch values reaching as high as 15 % for =0.5 at stretch value of 5. While this exercise attests to the accuracy of the homogenization scheme at low stretch values, it suggests that mesh refinement is crucial for retaining accuracy at high level of stretch. This will be accounted for automatically in the QC application by turning on dynamic adaptivity and enforcing a refinement criterion based on the local stretch level as explained in Section 4.
5.2 Nonuniform deformation gradients– Effect of refinement:
Here, we run two tests where the deformation gradient is not uniform over the domain. In the first test, we model a pristine sample with dimensions 64x64 where the bottom edge is restrained from movement in both directions and the top edge is clamped and is pulled upward. The results shown in Fig.5 using only continuum elements, homogenized through the QC procedure, show that although the force vs. stretch shows good agreement with the full discrete model, up to a stretch value of 6, there is a discrepancy in the deformed shape specially near the edges. The error distribution plot also shows that the nodal displacement error is highest in the edge elements, closest to the free lateral boundaries, where the network nodal displacements are interpolated from the edge nodes of the continuum elements. To address this problem, we designate these areas as high interest areas and model the network links in these regions explicitly. This leads to eliminating the concentrated errors near the edges. Although the repnodes ratio increased from 0.5% to around 5%, the deformed shapes and the force vs. stretch curves using both QC and fully discrete match better in this case. Furthermore, the error in the displacement distribution falls almost an order of magnitude throughout the domain. In practical applications of the QC method, the discovery of high interest areas and the change in the model resolution is done automatically using the adaptive mesh refinement algorithm.
In the second test, we use the same geometry and boundary conditions as in the first one, but the top edge now moves both upward and to the right such that the sample is subjected to both tension and shear. In this case, mesh refinement reduces the error in nodal displacement and leads to a better agreement in the force vs. stretch response between the QC and the full discrete models. The results shown in Fig.5 and Fig.6 illustrate the capability of the the QC method in bridging different scales, and also highlights the critical role of mesh adaptivity in reducing the errors specially in areas where the deformation gradient is not uniform. Furthermore, since the results in Fig.5, and Fig.6 suggest that the error in the displacement is usually high near the free surfaces, we have chosen to always use an initial QC mesh that is refined around the sample boundaries in the remaining cases presented in this paper.
5.3 Fracture test:
Here, we model crack propagation in a notched network using both fully discrete and QC models. The network dimensions are (500x250) chain units with around 250,000 degrees of freedom and the notch is 250 units long. The sample dimensions and boundary conditions are shown in the insert in Fig.7a. All links follow the worm like chain model, with forcestretch relation given by Eq.30. The chain length is drawn from a uniform random distribution with . A link is set to break if the value of the stretch exceeds 0.8. Automatic mesh adaptivty is turned on to allow the high interest area to evolve with the crack propagation.
The results in Fig.7 show that the crack path and the force vs. stretch curves are almost identical for both QC and fully discrete models. Specifically, the error in the peak force is around 1.0%, the error in the peak stretch is 1.4%, and the error in the total area under the curve is under 0.5%. The QC mesh evolution is also shown in the same figure. The first mesh Fig.7b has a rep. nodes ratio of only 2.1%, whereas the final mesh Fig.7d has 5% rep. nodes ratio. The significant reduction in the number of degrees of freedom lead to extreme computational saving where the run time of the QC method is of the order of 1 percent of time required for the fully discrete method. The computational savings are expected to increase even further as the sample size increases since most elements will be modeled as continuum elements in this case.
6 Results
In this section, we apply the QC method to investigate the effects of local topology and disorder on the network mechanical properties and its fracture characteristics. We show the advantage of the QC method in capturing those effects without the need to explicitly model every node of the network. As explained earlier, the discrete representation is limited to the zone around the crack whereas the regions away from the crack are homogenized consistently. This introduces a computationally efficient accurate representation of the fracture process.
6.1 Damage evolution in samples with different disorder parameter:
Figure.8 shows the damage evolution and the crack path for the single notched sample setup shown in Fig.7a with two different disorder parameter () values and , where is the coordination number defined as the number of connections a node may have to other nodes in the network, averaged over all the nodes.
For , which represents an almost uniform network, the crack propagates coherently perpendicular to the applied loading direction. The crack path is nearly horizontal, and the damage is localized along the crack path. For highly disordered networks with , although the dominant crack propagation direction is, on average, perpendicular to the applied loading direction, the crack path has more kinks as it propagates. Furthermore, the crack thickness, which may be taken as a basis for a material length scale in nonlocal theories such as phase field methods, is not set apriori. It emerges as part of the solution and shows nonuniformity as it tends to increase at some locations and decreases in others. Damage is not localized around the crack path, as the figure shows multiple local damage zones away from the crack tip. Some of these local damage spots coalesce with the main crack while some remain isolated. Fig.8 also shows the increase in the number of repnodes, represented by the blue dots, with crack propagation due to mesh adaptivity. The rate of increase in the number of repnodes for more disordered network is higher than that in the more uniform ones. This attests to the flexibility of the QC method in adapting to the complexity of the underlying network structure.
6.2 Force distribution in the vicinity of the crack tip:
To gain further insight into the processes controlling crack initiation and growth, we analyze the force distribution in the discrete network in the vicinity of the crack tip for different values of the coordination number and disorder parameter . As defined earlier, the coordination number is the average number of nodes directly connected to a given node and measures how densely connected the network is.
The results in Fig.9show how changing the local connectivity changes the force field around the crack altering conditions for the crack initiation, propagation, and growth direction. For example, for and , the stress field has the classical peanutlike distribution as expected for Mode I fracture in planar elasticity [35]. On the other hand, for and , the anisotropic network channel the forces in the links that are closely aligned with the loading direction. For , the force distribution shows an interesting pattern with some increase in the forces along diagonal directions running backward relative to the crack tip. For each coordination number, increasing the disorder in the chain length significantly changes the force field. In particular, the force field becomes more diffuse and loses its symmetry. Furthermore, increased disorder may lead to elevated force magnitudes in regions away from the crack tip leading to a more complex fracture pattern and potentially thicker cracks as was observed in Fig.8. Such complex force field distribution in highly anisotropic and disordered networks suggest the possibility of potential crack growth in unusual directions as has been observed recently as sideway cracks in silicone elastomers [36]. The QC methodology may enable systematic discovery of these new patterns.
6.3 Effect of disorder in the chain length:
We simulate the disorder in the network by drawing the chain length from a stochastic distribution. Fig.10 shows the response of a network with a single notch, the same setup shown in Fig.7a, when the available chain length is drawn from a random uniform distribution ranging between where is a predefined contour length. Fig.10ac show the effect of the chain length variability on the crack path. As the variability increases, the crack path becomes more meandering while propagating through the network. The QC mesh also shows higher refinement (high repnodes ratio) as the disorder increases. This is because as the range of variability in the chain length increases, some links with shorter contour length become critically stressed at locations not at the crack tip requiring further refinement in these regions.
To further investigate the effect of disorder, the average force vs. stretch curves are plotted for different disorder parameters in Fig.10d. Each curve is the average of ten different cases with the same value. The results show an increase in the stiffness and a decrease in the peak force as the disorder increases. The increase in stiffness is intrinsic in the nature of force stretch curve of the single worm like chain where the force increases rapidly as the stretch ratio approaches the stretch failure threshold. For more disordered networks, the increase in the force in the shorter links dominate over the reduction in force in the longer links. For the same reason, networks with higher variability in chain length tends to initiate failure at lower peak force, as they have a larger number of shorter links which fail faster. The results also show that the ductility (i.e. the ratio between stretch at total failure and stretch at peak force) increases with increasing . In less disordered networks, the links are more uniformly stressed and hence they fail rapidly once failure initiates. As the variability in the chain length increases, the network has more longer chains which helps the material system to survive longer under damage, this agrees with the experimental results reported in Bonn et al. [37]. Fig.10e shows that the average peak force decreases with increasing , whereas the average ductility increases. Furthermore, the uncertainty in both peak force and ductility, represented by error bars, increase with increasing the material disorder.
Figure.11 shows the effect of increasing the network disorder on the crack path of a network with symmetric and asymmetric double edge notches. As the disorder increases, crack path tends to be more nonuniform and thicker. In addition, a more extended spatial region is refined as the variability in the chain length increases, to capture the local failure of shorter chains in locations away from the crack tip. Furthermore, these results show how the two cracks emerging from the double notches interact with each other at different values. The crack path shape for the case of matches previously reported experimental and numerical observations [38, 39] .It also matches the numerical results for elastomeric materials fractured by crosslink failure using the phase field method [9]. However, our results suggests that local topology may have a profound effect on this crack path that may be challenging to capture using fully homogenized approaches. In particular, the curving of the cracks towards one another may completely disappear as the disorder increases.
6.4 Effect of coordination number Z:
In this section, we investigate the effect of network connectivity through changing he coordination number . Fig.12 shows the force vs. stretch for three different coordination numbers, , and . In constructing these networks, we keep the total weight of the network ,i.e. the sum of all polymer chain average contour lengths, as a constant. We also draw the contour length of the different chains from the same uniform random distribution with for the three cases. The results show that networks with higher coordination number has higher toughness and peak force compared to networks with lower coordination number. As the coordination number increases, more links contribute to load sharing and hence the network tends to become stiffer and sustain higher load. Furthermore, the case with higher coordination number shows a larger ductility ratio, defined as the ratio between stretch at failure and the stretch at peak force, but a slightly smaller maximum stretch. As the coordination number increases, there are, on average, more links connected per node. This enables a more efficient force redistribution among neighboring chains when one of the chains break off, leading to slower force drop and delayed damage accumulation. Similar results were obtained for nanocomposite gels with nanoparticle crosslinkers that increase the network coordination interactions [40]. As the coordination number decreases, the material becomes more compliant, and is capable of achieving slightly higher stretch values although with lower overall toughness and strength. Specifically, the area under the curve for is 1.5 times that for and the area under the curve for is twice that for .
6.5 Effect of crosslinking density:
To investigate the effect of increasing crosslinking density, defined as the number of network nodes per unit area, we simulate multiple samples with increased number of nodes and links in the the network while keeping the polymer concentration constant. This is achieved by adjusting the average chain contour length such that the sum of lengths of all chains remains constant. The results in Fig.12b shows that the peak force increases with increasing the crosslinking density up to a maximum value, then it starts decreasing. This non monotonic response is due to the reduction of the average chain length with increasing the crosslinker concentration. When the average chain length decreases, the network becomes stiffer and attains higher forces for a given value of stretch. However, as the crosslinker concentration further increases beyond a certain value, which depends on the chain length, the network becomes critically stressed as the initial stretch of the short chains becomes closer to their contour length. These short chains accumulate damage at a faster rate leading to a premature failure at lower global force levels. With increasing the available chain length, the peak force for a given concentration increases, and the critical crosslinker concentration corresponding to this peak force also increases.
7 Discussion:
We have introduced a new adaptive numerical algorithm for modeling damage and fracture in networked materials in general, and polymer networks in particular, using an extended version of the QuasiContinuum method that accounts for both material and geometric nonlinearities. The method allows transition between high resolution discrete representation at areas of high interest to exactly capture the local topology effect on fracture properties, and low resolution computationally homogenized representation at low interest zones to reduce the number of degrees of freedom and computational cost. Dynamic adaptivity in space and time is implemented to guarantee efficient transition between the two resolutions and accurate representation at different levels of stretch. The method enables accurate modeling of crack propagation without any prior constraint on the fracture energy while maintaining the influence of largescale elastic loading in the bulk. In this initial study, we have achieved two orders of magnitude speed up in computation, without parallelization, compared to a full discrete model and we expect the efficiency to further increase as the sample size increases.
While methods for modeling the nonlinear elastic response of polymeric materials are well established [41, 6] , modeling fracture and damage in soft polymeric materials continue to elicit challenges in science and engineering [42, 43, 44, 45, 46, 47, 48, 49, 50]. However, despite enabling useful insights, most of the existing techniques suffer from a number of limitations. For example, most continuum fracture models are insensitive to local topological or microstructural details at the network level. They also require input for critical parameters such as fracture energy which may not be easily measured in many soft materials or may exhibit size effects [51]. For phase field methods, and other nonlocal theories, a length scale is also needed for regularization of the fracture [52]. While in the case of phase field approaches, this length scale is introduced for numerical purposes and the framework should converge to Griffith formulation upon taking this length scale to the zero limit, it is not clear how this length scale is to be chosen in cases when the material has an intrinsic one arising from its microstructure. This compromises the predictive power of these modeling approaches. While fully discrete simulations provide a powerful alternative, as they capture the full details of the material microstructure, the computational cost of these approaches become prohibitive as the sample size increases.
The main advantage of the QC framework is that it enables efficient modeling of large scale problems without neglecting the effect of local topology and without assuming any critical macroscopic parameters such as fracture energy, or material length scale. Only a local failure criterion at the chain level is needed and this may be rigorously computed directly from explicit model of a single chain [32, 3, 9]. The QC method requires much less run time, and memory than a full discrete model while providing higher accuracy and insights into the fracture process than the continuum approaches. Furthermore, the QC has the flexibility of adding small scale physics and topology. This makes it suitable for taking advantage of recent advances in experimental techniques.
While most current experimental techniques do not allow direct visualization of the deformation and damage on the network micro scale, there have been a few exciting advances in recent years towards that goal. For example, the work of Ducrot et al. [53] where a special type of mechanophore molecules is embedded in the polymer network and enables the direct imaging of the spatiotemporal evolution of damage density due to bond breakage. Also the work of Poulain et al. [49] where high speed photography is used to track the onset and growth of cavitation and subsequent coalescence and crack growth. These recent advances in imaging and experimental measurements among others are probing deformation and failure mechanisms in soft materials with increasingly high resolution challenging some longheld ideas about these phenomena [54, 55]. The proposed QC approach is uniquely suited for taking advantage of these advances and making use of the increasingly available data at the microscale. Furthermore, the QC approach provides a rigorous tool for linking the micro and macroscales with minimal set of assumptions which may potentially help resolve several challenges in failure of hierarchical and heterogeneous polymeric materials.
It is worth to note that there exists a large literature on discrete models of fracture that abstract a continuum domain using a lattice representation, such as Lattice Discrete Particle Modeling [56], Graphbased Finite Element Analysis GraFEA [57], and Virtual bond model [58]. A major difference between the QC approach and these methods is that the discrete component of the QC model is physical and represents an actual networked material with a well defined failure criterion at the discrete elements scale. Thus, the QC method does not suffer from some limitations that exist in these other lattice models such as mesh sensitivity, possible ambiguity in defining the fracture criteria upon abstraction of the continuum scale to the lattice scale, or the necessity of introducing nonlocal effects [52] or phenomenological length scales. In particular, the crack path computed from the QC model is not mesh dependent, addressing an important challenge in modeling fracture of inhomogeneous materials [59].
Our results for fracture of disordered networks suggest that as the disorder increases, the expected value of the peak force decreases but the ductility ratio and energy dissipation increase. It is also observed that the crack path becomes nonuniform with increased disorder, as the fracture mechanism includes numerous links failing in the network away from the crack tip zone and later coalescing into a thick crack path (Fig.8b). Meandering crack paths suggest a more ductile response, which is observed in Fig.10 d,e and increases fracture toughness compared to straight cracks. Such rough cracks are observed extensively in experiments [60, 61, 62] and are usually attributed to dynamic instabilities. Here we show that they may also emerge naturally due to the increased level of disorder in the network structure. Transition in the failure mode from brittle to ductile with increased disorder has been reported previously using fully discrete models on smaller samples [63]. The QC method provides a flexible tool for reassessing these phenomena at larger scales and with less sensitivity to sample boundaries.
The current formulation of the QC implementation has a few limitations. For example, we have assumed an affine projection approximation to couple micro and macroscales in the homogenized region of the domain. While we have shown that the homogenization error is within acceptable limits for the cases tested, this affine approximation has some disadvantages. For instance, the accuracy was shown to deteriorate as the disorder increases unless mesh adaptivity is turned on. In addition, we only considered networks with axial interactions. Some polymer chain models have a contribution to their potential energy originating from the rotational degrees of freedom which introduces a nonaffine component in the deformation field [64]. While any forcestretch law may be readily used in the current formulation in place of the wormlike chain model, future work will explore extension of the homogenization algorithm to account for nonaffine deformation.
In this initial study we opted to evaluate the stress and material tangent tensors from the contributions of all links underlying each continuum elements. It is possible, however, to take a representative volume element (RVE) and estimate the stress and tangent tensors from the contribution of links within this volume. The general form of the tensors would be the same, but only links within this RVE will be accounted for. While this homogenization approach leads to increased computational savings, it also has some drawbacks. For example, it is not clear apriori how large the RVE should be in the case of disordered networks. Also the approximation based on an RVE will introduce additional errors. The pros and cons of this alternative methodology will be investigated further in future work.
Viscoelasticity is a critical feature in the mechanical response of polymeric materials such as rubber and gel. In our prior work [32, 3], we have modeled discrete polymer networks with viscoelastic response emerging from two major contributions: (1) breakage and formation of end bonds and crosslinkers, and (2) vsicous drag arising from the relative motion of the deforming polymer network and infused fluids. Extensive work has also been done previously by other groups for experimental characterization of the elastic and loss mdodulii [65, 66]. Continuum models based on phenomenological as well as homogenization approaches also exist [67, 68]. The QC method is capable of incorporating the viscoelastic response. The key point is to map the viscoelastic behavior from the microscale (network level) to the macroscale (QC level) using a proper homogenization rule. As a starting point, one may assume a damping force in each link that is proportional to the link stretch rate in addition to damping force at each network node proportional to its velocity. Then this damping matrix may be homogenized in a similar way as the stiffness component. The first term will result in a damping matrix that has the same structure as the stiffness matrix whereas the latter will result in a diagonal damping matrix. We will then use predictorcorrector schemes to solve the resulting linearized quasidynamic system at each time step. This approach would lead to the emergence of viscoelasticity on the continuum scale.
In this initial work, we considered polymer networks in two dimensions. However, the method may be easily extended to model networks in the 3D domain. Unlike in atomistic systems, where the QC computational cost may exponentially increase from 2D to 3D due to the presence of long range interactions, the extension in lattice systems is straightforward since physical links introduce interaction between directly connected nodes only. Furthermore, extension to 3D will enable incorporation of additional physical mechanisms such as entanglement, which may be modeled at the discrete level as an additional type of connections between the links, such as ring sliding connections [69] .This may then be homogenized using similar methodology to the one presented in this work and passed on to the continuum scale.
To summarize, the QC framework enables significant insight into the role of microstructure on macroscopic fracture response in networked materials. The unique potential of the method lies in its predictive power which qualifies it to play the role of a highly controlled virtual experimental setup. This will not only lead to identification of fundamental processing controlling fracture and deformation in a variety of networks, including polymers, but may also enable the discovery of new material designs with unusual fracture response.
8 Conclusions:
Fracture is a highly nonlinear problem in which macroscopic response depends on the details of microscopic physics. Here we have introduced an adaptive model for analysis and design of networked materials that starts from fundamental link scale mechanics and consistently computes macroscopic response. The method predicts the material response with both material and geometric nonlinearities. We apply this framework to investigate polymer networks where indivdiual polymer chains follow the wormlike chain model. Our main conclusions are summarized as follows:

QC method provides a promising platform for investigating fracture in disordered networked materials that combines advantages from both discrete and continuum approaches.

The QC method for fracture in polymer networks only requires a failure criterion on the chain level. It does not require information about critical energy release rate or a material length scale. The QC discover these quantities as part of the solution. This predictive power and scale bridging capabilities give the QC framework an advantage over other existing methdologies.

Adaptivity allows the QC mesh to evolve as the crack grows, increasing accuracy and computational savings. For the cases presented here, we have achieved two orders of magnitudes in computational savings compared to the full discrete models. Ongoing work with parallel implementation will enable bridging a wider range of scales with further computational savings.

Crack path is controlled by the network topology and stochastic properties of the chains. Thus, QC method illuminates new pathways for design through direct control of the microstructure.

The details of the network local topology have a significant effect on the local force distribution around the crack tip. Disorder makes the crack tip force field more diffused and also channel forces in directions dominated by network anisortropy.

Networks with higher disorder in the geometry or/and chain length have higher stiffness but lower average peak force. Disorder also enhances the softening behavior as disordered networks fails in a more gradual manner compared to uniform networks.

For networks with the same polymer chain concentration, increasing the coordination number leads to stiffer networks with higher peak force, lower stretch at failure, and more gradual softening behavior.

The peak force has a nonmonotonic dependance on the crosslinker concentration where it increases first as the crosslinker concentration increases but then decreases if the concentration exceeds a threshold value. Networks with longer chains have higher threshold value than those with shorter ones.
Acknowledgment
This work is supported by the Center for Geologic Storage of CO2, an Energy Frontier Research Center funded by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES), under Award No. DESC0C12504, and the National Science Foundation CAREER Award (Grant No. 1753249).
References
 Creton and Ciccotti [2016] C. Creton, M. Ciccotti, Fracture and adhesion of soft materials: a review, Reports on Progress in Physics 79 (2016) 046601.
 Everaers and Kremer [2009] R. Everaers, K. Kremer, Topology effects in model polymer networks, AIP Conference Proceedings 615 (2009) 615–626.
 Kothari et al. [2018] K. Kothari, Y. Hu, S. Gupta, A. Elbanna, Mechanical Response of TwoDimensional Polymer Networks: Role of Topology, Rate Dependence, and Damage Accumulation, Journal of Applied Mechanics 85 (2018) 031008.
 Gent [1996] A. N. Gent, A New Constitutive Relation for Rubber, Rubber Chemistry and Technology 69 (1996) 59–61.
 Rivlin [1948] R. S. Rivlin, Large Elastic Deformations of Isotropic Materials. IV. Further Developments of the General Theory, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 241 (1948) 379–397.
 M. Arruda and C. Boyce [1993] E. M. Arruda, M. C. Boyce, A threedimensional constitutive model for the large stretch behavior of rubber elastic materials, Journal of the Mechanics and Physics of Solids 41 (1993) 389–412.
 Hu et al. [2010] Y. Hu, X. Zhao, J. J. Vlassak, Z. Suo, Using indentation to characterize the poroelasticity of gels, Applied Physics Letters 96 (2010).
 Chester and Anand [2011] S. A. Chester, L. Anand, A thermomechanically coupled theory for fluid permeation in elastomeric materials: Application to thermally responsive gels, Journal of the Mechanics and Physics of Solids 59 (2011) 1978–2006.
 Mao and Anand [2018] Y. Mao, L. Anand, Fracture of Elastomeric Materials by Crosslink Failure, Journal of Applied Mechanics 85 (2018) 081008.
 Miehe and Schänzel [2014] C. Miehe, L. M. Schänzel, Phase field modeling of fracture in rubbery polymers. Part I: Finite elasticity coupled with brittle failure, Journal of the Mechanics and Physics of Solids 65 (2014) 93–113.
 Kumar et al. [2018] A. Kumar, G. A. Francfort, O. LopezPamies, Fracture and healing of elastomers: A phasetransition theory and numerical implementation, Journal of the Mechanics and Physics of Solids 112 (2018) 523–551.
 Spatschek et al. [2011] R. Spatschek, E. Brener, A. Karma, Phase field modeling of crack propagation, Philosophical Magazine 91 (2011) 75–95.
 Griffith [1921] A. A. Griffith, The Phenomena of Rupture and Flow in Solids, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 221 (1921) 163–198.
 Zhang et al. [2017] J. Zhang, S. O. Poulsen, J. W. Gibbs, P. W. Voorhees, H. F. Poulsen, Determining material parameters using phasefield simulations and experiments, Acta Materialia 129 (2017) 229–238.
 Buehler [2008] M. J. Buehler, Atomistic Modeling of Mechanical Failure, 2008.
 Deutsch and Binder [1991] H. P. Deutsch, K. Binder, Interdiffusion and selfdiffusion in polymer mixtures: A Monte Carlo study, The Journal of Chemical Physics 94 (1991) 2294–2304.
 Carmesin and Kremer [1988] C. Carmesin, K. Kremer, The Bond Fluctuation Method: A New Effective Algorithm for the Dynamics of Polymers in All Spatial Dimensions, Macromolecules 21 (1988) 2819–2823.
 Tadmor et al. [1996] E. B. Tadmor, M. Ortiz, R. Phillips, Quasicontinuum analysis of defects in solids, Philosophical Magazine A: Physics of Condensed Matter, Structure, Defects and Mechanical Properties 73 (1996) 1529–1563.
 Miller et al. [2002] R. E. Miller, E. B. Tadmor, R. E. Miller, The Quasicontinuum Method: Overview, applications and current directions, Journal of ComputerAided Materials Design 9 (2002) 203–239.
 Shenoy et al. [1999] V. B. Shenoy, R. Miller, E. B. Tadmor, D. Rodney, R. Phillips, M. Ortiz, An adaptive finite element approach to atomicscale mechanics  The quasicontinuum method, Journal of the Mechanics and Physics of Solids 47 (1999) 611–642.
 Tadmor et al. [1999] E. B. Tadmor, R. Miller, R. Phillips, M. Ortiz, Nanoindentation and incipient plasticity, Journal of Materials Research 14 (1999) 2233–2250.
 Knap and Ortiz [2001] J. Knap, M. Ortiz, An analysis of the quasicontinuum method, Journal of the Mechanics and Physics of Solids 49 (2001) 1899–1923.
 Miller et al. [1998] R. Miller, E. B. Tadmor, R. Phillips, M. Ortiz, Quasicontinuum simulation of fracture at the atomic scale, Modelling and Simulation in Materials Science and Engineering 6 (1998) 607–638.
 Mikeš and Jirásek [2017] K. Mikeš, M. Jirásek, Quasicontinuum method extended to irregular lattices, Computers and Structures 192 (2017) 50–70.
 Beex et al. [2014] L. A. Beex, P. Kerfriden, T. Rabczuk, S. P. Bordas, Quasicontinuumbased multiscale approaches for platelike beam lattices experiencing inplane and outofplane deformation, Computer Methods in Applied Mechanics and Engineering 279 (2014) 348–378.
 Phlipot and Kochmann [2018] G. P. Phlipot, D. M. Kochmann, A quasicontinuum theory for the nonlinear mechanical response of general periodic truss lattices, Journal of the Mechanics and Physics of Solids 124 (2018) 758–780.
 Hill [1963] R. Hill, Elastic properties of reinforced solids: Some theoretical principles, Journal of the Mechanics and Physics of Solids 11 (1963) 357–372.
 Suquet [1985] P. M. Suquet, Local and global aspects in the mathematical theory of plasticity., In Plasticity Today: Modelling, Methods and Applications (1985) 279–310.
 Geers et al. [2017] M. G. D. Geers, V. G. Kouznetsova, K. Matouš, J. Yvonnet, Homogenization Methods and Multiscale Modeling: Nonlinear Problems, in: Encyclopedia of Computational Mechanics Second Edition, John Wiley & Sons, Ltd, Chichester, UK, 2017, pp. 1–34.
 Marko and Siggia [1995] J. F. Marko, E. D. Siggia, Stretching DNA, Macromolecules 28 (1995) 8759–8770.
 Marko [1998] J. F. Marko, DNA under high tension: Overstretching, undertwisting, and relaxation dynamics, Physical Review E  Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 57 (1998) 2134–2149.
 Lieou et al. [2013] C. K. Lieou, A. E. Elbanna, J. M. Carlson, Sacrificial bonds and hidden length in biomaterials: A kinetic constitutive description of strength and toughness in bone, Physical Review E  Statistical, Nonlinear, and Soft Matter Physics 88 (2013) 1–10.
 Rokoš et al. [2017] O. Rokoš, R. H. Peerlings, J. Zeman, eXtended variational quasicontinuum methodology for lattice networks with damage and crack propagation, Computer Methods in Applied Mechanics and Engineering 320 (2017) 769–792.
 Rivara [1997] M. C. Rivara, New longestedge algorithms for the refinement and/or improvement of unstructured triangulations, International Journal for Numerical Methods in Engineering 40 (1997) 3313–3324.
 Anderson [2017] T. L. Anderson, Fracture Mechanics: Fundamentals and Applications, Fourth Edition, CRC Press, 2017.
 Lee and Pharr [2019] S. Lee, M. Pharr, Sideways and stable crack propagation in a silicone elastomer, Proceedings of the National Academy of Sciences (2019) 201820424.
 Bonn et al. [2014] D. Bonn, H. Kellay, K. Bendjemiaa, J. Meunier, Delayed Fracture of an Inhomogeneous Soft Solid, Science 280 (2014) 265–267.
 Ghelichi and Kamrin [2015] R. Ghelichi, K. Kamrin, Modeling growth paths of interacting crack pairs in elastic media, Soft Matter 11 (2015) 7995–8012.
 Fender et al. [2010] M. L. Fender, F. Lechenault, K. E. Daniels, Universal shapes formed by two interacting cracks, Physical Review Letters 105 (2010) 1–4.
 Wang and Gao [2016] Q. Wang, Z. Gao, A constitutive model of nanocomposite hydrogels with nanoparticle crosslinkers, Journal of the Mechanics and Physics of Solids 94 (2016) 127–147.
 Boyce and Arruda [2000] M. C. Boyce, E. M. Arruda, Constitutive Models of Rubber Elasticity: A Review, Rubber Chemistry and Technology 73 (2000) 504–523.
 Qi et al. [2018] Y. Qi, J. Caillard, R. Long, Fracture toughness of soft materials with rateindependent hysteresis, Journal of the Mechanics and Physics of Solids 118 (2018) 341–364.
 Zhao [2012] X. Zhao, A theory for large deformation and damage of interpenetrating polymer networks, Journal of the Mechanics and Physics of Solids 60 (2012) 319–332.
 Nguyen et al. [2005] T. D. Nguyen, S. Govindjee, P. A. Klein, H. Gao, A material force method for inelastic fracture mechanics, Journal of the Mechanics and Physics of Solids 53 (2005) 91–121.
 Guo et al. [2018] J. Guo, M. Liu, A. T. Zehnder, J. Zhao, T. Narita, C. Creton, C. Y. Hui, Fracture mechanics of a selfhealing hydrogel with covalent and physical crosslinks: A numerical study, Journal of the Mechanics and Physics of Solids 120 (2018) 79–95.
 Yu et al. [2018] K. Yu, A. Xin, Q. Wang, Mechanics of selfhealing polymer networks crosslinked by dynamic bonds, Journal of the Mechanics and Physics of Solids 121 (2018) 409–431.
 Talamini et al. [2018] B. Talamini, Y. Mao, L. Anand, Progressive damage and rupture in polymers, Journal of the Mechanics and Physics of Solids 111 (2018) 434–457.
 Mao and Anand [2018] Y. Mao, L. Anand, A theory for fracture of polymeric gels, Journal of the Mechanics and Physics of Solids 115 (2018) 30–53.
 Poulain et al. [2017] X. Poulain, V. Lefèvre, O. LopezPamies, K. RaviChandar, Damage in elastomers: nucleation and growth of cavities, microcracks, and macrocracks, International Journal of Fracture 205 (2017).
 Lavoie et al. [2019] S. R. Lavoie, P. Millereau, C. Creton, R. Long, T. Tang, A continuum model for progressive damage in tough multinetwork elastomers, Journal of the Mechanics and Physics of Solids 125 (2019) 523–549.
 Zhang et al. [2019] B. Zhang, C. S. Shiang, S. J. Yang, S. B. Hutchens, YShaped Cutting for the Systematic Characterization of Cutting and Tearing, Experimental Mechanics (2019).
 PijaudierCabot and Bazant [1988] G. PijaudierCabot, Z. P. Bazant, Nonlocal damage theory, Journal of Engineering Mechanics 113 (1988) 1512–1533.
 Ducrot et al. [2014] E. Ducrot, Y. Chen, M. Bulters, R. P. Sijbesma, C. Creton, Toughening elastomers with sacrificial bonds and watching them break, Science 344 (2014) 186–189.
 Lefèvre et al. [2015] V. Lefèvre, K. RaviChandar, O. LopezPamies, Cavitation in rubber: an elastic instability or a fracture phenomenon?, International Journal of Fracture 192 (2015) 1–23.
 Millereau [2019] P. C. C. C. Y. Millereau, Detection of Molecular Fracture in Elastomers with Mechanophores, in: Bulletin of the American Physical Society, p. 1.
 Cusatis et al. [2011] G. Cusatis, D. Pelessone, A. Mencarelli, Lattice Discrete Particle Model (LDPM) for failure behavior of concrete. I: Theory, Cement and Concrete Composites 33 (2011) 881–890.
 Khodabakhshi et al. [2016] P. Khodabakhshi, J. N. Reddy, A. Srinivasa, GraFEA: a graphbased finite element approach for the study of damage and fracture in brittle materials, Meccanica 51 (2016) 3129–3147.
 Zhang [2013] Z. Zhang, Discretized virtual internal bond model for nonlinear elasticity, International Journal of Solids and Structures 50 (2013) 3618–3625.
 Bazant [2010] Z. P. Bazant, Can MultiscaleMultiphysics Methods Predict Softening Damage and Structural Failure?, International Journal for Multiscale Computational Engineering 8 (2010) 61–67.
 RaviChandar and Knauss [1984] K. RaviChandar, W. G. Knauss, An experimental investigation into dynamic fracture: III. On steadystate crack propagation and crack branching, International Journal of Fracture 26 (1984) 141–154.
 Carlson [2018] E. K. Carlson, Watch tiny cracks travel in 3D, EOS 99 (2018).
 Steinhardt and Rubinstein [2015] W. Steinhardt, S. Rubinstein, High Speed Strain Measurements Surrounding Hydraulic Fracture in Brittle Hydrogel, in: APS Division of Fluid Dynamics Meeting Abstracts, p. A1.004.
 Driscoll et al. [2016] M. M. Driscoll, B. G.g. Chen, T. H. Beuman, S. Ulrich, S. R. Nagel, V. Vitelli, The role of rigidity in controlling material failure, Proceedings of the National Academy of Sciences 113 (2016) 10813–10817.
 Broedersz and Mackintosh [2014] C. P. Broedersz, F. C. Mackintosh, Modeling semiflexible polymer networks, Reviews of Modern Physics 86 (2014) 995–1036.
 Li et al. [2012] J. Li, Y. Hu, J. J. Vlassak, Z. Suo, Experimental determination of equations of state for ideal elastomeric gels, Soft Matter 8 (2012) 8121–8128.
 Nguyen et al. [2010] T. D. Nguyen, C. M. Yakacki, P. D. Brahmbhatt, M. L. Chambers, Modeling the relaxation mechanisms of amorphous shape memory polymers, Advanced Materials 22 (2010) 3411–3423.
 Kumar and LopezPamies [2016] A. Kumar, O. LopezPamies, On the twopotential constitutive modeling of rubber viscoelastic materials, Comptes Rendus  Mecanique 344 (2016) 102–112.
 Zhang and OstojaStarzewski [2015] J. Zhang, M. OstojaStarzewski, Mesoscale bounds in viscoelasticity of random composites, Mechanics Research Communications 68 (2015) 98–104.
 Noda et al. [2014] Y. Noda, Y. Hayashi, K. Ito, From topological gels to slidering materials, Journal of Applied Polymer Science 131 (2014) 1–9.