Computation of Onset and Growth of Delamination in Double Cantilever beam Specimens Subjected to Fatigue Loading

In this article, the delamination onset and growth behavior of double cantilever beam (DCB) specimens has been presented. The modeling of a debonded region using master and slave surface technique for DCB specimens is done in ABAQUS CAE. The analysis of DCB specimens comprising of fatigue cyclic load has been done in ABAQUS. An onset and Paris delamination growth regimes are plotted. The growth regime being linear in log-log scale, the prediction of constants of this regime has been obtained using the polyfit command in the MATLAB environment. To obtain these constants has been explained in this article. Comparison of experimental and analytical results is shown for delamination growth. The strain energy release rate values for threshold and critical are indicated on the graphs. The number of cycles for delamination onset and growth has been tabulated for various load cases.

${\delta }_{\mathrm{max}}$                 Maximum displacement
${\delta }_{\mathrm{min}}$                 Minimum displacement
${A}_{o},\text{\hspace{0.17em}}b$                Cyclic loading constants
$\omega$                Circular frequency
t               Incremental time
${t}_{o}$                 Starting time
a                 Delamination length
${a}_{o}$                 Initial delamination length
$\Delta a$                 Element length
${G}_{\mathrm{Im}ax}$                 Maximum strain energy release rate in mode-I
${G}_{Ic}$                 Critical equivalent strain energy release rate in mode-I
${G}_{th}$                 Threshold strain energy release rate
${N}_{o}$                 Number of cycles for delamination onset
${N}_{G}$                Number of cycles for delamination growth
${N}_{T}$                  Total number of cycles
$da}{dN}$                 Delamination growth per loading cycle
${m}_{0},\text{\hspace{0.17em}}{m}_{1}$                 Power law constants
${c}_{1},\text{\hspace{0.17em}}\text{\hspace{0.17em}}{c}_{2}$                 Material constants
C, m                 Paris law constants

Delamination is the most common modes of failure in layered composites structures. Delamination is created when two layers de-bond from each other1. The analysis of delamination is commonly divided into the study of an initiation and growth under load/displacement control. To characterize an onset and growth of delamination, the use of fracture mechanics has become common practice over the past three decades2. The virtual crack closure technique (VCCT) method is based on the hypothesis that when a delamination extends by a small amount, the energy released in the process is equal to the work required to close the delamination to its original length3. The VCCT is widely used for computing energy release rates based on results from continuum finite element analysis to predict the mode of separation especially in the mixed-mode fracture criterion4,5. This method has been applied to calculate the individual modes of strain energy release rate (SERR) from the stresses and displacements6. As new approaches for analyzing composite delamination are incorporated into FE codes, the need for comparison and benchmarking becomes important7-10.

Understanding delamination growth behavior under cyclic/fatigue loading is necessary for damage tolerant design and reliability assessment of composite structures. As these structures often operate under goes fatigue loading. The behavior of composite material is more complicated than the metals. In a composite material, damage starts very early and the extent of the damage zone grows steadily while the damage type in these zones may change. To study the fatigue phenomena in composite materials Fracture mechanics is commonly used11,12. Fracture mechanics relate the variation of SERR with the crack growth. Fatigue loading is related to sinusoidal stress cycles of constant/variable amplitude loading. Delamination growth rate is defined as the delamination extension per number of cycles. The correlation of the delamination growth rate with the amplitude of energy release rate is most commonly represented in a log-log diagram. The curve can be divided into three regions according to the curve shape. In region I there is a threshold value, Gth, below which crack remain dormant or additional crack growth is negligible. Above this value, the crack growth increases relatively quickly. Region II defines a stable crack growth generally characterized by a linear part of the curve in the log-log plot. Finally, in Region III the crack growth rate rises asymptotically that corresponds to the critical fracture toughness values, Gc, where static fracture is achieved. To describe the delamination growth rate Paris law is widely used13,14. Paris law describes the linear portion (Region II) of the curve.

In this study, the computation of an onset and growth of delamination in a double cantilever beam (DCB) specimen has been presented. The DCB Specimens were modeled using two-dimensional (2D) plane strain elements. Experimental results were taken from Saponara15, et al. for a DCB specimen. Plotted the complete delamination growth curve, from this obtained the linear part to this fitting of the polynomial. From the fit obtained the Paris constants C and m for the delamination growth region II.

2.1 Specimen Dimensions

The double cantilever beam (DCB) with isometric view is shown in Fig. 1, has the following dimensions: B = 25.0 mm, 2h = 3.0 mm, 2L = 150.0 mm, a0 = 30.5 mm. The DCB lay up consist of 24 layers of unidirectional (00), layers, each layer thickness = 0.125.

2.2   Material Properties

The DCB specimens consist of T300/1076 unidirectional graphite/epoxy with a unidirectional layup. The material properties are given in Table 1(material properties) and fracture toughness values and the parameters required for delamination onset and growth are given in Table 2 (fracture parameter).

Table 1. Material properties considered for the model.

Table 2. Fracture parameters.

2.3  Finite Elements Mesh Details

A two-dimensional finite element model of a DCB specimen is shown in Fig. 2(a). The specimen was modeled with solid plane strain elements (CPE4) and element length is Δa = 0.5 mm, detail numbers are indicated in Table 3. The models were created as upper and lower part of the specimens with identical nodal point coordinates in the plane of delamination. Two surfaces (top and bottom) were defined to identify the contact area in the plane of delamination as shown in the zoom view in Fig. 2(b) of Fig. 2(a). Additionally, a node set was created to define the intact (bonded nodes) region.

Table 3. Mesh details for DCB model.

3.1  Static Analysis

The numerical computation has been performed using ABAQUS finite element software package16,17. In this article to avoid unnecessary complications, experimental anomalies such as fiber bridging18,19 were not addressed. Static analyses were performed to ensure that the model geometry and material input data are correct and that the model responds as expected. This step can be used to check whether the peak loading leads to static delamination propagation or not. If the delamination does not propagate, increase the sample static load. Based on the results it was assumed that the geometry, maximum applied displacements and material input were defined correctly and model adequately represented the benchmark case.

First, models simulating specimens with seven different delamination lengths (30 mm ≤ a0 ≤ 40 mm) were analyzed. For each delamination length modeled, the reaction loads P at the location of the applied displacement were calculated and plotted versus the applied opening displacement δ/2 as shown in Fig. 3.

Figure 3.Load-displacement behavior for different delamination lengths.

The critical load Pcrit and δcrit/2 critical displacement were calculated for each delamination length modeled and the results were included in the load/displacement plots. The results indicate that, with increasing delamination length, lesser load magnitude is required to extend the delamination. This means that the DCB specimen exhibits unstable delamination propagation under load control. Therefore, prescribed opening displacement δ was applied in the analysis instead of nodal point loads P to avoid problems with numerical stability of the analysis.

3.2  Fatigue Delamination Onset

The number of cycles to delamination onset, NO, may be obtained by solving Eqn. (1). An onset curve is the power law fit of the experimental data obtained from a DCB test using the respective standard for delamination growth onset.

$G={m}_{0}\cdot {N}_{O}^{{m}_{1}}$         (1)

${N}_{O}={\left(\frac{1}{{m}_{0}}\right)}^{\frac{1}{{m}_{1}}}\cdot {G}^{\frac{1}{{m}_{1}}}⇒{N}_{O}={c}_{1}\cdot {G}^{{c}_{2}}$          (2)

where c1 and c2 are material constants evaluated by experimentally.

For the cyclic loading of the specimen, guidance was taken from a draft standard designed to determine mode I fatigue delamination propagation20. In the draft document, it is recommended to start the test at a maximum displacement,δmax , which causes the energy release rate at the front, GImax, to reach initially about 80 % of GIc.

$\frac{{G}_{\mathrm{Im}ax}}{{G}_{Ic}}=0.8$          (3)

In the ASTM document21, it is suggested to use a load ratio R = 0.1 for testing. The corresponding minimum load, Pmin, and minimum displacement,δmin were calculated.

$R=\frac{{P}_{\mathrm{min}}}{{P}_{\mathrm{max}}}=\frac{{\delta }_{\mathrm{min}}}{{\delta }_{\mathrm{max}}}=0.1⇒{\delta }_{\mathrm{min}}=0.1\cdot {\delta }_{\mathrm{max}}$        (4)

Further, it was suggested in the ASTM document to use a frequency f =10 Hz for testing. The applied displacement δ taken from22 is represented as a function of time, t.

$\delta =\left[{a}_{0}+{b}_{1}\cdot \mathrm{sin}\omega \left(t-{t}_{0}\right)\right]\cdot {\delta }_{\mathrm{max}}$         (5)

where δmax = 0.67 mm is the maximum displacement. The constants a0 = 0.55, b1 =0.45, the circular frequency $\omega$ = 20 π = 62.832 and the starting time t0 = 0.025 are calculated from load ratio R = 0.1 and the frequency f = 10 Hz for testing. The sensitivity of the load ratio for mode-I DCB is shown in reference23. The resulting Eqn. (5) to calculate the applied displacement δ is being applied.

In a second step of analysis, a single load cycle was analyzed, starting at the previously applied maximum displacement. This step was performed to check that the amplitude input was defined correctly and resulted in the desired periodic cyclic loading during the analysis. Based on a comparison with the desired fatigue loading it is assumed that the amplitude input was defined correctly and the increments are small enough to adequately represent the desired periodic fatigue cyclic defined by Eqn. (5).

3.3  Fatigue Delamination Growth

The number of cycles during stable delamination growth, NG, can be obtained from the fatigue delamination propagation relationship (Paris law) plotted in Fig. 4. The delamination growth rate can be expressed as a power law function.

$\frac{da}{dN}=C\cdot {G}_{\mathrm{max}}^{m}$        (6)

where da/dN is the increase in delamination length per cycle and Gmax is the maximum energy release rate at the front at peak loading. The constant C and exponent m are obtained by fitting the curve to the experimental/computational data obtained from DCB tests. The critical SERR or fracture toughness, GIc, was included in the plot of Fig. 4. Once a delamination length is reached where the SERR less than or equal to (≤) cutoff value, Gth, below which delamination growth is dormant.
The number of cycles during stable delamination growth can be obtained by solving Eqn. (2) for NG

${N}_{G}=\int dN=\int \frac{1}{C}{G}_{\mathrm{max}}^{-m}.da$            (7)

For practical applications, Eqn. (2) can be replaced by an incremental equivalent expression.

$\frac{\Delta a}{\Delta N}=C\cdot {G}_{\mathrm{max}}^{m}$        (8)

For the current study, increments of Δa=0.5 mm were chosen. Starting at the initial delamination length a0=30.5 mm, the SERR Gimax were obtained for each increment, i, from the curve fit plotted in Fig. 5. These energy release rate values were then used to obtain the increase in delamination length per cycle or growth rate Δa/ΔN from the Paris law. The number of cycles during stable delamination growth, NG, was calculated by summing the increments ΔNi.

${N}_{G}=\sum _{i=1}^{k}\Delta {N}_{i}=\sum _{i=1}^{k}\frac{1}{C}{G}_{i,\mathrm{max}}^{-m}\cdot \Delta a⇒a={a}_{0}+\sum _{i=1}^{k}\Delta a={a}_{0}+k\cdot \Delta a$        (9)

where k is the number of increments. The resulting delamination length, a, was calculated by adding the incremental lengths Δa to the initial length a0.

Figure 4. Delamination growth rate (Paris law) for T300-914C.

3.4   Polynomial Fit for Constants Determination

The MATLAB24 has an easy command for finding the slope and intercept for a straight line on a log domain. The polynomial fit (polyfit or cftool) command render two numbers, the first of these numbers, is the slope and the second number is y-intercept. As the plot is on log-log scale, either loge or log or ln all these three are same. Simply log means natural logarithm. To convert a natural logarithm to base10 logarithm divide by the conversion factor of 2.303. To the obtained the constants numbers raising them to 10 power equally on both sides of an equation and using laws of exponent. Or using the tools, fit the power curve directly either in an Excel25 (Microsoft) or the MATLAB software. There is no big rigmarole about the polynomial fit procedure. The crack growth data along with a linear plot is shown in Fig. 5. Fig. 5 is the linear portion of the Fig. 4. The linear Eqn. (y = mx + C) should not be compared with power Eqn. (y = C.xm) to obtain the values of Paris exponent m and coefficient C as in Sahoo26, et al.

Figure 5. Poly fit for power equation on a log-log plot.

3.5   Combined Fatigue Delamination Onset and Growth

For the combined case of delamination onset and growth, the total life, NT, may be expressed as
NT = NO + NG      (10)
where NO, and NG are the number of cycles to delamination onset and growth respectively. For the initial NO cycles, the delamination length remains constant, followed by a growth up to NG cycles, the delamination length increases following the Paris law. For different displacement amplitudes, the corresponding NO and NG cycles are indicated in Table 4. As the applied amplitude of displacement decreases there is increase in the number of cycles for delamination initiation/onset and decrease in growth. Figure 6 shows a typical plot for a displacement magnitude of 0.67, other plots are not shown here for brevity.

Table 4. Delamination onset and growth cycles for varying displacement applied.

Figure 6. Delamination onset and growth behavior for T300/ 914C material.

In this study, modeled the DCB specimen using ABAQUS CAE. The delamination, region has been modeled using master and slave surface bonding technique. DCB specimen subjected to fatigue loading of constant amplitude.

Plotted the complete delamination growth curve from this obtained the linear part to this fitting of the polynomial is described. From the fit obtained the Paris constants C (6.59E6) and m (10.78) for the delamination growth Region-II.

It is found that the delamination initiation occurred at 935 cycles for maximum displacement amplitude of 0.67 mm. As the applied amplitude of displacement decreases there is increase in the number of cycles for delamination initiation. Below 0.5 amplitude of displacement there is no delamination initiation and growth. As there is decrease in magnitude of load/displacement (0.75, 0.67, and 0.5) there is increase in the number of cycles for delamination onset (35, 935, and 1.2x106) and decrease in number of cycles for growth (5.9x106, 5.5x106, and 4.3x106) respectively.

The authors thank CSIR for granting the Supra Institutional Projects SIP-STTD-08, NAL. This work has been part of this project work. The authors express gratitude to the Head, STTD, Dr Satish Chandra. We also thank Padmashree Prof. B. Dattaguru, Mr. M. Subba Rao, and Dr Gangan Pratap for their moral support in carrying out this work. Dr Ranganath V.R., STTD, NAL, has been acknowledged immensely for the number of discussion and suggestions rendered. Authors also thank all those who are indirectly supported in making this work possible.

1. O’Brien, T.K. Fracture mechanics of composite delamination. ASM Handbook, vol. 21, Composites: ASM International, 2001, pp. 241–245.

2. O’Brien, T.K. Characterization of delamination onset and growth in a composite laminate. NASA Technical Report No. NASA-TM-81940, 1981.

3. Anderson, T. Fracture mechanics: Fundamentals and applications. Ed. 2. CRC Press, 1994.

4. Krueger, R. Virtual crack closure technique: History, approach and applications. NASA/CR-2002-211628, ICASE Report No. 2002-10, 2002.

5. Benzeggagh, M.L. & Kenane, M. Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus. Compos. Sci. Technol., 1996, 56(4), 439-49.

6. Krueger, R. & Dirk, Goetze. Influence of finite element software on energy release rates computed using the virtual crack closure technique. NASA/CR-2006-214523, NIA Report No. 2006-06, 2006.

7. Krueger, R. An approach to assess delamination propagation simulation capabilities in commercial finite element codes. NASA/TM-2008-215123, 2008.

8. Adrian, C. Orifici & Ronald, Krueger. Assessment of static delamination propagation capabilities in commercial finite element codes using benchmark analysis. NASA/CR-2010-216709, NIA Report No. 2010-03, 2010.

9. Adrian, C. Orifici & Ronald, Krueger. Benchmark assessment of automated delamination propagation capabilities in finite element codes for static loading Finite Elem. Anal. Des., 2012, 54, 28–36. doi:10.1016/j.finel.2012.01.006.

10. ABAQUS example problems manual. ABAQUS® Standard, Version 6.10, DSS Simulia, 2010.

11. Sham Prasad, M.S.; Venkatesha, C.S. & Jayaraju, T. Experimental methods of determining fracture toughness of fiber reinforced polymer composites under various loading conditions. J. Min. Mater. Charact. Eng., 2011, 10(13) 1263-75.

12. Pradeep, K.R. & Nageswara Rao, B. Interface fracture assessment on sandwich DCB specimens. J. Reinf. Plast. Comp., 2010, 29(13) 1963-77.

13. Pugno, N.; Ciavarella, M.; Cornetti, P. & Carpinteri, A. A generalized Paris’ law for fatigue crack growth. J. Mech. Phys. Solids, 2006, 54(7), 1333–49.

14. Nguyen, O.; Repetto, E.A.; Ortiz, M. & Radovitzky,  R.A. A cohesive model of fatigue crack growth. Int. J. Fracture, 2001, 110(4) 351-69.

15. Saponara, Valeria La; Muliana, H.; Hai-Ali, Rami & Kardomateas, G.A. Experimental and numerical analysis of delamination growth in double cantilever laminated beams. Eng. Fract. Mech., 2002, 69, 687-99.

16. Abaqus Analysis User’s Manual, ABAQUS® Standard, Version 6.10, DSS Simulia, 2010.

17. ABAQUS theory manual. ABAQUS® Standard, Version 6.10, DSS Simulia, 2010.

18. Lin, Ye & Friedrich, K. Fibre bridging in double cantilever beam specimens and its effect on mode I interlaminar fracture toughness. J. Mater. Sci. Lett., 1992, 11(22) 1537-39.

19. Johnson, W.S. & Mangalgiri, P.D. Investigation of fiber bridging in double cantilever beam specimens. J. Compos. Tech. Res., 1987, 9(1), 10-13.

20. Konig, M.; Krueger, R.; Kussmaul, K.; Von Alberti, M. & Gadke, M. Characterizing static and fatigue interlaminar fracture behavior of a first generation graphite/epoxy composite. In the 13th Composite materials: Testing and design, Vol. 13, ASTM STP 1242. edited by S.J. Hooper, 1997, 60-81.

21. ASTM. Standard test method for mode I fatigue delamination growth onset of unidirectional fiber-reinforced polymer matrix composites. ASTM D6115-97(2011). www.astm.org. [Accessed on 01 March 2013].

22. Krueger, R. Development of a benchmark example for delamination fatigue growth prediction, NASA/CR-216723, NIA Report No. 2010-04, 2010.

23. Turon, A.; Costa, J.; Camanho & Davila, C.G. Simulation of delamination in composites under high-cycle fatigue. Composites Pt. A- Appl. S., 2007, 389(11), 2270–82.

24. http://www.mathworks.in/help/matlab/index.html, [Accessed on 28 February2013].

25. http://chandoo.org/wp/category/excel/ [Accessed on 02 March2013].

26. Sahoo, P.K.; Dattaguru, B.; Manjunatha, C.M. & Murthy, C.R.L. Fatigue de-bond growth in adhesively bonded single lap joints. Sadhana, 2012, 37(1),79-88.

 Dr Krishna Lok Singh Dr Krishna Lok Singhdid his BE from University Visvesvaraya College of Engineering and ME from M.S. Ramaiah Institute of Technology, of Bangalore University, in 1993 and 1996, respectively. He obtained his PhD (Aerospace Engineering) from the Indian Institute of Science in 2003. Currently, he is a senior scientist in the Structural Technologies Division of the National Aerospace Laboratories, Bangalore. His primary fields of interest are : Computational methods, fatigue of structures and materials and health monitoring and management of aerospace vehicles. Mr Madhu K.S. Mr Madhu K.S.did his BE (Mechanical Engineering) under Vivesvaraya Technological University, Belgaum, in 2011. Presently he is working as Mechanical Engineer in Aero Structures Group of Safran Engineering Services India Pvt. Ltd. Bangalore. Wg. Cdr. Mallikarjun Vaggar (Retd) Wg. Cdr. Mallikarjun Vaggar (Retd)obtained MTech (Mechanical Engineering) from Institute of Armament Technology, Pune University, 1990. Currently he is pursuing his PhD, in Mechanical Engineering, Vivesvaraya Technological University (VTU) Belgaum. He also holds the position of Professor at Aeronautical Engineering Department of Dayananda Sagar College of Engineering, Bangalore.