Experimental-Numerical Investigation of

Experimental-Numerical Investigation of Delamination Buckling and Growth

R. Krüger, C. Hänsel, M. König ISD-Report 96/3, November 1996

ISD-Report 96/3

1

Abstract The growth of delaminations in CFRE (Carbon Fibre Reinforced Epoxy) specimens during R=-1 fatigue loading has been studied. Artificial circular delaminations have been introduced during manufacturing to simulate a pre-damaged structure and to cause delamination growth. Criteria based on fracture mechanics are used to describe the delamination failure. Predicting delamination growth with this approach requires the distribution of the local energy release rate along the delamination front. For obtaining this energy release rate distribution the virtual crack closure method was found to be most favourable for three-dimensional finite element analysis, as the separation of the total energy release rate into the contributing modes is inherent to the method and only one complete finite element analysis is necessary. Plots of measured delamination progression per load cycle versus computed energy release rates have been included in a Paris Law diagram obtained experimentally for the unidirectional interface. The results indicate that Paris Law parameters determined for a unidirectional interface are not valid for delamination propagation that occurs at interfaces between layers of dissimilar orientations. The corresponding Paris Law line seems to be less steep as for the the unidirectional case. This is essential when considering a damage tolerance approach for the design and operation of CFRE structures. However, the growth law for the interface considered should be confirmed by material characterization tests to support the current results.

ISD-Report 96/3

1

2

Objectives and Background

Motivated by the increasing use of composites in primary structural components, research has been focussed on the disbond of two adjacent fibre reinforced layers of a laminate which is a prevalent state of damage, commonly known as delamination. A recent survey on problems concerning composite parts of civil aircraft shows that delamination, mainly caused by impact, presents 60% of all damage observed [1]. In a 1991 IATA survey, air carriers reported that about 40% of all damage comes from platform or ground handling and maintenance [2]. Up to now failure caused by delamination is prevented by using empirically determined design criteria – based on maximum allowable strains (0,3%) – during layout and construction of components made of fibre reinforced materials [3, 4, 5]. For an optimal utilisation of the potential offered by those materials, as well as for the determination of inspection intervals however, it is essential to be able to predict delamination growth. In this investigation criteria based on fracture mechanics are used to describe the delamination failure. Propagation therefore is to be expected when a function of the mixed mode energy release rates GI , GII , GIII along the delamination front locally exceeds a certain value. This value can be regarded as a property of the interface and depends on the material and on the ply orientations of the layers adjacent to the plane of delamination. The goal of the investigation presented is to obtain information on the dependence of delamination growth on the mixed mode energy release rates for cyclic loading in real structures using a combined experimental-numerical procedure to verify the approach considered. For the experimental determination of critical energy release rates several simple test specimens have been developed, such as the double cantilever beam (DCB) specimen for determination of G I c , the end notched flexure (ENF) specimen and transverse crack tension (TCT) specimen for GII c . Following a global/local testing and analysis approach as shown in figure 1, failure criteria need to be verified and improved at each stage, i.e. from the level of material characterization via coupon and sub-element level up to the substructure and full structure level [6]. Therefore more complex specimens as that shown in figure 2 were used on a sub-element level. Looking at the fracture toughness data for untoughened carbon/epoxy materials we notice that the mode II values are about three to eight times higher than the corresponding mode I values. When investigating failure mechanisms, it is therefore necessary to consider the individual mode contributions of the energy release rate along a delamination contour. A simple analytical evaluation of the test data (as e.g. beam theory for DCB and ENF specimens) is no longer possible for the specimens considered. Computational methods based on three-dimensional finite element modelling are a meaningful and efficient tool, by which the energy release rate at delamination growth along the entire delamination front can be evaluated. The specimen shown in figure 2, made of unidirectionally reinforced prepreg material and containing an artificial circular delamination was designed to simulate a structural com-

ISD-Report 96/3

3

ponent subject to tension-compression loading which due to an impact has been damaged near the surface. Experiments with this specimen type will provide input data to set up a predictive model for onset and growth of a delamination located between layers of different ply orientations.

2

Experimental Programme

A rather extensive experimental program has been the first intermediate step in this study. The specimens - as shown in figure 2 - are cut from prefabricated plates made of prepreg = = = = material (Ciba Geigy T300/914C) with a stacking sequence of = == = = = . Artificial circular delaminations have been introduced at interface 2/3 during manufacturing by embedding two stacked foils, cut from a 20 m thick release film, to simulate a pre-damaged structure. The two stacked foils guarantee layer separation. An antibuckling guide with rectangular windows on both sides of the specimen is clamped to the test specimen, suppressing its global buckling. The loading , with a stress amplitude of 180 N/mm 2 and 220 N/mm2 for frequency is 10 Hz at R an initial delamination diameter of 20 mm and 10 mm, respectively. These stress amplitudes have shown to yield stable delamination growth under cyclic loading.

85 0

45

5 + 45

[ 5 +45 5 45 0

5]

= 1

2.1

Determination of Delamination Contours Using Moir e´ Technique

The principle of this optical in-situ technique for measuring out-of-plane displacements of a surface is illustrated in figure 3. When two repetitious patterns deform relative to each other, a Moiré effect is observed, resulting in a fringe pattern. The information necessary to obtain the surface deformation is contained in the Moiré picture. However, using this method, only that part of the displacement vector which is normal to the surface will be captured. By calibrating the Moiré fringes, it is possible to relate the fringes to out-ofplane displacements. During measurement, the object surface is illuminated with a grid pattern produced by a Michelson interferometer, which is monitored and digitized by a CCD video camera. Owing to the specific frame time (40 ms) of the CCD camera, it is necessary to stop dynamic loading and to monitor the test specimen at two static load levels. The difference of the phase patterns for the loaded and the unloaded specimen results in a phase distribution representing the buckling shape of the delaminated region. The evaluation is carried out by the phase-shift method. The phase-shift method utilizes the relationship between the brightness of the pixel and the induced phase shift to extract the underlying phase distribution of a fringe pattern from three recorded interferograms Æ of the optical path difference. By using a piezo separated by discrete phase shifts of controlled reflector within the Michelson interferometer, the phase shift is carried out for each point of the image. The resulting intensity distribution is the modulated interfer-

120

2

ISD-Report 96/3

4

ence phase. From the interference phase distribution, the out-of-plane displacements of the deformed shape are computed. These out-of-plane displacements result in an intensity distribution which can be figured in grey levels of constant displacements in the range of 0 - 255 (figure 4). Using additional image processing software, the size and shape of the delaminated sublaminate is determined from this information, yielding the delamination contours [7]. An example of delamination contours, obtained by this procedure from four different specimens, is presented in figure 5.

2.2

X-Ray Photographies and Ultrasonic Scans

A small number of specimens was subjected to X-ray radiography and ultrasonic scanning inspection, which was carried out at the German Aerospace Research Establishment (DLR), Institute for Structural Mechanics in Braunschweig. Matrix cracks parallel to the fibres are observed along the edges of the artificial delamination before any significant delamination propagation becomes apparent (see figures 6 through 8 for a specimen with a delamination of 10 mm diameter). These matrix cracks prevent the delamination from growing perpendicular to the loading direction as can be seen in the X-ray photographies after 200,000, 500,000 and 1 million load cycles (figures 9 through 12) as well as in corresponding images obtained from ultrasconic scanning after 500,000 and 1 million load cycles (figures 13 and 14). A similar observation has been described by O’Brien in a recent publication [8].

2.3

Determination of Smoothed Delamination Contours

The delamination contours for specimens with an initial delamination diameter of 10 mm obtained using Moiré technique, as described in the previous section, are shown in figure 15. As indicated above, delamination growth is only observed in load direction. Due to this observation it was assumed that the fronts can be approximated by elliptical contours the minor axes of which remain constant and equal to the radius of the original artificial delamination. Using the image processing software mentioned above, the delaminated area is determined from the information obtained by the Moiré inspection. Keeping the minor axes constant (5 mm), the major axes of the corresponding ellipse covering the equivalent delaminated area is easily calculated. Subtracting the initial radius (5 mm) from the calculated major axes yields the radial growth of the delamination for a particular specimen after a certain number of load cycles N . This information is plotted in figure 16 for five specimens. For comparison radial growth directly measured from X-ray photographies and images obtained from ultrasonic scanning has been included in the plot of figure 16. The values obtained using these more accurate methods still show significant scatter. We therefore assume that the scatter observed for the values determined from the Moiré measurements

ISD-Report 96/3

5

is acceptable. Using a fit through the Moiré data the radial growth for a given number of load cycles N can be determined and the corresponding smoothed elliptical fronts can be found. This procedure has been used for obtaining the fronts s 1, s2 , s3 , s4 ,s5 and s6 after 100,000, 200,000, 300,000, 400,000, 500,000 and 1 million load cycles, which are shown in figure 17. These six fronts and the initial artificial front s 0 were used as input to the corresponding finite element models. In addition, buckling displacements w of the delaminated sublaminate and the base laminate have been obtained from the measured data and are displayed in figure 18. The diagram shows the lateral displacements w of the center of the delaminated sublaminate and the center of the base laminate as a function of the compression load [9]. This was done for a later examination of the correspondence between the measured and computed postbuckling deformations as described in section 4.2. In the case of an artificial initial delamination of 20 mm diameter, the picture of the delamination growth is significantly different, as can be seen from the figures 19 and 20. After appearance of the matrix cracks (as described for the 10 mm initial delamination) this delamination grows into a more rectangular form (figure 19). During further cycling the matrix cracks cause a complete separation of the sublaminate over the width of the initial delamination and this open delamination subsequently grows in length direction of the specimen (see figure 20). Hence, for the analysis of delamination growth of the specimen with 20 mm artificial initial delamination, this splitting of the sublaminate has to be taken into account. The analysis of this situation will be part of a future study.

3

Analytical Tools

3.1

A Continuum-Based 3D-Shell Element

In order to achieve correct results for mode separation when computing energy release rates, the use of three-dimensional FE-models is required. Due to the extensive computation times already noticeable for 3D-models of simple DCB- and ENF-specimens [10, 11], the development of a layered volume element using a continuum-based three-dimensional shell theory [12] has been found to be necessary as

the computation of the complete load path using a layer of brick elements for each ply of the specimen will be extremely computer time consuming the standard isoparametric eight- as well as twenty-noded volume elements have the tendency to lock for small element thickness to element length ratios, leading to an unnaturally stiff behaviour of the structure during computation

ISD-Report 96/3

6

a volume-type element with eight nodes is necessary to assure complete compatibility with the contactor and target elements that are used to avoid structural overlapping in the vincinity of the crack front.

Several orthotropic layers of different orientations may be included in the developed layered 3D-shell element. An extended three-dimensional ABD-Matrix has to be supplied by the user as an input. The individual components of the matrix are calculated as for classical laminated plate theory [13] with

A =

m X i=1

Q [z i 0 (i)

z (i

( )

1)

]

m X 1 B = 2 Q [(z i ) (z i ) ] i m X 1 D = 3 Q [(z i ) (z i ) ] i where A denotes the extended membrane stiffness, D the flexural stiffness, B the couthe off-axis three-dimensional stiffness of the i-th ply and z i and pling stiffness, Q 0

(i)

( ) 2

(

1) 2

0

(i)

( ) 3

(

1) 3

=1

=1

0

z

(i)

( )

(i 1)

the distances of the surfaces of the i-th ply from the element mid plane as shown in figure 21. We thus obtain the three-dimensional relationship

"

Nij Mij

#

=

"

Aij Bij Bij Dij

# "

"ij ij

#

i; j = 1; 6

where the Nij and Mij denote the force and moment resultants, the A ij , Bij and Dij the stiffness coefficients and the "ij and ij the strains and curvatures respectively.

3.2

Virtual Crack Closure Method

The most significant step for the current approach is the accurate computation of the distribution of the mixed mode energy release rates along arbitrarily shaped delamination fronts, which means that a classical analytical evaluation of energy release rates - as possible for geometrically simple cases - is not applicable. Thus, computational methods based on three-dimensional finite element (FE) modeling are utilized. It has been found that the virtual crack closure method is most favourable for the computation of energy release rates because the separation of the total energy release rate into the contributions by the different crack opening modes is possible in a straight forward manner [14, 15]. When using this method only one FE-computation is necessary for a given delamination front which is beneficial especially when solving large geometrically non-linear problems.

ISD-Report 96/3

7

When using eight-noded elements - like the 3D-shell element described in the previous section - the procedure for the computation of the energy release rates can be illustrated by figures 22 and 23. The crack is virtually closed along the complete delamination front over the distance a 1. The energy release rate which is computed for the area A is assigned to node 0 on the crack front. This is mainly done to obtain the value at a nodal point of the FE-mesh which simplifies data management for postprocessing. For a specimen containing an arbitrarily shaped delamination contour, which is automatically meshed by commercially available software however, it is not useful to assume an equal element area A along the entire front. Therefore, in the most general case where the element widths bi as well as the element lengths ai may vary, the relative displacements u1; v1 and w1 computed at node behind the crack front for an element length a1 have to be corrected to fit to the forces X 1 ; Y1 and Z1 computed directly at the front (node 0 ) for an element length a 2. This may be done by taking into account the shape functions of the elements. For volume elements with eight nodes the displacements vary linearly along the edges and we obtain the desired values which correspond to the computed forces at point 0 by linear interpolation ( a 1 > a2) or linear extrapolation ( a1 < a2) of the computed displacements at point . Further, allowing a variation of the element widths bi requires an adjusted calculation of the adjacent element surfaces according to Ai 21 a2 bi.

1

1

1

0

0

1

0

1

=

We thus obtain

GI GII GIII

= 12 A = 12 A = 12 A

1

1

1

1 + A 1 + A 1 + A

a v a

Y1

1

0

2

X1

0

2

Z1

0

2

2

a u a a w a 1

1

1 2

1 2

1

:

Other techniques for virtually closing the crack along the crack front are mentioned in [16, 17].

3.3

A Contact Processor to Avoid Interpenetration

Local contact occuring in the delaminated area is a phenomenon that has to be considered in each model used for investigating the delamination behaviour. In the zone of contact the delaminated sublaminate is locally supported by the base laminate, influencing the deformation behaviour of the sublaminate. Using classical linear elastic finite element analysis this phenomenon can not be simulated, i.e. an interpenetration of the different layers may not be prevented. Therefore a non-linear analysis becomes necessary if the contact problem has to be taken into account.

ISD-Report 96/3

8

The contact processor implemented in the NOVA software accounts for frictionless contact of deformable bodies by applying the penalty method [18]. One of the bodies serves as a so called contactor, the other as the target as shown in figure 24. The contact problem is described using two independent nets, one applied to the surface of the three dimensional mesh representing the body of the contactor the other to the surface of the corresponding mesh of the target. This is done to keep the description of the surfaces involved fairly simple. It needs to be mentioned that these surface elements are not structural elements in the common sense as they do not introduce additional degrees of freedom. Therefore these surface elements are initialized and handled separately. Although the setup is arbitrary, it is suggested – for mere numerical reasons – to chose the softer body to serve as the contactor. The nodal points connected to the mesh of the contactor are called dependent nodes, those belonging to the mesh of the target are called guiding nodes. A dependent node can touch a target element which is defined by the guiding nodes along its edges. Using the penalty method the amount of contact is directely controlled by the penalty factor. This allows for a user control of the process from the soft side, which is advantageous for many applications. On the other hand this may also lead to numerical problems if the penalty factor has been chosen too large. In most cases it is helpful to start with a small penalty factor which then has to be increased once equilibrium is reached. It has to be considered, however, that the state of equilibrium which can be attained depends on the penalty factor chosen. It is possible that the iteration stops due to the fact that the penalty factor chosen is too small, thus avoiding further iterations to reach a smaller error. Looking at figure 24 the kinematics of the problem can easily be explained. Three dependent nodes have penetrated the target violating the contact conditions and thereby creating an active contact element during the finite element computation. An active element consists of the dependent node and the target element involved. The computed stiffness matrix provides the relationship between the contacting bodies and the resulting contact forces. It has to be considered that the contact elements should be kinematically compatible with those elements they are attached to, thus transmitting the forces of the base elements in a kinematically compatible manner. The contact algorithm has to push the dependent nodes onto the target surface following a vector normal to this surface. The intersection of this vector with the surface yields the physical location C of the contact. The length of the vector is equal to the distance. In a numerical analysis contact can only be achieved within certain limits, i.e. the dependent node will always slightly penetrate, even if equilibrium is reached. Therefore this distance is important when checking for the contact conditions. Looking at figure 24 we notice that a dependent node is free if the distance is greater than zero and is in contact or penetrates once this distanace becomes negative. Defining the maximum allowable penetration as gap the kinematics of the contact problem can be described as follows:

distance > 0: the dependent node is free

gap < distance < 0: the dependent node is in contact

ISD-Report 96/3

9

distance < gap: the dependent node has violated the contact condition and has penetrated. It therefore needs to be brought back to the target surface.

The maximum allowable penetration depth gap is a parameter, the value of which is to be set by the individual user. The contact processor implemented in the NOVA software package allows each of the dependent nodes to get into contact with each of the target elements, therefore the size and required space of the problem to be analyzed depends on the maximum difference in nodal point numbers of the surface meshes involved. In order to save space and computation time it is therefore adviseable not to extend the area meshed beyond the zone of expected contact and further to carefully chose the numbering of the nodal points.

4

Numerical Simulation of Delamination Buckling and Growth

To be in a position to gain better insight into the mechanism of delamination growth, a close correlation of the experimentally observed behaviour of the specimens and the results of its numerical model counterparts is essential. Hence, an accurate representation of the postbuckling behaviour of the delamination specimens is required. Reliable energy release rate distributions, which will be used to formulate a growth criterion, can then be computed. The numerical simulation in this investigation is entirely focussed on the analysis of the specimen containing an artificial delamination of 10 mm diameter. For the specimen containing an artificial delamination of 20 mm diameter only the initial circular front has been investigated. Results of the numerical simulation of delamination buckling are published in [19].

4.1

Material and Geometrical Data

The following material data were considered for the fibre matrix combination Ciba Geigy T300/914C (notation according to Tsai [13] with index 1 denoting the fibre direction):

E1 E2 = E3 E5 = E6 = G31 = G12 21 23

=

130000 N/mm2

=

9500 N/mm2

=

4300 N/mm2

= =

0:32 0:436

ISD-Report 96/3

10

The geometrical data were: Section modelled with finite elements: 55 mm 40 mm (equal to size of cutout of antibuckling device, see figure 25) = = = = = = = = = Layup: = Thickness: t mm mm embedded between the second and third Diameter of artificial delamination: d layer.

[ 5 + 45 =2

4.2

5

45 0

85 0 45 = 10

5 + 45

5]

Computation of Delamination Buckling Behaviour

In a preliminary study quadrilateral layered plate elements (with four nodes) which allow first order transverse shear deformation have been used to keep computing times low. The investigations during this phase of the project focussed on the boundary conditions to be applied along the modeled section of the specimen, the influence of the delamination size modeled and the global performance of the nonlinear FE-code used [20, 21]. Only the region inside the window of the antibuckling guide was modeled by two subplates with the nodes located at the plane of delamination over the entire region modeled. It was found that the boundary conditions applied to the edges of the region modeled have no significant influence on the computed deformations, provided we focus our interest on the region of stable delamination growth (e.g. load levels of 220-240 N/mm 2). This has been checked by applying first a clamped and then a simply supported condition. The same holds for the value of the linear buckling load for sublaminate buckling [9, 22]. Due to the existence of plies, models of the entire window section are necessary for accurately computing mixed mode energy release rates along the delamination fronts, even in the case of geometrically symmetric specimens. Using the newly developed layered 3D shell element with eight nodes [12], the top two layers of the sublaminate layer have been modeled by one element over the thickness of each layer. Also the of the base laminate which is adjacent to the plane of delamination has been modeled by one element over the thickness. The remaining 13 layers of the base laminate have been = = and = = = = layup joined in two elements of = respectively. A simply supported condition has been applied as boundary condition along the edges of the modeled section. This boundary condition represents the edges of the antibuckling guide. As already mentioned above the boundary conditions chosen do not have a significant influence on the buckling behaviour of the sublaminate. The geometry of the specimen and the section modeled with finite elements are shown in figure 25.

[ 5]

[ 5

45 0

85]

[0

45

[+45]

5 + 45

5]

A comparison of the computed buckling displacements w with those measured during a quasi-static test and during the first load cycle of the fatigue tests (see figure 18) shows that the computed displacements for a 10 mm initial delamination are considerably smaller than those which were determined experimentally (see figure 26). However, it was found

ISD-Report 96/3

11

by X-ray radiography that it is probable that the effective delamination width is larger than the width of the embedded release film (up to 11 mm). This results from relative movement of the two release foils, which can occur during the manufacturing process of the plates. Additional 3D analyses with an assumed larger delamination diameter were carried out by opening one ring of elements in the existing model, which results in a delamination diameter of 10.6 mm. A significant drop in the computed linear buckling load was observed. The obtained displacement values of the sublaminate lie well within the scatter band of the experiments for the load level of interest. However, in both cases a perfect buckling behaviour is computed, while the experiments clearly show an imperfection from the beginning of the load history, which partly results from the embedded 40 m thick double release foil used to create the artificial delamination. This fact was taken into account in further analyses where the foil was simulated by a simple coordinate shift of the nodal points in the region of the delaminated sublaminate. Results for a 10.0 mm and 10.6 mm delamination are included in figure 26. The model with an increased delamination size and coordinate offset, representing the embedded double foil, yields the best results compared with the experiment. The remaining discrepancy may be accounted for by a rather loose fitting of the specimens into the antibuckling device. The deformed geometry of a specimen with the 10 mm initial delamination is shown in : N/mm 2. This compressive load is figure 27 for an applied compression load of close to the bifurcation point for global buckling as computed by linear buckling analysis. The buckling of the sublaminate is shown in detail in figure 28 for an applied compression : N/mm2 which corresponds to the load level as applied for the fatigue load of experiments.

= 280 0

= 220 0

4.3

Computation of Energy Release Rates

The most significant step in the current investigation is the accurate computation of energy release rates along arbitrarily shaped delamination fronts. Therefore, several numerical techniques based on the finite element method were investigated using two- and three-dimensional models of DCB and ENF specimens. Two methods - the virtual crack extension (VCE) [23, 24] and the virtual crack closure method (also called modified crack closure (MCC) method) [14, 15, 16] - were found to be particularly well-suited for threedimensional finite element analysis. In both methods, only one complete finite element solution is necessary for the calculation of energy release rates along a given delamination front. Using the virtual crack closure method is advantageous for this study as it provides the possibility to easily split up the total energy release rate G T into the contributing modes GI ; GII and GIII . Investigations employing DCB and ENF specimens assured the reliability of the techniques used for computation of energy release rate distributions along straight and curved delamination fronts [10, 11]. In addition, convergence studies for ENF and SLB specimens [25] did not show the reported non-convergence [26, 27] of

ISD-Report 96/3

12

the virtual crack closure method associated with the oscillatory singularity. This may be due to a very small bimaterial mismatch of the considered interfaces [28]. The issue will be investigated further. It has been found that the assumption of a slightly larger initial delamination diameter which results in a considerably larger post buckling displacement as shown in figure 26 also has a significant effect on the computed energy release rates. A slight modification of the initial delamination diameter, (i.e. 10.6 mm instead of nominally 10.0 mm) results in a 53% increase in computed total energy release rate as shown in figure 29. Results of analyses without and with the release foil (thickness 40 m) which have been included in figure 29 show an increase in computed energy release rates up to 290% for the model with foil. Therefore the release foil needs to be included in the model [29, 19]. Computed mixed mode energy release rates for a model which accounts for the release foil and assumes a circular delamination of 10.6 mm diameter are shown in figure 30, for the maximum compression load (220.0 N/mm 2) applied in the fatigue experiments. The computed energy release rates in the phase of tension loading are negligible, which confirms previous assumptions that delamination growth occurs in the compression phase GII GIII and GT , we expect of the test. Looking at the distribution of G I , Gshear the delamination to start growing into load direction, whereas the potential to grow perpendicular to the load direction seems to be very low, as computed results in this area are significantly below the threshold values G I th and GII th (see section 5 on Paris Law). This is in agreement with experimental results shown in section 2.2. It is, however, different from the observations in [30] and [31] where an isotropic material behaviour and a homogeneous quasi-isotropic laminate were assumed, respectively. This indicates that delamination growth depends on the stacking sequence of the delaminated sublaminate and the base laminate.

=

+

Figures 31 through 53 show details of the FE-meshes, the deformed geometries, contact forces and the distributions of the computed energy release rates along the propagating delamination for fronts s 1 ; s2; s3 ; s4; s5 and s6 respectively. The distribution of the energy release rate along the front s 1 (i.e. after 100,000 load cycles) as plotted in figure 33, is significantly different from the distribution computed for the initial circular front. The global maximum of the total energy release rate G T is now found perpendicular to the loading direction (s ; and s ; ) with a second maximum occuring in loading direction. The values for interlaminar shear failure also reach their maximum perpendicular to the load direction (s ; und s ; ) and a second maximum for s and s : . This also holds for the GI distribution where again the maximum is reached perpendicular to the loading direction, in loading direction however it ceases to zero. Due to this distribution the delamination propagation would be expected to occur primarily perpendicular to the load direction, this however is not observed as discussed in section 2.2.

0 25

=05

0 75

0 22

0 72

=0

Looking at the detail of the deformed geometries (figures 32, 35, 39, 43,47, 51) a signifi-

ISD-Report 96/3

13

cant crack opening is observed perpendicular to the loading direction while the opening in load direction vanishes for higher numbers of load cycles, with the delaminated sublaminat touching the surface of the base laminat. Interpenetration of the layers is prevented by using a contact processor that utilizes a contactor target concept applying the penalty method as described in section 3.3. Contact occurs first after 200,000 load cycles as indicated by the computed contact forces shown in figure 36. With an increasing number of load cycles, i.e. increasing delamination size, the contact zone becomes larger as shown by the corresponding contact forces displayed in figures 40, 44, 48 and 52. In the region of contact GI along the front ceases to zero as shown in the plots of the energy release rate distribution (figures 33, 37, 41, 45, 49 and 53). Perpendicular to the load the values increase while in load direction the energy release rate drops during delamination propagation. This is shown also in figure 54 where the total energy release rate GT along all fronts investigated is presented. The crack opening observed perpendicular to the load direction is associated with a high mode I contribution as shown in figure 55 causing a severe bending of the [+5/-5] plies of the delaminated sublaminate. This bending perpendicular to the loading direction leads to the observed matrix cracking parallel to the fibres as discussed in section 2.2. The current investigation does not include the development (formation) of these cracks, they may however be included by either modeling the discrete crack itself or by properly adjusting the bending stiffness of the cracked plies. Two different failure mechanisms – delamination propagation and matrix cracking – are observed, occuring simultaneously while this specimen is subjected to tension-compression fatigue loading. The development of the matrix cracks has not been included in the current computation. When moving along an elliptically shaped delamination contour, high energy release rates are being computed perpendicular to the loading direction as expected. Delamination growth however is prevented by the formation of matrix cracks early in the load history. For the following investigation on delamination propagation it is assumed that these distrubances do not significantly influence the results obtained in loading direction and that the computed energy release rates in this section of the front may be used for consideration of the growth behaviour of delamination propagation between [-5/45] plies.

5

Comparison with Paris Law

Mode I and mode II failure under cyclic loading have been characterized for the considered material in [32]. For mode I characterization unidirectional DCB specimens were tested under R=0.1 loading. The results are shown in figure 56. Commonly known unidirectional ENF specimens and TCT (Transverse Crack Tension) specimens developed at DLR were used for mode II characterization under R=0.1 loading, yielding the information shown in figure 57. The MMB fixture has been employed to measure delamination growth under mixed mode fatigue R=0.1 loading conditions [33]. The results obtained : are plotted in figure 58. for GI =GII

= 0 25

ISD-Report 96/3

14

For the results of the tests, a straight line can be obtained by regression that fits very well to a Paris Law

da dN

= c Gn ;

approximating the delamination growth rate da=dN as a function of the maximum cyclic energy release rate G. The complete characterization of delamination growth under cyclic loading should include the threshold G T th below which the delamination does not develop and a critical value GT c above which the growth is unstable. These values as well as the Paris Law parameters c and n are given in Table 1. There is quite a difference between Paris Law lines for mode I and mode II, indicating that the delamaintion growth rate in general depends on the mixed mode ratio G I =GII . This is consistent with the differences observed for fracture toughness [33]. For delaminations propagating in a 0 Æ/0Æ interface, measured progression fits well into a Paris Law as obtained experimentally using simple specimens for material characterization. These results suggest that growth prediction based on Paris Law is an approach to be considered [34, 35, 36, 37]. Therefore the delamination growth a was measured between two consecutive fronts along the fronts investigated (s 1 – s6 ). Resulting delamination growth rates a= N (N number of load cycles) have been plotted versus the computed total energy release rates of the corresponding front and included in the Paris Law diagram for the unidirectional interface, as a growth law for a -5 Æ/45Æ interface is currently not available. Looking at the resulting plot in figure 59, a significant difference in the growth behaviour for the -5 Æ/45Æ interface and the 0Æ /0Æ interface becomes obvious. The results indicate that Paris Law parameters determined for a unidirectional interface are not valid if delamination growth occurs between plies of dissimilar orientations. This is an interesting result with respect to an application of the concept of damage tolerance for CFRE structures.

TABLE 1 – Fatigue test and Paris Law data for T300/914C Specimen

GI =GII GII =GT

GT th

GT c

c1

n1

0.406

4.344

References

[kJ/m2] [kJ/m2] ENF

0.0

1.0

TCT MMB DCB 1

0.25

1

0.0609

0.450

0.8

0.067

0.37

0.0

0.068

0.149

Paris Law parameters (da=dN

= c Gn )

[32, 33] 58.613

2:44 10

6

6.3695

[33]

10.610

[32, 33]

ISD-Report 96/3

6

15

Conclusions and Outlook

The growth of delaminations embedded near the surface of laminated CFRE specimens during tension-compression fatigue loading has been investigated. From the buckling deformations, which have been determined experimentally via Moiré technique, delamination contours have been determined. Obtained delamination fronts are in satisfactory agreement with those extracted from corresponding X-ray photographies and ultrasonic scans. Matrix cracks prallel to the fibres have been detected in the top layers, preventing delamination propagation perpendicular to the load direction. The postbuckling state of the sublaminate in the compression phase of the test, which is driving the delamination growth, has been computed by three-dimensional (layered 3D shell) finite element models. Computations show that slight variations in delamination size, applied loading and the consideration of the release foil in the model have a significant influence on the buckling behaviour of the sublaminate and the corresponding energy release rates, while the boundary condition of the window section that has been modeled only affects the buckling behaviour of the entire laminate. The assumption of a slightly larger initial delamination diameter (10.6 mm instead of 10.0 mm), simulating the as manufactured specimens, yields a computed load deflection behaviour that is in satisfactory agreement with the experiments. Local mixed mode energy release rates have been computed along a circular initial delamination front and along smoothed elliptical fronts after 100,000, 200,000, 300,000, 400,000, 500,000 and 1 million load cycles using three-dimensional models, which employ the virtual crack closure method. Along the elliptically shaped delamination contours high energy release rates are being computed perpendicular to the loading direction as expected. Delamination propagation however is observed in loading direction while the formation of matrix cracks early in the load history prevent growth in the direction perpendicular. The current model does not account for the formation of these matrix cracks but it is assumed that the resulting disturbances are negligible and do not effect the computed energy release rates in loading direction which are used to characterize the growth behaviour of a -5 Æ/45Æ interface. Although delamination propagation occurs at an interface between layers of dissimilar orientations, plots of measured delamination progression per load cycle versus computed energy release rates have been included in a Paris Law diagram obtained experimentally for the unidirectional interface. This was done, as a growth law for the interface investigated is currently not available. Looking at the results, a significant difference in the growth behaviour becomes obvious, indicating that Paris Law parameters determined for a unidirectional interface are not valid in this case. This is essential when considering the application of a damage tolerance approach to CFRE strctures. Following a global/local testing and analysis approach, failure criteria need to be verified and improved at each stage, i.e. from the level of material characterization via coupon and

ISD-Report 96/3

16

sub-element level up to the sub-structure and full structure level and vice versa. Therefore the current results need to be verified on the level of material characterization using DCB, ENF and MMB tests yielding the growth law for the interface of interest. This work points to the need for computations of mixed mode energy release rates along measured arbitrarily shaped delamination fronts and the derivation of a generalized crack driving force G? taking into account mixed mode conditions, the orientations of the plies adjacent to the plane of delamination, and the angle of the delamination front with respect to the ply orientations. Finally, the investigation of delamination growth with respect to the number of applied load cycles will result in a delamination progression versus G ? relation, which is expected to fit into Paris Law.

Acknowledgement This work has been supported by the Deutsche Forschungsgemeinschaft (Az.: Kr 668/171, Kr 668/17-2). The authors would like to thank Dr. Michael G”adke, Dr. W. Hilger, H.C. Goetting and B. Friederichs for the specimen manufacturing, the X-ray radiography and the ultrasonic scanning of the specimens, carried out at the German Aerospace Research Establishment (DLR), Institute for Structural Mechanics in Braunschweig.

ISD-Report 96/3

17

References [1] A.G. Miller, D.T. Lovell, and J.C. Seferis. The evolution of an aerospace material: Influence of design, manufacturing and in-service performance. Composite Structures, 27:193–206, 1994. [2] V.P. McConnel. Getting a fix on repair. May/June:19–24, 1994.

High-Performance Composites,

[3] A.C. Garg. Delamination - A damage mode in composite structures. Eng. Fracture Mech., 29(5):557–584, 1988. [4] H. Eggers, H.C. Goetting, and H. Bäumel. Synergism between layer cracking and delaminations in MD-laminates of CFRE. In AGARD Structures and Material Panel, Patras, Greece, 1992. [5] W.G.J. Hart and R.H.W.M. Frijns. Delamination growth in improved carbon composites under constant amplitude fatigue loading. Technical Report NLR TP 89008 U, National Aerospace Laboratory NLR, 1989. [6] R.H. Martin. Local fracture mechanics analysis of stringer pull-off and delamination in a post-buckled compression panel. In A. Poursartip and K. Street, editors, The Tenth International Conference on Composite Materials, Vol. I, pages 253–260. Woodhead Publishing Ltd., 1995. ISBN 1-85573-222-1. [7] C. Hänsel and K. Eberle. Measuring propagation of delaminations in CFRPlaminates by Moiré technique. In Composites Testing and Standardisation, ECCMCTS, Amsterdam, pages 417–424. EACM, European Association for Composite Materials, 1992. ISBN 2-9506577-1-0. [8] T.K. O’Brien. Local delamination in laminates with angle ply matrix cracks, part I: Delamination fracture analysis and fatigue characterization. In Composite Materials: Fatigue and Fracture, Fourth Volume, ASTM STP 1156, pages 507–538. American Society for Testing and Materials, 1993. [9] J. Albinger, M. König, and C. Hänsel. Berechnung des Ausbeulens einer CFKDruckprobe mit einseitiger oberflächennaher Delamination und Vergleich mit dem Experiment. ISD-Bericht 92/7, Institut für Statik und Dynamik der Luft- und Raumfahrtkonstruktionen, Universität Stuttgart, 1992. [10] S. Fan, R. Krüger, and M. König. Anwendung numerischer Methoden der Bruchmechanik auf die Delamination von CFK-Laminaten. ISD-Bericht 92/11, Institut für Statik und Dynamik der Luft- und Raumfahrtkonstruktionen, Universität Stuttgart, 1992.

ISD-Report 96/3

18

[11] R. Krüger, M. König, and T. Schneider. Computation of local energy release rates along straight and curved delamination fronts of unidirectionally laminated DCBand ENF-specimens, AIAA-93-1457-CP. In The 34th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, La Jolla, California, pages 1332–1342, 1993. [12] H. Parisch. A continuum-based shell theory for non-linear applications. Int. J. Num. Meth. Eng., 38:1855–1883, 1995. [13] S.W. Tsai. Composite Design. Think Composites, 4th edition, 1988. ISBN 09618090-2-7. [14] E.F. Rybicki and M.F. Kanninen. A finite element calculation of stress intensity factors by a modified crack closure integral. Eng. Fracture Mech., 9:931–938, 1977. [15] F.G. Buchholz, H. Grebner, K.H. Dreyer, and H. Krome. 2D- and 3D- applications of the improved and generalized modified crack closure integral method. In S.N. Atluri and G. Yagawa, editors, Computational Mechanics ’88. Springer Verlag, 1988. [16] I.S. Raju, K.N. Shivakumar, and J.H. Crews Jr. Three-dimensional elastic analysis of a composite double cantilever beam specimen. AIAA J., 26:1493–1498, 1988. [17] K.N. Shivakumar, P.W. Tan, and J.C. Newman Jr. A virtual crack-closure technique for calculating stress intensity factors for cracked three dimensional bodies. Int. J. Fracture, 36:R43–R50, 1988. [18] H. Parisch. A consistent tangent stiffness matrix for three-dimensional non-linear contact analysis. Int. J. Num. Meth. Eng., 28:1803–1812, 1989. [19] R. Krüger, S. Rinderknecht, C. Hänsel, and M. König. Computational Structural Analysis and Testing: An Approach to Understand Delamination Growth. In E. A. Armanios, editor, Interlaminar Fracture of Composites, pages 181–202. Key Eng. Mater., Vols. 120-121, Transtec Publications Ltd., 1996. ISSN 1013-9826. [20] H. Parisch. NOVA User Manual. Institute for Statics and Dynamics of Aerospace Structures, University of Stuttgart, 1991. [21] H. Parisch. An investigation of a finite rotation four node assumed strain shell element. Int. J. Num. Meth. Eng., 31:127–150, 1991. [22] M. König, J. Albinger, and C. Hänsel. Delamination Buckling: Numerical Simulation of Experiments. In A. Miravete, editor, ICCM-9, Madrid, Vol. VI, pages 535–542. Woodhead Publishing Ltd., 1993. ISBN 1-85573-140-1. [23] T.K. Hellen. On the method of virtual crack extension. Int. J. Num. Meth. Eng., 9:187–207, 1975.

ISD-Report 96/3

19

[24] D.M. Parks. A stiffness derivative finite element technique for determination of crack tip stress intensity factors. Int. J. Fracture, 10:487–502, 1974. [25] R. Krüger. Three Dimensional Finite Element Analysis of Multidirectional Composite DCB, SLB and ENF Specimens. ISD-Report No. 94/2, Institute for Statics and Dynamics of Aerospace Structures, University of Stuttgart, 1994. [26] I.S. Raju, J.H. Crews, and M.A. Aminpour. Convergence of Strain Energy Release Rate Components for Edge-Delaminated Composite Laminates. Eng. Fracture Mech., 30(3):383–396, 1988. [27] C. Hwu and J. Hu. Stress Intensity Factors and Energy Release Rates of Delaminations in Composite Laminates. Eng. Fracture Mech., 42(6):977–988, 1992. [28] H. Gao, M. Abbudi, and D.M. Barnett. Interfacial Crack-Tip Field in Anisotropic Elastic Solids. J. Mech. Phys. Solids, 40(2):393–416, 1992. [29] J. Albinger, C. Hänsel, M. König, R. Krüger, H. Parisch, and S. Rinderknecht. Ein kombiniertes experimentelles und numerisches Verfahren zur Bestimmung der Energiefreisetzungsraten beim Delaminationsfortschritt in CFK. In D. Dinkler and U. Hänle, editors, Numerische und experimentelle Methoden in der Statik und Dynamik, pages 265–289. Institut für Statik und Dynamik der Luft- und Raumfahrtkonstruktionen, Universität Stuttgart, 1994. ISBN 3-930683-00-8. [30] K.-F. Nilsson and B. Stor˚akers. On interface crack growth in composite plates. Journal of Applied Mechanics, 59:530–538, 1992. [31] J.D. Whitcomb. Three-dimensional analysis of a postbuckled embedded delamination. J. Composite Mat., 23:862–889, September 1989. [32] R. Prinz and M. Gädke. Characterization of interlaminar mode I and mode II fracture in CFRP laminates. In Proc. Int. Conf.: Spacecraft Structures and Mechanical Testing, pages 97–102. ESA SP-321, 1991. [33] M. König, R. Krüger, K. Kussmaul, M. v. Alberti, and M. Gädke. Characterizing Static and Fatigue Interlaminar Fracture Behaviour of a First Generation Graphite/Epoxy Composite. To appear in ASTM STP 1242, 1996. [34] M. König and R. Krüger. Computation of Energy Release Rates: A Tool for Predicting Delamination Growth in Carbon Fibre Reinforced Epoxy Laminates. In E. Oñate, D.R.J Owen, and E. Hinton, editors, Computational Plasticity. – Fundamentals and Applications – Proceedings of 4th International Conference,Barcelona, pages 1167–1178. Pineridge Press, 1995. ISBN 0-906674-85-9.

ISD-Report 96/3

20

[35] M. König, R. Krüger, and M. Gädke. Prediction of Delamination Growth Under Cyclic Loading Using Fracture Mechanics. In A.H. Cardon, H. Fukuda, and K. Reifsnider, editors, Progress in Durability Analysis of Composite Systems, Proceedings of the international conference DURACOSYS 95, Brussels, pages 45–52. A.A. Balkema Publishers, 1996. ISBN 90-5410-809-6. [36] R. Krüger and M. König. Prediction of Delamination Growth Under Cyclic Loading, To appear in. In E.A. Armanios, editor, Composite Materials: Fatigue and Fracture - 6th Vol., ASTM STP 1285. American Society for Testing and Materials, 1997. [37] R. Krüger, M. König, and M. Gädke. Predicting Delamination Growth under Cyclic Loading: An Approach Using Computational Structural Analysis and Testing. In A. Poursartip and K. Street, editors, The Tenth International Conference on Composite Materials, Vol. I, pages 561–568. Woodhead Publishing Ltd., 1995. ISBN 1-85573-222-1.

ISD-Report 96/3

21

Figures

Structure

Sub-Structure

Specimen D

Specimen B

Specimen A,C

Material Characterization

Coupon

Figure 1: Global/local testing and analysis approach

ISD-Report 96/3

22

87.0 150.0

1.0 2.0 87.0 60.0

+5 -5 +45 +5 -5 -45 0 +85 -85 0 -45 -5 +5 +45 -5 +5

Figure 2: Specimen investigated

ISD-Report 96/3

23

CCD-Camera

Michelson− Interferometer Figure 3: Principle of measurement by Moiré technique

10.0 mm

Figure 4: Example of grey level distribution of out-of-plane (buckling) displacements

ISD-Report 96/3

24

8

6

4

2

0

-2

-4

-6

-8 -8

-6

-4

-2

0

2

4

6

8

10.0 mm

Figure 5: Example of delamination contours, obtained from four specimens

ISD-Report 96/3

25

matrix crack

artificial delamination between plies 2 and 3 Figure 6: X-ray photography of specimen 1069/4 after 10,000 load cycles (DLR Braunschweig) – Surface application of contrasting agent

matrix crack

artificial delamination between plies 2 and 3 Figure 7: X-ray photography of specimen 1069/4 after 10,000 load cycles (DLR Braunschweig) – Internal application of contrasting agent

ISD-Report 96/3

26

Figure 8: D-Scan of specimen 1069/4 after 10,000 load cycles (DLR Braunschweig)

ISD-Report 96/3

27

end tab

artificial delamination between plies 2 and 3

growing delamination between plies 2 and 3

end tab

Figure 9: X-ray photography of specimen 1071/1 after delamination growth (DLR Braunschweig)

ISD-Report 96/3

28

artificial delamination between plies 2 and 3 matrix crack

growing delamination between plies 2 and 3 after 200,000 loadcycles Figure 10: X-ray photography of specimen 1071/1 after 200,000 load cycles (DLR Braunschweig)


growing delamination between plies 2 and 3 after 500,000 loadcycles Figure 11: X-ray photography of specimen 1071/1 after 500,000 load cycles (DLR Braunschweig)

ISD-Report 96/3

29


growing delamination between plies 2 and 3 after 1,000,000 loadcycles

Figure 12: X-ray photography of specimen 1069/3 after 1,000,000 load cycles (DLR Braunschweig)

ISD-Report 96/3

30

Figure 13: C-Scan of specimen 1071/2 after 500,000 load cycles (DLR Braunschweig)

ISD-Report 96/3

31

Figure 14: D-Scan of specimen 1070/2 after 1,000,000 load cycles(DLR Braunschweig)

ISD-Report 96/3

32

8

8

6

6

4

4

2

2

0

0

-2

-2

-4

-4

-6

-6

-8

-8 -8

-6

-4

-2

0

2

4

N=0

6

8

8

8

6

6

4

4

2

2

0

0

-2

-2

-4

-4

-6

-6

-8

-8

-6

-4

-2

0

-8

-6

-4

-2

0

-8

-6

-4

-2

-8

-6

-4

-2

2

4

6

N=50,000

8

-8 -8

-6

-4

-2

0

2

4

6

8

N=100,000

8

8

6

6

4

4

2

2

0

0

-2

-2

-4

-4

-6

-6

-8

2

4

6

8

0

2

4

6

8

0

2

4

6

8

N=200,000

-8 -8

-6

-4

-2

0

2

4

6

8

N=300,000 8

8

6

6

4

4

2

2

0

0

-2

-2

-4

-4

-6

-6

-8

N=400,000

-8 -8

-6

-4

-2

0

2

4

6

8

N=500,000

N=1,000,000

Figure 15: Delamination contours determined by Moiré technique

ISD-Report 96/3

Figure 16: Experimentally determined delamination growth

33

ISD-Report 96/3

34

s

Figure 17: Delamination contours (smoothed)

Figure 18: Experimentally determined displacements for specimens with circular delaminations of 10 mm and 20 mm diameter

ISD-Report 96/3

35

Figure 19: C-scan of specimen containing 20 mm artificial delamination (DLR Braunschweig) –Initial delamination growth

Figure 20: C-scan of specimen containing 20 mm artificial delamination (DLR Braunschweig) – Advanced delamination growth

ISD-Report 96/3

36

ζ 1

3

4

2

7

5

8 zi layer i

6

ξ

η

Figure 21: 3D-shell element ∆ a1

a

length ∆ a 2 crack closed

y, v

Y1’ ∆ v1 X1’ ∆ u1

X1’

x, u

Y1’

Figure 22: Virtual crack closure method, cut through plane of delamination

ISD-Report 96/3

37

crack front cracked area

uncracked area

∆A area virtually closed

1

1’ x,u,X

b1

∆ A1 ∆ A2 area

b2

z,w,Z 1 1’

virtually closed x,u,

X

z,w

,Z

forces displacements

∆a1 ∆ a2

Figure 23: Virtual crack closure method, top view of the plane of delamination

ISD-Report 96/3

38

contactor contactor element target element

contact

C

gap target penetration guiding node dependent node Figure 24: Contact analysis

87.0 150.0 87.0

60.0

Figure 25: Specimen and section modelled

ISD-Report 96/3

Figure 26: Comparison of experimentally determined displacements with FE results

39

ISD-Report 96/3

40

Figure 27: Deformed geometry for initial delamination front s 0

5.0 mm

Figure 28: Detail of deformed geometry for initial delamination front s 0

ISD-Report 96/3

41

Figure 29: Distribution of energy release rates along initial delamination front s 0

Figure 30: Distribution of energy release rates along initial delamination front 10.6 mm diameter, including foil

s0

with

ISD-Report 96/3

42

Figure 31: Detail of FE-mesh for delamination front s 1

Figure 32: Detail of deformed geometry for delamination front s 1

ISD-Report 96/3

Figure 33: Distribution of energy release rates along delamination front s 1

43

ISD-Report 96/3

44



ISD-Report 96/3

45

Figure 36: Contact force along delamination front s 2


ISD-Report 96/3

46



ISD-Report 96/3

47



ISD-Report 96/3

48



ISD-Report 96/3

49



ISD-Report 96/3

50



ISD-Report 96/3

51



ISD-Report 96/3

52

Figure 50: Detail of FE-mesh for fdelamination ront s 6


ISD-Report 96/3

53



ISD-Report 96/3

54

Figure 54: Distribution of total energy release rates G T along delamination fronts s 1 – s6

Figure 55: Distribution of mixed mode ratio along delamination fronts s 1 – s6

ISD-Report 96/3

Figure 56: Paris Law for mode I, obtained using DCB specimen (R=0.1)

Figure 57: Paris Law for mode II, obtained using ENF and TCT specimens (R=0.1)

55

ISD-Report 96/3

Figure 58: Paris Law for G I =GI I = 0:25, obtained using MMB specimen (R=0.1)

Figure 59: Results obtained in comparison to Paris Law for the 0 Æ /0Æ interface

56