Modeling Brittle Fractures in Epoxy Nanocomposites Using Extended Finite Element and Cohesive Zone Surface Methods

Linear elastic fracture modeling coupled with empirical material tensile data result in good quantitative agreement with the experimental determination of mode I fracture for both brittle and toughened epoxy nanocomposites. The nanocomposites are comprised of diglycidyl ether of bisphenol A cured with Jeffamine D-230 and some were filled with core-shell rubber nanoparticles of varying concentrations. The quasi-static single-edge notched bending (SENB) test is modeled using both the surface-based cohesive zone (CZS) and extended finite element methods (XFEM) implemented in the Abaqus software. For each material considered, the critical load predicted by the simulated SENB test is used to calculate the mode I fracture toughness. Damage initiates in these models when nodes at the simulated crack tip attain the experimentally measured yield stress. Prediction of fracture processes using a generalized truncated linear traction–separation law (TSL) was significantly improved by considering the case of a linear softening function. There are no adjustable parameters in the XFEM model. The CZS model requires only optimization of the element displacement at the fracture parameter. Thus, these continuum methods describe these materials in mode I fracture with a minimum number of independent parameters.


Introduction
Material ensembles such as adhered layers, fiber-reinforced polymer composites (FR-PCs), and multi-layer coatings serve a wide range of industries including: construction, military, transportation, and aerospace [1,2]. In such industries, epoxy resins demonstrate remarkable utility due to their mechanical properties, chemical resistance, and sustainability. Epoxy resins are environmentally benign, a property that is attributed to curing processes that render a thermo-mechanically stable material. However, these curing processes yield potentially brittle materials, which is exemplified by the low strain-to-failure capacity of FRPCs [3,4]. There is hence an immediate need to develop efficient damage prediction tools for epoxy resins to guide the formulations of tough composites.
The field of theoretical continuum mechanics provides reliable techniques for simulating fracture processes. We assume that a linear elastic fracture mechanics (LEFM)-based approach is applicable for epoxy resins because the calculated size of their crack tip plastic zones was found to be small compared to common specimen dimensions and it is confined to a region close to the crack tip [2]. This study reports the use of two methods to model fracture processes available in the commercial software Abaqus (3DS-SIMULIA) [5]. The extended finite element method (XFEM) circumvents the need for mesh refinement in fracture process zones and conveniently inherits model definitions from conventional finite element methods [6]. XFEM provides an efficient framework for failure analysis when either an accurate description of in situ failure is available or a user-defined criterion is implemented to tame the over-prediction of failure conditions [7]. Cohesive zone models are also commonly used to model fractures in either an element-based or surface-based (CZS) implementation. Compared to XFEM, the cohesive zone methods are disadvantaged by having more stringent mesh dependencies for attaining convergence and by requiring predefined crack paths, which are difficult to ascertain in many realistic systems. In XFEM and CZS simulations, the cohesive traction-separation law (TSL) functions as the constitutive relationship that governs the crack initiation and propagation processes. Conveniently, TSLs can be customized to model a wide spectrum of constitutive equations that describe different classes of materials, such as glassy poly(methyl-methacrylate), concrete, and steel [8]. Recent studies using XFEM to model the failure of composite structures demonstrate the importance of selecting damage initiation and damage evolution criterion that are appropriate to the materials and systems studied, such as transverse plies of cross-ply carbon fiber/epoxy composites [9], glass fiber/polyvinylesters [7], and heavily crosslinked epoxy [4]. To the best of our knowledge, these methods have not been applied to the materials considered in this study.
We developed the methodologies herein with the intention of modelling the fracture processes exhibited by epoxy resin nanocomposites studied by Bain et al. [2]. We use the XFEM and CZS techniques to model the mode I fracture exhibited by samples of epoxy resins during the quasi-static single-edge notched bending (SENB) test. The resins are composed of diglycidyl ether of bisphenol A (Hexion Inc.) that has been cured with Jeffamine D-230 (Huntsman Chemical), a diamine propylene glycol oligomer with a molecular weight of 230 Da, and filled with varying concentrations of core-shell rubber nanoparticles MX-125 or MX-257 (Kaneka USA). The materials studied include two unfilled D230 epoxies, namely D230-1 and D230-2, and filled nanocomposites D230 + 1 wt% MX-125, D230 + 5 wt% MX-125, D230 + 1 wt% MX-257, and D230 + 5 wt% MX-257. The chemical structures of the two unfilled D230 materials are identical, except that the D230-1 was processed at a higher curing temperature and therefore has a higher glass transition temperature (T g = 99 • C) compared to the D230-2 (T g = 96 • C).
The predicted load-displacement behavior and mode I fracture toughness show satisfactory agreement with the experimental results for all simulations presented. The proposed methodologies are advantaged by circumventing user-defined subroutines and excessive non-physical hyperparameters. The fracture processes analyzed are induced by tensile stresses at the crack tip and hence uniaxial tensile data are sufficient for parameterizing the fracture model. Our models require three physical parameters, while the CZS model uses an additional fitting parameter, to accurately predict trends in load-displacement relationships and in the fracture toughness among the considered epoxy resins. All parameters are discussed in detail in Section 2. These simple LEFM-based methodologies can be generalized to predict the same trends among other classes of materials while relying on limited experimental data.

Fracture Model and System Layout
To model fracture in epoxy nanocomposites, we used constitutive models that include parameters formulated from the results of experimental SENB and uniaxial tensile tests. We conducted these tests according to the procedures adopted by Bain et al. [2]. Nominal stressstrain curves were determined for the epoxy resins from tensile tests that were conducted and post-processed according to ASTM standard D638-14 (see Figure S1, Supplementary Materials) [10]. However, we assumed that all materials behave as linear elastic solids and account for both plasticity and ductility by using TSLs to predict material damage and the separation between finite elements. Fracture toughness quantifies the resistance of a sample to fracture. Determining the fracture toughness of a sample subjected to the SENB test requires measuring the load applied to the moving center pin of the three-point bending test apparatus at the critical point. The critical point of a SENB test is indicated at the maximum of the loaddisplacement curve, at which unstable fracture initiates upon further displacement of the center pin. The critical stress intensity factor K IC , also known as the fracture toughness, is calculated as [11] where P is the load applied at the critical point, B is the material depth or thickness, W is the material width, and a is the crack length. For materials obeying LEFM, the crack length is assumed to remain constant until the point of unstable fracture. The function f a W is defined according to ASTM standard D5045-14 [11].
One goal of this study is to validate using LEFM assumptions in modelling the quasistatic SENB test. Figure 1 presents the configuration of a SENB test simulated using the Abaqus software. The dimensions of all modeled samples comply with the guidelines for the SENB test for plastic materials established by ASTM standard D5045-14 [11]. We performed simulations for six epoxy resin samples possessing different initial crack length a. Note that the initial crack of a sample is modeled as a seam line at which the adjacent material surfaces are unbound to each other. The dimensions of all samples modeled in this study are presented in Table 1, except for the width W, which was assumed to be 12.7 mm for each sample. This is a reasonable assumption because the experimental measurements differ from this dimension by 1% at most. Note that Figure 1 does not show the depth of the sample because all samples were modeled with an implicitly in-plane depth B. Furthermore, these samples' thicknesses are at least 50x the theoretical plastic zone radius size and thus the crack tip region is in a state of plane strain deformation [2,11]. Therefore, our assumption that the samples are in plane strain deformation during the SENB test is justified.  In this study, the tensile yield stress is used to inform the criterion for initiating fracture between finite elements. This is justified by considering that the yield stress of a material represents the stress at which significant plastic deformation occurs and that concentrated regions of both stress and plastic strain form ahead of a crack tip when a load is applied in the SENB test. The yield stress of each material is assumed to be the local maximum stress from the experimental stress-strain curves. The resulting yield stresses are listed in Table 2.
In the case of unfilled D230, the reported yield stress is the average of the measured yield stresses of the D230-1 and D230-2 materials. For each material the yield stress reported in Table 2 is included as a parameter in the TSLs. Fracture data for epoxy resins obtained from experimental SENB tests were used to formulate criteria for predicting fracture initiation and crack propagation. Poisson's ratio and an effective elastic modulus were the physical parameters used to define the materials studied. A Poisson's ratio of 0.35 was assumed for every material to be consistent with experiments conducted by Bain et al. [2]. For each epoxy system, the experimentally determined mode I fracture toughness K IC and critical strain energy release rate G IC (Table 3) were used to determine [11] the effective elastic modulus E eff . This definition of E eff sufficiently describes brittle isotropic solids under plane strain deformation [12]. G IC informs the crack propagation (damage evolution) because it is the rate of energy change corresponding with the increase of surface area that occurs during mode I crack growth. Note that Equation (2) provides a lower bound for E eff because K IC represents a lower limiting value of the mode I fracture toughness [11]. The experimental fracture test results and calculated E eff are provided for each material in Table 3. These physical parameters were used to model the SENB test for each material considered in this study.

Surface-Based Cohesive Zone Model
The cohesive zone model incorporates discontinuities in materials by forming a cohesive damage zone that allows for the modelling and analysis of interfacial fracture [13]. The CZS method is often used to study delamination of surfaces [14], unlike the elementbased cohesive zone model approach frequently used to study adhesive layers [15]. In our CZS simulations, we define contact pairs between material interfaces to represent the cohesive zone, as performed by Li et al. [16]. Damage evolution (crack propagation) is governed by a TSL that represents the functional relationship between the traction T r and displacement between nodes at the crack tip ∆. Fracture initiates after stresses at elements satisfy a pre-defined criterion and the progress in separation between elements is governed by a softening function. The quadratic nominal stress (QUADS) and maximum principal stress (MAXPS) fracture initiation criterion were used for CZS and XFEM simulations, respectively. A softening function can be selected to better simulate the expected brittle or ductile type of fracture process [17]. Furthermore, the fracture processes exhibited by similar materials tend to be well-described by TSLs using the same class of softening functions. Our work adapts a generalization of trapezoidal softening functions, which are often used to model ductile materials [18]. Commonly used softening functions are defined in terms of hyperparameters that introduce additional fitting capabilities.
This study assesses the applicability of a truncated linear softening function in modeling the brittle fracture of epoxy nanocomposites. This approach has been reviewed previously for other materials [8]. Considering mode I fracture was being modelled, the area enclosed by the softening function and loading portion of the TSL (Figure 2) was equivalent to the critical strain energy release rate G IC [11]. Furthermore, the interface stiffness K determines the initial slope of the TSL. We explored the proper selection of values for the stiffness parameter K as well as its connection to the TSL [19]. Additionally, a softening function can be adjusted to model the type of material in question more accurately [8,19]. The general form of the truncated linear TSL is defined as, where the hyperparameter c is a geometrical descriptor that determines the point at which the softening function is truncated, ∆ T is the displacement at fracture initiation T K , and ∆ c is the displacement at complete fracture. At the region of TSL in which no damage is done to the material (0 < ∆ < ∆ T ), the displacement between nodes is miniscule and entirely reversible. The peak traction T is used to define the QUADS parameter in Abaqus simulations and thus it informed the damage initiation criterion of our CZS model. The QUADS parameter is calculated as the sum of quadratic ratios between the measured stresses and peak nominal stresses in three dimensions, which includes tractions normal to the crack interface and in two shear directions. We assigned the tensile yield stress as the peak nominal stress of each dimension. The area enclosed by the TSL is constrained to be the value G C [8] and thus for any valid value of c (0 < c ≤ 1), the TSL needs to be adjusted to uphold this constraint. Furthermore, the interfacial strength T (i.e., the maximum traction value of the TSL) is assumed to be equivalent to the tensile yield stress of a material. Thus, the displacement ∆ c can be calculated using c, stiffness K, and traction T as For the truncated linear softening function, the damage variable D is calculated as: where the input data are included in Abaqus simulations as tabular values of D as a function of ∆ − ∆ T . Multiple displacement values between ∆ T and ∆ c with corresponding D values are necessary to include in the definition of the TSL to obtain well-converged results. Elements along the crack that are separated such that ∆ T < ∆ < ∆ c are partially damaged (0 < D < 1). Displacements less than ∆ T indicate that the relevant elements are undamaged (D = 0), whereas displacements greater than or equal to ∆ c indicate that the relevant elements are separated and completely damaged (D = 1). Elements along the initial crack do not experience cohesive interactions and are considered completely damaged.
In implementing the truncated linear TSL, we observed that the hyperparameter of c = 1 led to the most accurate prediction of fracture data for all the materials considered. We have decided to report only the results obtained from this case because we are primarily interested in the ability of these methods to accurately model the fracture processes of epoxy resins. Thus, we must emphasize that a linear TSL was used in all simulations for the data reported in this study. Table 4 provides all parameters that define the damage initiation and evolution criteria for the CZS simulations conducted. Note that in the case of CZS simulations, the high stress concentrations surrounding the crack tip (the fracture process zone) necessitate fine mesh resolutions to accurately model fracture processes. It is observed that the quality of such simulations depends on the number of mesh elements within the fracture process zone [19]. The stiffnesses K were selected to approximately fit the experimental fracture data for each material.

Extended Finite Element Method
The XFEM method adds enrichment functions to the conventional finite element mesh to allow discontinuities in the displacement field u(the crack domain Γ c ) and singularities that occur at crack tips [4,5]. The discretized displacement field u h in XFEM is thus written as where x denotes material coordinates, I is the set of all nodes in the mesh, J is the set of crack tip nodes, K is the set containing related enrichment nodes, u i is the displacement field degrees of freedom at the node i, N n (x) is the shape function associated with the node n, H(x) is the step function, and b j is the nodal enrichment degrees of freedom which represent the amount of displacement in the crack line. The last expression in Equation (6) is the crack enrichment term, which contains the crack tip asymptotic functions F l (x), and the crack tip-enriched nodal degrees of freedom c l k , which represent the displacement in the crack tip [20].
In our XFEM simulations, the fracture process was governed by a linear TSL. For each material, the stiffness K was equivalent to the effective elastic modulus (Table 3) divided by 1 mm, which is a default setting available in the Abaqus software for XFEM simulations. The critical strain energy release rate was calculated using the Benzeggagh and Kenane approach [21], while the energy released per unit area during pure mode I fracture was equivalent to G IC ( Table 3). The tensile yield stress, converted to true stress, was assigned as the MAXPS fracture initiation criterion parameter for all XFEM simulations. The progress of a fracture was monitored using the damage variable STATUSXFEM in Abaqus simulations, which changes from 0 to 1. A STATUSXFEM value of 1 indicates a complete separation of crack face elements.

Simulation Parameters
Significant computational model parameters including details on meshing conditions are listed in Table 5. In both the XFEM and CZS models, we assigned a uniform mesh to the bulk material, and in the CZS model, we included a fine mesh along the predefined crack path. In all simulations, we defined the crack as a line segment with an infinitesimally small width and the initial length was observed in the experimental SENB test. Note that we simulated a version of the SENB test that is designed such that the bottom pins are rotatable with as minimal friction as possible [11]. We performed simulations with and without rotatable pins, and did not observe significant differences in results. We have presented data obtained from modelling the two adjacent pins as explicitly rotatable because this feature introduced additional control over the emulation of experimental conditions. Furthermore, the sliding contact between all pins and material surfaces had a friction coefficient of 0.1 and the punch (top pin) was modeled as a rigid analytical surface. We also observed that the prediction of the critical point was insensitive to the material density. This implies that the mode I fracture toughness predicted for each material is independent of material density and thus the assumption that every material possessed the same density is acceptable. Hence, we assumed that the density of each material was equivalent to the experimentally measured density of unfilled D230 resin measured as 1.156 × 10 −9 tonne/mm 3 (1.156 g/cm 3 ). Square mesh elements were used in each simulation. Mesh elements in the bulk material had edge lengths of 0.05 mm in the XFEM simulations and 0.2 mm in the CZS simulations. Extensive validation of the mesh design was performed for the XFEM model, documenting that the selected mesh yields converged results (see Figures S2 and S3). The CZS simulations included a sub-laminate layer containing smaller mesh elements, as depicted in Figure 3. This adjustment was required to obtain convergent CZS simulations and better resolution for monitoring stress distributions near the crack tip during the fracture process. We found that it was necessary for the sub-laminate layers to be wider than the length of the plastic zone ahead of the crack tip predicted under plain strain deformation in the Tresca yielding criteria [22] for unfilled D230: where σ ys is the tensile yield stress ( Table 2) and K IC is the mode I fracture toughness (Table 3) of unfilled D230. The predicted plastic zone length for unfilled D230 was 47 µm and hence our CZS simulations contained 50 µm sub-laminate layers along each side of the crack. Using the predicted plastic zone of unfilled D230 as a constraint for the sub-laminate layer width is justified because our methodologies assume that all materials are brittle linear elastic materials. The square mesh elements within the sub-laminate layers were assigned edge lengths of 0.005 mm, which was found to be sufficiently small for convergent simulations. Furthermore, we found that smaller sub-laminate mesh elements led to a numerical stability that was exhibited as hour-glassing.

Results and Discussion
The XFEM and CZS methods were used to simulate the brittle mode I fracture of six epoxy nanocomposites that occurs during the quasi-static SENB test. From each simulation, a load-displacement plot was generated and compared with experimental results. We identified the critical load that occur during the SENB test when there is complete separation between elements at the crack tip and where unstable fracture initiates. The critical point is identified when the parameter STATUSXFEM attains a value of 1 during XFEM simulations or, in the case of the CZS simulations, when the separation between elements at the crack tip become larger than the ∆ c value listed in Table 4. The load applied to the center pin at the critical point is used to calculate the mode I fracture toughness K IC according to Equation (1). The maximum load and corresponding load-line displacement (LLD) are identified and labeled as P max and d c , respectively, to compare with experimental results. The results are summarized in Table 6 and load-displacement curves predicted using XFEM and CZS are shown in Figures 4-6 with the experimental data. For each material, these simulations approximately predict the critical point of the SENB test and the mode I fracture toughness K IC . One can see that the simulations follow the trends observed in the experiments, as shown in Figures 7-9.       These simulations successfully predicted that brittle fracture processes occur in each material studied when subjected to the SENB test. The proposed methodologies quantitatively predict the critical point exhibited in the SENB test. Therefore, a linear elastic constitutive model that accounts for fracture using a linear TSL for modelling fracture is sufficient for simulating brittle fracture processes exhibited by epoxy resins. Furthermore, the yield stress determined from the experimental uniaxial tensile tests suffices for defining the criterion for the initiation of brittle fracture. Overall, there is good quantitative correlation between the experimental results and both XFEM and CZS predictions for both the unfilled and rubber-filled epoxy materials. This approach does not describe complex micromechanical toughening mechanisms and plastic deformations at the crack tip, but rather offers a sufficient approximation at a continuum level to predict fracture for the class of epoxies studied.
Error propagation analysis indicates that the experimental data were adequately reproduced by the XFEM and CZS methods. Percent errors between the simulated and experimental results presented in Table 6 were less than 10% in all cases other than the XFEM simulation for D230 + 5 wt% MX-257, which attained values up to 11.6%. In this case, the expected uncertainty for predictions was approximately 12%, which is the quadrature sum of the relative standard deviations of the experimentally determined tensile yield stress and elastic modulus. Our model inherits this uncertainty due to incorporating the tensile yield stress in the fracture initiation criteria and due to using both the experimentally determined mode I fracture toughness K IC and critical strain energy release rate G IC in Equation (1) to calculate the effective elastic modulus E eff [11]. Note that this argument assumes that the predicted peak load, critical displacement, and mode I fracture toughness change linearly with respect to the tensile yield stress and elastic modulus, which is reasonable for small deviations. In general, our predictions are within 2.4% of the uncertainties derived from the standard deviations of the experimentally determined yield stresses and elastic moduli. Note that discrepancies between our CZS and XFEM-predicted results are not attributed to the different degrees of convergence because our simulations are fully converged with respect to the meshes used. The XFEM and CZS models also have different degrees of freedom to fit experimental data. In the CZS simulations, the stiffness parameter K was adjusted to enable closer agreement with experimental results. In the XFEM model, the stiffness parameter was equal to the effective elastic modulus and was not adjusted to fit the experimental data. This difference between these models contributes to the discrepancies of their predictions. Overall, our models sufficiently model the SENB test and predict the corresponding fracture properties for each material considered.
Our models tended to predict the critical loads more accurately than the critical LLDs. On average, the error of the predicted critical loads was approximately 4.5 times smaller than the error of the predicted critical LLDs. In our simulations the pins are treated as rigid, non-deformable objects. Therefore, the system compliance observed in experimental SENB tests, which additively contributes to the measured LLDs, did not affect our simulations [11]. Our models also did not account for the formation of biaxial stress states and material deformations that formed in the crack plane, i.e. the face formed by the crack length and width [23]. This resulted from our assumption that each sample exhibits a 2D homogeneous plane strain condition. SENB tests modelled under this assumption tend to overestimate the material stiffness [24]. Consequently, our underestimation of critical LLDs is attributed to the assumption that the samples exhibit the 2D homogeneous plane strain condition.
We calculated principal stress contours and monitored the separation between the elements at the initial crack tip. Figure 10 includes results obtained from analysis of the CZS simulation for D230 + 5 wt% MX-257. Figure 10a shows a load-displacement curve with superimposed contours of principal stress that enclose the fracture process zone. Prior to the initiation of the unstable fracture, the fracture process zone increased in size with the applied load, reaching a maximum value of 257.0 N ( Table 6). The size of the fracture process zone remained approximately constant as the crack continued to progress. Although we did not explicitly employ plasticity models in the simulations, the formation of the fracture process zone and plastic zone that formed around the crack tip correspond to the increased toughness of the nanocomposites. We also found that the maximum value of the principal stress was attained close to the critical point of the load-displacement curve. The separation between the crack faces at the initial crack tip increased with the applied load, but the crack did not propagate until the critical point was reached. The critical point is reached when the distance between elements at the initial crack becomes larger than the ∆ c value of 0.0791 mm (Figure 10b). The corresponding load of 256.9 N was very close to the observed maximum load of 257.0 N. Using the value of 256.9 N to calculate the K IC value according to Equation (1) accurately predicts the experimental result, as shown in Table 6. The separation between elements at the crack tip increased non-linearly ( Figure 10b) and fracture occurred at d c = 0.598 mm ( Table 6). We also calculated the size of the fracture process zone at the critical point for each polymer studied. We observed that the size of the fracture process zone increased with the toughness of the polymer nanocomposite system. Detailed analysis and quantitative assessment of the fracture process zone and plastic zone would require a much denser mesh than employed in this study. While the current approach effectively replicates the increased toughness of nanocomposites, which is well known to occur due to complex micromechanical plastic deformation processes near the crack tip [25], these models do not address this phenomenon directly. A more detailed approach could entail using a hyper-elastic-based constitutive model that predicts plastic deformations at the crack tip, particularly for rubber-filled materials [26]. One could also directly account for the plasticity in the constitutive models used for these materials. Material plasticity could be indirectly treated by accounting for plastic strains that occur during the uniaxial tensile test [27]. The deformed cross-sectional areas of tensile test samples, which are required to exactly calculate the true stress, are unavailable and thus stress values beyond the ultimate yield stress could be extrapolated [28]. Measurements of the transverse strain in one or two dimensions using high resolution digital image correlation of uniaxial tensile tests would allow for the more direct measurement of Poisson's ratio and would enable the use of complex empirical plasticity models [29]. The accuracy of our predictions could also be improved by optimizing the yield stress, which is used to define the fracture initiation criterion. Optimizing the offset that is used to determine the yield stress of each sample from the experimental tensile data considered would account for the varying degrees of plasticity among the materials considered [10,30]. Note that the purpose of this study is to validate continuum methods using a minimal number of fitted parameters and hence we did not pursue this additional optimization.

Conclusions
This study demonstrates two LEFM-based methodologies in modelling the brittle mode I fracture of epoxy resin nanocomposites observed during the quasi-static SENB test. The XFEM and CZS methods are comparably robust and successful in modelling the fracture processes exhibited by the epoxy resins considered. These methods accurately predict the critical point of the SENB test, which enables the accurate determination of mode I fracture toughnesses K IC . This overall agreement between the predicted and experimental values justifies applying the present methods in fracture modeling of both brittle and toughened epoxy resins.
We have proposed methodologies that are well advantaged by their simple parameterization and reliance on minimal independent empirical data. Damage initiates upon stresses at the crack tip, attaining the yield stress, which is determined from experimental uniaxial tensile tests. This approximation is sufficient for all epoxy systems considered in this study. A general form of truncated linear TSL was considered and the best results were obtained using an entirely linear formulation of the damage model. There are no adjustable parameters of the XFEM model and the CZS model requires optimization of only the stiffness parameter K, which defines the displacement between elements at complete fracture ∆ c . The demonstrated predictive capability of the CZS model asserts that this parameterization is consistent. The displacement ∆ c was determined through optimization of the TSL and was used to predict the onset of complete fracture between nodes at the simulated crack tip. The corresponding critical load was used in the calculations of the fracture toughness K IC . Calculating stress contours at the crack tip enables monitoring the progression of the unstable brittle fracture process.
This research study is concerned with the validation of continuum methods for modelling fracture properties of epoxy resins that depend on minimal experimental parameters. Neither method used in this study required extensive customization or detailed plasticity modeling. However, accounting for toughening mechanisms at the microscale was not undertaken in this work [31,32] and is left for further developments. Future work may also involve refining the TSL by either further optimizing the softening function or by using an established inverse method to fit a TSL to experimental data [33]. Our simulations also do not account for complex phenomena such as the occurrence of modes of plastic deformation, crack tip-blunting, and explicit dependencies on both crack tip position and velocity.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/polym13193387/s1, Figure S1: Engineering stress-strain curves determined from uniaxial tensile tests for each material, Figure S2: Convergence of critical load in XFEM simulations with respect to mesh size, Figure S3: Convergence of critical load-line displacement in XFEM simulations with respect to mesh size.

Data Availability Statement:
The data presented in this study are available upon request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.