Determination of the critical load and energy release rate in mode II delamination using a meshfree method

1 Introduction

Delamination is one of the most common types of damage in composites because of their relatively weak interlaminar strengths. Delamination may appear under various conditions such as transverse concentrated loads. This damage mode is especially important for the structural integrity of composite structures since it is difficult to detect during inspection. It has a detrimental influence on both stiffness and strength of structural components. In order to use the full capacity of a composite structure, it is important that the initiation of delamination cracks and the resulting reduction in laminate strength can be correctly modeled. Therefore, a significant research effort has been devoted to studying the laws governing the onset and evolution of delamination damage [6]. It is interesting to predict whether the delamination will grow as the structure is loaded.

Delamination of laminated composite structures has been widely examined, from both an experimental and a numerical standpoint, due to local failure and, sometimes, to sudden structural collapse. Therefore, knowledge of delamination modes and growth rate during the operating life of a composite structure is of primary concern. Several studies aimed at formulating a criterion for delamination onset and growth prediction have been published, and experimental testing was performed for validating theoretical models. An opening mode (mode I) and a shearing mode (mode II) are the main modes of fracture in composites. There are existing standard tests for determining the critical strain energy release rates associated with each mode. The criteria used to predict delamination propagation under mode loading conditions are usually established in terms of the components of the energy release rate (G_I, G_II, G_III) and fracture toughness (G_IC, G_IIC, G_IIIC). The strain energy release rate (G) is often appraised because of an experimentally measurable quantity and it has been mathematically well defined. Delamination propagation is predicted when the energy release rate is equal to the corresponding fracture toughness of the material. Therefore, several fracture test methods have been developed to measure the critical strain energy release rates for pure and mixed-mode loadings.

Numerical simulation of the delamination is based on either fracture mechanics or damage mechanics. Linear elastic fracture mechanics (LEFM) is used if the material nonlinearities can be neglected. Delamination growth of one or more cracks is predicted using LEFM. The assumption is that the delamination propagates when the associated energy release rate is greater than or equal to a critical value. The application of this method is easy since the crack front propagates in one dimension. When more than one crack propagates simultaneously, the mechanical behavior of the interface which includes the cohesive zone model is modeled on the basis of damage mechanics. The interface-traction relationship is based on a simple bilinear one-dimensional law. The area of the triangle delimited by the traction/relative displacement curve is equated to the critical energy release rate [1].

Computational methods such as the finite differences method (FDM), the finite element method (FEM) or the finite volume method (FVM) all make use of a mesh, i.e., a set of nodes with node-to-node connectivity, which has gained much attention over the last 10 years. Since these methods have been successfully used to solve engineering problems for decades, some limitations of these kinds of methods are becoming more and more obvious: mesh generation is very time consuming and, in problems with moving boundaries, a drastic loss of accuracy is observed due to the distortion of the elements. Therefore there is a need to develop meshfree (grid free, mesh less) methods that do not suffer from the problems described above. Meshfree methods and application of composite fracture continue to be used in in-depth research. Meshfree methods are very widely used for modeling problems with moving boundaries such as crack growth in solids, because changing discontinuities can be displayed very easily by modifying only the weighting functions, which is a well-known advantage compared to the classical finite element method. In the meshfree method, discontinuities are not restricted to element boundaries or even no elements are necessary at all. Only nodes are used for approximation [3, 9, 14].

Meshfree methods have existed for a number of years, and the element-free Galerkin (EFGM) was first proposed for the simulation of fracture 15 years ago by Belytschko et al. [4, 5]. Several review articles [2, 7, 8, 10, 11], which cover the classification, application and computer implementation aspects of the meshfree method, have been reported in the literature [12–22].

The aim of this study was to present a meshfree method for the simulation of the delamination and failure of composite materials based on the cohesive segments model. With the use of the partition of unity of moving least-squares (MLS) shape functions, the discontinuities at the cohesive segments can be approximated with the additional degrees of freedom of nodes. An iterative solution scheme between the continuous and discontinuous fields is presented to solve the mode II delamination growth problem. An experimental investigation was used to verify the meshfree method’s results. The experimental study used end-notched flexure (ENF) specimens made of carbon/epoxy laminate (AS4/3501-6) which consists of 10 plies in [0]₁₀ and [0/90/0/90/0]_s lay-up with delamination inserted in the middle of the laminate.

The problem was solved for [0]₁₀, [0/45/-45/90/0]_s, [0/90/0/90/0]_s, [0/90/0/90/30]_s, [0/90/0/90/45]_s and [0/90/0/90/60]_s laminates with mid-plane delaminations, and the results are verified for different composite materials. The critical fracture force, which can be experimentally measured, was used to calculate the mode II delamination fracture toughness for the carbon/epoxy laminate. In addition, values of the integral for 219 (11×19) and 253 (11×23) background meshes with equivalent interval sizes were compared.

2 The meshfree method

A number of different meshfree methods have been developed. One of the most widely used is the EFGM. This method implements the connectivity between the nodes and sets completely using approximation functions. The EFGM uses the moving least-squares (MLS) approximants in order to obtain test and trial functions; background cells are utilized for numerical integration in the weak form; and essential boundary conditions are implemented using Lagrange multipliers. The weak form of Galerkin was used to develop the discretized system of equations.

The meshfree method is described in a broad sense as a method where nodes are not required to be interconnected. A distribution of nodes is generated to define the domain and its boundaries, and the solution of the problem is sought at these nodes in the EFGM. However, cells are also prescribed over the domain and used for numerical integration. In the EFGM, the moving least squares (MLS) scheme is utilized for interpolation, and a mesh of background cells, with no required connection to the nodal discretization, is used for the purpose of numerical integration [11].

3 Cohesive zone model approach

The cohesive zone model approach is one of the most commonly used tools to investigate interfacial fracture. This approach gives the physical explanations of the failure process and it can be used for both damage tolerance and strength analysis. The cohesive zone model approach is defined by a bilinear constitutive law. In the cohesive zone model, relative displacements and stresses are defined as failure and they do not contain crack tip stress singularities. The fracture process region chosen is a narrow band of vanishing thickness ahead of a crack tip. Bonding of the surfaces of the region is caused by cohesive traction which follows a cohesive constitutive law. There are an initial high positive stiffness and a negative tangent stiffness representing softening. The constitutive law is a bilinear softening model which is regulated by both material strength parameters and fracture mechanics parameters. If the cohesive traction disappears, crack growth takes place. A critical strain energy release rate is equal to the amount of work done per unit area of crack surface.

Consider a domain Ω containing a material discontinuity Г_c, which divides the domain Ω into two parts: Ω₁ and, under crack, Ω₂. The stress field inside the domain σ is related to the external loading and to the closing tractions in the material discontinuity through the equilibrium equations. The solid can be loaded by body force b and surface traction (force) t in any distributed fashion in the volume of the solid and the traction τ(v) in the boundary Г_c. The traction τ depends on the displacement jump v on the cohesive segment Г_c.

The bilinear law is the most commonly utilized cohesive law because of its simplicity. The bilinear constitutive relationship is shown in Figure 2 for the softening traction-displacement curves of the interface material model for the cases of mode II delamination. The fracture toughness G_C, cohesive surface tractions τ and the cracking relative displacements v are input material properties. The onset displacement jump and final displacement jump are defined as v₀ and v_f, respectively. The area under the traction-displacement jump relation is defined as the fracture toughness or the critical fracture energy.

The bilinear stress relative displacement curve is composed of three main parts: elastic part, softening part and decohesion part. The constitutive equations are as follows:

where K₀ is the interface stiffness, D is a damage parameter whose initial value is 0. D starts growing when v>v₀ and reaches the value when v>v_f. The damage parameter can be computed with the following relations for the elastic part, softening part and decohesion part, respectively:

If the traction across the interface increases to the maximum, this part is defined as the elastic part. In the softening part, the traction across the interface decreases to zero. Also, two layers begin to separate from each other in the softening part. If there is no bond between the two layers, the traction across the interface is null. This case is defined as the decohesion part.

Damage initiation is related to the interfacial strength τ_m, i.e., the maximum traction on the traction-displacement jump relation. When v₀ is reached, the stress is equal to the interfacial tensile strength τ_m, the maximum stress level possible. For higher relative displacements, the interface accumulates damage and its ability to sustain stress decreases progressively. Once v exceeds v_f, the interface is fully debonded and it is no longer able to support any stress. If the load were removed after v₀ has been exceeded but before v_f has been reached, the model would unload to the origin. The slope of the constitutive equation before damage initiation K₀ is referred to as the interface stiffness [3, 19, 20].

The criteria used to predict delamination propagation are usually established in terms of the components of the energy release rate and fracture toughness. It is supposed that when the energy release rate G exceeds the critical value, delamination of the critical energy release rate G_C grows.

As illustrated in Figure 1, the shaded area under the τ-v curve is the energy dissipated per unit area; it is defined as the strain energy release rate and the two following relations exist among these parameters:

Figure 1

Bilinear interface traction for mode II delamination.

As shown in Figures 1 and 2, v30 and v3f are expressed as the displacement jumps corresponding to delamination onset in mode II and in shear mode, respectively, while N30 and N3f are defined as the mode interlaminar strengths. The following relations exist:

Figure 2

Cohesive zone ahead of delamination tip.

The opening values v_i are determined by equating the area under the curve to the critical fracture energies such that

The work done per unit area upon complete interface degradation is given by the integration of cohesive tractions as functions of the displacement jump,

where G_IIC, τ₁₂ and v3f are the critical strain energy release rate, the mode traction and the final relative displacements, respectively, corresponding to complete decohesion. A failure criterion for mode II delamination can be written when the following condition is fulfilled:

The existing literature [2] expression utilizes the interface as a resin-rich zone of small thickness e_i and has proposed to define stiffness as:

where it is assumed that E₃=E₂ and G₁₃=G₁₂. E₂ and E₃ are Young’s modulus in the transverse and thickness direction, respectively.

The stress-relative displacement relationship and the displacement jump can be formulated, respectively, as [3]:

where Φ is the shape function belonging to sub-domain Ω₁ and Ω₂:

4 Variational principle

The EFGM is formulated in the space coordinates x=[x, y] on the domain Ω bounded by Г. The equilibrium equation can be written at a material point x as follows:

where σ is the stress tensor, which corresponds to the displacement field u, and b is the body force vector. The boundary conditions are written as σ·n_c=τ(v) on Γ_c, σ⋅nt=t¯ on Γt and u=u¯ on Γu, where t¯ and u¯ are defined as traction on the traction (natural) boundaries and as displacement on the displacement (essential) boundaries; n_t is the normal unity vector of the boundary Г_t where the traction t¯ is prescribed; n_c is the normal unity vector of the boundary Г_c where the traction τ(v) is prescribed; and Г_u is the boundary where the displacement u¯ is imposed.

The boundary conditions in the EFGM cannot be performed as in classical computational techniques because the properties of Kronecker delta are not satisfied. Therefore special methods are required such as the Lagrange multipliers method, coupling with FEM, and penalty methods. In this study, Lagrange multipliers were used. The weak form is written with Lagrange multipliers, with essential boundary conditions [3]:

By taking a sufficiently high order of quadrature, the weak form is integrated numerically over the domain, as a function of the size of the background cells and of the nodal distribution.

Λ is a vector that collects the nodal Lagrange multipliers for all field nodes on essential boundaries:

If we rewrite Eq. (24), by defining K and F as the global stiffness matrix and the global force vector, respectively, we obtain

where K is called the nodal stiffness matrix. The terms in the stiffness matrix are defined by the following equation [3]:

in the delaminated part:

The global force F is

and

Since both δU and δΛ are arbitrary, this equation can be satisfied only if [3]

If it is rewritten, Eq. (32) becomes

then the solution to the whole problem can be obtained by the iterative solution of Eq. (33).

In the conventional EFGM, the stiffness matrix and load factor are computed by the Gaussian quadrature method with a background mesh.

Q and G are defined as a global vector form and a global matrix form, respectively:

To illustrate the applicability of the meshfree method, the delamination under mode II loading is simulated.

5 Numerical implementation: delamination of composite specimens

In this section, the validity and the accuracy of the meshfree EFGM are shown by determining the critical load and energy release rate in mode II delamination. To demonstrate the usefulness of the proposed method, the delamination problem is numerically analyzed by the EFG method. The fracture modes studied are mode II with an ENF test.

For the analysis of the critical load and energy release rate, some codes were developed in MATLAB. Four-point Gaussian quadrature method is employed in the EFG method. Linear basis functions in two dimensions are used. In addition, the cubic spline is taken as a weight function.

Load-energy release rate diagrams obtained from the meshfree method are shown in Figures 3–8 for different composite materials: e-glass/epoxy, s-glass/epoxy, Kevlar/epoxy (Aramid149/epoxy), carbon/epoxy (AS4/3501-6), carbon/PEEK (AS4/APC2) and carbon/epoxy (IM6/SC1081), with different orientation angles [13].

Figure 3

Energy release rate-load diagram for e-glass/epoxy obtained from meshfree method.

Figure 4

Energy release rate-load diagram for s-glass/epoxy obtained from meshfree method.

Figure 5

Energy release rate-load diagram for Kevlar/epoxy (Aramid149/epoxy) obtained from meshfree method.

Figure 6

Energy release rate-load diagram for carbon/epoxy (AS4/3501-6) obtained from meshfree method.

Figure 7

Energy release rate-load diagram for carbon/PEEK (AS4/APC2) obtained from meshfree method.

Figure 8

Energy release rate-load diagram for carbon/epoxy (IM6/SC1081) obtained from meshfree method.

Comparison of meshfree method results of the carbon/epoxy laminate (AS4/3501-6) with [0°]₁₀ and [0°/90°/0°/90°/0°]_s angle plies are shown in Table 1 for values of the integral for 11×19 and 11×23 background meshes.

6 Comparisons with experimental studies

To demonstrate the usefulness of the proposed method, delamination problem was investigated by analyzing the mode II delamination test. The ENF test was used to measure mode II delamination toughness. The most commonly used of these is the ENF test as it may be easily performed on thin specimens. The load measured at the onset of delamination from the precrack was used to evaluate the critical energy release rate. Tests were conducted using the ENF configuration (Figure 9), in which a specimen with an initial mid-plane delamination at one end is loaded in a three-point bending arrangement. Although there is no universally accepted standard for mode II testing, the ENF is commonly used. End-notched flexural specimens, end-loaded split (ELS) or short-length specimens are possible test configurations to quantify delamination in mode II [22].

Figure 9

Specimen for end-notched flexure (ENF) specimen.

The ENF specimen was fabricated with carbon/epoxy (AS4/3501–6), a high-strength carbon pre-impregnated with epoxy resin-based repairing materials. The delamination cracks were simulated by Teflon films. The specimen was 100 mm long, 40 mm wide and 2 mm thick, and it had an initial crack length of 30 mm (Figure 10). Simulation of the ENF specimen was performed for two different laminate lay-ups: Type 1 [0]₁₀ and Type 2 [0°/90°/0°/90°/0°]_s.

Figure 10

Specimen dimension.

The single-ply thickness was 0.2 mm. The mechanical properties of the unidirectional lamina are reported as follows: E₁₁=146.7 GPa, E₂₂=7.8 GPa, G₁₂=G₁₃=4.8 GPa and v₁₂=0.29 [23].

Fracture testing was performed using universal testing machine (Shimadzu, Autograph AG-X 100 series) with a 100 kN load cell. A test sample as shown in Figure 11 was loaded with a speed of 10 mm/min. The incremental test method for mode II fracture testing which requires loading-crack propagation-unloading sequences was followed. The specimens were loaded until crack propagation occurred.

Figure 11

End-notched flexure (ENF) specimen for mode II testing mounted in the fixture.

Typical load versus deflection curves (P-δ) obtained from testing the carbon/epoxy (AS4/3501-6) specimens using the ENF test are shown in Figures 12 and 13. As shown in Figure 13, when the displacement reaches about 3.4 mm, the discontinuity begins to propagate along the horizontal straight line. As the cohesive traction decreases with the increasing separation, the load-displacement curve rapidly declines as the crack propagates.

Figure 12

Load-displacement curve for [0°]₁₀ carbon/epoxy laminate in ENF test.

Figure 13

Load-displacement curve for [0°/90°/0°/90°/0°]_s carbon/epoxy laminate in ENF test.

In the present study, the existing literature expression was utilized to calculate the critical strain-energy release rate for the composites [7]:

In Eq. (36), P, a and h are defined as the applied load, the crack length and the thickness of each beam in the delaminated region, respectively. E₁ and G₃₁ are the equivalent elastic flexural modulus and the equivalent transverse shear modulus of the laminate, respectively. To calculate G_IIC from Eq. (36), first the critical load must be determined from the experiment. The critical load and energy release rate are given in Table 2 for carbon/epoxy laminate (AS4/3501–6).

The percentage relative errors of the critical load and energy release rate are defined as follows:

where E and EFG are the critical load obtained from the experiment and the meshfree EFG method, respectively.

7 Results and discussion

In this paper, critical load and energy release rate for mode II testing of a delaminated composite for various orientation angles and materials were obtained and compared. In the experimental method, the ENF specimen was used. The obtained results from experimental method were supported by the meshfree and finite element numerical methods. A comparison of the critical loads and energy release rate is given in Tables 3 and 4.

In Tables 3 and 4, the relative error percentages of the experimental data and the meshfree EFG are evaluated using both [0]₁₀ and [0/90/0/90/0]_s lay-up. The relative error percentages of the critical load and energy release rate are 5.86 and 12.72 for carbon/epoxy (AS4/3501-6) with the [0]₁₀ lay-up. Similarly, it is seen that the relative error percentages of the critical load and energy release rate are 7.74 and 16.69 with the [0/90/0/90/0]_s lay-up. The meshfree EFG results for the [0]₁₀ lay-up show good convergence and agreement with experimental data. The meshfree EFG method for the [0]₁₀ lay-up is computationally more effective than that for the [0/90/0/90/0]_s lay-up.

In addition to study the effects of background meshes, the critical load and the energy release rate for two different background meshes with equivalent interval sizes were calculated in the meshfree EFGM. The critical load and the energy release rate are presented using angle-ply [0]₁₀, [0/45/-45/90/0]_s, [0/90/0/90/0]_s, [0/90/0/90/30]_s, [0/90/0/90/45]_s and [0/90/0/90/60]_s laminates with mid-plane delamination with six different composite materials. A comparison of the critical loads and the energy release rate for values of the integral for 11×19 and 11×23 background meshes is given in Tables 5–12.

As can be easily observed, the maximum critical load was obtained for e-glass/epoxy with the [0/90/0/90/60]_s angle ply for 11×19 and 11×23 background meshes. The critical buckling loads with 11×19 and 11×23 background meshes were 1190 and 1215 N, respectively, for the [0/90/0/90/60]_s angle ply. A similar study of mesh size effect was repeated for the critical energy release rate. The critical energy release rate calculated using a mesh size of 11×19 and 11×23 was 2902.05 and 3117.57 N/m, respectively, for e-glass/epoxy with the [0/90/0/90/60]_s angle ply. Tables 5–8 show that the maximum critical load and the energy release rate predicted by the meshfree EFGM converge when increasing the background meshes.

As shown in Table 9, it is seen that D/A₃₂ affects the critical load and energy release rate values directly. The maximum and the minimum critical loads and energy release rates were obtained for the [0/90/0/90/60]_s and [0/45/-45/90/0]_s orientation angles, respectively.

8 Conclusions

In this paper, both numerical and experimental investigations are presented for the delamination under mode II loading, detailing the derivation of the formulations in numerical simulations of fracture. Based on the results of this study, the following conclusions can be drawn:

• It has been shown that the mode II delamination behavior of fiber-reinforced laminated composite beams could be analyzed by using a meshfree EFG method. Even though the node distribution used in this study was selected as uniform, node scattering could be changed to take into account high stress variation around the crack tip. This will help to improve the accuracy of the analysis.

• The experimental results and closed form solution were used together to calculate the mode II energy release rate for only one composite material having two different lay-ups. The critical load and energy release rate results of meshfree EFG and the FEM were compared with each other as a verification of the EFGM. Results obtained from the meshfree method showed very good agreement with experimental evidences for single-mode delamination under mode II loading. Therefore, the meshfree method is capable of accurately predicting delamination onset.

• In the meshfree method, critical loads and energy release rates are increased by an increase in the number of nodes. That is because the node density in the analysis increased everywhere including in the crack tip region. A relatively fine background mesh helps the analysis to estimate the critical load and energy release rate more accurately. Although the uniform mesh density was used in the analysis, for better performance varying the node density can be used such that more nodes may be placed around the crack tip while fewer nodes are used in the other regions.

• The maximum and the minimum critical loads and energy release rates were obtained for the [0/90/0/90/60]_s and [0/45/-45/90/0]_s orientation angles, respectively (Table 13). Since mode II energy release rate was calculated by using Eq. (17), where out-of-plane stress is given by

  (38)τγ=DA32-1Ny (38)

• the value of the expression D/A₃₂ affects the critical load and energy release rate values directly.

• Results obtained from the meshfree method showed very good agreement with experimental data for single-mode delamination under mode II loading. The results presented will help in the implementation of numerical analysis techniques that protect numerical accuracy while minimizing computational expense. Consequently, the meshfree method could be used effectively to study delamination growth in composite laminates and it is especially suitable for the simulation of complex delamination patterns that are difficult to model using other numerical methods.