Laurent Brochard
email: laurent.brochard@enpc.fr

György Hantal

Hadrien Laubie

Franz-Joseph Ulm

Roland J M Pellenq

Capturing material toughness by molecular simulation Accounting for large yielding effects and limits

Keywords: Molecular simulation, LEFM, large yielding, cohesive zone

The inherent computational cost of molecular simulations limits their use to the study of nanometric systems with potentially strong size effects. In the case of fracture mechanics, size effects due to yielding at the crack tip can affect strongly the mechanical response of small systems. In this paper we consider two examples: a silica crystal for which yielding is limited to a few atoms at the crack tip, and a nanoporous polymer for which the process zone is about one order of magnitude larger. We perform molecular simulations of fracture of those materials and investigate in particular the system and crack size effects. The simulated systems are periodic with an initial crack. Quasi-static loading is achieved by increasing the system size in the direction orthogonal to the crack while maintaining a constant temperature. As expected, the behaviors of the two materials are significantly different. We show

Introduction

A common issue in mechanical engineering is the material fracture failure, which has to be accounted for carefully for the safety of structures. It is well-known in fracture mechanics that the presence of cracks in a material induces stress concentrations that greatly affect the maximum loading capacity of the material. According to Linear Elastic Fracture Mechanics (LEFM), the stress field is singular at crack tips and its asymptotic form when approaching a tip is characterized by the stress intensity factor. LEFM predicts that a crack propagates when the stress intensity factor exceeds a critical value, called toughness. LEFM was very successful at predicting the failure of mechanical structures, even though the physical existence of infinite stresses at the tip is of course questionable and one expects yielding of the material close to the tip. Actually, as long as the structure considered is much larger than the process zone (small scale yielding assumption), the approximation of LEFM is reasonably valid and this theory predicts correctly the mechanical behavior [START_REF] Anderson | Fracture mechanics: fundamentals and applications[END_REF]. However, when the size of the process zone is no more negligible compared to the size of the structure, LEFM fails to predict the structure's failure. Theories were developed to capture the effect of yielding at the crack tip, e.g., HRR singularity [START_REF] Hutchinson | Singular behaviour at the end of a tensile crack in a hardening material[END_REF] or cohesive zone models [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF]Barenblatt, 1962), in which the criterion of crack propagation is no more the toughness (e.g., critical J-integral or the crack tip opening displacement). Those theories are equivalent to usual LEFM under the assumption of small scale yielding, and all the failure criteria are thus related.

The size of the process zone varies widely depending on the material, from Angströms for cleavage fracture of monocrystals to centimeters for concrete and rocks [START_REF] Bazant | Size effect in blunt fracture: concrete, rock, metal[END_REF]. Any measurement of the toughness of a sample for which the process zone size is not negligible should account for yielding effects. This is particularly true when considering molecular simulations approaches since the studied system sizes barely exceed a few tens of nanometers. More fundamentally, a question that arises at the nanoscale is whether continuum fracture mechanics still applies at this scale. For instance, in their early work, [START_REF] Thomson | Lattice trapping of fracture cracks[END_REF] have shown that the discrete nature of matter at the atomic scale leads to a lattice trapping of the crack, an effect unknown in continuum fracture mechanics. Likewise, a great deal of research efforts was devoted to various atomic scale peculiarities of fracture mechanics: competition between crack propagation and dislocation emissions [START_REF] Rice | Ductile versus brittle behaviour of crystals[END_REF][START_REF] Decelis | Molecular dynamics simulation of crack tip processes in alpha-iron and copper[END_REF][START_REF] Cheung | A molecular-dynamics simulation of crack-tip extension: the brittle-to-ductile transition[END_REF], role of inter-atomic potentials [START_REF] Sinclair | The influence of the interatomic force law and of kinks on the propagation of brittle cracks[END_REF][START_REF] Holian | Fracture simulations using large-scale molecular dynamics[END_REF][START_REF] Marder | Effects of atoms on brittle fracture[END_REF][START_REF] Buehler | Dynamical fracture instabilities due to local hyperelasticity at crack tips[END_REF], role of phonons [START_REF] Holian | Fracture simulations using large-scale molecular dynamics[END_REF][START_REF] Zhou | Dynamic crack processes via molecular dynamics[END_REF][START_REF] Gumbsch | Molecular dynamics investigation of dynamic crack stability[END_REF], crack velocity [START_REF] Marder | Origin of crack tip instabilities[END_REF][START_REF] Buehler | Dynamical fracture instabilities due to local hyperelasticity at crack tips[END_REF], dynamic instability [START_REF] Marder | Origin of crack tip instabilities[END_REF][START_REF] Abraham | Large-scale simulations of brittle and ductile failure in fcc crystals[END_REF][START_REF] Buehler | Dynamical fracture instabilities due to local hyperelasticity at crack tips[END_REF][START_REF] Kermode | Low-speed fracture instabilities in a brittle crystal[END_REF], effects of crystal orientation and grain boundaries [START_REF] Miller | Quasicontinuum simulation of fracture at the atomic scale[END_REF][START_REF] Abraham | Large-scale simulations of brittle and ductile failure in fcc crystals[END_REF][START_REF] Pérez | Directional anisotropy in the cleavage fracture of silicon[END_REF], effect of chemical environment and impurities [START_REF] Lawn | Physics of fracture[END_REF][START_REF] Kermode | Macroscopic scattering of cracks initiated at single impurity atoms[END_REF]. A long standing issue in atomic scale studies is the size effect that may arise at such small scales. [START_REF] Holland | Cracks and atoms[END_REF] designed a molecular simulation approach and analysis for silicon that is scale insensitive, provided the simulated system is large enough (strip at least 80 atoms high). Bouch-binder et al (2010) discussed in detail various origins of scale effects that cause deviations from LEFM theory, such as non-linear elasticity and irreversibility at crack tip, potential energy corrugation at the atomic scale, energy flow near the tip. In all cases, LEFM limitations arise from the assumption of continuum linear elasticity at the crack tip. Apart from the intrinsic limitations of continuum theories, material properties are also size dependent when considering nanoscale systems. In particular, the high resistance to fracture of tough materials involves processes that range from breaking of atomic bonds to the nucleation of voids, crazing, crystal rotations or microcracking on the micron and larger scales [START_REF] Bazant | Size effect in blunt fracture: concrete, rock, metal[END_REF][START_REF] Garrison | Ductile fracture[END_REF][START_REF] Anderson | Fracture mechanics: fundamentals and applications[END_REF][START_REF] Ward | Infinite row of parallel cracks in a strip[END_REF]. The toughness of systems smaller than the range of those processes are size dependent. For brittle crystals with process zones a few atoms large, size effects appear for systems a few nanometer large (Nazmus [START_REF] Sakib | On the sizedependent critical stress intensity factor of confined brittle nanofilms[END_REF]; whereas size effects are at the scale of civil engineering structures for building materials [START_REF] Bazant | Size effect in blunt fracture: concrete, rock, metal[END_REF]. The toughness can also depend on the crack size because of surface energy changes with strain, lattice trapping, and non-linear mechanical behaviors [START_REF] Mattoni | Atomic scale origin of crack resistance in brittle fracture[END_REF][START_REF] Zhang | Atomistic and multiscale analyses of brittle fracture in crystal lattices[END_REF]. That is, the scale at which the effect of crack size appears is very dependent on the nature of the material.

Although molecular simulation approaches can capture the atomic details of fracture mechanisms, these simulations do not exceed a few billions of atoms with supercomputers and highly efficient implementation [START_REF] Abraham | Simulating materials failure by using up to one billion atoms and the world's fastest computer: Work-hardening[END_REF], which amounts to structures 0.1µm large. Accordingly, effects due to small sizes are a recurring issue when considering molecular simulation. For fracture studies, the assumption of small scale yielding often fails which is responsible for discrepancies between molecular simulation results and LEFM. In this work, we present an approach to analyze molecular simulation results in order to capture the fracture toughness of materials while accounting for large yielding. This approach is based on cohesive zone models, developed originally by [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] and Barenblatt (1962) to model plastic effect and suppress unphysical stress singularity at crack tips. Cohesive zone models are not the only approach to circumvent the stress singularity of LEFM. Others exists such as the (equivalent) dislocation-based approach [START_REF] Bilby | The spread of plastic yield from a notch[END_REF], the non-local theory [START_REF] Kröner | Elasticity theory of materials with long range cohesive forces[END_REF] and some more recent theories [START_REF] Pugno | Quantized fracture mechanics[END_REF][START_REF] Oh | A theory of fracture based upon an extension of continuum mechanics to the nanoscale[END_REF]. Since cohesive zone models are known to be inadequate for the modeling of dynamic fracture because of their inherent mathematical formulation [START_REF] Langer | Critical examination of cohesive-zone models in the theory of dynamic fracture[END_REF], we limit ourselves to quasi-static evolution in this work and disregard crack dynamics. The DBCS analysis ('Dugdale-Bilby-Cottrell-Swinden') of [START_REF] Mattoni | Atomic scale origin of crack resistance in brittle fracture[END_REF] is relevant to the approach developed in this paper. While Mattoni et al limited their study to the effect of crack size and analyzed only the critical stress at failure, we consider in this work the system size effect and analyze the complete loading curve.

As case studies, two materials are considered for molecular simulations of fracture. The first one, a silica crystal, is a brittle material for which the process zone is expected to be a few atoms large. The second one, a molecular reconstruction of a saccharose coke, is a highly cross-linked nanoporous polymer which exhibits a significantly larger process zone. Accordingly, the yielding effects should be limited in the first material, but significant in the second, thus providing two complementary cases. For each material, we investigate the effect of system and crack sizes on the fracture behavior, e.g., the critical stress at failure. The molecular simulation results serve as a basis to challenge theories. We identify LEFM strong shortcomings as soon as yielding becomes significant, whereas cohesive zone models proves much more versatile and universal with respect to yielding effect.

In part 2 of this paper, we present the molecular simulations of the two materials considered and illustrate some typical results. In part 3, we propose an analysis that captures the yielding effects by mean of cohesive zone models, and apply it to the two case studies.

Molecular simulations

Simulation details

The design of a molecular simulation is key when considering the simulation of fracture. In particular, major fracture properties such as brittle to ductile transition, crack tip blunting, or crack velocity are very sensitive to small details of the simulation such as the inflection of the attractive part of the inter-atomic potential [START_REF] Belytschko | Atomistic simulations of nanotube fracture[END_REF], the way to handle emitted phonon and dislocations at boundaries [START_REF] Cheung | A molecular-dynamics simulation of crack-tip extension: the brittle-to-ductile transition[END_REF], and the loading technique [START_REF] Holian | Fracture simulations using large-scale molecular dynamics[END_REF]. With the restriction to quasi-static behaviors, we disregard crack dynamics, and, accordingly, precautions regarding phonon propagation or loading technique are not critical in this work. However, the accuracy of the inter-atomic potential, the choice of boundary conditions and the control of dissipated energy (temperature) are still important issues. Regarding the inter-atomic potential, we use the ReaxFF potential [START_REF] Van Duin | Reaxff: a reactive force field for hydrocarbons[END_REF]. ReaxFF is a reactive potential trained on experimental data or quantum calculations results. It has the ability to adapt to various atomic configuration, e.g., for carbon sp, sp2 or sp3 bonds. Therefore, it is able to predict the complex molecular rearrangements that occur in fracture simulations (bond breaking, bond formation, charge transfer etc.). The original version of ReaxFF was parameterized for the simulation of polymers and organic matter. In this work, we use a parameterization of ReaxFF extended to the simulation of silica [START_REF] Van Duin | Reaxffsio reactive force field for silicon and silicon oxide systems[END_REF][START_REF] Chenoweth | Simulations on the thermal decomposition of a poly (dimethylsiloxane) polymer using the reaxff reactive force field[END_REF] which is suitable for the two materials we consider in this work. Previous studies on the fracture of silicon have used ReaxFF successfully (Buehler et al, 2006[START_REF] Buehler | Threshold crack speed controls dynamical fracture of silicon single crystals[END_REF], and the potential formulation does not seem to exhibit spurious features that would strongly affect the fracture behavior such as the force peak in the REBO reactive potential [START_REF] Belytschko | Atomistic simulations of nanotube fracture[END_REF]. In our simulations, the fracture behavior simulated with reaxFF is consistent with the expected materials behaviors: brittle fracture but without cleavage for the silica [START_REF] Swiler | Molecular dynamics study of brittle fracture in silica glass and cristobalite[END_REF][START_REF] Swiler | Molecular dynamics study of brittle fracture in silica glass and cristobalite[END_REF], and chain elongation and scission of the nanoporous carbon [START_REF] Rottler | Fracture in glassy polymers: a molecular modeling perspective[END_REF].

In the literature, various boundary conditions are used to model fracture, from clamped boundaries (fixed displacement) [START_REF] Abraham | Large-scale simulations of brittle and ductile failure in fcc crystals[END_REF] to prescribed forces on boundary atoms (fixed stress often coupled with continuum mechanics) [START_REF] Decelis | Molecular dynamics simulation of crack tip processes in alpha-iron and copper[END_REF] and atoms insertions/removals to follow steady state dynamics [START_REF] Holland | Ideal brittle fracture of silicon studied with molecular dynamics[END_REF]. A recurring issue is that artificial manipulation of the atom dynamics at the boundary (either by imposing positions or by imposing forces coming from a 'smooth' continuum field) leads to unphysical molecular dynamics in the vicinity of the boundary. When simulating materials near 0K, which is most often the case in the literature, this issue is of little importance. But in the present work, we consider finite temperatures (a key parameter especially for the polymer fracture behavior) and large scale yielding potentially affected by the boundary. Accordingly, spurious effects due to artificial manipulation of boundary atoms are an important issue in this work.

To avoid it, we use periodic boundary conditions and control the periodic box size to impose a mechanical loading to the system (see Fig. 1). Periodic boundary conditions mimic an infinite body, do not require to artificially manipulate boundary atoms and avoid surface tensions that arise in non-periodic systems. However, the behavior of a periodic body is not always representative of that of an infinite body (e.g., discrete number of phonon modes). Regarding fracture simulation a recognized issue is the interactions between a crack and its periodic replicas: from a dynamical point of view, the vibrational waves emitted from a crack propagate to the replicas and alter the crack propagation; from a static point of view, the stress intensity at a crack tip is affected by the presence of periodic replicas of the crack. Since we disregard crack dynamics in this work, we are only concerned with the second issue. We properly address this issue in the analysis part.

As mentioned, we perform molecular simulations at room temperature (300 K). Temperature is associated with atom's velocities at the molecular scale, also referred to as 'thermal agitation', which is known to affect the fracture behavior. For instance, the velocity gap phenomenon vanishes upon temperature increase because thermal agitation, on the order of the atomic corrugation of the energy landscape, overcomes the lattice trapping effect [START_REF] Holland | Cracks and atoms[END_REF]. Likewise, temperature increase can induce brittle to ductile transition by favoring blunting, dislocation nucleation, or creep [START_REF] Cheung | A molecular-dynamics simulation of crack-tip extension: the brittle-to-ductile transition[END_REF]. In the present work, temperature is an essential parameter since it strongly affects the mechanics and failure properties of polymers [START_REF] Rottler | Fracture in glassy polymers: a molecular modeling perspective[END_REF][START_REF] Ward | Infinite row of parallel cracks in a strip[END_REF]. As an illustration, the strength of the nanoporous polymer we study is two times smaller at room temperature than at 0 K. We use a Nose-Hoover thermostat to impose a temperature of 300 K in our molecular simulation, with a damping parameter of 10 fs (time step of integration 0.1 fs) [START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF]. Since the energy dissipation during fracture generates heat, we systematically wait for thermal equilibration after each loading step before acquiring physical observables. Doing so our results can be considered quasi-static.

Simulation of silica

The first material we consider is α-cristobalite, a silica crystal. Previous molecular simulation studies have in- vestigated the failure behavior of this crystal [START_REF] Swiler | Molecular dynamics study of brittle fracture in silica glass and cristobalite[END_REF][START_REF] Swiler | Molecular dynamics study of brittle fracture in silica glass and cristobalite[END_REF]. While α-cristobalite is known to exhibit a brittle behavior, cracking does not occur along a specific cleavage plane. Instead, fracture surfaces exhibit fractal character which is typical of silicon oxides [START_REF] Nakano | Growth of pore interfaces and roughness of fracture surfaces in porous silica: Million particle molecular-dynamics simulations[END_REF]. We performed molecular simulations of fracture of α-cristobalite by initiating an elliptic crack (aspect ratio 1/5) in the (001) direction of the crystal and by increasing progressively the box size in the direction orthogonal to the cracks.

To prepare the initial atomic configurations, we considered the unit cell from Downs and Palmer (1994) (4.9717 × 4.9717 × 6.9223 Å3 ). We considered a total of 24 different situations with various system sizes and initial crack sizes, as listed in Table 1. We display in Figure 2 some of the systems simulated. The smallest system considered is 29.8 × 20.8 × 5.0 Å3 large (209 atoms), and the largest 99.4 × 83.1 × 19.9 Å3 large (11 324 atoms). The smallest crack considered is 7.5 Å long (18 atoms removed) and the largest 50 Å long (548 atoms removed).

To load the system we performed small increments of length of the system in the direction orthogonal to the crack. In a first series of simulations, we increased the length by 1% every 10 ps, i.e., an equivalent strain rate of 0.001 ps -1 . To avoid rate effects, this strain rate is two orders of magnitude smaller than the strain rate at which the failure of silica becomes strain rate dependent [START_REF] Swiler | Molecular dynamics study of brittle fracture in silica glass and cristobalite[END_REF]. The average tensile stress in the vertical direction (Σ) was computed every 1 fs after a period of relaxation of the system of 5 ps after each increment. Σ can be interpreted as the average stress on the top and bottom boundary of the periodic cell (Fig. 1) and was computed with the virial equation [START_REF] Allen | Computer simulation of liquids[END_REF]. In a second series of simulations, we optimized the loading procedure to reduce the computation time. We considered larger increments (2.5%) performed more frequently (every 2 ps), corresponding to a strain rate of 0.0125 ps -1 , still small enough to prevent any rate effects. Σ was computed every 1 fs after a period of equilibration of 1 ps after each increment. These values were chosen as a good compromise 22.90 † Loading procedure: increments of Ly of 1% every 10 ps and computation of Σ every 1 fs after a period of equilibration of 5 ps. ‡ Loading procedure: increments of Ly of 2.5% every 2 ps and computation of Σ every 1 fs after a period of equilibration of 1 ps. between computational cost and accuracy. All the simulations were performed with the LAMMPS software [START_REF] Plimpton | Fast parallel algorithms for shortrange molecular dynamics[END_REF] (http://lammps.sandia.gov).

In Table 1, we list the critical stress at the onset of failure for all studied systems. In Figure 3 the full loading curve for the three systems displayed in Figure 2 are shown. The chart represents the evolution of the average stress Σ as a function of the strain L y /L y0 -1 during the loading procedure up to complete failure (L y0 is the initial size of the system). The non-zero stress in the initial state is due to the surface tension induced by the initiated crack and to the reaxFF potential for which the ground state of the material is not exactly that of the initial atomic structure (unit cell from [START_REF] Downs | The pressure behavior of a-cristobalite[END_REF]). We also display on the same figure a few snapshots of the molecular structure for system H that illustrate the different steps of failure. As expected, the mechanical behavior of α-cristobalite is nearly linear elastic up to the onset of crack propagation. The material failure is brittle with unstable cracking through the whole periodic cell under displacement-controlled loading, leaving behind a rough crack surface. As expected, the critical stress and strain at the onset of failure depend on the crack and system size. In particular, smaller cracks lead to higher critical stresses. The highest critical stress of 25 GPa is obtained for the systems with the smallest crack only 6 atoms long (systems R to X). This value however is still much lower than the yield stress of αcristobalite (35 GPa, estimated from molecular simulations of bulk α-cristobalite). Therefore, the process zone remains small, typically the size of a few atoms.

Simulation of nanoporous polymer

The second material we consider is CS1000, a molecular reconstruction of a saccharose-based heat-activated carbon (CS) obtained by Hybrid Reverse Monte Carlo method [START_REF] Jain | Molecular modeling of porous carbons using the hybrid reverse monte carlo method[END_REF]. CS1000 is a glassy nanoporous material composed of a highly cross-linked polymeric network (Fig. 4). Fracture mechanisms of polymers are peculiar because of their underlying molecular structure made of chains: chain elongation in the process zone leads to the formation of fibrils bridging the crack faces separated by voids (crazing phenomenon) [START_REF] Kramer | Fundamental processes of craze growth and fracture[END_REF][START_REF] Ward | Infinite row of parallel cracks in a strip[END_REF]. The crazing phenomenon is a strong toughening mechanism: glassy polymers can exhibit fracture energy release rates orders of magnitude larger than their surface energy. Fibrils are typically a few nanometers large and can be studied by molecular simulations; see, for instance, [START_REF] Rottler | Cracks and crazes: on calculating the macroscopic fracture energy of glassy polymers from molecular simulations[END_REF]. In the particular case of CS1000, the high degree of cross-linking strongly limits crazing, since important cross-linking is known to increase the brittleness of polymers [START_REF] Sauer | Effect of molecular variables on crazing and fatigue of polymers[END_REF]. With a high degree of cross-linking, the linear chain structure of polymers is turned into a three-dimensional network which strongly limits the chain mobility, and, therefore, reduces the plastic behavior while increasing the stiffness. Thus, a highly cross-linked polymer exhibits lit-tle or no crazing and chain scission is the main failure mechanism. By molecular simulation of bulk CS1000, we estimated the yield strength of the flaw-less matrix to be σ Y S = 19 GPa (high value because of the absence of flaws). We anticipate a toughness of CS1000 in the range of toughness commonly observed for glassy solid polymers, i.e., from 0.5 to 5 MPa. √ m [START_REF] Ashby | Materials Selection in Mechanical Design (3rd Edition). Butterworth-Heinemann Barenblatt GI[END_REF]. Assuming a toughness of K Ic = 2 MPa.

√ m, we anticipate a process zone size of r p = 4 nm at the crack tip (following the Dugdale-Barenblatt estimate [START_REF] Anderson | Fracture mechanics: fundamentals and applications[END_REF]:

r p = π 8 K Ic σ Y S 2
). This is a small process zone with respect to polymers in general since process zones associated with crazing can be tens of micrometers large. While the CS1000 process zone is limited in size, it is still about one order of magnitude larger than in the case of silica, which should strongly affect the fracture behavior in our simulation with significant size effects.

An additional complexity of polymer mechanics is the dependency on both temperature and loading rate. When temperature is raised or when loading rate is decreased, amorphous polymers change from glass-like to rubber-like behavior [START_REF] Ward | Infinite row of parallel cracks in a strip[END_REF]. In the glassy state, atoms vibrate close to their equilibrium position, whereas in the rubbery state, polymer chains are quite flexible and can adopt a variety of conformations. As a consequence, rubbery polymers, dominated by entropy, are compliant and viscoelastic; while glassy polymers are stiff and elastic. The transition from one regime to another is quite sharp around the glass transition temperature (for a given loading rate). According to the time-temperature equivalence, the viscoelastic behaviour at a temperature is related to the behavior at another temperature by a change in the loading rate only. Following this equivalence, a single master curve fully characterizes the dependence of a polymer mechanical property on both temperature and loading rate. These concepts apply to polymer mechanical failure as well, and master curves for failure properties can be drawn [START_REF] Smith | Dependence of the ultimate properties of a gr-s rubber on strain rate and temperature[END_REF]. In the present case, at room temperature (300 K), CS1000 is well below its glass transition temperature since the heat activation of saccharose at 1000 K removed the most flexible part of the polymer. In addition, the loading rate we consider in this study is high (0.0125 ps -1 ) because of the time scale accessible by molecular simulation. Interestingly, at low temperatures / high loading rates, as is the case here, the time-temperature equivalence master curves are nearly constant, i.e., the properties are temperature and rate independent. Nevertheless, in this work, to avoid any mis-interpretation due to temperature and rate effects, we applied the same temperature and the same loading rate to all CS1000 systems. For one of the systems, we repeated the calculation with a loading rate about one order of magnitude smaller to investigate if significant rate effect could occur.

A series of 15 systems of various sizes constructed from the CS1000 model, with an elliptic crack of various length, but constant aspect ratio (1/5), was considered (detailed characteristics of the systems are listed in Table 2). Three of these systems are displayed in Figure 4. The systems were loaded by increasing their length orthogonal to the crack by 2.5% every 2 ps (strain rate of 0.0125 ps -1 ). The average tensile stress (Σ) in the direction orthogonal to the crack was computed every 1 fs after a period of equilibration of 1 ps after each increment. In addition, we repeated the molecular simulation for system C' but with a much smaller loading rate of 0.001 ps -1 (System P'). The systems were loaded until complete failure was achieved. We report in Table 2 the critical stress at the onset of failure for all the systems simulated, and display in Figure 5 the average tensile stress Σ as a function of the strain, L y /L y0 -1, for the three systems displayed in Figure 4. In addition, we include in Figure 5 a few snapshots of molecular configuration for system N' to illustrate the steps of CS1000 failure. Contrary to α-cristobalite, CS1000 failure is stable upon displacement-controlled loading. One can clearly observe significant yield at the crack tips. As for α-cristobalite, the critical stress and strain at the onset of failure depend on the crack and system size. That is, as expected, the smaller the crack the larger the critical stress. But, in the case of CS1000, the maximum critical stress (15.53 GPa, system A') is quite close to the yield strength (19 GPa). in fact, for this system, the yielding of the crack plane is almost complete since the uncracked area for this system represents 80% of the cross section in the crack plane.

3 Mechanical analysis

LEFM analysis

We first consider the strict application of LEFM to the present situation. LEFM assumes linear elasticity even 12.79 † Loading procedure: increments of Ly of 2.5% every 2 ps and computation of Σ every 1 fs after a period of equilibration of 1 ps. ‡ Loading procedure: increments of Ly of 1% every 10 ps and computation of Σ every 1 fs after a period of equilibration of 5 ps. at the crack tips, a consequence of which is the stress singularity at the tips (σ ∝ 1/ √ r, where r is the distance from the crack tip). The stress intensity K I in mode I loading, as is the case here, is defined as K I = lim x→0 + √ 2πxσ yy (y = 0) for a crack orthogonal to the y direction and for a crack tip located at (x = 0; y = 0). The stress intensity characterizes the singularity of the stress at the crack tip and depends on the geometry and loading conditions of the system. For instance, for a crack of size 2a in an infinite body subjected to a remote stress Σ orthogonal to the crack, the stress intensity is

K I = Σ √ πa.
Failure occurs when the stress intensity reaches the fracture toughness, K Ic . Therefore, for the finite crack in an infinite body, the critical stress at failure is a function of the crack size Σ cr = K Ic / √ πa. For the periodic geometry we consider in this work (Fig. 1), the situation is more complex, since the stress intensity at a crack tip is affected by the presence of the periodic replicas of the crack. By dimensional analysis, the dimensionless reduced stress intensity K I / (Σ √ πa) also depends on the dimensionless quantities characterizing the periodic geometry, i.e., L y /L x and 2a/L x :

K I Σ √ πa = C K I 2a L x , L y L x (1)
This specific periodic geometry has been studied within the framework of LEFM [START_REF] Ward | Infinite row of parallel cracks in a strip[END_REF][START_REF] Isida | Rectangular plates, strips and wide plates containing internal cracks under various boundary conditions[END_REF][START_REF] Karihaloo | Doubly periodic arrays of bridged cracks and short fibrereinforced cementitious composites[END_REF][START_REF] Karihaloo | On the solution of doubly periodic array of cracks[END_REF]. While there exists no analytical formulation of C K I (2a/L x , L y /L x ), one can quantify it numerically. In this work, we followed the numerical approach proposed by [START_REF] Karihaloo | Doubly periodic arrays of bridged cracks and short fibrereinforced cementitious composites[END_REF], which consists in a pseudo-traction method based on the solution of a single periodic array of collinear cracks [START_REF] Tada | The Stress Analysis of Cracks Handbook[END_REF]. The problem of an infinite body with a doubly periodic array of crack submitted to a remote loading Σ is first decomposed into a homogeneous crack-free problem submitted to the remote loading and a subsidiary problem in which the system is submitted to the loading on the crack faces instead. The subsidiary problem is further decomposed into an infinite number of identical sub-problems where a body with a single row of periodic collinear crack is submitted to pseudo-tractions σ p on the crack faces. For the second decomposition to be valid, the pseudo-tractions must satisfy the consistency equation:

σ p (x) -2 +∞ j=1 a 0 K σj (x, u) σ p (u) du = Σ (2)
where the kernel K σj (x, u) represents the stress induced by the pseudo-traction σ p (u) in the sub-problem with collinear periodic cracks at y = -jL y at the location of the pseudo-traction σ p (x) in the sub-problem with the collinear periodic cracks at y = 0. It is readily obtained from the analytic solution for an infinite body with periodic collinear cracks [START_REF] Tada | The Stress Analysis of Cracks Handbook[END_REF]: (3) [START_REF] Karihaloo | Doubly periodic arrays of bridged cracks and short fibrereinforced cementitious composites[END_REF] solved Equation 2by means of a Gauss-Legendre quadrature method in which pseudotractions σ p (x) are discretized at the integration points. Once the pseudo-tractions are known, one can easily access any quantity of interest. For instance, the stress intensity results from the pseudo tractions of a single sub-problem (only one sub-problem leads to a singular stress):

K σj (x, u) = Re (Z (x + i (jL y ) , u))
K I = a 0 K K I (u) σ p (u) du (4)
with

K K I (u) = 2 L x cos πu Lx L x tan πa Lx sin πa Lx 2 -sin πu Lx 2 (5)
Likewise, the deformation of the periodic cell is obtained by the superposition of the deformations in the homogeneous problem and in all the sub-problems. The displacement ∆ of the periodic boundary at y = L y /2 is:

∆ = L y 2 Σ E + +∞ j=-∞ a 0 K ∆j (u) σ p (u) du (6)
where

K ∆j (u) = (2/E ) Im Z i j + 1 2 L y , u -Z (i (jL y ) , u
) and Z = dZ/dz. E is the elastic modulus in plane strain, which is the case here: E = E/ 1 -ν 2 , with E the Young's modulus and ν the Poisson's ratio.

We applied the approach of [START_REF] Karihaloo | Doubly periodic arrays of bridged cracks and short fibrereinforced cementitious composites[END_REF] to account for the interactions between periodic cracks. The values of the correction factor C K I in Equation 1 are generally greater than 1, i.e., the influence of the periodic replicas of a crack leads to higher stress intensities than for an infinite body; up to 20% increase for the geometries we consider in this work. But, for some geometries, e.g. L y /L x < 0.6 and 2a/L x < 0.45, the stress intensity is moderately lower than in the infinite body (C K I < 1) which leads to a toughening of the material. Thus, the effect of the doubly periodic array of crack can be counter-intuitive and needs to be accounted for carefully.

Like the stress intensity, the deformation of the periodic cell is affected by the presence of periodic cracks. The ratio between the dimensionless deformation ∆ * = 2∆/L y and dimensionless stress Σ * = Σ/E is a function of the dimensionless quantities L y /L x and 2a/L x characterizing the periodic geometry, and can be computed from Equation 6:

2∆/L y Σ/E = C ∆ 2a L x , L y L x (7)
With these corrections, one can fully characterize the failure behavior as predicted by LEFM. Before failure, the mechanical behavior is linear elastic and follows Equation 7. Failure occurs when the stress intensity reaches the toughness. The critical stress at the onset of failure is

Σ cr = K Ic √ πaC K I 2a Lx , Ly Lx (8)
Then, in the quasi-static limit, the stress intensity is equal to the toughness all along the crack propagation. Therefore, for any crack length larger than the initial length a > a 0 , the stress can be obtained from Equation 8 and the corresponding strain from Equation 7. We summarize the loading and failure of the periodic system as follows:

Before failure, for a = a 0

∆ * K * Ic = C ∆ 2a0 Lx , Ly Lx Σ * K * Ic with Σ * K * Ic < C K I 2a0 Lx , Ly Lx -1
Then, for a > a 0

∆ * K * Ic = C ∆ 2a Lx , Ly Lx Σ * K * Ic with Σ * K * Ic = a a0 C K I 2a Lx , Ly Lx -1 (9)
where we introduced the dimensionless reduced strain ∆ * /K * Ic and the dimensionless reduced stress

Σ * /K * Ic with K * Ic = K Ic / E √ πa 0 .
Note that after failure the control variable is the crack length a that increases from a 0 to L x /2. Figure 6 displays the mechanical behavior of the periodic system as predicted by LEFM. It highlights that the behavior depends on the dimension ratios 2a/L x and L y /L x . The behavior is linear elastic until the onset of failure, with an elasticity depending on the geometric ratios. Then, failure is almost always unstable upon both stress-controlled loading and displacementcontrolled loading. That is, we almost always have Fig. 6 Loading curves as predicted by LEFM for various crack sizes and periodic cell elongations. We display the dimensionless reduced stress Σ * /K * Ic as a function of the dimensionless reduced strain ∆ * /K * Ic . The behavior is linear elastic until

K I = K Ic , i.e., Σ * /K * Ic = C K I a 0 L x , L y L x
. Then, cracking occurs until complete propagation through the periodic box. 0 and d∆ * da < 0. This is a peculiarity of the periodic system. In the periodic geometry, when the crack propagates, the tips come closer to their periodic replicas which enhances stress intensity and favors instability. Failure is stable upon displacement loading if ∂K I ∂a ∆ < 0. For the periodic system, this criterion is verified only for quite elongated systems (L y /L x small) and the range of crack length for which failure is stable is limited. The systems we considered are all close to a square geometry. Therefore, we expect unstable failure in our case.

The critical stress at failure (Eq. 8) also differs from the case of an infinite body, because of the correction term C K I . Equation 8 can be rewritten in the following dimensionless form:

Σ * cr G * c = π 2a L x C K I 2a L x , L y L x -1 (10)
where we introduced the dimensionless reduced failure stress

Σ * cr / G * c , with Σ * cr = Σ cr /E and G * c = G c / Lx 2 E . G c
is the critical energy release rate and is related to the toughness K Ic according to Irwin's formula [START_REF] Anderson | Fracture mechanics: fundamentals and applications[END_REF]:

G c = (K Ic )
2 /E . We display Fig. 7 LEFM prediction of the critical stress at the onset of failure. We display the dimensionless reduced failure stress Σ * cr / G * c as a function of the crack size 2a/L x for different elongations L y /L x of the periodic cell. We also display the case of a finite crack in an infinite body:

Σ * cr / G * c = 1/ π 2a L x
. When 2a/L x → 0 one recovers the critical stress of an infinite body. In contrast, significant deviations appears for larger crack sizes, especially for 2a/L x > 0.5.

in Figure 7 the dependence of Σ * cr / G * c on the geometric ratios 2a/L x and L y /L x . For small crack sizes (2a/L x → 0), we recover the critical stress of a finite crack in an infinite body. However, for moderate and large crack sizes (2a/L x > 0.2), the critical stress of a periodic system deviates from that of an infinite system. Particularly large deviations appear for 2a/L x > 0.5. In the limit of 2a/L x → 1 the crack tips merge at the periodic boundary, and, accordingly, the critical stress converges to 0. For the systems we simulated, the ratio 2a 0 /L x ranges from 0.13 to 0.60. Therefore, accounting for the correction C K I is necessary in our work to analyze the critical stresses.

Application of the LEFM analysis to molecular simulation results

We apply the LEFM analysis to the molecular simulation results of α-cristobalite and CS1000. Regarding loading curve shape and failure stability, it is clear that CS1000 results are in contradiction with LEFM prediction: CS1000 failure is stable upon displacement loading, and the loading curves (Fig. 5) exhibit significant ductility with an ultimate strain at the end of failure about twice the critical strain at the onset of failure. In contrast, the α-cristobalite results (Fig. 3) are consistent with LEFM: the behavior is almost linear and failure is unstable since stresses fall to zero right after the onset of failure. In our molecular simulations, it is not possible to follow the instable branches that appear in Figure 6. Instead, the stress drops to zero at failure.

A quantitative comparison with LEFM theory can be achieved by plotting the critical stress at failure as a function of the effective crack half length

a ef f = a C K I 2a Lx ,

Ly Lx

2

. If LEFM theory was respected, following Equation 8, we would observe that the critical stress is proportional to the inverse of the square root of the effective crack length: Σ cr = K Ic / √ πa ef f . By fitting the curve, one can recover the fracture toughness K Ic characterizing the material. We display in Figure 8 the critical stresses obtained for molecular simulations as a function of a ef f for all simulated systems. On the same Figure we display the fitted inverse square root dependency according to LEFM theory. For αcristobalite, the results seem consistent with LEFM.

In the Figure, we distinguish between the two loading rates (see Table 1). We do not observe any rate effects, but it appears that the long simulation procedure (smallest rate) provides more accurate results than the short procedure (largest rate) for which the results are significantly more dispersed. Despite the dispersion, the results are quite consistent with LEFM theory. The fitting leads to a fracture toughness K Ic = 0.87 MPa. √ m. This value is consistent with experimental values of toughness of silica. Indeed, [START_REF] Lucas | Determining fracture toughness of vitreous silica glass[END_REF] reported K Ic = 0.82 ± 0.07 MPa.m 1/2 . In contrast to α-cristobalite, the results of CS1000 exhibit a rather linear trend with the effective crack length, which is in contradiction with the inverse square root dependency predicted by LEFM. Any tentative application of LEFM to the results for CS1000 leads to serious discrepancies between theory and observation. A tentative application is displayed in Figure 8 with large discrepancies at low crack lengths; the corresponding toughness is the same as for α-cristobalite (K Ic = 0.87 MPa.

√ m).

As an alternative to critical stress, one can focus on the energy released during failure. The critical energy release rate G c is the energy released during the fracture per unit area of crack surface created. Fracture occurs when the system can release enough mechanical energy P upon cracking, with G c defining this threshold. Formally, G c is the variation of mechanical energy -dP due to a small advance of the crack (increment of crack area dA): G c = -∂P ∂A T,loading . Together, G c and P are the two forms of energy for a system exposed to mechanical loading whose response is limited to deformation or fracture. In isothermal conditions, as is the case here, the balance of Helmholtz free energy F is: dF = dP + G c dA. The average stress Σ is the derivative of the Helmholtz free energy F with respect to the periodic cell size L y : Σ = 1 LxLz ∂F ∂Ly T,Lx,Lz . Accordingly, integrating Σ with respect to L y provides a way to calculate the variation of free energy. In particular the integration of the complete loading curves (Fig. 6) leads to the critical energy release rate. Indeed, at the beginning of the loading curve, the system is unstressed, i.e., P = 0; at the end of the loading curve, the system is also unstressed and P = 0. Therefore, all the mechanical energy accumulated during the loading was converted into fracture energy at the end of the loading. Accordingly, the change of Helmholtz free energy calculated from the integration L x L z loading ΣdL y is equal to the total increase in fracture energy G c ∆A, where ∆A = L z (L x -2a 0 ) is the total area of crack created:

G c = L x L z ∆A loading ΣdL y (11)
The critical energy release rate G c , as obtained from Equation 11, can be converted into fracture toughness K Ic following Irwin's formula [START_REF] Anderson | Fracture mechanics: fundamentals and applications[END_REF]:

G c = (K Ic )
2 /E . We applied this thermodynamic integration to compute the toughness of α-cristobalite and CS1000 from the molecular simulation results. For αcristobalite, we obtained an average toughness of 1.13± 0.09 MPa. √ m, which is somewhat higher than the value of 0.87 MPa.

√ m derived from the critical stress. This discrepancy arises from the fact that integrating the curves of Figure 3 only provide an overestimation of the critical energy release rate because of instability. Indeed, proper integration should follow the unstable branch that recede to smaller strains when the crack advances (such as in Figure 6). But, our molecular simulations are unable to follow the unstable branch and instead the stress drops to zero. As a consequence, the integration of stress (area delimited by the loading curve) overestimates the critical energy release rates. For CS1000, the same restriction does not apply since failure is stable. The integration of the loading curves for CS1000 leads to an average toughness of 1.49 ± 0.14 MPa.

√ m, i.e., almost twice the value estimates from the critical stress and LEFM. With this toughness, the critical stresses predicted by LEFM would be twice that obtained by molecular simulations. This further confirms that LEFM is invalid for CS1000.

To illustrate the application of LEFM to α-cristobalite, we display in Figure 9 the loading curves from molecular simulations along with that derived from LEFM, for the three systems of Figure 2. LEFM captures satisfactorily the elastic behavior and predicts the onset of failure with a reasonable accuracy given the inherent variabilities of the simulations (in particular thermal agitation and variable crack tip initial configuration). The failure branches differ significantly as already discussed above. Figure 9 illustrates well the overestimation of the toughness when performing a thermodynamic integration on the molecular simulation results.

Cohesive zone analysis

LEFM does not capture the behavior of CS1000 in our molecular simulation. This is due to the important yielding at the scale of our molecular simulations. LEFM is valid only in the limit of small scale yielding, which is clearly not the case for the CS1000 systems we simulated. As a consequence, CS1000 response is ductile without instability at failure and the critical stresses obtained are much smaller than expected from LEFM because of the inherent physical limit of CS1000 yield stress. We propose here to introduce a cohesive zone approach in the mechanical modeling of the periodic system in order to capture the behavior of CS1000. Cohesive zones are a simple modeling of the yield zone at crack tips originally introduced by [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] and Barenblatt (1962) to suppress the unphysical stress singularities at crack tips. The approach is illustrated in Figure 10. Behind the crack tips, the crack faces are cohesive, i.e., an attractive force p between the faces opposes the opening of the crack. In the original formulation, the cohesive force is constant (plastic stress) and vanishes at a critical opening δ cr of the crack. More generally, the cohesive force p is a function of the opening δ, called the cohesive law. The negative stress intensity at the crack tip due to the cohesive forces cancels out the stress intensity from the external loading so that stress singularity is suppressed. This condition is used to determine the length of the cohesive zone. Failure occurs when the opening at the end of the cohesive zone exceeds the critical opening at which the cohesive force vanishes. Application of the J-integral to a contour following the boundary of the cohesive zone shows that the critical energy release rate is obtained by integration of the cohesive law [START_REF] Anderson | Fracture mechanics: fundamentals and applications[END_REF]:

G c = δcr 0 p (δ) dδ (12)
Thus, one can easily relate the cohesive zone model to usual LEFM. For small scale yielding, the cohesive zone is extremely small compared to the dimensions of the system, in particular the size of the crack. In Cohesive zones develop behind crack crack tips in which a force is transmitted across the crack, depending on the crack opening. The cohesive zone suppresses the stress singularity at the crack tip. At large openings, the cohesive forces vanish and the crack is unbridged. Loading proceeds in three phases: 1) the cohesive zone grows until the opening at the end of the zone reaches the critical opening, 2) the unbridged crack propagates, 3) the cohesive zones coalesce at the periodic boundaries. In the third phase, the mechanical problem is that of an infinite plate submitted to periodic traction loads on its faces.

that case, cohesive zone models and LEFM are strictly equivalent and the two approaches are related by Equation 12. When yielding becomes significant, cohesive zone predictions differ from that of LEFM. In particular, the stresses cannot exceed the yield stress, the main physical contradiction of LEFM.

Even though cohesive zone models are an idealization of the process zone, they are quite versatile and have been used to model very different failure processes from metal plasticity [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] to biocomposites damage [START_REF] Gao | Application of fracture mechanics concepts to hierarchical biomechanics of bone and bone-like materials[END_REF] and fracture in concrete [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF]. We propose here to apply it to the periodic geometry of our molecular simulations and investigate if it can capture the behavior of CS1000.

The cohesive zone modeling introduces bridging forces between the crack faces that depend on the opening of the crack. The loading of the system takes place as follows (see Fig. 10). Prior to loading, in the initial state, there is no cohesive zone. Then, at the beginning of the loading, the crack tip starts to move ahead and the cohesive zone develops. When the opening at the end of the cohesive zone exceeds the maximum opening of the cohesive law, the unbridged crack starts to move ahead as well. Finally, when the crack tip coalesce at the periodic boundary, the cohesive zones keep transmitting stress across the crack until the unbridged crack coalesce. For the first regime before crack tip coalescence, the mechanical problem can be solved by introducing bridging forces in the periodic crack problem treated in the LEFM approach. For the second regime after coalescence, the mechanical problem is that of an elastic plate submitted to symmetric periodic tractions on its faces. This second mechanical problem can be solved in the Fourier space. We present both solutions hereafter.

We first consider the system before crack tip coalescence. The cohesive zone introduces bridging forces p behind the crack tip as a function of the opening δ (x) at position x. These bridging forces modify the pseudotraction consistency Equation 2 as follows:

σ p (x) -2 +∞ j=1 a 0 K σj (x, u) σ p (u) du + p (δ (x)) = Σ (13)
The crack opening δ (x) can be computed from the pseudo-tractions. Since the displacement discontinuity that leads to the crack opening is due to one of the subproblems only, we can use of the analytic formulation for a row of periodic collinear cracks [START_REF] Tada | The Stress Analysis of Cracks Handbook[END_REF]. We have:

δ (x) = a 0 K δ (x, u) σ p (u) du (14)
with

K δ (x, u) =        8 πE tanh -1 1-(cos( πa Lx )/ cos( πu Lx )) 2 1-(cos( πa Lx )/ cos( πx Lx )) 2 if |x| ≤ u 8 πE coth -1 1-(cos( πa Lx )/ cos( πu Lx )) 2 1-(cos( πa Lx )/ cos( πx Lx )) 2 if u < |x| ≤ a (15)
Accordingly, the cohesive force distribution p (δ (x)) in Equation 13is a function of the pseudo-tractions σ p (x). Solving for the pseudo-tractions in Equation 13, one can then access any quantity of interest (displacement / deformation, stress intensity etc.). As before, we followed the approach of [START_REF] Karihaloo | Doubly periodic arrays of bridged cracks and short fibrereinforced cementitious composites[END_REF] (Gauss-Legendre quadrature method with σ p (x) discretized at the integration points). However, the problem is no more linear because of the cohesive law p (δ). We used an iterative Newton-Raphson method to solve the nonlinear problem. In addition, special attention must be paid to discretization with respect to the discontinuities of bridging forces to avoid integration and convergence errors. At each loading step, the length of the cohesive zone is adjusted so that the total stress intensity at the crack tip is zero.

After crack tip coalescence at the periodic boundaries, the mechanical problem becomes that of stacked periodic plates separated by cohesive zones which exert a load function of the distance between consecutive plates (Fig. 10). By symmetry, the stress and strain fields are even about the center of the unbridged crack. To solve this new mechanical problem, we consider the solution of [START_REF] Sneddon | Fourier Transforms[END_REF] for an infinite plate loaded by even forces on its faces. The solution of Sneddon consists in solving the mechanical problem in the Fourier space. Applying a cosine Fourier transform to the equations of mechanical equilibrium in the absence of body forces gives:

∂ 2 ∂y 2 -ξ 2 G = 0 (16) with G (ξ, y) = ∞ 0 χ (x, y) cos (ξx) dx
where χ is the Airy stress function. The function G replaces χ in the Fourier space and fully determines the stress, strain and displacement fields. In particular, the normal component of the stress in the vertical direction, the shear stress and the displacement in the vertical direction are:

σ y (x, y) = - 2 π ∞ 0 ξ 2 G cos (ξx) dξ (17)
τ xy (x, y) = 2 π ∞ 0 ξ ∂G ∂y sin (ξx) dξ (18) u y (x, y) = 2 π ∞ 0 1 + ν E 1 -ν ξ 2 ∂ 3 G ∂y 3 (19) -(2 -ν) ∂G ∂y cos (ξx) dξ
In the particular case of a plate, the solutions of Equation 16 take the form: Equation 22 gives the deformed configuration of a linear elastic plate under the action of an even symmetric periodic load on its faces. Owing to the cohesive law, the load p is a function of the opening δ of the crack, that is of the distance between two consecutive plates. The opening δ and the displacement of the boundaries u y (±b) are related according to:

G (ξ,
δ (x) = δ 0 -2 • u f aces y (x) -u f aces y (0)
, where δ 0 is a parameter representing the opening at x = 0. The mechanical problem can be formulated in terms of a consistency equation that must be satisfied by the opening δ:

δ (x) = δ 0 - 16 πE × ∞ k=0 sinh πk Ly Lx 2 2πk Ly Lx + sinh 2πk Ly Lx p k cos 2πk Lx x -1 2πk Lx (23)
where

p k = 2 Lx Lx 2 0 p (δ (x)) cos 2πk
Lx x dx and p (δ (x)) is the cohesive law. The opening function δ (x), which is the solution of Equation 23, fully characterizes the mechanics of the system since one can then derive from it the full stress, strain and displacement fields. Except for linear cohesive laws, the problem of Equation 23 is non linear and no simple analytic solution exists. In our work, we solve Equation 23 by adopting a similar approach as was done previously for computing the doubly periodic crack problem with cohesive zones. That is, we use Gauss-Legendre method for integration of the Fourier coefficients p k while the opening δ is discretized at the integration points. We solve Equation 23 iteratively with a Newton-Raphson method. Doing so, we can solve the mechanical problem for a wide variety of cohesive laws, but the discretization must be adapted carefully from one law to another to avoid integration and convergence errors (due to discontinuities in the law for instance).

The relationship between the dimensionless strain ∆ * = 2∆/L y and stress Σ * = Σ/E depends not only on the geometric ratios 2a/L x and L y /L x as in Equation 7, but also on the cohesive law. As an example, we display in Figure 11 the case of a plastic cohesive law, i.e., p (δ) = σ Y S for δ < δ cr with σ Y S the yield stress of the material. In Figure 11, we show how the loading curve changes with the dimensionless ratio δ * cr /σ * Y S = (2δ cr /L x ) / (σ Y S /E ) which quantifies the degree of ductility of the plastic cohesive law. Alternatively, we can also relate this quantity to the ratio between the critical stress for an infinite body (K * Ic ) and the yield stress (σ * Y S ) as follows:

δ * cr /σ * Y S = (2πa/L x )•(K * Ic /σ * Y S )
2 . Accordingly, we expect a significant difference with LEFM at large values of δ * cr /σ * Y S . This is well illustrated in Figure 11. For δ * cr /σ * Y S → 0, the plastic cohesive zone model predicts the same loading curve as LEFM. But with increasing ductility ratio, the cohesive zone model introduces non linearity and ultimately a plastic plateau, so that the shape of the loading curve significantly differ from the LEFM prediction.

In Figure 11, we also show how the ductility ratio δ * cr /σ * Y S affects the critical stress at failure. We display the dimensionless reduced failure stress Σ * cr / G * c = Σ * cr / δ * cr σ * Y S as a function of the crack size 2a/L x for L y /L x = 1. The cohesive zone approach suppresses the stress singularity, and, accordingly, at low 2a/L x , the critical stress converges to a finite value Σ cr = σ Y S , i.e., Σ * cr / G * c = 1/ δ * cr /σ * Y S . As for the loading curve, the predictions of the cohesive zone model converge to LEFM predictions for δ * cr /σ * Y S → 0. But, significant discrepancy appears with increasing ductility. In particular, the usual trend in 1/ √ a gradually evolves toward a linear trend. The limit case is when the cohesive zone extends over the whole uncracked body:

Σ cr σ Y S = 1 - 2a L x ⇒ Σ * cr G * c = 1 - 2a L x 1 δ * cr /σ * Y S (24)
3.4 Application of the cohesive zone analysis to molecular simulation results

We have shown in the previous section that cohesive zone models lead to significant changes in the loading curve and critical stress of the periodic system. The plastic cohesive zone illustrated in Figure 11 is an example, but it does not captures the shape of the loading curve we obtained for the nanoporous polymer (Fig. Fig. 11 Loading curves and critical as predicted with a plastic cohesive zone model for various ductility ratios δ * cr /σ * Y S . For the loading curves (top), we display the dimensionless reduced stress Σ * /K * Ic as a function of the dimensionless reduced strain ∆ * /K * Ic (with the geometrical ratios 2a/L x = 0.5 and L y /L x = 1). The loading curve from LEFM is displayed for comparison. The cohesive zone model introduces non linearity and the loading curves deviates from LEFM predictions. Deviation from LEFM increases with ductility and a plateau appears that is characteristic of the plastic cohesive law. For the critical stress (bottom), we display the dimensionless reduced failure stress Σ * cr / G * c as a function of the crack size 2a/L x (with the geometric ratio L y /L x = 1). For comparison, we also display the prediction of LEFM for the same periodic geometry and for a finite crack in an infinite body. The usual trend in 1/ √ a gradually evolves toward a linear trend with increasing ductility. 5). To capture the behavior of the nanoporous polymer, we considered a linear decreasing cohesive law instead: p (δ) = σ Y S (1 -δ/δ cr ) for δ < δ cr . We display in Figure 12 the loading curves and critical stress trend predicted with this cohesive law for various ductility ratios δ * cr /σ * Y S . The quantities displayed in Figure 12 are the same as in Figure 11. Note, however, that we consider a higher range of values for the ductility ratio δ * cr /σ * Y S since the linear cohesive law dissipates two times less energy than the plastic cohesive law (G * c = δ * cr σ * Y S /2). As a consequence, based on energy equivalence, the ductility ratios of Figure 12 are scaled by a factor of two compared to the ratios in Figure 11, whence the higher values. 11. The loading curves deviates from LEFM prediction with increasing ductility and exhibit a linear decreasing branch at large ductility characteristic of the linear cohesive law. The usual trend of the critical stress in 1/ √ a evolves toward a linear trend at high ductility.

In Figure 12, the qualitative shape of the loading curves at high ductility ratios are quite consistent with the molecular simulation results for CS1000 (Fig. 5). The trend of the critical stress is similar to that obtained with a plastic cohesive law: in the limit δ * cr /σ * Y S → 0 the critical stress converges to the LEFM prediction, while at large ductility, the critical stress depends linearly on the crack length, with a limit case:

Σ cr σ Y S = 1 - 2a L x ⇒ Σ * cr G * c = 1 - 2a L x 2 δ * cr /σ * Y S (25)
Equation 25 differs from Equation 24 by a factor of 2 which arises from the difference in energy release with the plastic cohesive law. We evaluated the trend of the critical stress from the results of molecular simulation of CS1000. We display in Figure 13 the critical stresses for the various CS1000 systems simulated. Figure 13 shows that the critical stress obtained by molecular simulation follow a linear trend in a limited range of crack length (0.2 < 2a 0 /L x < 0.6). Note that the systems considered in Figure 13 have different geometrical ratios L y /L x which is not accounted for in this analysis. Linear trend Systems A' to P' Bulk CS1000

13 Trend between the critical stress Σ cr and the reduced crack length 2a 0 /L x for all the CS1000 systems considered.

A major consequence is that the dimensionless ductility ratio δ * cr /σ * Y S differ from one system to another. But, in raw units (case of Fig. 13) the effect of L y /L x vanishes with increasing ductility ratio as the critical stress approaches the limit case (Eq. 25). The range of the linear trend (0.2 < 2a 0 /L x < 0.6) is not wide enough to characterize the degree of ductility of the cohesive law. Assuming that the linear relation applies (Eq. 25), the y-intercept of the linear trend in Figure 13 is the yield stress σ Y S . Doing so, we obtain σ Y S = 20 GPa, which compares well with the yield stress estimated by molecular simulation of bulk CS1000 (19.4 GPa). Accordingly, the linear trend is valid up to small crack length which is characteristic of a high ductility ratio (δ * cr /σ * Y S > 2). The analysis of the loading curve shape and critical stress confirms that the linear cohesive zone model is well adapted to model the behavior of CS1000.We calibrated the value of δ cr to capture the loading curves from molecular simulations. We found that with δ cr = 18.5 Å the cohesive zone model captures reasonably well the mechanical responses of all the CS1000 systems considered. As an illustration, we display in Figure 14 the loading curves for the three CS1000 systems of Figure 4 along with the predictions from the linear cohesive zone model. In the framework of the cohesive model, the energy dissipated during fracture is concentrated in the cohesive zone and the critical energy release rate G c can be obtained by integration of the cohesive law (Eq. 12). Thus, with the linear cohesive law, we have: G c = δ cr σ Y S /2 = 18.5 N/m. Alternatively, we can compute G c by integrating the loading curves from molecular simulation (Eq. 11). Unlike the case of α-cristobalite, this estimation is valid here since crack propagation is stable for CS1000. We obtain G c = 19.5 ± 3.2 N/m which is consistent with the value obtained from the cohesive law. In addition, no particular correlation is found between the energy re- lease rates and the system size. If the process zone size was larger than the system size, we would observe that G c increases with the system size up to an asymptotic value which is the relevant one at larger scales. Here, the high degree of cross-linking of CS1000 strongly limits the process zone size to the nanometer scale.

Conclusion

A careful analysis of molecular simulations of material fracture failure, in which we pay a special attention to the size effects, shows that usual LEFM predictions are still valid if the process zone size is small at the scale of molecular simulations. We illustrate this situation for a silica crystal, α-cristobalite, for which the process zone is only a few atoms large. Provided that boundary conditions are accounted for (periodic boundaries in our work), LEFM captures correctly the global material be-havior (remote stress and average strain). In particular, one can directly estimate material fracture properties, such as toughness, from molecular simulation. In contrast, when the process zone occupies a significant fraction of the molecular system, LEFM fails to predict the mechanical behavior. We illustrate this second situation with the example of a nanoporous polymer, CS1000, for which the process zone is about one order of magnitude larger than that of α-cristobalite. LEFM fails to predict the behavior of CS1000 both quantitatively and qualitatively: the LEFM scaling between critical stress and crack size (inverse square root) becomes linear, the stability of CS1000 failure does not follow LEFM analysis, and LEFM estimates of toughness from critical stress is not consistent with the total energy released during failure. To overcome the inherent limitations of LEFM in presence of large scale yielding, we develop an alternative analysis which includes cohesive zone models at the crack tips. Provided the cohesive law is well adapted, we recovered all the features of the mechanical behavior of CS1000 within this framework. From the cohesive law, one can then easily recover the macroscopic toughness and anticipate the mechanical behavior at larger scales. Thus, our analysis with cohesive zone models offers the possibility to upscale fracture properties from molecular simulation even with large yielding at the scale of the molecular system. It should be noted however that, the approach is limited by the size of the molecular system simulated since the process zone cannot extend further.

The process zone size is on the order of the square of the ratio between toughness and yield stress. Thus, estimating this ratio is critical to know whether or not the scale of molecular simulations is relevant, and to decide what kind of analysis to perform.

Fig. 1

 1 Fig. 1 Schematic representation of the simulated systems.We consider periodic boundary condition and increase the box size to load the system until complete failure. The overall stress Σ is computed.

Fig. 2

 2 Fig. 2 Three of the systems simulated for α-cristobalite.

Fig. 3

 3 Fig.3Results of the molecular simulations for three of the α-cristobalite systems studied. The average tensile stress Σ in the loading direction is represented as a function of the strain imposed to the systems. The snapshots of molecular configurations illustrate the case of system H.

Fig. 4

 4 Fig. 4 Three of the systems simulated for CS1000.

Fig. 5

 5 Fig.5Results of the molecular simulations for three of the CS1000 systems studied. The average tensile stress Σ in the loading direction is represented as a function of the strain imposed to the systems. The snapshots of molecular configurations illustrate the case of system N'.

  +jL y Im (Z (x + i (jL y ) , u)), with Z = dZ dz and Z the Westergaard stress function below for collinear periodic cracks problem:

Fig. 8

 8 Fig.8Critical stress Σ cr obtained from molecular simulation for α-cristobalite (top) and CS1000 (bottom) as a function of the effective half length of the crack a ef f . The results for α-cristobalite respect LEFM theory, whereas the results for CS1000 differ significantly from LEFM at low crack lengths.

Fig. 9

 9 Fig. 9 Loading curves from molecular simulations compared with LEFM predictions for the three systems of Figure 2.

Fig. 10

 10 Fig.10Schematic representation of the cohesive zone model. Cohesive zones develop behind crack crack tips in which a force is transmitted across the crack, depending on the crack opening. The cohesive zone suppresses the stress singularity at the crack tip. At large openings, the cohesive forces vanish and the crack is unbridged. Loading proceeds in three phases: 1) the cohesive zone grows until the opening at the end of the zone reaches the critical opening, 2) the unbridged crack propagates, 3) the cohesive zones coalesce at the periodic boundaries. In the third phase, the mechanical problem is that of an infinite plate submitted to periodic traction loads on its faces.

  y) = (A + Bξy) cosh (ξy) + (C + Dξy) sinh (ξy)(20) where A, B, C and D are constants that depend on the boundary conditions. We consider a plate submitted to an even vertical traction load p (x) on its upper and lower faces (y = L y /2 and y = -L y /2, considering the middle of the plate as the origin). Introducing these boundary conditions in Equations 17 and 18 leads to the following values for the constants: A = x) cos (ξx) dx is the cosine Fourier transform of the load p. The inverse transform is p (x) = 2 π ∞ 0 p (ξ) cos (ξx) dξ. Introducing the values of the constants into the expression of G (Eq. 16), we obtain the expression of the vertical displacement u f aces y (x) = u y (x, y = b) = -u y (x, y = -b) (Eq. 19) on the upper and lower faces: , the load p is a periodic function of period L x . Therefore, the cosine Fourier transforms can be replaced with Fourier series with the Fourier coef-

Fig. 14

 14 Fig. 14 Loading curves from molecular simulations compared with the cohesive model predictions for the three systems of Figure 4.

Table 1

 1 Systems considered for simulation of fracture in the (001) plane of α-cristobalite

	System	Lx ( Å)	Ly ( Å)	Lz ( Å)	Initial crack size ( Å)	Critical stress Σcr (GPa)
	A †	49.72 55.38 14.92	10	19.97
	B †	49.72 55.38 14.92	20	15.73
	C †	49.72 55.38 14.92	30	10.84
	D †	74.58 69.22 24.86	20	14.99
	E †	74.58 69.22 24.86	30	11.90
	F †	74.58 69.22 24.86	40	9.77
	G †	99.43 83.07 19.89	30	12.99
	H †	99.43 83.07 19.89	40	10.84
	I †	99.43 83.07 19.89	50	7.65
	J †	99.43 69.22 4.97	30	12.21
	K ‡	99.43 69.22 4.97	30	15.92
	L ‡	89.49 62.30 4.97	30	11.20
	M ‡	74.58 55.38 4.97	30	10.49
	N ‡	59.66 41.53 4.97	20	14.24
	O ‡	49.72 34.61 4.97	20	17.28
	P ‡	39.77 27.69 4.97	15	15.78
	Q ‡	29.83 20.77 4.97	12	20.77
	R ‡	29.83 20.77 14.92	7.46	25.35
	S ‡	29.83 27.69 14.92	7.46	24.66
	T ‡	29.83 34.61 14.92	7.46	25.40
	U ‡	29.83 41.53 14.92	7.46	24.64
	V ‡	44.75 34.61 14.92	7.46	20.17
	W ‡	59.66 27.69 14.92	7.46	24.06
	X

‡ 59.66 41.53 14.92 7.46

Table 2

 2 Systems considered for simulations of fracture in CS1000

	System	Lx ( Å)	Ly ( Å)	Lz ( Å)	Initial crack size ( Å)	Critical stress Σcr (GPa)
	A' †	50	50	25	10	15.53
	B' †	50	50	25	15	14.15
	C' †	50	50	25	20	13.83
	D' †	50	50	25	25	10.89
	E' †	50	50	25	30	9.44
	F' †	75	75	25	20	13.52
	G' †	75	75	25	25	13.49
	H' †	75	75	25	30	11.43
	I' †	75	75	25	35	10.63
	J' †	75	75	25	40	8.94
	K' †	100	75	25	30	12.41
	L' †	100	75	25	35	12.57
	M' †	100	75	25	40	11.99
	N' †	100	75	25	45	9.59
	O' †	100	75	25	50	9.31
	P' ‡	50	50	25	20

  Fig. 12 Loading curves and critical as predicted with a linear cohesive zone model for various ductility ratios δ * cr /σ * Y S . The notations are the same as in Figure

	Σ * /K * Ic	0.4 0.6 0.8 1.0	δ * cr /σ * YS 0.5 1 1.5 2 2.5 3		LEFM
		0.2
		0.0 0.0	0.2	0.4	0.6	0.8 ∆ * /K * Ic	1.0	1.2	1.4
	Σ * cr / G * c	1.0 1.5 2.0 2.5					Infinite body LEFM	δ * cr /σ * YS 0.4 0.8 1.2 1.6 2 2.4 2.8
		0.5
		0.0 0.0	0.2		0.4	2a/L x	0.6	0.8

Acknowledgements Funding for this work through the X-Shale project enabled through MIT's Energy Initiative, with sponsorship provided by Shell and Schlumberger, is greatly appreciated.