An Adaptive Quasi-Continuum Approach for Modeling Fracture in Networked Materials: Application to Modeling of Polymer Networks

An A daptiv e Quasi-Con tin uum Approac h for Mo deling F racture in Net w ork ed Materials: Application to Mo deling of P olymer Net w orks Ahmed Ghareeb 1 and Ahmed Elbanna ∗ 1 1 Departmen t of Civil and Environmen tal Engineering, Univ ersit y of Illinois at Urbana-Champaign, Illinois, USA Ma y 2019 Abstract Materials with net work-lik e microstructure, including p olymers, are the bac kb one for man y natural and h uman-made materials such as gels, biological tissues, metamaterials, and rubb ers. F racture pro cesses in these net w orked materials are in trinsically multiscale, and it is compu- tationally prohibitiv e to adopt a fully discrete approac h for large scale systems. T o ov ercome suc h a challenge, we introduce an adaptiv e numerical algorithm for modeling fracture in this class of materials, with a primary application to p olymer netw orks, using an extended v ersion of the Quasi-Contin uum metho d that accounts for b oth material and geometric nonlineari- ties. In regions of high interest, for example near crack tips, explicit representation of the lo cal top ology is retained where eac h p olymer c hain is idealized using the worm like chain mo del. A w ay from these imperfections, the degrees of freedom are limited to a fraction of the net work no des and the net work structure is computationally homogenized, using the micro- macro energy consistency condition, to yield an anisotropic material tensor consistent with the underlying net work structure. A nonlinear ﬁnite elemen t framework including b oth mate- rial and geometric nonlinearities is used to solv e the system where dynamic adaptivit y allo ws transition b et ween the contin uum and discrete scales. The metho d enables accurate mo delling of crack propagation without a priori constraint on the fracture energy while main taining the inﬂuence of large-scale elastic loading in the bulk. W e demonstrate the accuracy and eﬃciency of the method b y applying it to study the fracture in diﬀerent examples of netw ork structures. W e further use the method to inv estigate the eﬀects of net work topology and disorder on its fracture c haracteristics. W e discuss the implications of our metho d for m ultiscale analysis of fracture in net work ed material as they arise in diﬀerent applications in biology and engineering. Keywor ds: Quasi-Con tinuum Metho d, Polymer Netw orks, F racture 1 In tro duction P olymers are the building blo cks for many natural and artiﬁcial materials. The load transfer system in p olymeric materials ma y b e abstracted as a complex netw ork of cross-link ed nonlinear c hains [ 1 ]. The top ological prop erties of p olymer netw orks, suc h as lo cal connectivity , cross-linking density , b ond types, and sacriﬁcial b onds and hidden length, ma y dramatically aﬀect their mec hanical resp onse and fracture p roprieties [ 2 , 3 ]. Most curren t exp erimen tal techniques do not allo w direct visualization of the deformation and damage on the net work micro scale, thus, the relation b et ween the topology and mechanical b eha vior of these materials remains elusiv e. Multiscale n umerical mo dels are th us needed to bridge this gap [ 2 ]. Mo deling of soft p olymeric materials has b een mainly approached using contin uum theories, includ- ing linear elasticit y , h yp erelasticit y , visco elasticit y , viscoplasticity and p oro elasticit y [ 4 , 5 , 6 , 7 , 8 ]. ∗ Corresponding author: Ahmed Elbanna, elbanna2@illinois.edu 1 Con tinuum represen tation of materials has b een used in traditional solid mec hanics, where the material is assumed to b e a contin uous medium with its resp onse describ ed mathematically by a set of constitutive laws. Through exp erimen ts, the parameters of the constitutiv e la ws ma y b e determined and then utilized to solve the b oundary v alue problems using either analytical or com- putational approac hes. These tec hniques enabled signiﬁcan t progress in modeling the elasticity and distributed damage b eha vior of p olymeric materials with and without infused ﬂuids. Ho wev er, these contin uum metho ds are usually indep enden t of the size scale of the netw ork structures and they ma y not capture lo calized phenomena without enric hment, whic h in many cases is either ex- p ensiv e or is based on idealized models of microstructure. In particular, when it comes to fracture, a phenomenon that dep ends critically on the lo cal conditions in the vicinity of the propagating crac ks, contin uum metho ds are not suited for predicting the eﬀect of the netw ork lo cal top ology on its fracture resistance. F urthermore, con tin uum metho ds may also fall short in guiding the net work geometric design for enhanced fracture toughness. An exception in the recen t y ears has b een the emergence of phase ﬁeld models whic h is increasingly used to simulate fracture in a wide class of materials including gels and rubb ers [ 9 , 10 , 11 ]. The metho d may b e used to simulate curved crac k paths, crac k kinking and branc hing, and crac k fron t segmen tation. Phase ﬁled metho ds are based on gradient enhanced damage mo dels where a parameter is in tro duced to distinguish betw een fully intact material and fully damaged material in a gradual manner [ 12 ]. The width of the damage zone dep ends on a length scale that is introduced primarily for regularization purposes but has b een shown theoretically that as this length scale go es to zero, Griﬃth fracture criterion [ 13 ] is recov ered. The damage v ariable changes from 1 (fully intact) to 0 (F ully damaged) ov er this length scale. This length scale must b e carefully c hosen to regularize the strain-softening b ehavior during the fracture pro cess, and to a void mesh dep endency-related issues during ﬁnite elemen t simulations. The choice of this lengths scale for materials with microstructure, whic h p ossess physical length scales of their own, is not clear apriori and is left, in most cases, as a tunable parameter to ﬁt exp erimen tal observ ations [ 14 ]. Recen tly , thanks to adv ances in imaging, computational p o wer, and nano/micro fabrication tech- nologies, the mec hanics of materials tak es on another challenge of mo deling fundamen tal material b eha vior starting from the atomic scale [ 15 ]. F or the accurate represen tation of the material re- sp onse at this scale, the material m ust b e mo deled using discrete simulations. Discrete simulations ha ve been widely used to mo del polymeric materials as net works of irregular lattices with diﬀeren t link models that incorporate the polymer chain physics [ 16 ]. Adapting these discrete models where eac h p olymer c hain is mo deled explicitly , it is p ossible to accurately describ e lo calized phenomena suc h as fracture and ca vitation. F urthermore, a new synthesis paradigm has b een introduced in recen t years for realizing near ideal p olymer netw orks. These ideal cov alent net works, which are pro duced b y the cross-coupling of macro-molecular building blo c ks [ 17 , 18 ], are essentially free of en tanglements and hav e few top ological defects which leads to high strength [ 17 ]. Ideal netw orks pro vide a pathw a y for exploring an explicit connection b et ween net work design and its mec hanical resp onse [ 19 , 20 ], as w ell as a promising platform for designing of new lattice-like materials with tunable prop erties [ 21 ]. Motiv ated by these ideal netw orks, Alame and Brassart [ 20 ] studied the elasticit y of discrete near-ideal polymer netw orks, using a force extension model for polymers con- ﬁned to a surface, with diﬀerent co ordination num b ers. How ever, the study w as limited to small scale samples for computational considerations. In this pap er, W e introduce a new adaptive numerical algorithm for mo deling fracture in p olymer net works using an extended version of the Quasi-Contin uum method [ 22 , 23 , 24 ] that accoun ts for b oth the nonlinear elastic nature of the p olymer chains as well as the geometric nonlinearity asso ciated with their p oten tially large deformation. In regions of high interest where deformation gradien ts are non uniform, for example near crack tips, explicit representation of the lo cal top ology is retained where each p olymer chain is idealized using the worm like c hain model. A wa y from these imp erfections where the deformation gradient is suﬃcien tly uniform, the netw ork structure is computationally homogenized, using an equiv alence of the microscale and macroscale incremen tal w ork, to yield an anisotropic material tensor consistent with the underlying netw ork structure. Th us at any instant in time, only a fraction of the netw ork no des is solved. Dynamic adaptivity allo ws eﬃcient transition b et ween the t wo resolutions. Our goal is to develop a metho d that enables mo deling crack propagation without apriori constraint on the fracture energy or the need to assume phenomenological length scales, while maintaining the inﬂuence of large-scale elastic loading in the bulk. 2 The remainder of the paper is organized as follows. In Section 2 we presen t an o verview of the QC metho d. W e introduce the problem form ulation, the homogenization tec hnique, and the n umerical implemen tation algorithms used in this study in Sections 3 and 4. W e verify the numerical algo- rithm for mo deling b oth pristine and crac ked samples in Section 5. In Section 6 , w e demonstrate the eﬃciency of the metho d b y applying it to study the fracture of in mo del systems that are prohibitiv e for a full-ﬁeld discrete approach. W e further use the metho d to study the eﬀects of net work top ology on its fracture resistance. Finally , W e discuss the method implications for the analysis of netw ork ed material systems in Section 7 . 2 Quasi-Con tin uum metho d o v erview: T admor, Ortiz, and Phillips [ 22 ] prop osed the Quasi-contin uum (QC) metho d in 1996 to ov ercome the scale limitations of atomic sim ulations. The original QC metho d is a computational tec hnique to mo del an atomistic system without the need to treat all atoms in the domain explicitly . The QC pro vides a framework where degrees of freedom are judiciously eliminated and force or energy calculations are approximated at areas of low in terest, whereas exact discrete representation is used at the area of high in terest such as crac k tips or dislocation zones. The QC metho d also pro vides an adaptive framework for the fully resolved critical part to evolv e during the simulation to balance b et ween computational cost and error estimates. The QC metho d ﬁnds the deformation ﬁeld that minimizes the system potential energy while ac hieving the follo wing three criteria [ 23 ]: 1. The num b er of degrees of freedom is reduced by limiting the degrees of freedoms to a small fraction of the no des, called representativ e nodes (repno des), but the full discrete description is retained in high interest regions, around crac k tip for example. Hence, the computational cost of solving the system is signiﬁcantly reduced. 2. 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 signiﬁcantly reduced. 3. The fully resolv ed critical regions can evolv e with the deformation, during the simulation. Th us, the metho d ac hieves eﬃciency b y using adaptivit y to reduce the fully resolved domain as muc h as p ossible. Figure 1: A schematic for domain discretization for a notched sample with a lattice micro-structure using the Quasi-Contin uum metho d: (a) the underlying netw ork, and (b) the QC mesh with full net work resolution around the crack, and homogenized material aw ay from it. The degrees of freedom at the low in terest area are limited to the repno des. The QC method applied to atomistic systems has been successfully used to in vestigate man y lo calized phenomena, such as Nano-indentation, crac k-tip deformation, deformation and fracture of grain b oundaries, and dislo cation interactions [ 25 , 26 , 27 ]. 3 While the QC method was originally introduced for atomic simulations, it ma y b e extended to any discrete systems including net works or lattices. Beex et al. extended the QC metho d to regular discrete lattice netw orks with short-range nearest-neigh b or interactions for b oth conserv ativ e [ 28 ], and non-conserv ative lattice mo dels such as dissipation, b ond failure, and ﬁb er sliding [ 29 , 30 ], and used the metho d to model failure in electronic textiles [ 31 ]. The metho d was also extended to uniform linear elastic b eam-lik e lattices exp eriencing planer b ending [ 32 ] and b oth in plane and out of plane deformation [ 33 ]. Rok os et al. used energy-based v ariational QC framew ork for mo deling regular net works with linear lattices undergoing localized damage [ 34 , 35 ]. Mik es and Jirasek extended the quasi con tinuum metho d to irregular linear elastic lattices with axial interactions and small deformations [ 36 ]. They applied the metho d to sim ulate nanotextile structures [ 37 ], and mo del the microstructure based on the anisotropic microplane model [ 38 ]. Phlip ot et al. [ 39 ] used the metho d to model linear elastic p eriodic b eam lattices with co-rotational framework. In addition to lattice structures, the QC metho dology has b een extended for simulation of gran ular media where the discrete elemen t method is used at areas of high in terest while homogenization is utilized elsewhere for computational eﬃciency [ 40 , 41 ]. The ﬂexibilit y of the framework mak es it a suitable candidate for applications to irregular p olymer netw orks with both material and geometric nonlinearities, which is the focus of this study . 3 Quasi-Con tin uum for irregular p olymer net w orks: In this section, we discuss the problem setup and the ﬁnite element formulation for b oth the discrete part and the homogenized part of the netw ork mo del. F or the discrete part, w e use a nonlinear 1D ﬁnite element with ﬁnite deformation and rotation. F or the homogenized part, we use 2D nonlinear triangular ﬁnite element with ﬁnite deformations. The material properties for the 2D elements are deriv ed by applying homogenization rule based on the micro-macro energy consistency condition that conserves the v ariation of lo cal work b et ween the macro (2D homogenized elements) and micro (underlying netw ork) scales. 3.1 Finite Element F ormulation: 2D homogenized elemen ts The strong form for the b oundary v alue problem in the absence of inertia eﬀects is giv en by: ∂ P ik ∂ X k + ρ 0 b i = 0 (1) Where P is the 1 st Piola-Kirc hhoﬀ stress tensor, X is the material p oin t coordinate in the reference conﬁguration, b is the b o dy force, and ρ 0 is the material density at the reference conﬁguration. The b oundary conditions are giv en by: P ij N j = ¯ t i on Γ o t i and u i = ¯ u i on Γ o u i (2) Where Γ o t i ∩ Γ o u i = φ , Γ o t i ∪ Γ o u i = Γ o Where ¯ t i is the applied traction on b oundary Γ o t i , N is the unit v ector normal to the b oundary surface in the reference conﬁguration, ¯ u i is the prescrib ed displacemen t on b oundary Γ o u i , and Γ o is the b oundary surface. The weak form in a total Lagrangian formulation may b e written as: Z Ω o δ F ik P ik d Ω o = Z Ω o δ u i ρ 0 b i d Ω o + Z Γ o t i δ u i ¯ t i Γ o (3) Where F is the deformation gradient, and δ u ∈ u o where u o the space of kinematically admissible displacemen ts with the requirement that the displacemen ts v anish on displacement boundaries. The ﬁnite elemen t appro ximation is deriv ed in terms of the ob jective 2 nd Piola-Kirc hhoﬀ stress 4 tensor S . After lineraization and making use of V oigt notation, the in ternal force vector and material tangent matrix are giv en by: f int I = Z Ω o B T oI S d Ω o and K I J = K mat I J + K g eo I J (4) Where K mat I J = Z Ω o B T oI [ C S E ] B oJ d Ω o and K g eo I J = I Z Ω o B T oI SB oJ d Ω o (5) Where B o is the strain-displacemen t matrix. The goal of the homogenization is to deriv e the macroscopic stress S and material tangen t C S E tensors that map the micro-scale behavior with suﬃcien t accuracy . These quantities are then pro vided to the nonlinear ﬁnite element framework for solving the system. 3.2 Homogenization: W e replace the netw ork links, a wa y from the high interest zones, b y contin uum elements with material prop erties obtained by homogenization of the discrete netw ork. This results in a signiﬁ- can t n umber of netw ork links b eing remo ved from the assem bly pro cedure, leading to signiﬁcant computational sa vings. This approac h was pioneered b y Mik es and Jirasek [ 36 ] to mo del irregular linear elastic lattices with axial interactions and small deformations. In this study , we extend the form ulation to nonlinear links with ﬁnite deformations. 3.2.1 Micro-macro energy consistency condition: W e establish the micro-macro scale transition relation based on the Hill-Mandel condition [ 42 , 43 ]. This condition requires the volume av erage of the v ariation of work on the micro lev el to b e equal to the v ariation of lo cal w ork on the macroscale: δ W M acro = δ W micro (6) In terms of macro deformation gradien t tensor and the ﬁrst Piola-Kirc hhoﬀ stress tensor, the condition reads [ 44 ]: P M : δ F M = 1 V o N p X p =1  ~ f p .δ ~ l p  (7) Where ~ f p and ~ l p are the force and the length of eac h link in the curren t conﬁguration, N p is the n umber of links. Using δ ~ l p = ~ L p δ F T , where ~ L p is the link length in the original conﬁguration, Eq. 7 b ecomes: P M : δ F M = 1 V o N p X p =1  ~ f p ⊗ ~ L p  : δ F m (8) Where the subscript M denotes the macro level tensors, and m denotes the micro level tensors. Using V oigt assumption (i.e F m = F M ). The homogenized 1 st Piola-Kirc hhoﬀ stress is giv en by: P M = 1 V o N p X p =1  ~ f p ⊗ ~ L p  (9) T o calculate the tangen t tensor, w e start by taking the v ariation of P M : δ P M = 1 V o N p X p =1 ∂ ~ f p ∂ ~ l p .δ ~ l p ! ⊗ ~ L p (10) Substituting δ ~ l p = δ F ~ L p and rearranging: 5 δ P M = 1 V o N p X p =1 ~ L p ⊗ ∂ ~ f p ∂ ~ l p ⊗ ~ L p ! LT : δ F m (11) Where the subscript L T denotes the transpose of the left tw o indices. Again, using V oigt assump- tion, it follows that: C P F = 1 V o N p X p =1 ~ L p ⊗ ∂ ~ f p ∂ ~ l p ⊗ ~ L p ! LT (12) Where C P F is the homogenized 4 th order ﬁrst elasticit y tensor. It must be noted that the homog- enized stress tensor is a tw o-p oin t tensor (i.e not symmetric), and the elasticity tensor has only ma jor symmetry . It is more con venien t to use the ob jectiv e symmetric 2 nd Piola-Kirc hhoﬀ Stress tensor and the corresp onding elasticity tensor. By deﬁnition, the 2 nd Piola-Kirc hhoﬀ Stress tensor is given by: S M = F − 1 M P M (13) Using indicial notation and substituting for P M from Eq. 9 S M ij = F M ik − 1 1 V o N p X p =1  f p k L p j  (14) After using V oigt assumption: S M ij = 1 V o N p X p =1  F p ik − 1 f p k L p j  (15) W e note that the quan tity F p ik − 1 f p k is the (ﬁctitious) force vector corresp onding to the PK2 stress. T o derive the elasticity tensor relating the PK2 tensor and the Green strain, we start by v arying the relation b et ween PK1 and PK2: δ S M ij = F M im − 1 δ P M mj | {z } T erm 1 + δ F M iu − 1 P M uj | {z } T erm 2 (16) After substituting from Eq. 12 and using V oigt assumption, the ﬁrst term may then b e written as: (w e drop the subscript p for conv enience) F M im − 1 δ P M mj = 1 V o N p X p =1  F − 1 im L j ∂ f m ∂ l n L l  : δ F nl (17) Simplifying and using δ E kl = F nk δ F nl , δ E kl = δ E lk : F M im − 1 δ P M mj = 1 V o N p X p =1  L j F − 1 im ∂ f m ∂ l n F − 1 kn L l  : δ E kl (18) The second term is written as: δ F M iu − 1 P M uj = ∂ F − 1 iu ∂ F nk δ F nk P uj (19) After simpliﬁcation and some indices manipulation: δ F M iu − 1 P M uj = − F − 1 in F − 1 ku P uj : δ F nk (20) δ F M iu − 1 P M uj = − F − 1 in S kj : δ F nk (21) 6 δ F M iu − 1 P M uj = − F − 1 in F − 1 ln S kj : δ E lk (22) Hence: C S E ij kl = 1 V o N p X p =1  L j F − 1 im ∂ f m ∂ l n F − 1 kn L l  − F − 1 in F − 1 ln S j k (23) Where C S E is the homogenized 4 th order second elasticity tensor. 3.3 Assem bling the tensors: Since the links in the netw ork hav e axial deformations only , the deformation gradient for eac h link is written as: F ij = λn i N j (24) Where λ = l /L is the stretch, n i is the unit v ector in the current conﬁguration, and N i is the unit v ector in the original conﬁguration. Similarly: F − 1 ij = 1 λ N i n j (25) The PK1 stress tensor is then written as P ij = 1 V o N p X p =1 ( f L ) n i N j (26) And the PK2 stress tensor: S ij = 1 V o N p X p =1  1 λ f L  N i N j (27) Whereas the material tensors are given b y C P F ij kl = 1 V o N p X p =1  k L 2  n i N j n k N l (28) C S E ij kl = 1 V o N p X p =1  1 λ 2 k L 2 − 1 λ 3 f L  N i N j N k N l (29) Here k = ∂ f /∂ l for the link, f , and L are the magnitudes of the link force and reference length resp ectiv ely . The stress tensor S is symmetric, and the elasticit y tensor C S E has b oth ma jor and minor symmetries. 3.4 P olymer chain constitutiv e law: In this study , we mo del eac h link in the p olymer net work using a nonlinear elastic force elongation relation give by the w orm like c hain mo del. The force-elongation relation is given b y [ 45 ]: f = k B T b " 1 4  1 − x L c  − 2 − 1 4 + x L c # (30) Where f is the force, x is the end-to-end distance, b is the p ersistence length, k B is Boltzmann’ constan t, T is the temp erature, and L c is the contour length of the polymer c hain. W e assume the net work is force balanced at the reference conﬁguration. The chain stiﬀness is given by the slop e of the tangent of the force elongation curve and it tak es the form k = d f /dx = k B T b " 1 2 L c  1 − x L c  − 3 + 1 L c # (31) 7 These expressions are used in Eq. 27 and Eq. 29 to compute the homogenized stress and elasticity tensor. As the end-to-end distance x approac hes the con tour length L c , the p olymer chain resp onse b ecomes highly nonlinear as b oth the force and the stiﬀness v alues go to inﬁnity . This signals a limitation of this constitutiv e description which may b e circum ven ted b y accoun ting for strain energy of the c hain in addition to its en tropic energy [ 46 ]. Suc h correction will b e in v estigated further in future w ork. In this initial study , such extreme limit is not probed as w e adopt a maxim um stretch failure criterion for the links. That is, we assume the link w ould fail if its stretc h relativ e to the av ailable con tour length λ p = l/L c exceeds a threshold v alue that is predeﬁned. Other failure criterion based on transition state theory for b ond break age [ 47 ] may also b e used. 4 Numerical Implemen tation: T o couple the micro- and macro-scales, the net work nodal displacements u n are related to the represen tative no des (repno des) displacemen t u R , see Fig. 3 -c, through the interpolation functions: u n = N rep X i =1 φ R u R (32) where φ R is the interpolation function associated with repnode R and N rep is the num b er of repno des. Using the compact supp ort of the ﬁnite element shap e functions, the displacement of eac h no de is determined from the sum o ver the three vertices of the triangle con taining this no de. In this initial study , w e use linear shap e functions whic h turns out to p erform satisfactorily as will b e discussed shortly . After ev aluating the netw ork no dal displacemen ts, the stress and material tangen t tensors are assem bled as follo ws: Algorithm 1: Homogenization of QC elements Calculate the netw ork no dal displacemen ts from the last QC solution through in terp olation functions given in Eq. 32 ; for e ach QC element Ω i do for e ach link inside element Ω i do calculate link force f from Eq. 30 ; calculate link stiﬀness k from Eq. 31 ; calculate link stretch λ = l /L ; up date S using Eq. 27 ; up date C S E using Eq. 30 ; end end Since we use the total Lagrangian formulation, the reference conﬁguration quantities, suc h as the reference length L , and the reference unit v ector N , are calculated at the system initialization and are stored instead of calculating them each time step. This leads to further computational sa vings. As discussed in Section 2 , QC metho dologies normally use tw o reduction steps: limiting the degrees of freedom to a small fraction of the no des and sampling of the lattice interactions for eﬃcien t energy summation [ 29 , 35 ]. How ever, in the current work, due to the net work irregularit y , w e accoun t for the con tribution of all links in eac h con tin uum elemen t in the homogenization pro cedure [ 36 ]. W e discuss p ossible alternative approaches in Section 7 . 4.1 Automatic mesh adaptivit y: Automatic Mesh Adaptivit y is a critical ingredient for the eﬃcient implementation of the QC metho d as it enables the QC mesh to ev olve dynamically based on the crac k propagation path. In our proposed approach, w e mark a 2D ﬁnite element Ω i for reﬁnemen t if it meets 1 of the following criteria: 1. Link stretch criterion: λ p > λ th p ∈ Ω i (33) 8 Where λ p = l /L c is the stretch relativ e to the p olymer chain con tour length of any net work link inside the QC element Ω i , and λ th is a threshold v alue. The threshold v alue is chosen as a fraction of the failure stretch. 2. Deformation gradient error criterion [ 23 ]: ε i =  1 V Ω i Z Ω i ( ¯ F − F i ) : ( ¯ F − F i ) dV  1 2 > ε th (34) Where ε i is a scalar measure that quan tiﬁes the error introduced in to the solution by the curren t density of representativ e no des, V Ω i is the volume of element Ω i , F i is the QC solution for the deformation gradient in element Ω i , ¯ F is the L 2 -pro jection of the QC solution for F giv en b y ¯ F = φ F av g where φ is the shap e function arra y , and F av g is computed b y av eraging the deformation gradien t in each elemen t sharing a giv en repnode, and ε th is a speciﬁed threshold v alue. W e start by constructing right-triangulated initial mesh to reduce summation errors as p er [ 48 ]. Eac h QC mesh p oin t is then mov ed to the nearest netw ork no de. All elements marked for reﬁnement are added to a set E ref . F or eac h element inside E ref . W e use the standard Rivera algorithm [ 49 ], whic h conserves non-degeneracy , conformit y , and smo othness. This algorithm tracks the longest- edge propagation path (LEPP) asso ciated with a triangular element Ω i in a backw ard manner. The LEPP of Ω i in a conforming triangular mesh is deﬁned as an ordered list of triangles Ω i , Ω i +1 ,..., Ω n suc h that Ω i +1 is a neighbour b y the longest e dge of Ω i for each i = 1 , 2 , ..., n , the reﬁnemen t scheme is shown in Fig. 2 . The algorithm, whic h is summarized in Algorithm 2, is then used for each Ω i ∈ E ref and rep eated till Ω i is fully resolved. Figure 2: Mesh reﬁnemen t steps: (a) the initial mesh where Ω 1 is marked for reﬁnemen t, and (b) the ﬁnal mesh with the newly added edges shown as dashed lines. The LEPP for Ω 1 is trac ked and triangles are bisected in backw ard manner till Ω 1 is reﬁned. Algorithm 2: Mesh adaptivit y Algorithm for Ω i ∈ E ref do while Ω i isn ’t bise cte d do Find the LEPP of Ω i ; if The longest e dge of Ω n , the last triangle in LEPP of Ω i , is on Γ o then Bisect Ω n ; else Bisect the last tw o triangles in LEPP of Ω i , ( Ω n , Ω n − 1 ); end end Mo ve eac h rep no de to the nearest netw ork no de; if Ω i is ful ly r esolve d then T ransform from QC elemen t into discrete elements; end end After mesh reﬁnement, each represen tative no de is mov ed to the nearest netw ork no de. All fully reﬁned QC elemen ts (i.e. elemen ts that contain only 3 netw ork nodes) are transformed from 2D 9 triangular elements to 1D links consisten t with the underlying netw ork structure. This step also prev ents ill-conditioning of the stiﬀness matrix of the QC system [ 33 ]. 4.2 Time adaptivity: Time adaptivity is required for computational eﬃciency by allowing larger time steps, when p os- sible, without compromising accuracy . In addition, time adaptivit y is also critical for the accurate mo deling of fracture in high interest zones where small time steps are needed. Using a constan t time step may b e ov erly restrictive if a small time step is used throughout and may lead to inaccu- rate results if a large time step is p erp etually maintained. Time adaptivit y balances computational cost and error. The sim ulation starts with a stable large time step, that accurately and eﬃcien tly compute the nonlinear elastic deformation. The algorithm has tw o time reﬁnement criteria whic h are usually triggered when the discrete elemen ts approach failure. The ﬁrst one is to ensure con- v ergence where the time step is reduced, if needed, to reach the stable time step required for the solution conv ergence. The second criterion is to ensure accuracy in predicting the crack propa- gation as large time steps ma y lead to ov erestimating the n umber of failed links. Once failure initiates, the time step is reduced progressiv ely such that further reduction in the step size won’t aﬀect the num ber of elemen ts failing at the curren t time step. In the curren t study , the time step is the same for both the homogenized and fully discrete domains. W e disc uss possible alternative approac hes in Section 7 . 4.3 Nonlinear ﬁnite elemen t framework: The algorithm used for adaptiv e mo deling of fracture in p olymer net w orks is given in Algorithm 3 . The homogenization is p erformed on the ﬂy to pro vide the stress and the material tangent tensor eac h time step. Algorithm 3: A daptive QC based ﬁnite elemen t framework Initialize the system, construct the initial mesh with required information; for step = 1,2,......,n do Calculate the homogenized material consistent tangen t and stress tensor using Algorithm ( 1 ); Apply the step incremental boundary conditions at step n; Assem ble the system in ternal force vector and tangent matrix and solv e for the displacemen t; F or each 2D element chec k the mesh adaptivity conditions, add all elements marked for reﬁnemen t to a set E ref ; if E ref is not empty then Reﬁne the current mesh as p er Algorithm 2 , up date system information; Rep eat the curren t time step to balance the new mesh system; else Store the relev an t outputs for the curren t time step; Pro ceed to the next time step; end end 5 V eriﬁcation F or veriﬁcation of the prop osed homogenization tec hnique, we run multiple tests for b oth pristine and crack ed samples. In each case, we v erify the QC results by comparing it to the results of the fully discrete mo del and quantify the error in terms of the measured reaction normalized force at diﬀeren t levels of stretc h. W e deﬁne the normalized force as f b/k B t (Eq. 30 ). 10 5.1 Uniform deformation gradien t tests: F or the ﬁrst test, we consider a regular netw ork geometry sho wn in Fig. 3 -a with uniform material prop erties. All the netw ork links in the sample are homogenized using the pro cedure outlined in Section 4 . W e impose an aﬃne deformation gradient, considering b oth uni-axial and bi-axial tension loading cases, up to a stretch v alue of 6. The QC results match exactly the fully discrete mo del in terms of the reaction force and the total strain energy v alues. This v eriﬁes that the homogenization sc heme reco vers the exact response if the microscale deformation is aﬃne and the net work structure is regular. Figure 3: The pristine sample used for veriﬁcation: (a) uniform netw ork (in black) and the corre- sp onding QC mesh (in blue), (b) non uniform net work where the nodal location is shifted randomly from the uniform case, the maximum shift is controlled by the disorder parameter µ g , and (c) a QC elemen t showing the net work nodes inside this element and the repnodes used for interpolation. Next, to study the eﬀect of disorder in the geometric and material prop erties on the p erformance of the QC metho d, w e sim ulate th e resp onse of a non uniform netw ork, as shown in Fig. 3 -b and the zoomed in view in Fig. 3 -c, under uni-axial tension. W e introduce non-uniformit y in t wo wa ys. First, we shift the p osition of the netw ork no des randomly with resp ect to a reference perfectly ordered lattice. The maxim um shift v alues in x and y directions are µ g L x and µ g L y , respectively , where µ g is the geometric nonconformity parameter. Second, we dra w the p olymer c hain con tour length L c randomly from a uniform distribution ranging betw een [(1 − µ ) L c , (1 + µ ) L c ] , where µ is the contour length nonconformity parameter. F or eac h non-uniformit y parameter v alue, ten randomly generated meshes are simulated and the a verage results are plotted in Fig. 4 . W e deﬁne the relative error in the force as: E RR f = | P k R D k − P k R QC k | P k R D k (35) Where P k R k is the sum of reactions in the loading direction, D refers to the corresponding fully discrete simulation and QC refers to the QC simulations containing homogenized contin uum elemen ts only . Fig. 4 -a shows that as the disorder increases the error p ercen tage from the homogenization sc heme increases. Ho wev er, this error decreases as the fraction of repno des in the sample increases. This is further conﬁrmed in Fig. 4 -b whic h sho ws that v ariation of error as a function of v ariability in the chain lengths at tw o v alues of repnodes percentages. A t low v alues of disorder, µ less than 0.2, b oth the homogenized and discrete models give identical results (the error is almost zero). As the disorder increases, the error increases. How ever, it remains less than 8 % and it decreases with mesh reﬁnement. Fig. 4 -c shows that the error from the homogenization dep ends on the stretch lev el. F or lo w v alues of disorder, the error is negligible even at high stretch v alues (up to 6). As the disorder increases, the error is relatively small (few p ercen ts) at low stretch v alues but it rapidly increases at high stretch v alues reaching as high as 15 % for µ =0.5 at stretc h v alue of 5. While this exercise attests to the accuracy of the homogenization scheme at low stretch v alues, it suggests that mesh reﬁnement is crucial for retaining accuracy at high lev el of stretch. This will b e accoun ted for automatically in the QC application b y turning on dynamic adaptivit y and enforcing a reﬁnement criterion based on the lo cal stretch level as explained in Section 4 . 11 Figure 4: Homogenization error of the QC method applied to disordered netw orks: (a) the relativ e error in the force at λ = 1 vs. the ratio of the repno des to the netw ork no des for 3 diﬀeren t disorder parameter µ = 0 . 05 , 0 . 25 , and 0 . 5 , (b) the relative error in the force at λ = 1 vs. µ for tw o v alues of repno des ratio, and (c) the evolution of the error vs. stretch for µ = 0 . 05 , 0 . 25 , and 0 . 5 and repno des ratio = 7.5%. 5.2 Non-uniform deformation gradien ts– Eﬀect of reﬁnemen t: Here, w e run tw o tests where the deformation gradient is not uniform ov er the domain. In the ﬁrst test, w e model a pristine sample with dimensions 64x64 where the b ottom edge is restrained from mo vemen t in b oth directions and the top edge is clamp ed and is pulled up ward. The results sho wn in Fig. 5 using only contin uum elemen ts, homogenized through the QC pro cedure, sho w that although the force vs. stretch shows go o d agreement with the full discrete model, up to a stretc h v alue of 6, there is a discrepancy in the deformed shape specially near the edges. The error distribution plot also sho ws that the no dal displacemen t error is highest in the edge elements, closest to the free lateral b oundaries, where the netw ork no dal displacements are in terp olated from the edge no des of the contin uum elements. T o address this problem, we designate these areas as high in terest areas and mo del the net w ork links in these regions explicitly . This leads to eliminating the concentrated errors near the edges. Although the repno des ratio increased from 0.5% to around 5%, the deformed shap es and the force vs. stretc h curves using b oth QC and fully discrete matc h b etter in this case. F urthermore, the error in the displacement distribution falls almost an order of magnitude throughout the domain. In practical applications of the QC metho d, the discov ery of high interest areas and the change in the mo del resolution is done automatically using the adaptiv e mesh reﬁnement algorithm. In the second test, w e use the same geometry and b oundary conditions as in the ﬁrst one, but the top edge now mov es b oth upw ard and to the right such that the sample is sub jected to b oth tension and shear. In this case, mesh reﬁnemen t reduces the error in no dal displacemen t and leads to a b etter agreemen t in the force vs. stretch resp onse betw een the QC and the full discrete mo dels. The results shown in Fig. 5 and Fig. 6 illustrate the capability of the the QC metho d in bridging diﬀerent scales, and also highligh ts the critical role of mesh adaptivity in reducing the errors sp ecially in areas where the deformation gradien t is not uniform. F urthermore, since the results in Fig. 5 , and Fig. 6 suggest that the error in the displacemen t is usually high near the free surfaces, we hav e chosen to alw ays use an initial QC mesh that is reﬁned around the sample b oundaries in the remaining cases presented in this pap er. 12 Figure 5: The eﬀect of reﬁnemen t on the QC metho d applied to pristine netw orks under tensile loading: The net work and the QC mesh, deformed shap e, relative error in the no dal displacement, and force vs. stretc h from b oth discrete and QC results for (a) QC mesh where no high in terest area (i.e. no discrete links) is deﬁned, and (b) QC mesh where the elements adjacen t to the free v ertical edges are reﬁned to the discrete representation. Figure 6: The eﬀect of reﬁnement on the QC method applied to pristine netw orks under com bined tensile and shear loading: The netw ork and the QC mesh, deformed shap e, relativ e error in the no dal displacement and force vs. stretc h from both discrete and QC results for (a) QC mesh where no high interest area (i.e. no discrete links) is deﬁned, and (b) QC mesh where the elements adjacen t to the free edges are reﬁned to the discrete represen tation. 5.3 F racture test: Here, w e mo del crack propagation in a notched netw ork using both fully discrete and QC models. The netw ork dimensions are (500x250) chain units with around 250,000 degrees of freedom and the notch is 250 units long. The sample dimensions and b oundary conditions are shown in the insert in Fig. 7 -a. All links follow the w orm like chain mo del, with force-stretch relation given by Eq. 30 . The chain length is drawn from a uniform random distribution with µ = 0 . 05 . A link is set to break if the v alue of the stretch λ p exceeds 0.8. Automatic mesh adaptivty is turned on to allo w the high in terest area to ev olve with the crac k propagation. 13 Figure 7: V eriﬁcation of the QC metho d for a fracture test: (a) the normalized force vs. stretc h from QC and discrete results, the insert sho ws the geometry of the sample, (b,c,d) mesh ev olution and crack path from the QC results, and (d) crack path from the discrete mo del results. The results in Fig. 7 sho w that the crack path and the force vs. stretc h curv es are almost identical for b oth QC and fully discrete mo dels. Sp eciﬁcally , the error in the p eak force is around 1.0%, the error in the peak stretch is 1.4%, and the error in the total area under the curv e is under 0.5%. The QC mesh evolution is also sho wn in the same ﬁgure. The ﬁrst mesh Fig. 7 -b has a rep. no des ratio of only 2.1%, whereas the ﬁnal mesh Fig. 7 -d has 5% rep. no des ratio. This signiﬁcan t reduction in the num b er of degrees of freedom lead to extreme computational saving. 5.4 Crac k initiation: In this sec tion, we verify the capability of the prop osed QC methodology to capture crack initiation in an in tact sample and subsequen t crack propagation. F ormulating a criterion for crack nucleation in con tinuum mo dels ma y b e c hallenging since initiation of fracture in an in tact sample ma y b e based on strength rather than energetic consideration. As the QC blends both discrete and con tinuum representations, it does not suﬀer from such limitation. Fig. 8 shows the comparison b et w een the results of a fully discrete simulation and the corresp onding QC mesh for a pristine sample with dimensions (250x250) where the links in the zone shown in Fig. 8 -a hav e b een w eakened b y lo wering the failure stretc h to 50% of that of the netw ork links. The QC metho d captures the crac k initiation thanks to the link stretc h based mesh reﬁnemen t criterion. The force-stretc h relation for b oth cases are nearly identical with errors in p eak force and stretc h b elo w 0.5%. 14 Figure 8: Crack initiation and propagation due to a zone with w eak links: (a) The normalized force vs. stretc h from QC and discrete results, (b) crack path from fully discrete simulation results at p oin ts 1, 2, and 3 resp ectiv ely , and (c) mesh evolution and crack path from the QC results at p oin ts 1, 2, and 3 resp ectiv ely . 5.5 Computational time analysis: In this section, we compare the run time for the prop osed QC approach vs. that of fully discrete samples for the same geometry in Fig. 7 and diﬀerent sample sizes. F or solving the system, w e use conjugate gradient solver with incomplete Cholesky preconditioning [ 50 ] whic h is suitable for symmetric p ositiv e deﬁnite matrices. The computational times of an av erage step and the total time of a simulation of 500 steps are shown in T able 1 . The time p ortions consumed in solving the system of equations at each step, including all iterations, are shown separately . Although QC requires more time at the initialization to construct the course mesh and assign links to QC elemen ts, the total simulation time signiﬁcan tly reduces compared to the fully discrete sim ulations. The reduction is due to solving smaller system each time step, the time consumed in solving the system is reduced b y nearly tw o orders of magnitude in the cases shown in T able 1 . In addition, the time consumed in the homogenization is less than the time required to calculate the local stiﬀness matrix for eac h link and assem ble the global stiﬀness matrix in the fully discrete simulation. Thus, the total simulation time is reduced up to 15 times for the cases shown. T able 1 also shows that as the sample size increases, the ratio of the representativ e no des decreases and the sp eed up ratio b et w een QC and fully discrete sim ulation increases which suggests that computational adv antage of the method is even larger for larger scale problems. It is also worth noting that homogenization step is done indep enden tly for each QC element, th us it may be easily parallelized with almost ideal exp ected speed-up [ 36 ]. T able 1: Computational times consumed by fully discrete vs QC metho d for diﬀeren t sample sizes 6 Results In this section, we apply the QC metho d to in vestigate the eﬀects of lo cal top ology and disorder on the netw ork mechanical prop erties and its fracture c haracteristics. W e show the adv an tage of the QC metho d in capturing those eﬀects without the need to explicitly mo del every no de of the netw ork. As explained earlier, the discrete representation is limited to the zone around the 15 crac k whereas the regions a wa y from the crac k are homogenized consistently . This introduces a computationally eﬃcient accurate representation of the fracture pro cess. 6.1 Damage evolution in samples with diﬀeren t disorder parameter: Figure. 9 shows the damage evolution and the crack path for the single notched sample setup shown in Fig. 7 -a with t wo diﬀeren t disorder parameter ( µ ) v alues and Z ≈ 8 , where Z is the co ordination n umber deﬁned as the num b er of connections a node may hav e to other no des in the netw ork, a veraged o ver all the no des. Figure 9: Damage ev olution in a single-notc hed sample, the sample geometry is shown in the insert of Fig. 7 -a: (a) net work with lo w disorder, and (b) net w ork with high disorder. The damaged lattices are marked by red color, whereas the repno des are shown as blue dots. F or µ = 0 . 05 , whic h represents an almost uniform net work, the crack propagates coheren tly per- p endicular to the applied loading direction. The crac k path is nearly horizontal, and the damage is lo calized along the crack path. F or highly disordered netw orks with µ = 0 . 50 , although the dom- inan t crac k propagation direction is , on av erage, perp endicular to the applied loading direction, the crac k path has more kinks as it propagates. F urthermore, the crack thickness, which may b e tak en as a basis for a material length scale in nonlo cal theories such as phase ﬁeld metho ds, is not set apriori. It emerges as part of the solution and shows non-uniformity as it tends to increase at some lo cations and decreases in others. Damage is not lo calized around the crack path, as the ﬁgure shows m ultiple local damage zones aw a y from the crack tip. Some of these lo cal damage sp ots coalesce with the main crack while some remain isolated. Fig. 9 also shows the increase in the n umber of repno des, represented b y the blue dots, with crac k propagation due to mesh adaptivit y . The rate of increase in the num ber of repno des for more disordered netw ork is higher than that in the more uniform ones. This attests to the ﬂexibility of the QC metho d in adapting to the complexit y of the underlying netw ork structure. 6.2 F orce distribution in the vicinit y of the crack tip: T o gain further insigh t into the pro cesses controlling crack initiation and growth, w e analyze the force distribution in the discrete netw ork in the vicinit y of the crack tip for diﬀeren t v alues of the co ordination num ber Z and disorder parameter µ . As deﬁned earlier, the co ordination num ber is the a v erage n umber of no des directly connected to a giv en no de and measures ho w densely connected the netw ork is. 16 Figure 10: The force distribution in the net work links around the crack tip: (a,b,c) net works with Z ≈ 8 and µ = 0 . 05 , 0 . 25 , and 0 . 5 resp ectiv ely , (d,e,f ) netw orks with Z ≈ 6 and µ = 0 . 05 , 0 . 25 , and 0 . 5 respectively , and (g,h,i) netw orks with Z ≈ 4 and µ = 0 . 05 , 0 . 25 , and 0 . 5 resp ectiv ely . The results in Fig. 10 show how c hanging the lo cal connectivity changes the force ﬁeld around the crac k altering conditions for the crack initiation, propagation, and growth direction. F or example, for Z ≈ 8 and µ = 5% , the stress ﬁeld has the classical p ean ut-like distribution as exp ected for Mo de I fracture in planar elasticit y [ 51 ]. On the other hand, for Z ≈ 4 and µ = 5% , the anisotropic net work channel the forces in the links that are closely aligned with the loading direction. F or Z ≈ 6 , the force distribution shows an interesting pattern with some increase in the forces along diagonal directions running bac kward relativ e to the crack tip. F or eac h co ordination n umber, increasing the disorder in the chain length signiﬁcantly changes the force ﬁeld. In particular, the force ﬁeld b ecomes more diﬀuse and loses its symmetry . F urthermore, increased disorder ma y lead to elev ated force magnitudes in regions aw a y from the crac k tip leading to a more complex fracture pattern and p oten tially thick er cracks as was observ ed in Fig. 9 . Suc h complex force ﬁeld distribution in highly anisotropic and disordered net w orks suggest the p ossibilit y of p oten tial crac k gro wth in unusual directions as has b een observ ed recently as sidew ay cracks in silicone elastomers [ 52 ]. The QC metho dology may enable systematic disco very of these new patterns. 6.3 Eﬀect of disorder in the c hain length: W e sim ulate the disorder in the net work by drawing the c hain length from a sto c hastic distribution. Fig. 11 shows the resp onse of a netw ork with a single notch, the same setup sho wn in Fig. 7 -a, when the a v ailable c hain length is drawn from a random uniform distribution ranging b et ween [(1 − µ ) L c , (1 + µ ) L c ] where L c is a predeﬁned contour length. Fig. 11 -a-c show the eﬀect of the 17 c hain length v ariabilit y on the crack path. As the v ariability increases, the crack path b ecomes more meandering while propagating through the net work. The QC mesh also shows higher reﬁnemen t (high repnodes ratio) as the disorder increases. This is because as the range of v ariability in the c hain length increases, some links with shorter conto ur length b ecome critically stressed at lo cations not at the crack tip requiring further reﬁnemen t in these regions. Figure 11: The eﬀect of net work disorder on the mechanical resp onse of single-notched samples. The sample geometry is sho wn in the insert of Fig. 7 -a: (a,b,c) crack path for single notched sample and 3 diﬀerent disorder parameter µ = 0 . 05 , 0 . 25 , and 0 . 5 resp ectiv ely , (d) the force vs. stretc h relation for the three disorder parameter v alues, and (e) the eﬀect of µ on the peak force and ductilit y of the net works, the uncertaint y in these v alues d ue to the disorder is indicated b y the error bars. T o further in vestigate the eﬀect of disorder, the a v erage force vs. stretch curves are plotted for diﬀeren t disorder parameters µ in Fig. 11 -d. Eac h curv e is the av erage of num b er of randomly generated net works with the same disorder parameter. The num b er of netw orks was chosen such that the change in b oth the mean and the standard deviation of the p eak force distribution due to adding more net w orks falls b elo w 1%. W e performed this analysis for the highest disorder parameter which required 60 netw orks, then the same n umber of net works was used for all cases. More details ab out this pro cedure may b e found in App endix A. The results show an increase in the stiﬀness and a decrease in the p eak force as the disorder increases. The increase in stiﬀness is in trinsic in the nature of force stretch curv e of the single worm lik e chain where the force increases rapidly as the stretc h ratio approaches the stretc h failure threshold. F or more disordered net works, the increase in the force in the shorter links dominate ov er the reduction in force in the longer links. F or the same reason, net w orks with higher v ariability in chain length tends to initiate failure at lo wer p eak force, as they hav e a larger n umber of shorter links which fail faster. The results also sho w that the ductility (i.e. the ratio b et ween stretch at total failure and stretc h at p eak force) increases with increasing µ . In less disordered net works, the links are more uniformly stressed and hence they fail rapidly once failure initiates. As the v ariability in the chain length increases, the net work has more longer c hains which helps the material system to surviv e longer under damage, this agrees with the exp erimen tal results rep orted in Bonn et al. [ 53 ]. Fig. 11 -e sho ws that the av erage p eak force decreases with increasing µ , whereas the a verage ductility increases. F urthermore, the uncertaint y in b oth p eak force and ductilit y , represented b y error bars, increase with increasing the material disorder. Figure. 12 shows the eﬀect of increasing the net work disorder on the crack path of a netw ork with 18 Figure 12: The eﬀect of net w ork disorder on the crack path of double-notc hed samples: (a) the geometry of symmetric double notc hed sample, (b,c,d) crack path for symmetric double notched sample and 3 diﬀerent disorder parameter µ = 0 . 05 , 0 . 25 , and 0 . 5 resp ectively , (e) the geometry of asymmetric double notc hed sample, and (f,g,h) crack path for asymmetric double notched sample and 3 diﬀerent disorder parameter µ = 0 . 05 , 0 . 25 , and 0 . 5 resp ectiv ely . symmetric and asymmetric double edge notc hes. As the disorder increases, crack path tends to b e more nonuniform and thick er. In addition, a more extended spatial region is reﬁned as the v ariability in the c hain length increases, to capture the lo cal failure of shorter chains in lo cations a wa y from the crack tip. F urthermore, these results sho w how the t wo cracks emerging from the double notches in teract with each other at diﬀeren t µ v alues. The crac k path shap e for the case of µ = 0 . 05 was veriﬁed by running fully discrete sim ulations for the same geometry . The fracture pattern agrees qualitativ ely with previously reported exp erimen tal and n umerical observ ations [ 54 , 55 ] .It also agrees qualitativ ely with the numerical results for elastomeric materials fractured b y cross-link failure using the phase ﬁeld metho d [ 9 ]. Ho w ever, our results suggests that local top ology ma y hav e a profound eﬀect on this crack path that may be c hallenging to capture using fully homogenized approaches. In particular, the curving of the cracks tow ards one another ma y completely disapp ear as the disorder increases. 6.4 Eﬀect of co ordination n um b er Z: In this section, we inv estigate the eﬀect of net work connectivit y through changing he co ordination n umber Z . Fig. 13 sho ws the force vs. stretch for three diﬀeren t co ordination num b ers, Z ≈ 4 , Z ≈ 6 and Z ≈ 8 . In constructing these net works, w e keep the total w eight of the netw ork ,i.e. the 19 sum of all p olymer chain av erage con tour lengths, as a constant. W e also dra w the con tour length of the diﬀerent c hains from the same uniform random distribution with µ = 0 . 25 for the three cases. The results sho w that net works with higher coordination num b er has higher toughness and p eak force compared to net w orks with low er co ordination n umber. As the co ordination n umber increases, more links con tribute to load sharing and hence the net work tends to become stiﬀer and sustain higher load. Similar conclusion was recently rep orted b y Alame and Brassart [ 20 ]. F urthermore, the case with higher coordination n umber shows a larger ductility ratio, deﬁned as the ratio betw een stretc h at failure and the s tretc h at p eak force, but a slightly smaller maximum stretc h. As the co ordination n umber increases, there are, on av erage, more links connected per no de. This enables a more eﬃcient force redistribution among neighboring chains when one of the c hains break oﬀ, leading to slow er force drop and delay ed damage accumulation. Similar results w ere obtained for nano-composite gels with nanoparticle cross-linkers that increase the netw ork co ordination in teractions [ 56 ]. As the co ordination n umber decreases, the material b ecomes more complian t, and is capable of ac hieving sligh tly higher stretch v alues although with lo w er o v erall toughness and strength. Sp eciﬁcally , the area under the curve for Z ≈ 6 is 1.5 times that for Z ≈ 4 and the area under the curve for Z ≈ 8 is twice that for Z ≈ 4 . Figure 13: Eﬀect of net w ork top ology on its mec hanical response: (a) the eﬀect of the co ordination n umber Z on the force vs. stretc h b eha vior, and (b) the eﬀect of cross-linking densit y , measured relativ e to a reference mesh, on the net work p eak force for three v alues of p olymer chain av ailable con tour length L c . 6.5 Eﬀect of cross-linking densit y: T o inv estigate the eﬀect of increasing cross-linking density , deﬁned as the n umber of net work no des p er unit area, w e simulate multiple samples with increased num b er of no des and links in the the net work while k eeping the p olymer concen tration constan t. This is ac hieved by adjusting the a verage chain contour length such that the sum of lengths of all c hains remains constan t. The results in Fig. 13 -b sho ws that the p eak force increases with increasing the cross-linking density up to a maximum v alue, then it starts decreasing. This non monotonic response is due to the reduction of the av erage chain length with increasing the cross-linker concentration. When the a verage chain length decreases, the netw ork b ecomes stiﬀer and attains higher forces for a given v alue of stretch. How ev er, as the c ross-link er concentration further increases beyond a certain v alue, which dep ends on the chain length, the netw ork b ecomes critically stressed as the initial stretc h of the short chains becomes closer to their contour length. These short c hains accumulate damage at a faster rate leading to a premature failure at low er global force lev els. With increasing the a v ailable c hain length, the peak force for a giv en concen tration increases, and the critical cross-link er concentration corresp onding to this p eak force also increases. 20 7 Discussion: F racture is a highly nonlinear problem in whic h macroscopic resp onse dep ends on the details of microscopic physics. Here we hav e in tro duced an adaptive mo del based on the Quasi-Contin uum approac h for analysis and design of netw orked materials that starts from fundamental link scale mec hanics and consistently computes macroscopic resp onse. While this w ork builds on recent dev elopments of extending the quasi-con tinuum approac hes to lattice materials, it brings new inno v ations in b oth the application domain as well as the metho d formulation. F or example, on the application level, a ma jor contribution of this study is to bring the QC metho d into the discussion of mo deling fracture in polymeric and soft materials. This is a topic of intense con temp orary in terest whic h lends itself naturally to the QC framework but has been addressed in the solid mechanics comm unity , for the most part, using con tinuum mo dels or phase ﬁeld approaches that generally o verlook the underlying topological structure or ha ve limited predictive p o wer. F urthermore, a unique adv antage of the QC metho dology is that it allows capturing crack initiation and nucleation, a strength-based process, that other approac hes built on mo diﬁed Griﬃth criterion, ma y fail to capture in a straigh tforward wa y . On the metho d formulation, due to the highly nonlinear nature of the p olymer netw orks, we hav e prop osed a homogenization rule for axially loaded lattices that accoun ts for both material and geometric nonlinearities and implemented it in a nonlinear ﬁnite elemen t framew ork. Moreo ver, w e hav e developed an adaptivity sc heme that allo ws for transition b et w een homogenized material and net work lattices where the reﬁnement criteria dep end not only on the deformation gradient, as has b een done in earlier studies but also on link stretch. This is necessary , particularly for highly di sordered netw orks, and also to allow for triggering reﬁnement and crack nucleation at the lo cations of w eaker links in absence of a notch or pre-crac k While methods for mo deling the nonlinear elastic resp onse of p olymeric materials are well es- tablished [ 57 , 6 ] , modeling fracture and damage in soft p olymeric materials con tin ue to elicit c hallenges in science and engineering [ 58 , 59 , 60 , 61 , 62 , 63 , 64 , 65 , 66 ]. How ever, despite enabling useful insights, most of the existing techniques suﬀer from a n umber of limitations. F or example, most con tinuum fracture mo dels are insensitiv e to lo cal topological or micro-structural details at the net work level. They also require input for critical parameters such as fracture energy whic h ma y not b e easily measured in many soft materials or ma y exhibit size eﬀects [ 67 ]. F or phase ﬁeld metho ds, and other non-lo cal theories, a length scale is also needed for regularization of the frac- ture [ 68 ]. While in the case of phase ﬁeld approaches, this length scale is introduced for numerical purp oses and the framew ork should con verge to Griﬃth form ulation up on taking this length scale to the zero limit, it is not clear ho w this length scale is to b e c hosen in cases when the material has an in trinsic one arising from its micro-structure. This compromises the predictive p o w er of these mo deling approaches. While fully discrete sim ulations pro vide a p o werful alternative, as they capture the full details of the material microstructure, the computational cost of these approaches b ecome prohibitive as the sample size increases. The main adv antage of the QC framework is that it enables eﬃcien t mo deling of large scale prob- lems without neglecting the eﬀect of lo cal topology and without assuming any critical macroscopic parameters such as fracture energy , or material length scale. Only a lo cal failure criterion at the c hain level is needed and this may b e rigorously computed directly from explicit mo del of a single c hain [ 47 , 3 , 9 ]. The QC metho d requires m uch less run time, and memory than a full discrete mo del while providing higher accuracy and insights into the fracture pro cess than the contin uum approac hes. F urthermore, the QC has the ﬂexibility of adding small scale ph ysics and top ology . This makes it suitable for taking adv antage of recent adv ances in exp erimen tal techniques. While most curren t exp erimental techniques do not allo w direct visualization of the deformation and damage on the netw ork micro scale, there hav e b een a few exciting adv ances in recent years to wards that goal. F or example, the work of Ducrot et al. [ 69 ] where a sp ecial t yp e of mechanophore molecules is em b edded in the polymer net work and enables the direct imaging of the spatio- temp oral ev olution of damage density due to b ond break age. Also the work of P oulain et al. [ 65 ] where high sp eed photography is used to track the onset and gro wth of cavitation and subsequent coalescence and crack gro wth. These recent adv ances in imaging and exp erimen tal measurements among others are probing deformation and failure mec hanisms in soft materials with increasingly high resolution challenging some long-held ideas ab out these phenomena [ 70 , 71 ]. The prop osed QC approach is uniquely suited for taking adv antage of these adv ances and making use of the increasingly a v ailable data at the microscale. F urthermore, the QC approach provides a rigorous 21 to ol for linking the micro and macroscales with minimal set of assumptions whic h ma y p otentially help resolve several c hallenges in failure of hierarc hical and heterogeneous polymeric materials. It is worth to note that there exists a large literature on discrete mo dels of fracture that abstract a contin uum domain using a lattice represen tation, suc h as Lattice Discrete Particle Mo deling [ 72 ], Graph-based Finite Element Analysis GraFEA [ 73 ], and Virtual b ond mo del [ 74 ]. A ma jor diﬀerence b et w een the QC approach and these metho ds is that the discrete comp onent of the QC mo del is physical and represents an actual netw orked material with a well deﬁned failure criterion at the discrete elemen ts scale. Thus, the QC metho d do es not suﬀer from some limitations that exist in these other lattice mo dels such as mesh sensitivity , p ossible ambiguit y in deﬁning the fracture criteria up on abstraction of the contin uum scale to the lattice scale, or the necessit y of in tro ducing non-lo cal eﬀects [ 68 ] or phenomenological length scales. In particular, the crack path computed from the QC mo del is not mesh dep enden t, addressing an important challenge in modeling fracture of inhomogeneous materials [ 75 ]. Our results for fracture of dis ordered net w orks suggest that as the disorder increases, the expected v alue of the p eak force decreases but the ductility ratio and energy dissipation increase. It is also observ ed that the crack path b ecomes non uniform with increased disorder, as the fracture mec hanism includes numerous links failing in the netw ork aw a y from the crack tip zone and later coalescing in to a thick crac k path (Fig. 9 -b). Meandering crac k paths suggest a more ductile resp onse, which is observed in Fig. 11 -d,e and increases fracture toughness compared to straight crac ks. Suc h rough crac ks are observ ed extensively in exp eriments [ 76 , 77 , 78 ] and are usually attributed to dynamic instabilities. Here we show that they ma y also emerge naturally due to the increased lev el of disorder in the netw ork structure. T ransition in the failure mo de from brittle to ductile with increased disorder has been reported previously using fully discrete mo dels on smaller samples [ 79 ]. The QC metho d pro vides a ﬂexible to ol for reassessing these phenomena at larger scales and with less sensitivity to sample b oundaries. The current formulation of the QC implementation has a few limitations. F or example, w e ha ve assumed an aﬃne pro jection approximation to couple micro and macro-scales in the homogenized region of the domain. While w e hav e shown that the homogenization error is within acceptable limits for the cases tested, this aﬃne approximation has some disadv antages. F or instance, the accuracy was shown to deteriorate as the disorder increases unless mesh adaptivity is turned on. In addition, w e only considered net works with axial in teractions. Some p olymer chain models hav e a contribution to their p oten tial energy originating from the rotational degrees of freedom which in tro duces a non-aﬃne comp onent in the deformation ﬁeld [ 16 ]. While an y force-stretch la w ma y b e readily used in the current form ulation in place of the worm-lik e chain mo del, future w ork will explore extension of the homogenization algorithm to accoun t for non-aﬃne deformation. In this initial study we opted to ev aluate the stress and material tangen t tensors from the con tri- butions of all links underlying eac h contin uum elements. It is p ossible, how ev er, to determine a represen tative volume element (R VE) and estimate the stress and tangent tensors from the con- tribution of links within this v olume. The general form of the tensors would be the same, but only links within this R VE will b e accounted for. While this homogenization approac h leads to increased computational sa vings, it also has some drawbac ks. F or example, it is not clear apriori ho w large the R VE should be in the case of disordered netw orks. Also the approximation based on an R VE will in tro duce additional errors. The pros and cons of this alternative metho dology will b e inv estigated further in future w ork. Automatic mesh adaptivity leads to signiﬁcant computational savings by reducing the initial fully resolv ed zone. Ho wev er, as the crac k propagates, the ratio of repno des increases, especially in cases of highly non-uniform netw orks, leading to reducing the computational savings. It is p ossible to implement coarse-graining behind the crac k pro cess zone to reduce the num b er or repno des when the full resolution of a speciﬁc zone is not further needed. QC with coarse graining was accomplished by Rokos et al [ 48 ] were signiﬁcant computational savings were rep orted. In addition, time adaptivity allows taking larger time steps, when p ossible, without compromising accuracy . In the current study , we use the same time step for the whole domain. The choice of time step is based on the s mallest stable v alue p ossible in the global assem bly . How ever, it is p ossible to adopt an asynchronous time stepping approac h where diﬀeren t time steps for the fully discrete and homogenized domains are used and couple them in a parallel framework [ 80 ]. This might lead to signiﬁcan t time sa vings as the fully discrete domain usually requires small time steps to resolv e 22 the fracture accurately and will b e explored further in a future study . Visco elasticit y is a critical feature in the mechanical resp onse of p olymeric materials such as rubb er and gel. In our prior w ork [ 47 , 3 ], we hav e mo deled discrete p olymer net works with visco elastic resp onse emerging from t w o ma jor contributions: (1) break age and formation of end b onds and cross-link ers, and (2) vsicous drag arising from the relativ e motion of the deforming polymer net work and infused ﬂuids. Extensiv e work has also b een done previously by other groups for exp erimen tal characterization of the elastic and loss mdo dulii [ 81 , 82 ]. Con tinuum mo dels based on phenomenological as well as homogenization approac hes also exist [ 83 , 84 ]. The QC metho d is capable of incorp orating the visco elastic resp onse. The key p oin t is to map the visco elastic b eha vior from the microscale (netw ork level) to the macroscale (QC level) using a prop er homogenization rule. As a starting p oint, one may assume a damping force in eac h link that is prop ortional to the link stretc h rate in addition to damping force at each netw ork no de prop ortional to its v elo cit y . Then this damping matrix ma y b e homogenized in a similar wa y as the stiﬀness comp onen t. The ﬁrst term will result in a damping matrix that has the same structure as the stiﬀness matrix whereas the latter will result in a diagonal damping matrix. W e will then use predictor-corrector sc hemes to solve the resulting lin earized quasi-dynamic system at each time step. This approach w ould lead to the emergence of visco elasticit y on the contin uum scale. In this initial w ork, we considered polymer netw orks in t wo dimensions. Ho wev er, the metho d ma y b e easily extended to mo del net works in the 3D domain. Unlik e in atomistic systems, where the QC computational cost ma y exponentially increase from 2D to 3D due to the presence of long range in teractions, the extension in lattice systems is straightforw ard since ph ysical links introduce in teraction b et ween directly connected nodes only . F urthermore, extension to 3D will enable incorp oration of additional physical mechanisms suc h as en tanglement, whic h may b e mo deled at the discrete level as an additional type of connections b et ween the links, suc h as ring sliding connections [ 85 ] .This ma y then b e homogenized using similar metho dology to the one presen ted in this work and passed on to the contin uum scale. T o summarize, the QC framew ork enables signiﬁcant insight in to the role of microstructure on macroscopic fracture resp onse in net work ed materials. The unique p oten tial of the metho d lies in its predictive p o wer which qualiﬁes it to pla y the role of a highly con trolled virtual experimental setup. This will not only lead to iden tiﬁcation of fundamen tal processing con trolling fracture and deformation in a v ariet y of net works, including p olymers, but ma y also enable the discov ery of new material designs with unusual fracture response. 8 Conclusions: Our main conclusions are summarized as follows: 1. QC metho d pro vides a promising platform for inv estigating fracture in disordered net work ed materials that combines adv antages from b oth discrete and con tinuum approac hes. 2. The QC metho d for fracture in polymer net works only requires a failure criterion on the c hain lev el. It do es not require information ab out critical energy release rate or a material length scale. The QC discov er these quantities as part of the solution. This predictiv e p o wer and scale bridging capabilities giv e the QC framew ork an adv antage o v er other existing methdologies. 3. A daptivity allows the QC mesh to ev olve as the crack grows, increasing accuracy and com- putational savings. F or the cases presented here, we ha ve achiev ed up to 15 times sp eed up compared to the full discrete models. Ongoing work with parallel implemen tation will enable bridging a wider range of scales with further computational savings. 4. Crac k path is con trolled b y the net w ork top ology and sto c hastic prop erties of the c hains. Th us, QC metho d illuminates new pathw a ys for design through direct control of the mi- crostructure. 5. The details of the netw ork lo cal top ology hav e a signiﬁcant eﬀect on the local force distri- bution around the crack tip. Disorder mak es the crack tip force ﬁeld more diﬀused and also c hannel forces in directions dominated by netw ork anisortropy . 23 6. Net works with higher disorder in the geometry or/and c hain length hav e higher stiﬀness but low er a verage p eak force. Disorder also enhances the softening b ehavior as disordered net works fails in a more gradual manner compared to uniform net works. 7. F or netw orks with the same p olymer chain concentration, increasing the co ordination num b er leads to stiﬀer netw orks with higher peak force, lo wer stretch at failure, and more gradual softening b eha vior. 8. The p eak force has a non-monotonic dep endance on the cross-linker concentration where it increases ﬁrst as the cross-link er concen tration increases but then decreases if the concen- tration exceeds a threshold v alue. Netw orks with longer c hains hav e higher threshold v alue than those with shorter ones. A c kno wledgmen t This w ork is supp orted b y the Center for Geologic Storage of CO2, an Energy F rontier Research Cen ter funded by the U.S. Departmen t of Energy (DOE), Oﬃce of Science, Basic Energy Sciences (BES), under A ward No. DE-SC0C12504, and the National Science F oundation CAREER A ward (Gran t No. 1753249). 24 App endix A: Eﬀect of disorder on the force stretch relation W e p erform a simple statistical analysis to determine the minimum num b er of samples enough to dra w some representativ e conc lusions ab out the eﬀect of netw ork disorder on the force-stretch relation. W e start with the highest disorder parameter ( µ =0.50), and w e determine the num b er of sim ulations based on the saturation of the mean and the standard deviation of the p eak force dis- tribution where the rate of change of b oth the mean the standard deviation after adding additional runs falls below 1%. F or the sp eciﬁc parameters we use, saturation is achiev ed after around 60 random netw orks. W e used this n umber for all other disorder parameters for consistency although lo wer n umber is suﬃcient in these cases. Fig. A.1 a-c show the normalized force vs stretch curves for the 60 runs along with the a verage normalized force-stretc h curv e for each disorder parameter. Fig. A.1 d,e sho w the ev olution of the mean and standard deviation with the num b er of runs. The a verage plots are shown in the man uscript. Figure A.1: The results from the statistical analysis of disorder net works: (a-c) Normalized force vs, stretch curves for all simulated netw orks with disorder parameters of µ = 0 . 05 , 0 . 25 , and 0 . 5 resp ectiv ely . The a v erage curv es are also plotted for each case, and (d,e) The evolution of the mean and the standard deviation of the peak force with the n umber of sim ulated netw orks for the three diﬀeren t studied disorder parameters. 25 References [1] C. Creton, M. Ciccotti, F racture and adhesion of soft materials: a review, Rep orts on Progress in Physics 79 (2016) 046601. [2] R. Ev eraers, K. Kremer, T op ology eﬀects in model p olymer netw orks, AIP Conference Pro ceedings 615 (2009) 615–626. [3] K. K othari, Y. Hu, S. Gupta, A. Elbanna, Mechanical Resp onse of T wo-Dimensional P olymer Net works: Role of T op ology , Rate Dep endence, and Damage Accum ulation, Journal of Applied Mec hanics 85 (2018) 031008. [4] A. N. Gent, A New Constitutiv e Relation for Rubb er, Rubber Chemistry and T echnology 69 (1996) 59–61. [5] R. S. Rivlin, Large Elastic Deformations of Isotropic Materials. IV. F urther Dev elopments of the General Theory, Philosophical T ransactions of the Roy al Society A: Mathematical, Ph ysical and Engineering Sciences 241 (1948) 379–397. [6] E. M. Arruda, M. C. Boyce, A three-dimensional constitutiv e mo del for the large stretch b eha vior of rubber elastic materials, Journal of the Mec hanics and Physics of Solids 41 (1993) 389–412. [7] Y. Hu, X. Zhao, J. J. Vlassak, Z. Suo, Using indentation to c haracterize the p oroelasticity of gels, Applied Physics Letters 96 (2010). [8] S. A. Chester, L. Anand, A thermo-mec hanically coupled theory for ﬂuid p ermeation in elastomeric materials: Application to thermally resp onsiv e gels, Journal of the Mechanics and Ph ysics of Solids 59 (2011) 1978–2006. [9] Y. Mao, L. Anand, F racture of Elastomeric Materials by Crosslink F ailure, Journal of Applied Mec hanics 85 (2018) 081008. [10] C. Miehe, L. M. Schänzel, Phase ﬁeld modeling of fracture in rubb ery polymers. Part I: Finite elasticit y coupled with brittle failure, Journal of the Mechanics and Physics of Solids 65 (2014) 93–113. [11] A. Kumar, G. A. F rancfort, O. Lop ez-P amies, F racture and healing of elastomers: A phase- transition theory and numerical implementation, Journal of the Mechanics and Physics of Solids 112 (2018) 523–551. [12] R. Spatschek, E. Brener, A. Karma, Phase ﬁeld mo deling of crack propagation, Philosophical Magazine 91 (2011) 75–95. [13] A. A. Griﬃth, The Phenomena of Rupture and Flo w in Solids, Philosophical T ransactions of the Roy al So ciet y A: Mathematical, Physical and Engineering Sciences 221 (1921) 163–198. [14] J. Zhang, S. O. P oulsen, J. W. Gibbs, P . W. V o orhees, H. F. Poulsen, Determining material parameters using phase-ﬁeld sim ulations and exp eriments, A cta Materialia 129 (2017) 229– 238. [15] M. J. Buehler, Atomistic Modeling of Mechanical F ailure, 2008. [16] C. P . Bro edersz, F. C. Mac kintosh, Mo deling semiﬂexible polymer net works, Reviews of Mo dern Physics 86 (2014) 995–1036. [17] T. Sak ai, T. Matsunaga, Y. Y amamoto, C. Ito, R. Y oshida, S. Suzuki, N. Sasaki, M. Shiba yama, U. I. Ch ung, Design and fabrication of a high-strength hydrogel with ide- ally homogeneous net work structure from tetrahedron-lik e macromonomers, Macromolecules 41 (2008) 5379–5384. [18] T. Sak ai, Y. Ak agi, T. Matsunaga, M. Kurak azu, U. I. Chung, M. Shibay ama, Highly elastic and deformable h ydrogel formed from tetra-arm p olymers, Macromolecular Rapid Comm uni- cations 31 (2010) 1954–1959. 26 [19] Y. Ak agi, T. Katashima, H. Sakurai, U. I. Ch ung, T. Sak ai, Ultimate elongation of p olymer gels with controlled netw ork structure, RSC A dv ances 3 (2013) 13251–13258. [20] G. Alamé, L. Brassart, Relative contributions of c hain densit y and top ology to the elasticity of tw o-dimensional p olymer netw orks, Soft Matter 15 (2019) 5703–5713. [21] H. Kamata, Y. Ak agi, Y. Kay asuga-Kariy a, U. I. Chung, T. Sak ai, "Nonswellable" hydrogel without mechanical hysteresis, Science 343 (2014) 873–875. [22] E. B. T admor, M. Ortiz, R. Phillips, Quasicontin uum analysis of defects in solids, Philosoph- ical Magazine A: Physics of Condensed Matter, Structure, Defects and Mechanical Prop erties 73 (1996) 1529–1563. [23] R. E. Miller, E. B. T admor, R. E. Miller, The Quasicon tinuum Metho d: Ov erview, applications and current directions, Journal of Computer-Aided Materials Design 9 (2002) 203–239. [24] V. B. Sheno y , R. Miller, E. B. T admor, D. Ro dney , R. Phillips, M. Ortiz, An adaptiv e ﬁnite elemen t approac h to atomic-scale mechanics - The quasicontin uum metho d, Journal of the Mec hanics and Ph ysics of Solids 47 (1999) 611–642. [25] E. B. T admor, R. Miller, R. Phillips, M. Ortiz, Nanoindentation and incipient plasticit y, Journal of Materials Research 14 (1999) 2233–2250. [26] J. Knap, M. Ortiz, An analysis of the quasicontin uum metho d, Journal of the Mechanics and Ph ysics of Solids 49 (2001) 1899–1923. [27] R. Miller, E. B. T admor, R. Phillips, M. Ortiz, Quasicontin uum sim ulation of fracture at the atomic scale, Mo delling and Simulation in Materials Science and Engineering 6 (1998) 607–638. [28] L. A. A. Beex, R. H. J. P eerlings, M. G. D. Geers, A quasicontin uum methodology for m ultiscale analyses of discrete microstructural mo dels, International Journal for Numerical Metho ds in Engineering 87 (2011) 701–718. [29] L. A. Beex, R. H. P eerlings, M. G. Geers, A multiscale quasicontin uum metho d for dissipative lattice mo dels and discrete net works, Journal of the Mechanics and Ph ysics of Solids 64 (2014) 154–169. [30] L. A. Beex, R. H. P eerlings, M. G. Geers, A multiscale quasicon tinuum metho d for lattice mo dels with b ond failure and ﬁb er sliding, Computer Metho ds in Applied Mec hanics and Engineering 269 (2014) 108–122. [31] L. A. Beex, R. H. Peerlings, K. V an Os, M. G. Geers, The mechanical reliability of an elec- tronic textile in vestigated using the virtual-pow er-based quasicon tinuum metho d, Mechanics of Materials 80 (2015) 52–66. [32] L. A. Beex, O. Rokoš, J. Zeman, S. P . Bordas, Higher-order quasicontin uum metho ds for elastic and dissipative lattice mo dels: Uniaxial deformation and pure bending, GAMM Mitteilungen 38 (2015) 344–368. [33] L. A. Beex, P . Kerfriden, T. Rabczuk, S. P . Bordas, Quasicontin uum-based multiscale ap- proac hes for plate-lik e b eam lattices exp eriencing in-plane and out-of-plane deformation, Com- puter Metho ds in Applied Mec hanics and Engineering 279 (2014) 348–378. [34] O. Rok oš, L. A. Beex, J. Zeman, R. H. P eerlings, A v ariational form ulation of dissipativ e quasicon tinuum metho ds, International Journal of Solids and Structures 102-103 (2016) 214– 229. [35] O. Rok oš, R. H. Peerlings, J. Zeman, L. A. Beex, An adaptive v ariational Quasicontin uum metho dology for lattice net works with localized damage, In ternational Journal for Numerical Metho ds in Engineering 112 (2017) 174–200. [36] K. Mikeš, M. Jirásek, Quasicontin uum metho d extended to irregular lattices, Computers and Structures 192 (2017) 50–70. 27 [37] K. Mikeš, M. Jirásek, Quasicon tinuum sim ulation of nanotextile based on the microplane mo del, Key Engineering Materials 714 (2016) 143–147. [38] K. Mikeš, M. Jirásek, Quasicon tinuum Metho d Combined with Anisotropic Microplane Mo del, A dv anced Materials Research 1144 (2017) 142–147. [39] G. P . Phlipot, D. M. Kochmann, A quasicon tinuum theory for the nonlinear mechanical resp onse of general p eriodic truss lattices, Journal of the Mec hanics and Ph ysics of Solids 124 (2018) 758–780. [40] C. W ellmann, P . W riggers, A t wo-scale mo del of granular materials, Computer Metho ds in Applied Mechanics and Engineering 205-208 (2012) 46–58. [41] Y. Y ue, B. Smith, P . Y. Chen, M. Chanthara yukhonthorn, K. Kamrin, E. Grinspun, Hybrid grains: Adaptiv e coupling of discrete and contin uum simulations of granular media, SIG- GRAPH Asia 2018 T echnical Papers, SIGGRAPH Asia 2018 37 (2018). [42] R. Hill, Elastic properties of reinforced solids: Some theoretical principles, Journal of the Mec hanics and Ph ysics of Solids 11 (1963) 357–372. [43] P . M. Suquet, Lo cal and global asp ects in the mathematical theory of plasticity ., In Plasticity T o day: Modelling, Metho ds and Applications (1985) 279–310. [44] M. G. D. Geers, V. G. Kouznetso v a, K. Matouš, J. Y v onnet, Homogenization Metho ds and Multiscale Mo deling: Nonlinear Problems, in: Encyclop edia of Computational Mec hanics Second Edition, John Wiley & Sons, Ltd, Chic hester, UK, 2017, pp. 1–34. [45] J. F. Marko, E. D. Siggia, Stretching DNA, Macromolecules 28 (1995) 8759–8770. [46] J. F. Marko, DNA under high tension: Ov erstretching, undertwisting, and relaxation dynam- ics, Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary T opics 57 (1998) 2134–2149. [47] C. K. Lieou, A. E. Elbanna, J. M. Carlson, Sacriﬁcial b onds and hidden length in biomaterials: A kinetic constitutive description of strength and toughness in b one, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 88 (2013) 1–10. [48] O. Rok oš, R. H. Peerlings, J. Zeman, eXtended v ariational quasicon tin uum methodology for lattice net works with damage and crac k propagation, Computer Methods in Applied Mec hanics and Engineering 320 (2017) 769–792. [49] M. C. Riv ara, New longest-edge algorithms for the reﬁnement and/or improv emen t of un- structured triangulations, International Journal for Numerical Metho ds in Engineering 40 (1997) 3313–3324. [50] M. T. Heath, SCIENTIFIC COMPUTING An Introductory Surv ey , Philosophical transac- tions. Series A, Mathematical, physical, and engineering sciences 363 (1997) 1707–13. [51] T. L. Anderson, F racture Mechanics: F undamentals and Applications, F ourth Edition, CRC Press, 2017. [52] S. Lee, M. Pharr, Sidew ays and stable crack propagation in a silicone elastomer, Pro ceedings of the National Academ y of Sciences (2019) 201820424. [53] D. Bonn, H. Kellay , K. Ben-djemiaa, J. Meunier, Delay ed F racture of an Inhomogeneous Soft Solid, Science 280 (2014) 265–267. [54] R. Ghelic hi, K. Kamrin, Mo deling growth paths of interacting crack pairs in elastic media, Soft Matter 11 (2015) 7995–8012. [55] M. L. F ender, F. Lechenault, K. E. Daniels, Universal shap es formed by tw o interacting cracks, Ph ysical Review Letters 105 (2010) 1–4. [56] Q. W ang, Z. Gao, A constitutive mo del of nano composite hydrogels with nanoparticle crosslink ers, Journal of the Mec hanics and Physics of Solids 94 (2016) 127–147. 28 [57] M. C. Boyce, E. M. Arruda, Constitutive Mo dels of Rubber Elasticity: A Review, Rubb er Chemistry and T ec hnology 73 (2000) 504–523. [58] Y. Qi, J. Caillard, R. Long, F racture toughness of soft materials with rate-independent h ysteresis, Journal of the Mec hanics and Physics of Solids 118 (2018) 341–364. [59] X. Zhao, A theory for large deformation and damage of interpenetrating p olymer netw orks, Journal of the Mechanics and Ph ysics of Solids 60 (2012) 319–332. [60] T. D. Nguyen, S. Govindjee, P . A. Klein, H. Gao, A material force metho d for inelastic fracture mechanics, Journal of the Mec hanics and Physics of Solids 53 (2005) 91–121. [61] J. Guo, M. Liu, A. T. Zehnder, J. Zhao, T. Narita, C. Creton, C. Y. Hui, F racture mechanics of a self-healing hydrogel with cov alent and physical crosslinks: A numerical study, Journal of the Mechanics and Ph ysics of Solids 120 (2018) 79–95. [62] K. Y u, A. Xin, Q. W ang, Mec hanics of self-healing p olymer net works crosslinked by dynamic b onds, Journal of the Mechanics and Physics of Solids 121 (2018) 409–431. [63] B. T alamini, Y. Mao, L. Anand, Progressive damage and rupture in p olymers, Journal of the Mec hanics and Ph ysics of Solids 111 (2018) 434–457. [64] Y. Mao, L. Anand, A theory for fracture of polymeric gels, Journal of the Mec hanics and Ph ysics of Solids 115 (2018) 30–53. [65] X. Poulain, V. Lefèvre, O. Lop ez-Pamies, K. Ravi-Chandar, Damage in elastomers: nucleation and growth of ca vities, micro-crac ks, and macro-cracks, International Journal of F racture 205 (2017). [66] S. R. Lav oie, P . Millereau, C. Creton, R. Long, T. T ang, A contin uum mo del for progressive damage in tough m ultinetw ork elastomers, Journal of the Mechanics and Physics of Solids 125 (2019) 523–549. [67] B. Zhang, C. S. Shiang, S. J. Y ang, S. B. Hutchens, Y-Shap ed Cutting for the Systematic Characterization of Cutting and T earing, Exp erimen tal Mechanics (2019). [68] G. Pijaudier-Cab ot, Z. P . Bazant, Nonlo cal damage theory, Journal of Engineering Mec hanics 113 (1988) 1512–1533. [69] E. Ducrot, Y. Chen, M. Bulters, R. P . Sijbesma, C. Creton, T oughening elastomers with sacriﬁcial b onds and w atching them break, Science 344 (2014) 186–189. [70] V. Lefèvre, K. Ravi-Chandar, O. Lop ez-P amies, Ca vitation in rubb er: an elastic instability or a fracture phenomenon?, International Journal of F racture 192 (2015) 1–23. [71] P . C. C. C. Y. Millereau, Detection of Molecular F racture in Elastomers with Mechanophores, in: Bulletin of the American Ph ysical So ciet y , p. 1. [72] G. Cusatis, D. Pelessone, A. Mencarelli, Lattice Discrete Particle Mo del (LDPM) for failure b eha vior of concrete. I: Theory, Cement and Concrete Comp osites 33 (2011) 881–890. [73] P . Khodabakhshi, J. N. Reddy , A. Sriniv asa, GraFEA: a graph-based ﬁnite elemen t approac h for the study of damage and fracture in brittle materials, Meccanica 51 (2016) 3129–3147. [74] Z. Zhang, Discretized virtual in ternal b ond model for nonlinear elasticit y, In ternational Journal of Solids and Structures 50 (2013) 3618–3625. [75] Z. P . Bazan t, Can Multiscale-Multiph ysics Methods Predict Softening Damage and Structural F ailure?, In ternational Journal for Multiscale Computational Engineering 8 (2010) 61–67. [76] K. Ravi-Chandar, W. G. Knauss, An exp erimental inv estigation into dynamic fracture: I I I. On steady-state crack propagation and crac k branc hing, In ternational Journal of F racture 26 (1984) 141–154. [77] E. K. Carlson, W atch tiny crac ks trav el in 3-D, EOS 99 (2018). 29 [78] W. Steinhardt, S. Rubinstein, High Sp eed Strain Measurements Surrounding Hydraulic F rac- ture in Brittle Hydrogel, in: APS Division of Fluid Dynamics Meeting Abstracts, p. A1.004. [79] M. M. Driscoll, B. G.-g. Chen, T. H. Beuman, S. Ulrich, S. R. Nagel, V. Vitelli, The role of rigidit y in controlling material failure, Proceedings of the National Academ y of Sciences 113 (2016) 10813–10817. [80] M. Astorino, F. Chouly , A. Quarteroni, A time-parallel framew ork for coupling ﬁnite element and lattice b oltzmann methods, Applied Mathematics Researc h eXpress 2016 (2016) 24–67. [81] 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. [82] T. D. Nguyen, C. M. Y ak ac ki, P . D. Brahmbhatt, M. L. Chambers, Mo deling the relaxation mec hanisms of amorphous shap e memory p olymers, Adv anced Materials 22 (2010) 3411–3423. [83] A. Kumar, O. Lop ez-P amies, On the tw o-p oten tial constitutive mo deling of rubb er visco elastic materials, Comptes Rendus - Mecanique 344 (2016) 102–112. [84] J. Zhang, M. Osto ja-Starzewski, Mesoscale b ounds in visco elasticit y of random comp osites, Mec hanics Research Comm unications 68 (2015) 98–104. [85] Y. No da, Y. Ha yashi, K. Ito, F rom top ological gels to slide-ring materials, Journal of Applied P olymer Science 131 (2014) 1–9. 30

An Adaptive Quasi-Continuum Approach for Modeling Fracture in Networked Materials: Application to Modeling of Polymer Networks

Original Paper

Comments & Academic Discussion

Leave a Comment

Original Paper

Related Papers

Comments & Academic Discussion

Leave a Comment