Next Article in Journal
Advanced Sandwich Composite Cores for Patient Support in Advanced Clinical Imaging and Oncology Treatment
Next Article in Special Issue
The Elastic Share of Inelastic Stress–Strain Paths of Woven Fabrics
Previous Article in Journal
Daylight-Active Cellulose Nanocrystals Containing Anthraquinone Structures
Previous Article in Special Issue
Modelling Inhomogeneity of Veneer Laminates with a Finite Element Mapping Method Based on Arbitrary Grayscale Images
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
School of Engineering, University of British Columbia, Kelowna, BC V1V 1V7, Canada
2
Materials and Manufacturing Research Institute, University of British Columbia, Kelowna, BC V1V 1V7, Canada
*
Author to whom correspondence should be addressed.
Materials 2020, 13(16), 3548; https://doi.org/10.3390/ma13163548
Submission received: 2 July 2020 / Revised: 5 August 2020 / Accepted: 9 August 2020 / Published: 12 August 2020

Abstract

:
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.

1. 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 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.

2. 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).
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 I I C , (Figure 2b) was calculated based on the compliance method using the equation below [13]:
G I I C = 9 a 2 P 2 2 w 2 t 3 E longitudinal
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.
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.

3. 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.

3.1. 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]:
u ( X ) = i n i N a l l ϕ i ( X ) u i ord + j n j N f ϕ j ( X ) ψ ( X ) u j enr
where ϕ i ( X ) is the conventional shape function, ψ ( X ) is the general enrichment function, N a l l is the finite element mesh nodes, N f is the enriched nodes of the mesh, u i ord is the classic degrees of freedom at the node i, and u j enr is the additional enriched degree of freedom at the enriched node j .
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:
H ( X ) = { + 1 , i f   ( X X * ) . e n > 0 1 , o t h e r w i s e

3.2. 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 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]:
K T = K M a t + K G e o = Γ V B ¯ T D S e p B ¯ d Γ V + Γ V G s T M S G s d Γ V
where K M a t , K G e o , B ¯ , D s e p , 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:
B ¯ = [ B ¯ u ord B ¯ u enr ]
G s = [ G s u ord G s u enr ]
where B ¯ u enr and G s u enr are defined as follows (H is the Heaviside function) [32]:
B ¯ u enr = [ ( N i H ) X x X ( N i H ) X y X ( N i H ) X z X ( N i H ) Y x Y ( N i H ) Y y Y ( N i H ) Y z Y ( N i H ) Z x Z ( N i H ) Z y Z ( N i H ) Z z Z ( N i H ) Y x X + ( N i H ) X x Y ( N i H ) Y y X + ( N i H ) X y Y ( N i H ) Y z X + ( N i H ) X z Y ( N i H ) Z x Y + ( N i H ) Y x Z ( N i H ) Z y Y + ( N i H ) Y y Z ( N i H ) Z z Y + ( N i H ) Y z Z ( N i H ) Z x X + ( N i H ) X x Z ( N i H ) Z y X + ( N i H ) X y Z ( N i H ) Z z X + ( N i H ) X z Z ]
G s u enr = [ ( N i H ) X 0 0 0 ( N i H ) X 0 0 0 ( N i H ) X ( N i H ) Y 0 0 0 ( N i H ) Y 0 0 0 ( N i H ) Y ( N i H ) Z 0 0 0 ( N i H ) Z 0 0 0 ( N i H ) Z ]
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]:
D S e p = [ K ¯ 11 0 0 0 K ¯ 22 0 0 0 K ¯ 33 ]
where K ¯ i i 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]:
δ n 0 , P C o n t a c t 0 , ( δ n ) × ( P C o n t a c t ) = 0
where δ n is the crack opening displacement, and P C o n t a c t 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].

3.3. 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.
To model the cohesive behavior ahead of the crack tip, a cohesive transformation matrix ( B ¯ C o h ) is used to relate crack normal opening and sliding displacements, δ , to the displacement of a point at the crack interface, u k :
δ = B ¯ C o h   u k
In addition, the tractions on crack faces, T , is related to δ as defined below:
T = D ¯ Interface δ
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:
u ( X ) = i ( N i u i ord ) + j ( N j ( k N k H k H j ) u j enr ) = [ N 0 0 N enr ] [ u ord u enr ]
where
N j enr = N j ( k N k H k H j )
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:
N enr ( b o t t o m ) = N enr ( t o p )
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:
δ ¯ = u k ( t o p ) u k ( b o t t o m ) = [ N N enr ] [ u ord u enr ] [ N N enr ] [ u ord u enr ]
δ ¯ = 2 [ N N enr ] [ u ord u enr ]
In order to find the relative crack displacements in a global coordinate system, a simple transformation based on the normal and tangential directions, m i j , of the crack plane with respect to the global coordinate can be employed:
δ = [ m 11 m 12 m 13 m 21 m 22 m 23 m 31 m 32 m 33 ] [ δ ¯ X δ ¯ Y δ ¯ Z ] = 2 [ m 11 m 12 m 13 m 21 m 22 m 23 m 31 m 32 m 33 ] [ 0 N enr ] [ u ord u enr ] = B ¯ C o h   u k
Consequently, Equation (18) can be substituted into Equation (12) and used in the tangential stiffness formulation to introduce process zone properties within enriched elements:
K T = K M a t + K G e o + K C o h = Γ V B ¯ T D S e p B ¯ d Γ V + Γ V G T M S G d Γ V + Γ c ( B ¯ C o h ) T D ¯ I n t e r f a c e B ¯ C o h d Γ c
Finally, in order to evaluate the internal forces, one can simply employ Equation (20) as follows:
F int = Γ V B ¯ T σ d Γ V + Γ C N enr T f t d Γ C
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 103 N/mm3 to 107 N/mm3 was considered in trial simulations. It was observed that applying penalty stiffness lower than 104 N/mm3 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 106 N/mm3, for which the deflection of the specimen reduced unrealistically and prevented the physical crack opening behavior from being modeled. A value of 105 N/mm3 for penalty stiffness was eventually selected and used for subsequent simulations.

3.4. 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 I I C ) was picked to form a stochastic bilinear traction–separation law for the enriched elements in the cohesive zone (Figure 6):
G I I C = G I I C ( a v e ) + ( 1 ) R a n d 1 R a n d 2 × G I I C ( s t d )
where R a n d 1 is a random integer (odd or even to assign a random sign), and G I I C ( a v e ) and G I I C ( s t d )   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:
R a n d 2 w e i b u l l = [ 1 α 1 ln ( 1 R a n d 2 u n i f o r m ) ] 1 β 1
where α 1 > 0 is a shape parameter, β 1 > 0 is the scale parameter of distribution, and both are considered to be equal to 3.
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, T max :
G I I C = T max δ f 2
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].

4. 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.

5. 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.

Author Contributions

Conceptualization, A.S.M. and D.M.; methodology, A.S.M. and D.M.; software, D.M.; validation, D.M., A.S.M., and M.T.; data curation, D.M.; writing—original draft preparation, D.M.; writing—review and editing, A.S.M. and M.T.; visualization, D.M. and M.T.; supervision, A.S.M. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by The Natural Sciences and Engineering Research Council of Canada, Discovery Grant #2018-04575.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Prasad, M.S.; Venkatesha, C.S.; Jayaraju, T. Experimental methods of determining fracture toughness of fiber reinforced polymer composites under various loading conditions. J. Miner. Mater. Charact. Eng. 2011, 10, 1263–1275. [Google Scholar] [CrossRef]
  2. Gaggar, S.; Broutman, L.J. Fracture toughness of random glass fiber epoxy composites: An experimental investigation. Flaw Growth Fract. ASTM STP 1977, 631, 310–330. [Google Scholar]
  3. Mower, T.M.; Li, V.C. Fracture characterization of short fiber reinforced thermoset resin composites. Eng. Fract. Mech. 1987, 26, 593–603. [Google Scholar] [CrossRef]
  4. Russell, A.J. On the Measurement of Mode II Interlaminar Fracture Energies; Material Report 82; Defence Research Establishment Pacific: Victoria, BC, Canada, 1982. [Google Scholar]
  5. Murri, G.B.; Obrien, T. Interlaminar GIIC evaluation of toughened resin matrix composites using the end-notched flexure test. In Proceedings of the 26th AIAA/ASME/ASCE/AHS Structures, Structural Dynamics, and Materials Conference, Orlando, FL, USA, 15–17 April 1985. [Google Scholar]
  6. Harper, P.W.; Hallett, S.R. Cohesive zone length in numerical simulations of composite delamination. Eng. Fract. Mech. 2008, 75, 4774–4792. [Google Scholar] [CrossRef] [Green Version]
  7. Fan, C.; Jar, P.Y.B.; Cheng, J.R. Cohesive zone with continuum damage properties for simulation of delamination development in fibre composites and failure of adhesive joints. Eng. Fract. Mech. 2008, 75, 3866–3880. [Google Scholar] [CrossRef]
  8. Silberschmidt, V.V. Matrix cracking in cross-ply laminates: Effect of randomness. Compos. A Appl. Sci. Manuf. 2005, 36, 129–135. [Google Scholar] [CrossRef]
  9. Skordos, A.; Sutcliffe, M. Stochastic simulation of woven composites forming. Compos. Sci. Technol. 2008, 68, 283–296. [Google Scholar] [CrossRef]
  10. Ashcroft, I.A.; Khokar, Z.R.; Silberschmidt, V.V. Modelling the effect of micro-structural randomness on the mechanical response of composite laminates through the application of stochastic cohesive zone elements. Comput. Mater. Sci. 2012, 52, 95–100. [Google Scholar] [CrossRef]
  11. Jumel, J. Crack propagation along interface having randomly fluctuating mechanical properties during DCB test finite difference implementation–Evaluation of Gc distribution with effective crack length technique. Compos. Part B Eng. 2017, 116, 253–265. [Google Scholar] [CrossRef]
  12. ASTM Standard. D7905/D7905M–14. Standard Test Method for Determination of the Mode II Interlaminar Fracture Toughness of Unidirectional Fiber-Reinforced Polymer Matrix Composites; ASTM Standard: West Conshohocken, PA, USA, 2014. [Google Scholar]
  13. Park, S.J.; Seo, M.K. Modeling of Fiber–Matrix Interface in Composite Materials. In Interface Science and Technology; Elsevier: Amsterdam, The Netherlands, 2011; Volume 18, pp. 739–776. [Google Scholar]
  14. Motamedi, D. Nonlinear XFEM Modeling of Delamination in Fiber Reinforced Composites Considering Uncertain Fracture Properties and Effect of Fiber Bridging. Ph.D. Thesis, University of British Columbia, Vancouver, BC, Canada, 2013. [Google Scholar]
  15. Monaghan, J.J. Smoothed particle hydrodynamics. Rep. Prog. Phys. 2005, 68, 1703. [Google Scholar] [CrossRef]
  16. Belytschko, T.; Lu, Y.Y.; Gu, L. Element-free Galerkin methods. Int. J. Numer. Methods Eng. 1994, 37, 229–256. [Google Scholar] [CrossRef]
  17. Coates, R.T.; Schoenberg, M. Finite-difference modeling of faults and fractures. Geophysics 1995, 60, 1514–1526. [Google Scholar] [CrossRef] [Green Version]
  18. Zienkiewicz, O.C.; Taylor, R.L. The Finite Element Method: Solid Mechanics; Butterworth-Heinemann: Oxford, UK, 2000; Volume 2. [Google Scholar]
  19. Mi, Y.; Aliabadi, M.H. Dual boundary element method for three-dimensional fracture mechanics analysis. Eng. Anal. Bound. Elem. 1992, 10, 161–171. [Google Scholar] [CrossRef]
  20. Douillet-Grellier, T.; Jones, B.D.; Pramanik, R.; Pan, K.; Albaiz, A.; Williams, J.R. Mixed-mode fracture modeling with smoothed particle hydrodynamics. Comput. Geotech. 2016, 79, 73–85. [Google Scholar] [CrossRef]
  21. Liu, P.F.; Zheng, J.Y. Recent developments on damage modeling and finite element analysis for composite laminates: A review. Mater. Des. 2010, 31, 3825–3834. [Google Scholar] [CrossRef]
  22. Barenblatt, G.I. The formation of equilibrium cracks during brittle fracture. General ideas and hypotheses. Axially-symmetric cracks. J. Appl. Math. Mech. 1959, 23, 622–636. [Google Scholar] [CrossRef]
  23. Dugdale, D.S. Yielding of steel sheets containing slits. J. Mech. Phys. Solids 1960, 8, 100–104. [Google Scholar] [CrossRef]
  24. Miehe, C.; Welschinger, F.; Hofacker, M. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations. Int. J. Numer. Methods Eng. 2010, 83, 1273–1311. [Google Scholar] [CrossRef]
  25. Molnár, G.; Gravouil, A. 2D and 3D Abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture. Finite Elem. Anal. Des. 2017, 130, 27–38. [Google Scholar] [CrossRef] [Green Version]
  26. Leguillon, D. Strength or toughness? A criterion for crack onset at a notch. Eur. J. Mech. A Solids 2002, 21, 61–72. [Google Scholar] [CrossRef]
  27. Doitrand, A.; Martin, E.; Leguillon, D. Numerical implementation of the coupled criterion: Matched asymptotic and full finite element approaches. Finite Elem. Anal. Des. 2020, 168, 103344. [Google Scholar] [CrossRef]
  28. Melenk, J.M.; Babuška, I. The partition of unity finite element method: Basic theory and applications. Comput. Methods Appl. Mech. Eng. 1996, 139, 289–314. [Google Scholar] [CrossRef] [Green Version]
  29. Belytschko, T.; Black, T. Elastic crack growth in finite elements with minimal remeshing. Int. J. Numer. Methods Eng. 1999, 45, 601–620. [Google Scholar] [CrossRef]
  30. Moёs, N.; Dolbow, J.; Belytschko, T. A finite element method for crack growth without remeshing. Int. J. Numer. Methods Eng. 1999, 46, 131–150. [Google Scholar] [CrossRef]
  31. Motamedi, D.; Milani, A.S.; Komeili, M.; Bureau, M.N.; Thibault, F.; Trudel-Boucher, D. A stochastic XFEM model to study delamination in PPS/Glass UD composites: Effect of uncertain fracture properties. Appl. Compos. Mater. 2014, 21, 341–358. [Google Scholar] [CrossRef]
  32. Khoei, A.R.; Biabanaki, S.O.R.; Anahid, M. A Lagrangian extended finite element method in modeling large plasticity deformations and contact problems. Int. J. Mech. Sci. 2009, 51, 384–401. [Google Scholar] [CrossRef]
  33. Unger, J.F.; Eckardta, S.; Könke, C. Modelling of cohesive crack growth in concrete structures with the extended finite element method. Comput. Methods Appl. Mech. Eng. 2007, 196, 4087–4100. [Google Scholar] [CrossRef]
  34. Khokhar, Z.R.; Ashcroft, I.A.; Silberschmidt, V.V. Simulations of delamination in CFRP laminates: Effect of microstructural randomness. Comput. Mater. Sci. 2009, 46, 607–613. [Google Scholar] [CrossRef] [Green Version]
  35. Motamedi, D.; Milani, A.S. 3D nonlinear XFEM simulation of delamination in unidirectional composite laminates: A sensitivity analysis of modeling parameters. Open J. Compos. Mater. 2013, 3, 113–126. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Experimental setup of end-notched flexure (ENF) testing.
Figure 1. Experimental setup of end-notched flexure (ENF) testing.
Materials 13 03548 g001
Figure 2. ENF results with a constant crack length (43 mm). (a) Variation of critical load versus middle point deflection and (b) variation of fracture toughness versus middle point deflection.
Figure 2. ENF results with a constant crack length (43 mm). (a) Variation of critical load versus middle point deflection and (b) variation of fracture toughness versus middle point deflection.
Materials 13 03548 g002
Figure 3. An arbitrary finite element mesh with a discontinuity (circles represent the enriched nodes of the mesh).
Figure 3. An arbitrary finite element mesh with a discontinuity (circles represent the enriched nodes of the mesh).
Materials 13 03548 g003
Figure 4. Unit tangential and normal vectors for the Heaviside function and nearest point to X on the crack surface; X * .
Figure 4. Unit tangential and normal vectors for the Heaviside function and nearest point to X on the crack surface; X * .
Materials 13 03548 g004
Figure 5. Bilinear traction–separation law for modeling material degradation.
Figure 5. Bilinear traction–separation law for modeling material degradation.
Materials 13 03548 g005
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.
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.
Materials 13 03548 g006
Figure 7. Comparison of stochastic measured and predicted force–displacement values in ENF tests on the polyphenylene sulfide (PPS)/glass samples.
Figure 7. Comparison of stochastic measured and predicted force–displacement values in ENF tests on the polyphenylene sulfide (PPS)/glass samples.
Materials 13 03548 g007
Figure 8. Stages of delamination propagation within the ENF numerical model. (a) Stress concentration at onset of crack propagation. (b) Stress distribution during degradation in traction–separation law. 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.
Figure 8. Stages of delamination propagation within the ENF numerical model. (a) Stress concentration at onset of crack propagation. (b) Stress distribution during degradation in traction–separation law. 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.
Materials 13 03548 g008

Share and Cite

MDPI and ACS Style

Motamedi, D.; Takaffoli, M.; S. Milani, A. Nonlinear XFEM Modeling of Mode II Delamination in PPS/Glass Unidirectional Composites with Uncertain Fracture Properties. Materials 2020, 13, 3548. https://doi.org/10.3390/ma13163548

AMA Style

Motamedi D, Takaffoli M, S. Milani A. Nonlinear XFEM Modeling of Mode II Delamination in PPS/Glass Unidirectional Composites with Uncertain Fracture Properties. Materials. 2020; 13(16):3548. https://doi.org/10.3390/ma13163548

Chicago/Turabian Style

Motamedi, Damoon, Mahdi Takaffoli, and Abbas S. Milani. 2020. "Nonlinear XFEM Modeling of Mode II Delamination in PPS/Glass Unidirectional Composites with Uncertain Fracture Properties" Materials 13, no. 16: 3548. https://doi.org/10.3390/ma13163548

APA Style

Motamedi, D., Takaffoli, M., & S. Milani, A. (2020). Nonlinear XFEM Modeling of Mode II Delamination in PPS/Glass Unidirectional Composites with Uncertain Fracture Properties. Materials, 13(16), 3548. https://doi.org/10.3390/ma13163548

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop