Characterization of Delamination Crack in Multidirectional E-glass/epoxy Composite under Mode I Loading

In this study, mode I critical strain energy release rate (GIc) of unidirectional E-glass/epoxy was determined by double cantilever beam (DCB) test. Calculated GIc was used to initiate delamination in numerical models. The effect of stacking sequence and fiber orientation of the sublaminates and the effect of thickness on mode I delamination crack length and GIc distribution along the specimen width in 0 o//0o interface have been studied using ANSYS®. 3-D 8-node linear interface element INTER 205 is used to create a predefined crack path. To investigate the fiber orientation effect, composites with [+θ2, 902, -θ2, 02]s stacking sequences were modeled. Interlaminar fracture analyses were performed by Virtual Crack Closure Technique (VCCT). Experimental and numerical critical loads (Pcr) showed a good agreement. According to the results, while fiber orientation affects significantly the extended crack length and the GIc distribution along the specimen width, stacking sequence only affects the GIc distribution.


INTRODUCTION
Fiber-reinforced laminated composite materials are widely used in the aerospace industry due to their great potential in weight saving, high specific modulus and high specific strength.There withal the laminated composite materials made of brittle matrices are susceptible to interlaminar cracks (interlaminar mode of fracture -delamination) and to propagation of that cracks also.Especially low-velocity impact damage and micro-cracks formed during manufacturing, service or maintenance cause to delamination in laminated composite materials.
It is well known that the laminated composites generally have no reinforcement through the thickness [1,2] and the mechanical properties of fibers are superior to mechanical properties of matrix.On the other hand resin rich zones are excessively take place between the layers and so the delamination propagation is mostly matrix dependent interlaminar fracture phenomenon.In the cases of interlaminar shearing fracture and sliding fracture (mode II & III), frictional forces between the fractured surfaces affect the fracture toughness.Because of opening fracture (mode I) phenomenon is the friction free, interlaminar fracture toughness in pure mode I is slightly less than the other two mode and more critical.In pure mode I, interlaminar fracture toughness can be measured by the strain energy release rate (SERR) [3].While interlaminar fracture initiation is defined by critical SERR (G Ic ), delamination propagation is defined by SERR (G I ).To measure G I and G Ic of polymer matrix composites, DCB test is described by American Society of Testing Materials standard (ASTM D5528-01) [4].DCB specimen that has a precrack located at the middle interface, is subjected to two equal forces through reversed directions at the ends of the specimen [4].While applying load on DCB specimen with a rate of constant displacement to calculate the fracture energy release rate, it is concurrently needed to measure the interlaminar crack extension.
Several authors have focused on and conducted numerical and experimental investigations to define the influence of the stacking sequence on the interlaminar fracture initiation and the propagation by means of delamination resistance ( R-curve behavior) and G Ic distribution along the delamination tip [2,[5][6][7][8][9].Davidson et al. [8] represented that the energy release rate distribution along the delamination tip is a function of a non-dimensional ratio, D c (=D 2  12 / D 11 D 22 ), and bending -twisting coupling measurement B t (=|D 16 |/D 11 ) whereas D ij are bending stiffness's.To obtain symmetric strain energy release rate (SERR) distribution along the delamination tip, D c ratio should be less than 0.25 for DCB specimens with ±30 o //±30 o crack interfaces [8] and the value of B t as low as possible [7][8][9].
Shokrieh and Heidari-Rarani [5] studied the effect of stacking sequence on mode I delamination resistance of polymer matrix DCB specimen with 0 o //0 o crack interface and they reported that non-dimensional ratio D c may affect R-curve behavior of DCB specimens.Shetty et al. [10] concluded that the delamination propagation in 90 o //90 o interface exhibits highest fracture resistance, since the propagations in +45 o //-45 o and +60 o //-60 o interfaces are lower and similar to each other.Because of that the fiber bridging excessively takes place in 0 o //θ o interfaces, more uniform delamination propagation and lower G Ic values are conversely observed in 0 o //0 o interface of carbon-epoxy composite [6,7].
To simulate delamination propagation with Finite Element Method (FEM), VCCT and Cohesive Zone Method (CZM) are used.The VCCT is quite appropriate to analyze crack propagations in laminated composite materials with brittle matrix.The VCCT is gaining a great acceptance during the last decade as a Linear Elastic Fracture Mechanic (LEFM) based numerical method to characterize of mode I delamination growth [11,12].Zhang and Chryssanthopoulos [13] reported that VCCT can be used to simulate mode I delamination growth even though the technique exhibits significantly overestimated critical strain energy release rate.Bonhomme et al. [14] investigated mode I interlaminar fracture toughness of carbon/epoxy composite by using a two-step numerical method similar to the VCCT.Bonhomme et al. [15] also used their two-step numerical method to investigate mode I, mode II and mixed mode I & II interlaminar fracture toughness of carbon/epoxy laminate.They reported that a good agreement found between the experimental and the numerical results with a difference less than 10%.Beuth [16] presented a modified method by aid of VCCT to extract crack extension independent strain energy release rate calculation.
Previous researches are excessively dealed with the crack propagation and uniformity of the crack propagation in 0 o //θ o and +θ 0 //-θ 0 interfaces of the anti-symmetric carbon/epoxy composite laminates [6][7][8][9]11].But the delamination propagates more easily between 0 o //0 o interfaces than +θ 0 //-θ 0 interfaces.Just a few researchers studied on G Ic distribution along width and delamination propagation in 0 0 //0 0 interfaces of unidirectional laminated composites [5,9,10].But there is no study focused on the effect of multidirectional sublaminates to G Ic distribution along width and delamination propagation in 0 0 //0 0 interface.For this reason, the aim of this research is to investigate numerically effects of the fiber orientations and the stacking sequence on the overall delamination propagation and energy release rate distribution along the width in 0 o //0 o interface of E-glass/epoxy laminates with balanced [+θ 2 ,90 2 ,-θ 2 ,0 2 ] s stacking sequence.Initially, interlaminar fracture toughness of unidirectional E-glass/epoxy composite laminate was experimentally determined as a material property and then VCCT was used to perform numerical investigations on multidirectional models by employing the experimental result.

Manufacture of Specimens
The laminates with [0 o ] 16 stacking sequence are manufactured from 330 gr/m 2 unidirectional E-glass fabric and epoxy resin by aid of vacuum resin infusion transfer molding (VRITM) method.The fabric consists of 300 gr/m 2 E-glass weft, 27 gr/m 2 E-glass warp and 3 gr/m 2 polypropylene stitching yarns with the density of 10 stitches/inch.The resin system is prepared by mixture of bisphenol A based epoxy resin Hexion MGS L160 and a resin hardener Hexion MGS H160 in the mass ratio of 100:25.Fiber volume fraction (V f ) of the manufactured laminate is 47% and the density is 1.87 gr/cm 3 .Obtained laminate thickness is 4 mm.Experimentally determined mechanical properties of the laminate are listed in Table 1.While fabric stacking, mold release agent coated PET film with 12 µm thickness is located into the mid-plane to obtain a 50 mm precrack.Edge of the PET film is marked with a permanent marker before inserting into the mid-plane to detect the location of the precrack tip easily.The laminates are cured at the room temperature for 24 h under 1 atm pressure and followed by post-curing at 50 o C for 6 h.DCB specimens are cut accurately by computer numeric controlled water jet cutting machine in dimensions (width: w = 25 mm, length: L = 150 mm) like as suggested at ASTM D5528-01 [4].Edges of the specimens are coated with a thin layer of water-based typewriter correction fluid and then a ruler is marked on the coated edge starting from the precrack tip.So it is easy to detect the delamination growth initiation and to measure the delamination propagation.

Data Reduction
Due to simple beam theory, the energy release rate of a DCB specimen can be calculated by the following equation [4]: If the critical load (P cr ) and the critical displacement (δ cr ) values that corresponding to the crack initiation are inserted into Eq.1 used to determine SERR then G Ic can be calculated.However, since the arms of the DCB specimen are not perfectly built-in and the rotation may occur at the delamination front [17], a correction factor must be added into Eq.1 to prevent overestimated G Ic values.Three data reduction methods including three different correction factors (Δ, n, A 1 ) are proposed by ASTM D5528-01 [4] standard to determine G I and G Ic .These are: a) The Modified Beam Theory (MBT) Method [4]: ( ) b) The Compliance Calibration (CC) Method [4]: c) The Modified Compliance Calibration (MCC) Method [4]: C is the compliance of DCB specimen where calculated as δ/P.In this study, all of the three methods are used to determine G I and G Ic values.
Similarly, three methods are proposed in ASTM D5528-01 [4] to detect the point specified by P cr and δ cr .These are [4]; (1) The point of deviation from linearity (NL point) in the load-displacement curve (P-δ curve), (2) The point at which delamination onset is visually detected (VIS point) and (3) The point at which the compliance has increased by 5% (5%/max point).
Detailed knowledge about these methods is available in ASTM 5528-01 [4].

Virtual Crack Closure Technique
The method VCCT assumes that the energy needed to create the fractured surfaces is equal to the energy needed to close the fractured surfaces and the crack propagates just along a predefined path.Since the mechanical properties of laminae are superior to matrix, crack branching is hindered and so the delamination propagation is bordered into resin rich interface.As a result, VCCT is successfully applied to delamination fracture mechanics of the laminated composite materials.
The VCCT technique is based on two essential assumptions; • If there is a small amount of crack propagation, the released strain energy is equal to the energy needed to close the fractured surfaces to its original positions [18].
• If the crack propagates a small amount compared to the crack length, there are no noticeable changes at the stresses in the crack tip during crack propagation [14].
The VCCT uses the nodal forces at the nodes at the crack tip (F xi , F yi ) and the relative displacements at the previous nodes (u i , ν i ) [19] as shown in Figure 2. The SERR for mode I fracture is computed as: ( ) ( ) where Δa is crack propagation length between adjacent nodes.Element size at the precrack zone is selected 0.25 mm, 0.5 mm, 1 mm and 2 mm respectively to perform the mesh sensitivity analysis.Element size through the thickness which is equal to the ply thickness is kept 0.25 mm for the mesh sensitivity analysis.In the mid-plane (0 o //0 o interface), all of the models include initial delamination with a length of 50 mm and extended crack path from the delamination tip to model end.3-D 8-node linear interface element INTER 205 is used to create a predefined crack path.G Ic value is determined from experimental P-δ results and used in FE analysis to initiate the fracture and to separate the coincident nodes along the predefined crack path.
Each arm of the model is displaced contrary along y-axis until 20 mm crack opening displacement is obtained.Motions of the arms in x and z directions are set free.The numerical model is shown in Figure 3.

Experimental Results
Initially, DCB tests are performed and typical P-δ curve is drawn (Figure 4).As can be seen in Figure 4, NL point and 5%/max point are so close each other and thus in each tests P cr values and corresponding to δ cr values are almost identical while VIS point exhibits lower values.Because of that the edges of the specimen are fiber-bridging free, P cr value detected at VIS point is lower than the values at NL point and 5%/max point.Fiber bridging is also caused to variable crack extension rates (stick-slip behavior) at all DCB specimens.It was an inherent outcome of the excessive warp and the stitching yarns in unidirectional E-glass fabric.Similar observation was reported by Shokrieh et al [20].3. G Ic value obtained by MBT method is used for the numerical analysis since MBT method calculates more conservative G Ic value than the other two [4].Similar results comparing the methods are seen at experimental results.The highest G Ic is found by using CC method while the lowest is obtained by MBT method.According to the NL point, G Ic obtained by CC method is 5.27 % higher than G Ic obtained by MBT method.
G Ic values are calculated using P cr and δ cr values detected at the VIS point and 5%/max point (5%/max G Ic ) and listed in Table 4. G Ic value according to the VIS point (VIS G Ic ) is lower than that of the NL point.Since the VIS point method does not consider fiber bridging, P -δ values at the NL point are assumed as critical ones.Shokrieh et al. [20] also explained that NL point is a criterion for delamination onset.
Test method ASTM D5528-01 [4] indicates G Ic calculated by NL point (NL G Ic ) is lower than VIS G Ic and 5%/max G Ic for composites with tough matrix.For the composites with brittle matrix, NL G Ic is very close to VIS G Ic [4].Whereas VIS G Ic is evaluated as 18% lower than NL G Ic in this experimental study.Brittle matrix and excessive fiber bridging of vacuum infused E-glass/epoxy laminates may be the reason for the above contradiction.

Numerical Results
The comparisons of the critical displacements and the critical forces obtained from the experiments and the numerical analysis with different mesh sizes are shown at Table 5.According to the results, the mesh size is selected as 2 mm since the numerical critical force with Δa = 2 mm are closer to the experimental result.The comparison of the numerical and the experimental P-δ curves is given in Figure 6.It is seen that VCCT determines P cr value of the unidirectional fiber reinforced model (L0_4mm) is 10.15% higher than the experimental ones.On the other hand, corresponding δ cr value of the numerical model is 30.43%lower than the experimental displacement.This discrepancy was explained with linear elastic model of ANSYS programme by Bonhomme et al. [14].They said that the DCB specimen had non-linear property [14].LEFM also doesn't simulate likely damages taken place in DCB arms at the time of loading [21].Additionally VCCT assumes that DCB model is perfectly built in and ignores correction factors (Δ, n, A 1 ) which are specified in experimentally calculations of strain energy release rate.Fiber bridging phenomenon is not taken care by numerical model can also explain this mismatch.the flexural rigidity increase with the specimen thickness.The difference between P cr of L0_a and L0_c is nearly 60%. Figure 8 represents P-δ curves of the unidirectional and the multidirectional finite element models with 4 mm thickness.
Although the precracks of all the models located in 0 o //0 o interface, the cross ply model and all of the angle ply models needs more crack opening displacement than L0_b model to initiate cracking except for L15_a.Highest P cr occurred at the model L30 and the lowest occurred at the model L15_a.Figure 9 exhibits the stacking sequence effect on P cr and δ cr .Among the three L15 models, L15_c model has highest P cr and δ cr values and L15_b model has highest rigidity.It is also seen as a result that B t ratio is not proportional to symmetry of the G Ic distribution for E-glass/epoxy.Previous researchers represent the effect of D c and B t ratios just on G Ic distribution along the model width.For this reason, the relation between extended crack length and D c and B t ratios is investigated finally.Crack extensions and the difference in crack extension length along the width can be seen in Table 6.Skewed crack extension causes differences in the extended crack length at y=0 mm edge and y=25 mm edge except for the models L15_a and L15_b.The extended crack length in models L30, L45, L60 changes along the specimen width.On the other hand, the models L15_a and L15_b have uniform extended crack length although they have skewed delamination front.As a result, asymmetric or symmetric G Ic distribution does not affect the extended crack length.So it is concluded that extended crack length can not considered as a function of D c or B t ratios in this study.

CONCLUSIONS
In this study, characterization of delamination crack in multidirectional E-glass/epoxy composite under mode I loading is researched numerically.Initially G Ic of [0 16 ] unidirectional E-glass/epoxy was determined by DCB test and then experimentally determined G Ic was used to initiate the delamination in the numerical models.The effect of stacking sequence and fiber orientation of sublaminates and the effect of thickness on mode I delamination crack length and G Ic distribution along the specimen width in 0 o //0 o interface have been investigated.3-D models were used to simulate mode I fracture by VCCT.The results showed that; • VIS G Ic is found 18% lower than NL G Ic .But the test method ASTM D 5528 indicates NL G Ic is very close to VIS G Ic for brittle composites.Brittle matrix and excessive fiber bridging of vacuum infused E-glass/epoxy laminates may be the reason for this conclusion.
• While P cr value of the unidirectional fiber reinforced model (L0_b) is 10.15% higher than the experimental ones, experimental δ cr is 30.43% higher than the numerical δ cr .Linear elastic model chosen at ANSYS® [14], possible damages in displaced arms of DCB [21] and neglected correction factors (Δ, n, A1 ) by VCCT can explain this mismatch.
• D c ratios higher than 0.03 exhibits asymmetric G Ic distribution along width.Maximum skewness of G Ic distribution is seen for L45 model with maximum D c ratio.
• The model L15_c with [90 2 ,+15 2 ,-15 2 ,0 2 ] s stacking sequence exhibits asymmetric G Ic distribution in spite of its B t ratio is zero.Magnitude of the B t ratio is not a scale for the symmetry of G Ic distribution for E-glass/epoxy laminates.
• Extended crack length can not be considered as a function of D c or B t ratios.

Fig. 1 a
Fig.1 a) Schematic illustration of DCB specimen, b) Experimental setup 2.2 Test Procedure A universal testing machine (Shimadzu AGS-X 100kN) with a class 0.5 load cell is used to perform mode I interlaminar fracture toughness tests with the DCB specimens.Five [0 o ] 16 E-glass/epoxy DCB specimens are tested.The crosshead speed is chosen as 3 mm/min which is in the interval ( 1 -5 mm/min) advised by ASTM D5528-01 [4].A high precision digital camera that has a lens with 25 mm focal length is used to observe and to record the motion of the crack tip and so the delamination propagation.The load (P) and crosshead displacement (δ) values are recorded by the universal testing machine and the crack length (a) is recorded by visual observation.Tests are continued until the overall crack length is reached to 95 mm.

Fig. 2 2D
Fig.2 2D Delamination tip model for VCCT 3.2 Finite Element Simulation For the present analysis, finite element software ANSYS® with VCCT is used to determine the G Ic distribution along the specimen width at the precrack tip.Experimentally determined mechanical properties of the unidirectional E-glass/epoxy composite are used to develop finite element (FE) models.Nine E-glass/epoxy composites are modeled with different stacking sequences and fiber orientations.Fiber orientations and stacking sequences are chosen considering balance and symmetry of laminate.And three unidirectional E-glass/epoxy composites are modeled with different thickness.Codes, stacking sequences, D c (D 2 12 / D 11 D 22 ) and B t (= |D 16 |/D 11 ) ratios of the models are listed in Table 2.As the element type, SOLID 185 with eight nodes having three degrees of freedom in each node is used for 3-D modeling of the composite

Fig. 4
Fig.4 Experimental P-δ curve of DCB specimens SERR (G I ) of the specimens is calculated using visual inspection of the crack extension and acquired P-δ values.Calculated SERR (G I ) according to MBT, CC, MCC methods are shown in Figure 5. Experimental G Ic values calculated by MBT, CC and MCC methods at the NL point are listed in Table3.G Ic value obtained by MBT method is used for the numerical analysis since MBT method calculates more conservative G Ic value than the other two[4].Similar results comparing the methods are seen at experimental results.The highest G Ic is found by using CC method while the lowest is obtained by MBT method.According to the NL point, G Ic obtained by CC method is 5.27 % higher than G Ic obtained by MBT method.

Fig. 5
Fig.5 Characteristic R-curve behavior obtained by visual inspection

Fig. 6
Fig.6 Comparison of experimental and numerical P-δ curves P-δ curves of the unidirectional numerical models with different thickness are shown in Figure 7.It is seen that P cr and

Fig. 9 Fig. 10 G
Fig.9 Numerical P-δ curves of DCB specimens with different stacking sequences of L15 model

Fig. 11 GFig. 12 a
Fig.11 G Ic distribution at precrack tip along the width of L15 models with different stacking sequences

Table 1 .
Physical and mechanical properties of the unidirectional E-glass/epoxy

Table 2 .
D c and B t ratios of composite models with different stacking sequences

Table 3 .
G Ic of DCB specimens according to NL point

Table 4 .
G Ic of DCB specimens according to VIS point and %5/max point

Table 5 .
Comparison of experimental and numerical DCB test results

Table 2 and
[7]ures 10, 11 and 12show the effect of B t ratio on G Ic distribution along width that if the B t ratio is nonzero, asymmetric G For example, E 1 /E 2 ratio of E-glass/epoxy used in this study is 2.95 whereas E 1 /E 2 ratio of AS4/8552 carbon/epoxy is 16.79[7].
[8]distribution is inevitable.But in spite of B t ratio is zero for the model L15_c, full symmetric G Ic distribution could not obtained along the width.On the other hand, although B t ratio of the model L60 (0.0270) is approximately four times lower than B t ratio of the model L30 (0.1144), similar asymmetric G Ic distribution is obtained.But Sebaey et al.[7]and Davidson et al.[8]explained that symmetry of G Ic distribution depend on the decrease of B t ratio for AS4/8552 carbon/ epoxy and C12K/R6376 graphite/epoxy respectively.These differences between the results may be explained by elastic properties of chosen materials.

Table 6 .
Extended crack length of numerical models at edges and middle section of the