Nonlinear XFEM Modeling of Mode II Delamination in PPS/Glass Unidirectional Composites with Uncertain Fracture Properties

Initiation and propagation of cracks in composite materials can severely affect their global mechanical properties. Due to the lower strength of the interlaminar bonding compared to fibers and the matrix, delamination between plies is known to be one of the most common failure modes in these materials. It is therefore deemed necessary to gain more insight into this type of failure to guide the design of composite structures towards ensuring their robustness and reliability during service. In this work, delamination of interlaminar bonding in composite end-notched flexure (ENF) samples was modeled using a newly developed stochastic 3D extended finite element method (XFEM). The proposed numerical scheme, which also incorporates the cohesive zone model, was used to characterize the mode II delamination results obtained from ENF testing on polyphenylene sulfide (PPS)/glass unidirectional (UD) composites. The nonrepeatable material responses, often seen during fracture testing of UD composites, were well captured with the current numerical model, demonstrating its capacity to predict the stochastic fracture properties of composites under mode II loading conditions.


Introduction
Nowadays, fiber reinforced polymer (FRP) composites are widely used in various engineering applications, including aeronautical, marine, and automotive industries. These materials have high strength-to-weight ratios as well as good corrosion resistance and can be engineered based on the required strength or performance objectives of a given application. Although FRP composite structures have proven to provide numerous advantages, initiation and propagation of cracks in these materials can drastically affect their mechanical properties. The most common failure modes in these composites are classified as fiber breakage, fiber pull-out, matrix cracking, and interlaminar delamination. Among these, interlaminar delamination is perhaps the most common type and may occur because of weak bonding between composite layers. This failure mechanism can significantly reduce the structural stiffness of FRP structures and weaken their tensile or shear capacity under service loads [1].
To improve the mechanical performance of FRP composites in the presence of process-induced or loading-induced cracks, extensive studies on their fracture properties have been performed, both experimentally and numerically [2][3][4][5]. Gaggar and Broutman [2] utilized both single and double edge-notched tensile tests as well as a notched bend test to extract the critical stress intensity factors of these cracks in composite materials. Mower and Li [3] summarized the experimental results 2 of 13 from previous investigations and concluded that linear elastic fracture mechanics (LEFM) is not a valid approach for long-fiber composites and that a nonlinear material constitutive model is instead required to accurately characterize their fracture response. In an attempt to measure mode II fracture toughness, Russel [4] proposed end-notched flexure (ENF) testing. Murri and O'Brien [5] employed ENF samples to study the mode II fracture toughness of FRP composites while investigating the error associated with neglecting nonlinear terms in the calculation of strain energy release rate. ENF testing is currently one of the most recognized experimental methods for studying the mode II fracture of FRP composite materials.
Numerical models of crack initiation and propagation during ENF testing of composites have also been developed in the past (e.g., Harper and Hallett [6] and Fan et al. [7]) to facilitate the design of composite structures prior to prototyping and testing stages. However, these investigations do not typically account for the heterogeneous and stochastic properties of composites, which are commonly observed during material and product testing. There can be a variety of sources for such nonuniformity of material properties in composites. Examples include random distributions of fibers within samples, fiber penetration between layers, existence of voids within the matrix, human error in manufacturing process, and uneven heating or cooling of samples during molding. Adopting deterministic approaches and ignoring the spatial variability that exists in composites can introduce errors into large-scale simulations. Consequently, stochastic modeling of effective material properties appears to be essential for more precise assessment of the mechanical behavior of composites, especially during the prediction of critical failure loads and crack formation patterns [8,9].
Among recent stochastic modeling works in the field, Ashcroft et al. [10] introduced microstructure randomness in the fracture properties of carbon fiber reinforced polymer composite materials for finite element simulation of double cantilever beam (DCB) tests using interface cohesive elements. Nonuniformity and random distribution of material fracture properties were considered by means of uniform and Weibull distributions. Jumel [11] employed a finite difference numerical method to study crack initiation and propagation in DCB specimens with randomly fluctuating interface properties along the crack path. The effect of this variability at the microscopic level on the parameters measured at the macroscopic level was investigated. The results emphasized the need for further development of computational approaches that account for randomness and variability in material properties.
The present study was aimed at developing and examining an enhanced numerical approach for simulating fracture in unidirectional (UD) composite structures by considering both material and geometric nonlinearities along with stochastic fracture properties. An Abaqus user element subroutine was developed and linked to MATLAB to model under ENF testing (mode II fracture). Here, we adopted the extended finite element method (XFEM) and enhanced it with contact and cohesive zone modeling capabilities along with stochastic fracture properties to improve the simulation of crack propagation in the composite laminates. The cohesive zone model constituted a bilinear traction-separation law at the crack front, which enabled modelling of the fiber bridging mechanism in the process zone ahead of the crack tip during loading. Stochastic distributions of the fracture properties were captured within the bilinear traction-separation law. The numerical results were compared with a set of tests performed on polyphenylene sulfide (PPS)/glass UD composites.

Experimental
PPS/glass UD composite samples comprising 14 plies and with nominal dimensions of 250 mm length, 25 mm width, and 3 mm thickness were prepared for ENF testing. To introduce an edge crack with a nominal size of 43 mm in each sample, a polyimide Teflon sheet was placed between the middle layers of the stacked laminates before placing them in the oven for curing.
Cured specimens with preinserted delamination were put into a three-point bending test fixture and made to undergo mid-span deflection at a rate of 2 mm/min on the grips. At the onset of delamination extension, the force on the loading cell was recorded while the crack was permitted to propagate ( Figure 1). extension, the force on the loading cell was recorded while the crack was permitted to propagate ( Figure  1). During ENF testing, crack propagation occurs due to excessive shear deformation, thus exhibiting mode II failure with an abrupt nature. Such behavior makes it difficult to control the external load to achieve a target crack extension. Due to this limitation, only the first step of crack propagation during ENF testing was considered in the present study following ASTM D7905/D7905M [12].
The results from ENF testing on five FRP samples are depicted in Figure 2. The points in Figure  2a correspond to the load and the middle point deflection at which delamination started to propagate for each sample. Mode II fracture toughness, , (Figure 2b) was calculated based on the compliance method using the equation below [13]: where is the initial crack length, is the critical load at the onset of crack propagation, is the sample width, is the sample thickness, and longitudinal is the elastic modulus of the composite along the specimen length direction aligning with the fiber's orientation. The elastic properties of these samples in different directions have already been measured and reported in [14], where the modulus of elasticity and elongation were found to be 44058.68 ± 3460.23 MPa and 7.93% ± 1.40% along the longitudinal direction and 1814.72 ± 697.73 MPa and 0.08% ± 0.03% along the transverse direction, respectively.
The variation in the test results demonstrates that the material fracture properties change from one sample to another due to uneven distribution of fibers, pressure, heating, and cooling during manufacturing as well as human error. Such behavior of UD laminated composites embraces the necessity of considering the randomness of properties in simulation studies of these materials.

XFEM Modeling of Delamination
In many cases during the early stages of product development, the structure's dimensions and the test setup configurations make full-scale experimental evaluation cumbersome, increasing the demand for numerical analyses instead. A variety of numerical modeling techniques have been proposed in the past decades, and they can be categorized into mesh-free methods, such as smoothed particle hydrodynamics (SPH) [15], element-free Galerkin method (EFGM) [16], and finite difference method (FDM) [17], and mesh-based methods, such as finite element method (FEM) [18] and boundary element method (BEM) [19]. When it comes to modeling cracks, each method has its own advantages and disadvantages, and while the formulations of some techniques allow easier definition of cracks, e.g., SPH [20], others require modifications to their mathematical foundations so that analytical information about cracks can be included in the numerical space, e.g., FEM [21]. During ENF testing, crack propagation occurs due to excessive shear deformation, thus exhibiting mode II failure with an abrupt nature. Such behavior makes it difficult to control the external load to achieve a target crack extension. Due to this limitation, only the first step of crack propagation during ENF testing was considered in the present study following ASTM D7905/D7905M [12].
The results from ENF testing on five FRP samples are depicted in Figure 2. The points in Figure 2a correspond to the load and the middle point deflection at which delamination started to propagate for each sample. Mode II fracture toughness, G IIC , (Figure 2b) was calculated based on the compliance method using the equation below [13]: where a is the initial crack length, P is the critical load at the onset of crack propagation, w is the sample width, t is the sample thickness, and E longitudinal is the elastic modulus of the composite along the specimen length direction aligning with the fiber's orientation. The elastic properties of these samples in different directions have already been measured and reported in [14], where the modulus of elasticity and elongation were found to be 44058.68 ± 3460.23 MPa and 7.93% ± 1.40% along the longitudinal direction and 1814.72 ± 697.73 MPa and 0.08% ± 0.03% along the transverse direction, respectively.  Various modifications to FEM have been proposed to account for crack nucleation and propagation. In cohesive zone models [22,23], it is assumed that a fracture process zone exists ahead of the crack tip, which develops according to a traction-separation law. With this approach, the stress singularity at the crack tip is avoided; however, its major limitation is that it requires prior knowledge of the crack paths for incorporating cohesive elements. Phase-field models [24,25] implicitly track cracks in the computational domain by introducing an auxiliary scalar field variable to represent crack topology. A governing equation defines the interaction between displacement and auxiliary parameter fields to model crack initiation and propagation. Another method is the coupled criterion framework [26,27], which combines energy-based and stress-based conditions for crack nucleation to  The variation in the test results demonstrates that the material fracture properties change from one sample to another due to uneven distribution of fibers, pressure, heating, and cooling during manufacturing as well as human error. Such behavior of UD laminated composites embraces the necessity of considering the randomness of properties in simulation studies of these materials.

XFEM Modeling of Delamination
In many cases during the early stages of product development, the structure's dimensions and the test setup configurations make full-scale experimental evaluation cumbersome, increasing the demand for numerical analyses instead. A variety of numerical modeling techniques have been proposed in the past decades, and they can be categorized into mesh-free methods, such as smoothed particle hydrodynamics (SPH) [15], element-free Galerkin method (EFGM) [16], and finite difference method (FDM) [17], and mesh-based methods, such as finite element method (FEM) [18] and boundary element method (BEM) [19]. When it comes to modeling cracks, each method has its own advantages and disadvantages, and while the formulations of some techniques allow easier definition of cracks, e.g., SPH [20], others require modifications to their mathematical foundations so that analytical information about cracks can be included in the numerical space, e.g., FEM [21].
Various modifications to FEM have been proposed to account for crack nucleation and propagation. In cohesive zone models [22,23], it is assumed that a fracture process zone exists ahead of the crack tip, which develops according to a traction-separation law. With this approach, the stress singularity at the crack tip is avoided; however, its major limitation is that it requires prior knowledge of the crack paths for incorporating cohesive elements. Phase-field models [24,25] implicitly track cracks in the computational domain by introducing an auxiliary scalar field variable to represent crack topology. A governing equation defines the interaction between displacement and auxiliary parameter fields to model crack initiation and propagation. Another method is the coupled criterion framework [26,27], which combines energy-based and stress-based conditions for crack nucleation to improve the ability and efficiency of FEM in predicting crack initiation loading.
FEM can also be enhanced and used in modeling discontinuities by enriching its polynomial approximate functions, the so-called shape functions, with the partition of unity method (PUM) proposed by Melenk and Babuška [28]. This hybrid scheme is known as the extended finite element method (XFEM) [29,30], which offers a more accurate and efficient means to study the geometries and evolution of cracks in structures. In XFEM, similar to conventional FEM, the finite element mesh is generated regardless of the discontinuity locations. Then, specific search algorithms, such as the level-set or fast marching methods are utilized to identify the location of any discontinuity with respect to the existing mesh and apply enrichment on the affected elements. Next, additional auxiliary degrees of freedom are added to a group of nodes around the discontinuity. These degrees of freedom assist the model in capturing the displacement jumps caused by discontinuities.
Based on the XFEM modeling framework, we have previously developed a user-defined element in the finite element (FE) software Abaqus to simulate delamination in UD composites under mode I loading [31]. Here, in the next sections, we provide the fundamental formulations of our approach and demonstrate its efficacy to analyze mode II delamination failure.

Modeling Cracks in XFEM
Assume a discontinuity (a crack) within an arbitrary finite element mesh ( Figure 3). The displacement field of point X, u(X), inside the domain is described with two parts: one related to the conventional finite element approximation and the other associated with the XFEM enriched field defining the discontinuity [29]: where φ i (X) is the conventional shape function, ψ(X) is the general enrichment function, N all is the finite element mesh nodes, N f is the enriched nodes of the mesh, u ord i is the classic degrees of freedom at the node i, and u enr j is the additional enriched degree of freedom at the enriched node j.  In order to choose the enrichment function, i.e., ( ), any discontinuous function in the problem domain can be employed to estimate the displacement field approximation in the vicinity of the crack. A function that satisfies such a requirement is the Heaviside step function, ( ). It gains a value of +1 on one side of the crack and −1 on the other side and can be utilized when the crack propagation is modeled by multiple straight line segments. To find the Heaviside function value at each node of an element, tangential and normal vectors of the crack surface curve should be measured. If ( s e is the unit tangential vector, and is the out of plane vector), then using a scalar product between the distance vector of the element's nodes and the normal vector of the crack surface, the Heaviside function value can be calculated as follows:

Modeling Contact on Material Interfaces Using XFEM
During ENF testing, samples undergo large deformation and may experience plasticity ( Figure 1). We therefore adopted a technique proposed by Khoei et al. [32] that integrates large deformation and contact modeling into XFEM formulations. Here, we provide the fundamental formulations for this technique and invite the reader to refer to [32] for detailed development of the model.
For nonlinear problems, implementation of a finite element method results in a set of nonlinear algebraic equations, ( ) = , in which is the external load and is a nonlinear function of nodal displacements, (X). This set of nonlinear equations is solved through an incremental iterative technique in which the problem is divided into incremental load steps, and a linear set of equations In order to choose the enrichment function, i.e., ψ(X), any discontinuous function in the problem domain can be employed to estimate the displacement field approximation in the vicinity of the crack. A function that satisfies such a requirement is the Heaviside step function, H(X). It gains a value of +1 on one side of the crack and −1 on the other side and can be utilized when the crack propagation is modeled by multiple straight line segments. To find the Heaviside function value at each node of an element, tangential and normal vectors of the crack surface curve should be measured. If X * is the nearest point of a crack to the node X ( Figure 4) and e n is the unit normal vector of the crack at point X * where e s × e n = e z (e s is the unit tangential vector, and e z is the out of plane vector), then using a scalar product between the distance vector of the element's nodes and the normal vector of the crack surface, the Heaviside function value can be calculated as follows: Materials 2020, 13, x FOR PEER REVIEW 5 of 13 ( s e is the unit tangential vector, and is the out of plane vector), then using a scalar product between the distance vector of the element's nodes and the normal vector of the crack surface, the Heaviside function value can be calculated as follows:

Modeling Contact on Material Interfaces Using XFEM
During ENF testing, samples undergo large deformation and may experience plasticity (Figure 1). We therefore adopted a technique proposed by Khoei et al. [32] that integrates large deformation and contact modeling into XFEM formulations. Here, we provide the fundamental formulations for this technique and invite the reader to refer to [32] for detailed development of the model.
For nonlinear problems, implementation of a finite element method results in a set of nonlinear algebraic equations, ( ) = , in which is the external load and is a nonlinear function of nodal displacements, (X). This set of nonlinear equations is solved through an incremental iterative technique in which the problem is divided into incremental load steps, and a linear set of equations

Modeling Contact on Material Interfaces Using XFEM
During ENF testing, samples undergo large deformation and may experience plasticity (Figure 1). We therefore adopted a technique proposed by Khoei et al. [32] that integrates large deformation and contact modeling into XFEM formulations. Here, we provide the fundamental formulations for this technique and invite the reader to refer to [32] for detailed development of the model.
For nonlinear problems, implementation of a finite element method results in a set of nonlinear algebraic equations, R(u) = P, in which P is the external load and R is a nonlinear function of nodal displacements, u(X). This set of nonlinear equations is solved through an incremental iterative technique in which the problem is divided into incremental load steps, and a linear set of equations Materials 2020, 13, 3548 6 of 13 needs to be solved for each step, i.e., K T ∆u = ∆P, where K T is the tangential stiffness matrix. In our nonlinear XFEM model, K T is defined as follows [32]: where K Mat , K Geo , B, D ep s , G s , and M s are matrices associated with materials stiffness, geometrical stiffness, strain gradient, contact constitute behavior, Cartesian shape function derivatives, and rearranged second Piola-Kirchhoff stress, respectively. Γ V also denotes the element domain. In addition to the conventional FEM part of nodal vectors, B and G s also includes enriched elements as given below: where B u enr and G u enr s are defined as follows (H is the Heaviside function) [32]: where (X,Y,Z) and (x,y,z) refer to material and spatial coordinates, respectively. The contact constitutive matrix in Equation (4) is defined as follows [32]: where K ii is the penalty stiffness assigned to the local coordinates on the contact surfaces. K 11 provides the impenetrable characteristic to the normal direction of the crack plane, which follows the Kuhn-Tucker thresholds [33]: where δ n is the crack opening displacement, and P Contact contains the vector of contact forces. The remaining terms in the constitutive matrix, K 22 and K 33 , create the friction forces and prevent the contact surfaces from abrupt sliding. For these terms, standard static and dynamic friction laws can be applied to perform the analysis [32].

Cohesive Zone Implementation
As reported in previous studies [6,7], depending on the lay-up of the FRP composite, a large processing zone is expected to form during crack propagation. It is therefore necessary to incorporate cohesive crack modeling into XFEM formulations to analyze mode II delamination in FRP composites. Here, we have used a bilinear traction-separation law ( Figure 5) that relates traction on crack faces, T, to mode I and mode II nodal displacements. In the bilinear form, the traction at the interface increases linearly with the crack tip opening displacement (CTOD) to a limit value, T max . The interface element then experiences softening until a traction of zero is reached, which corresponds to complete debonding. The total area enclosed by the bilinear curve represents the fracture toughness of the material.

Cohesive Zone Implementation
As reported in previous studies [6,7], depending on the lay-up of the FRP composite, a large processing zone is expected to form during crack propagation. It is therefore necessary to incorporate cohesive crack modeling into XFEM formulations to analyze mode II delamination in FRP composites. Here, we have used a bilinear traction-separation law ( Figure 5) that relates traction on crack faces, , to mode I and mode II nodal displacements. In the bilinear form, the traction at the interface increases linearly with the crack tip opening displacement (CTOD) to a limit value, max . The interface element then experiences softening until a traction of zero is reached, which corresponds to complete debonding. The total area enclosed by the bilinear curve represents the fracture toughness of the material.
To model the cohesive behavior ahead of the crack tip, a cohesive transformation matrix ( Coh B ) is used to relate crack normal opening and sliding displacements, , to the displacement of a point at the crack interface, : In addition, the tractions on crack faces, , is related to as defined below: where Interface includes the cohesive interface material properties [6]. To model the cohesive behavior ahead of the crack tip, a cohesive transformation matrix (B Coh ) is used to relate crack normal opening and sliding displacements, δ, to the displacement of a point at the crack interface, u k : δ = B Coh u k (11) In addition, the tractions on crack faces, T, is related to δ as defined below: where D Interface includes the cohesive interface material properties [6]. The cohesive transformation matrix can be extracted by finding the displacement in an enriched element. The displacement vector of a point in the enriched element, u(X), is as follow: (13) where N enr The conventional finite element shape function's value remains constant for different points in the enriched element, while the enriched shape function's value demonstrates an odd function property with respect to the interface position: Thus, the global relative crack displacement, δ, can be described in the form of the displacement difference between two points above and beneath the crack surface: In order to find the relative crack displacements in a global coordinate system, a simple transformation based on the normal and tangential directions, m ij , of the crack plane with respect to the global coordinate can be employed: Consequently, Equation (18) can be substituted into Equation (12) and used in the tangential stiffness formulation to introduce process zone properties within enriched elements: Finally, in order to evaluate the internal forces, one can simply employ Equation (20) as follows: Similar to conventional application of cohesive zone models, a predefined crack path was utilized to model the cracking behavior. With regard to the parameters in the traction-separation law, T max was defined as 5.5 MPa based on the shear strength of PPS/glass composite samples [31] considering that the shear strength is 0.577 of tensile strength according to the maximum distortion energy theory. To identify an optimal value for k pen , a range of 10 3 N/mm 3 to 10 7 N/mm 3 was considered in trial simulations. It was observed that applying penalty stiffness lower than 10 4 N/mm 3 would result in extensive softening and a significant reduction in the peak load. Such behavior results from lower rigidity in the hardening region of the traction-separation law for the given material. On the other hand, excessive hardening was observed in simulations with k pen higher than 10 6 N/mm 3 , for which the deflection of the specimen reduced unrealistically and prevented the physical crack opening behavior from being modeled. A value of 10 5 N/mm 3 for penalty stiffness was eventually selected and used for subsequent simulations.

Stochastic Fracture Properties
Based on the experimental results given in Figure 2, we considered that the fracture properties of the tested UD samples had a stochastic nature as opposed to deterministic approaches where averaged values of experimental results are assumed for the estimation of material properties. The stochastic fracture properties of the materials were incorporated into our XFEM model through the procedure explained below.
A random number within the range of experimentally measured mode II fracture toughness values (G IIC ) was picked to form a stochastic bilinear traction-separation law for the enriched elements in the cohesive zone ( Figure 6): where Rand 1 is a random integer (odd or even to assign a random sign), and G IIC(ave) and G IIC(std) are average and standard deviations of fracture toughness, respectively. The second random number (Rand 2 ), corresponding to individual enriched elements in each stage of damage evolution, was taken from a uniform distribution with a range of 0 to 1. The randomly selected values were then converted to a Weibull two-parameter distribution between 0 and 1 via the following formula: where α 1 > 0 is a shape parameter, β 1 > 0 is the scale parameter of distribution, and both are considered to be equal to 3.
Materials 2020, 13, x FOR PEER REVIEW 9 of 13 A random number within the range of experimentally measured mode II fracture toughness values ( ) was picked to form a stochastic bilinear traction-separation law for the enriched elements in the cohesive zone ( Figure 6): where is a random integer (odd or even to assign a random sign), and ( ) and ( ) are average and standard deviations of fracture toughness, respectively. The second random number (Rand2), corresponding to individual enriched elements in each stage of damage evolution, was taken from a uniform distribution with a range of 0 to 1. The randomly selected values were then converted to a Weibull two-parameter distribution between 0 and 1 via the following formula: Rand Rand (22) where 1 α > 0 is a shape parameter, 1 β > 0 is the scale parameter of distribution, and both are considered to be equal to 3. Figure 6. Proposed stochastic bilinear traction-separation behavior (Rand2 is a random number taken from a two-parameter Weibull distribution; GIICL and GIICH correspond to the lower and upper limits of GIIC. It was assumed that the onset of crack propagation was mainly influenced by the resin itself, and given that the resin showed more deterministic material properties, Tmax was assumed to have no variation in the model. Because the fracture toughness distribution is considered to be a function of the crack length, a linear interpolation was utilized to extract GIIC(ave) for each specific crack length. For GIIC(std), it can be a constant or, in a more general form, it can be scaled with GIIC(ave), which in turn becomes a function of crack length. It should also be added that according to the bilinear traction-separation law, a direct relationship exists between the critical fracture toughness, GIIC, critical sliding displacement, f δ , and maximum interface shear strength, max T : Therefore, the obtained statistical distribution of the fracture toughness can be converted into the variation of failure crack sliding and/or maximum interface strength of material via Equation (23). Khokhar et al. [34] introduced randomness into their simulation by implementing a relationship between random fracture toughness and the maximum interface strength by keeping the failure crack sliding displacement constant. To improve the convergence of the numerical simulation, the present study assumed a constant value for the maximum interface shear strength ( max = 5.5 MPa), while the failure crack sliding displacement was randomly varied during damage evolution (see Figure 6). This approach relies on constant penalty stiffness and prevents the over-strengthening of the stiffness of Figure 6. Proposed stochastic bilinear traction-separation behavior (Rand 2 is a random number taken from a two-parameter Weibull distribution; G IICL and G IICH correspond to the lower and upper limits of G IIC . It was assumed that the onset of crack propagation was mainly influenced by the resin itself, and given that the resin showed more deterministic material properties, T max was assumed to have no variation in the model. Because the fracture toughness distribution is considered to be a function of the crack length, a linear interpolation was utilized to extract G IIC(ave) for each specific crack length. For G IIC(std) , it can be a constant or, in a more general form, it can be scaled with G IIC(ave) , which in turn becomes a function of crack length. It should also be added that according to the bilinear traction-separation law, a direct relationship exists between the critical fracture toughness, G IIC , critical sliding displacement, δ f , and maximum interface shear strength, T max : Therefore, the obtained statistical distribution of the fracture toughness can be converted into the variation of failure crack sliding and/or maximum interface strength of material via Equation (23). Khokhar et al. [34] introduced randomness into their simulation by implementing a relationship between random fracture toughness and the maximum interface strength by keeping the failure crack sliding displacement constant. To improve the convergence of the numerical simulation, the present study assumed a constant value for the maximum interface shear strength (T max = 5.5 MPa), while the failure crack sliding displacement was randomly varied during damage evolution (see Figure 6). This approach relies on constant penalty stiffness and prevents the over-strengthening of the stiffness of elements in the process zone. It also stands with the fact that the crack length extension in test specimens is a function of the CTOD and the energy release rate in front of the crack tip. Here, the process zone size was selected based on experimental results. Sensitivity analysis was performed with a range of process zones from 15 to 30 mm (length of approximately 15 to 30% of specimens) to confirm the assumption. A value of 25 mm for the cohesive zone length was used for stochastic simulations [35].

Results and Discussion
The above formulations were integrated into an Abaqus user-defined element subroutine, which is accessible in [14], coupled with MATLAB. The implicit solver of Abaqus was used for accurate evaluation of displacement field. Randomness was introduced into the analysis by means of the constant standard deviation method. In each simulation, a preassigned initial crack was considered in the specimen, and the stochastic fracture properties were assigned to the elements in the area near the crack front. Once fracture toughness value was assigned to a given element, it remained unchanged during the given simulation run. As the crack propagated, new elements would enter into the process zone and stochastic fracture properties would similarly be assigned to them.
The results for continuous delamination are given in Figure 7. The results of the stochastic simulations were in great agreement with the nonrepeatable experimental data obtained from ENF testing. The critical sliding displacement, δ f , was different in each simulation run, affecting the value of the maximum load at which the mid-plane crack started to propagate. However, the shape of the load-displacement curve was the same irrespective of the δ f value. Based on the trend in the simulation data, it can be concluded that the fiber bridging effect in ENF testing had minimal effect on the global characteristics of the response. Figure 8 shows the XFEM model contours under different stages of delamination during ENF testing. It should be noted that, from an application perspective, particularity in the context of composite-forming processes, the mode II loading, similar to ENF testing, would be more relevant due to the likelihood of sliding between layers of the laminate under the punch load.
Materials 2020, 13, x FOR PEER REVIEW 10 of 13 elements in the process zone. It also stands with the fact that the crack length extension in test specimens is a function of the CTOD and the energy release rate in front of the crack tip. Here, the process zone size was selected based on experimental results. Sensitivity analysis was performed with a range of process zones from 15 to 30 mm (length of approximately 15 to 30% of specimens) to confirm the assumption. A value of 25 mm for the cohesive zone length was used for stochastic simulations [35].

Results and Discussion
The above formulations were integrated into an Abaqus user-defined element subroutine, which is accessible in [14], coupled with MATLAB. The implicit solver of Abaqus was used for accurate evaluation of displacement field. Randomness was introduced into the analysis by means of the constant standard deviation method. In each simulation, a preassigned initial crack was considered in the specimen, and the stochastic fracture properties were assigned to the elements in the area near the crack front. Once fracture toughness value was assigned to a given element, it remained unchanged during the given simulation run. As the crack propagated, new elements would enter into the process zone and stochastic fracture properties would similarly be assigned to them.
The results for continuous delamination are given in Figure 7. The results of the stochastic simulations were in great agreement with the nonrepeatable experimental data obtained from ENF testing. The critical sliding displacement, f δ , was different in each simulation run, affecting the value of the maximum load at which the mid-plane crack started to propagate. However, the shape of the load-displacement curve was the same irrespective of the f δ value. Based on the trend in the simulation data, it can be concluded that the fiber bridging effect in ENF testing had minimal effect on the global characteristics of the response. Figure 8 shows the XFEM model contours under different stages of delamination during ENF testing. It should be noted that, from an application perspective, particularity in the context of composite-forming processes, the mode II loading, similar to ENF testing, would be more relevant due to the likelihood of sliding between layers of the laminate under the punch load.     Note that interface is considered here as an equivalent (macro) medium, dominated by the resin along with random fiber bridges; notice the stress concentration at the crack front.

Conclusion
In the present work, a framework is presented to numerically simulate the mode II fracture behavior of FRP composite materials. For this purpose, extending the capability of modeling delamination and contact interfaces in large deformation problems, a user-element subroutine was developed in Abaqus by incorporating nonlinear XFEM element properties with the cohesive zone model and contact formulation. In addition, the stochastic fracture properties of the composite samples were incorporated into the code to capture the randomness in properties evidenced in nonrepeatable ENF testing results of composite materials. The efficacy of the model was demonstrated by predicting the stochastic fracture behavior of PPS/glass UD laminates, which was in good agreement with the experimental data.   Note that interface is considered here as an equivalent (macro) medium, dominated by the resin along with random fiber bridges; notice the stress concentration at the crack front.

Conclusion
In the present work, a framework is presented to numerically simulate the mode II fracture behavior of FRP composite materials. For this purpose, extending the capability of modeling delamination and contact interfaces in large deformation problems, a user-element subroutine was developed in Abaqus by incorporating nonlinear XFEM element properties with the cohesive zone model and contact formulation. In addition, the stochastic fracture properties of the composite samples were incorporated into the code to capture the randomness in properties evidenced in nonrepeatable ENF testing results of composite materials. The efficacy of the model was demonstrated by predicting the stochastic fracture behavior of PPS/glass UD laminates, which was in good agreement with the experimental data.

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