Citation
Finite elements for free edge delaminations and delaminated anisotropic beams

Material Information

Title:
Finite elements for free edge delaminations and delaminated anisotropic beams
Creator:
Pinheiro, Marco Antonio Santos, 1947-
Publication Date:
Language:
English
Physical Description:
x, 131 leaves : ill. ; 29 cm.

Subjects

Subjects / Keywords:
Beams ( jstor )
Cantilever beams ( jstor )
Composite materials ( jstor )
Delamination ( jstor )
Laminates ( jstor )
Matrices ( jstor )
Residual stress ( jstor )
Resultants ( jstor )
Stiffness ( jstor )
Thermal stress ( jstor )
Aerospace Engineering, Mechanics and Engineering Science thesis Ph. D
Composite construction ( lcsh )
Composite materials -- Delamination ( lcsh )
Dissertations, Academic -- Aerospace Engineering, Mechanics and Engineering Science -- UF
Finite element method ( lcsh )
Genre:
bibliography ( marcgt )
non-fiction ( marcgt )

Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 1991.
Bibliography:
Includes bibliographical references (leaves 128-129)
General Note:
Typescript.
General Note:
Vita.
Statement of Responsibility:
by Marco Antonio Santos Pinheiro.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
The University of Florida George A. Smathers Libraries respect the intellectual property rights of others and do not claim any copyright interest in this item. This item may be protected by copyright but is made available here under a claim of fair use (17 U.S.C. §107) for non-profit research and educational purposes. Users of this work have responsibility for determining copyright status prior to reusing, publishing or reproducing this item for purposes other than what is allowed by fair use or other copyright exemptions. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder. The Smathers Libraries would like to learn more about this item and invite individuals or organizations to contact the RDS coordinator (ufdissertations@uflib.ufl.edu) with any additional information they can provide.
Resource Identifier:
027030612 ( ALEPH )
25605189 ( OCLC )

Downloads

This item has the following downloads:


Full Text













FINITE ELEMENTS FOR FREE EDGE DELAMINATIONS
AND DELAMINATED ANISOTROPIC BEAMS












BY.
MARCO ANTONIO SANTOS PINHEIRO




















A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA

1991



























In memory of my father-in-law

Fernando Pinto

still and forever a beacon to my family














ACKNOWLEDGMENTS



I would like to express my sincere gratitude to my advisor Dr. Bhavani V. Sankar for his continuous support, encouragement and constructive criticism during all my doctoral studies. His competence and understanding of his students always impressed me and are examples to be followed in my future life. Gratitude is extended to Professors U.H. Kurzweg, L.E. Malvern, M. Rao and C.T. Sun, members of my doctoral supervisory committee, for their assistance during my research. My appreciation goes to Dr. C.E. Taylor for his participation in the initial part of my work.

Special thanks are owed to the Brazilian Army for giving me the opportunity to come to the United States of America to pursue doctoral studies.

I would like to thank my friends at the Department of Aerospace Engineering, Mechanics and Engineering Science, whose friendship have been motivation for the completion of this work.

Finally, I am grateful to my wife, Fernanda, my daughter Simone, my son Luis Marco, my parents Jose and Nely, and my mother-in-law Geny. They have been always a source of encouragement, support and love.


iii















TABLE OF CONTENTS

Page

ACKNOWLEDGMENTS ........................................ iii

LIST OF TABLES ........................................... vi

LIST OF FIGURES .............. ............................ vii

ABSTRACT ............................................... ix

CHAPTERS

1 INTRODUCTION ........................................... 1

1.1 Motivation and Objectives ...................... 1
1.2 Literature Survey ................................. 4
1.3 Scope and Outline of this work .................. 8

2 AN OFFSET BEAM FINITE ELEMENT FOR MODELING FREE EDGE
DELAMINATIONS IN COMPOSITES ......................... 10

2.1 Introduction .................................... ..10
2.2 Laminated Beam .................................. 16
2.3 Fabrication Process .............................. 16
2.4 Free Edge Problem ...................... .... 22
2.5 Finite Element Solution .......................... 26
2.6 Finite Element Models .............................. 28
2.7 Strain Energy Release Rate ..................... 30
2.8 Stresses Ahead of the Crack ...................... 39
2.9 Numerical Results ............................... 43
2.9.1 Strain energy release rate ................. 44
2.9.2 Stresses ahead of the crack ............... 45
2.9.3 Conclusions ........... ............... 49

3 MODELING DELAMINATION IN COMPOSITE BEAMS AND APPLICATION TO TORSION PROBLEM ................. . 56

3.1 Introduction .................................... 56
3.2 Beam Equations ................................... ... 58
3.3 New Set of Force and Moment Resultants ........... 64
3.4 Equilibrium Equations ............................ 65
3.5 Finite Element Solutions ......................... 69
3.6 Strain Energy Release Rate ................... 74
3.7 Rigid and Gap Elements ........................... 75

iv










4 NUMERICAL RESULTS AND DISCUSSION ...................... 79

4.1 Numerical Results ................................ 79
4.2 Material Properties and Stiffness Coefficients ... 80 4.3 Cantilever Beam ................................... 81
4.3.1 Model with two sublaminates ................ 81
4.3.1.1 Beam under Force 1 (F1) ............. 81
4.3.1.2 Beam under Force 2 (F2) ............ 86
4.3.1.3 Beam under Force 3 (F3) ............ 86
4.3.1.4 Beam under Force 4 (F4) ............ 86
4.3.1.5 Beam under Force 5 (F5) ............ 87
4.3.1.6 Beam under Force 6 (F6) ............ 87
4.3.2 Model with top sublaminate ................. 89
4.3.2.1 Beam under Force 1 (Fl) ............ 89
4.3.2.2 Beam under Force 2 (F2) ............ 93
4.3.2.3 Beam under Force 3 (F3) ............ 95
4.3.2.4 Beam under Force 4 (F4) ............ 95
4.3.2.5 Beam under Force 5 (F5) ............ 97
4.3.2.6 Beam under Force 6 (F6) ............ 98
4.3.3 Comments on the finite element solution .... 98
4.4 Delaminated Cantilever Beam under Torsion ........ 102
4.4.1 Delamination at the middle-plane of the beam 105
4.4.2 Delamination position varying along the
thickness .................................. 108
4.5 Simply Supported Beam under Transverse Force ...... 112

5 SUMMARY AND FUTURE RESEARCH .......................... 116

5.1 Summary .......................................... 116
5.2 Future Research .................................. 118

APPENDIX ............................................... 119

REFERENCES .............................................. 128

BIOGRAPHICAL SKETCH .................................... 130
















V















LIST OF TABLES


Table Page

2.1 Comparison of energy release rate obtained
from different methods ............................. 46

2.2 Relative difference (%) between the results
from present method and other methods ............. 47

4.1 Stiffness coefficients ........................... 82

4.2 Displacements for model with 42 elements .......... 85

4.3 Displacements for model with 21 elements .......... 92

4.4 Strain energy release rate (G) for cantilever beam. 110

4.5 Strain energy release rate (G) for simply
supported beam ................................. 114


























vi















LIST OF FIGURES


Figure Page 1.2 Types of delamination .............................. 3
(a) Free edge
(b) Transverse delamination

2.1 Free edge delamination ............................ 11

2.2 Sublaminates ........................................ 13

2.3 Finite element model for strain energy
release rate ...................................... 15

2.4 Finite element model for peel stresses ............ 15

2.5 Values of h for top sublaminate .................... 19

2.6 Force and moment resultants (top sublaminate) .... 23 2.7 Forces and degrees of freedom ...................... 29
(a) Nodal and average forces
(b) Degrees of freedom

2.8 Zero-volume path for J-integral .................... 32

2.9 Crack-closure method .............................. 37

2.10 Model for spring constant calculation .............. 40

2.11 Stresses ahead of the crack tip for specimen A ... 50 2.12 Stresses ahead of the crack tip for specimen B ... 51 2.13 Stresses ahead of the crack tip for specimen C ... 52 2.14 Stresses ahead of the crack tip for specimen D ... 53 2.15 Influence of temperature variation ...................54

2.16 Influence of moisture level variation ............ 55




vii










3.1 Laminated composite beam ............................ 59

3.2 Force and moment resultants ......................... 63

3.3 Rigid elements: ..................................... 77
(a) Beam under torsion
(b) Interference on left side (front view)
(c) Interference on right side (front view)

4.1 Loading conditions at middle plane .................. 83
(a) Perspective view
(b) Front view

4.2 Model with 42 elements ............................. 84

4.3 Loading conditions at the bottom plane .............. 90
(a) Perspective view
(b) Front view

4.4 Model with 21 top elements .......................... 91

4.5 Beam under loading F4 ............................... 96

4.6 Angle of twist along the length (L/h=100) .......... 100

4.7 Angle of twist along the length (L/h= 10) ......... 101 4.8 Delaminated beam .......................... ........ 103
(a) Under torsion
(b) Separated by parts

4.9 Finite element model ................................ 104

4.10 Delamination close to the top plane ............... 109

4.11 Strain energy release rate for cantilever beam .... 111 4.12 Simply supported beam ............................. 113

4.13 Strain energy release rate for simply supported
beam .............................................. 115

A Stiffness matrix K ................................ 120









viii















Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy

FINITE ELEMENTS FOR FREE EDGE DELAMINATIONS
AND DELAMINATED ANISOTROPIC BEAMS By

Marco Antonio Santos Pinheiro December, 1991

Chairman: Dr. Bhavani V. Sankar Major Department: Aerospace Engineering, Mechanics, and Engineering Science

The major objective of this study is to develop a

technique to analyze delaminations in composite beams by using the finite element method.

A method is proposed to obtain the strain energy release rate and the stresses ahead of the crack for the specific case of free edge delaminations. A shear deformable beam finite element with nodes offset to either top or bottom side has been developed. The thermal effects originating from the curing process are included in the formulation of the problem.

Three different expressions to calculate the strain

energy release rate are derived, one in terms of the force and moment resultants in the elements surrounding the crack-tip, and the two others in terms of the forces transmitted by a rigid element placed at the tip of the crack. The stresses ahead of the crack are obtained by using discrete spring ix









elements in the uncracked portion of the beam. The results calculated are in good agreement with results obtained from other studies of the same problem. The advantage of using the present method is that it requires less computational effort, and is capable of modeling any type of geometry, loading, and boundary conditions.

An efficient beam finite element has been developed

based on a generalized beam theory for anisotropic beams, in which the equilibrium equations are assumed to be satisfied in an average sense over the width of the beam. This introduces a new set of force and moment resultants for the beam. The displacements in the beam are expressed in terms of three deflections, three rotations and one warping term. The effectiveness of the method is illustrated by solving problems where the beam theory solutions are known. The strain energy release rate is determined for the cases of delaminated cantilever beam under torsion and simply supported delaminated beam under transverse loading. Studies are performed to understand the influence of the position of the delamination on the strain energy release rate.













x















CHAPTER 1

INTRODUCTION



1.1 Motivation and Objectives

Structural materials can bear defects proceeding from the fabrication processes or accumulate damage during their service life. Frequently, such defects and damage are not visually observable and have to be detected by using nondestructive evaluation techniques (NDE).

In a variety of applications, such as aerospace and aircraft structures, offshore platforms, bridges and tall buildings, the replacement of the imperfect parts is not a cost-effective solution for the problem. The costs involved in keeping the equipment out of order, the difficulty of access, and other reasons suggest that there is a need for research and development of methods to predict the occurrence and extent of the damages.

Composite materials and structures, for example, are highly susceptible to damage from low-energy impacts. Incidents such as tool-drops, improper handling or hailstorms (aircraft structures) can cause internal damage and significant reduction in compressive strength of the material. Under subsequent transverse loadings, the damage will grow


1









2

causing a redistribution of the stresses in the plies of a laminate and consequently reduction in the material stiffness and fatigue life. However, depending on the subsequent loading conditions and the damage characteristics, the structure often can still be used despite the existence of damage.

Failure in a composite material often initiates in the form of a matrix crack or a delamination. The former refers to intralaminar failure or transverse crack, while the latter refers to interlaminar failure. Both types of failures are illustrated in Fig. 1.1.

Delamination, or separation of the layers, is a mode of failure unique to laminated materials, and it is the predominant failure mode in continuous fiber composites. Based on the location and direction of the growth, we can discern two types of delamination: the free edge delamination and the transverse delamination. In the edge delamination the growth initiates at the load-free edges of the laminate and propagates normal to the load direction. Under transverse delamination, also known as crack tip or local delamination, the growth starts from a matrix crack and propagates normal to the transverse crack from which it originates. The Fig. 1.2 shows both types of delamination. Since delamination is a very commonly observed failure mode in composite materials, its study is of great importance in the field of failure prediction for composite structures.














matrix crack delamination









(a) (b)







Figure 1.1 Failures in composite materials

(a) Matrix crack (b) Delamination








matrix crack

delamination





delamination

(a) (b)


Figure 1.2 Types of delamination

(a) Free edge (b) Transverse delamination









4

1.2 Literature Survey

In the past few years, several researchers have studied delamination occurrence and growth. Usually they predict delamination growth from the strain energy release rate (G), and try to develop analytical and numerical methods capable of modeling this phenomenon. A brief description of some of these studies is important for accurate comprehension of the present work.

Chan and Ochoa (1987) developed a two-dimensional finite element model to calculate interlaminar stresses and strain energy release rate in composites subjected to uniaxial tension, bending and torsion. The element is an isoparametric element with eight nodes and three degrees of freedom per node. The basic idea in their work is to represent the effect of longitudinal or twisting curvature as a "load effect." The study concentrates on establishing relationships of G,, GH and GIII as functions of crack length and stacking sequence.

Pagano (1974) studied the case of fibrous composite laminates under a uniform axial extension and developed an approximate method to define the distribution of the interlaminar normal stress along the central plane of a symmetric laminate. A simple expression for the interlaminar normal stress was derived and the results obtained by using this expression were compared with those obtained by using a finite difference approximation of the elasticity equations. The two methods agreed well, and both results show that









5

significant stresses exist only in narrow regions next to the free edges. The length of this region is approximately equal to the laminate thickness.

Whitney (1986) proposed an approximation to the interlaminar normal stress distribution ahead of the crack and the influence of specimen geometry on Mode I strain energy release rate for edge delamination specimens. The specimens were subjected to tension and compression loadings and were of the classes [03/903]s or [8/-82/8/902]s. Whitney used a higher order plate theory, which considers transverse shear deformation and thickness-stretching. The effect of thermal residual stresses due to lamination fabrication was considered in the formulation of the problem. The numerical results suggested significant influence of the residual thermal stresses on both the stress distribution ahead of the crack and the Mode I strain energy release rate. In general, the thermal residual stresses increase the stresses at crack tip and the Mode I strain energy release rate for specimens under tension, and has the inverse effect for specimens under compression. Again the results confirmed that the effect of normal stress ahead of the crack dissipates within a distance of approximately twice the thickness from the crack front.

Armanios and Badir (1989) used a shear deformation theory to calculate the interlaminar stresses and energy release rate associated with both Mode I and mixed-mode edge delaminations. The formulation used by Armanio and Badir does not consider









6

the thickness strain effect. The analysis includes the effects caused by the residual thermal stresses and moisture stresses. The classes of laminates considered are the same ones used by Whitney (1986). Results suggest that the residual thermal stresses tend to increase the interlaminar stresses and total energy release rate for both modes, while the moisture content tends to alleviate the thermal effects. The value of moisture content necessary to alleviate totally the interlaminar stresses and the strain energy release rate ranged between 0.70% and 0.80% in the examples considered, and it was found that this value has a weak dependency on both the stacking sequence and the mode considered.

Usman and Gandhi (1989) studied the effect of moisture and temperature on the vibrational characteristics of composite materials. The nonlinear vibration of moderatelythick laminated composite plates under hygro-thermal conditions are considered by incorporating the effects of transverse shear and rotary inertia. The hygro-thermal changes along the thickness of the composite were explicitly incorporated in the form of the deformation field. The response of the composite under hygro-thermal environments was modeled by creating one model formed by two parts, one anisotropic solid (composite) and one ideal fluid (moisture). The equations presented in this work are significantly different from the corresponding equations in classical laminated composite plate theory.









7

Hu et al.(1989) developed a novel beam element to study the crack propagation in beam types of specimens. Finite elements are used to model the specimens, and the beam element has its nodes offset to one side, either bottom or top side. This approach allows modeling crack propagation without renumbering the nodes as the crack propagates, and is convenient in modeling the partial contact and slip at a delamination interface. The bottom and top un-delaminated portions of the beam were connected either by foundation springs or rigid elements. For some specimens the strain energy release rate was computed by using a zero-volume path integral and a crack closure method, and the results were found to be in good agreement. In the present study the beam theory is extended to analyze delaminated composite beams.

Sankar (1991 b) developed a beam theory for laminated composite beams derived from the shear deformable plate theory. The equilibrium equations are assumed to be satisfied in an average sense over the width of the beam, which introduces a new set of force and moment resultants for the beam. The displacements in the beam are expressed in terms of three deflections, three rotations, and one warping term. The closed solutions for the cases of a cantilever beam subjected to end loads and a specially orthotropic laminate beam under torsional loading are presented. The solution considers the shear deformation and the edge effect. The results in terms of angle of twist agree well with available solutions.









8

1.3 Scope and Outline of this Work

Numerical and analytical methods have been used to express the complex stress state near free edges and to calculate the energy release rate during crack propagation. Among several numerical methods, the finite element method was selected to be used in the present work. Although requiring intensive computational effort, this method is very popular and provides accurate solutions. Also, this method is very convenient in modeling structural elements of variable cross section.

In the present work, the typical element used to model the structure is the beam element. This element not only is very common in mechanical structures but also is the most appropriate to model some standard test specimens.

The scope of this work can be described as to develop a method to analyze delaminations in anisotropic composite beams, using the finite element method to model the problem and selecting the strain energy release rate as a delamination parameter. The loading is considered independent of time.

This dissertation consists of five chapters. Chapter 2 is focused on free edge delamination. A method is proposed for modeling this type of delamination using offset beam finite elements. The effect of temperature variation during the curing process (residual stresses) is considered in the problem. Chapters 3 and 4 are devoted to presenting a new and









9

efficient method to analyze composite laminated beams, using a more complete and accurate beam finite element. The capability of this new finite element in modeling the problem was verified by solving the problem of a cantilever beam subjected to different end loading conditions, including torsion. Applications to delaminated beams under torsion loading and the simulation of a quasi-static impact test on a delaminated beam are also presented in these chapters. Finally, chapter 5 summarizes the principal conclusions of the present study, and presents suggestions for future studies in this field.














CHAPTER 2

AN OFFSET BEAM FINITE ELEMENT FOR MODELING
FREE EDGE DELAMINATIONS IN COMPOSITES


2.1 Introduction

The separation of plies or delamination is a predominant failure mode in laminated composite structures. It occurs when interlaminar stresses exceed the interlaminar strength of the material, or as a result of the fabrication process, environmental degradation and improper handling of the material. Common examples of such failure regions during the normal use of composite structures are those near cut-outs, holes and free edges. In these regions, the abrupt changes in the geometry or material properties tend to intensify the stresses when the laminate is subjected to some external loads.

Specifically, the occurrence of delamination in free edges has been receiving increasing attention of many investigators, in their efforts to understand and to prevent delaminations in composite structural designs. The reason is the presence of high interlaminar stresses, the "peeling" stress, in the neighborhood of a free boundary. The "free edge problem" usually is modeled by a plate subjected to a uniform axial strain ( eo ), as shown in Fig. 2.1.


10














e 11









x/


















e 0
free edge delamination (applied strain)



Figure 2.1 Free edge delamination









12

During the usual fabrication process, the composite material is subjected to temperatures of several hundred degrees Fahrenheit, and then it is cooled down to the room temperature. The composite is assumed to be stress-free at the curing temperature. The change in temperature causes unavoidable stresses to the composite material layers when cooled to room temperature. Usually these stresses are called residual stresses or thermal stresses, and they are considered in the formulation of the present problem.

A simple procedure is developed in this study to obtain an approximate solution for the peel stresses ahead of the crack tip and to calculate the strain energy release rate, for the case of free edge delaminations including thermal effects. The problem depicted in Figure 2.1 is modeled using a shear deformable beam finite element with offset nodes (Sankar, 1989). The three-dimensional problem is modeled as a plane problem in the x-z plane assuming that the stresses and strains do not vary along the y-direction. This will be true for most of the length of the plate in the y-direction. The delamination plane divides the beam into two parts, the top sublaminate and the bottom sublaminate (Figure 2.2), with nodes offset to the bottom and top, respectively. The delamination plane is assumed to be the weakest one, therefore the crack propagates along this plane, without crack branchings.








13









free edge delamination








e

top sublaminate





x



bottom sublaminate





Figure 2.2 Sublaminates









14

When calculating the energy release rate, the uncracked region of the beam is connected by common nodes, except the tip of the crack, where a rigid element is used to connect the top and bottom nodes (Fig. 2.3). This rigid element not only ensures the continuity of displacements and rotations at that section, but also brings in its formulation the generalized forces that are transmitted between the sublaminates. These forces are used in the calculation of the energy release rate. In evaluating the peel stresses ahead of the crack, the uncracked region of the beam is connected by spring elements (Fig. 2.4), stiffness of which has to be judiciously selected.

The strain energy release rate has been successfully adopted by researchers in characterizing delamination. In this study its value is obtained from three different methods (Sankar and Pinheiro, 1990). The first one uses in its formulation the forces transmitted by the crack-tip rigid element. In the second method, the result is obtained by using a method similar to the crack-closure method used in twodimensional fracture problems. Finally, composites are very brittle materials,so that, the plastic zone in the vicinity of the crack tip is negligible, and hence the linear elastic fracture mechanics can be applied in which case the J-integral is equivalent to the energy release rate. Therefore, the third method consists in obtaining the J-integral from the beam model. The force and moment resultants in the region surrounding the crack-tip are used in this calculation.








15

cracked region uncracked region


top sublaminate rigid element common nodes bottom sublaminate

Figure 2.3 Finite element model for strain energy release rate




cracked region uncracked region top offset element




spring element




bottom offset element


Figure 2.4 Finite element model for peel stresses









16

The results from the three different methods are compared and found to be in good agreement.



2.2 Laminated Beam

In this section a homogeneous orthotropic laminate subjected to a uniform initial strain E = e0 applied over a length 2L is analyzed. The thickness of the laminate is 2h and the width is 2b. The thickness 2h is much smaller than 2L and 2L much smaller than 2b (Fig. 2.1). Due to symmetry conditions, just one half of the laminate is considered here.

The problem is divided into two different phases. In the first phase no initial strain is applied to the specimen, which is subjected only to thermal stresses due to the curing process. In the second phase the cured specimen containing the edge delamination is subjected to the axial strain e = eo in addition to the residual stresses.

The equations derived in the next sections refer to a laminated beam situated above the delamination plane. The expressions for the laminated beam situated below the delamination plane differ only in the integration limits ( -h to 0, instead of 0 to h ).


2.3 Fabrication Process

During the composite fabrication process there are heating and cooling phases. In the cooling phase, the temperature drops from the stress-free temperature to the room












17
temperature. In this case, each layer influences the expansion or contraction of the other because of the difference in their coefficient of thermal expansion. A similar situation exists when the composite absorbs moisture during service. As a result of the temperature variation, each layer has residual stresses after the laminate fabrication.

The effects of the residual stresses on the laminate as a whole are calculated and expressed in terms of thermal forces and moments. In the next phase, these thermal forces and moments are added to the external forces and moments.

During the cooling process the in-plane and transverse displacements are assumed to be of the form:

u(x,z)- uo(x) + z*X(x) (2.1) v(y,z)- vo(y) + zJ,(y)
w(x, z) - 0 (x)

where uo and vo are the axial displacements of points on the x-y plane in the x and y directions respectively, and W4 and '4 are the rotations about the y and x axes respectively. It is assumed that the modification in the thickness of the laminate is negligible when compared with the modifications in the other two dimensions. Therefore, w is assumed to be constant through the thickness of the laminate (no thickness strain).












18

From Equations 2.1 the strains for the cooling process are derived as
exx- uO,x + Z x.x- cox + zxx ey,- Vo + Zly- e0 + zxy (2.2) Yzx" aj + Wx A state of plane stress normal to the z direction is assumed during the curing process. Therefore, the thermal stresses developed within the kth layer (Agarwal and Broutman, 1980) are



Sk Q16026Q66 Y (2.3)



where Qi, are the stiffness coefficients and can be found in Whitney (1987), and


(x- ax AT
Cy a, AT
C -" aXY A T (2.4) In the above equations a is the thermal expansion coefficient and AT is the change in temperature.

By definition, the force (N) and moment (M) acting on the top laminate are

h
(Nx, Mx) -fx,[1, z] dz
0
h
(NY, MY) -foY[, z] dz (2.5)
0












19

In the present case, there are no net resultants acting on the laminate. Therefore, considering the definition of force and moment resultants, we have for the whole laminate "0 " (a YkdZ k-1 4 kdz (2.6) and



( k faXzdz (2.7) where n is the total number of layers on the top sublaminate and (hk - hk.1) is equal to the thickness of the kth layer (Fig. 2.5).


layer number





h 0 3





Sn-n
n --------h


plane of delamination



Figure 2.5 Values of h for top sublaminate











20

The stiffness coefficients of the top sublaminate are defined as

h
[a1y,biy, dI - fij [l, z, z2] dz (2.8)
0

i,j = 1,2,...,6

In the above equation, a shear correction factor equal to 5/6 (Whitney, 1987) is introduced in the calculation of the terms aii,bi, and d1i, wherever the index i is greater then 3.

Substituting Eq. 2.3 in Eq. 2.6 and 2.7, and using the definition of the stiffness coefficients in Eq. 2.8, we obtain the following expressions for the force and moment resultants

("-) - (0) - (a: ::) (:":) bi (: (:1) - (a; (2.9)








aZ;) - O d+ a22 ( 22 (XO) (N) (2.10)



where NTx and MTx are the thermal force and moment in the x direction and so are NTy and MT in the y direction.











21
The thermal forces and moments are defined as

h
[N',MT - (S 11x + D12 y + D6 ) [1, z] dz (2.11)
0

h
[NT, MyT] (12 Cx + Q22 y + Q26 CX) [1,z] dz (2.12)
0

Defining the matrices


[L]- bn d) (2.13)



1[M2 b12 (2.14) b12 d12)

and
[N]- 22 b22 (2.15) Equations 2.9 and 2.10 take the form of


)- [L + I (O) - (2.16) and

()- [M () + IN(: -( (y (2.17)

From Eq. 2.17 we obtain



:- [N -1 (N [NI -1 [M] (ex (2.18) and, after substituting from Eq. 2.18 in Eq. 2.16, we obtain


[P- I( ) + I [N- (2.19)












22

where

[P]- [L] - [M] [NJ- [M


Eq. 2.19 defines the thermal force and moment, pT and MT respectively, as


- - [eM [ N1( 1( (2.20)



From now on, the thermal force and moment resultants defined in Eq. 2.20 are considered as external forces and moments acting on the composite specimen (Agarwal and Broutman, 1980). The finite element method is used to solve Eq. 2.20 and to obtain the solution for E0x(x) and Kx(x). For bottom element the derivation is the same, except for the limits of integration in Eq. 2.5, 2.8, 2.11 and 2.12.



2.4 Free Edge Problem

In defining the in-plane and transverse displacements for the free edge problem (second phase), Eq. 2.1 remains the same, except for displacement in the y direction, which takes the form of

v(y,z) - e0y (2.21)











23

where eo represents the applied strain in the y direction. Consequently, the strain E in Eq. 2.2 is replaced by e, - eo (2.22) Because the dimension 2b in Figure 2.1 is considered much greater than the dimension 2L, it is assumed that a state of plane strain parallel to the x-z plane exists in the second phase of the free edge problem. The stress-strain relations for the kth layer in this case are



a 1,Q Q22 0 *e ej (2.23)



The in-plane force and moment acting on the beam are shown in Fig. 2.6





V








2M M




Figure 2.6 Force and moment resultants (top sublaminate)












24

P, M and V are defined as

h
[PM, V]- f [OZOxx,Tx,,,] dz (2.24)
0

or, using the definitions given in Eq. 2.8 and 2.23, as



S55 b1 dl 0 eo + MI
0 ass Y x 0


where PP and MP are the force and moment induced by the applied strain eo. They are defined as

h
[PP,MP] - f, (1, z) dz - [a12,b12] e0 (2.26)
0

The thermal force and moment resultants defined in Eq. 2.20 are added to the problem in the form of external forces, and hence the problem takes the final form of


p p a ll by, O0 p
M - b d11 0 eo0 M (2.27)
V 0 0 ass Y xz 0











25

or


P all bl0 e" PT PP
d1 eoj - I -IMIf (2.28)
S0 a55, Yxz k,0 0 The last equation can be expressed in a more concise form as

F- Se - R (2.29) where S is the stiffness matrix, components of which are (al bl 0

[S] - lb1 d1 0 (2.30)
0 0 a55ss)


and R is the vector that represents the forces and moments induced by the applied strain and by the thermal effects



(R)- M - M (2.31) The vector F refers to the external forces acting on the beam and its components are


(F) - ( (2.32)












26
while e is the vector of deformations in the form of



(e)- Kx (2.33)


Finally, the Equation 2.27 can also be expressed as

e - c(F+R) (2.34) where c = s"' and is the compliance matrix.



2.5 Finite Element Solution

The Finite Element Method is used to solve Eq. 2.29, and a beam finite element with offset nodes is selected to model the problem. The derivation of this specific beam finite element was done by Sankar (1989). The finite element has two nodes and three degrees of freedom (u,W and w) at each node. The nodal forces and moment for the ith node are denoted by pi, mi and vi. In the present study a slightly different method is used to derive the stiffness matrix.

The average force and moment resultants (P, M and V) at the center of the element are given by ( p2 p9
(m2 - ml) / 2 (2.35) (V, - v() / 2












27

The deformations at the center of the element are expressed in terms of the nodal displacements as Ca, ' (u.-ul) /L
(e)- - /L (2.36) Yzx, T2 ) /2+ (w,-w) / where L is the element length.

The beam finite element forces and degrees of freedom are shown in Fig.2.7. By substituting from Eq. 2.33-2.34 into Eq. 2.29, the relations between nodal forces and displacements are obtained in the form of f - k q - R (2.37) where the vectors of nodal forces (f) and displacements (q) have the following components

(q) r. [u, l 11 W1, u212 2w2] (2.38)
(f) T- [P1, 1, v1p2,m, V2]

The vector R and its components are defined by Eq. 2.20,

2.26 and 2.31. The elements of the stiffness matrix k are











28

all b 0 all b 0 L L L L
by, d11 asL) ass b (_ d + asL) a L L 4 2 L L 4 2 Sa55 a55 0 a55 as5 2 L L L [K]
Sa1 b a11 bL L L L
b1_ d1 +_ a5sL) L as5 b, d a5sL a55 L L 4 L L L 4 2 0as55 a5 0 a55 a55 L L 2 L (2.39)


2.6 The Finite Element Models

Fig. 2.3 and 2.4 show the models used in modeling the free edge delaminations. They consist of two sublaminates above and below the delamination plane. The assumptions are that the delamination plane is the weakest plane in the material and that there is no frictional resistance against sliding of the sublaminates. The first assumption means that the crack will propagate along the delamination plane, without crack branching.

In the uncracked region, the top and bottom portions of the beam are connected by common nodes, when calculating the strain energy release rate (Fig. 2.3), and by spring elements, when predicting the stress ahead of the crack (Fig. 2.4).










29




V 2





Mv P
1





1p

I1 - 2 x
(a)



z

w
w1 2





1
u. 1 2 u2



(b)




Figure 2.7 Force and degrees of freedom
(a) Nodal and average forces (b) Degrees of freedom
(a) Nodal and average forces (b) Degrees of freedom











30

2.7 Strain Enerqy Release Rate

In the calculation of the strain energy release rate (G), one rigid element is placed in the tip of the crack connecting both sublaminates (Fig. 2.3). The rigid element has nine degrees of freedom. Apart from the three displacements at each node, the three generalized forces, Ft, transmitted by the rigid element are also considered as unknown degrees of freedom. The constraint equations corresponding to the rigid element are
Ff 03 03 13O
, - 03 03 I1

0 3 0z F (2.40) where F, and F2 are the nodal forces at nodes 1 and 2 and ql and q2 are the nodal displacements. The symbol 03 represents the 3x3 null matrix and 13 the 3x3 identity matrix.

Considering a zero volume path surrounding the crack tip as shown in Fig. 2.8, the expression for the J-integral relative to two-dimensional states can be found in Hellan (1984) as

au(
S- ( Wdx2-Pi ds) (2.41) r 1a












31

where

Pi - n Gji i-1,2

For the present case, and considering the path r very thin (zero-volume), the Eq.2.41 takes the form of


J- f dz - fnxxxU,xds- fnxrzxWxds (2.42)
r r r

where W is the strain energy density and is defined as


W - i (oXXeX + oyeY + zxYzx) (2.43)


The path r has four segments, one in each sublaminate behind and ahead of the crack tip (Fig. 2.8). The values of nx (unit vector normal to the path) and ds for each element are:

Elements 1 and 4: nx =-1; ds =-dz Elements 2 and 3: nx = 1; ds = dz

Considering the Element 1 in Fig.2.8 we have:

h h h
J1 - -f dz + foudz+ dzZ (2.44)
0 o 0

where the subscript stands for the element number.

The first term of the right hand side of Eq.2.44 is defined by using Eq.2.43 and the stress-strain relations given in Eq.2.23 as

h h (2.45)

fdz - 1e," + 2 Q12exe +22e02 +55 zx d
0 0









32



L1,4 - 2,3

z .e.e.ment.I.ement2..





telement4 element 3J

elements ahead elements behind of the crack tip of the crack tip





element 1 2 21 element 2 v 1 2V


element 4 47 J4 33 3 F3 element 3







Figure 2.8 Zero-volume path for J-integral












33

From the definition of the stiffness coefficients and Eq.

2.30 and 2.33 we obtain the following relation

h
eTSe - f(11exX2+ssyz2) dz (2.46)
0

Therefore, the Eq.2.45 can be modified and expressed as

h
Wdz- Se + a12ee + a22 e2 (2.47)
0

The second term of the right hand side of Eq.2.44 is modified in a similar way to what has been done for the first term.

h h
foju,xdz-f (D1iexx2+i12exxe,) dz (2.48)
0 0

or

h h
OXXuxdz- f(Uzzex 12 xx oD55 x 255 zx2) dz (2.49)
0 0

Finally,

h
foxxuxdz- e1 TSe1 + a12 e0exx- a55Yzx2 (2.50)












34

The third term of Eq. 2.42 is

h
rzxw,xdz- ,V i(2.51)
0

where V1 and qi are the shear force resultant (Eq. 2.24) and the rotation at the crack tip, respectively, for the element number 1.

Substituting Eq. 2.48, 2.50 and 2.51 into Eq. 2.44 we obtain


J- 1elTSel + Uc (2.52) where Uc is a energy term defined by


Uc - 1 a22 e02 - aszx2 (2.53) Similar expressions can be obtained for the other elements in the form of


2 -- - e2 rSe2 - U�

J3 - - e3 TSe3 - Uc (2.54) S-e4TSe4 + U

For linear elastic fracture, the J-integral is equivalent to the energy release rate (Hellan, 1984). This is valid for fractures in brittle composite materials, where the plastic zone in the vicinity of the crack is negligible. Hence,












35

G - J- J + J + J + J4- (e,'sel + e4TSe4 - e2 Se2 - e3Se) (2.55)

The Eq.2.55 can be calculated in a more convenient form,


G-1 (FiTc F + F4TcF4 - F2 ctF2 - F3TcbF3) (2.56) where the forces Fi (i=1,2,3,4) are calculated from the finite element solution, and so are the compliance matrices ct and cb from the input data.

The elements 2 and 3 are in the uncracked region, and so, they are under the same deformations. We thus obtain the relations,

e- e3
CtF2 - CbF3

Substituting the last equation, and also the equilibrium relations (F1 + F4 = F2 + F3), in Eq. 2.56 we obtain G-1 (FiTct - F4Tcb) (F1-F2) (2.59) or, using again the equilibrium relations,


G-1 (F-F2) T(ct+Cb) (F1-F2) (2.60)











36

Actually, (F1 - F2) is the force transmitted by the rigid element, Fr, between bottom and top crack tip nodes. Therefore, the strain energy release rate (G) can be expressed in a more convenient form as


G - fFT(Ct+Cb) F (2.61)


where Fr are considered unknown degrees of freedom and part of the finite element solution of the problem (Eq. 2.40).

One can also evaluate the value of the strain energy release rate using a method analogous to the virtual crack closure method in two-dimensional fracture problems ( Hellan, 1984). In this method, the energy necessary to "close" the crack is equal to the energy release rate.

In Fig. 2.9, as the crack propagates, the crack tip Nodes 1 and 2 have displacements A1 and A2, respectively, and the crack tip moves to Node 5. If the length of the elements behind the crack tip is sufficiently small, AL, then the crack opening, (Aq,-Aq2), will be approximately equal to the difference in displacements of Nodes 3 and 4,(Aq3-Aq4). Therefore, the problem is analyzed considering the displacements of Nodes 1 and 2 as the tip displacements of cantilever beams of length AL under applied forces Fr on each beam. For very short cantilever beams (thickness/length <<1)










37
AL




3 ----------2





---------cantilever beam
F

r










cantilever beam LF Figure 2.9 Crack-closure method












38

the tip deflections (qti) can be calculated using the expression

tip - L [c] F (2.62) where L and F are the beam length and the forces applied at the tip of the beam, respectively, and c is the compliance matrix. In the present case, Aq,- ALctFr (2.63) Aq - -ALcbF

or

L (c + cb) F. (2.64)


Substituting the last equation in Eq. 2.61 gives


G 1- FrT(A I Aq2) (2.65) 2AL

In this equation, the displacements in the tip nodes and the forces in the rigid element are part of the finite element solution.

In the present study, the strain energy release rate was calculated by three different methods. These methods are represented by Eq.2.56, 2.61 and 2.64, and the results obtained from the three equations are in very good agreement.











39

2.8 Stresses Ahead of the Crack

In evaluating the stresses ahead of the crack tip, the uncracked region was connected by discrete spring elements (Fig. 2.4). The stress ahead of the crack in the ith element was obtained by


i ALb (2.66)


where V and L are the vertical force in the spring element and the length of the element, respectively. The vertical forces are obtained from the Finite Element solution of the problem, and the element lengths are part of the input data.

The determination of the spring constant for the spring elements is of significant importance in appropriately modeling the problem. The spring constants are derived from the change in thickness of a beam subjected to transverse forces. The exact change in thickness is obtained from the two dimensional elasticity theory (Timoshenko, 1970). In the first step, just the top part of the beam is considered and it is represented in Fig. 2.10. From the elasticity theory we have, for this case


y 3q 1 y3- 2y+ 2 c3) (2.67) 4c3 E* 3 3








40









-------------- -- L











igue 2.10 defor spring constant Calculation
~ on












41

where E* is defined as


E*- C3 1-v12 21 (2.68) E EE2A


and


A -V 1 V1-v23v32 -31V13 -221V32V13 (2.69) E E2 E3


The vertical displacement (v) is obtained by integrating the deformation E in the y direction. The result is


v- eydy-- 3q ( -Cy + c3y) + C1 (2.70) S 4c3E 12 2 3



where C1 is a constant.

From Eq. 2.70 the displacements at the positions y=-c and y=0 are determined as 13 qc (2.71) 16 E*
v(0)- C1

The relative displacement between these two particular positions is


VR- v(-c)-v(O)- 13 qc (2.72) 16 E*











42

For a beam with unit area we have that P = q. In the present case, c is equal to half the thickness of the beam (c = h/2). Therefore Eq. 2.72 is rewritten as 32E*v (2.73) 13h

For a spring model we have

P- k VR (2.74) where k is the spring constant.

By comparing between Eqs. 2.73 and 2.74 we have the value for the spring constant per unit length considering just the top part of the element.


k- 32E* (2.75) 13h

Considering both parts, top and bottom, we have two springs placed in a series arrangement (Fig. 2.10). The spring constant equivalent (kq) to this kind of arrangement is given by


ktkb
k ktkb (2.76) kt+kb











43

where the subscripts t and b stand for top and bottom respectively. Substituting Eq. 2.80 in the last equation we obtain the final expression for the spring constant used in the present study

32 E*
k 32E* (2.77) e13 (ht+hb)


2.9 Numerical Results

Some results from the investigation of strain energy release rate and the interlaminar stresses in a free edge specimen can be found in Whitney (1987), (Pagano, 1989) and Raju (1986). These researchers computed results for the same composite specimen. The present study uses the same specimen to make the comparison easier.

The free-edge specimen width to thickness ratio (b/h) is 25 (Fig 2.11-2.16), which is typical for conventional edge delamination test (EDT) specimens. Considering a coordinate system where xi and X2 are parallel and transverse to the fibers respectively, and x3 is through the thickness, the ply properties have the following relations and values:

Ej/E2 = 14 E2/E3 = 1

G12/E2 = 0.533 G23/E2 = 0.323

12 = 0.3 V23 = 0.55
a = -9 E-07/OC a2 = a3 = 23 E-06/OC S= 0 2 = S = 5560 E-06/%












44

where Ei is the elastic moduli in the ith direction, Gij is the shear modulus in the xi-x plane and vil are the Poisson's ratios. Finally, the thermal expansion coefficients and the moisture swelling coefficient in the xi direction are denoted by ai and fi, respectively. These properties and relations among them are typical of current high performance graphite/epoxy unidirectional composites.

The specimens considered in the present study have the following lay-ups:

Specimen A: [903/03]s Specimen B: [-60/602/-60/02]s Specimen C: [-45/452/-45/02 s Specimen D: (02/-60/602/-60]s

A reasonable stress-free temperature for standard graphite/epoxy laminates is 1380C (2800F) (Whitney, 1987). Assuming 210C (700F) as the room temperature, the difference in temperature (AT) that causes residual stresses in the plies is equal to -1170C (2100F).

The critical axial strains (eo) experimentally determined in edge delamination tests are 0.3% for angle-ply laminates, and 1% for the special case of a bidirectional class of laminates (Whitney, 1986). Therefore, these values of strains are used in the numerical examples.



2.9.1 Strain energy release rate

The normalized strain energy release rate, G,, is given by












45

G
GN " e (2.78) E,h ea

In the present study the energy release rate was calculated using Eqs. 2.56, 2.61 and 2.65. The results obtained from these three methods were dn good agreement and they are refered to as "Present" in Table 2.1. Comparison of the results from the present study with available results is also present in Table 2.1. The relative differences between the results from the offset-node method and each of the other three approaches are shown in Table 2.2.

From Tables 2.1 and 2.2 it may seen that the present results compare very well with Raju (1986), whereas the comparison with the other two methods is not good in a few cases. It should be mentioned that of the three methods we compared with, Raju (1986) used a detailed finite element analysis.



2.9.2 Stress ahead of the crack

The interlaminar stresses ahead of the crack tip are evaluated using spring elements in the uncracked region of the specimen. The spring constant for these elements is calculated by Eq. 2.77 and found to be equal to 300x 109 N/m .











Table 2.1 Comparison of strain energy release rate values obtained from different methods.


LAMINATE (%) AT WHITNEY PAGANO RAJU PRESENT 0.3 0. 1.300 1.300 1.2506 1.3000 (-60/602,-60,02)s 0.3 -117. 1.700 2.000 1.9551 2.0580 0.3 0. 1.300 - 1.2700 1.3000 (02,-60,602,-60) 0.3 -117. 1.000 - 0.6700 0.7160 1.0 0. 0.016 0.016 0.0155 0.0156 (903,03) s 1.0 -117. 0.077 0.045 0.0605 0.0650 0.3 0. 0.770 0.770 0.7285 0.7440 (-45,452,-45,02)s 0.3 -117. 0.800 0.980 0.9194 0.8640









Table 2.2 Relative difference (%) between the results from present method and
other methods.
Relative difference=100 x (present method - other's) / other's


LAMINATE A (%) AT WHITNEY PAGANO RAJU 0.3 0. 0.0 0.0 4.0 (-60/602,-60,02), 0.3 -117. 21.6 2.9 5.3 0.3 0. 0.0 - 2.4 (02,-60,602,-60)s 0.3 -117. -28.4 - 6.9 1.0 0. 2.5 2.5 0.6


(903, 03)s 1.0 -117. -15.6 44.4 0.0 0.3 0. -3.4 -3.4 2.1


(-45,452,-45,02) s 0.3 -117. 8.0 -11.8 -6.0


" -.











48

Fig. 2.11-2.15 show the results for peel stresses for different laminate configurations, with and without thermal effects. From the results we find that the effect of residual thermal stresses is to cause increase in the maximum peel stress at the free edge. The effect is significant for specimen A ( [903/03] ) and less significant for specimen C ( [-45/452/-45/02]s ). In Fig. 2.11 the stress distribution is compared with Whitney (1987) for A= OOC and A= -1170C for the specimen B ( [-60/602/-60/02]s ). The agreement is excellent. For the bidirectional laminates, specimen A, the influence of temperature can be better appreciated in Fig.2.15, where the temperature variation goes from zero to -1200C in increments of -300C.

Besides thermal stresses, some studies on effect of moisture absorption on the free edge stresses were also performed. In fact, the analysis procedure is similar except that a AT has to be replaced by 8 AH, where 8 is the swelling coefficient and AH represents the variation in the moisture concentration. Fig. 2.16 shows the influence of the moisture level variation on specimen A. The moisture level was allowed to vary from 0% to 1.2% of the laminate weight, which reflects feasible conditions. This result agrees with those obtained by Armanios and Badir (1989). Unlike the thermal stresses, the residual moisture tends to decrease the interlaminar stress close to the crack front.












49

2.9.3 Conclusions

The residual thermal stress tends to increase the total energy release rate and the interlaminar stress at the tip of the crack, while the residual moisture alleviates these effects.

The results obtained from the offset nodes method are very close to those obtained by Whitney (1987), Pagano(1989), Raju(1986) and Armanios and Badir (1989) using different methods. Therefore, the present offset-node method is a good method to predict the energy release rate and stresses ahead of the crack tip in the free-edge problems. The present finite element method has the facility to model any type of geometry, loading and boundary conditions, and avoids renumbering the nodes and elements as the crack propagates.

The offset node method can be used in studying dynamic delamination, and can be extended to problems of sublaminate buckling and delaminations in composite plates.









0.7


0.6
2h




o 0.4
0.4 b/h=25
-, \A CT=-1170C e0= 1% b0.3- __ _T= oC e o Whitney's result : & T= -1170C
0.2- * Whitney's result : A T= OOC i Present result

0.1

0
z 0


-0.1


-0.2 I I I I I
0 0.4 0.8 1.2 1.6 2 2.4 x/h

Figure 2.11 Stresses ahead of the crack tip for specimen A 0









3

Zt
2.5 2h

ILp- b 2
2 A T=-117 C [-60/6 /-60/0s ST= 00 C b/h=25 e = 0.3%
1.5 0









05







x/h

Figure 2.12 Stresses ahead of the crack tip for specimen B















1.8 0



a [-45/4!V-45/%4j I- 1.2 AT=-111C blh=25 AT= 00oC e = 0.3%


0.8

w 0.6


E0.4 EO
s
0
Z 02


0

0.2



0 .4 0.8 12 1.6 2 2.4
x/h

Figure 2.13 Stresses ahead of the crack tip for specimen C








3

z

2.5 -h 2h
b x



2 /-60/64/-60]s ep A T=-1170 C b/h=25

1.5 A T= 00 C eb=-0.3%
N


1



E

0




-0.5
0 0.4 0.8 1.2 1.6 2 2.4 x/h

Figure 2.14 Stresses ahead of the crack tip for specimen D










0.6


0.5
2h
x

0.4 - AT= -120C b L\jT= -90 93 u(C ZS, T= - 60 C: b/h=25
0.3
A T= - 30 CC 0 AT= OC
0.2


� 0.1


Z 0



-0.1



-0.2
0 0.4 0.8 1.2 1.6 2 2.4
x/h
Figure 2.15 Influence of temperature variation









0.3


0.2


0.1






-0.1
L /H= 0.0% = -0.2 AH = 0.3% AH = 0.6% H = 0.9% 2h z -0.3 - H = 1.2 % 0 b x


-0.4 - 9/ b/h=25
e = 1%

-0.5 I I I I I I I I
0 0.4 0.8 1.2 1.6 2 2.4 x/h
Figure 2.16 Influence of moisture level variation















CHAPTER 3

MODELING DELAMINATION IN COMPOSITE BEAMS
AND APPLICATION TO TORSION PROBLEM


3.1 Introduction

The application of composite materials is increasing day by day. Not only traditional areas like aerospace and automobile structures, but also several other areas, such as robotics, marine industries and medical devices and prostheses, have shown that the application of composite materials is major factor in the improvement of their mechanical structures.

The advances in the manufacturing processes of composite materials, the development of new matrices and fibers, and the progress in methods for study of the mechanical behavior of composites have been made it possible for the composite materials to replace traditional materials in several mechanical structures. Usually composite materials are considered as a possible and adequate solution whenever the weight savings and the capability of "tailoring" the material in according to specific needs are more critical than the costs.

As the use of composite materials in structural applications increases, there is more need for structural


56









57

analysis. Unlike an isotropic beam, stretching a laminated composite beam can produce bending, shear and twisting in addition to the extensional deformation. The coupling among the modes of deformation is the major obstacle to precisely modeling the mechanical behavior of beam like composite structures.

In the present chapter, a new and efficient method to analyze delaminations in laminated composite beams is presented. The equations of equilibrium were derived by Sankar (1991) from the classical shear deformable laminated plate theory (Whitney, 1987). The equilibrium equations are assumed to be satisfied in an average sense over the width of the beam. From this assumption results a new set of force and moment resultants, which amplify the possibilities for modeling beams under different loading conditions, including torsion. The equilibrium equations are derived from the Minimum Potential Energy Principle.

In the present study a new offset beam finite element is developed for modeling the problem. Initially, the capability of the new finite element in modeling the problem is verified by solving the problem of a cantilever beam under different end loading conditions. The results are compared with those from beam theory solutions found in Timoshenko (1970), Reismann and Pawlik (1980), Whitney (1987) and Sankar (1991). Next, the strain energy release rate for a delaminated specially orthotropic cantilever beam under torsion is









58

calculated. The delamination is considered to be at the middle-plane of the beam. The results for G from the finite element method are compared with analytical solutions of the problem. The effect of position of delamination on the strain energy release rate values is presented. Finally, the problems of a specially orthotropic and a quasi-isotropic simply supported delaminated beam under a transverse force is considered. The strain energy release rate values for various positions of the delamination are computed by using the finite elements developed in the present study. The results will be useful in explaining delamination propagation in composite beams due to quasi-static impact.



3.2 Beam Equations

Considering the laminated beam shown in Fig. 3.1, the Shear Deformation Theory (Whitney, 1987) for laminated plates assumes the displacement field as

u(x, y, z - uo (x, y) +z* x(x, y)
v(x,y, z)-vo x,y) +z (x, y) (3.1)
w(x,y, z) - WO (, y)

where uo, vo and w0 are the displacements in the middle plane of the plane in the x, y and z directions, respectively, and so, are yx and (,y with respect to the rotations about the X and Y axes.

All five displacements and rotations in equation 3.1 can be expanded as a Taylor series, in the direction of the width









59


























L h Figure 3.1 Laminated composite beam









60

of the plate. If just the first order terms of this series are retained, the displacements and rotations can be expressed in the form of
u0 (x, y) - U(ix) +yF(x) v (x, y) - V(x) +yG(x) w (x, y) - W(x) + ye (x) (3.2) x (x, y) -(x) + yaO (x)
S(x, y) - (x) +y (x)

The displacement field is now represented in terms of ten functions of the x coordinate.

u(x, y, z) - U(x) +yF(x) + z (x) +yza (x) v(x, y, z) -V(x) +yG (x) + zY (x) +yz P (x)
w(x, y, z) - W(x) + ye (x)

From the beam theory, we assume that there is no deformation in cross sections normal to the x direction (ezz=0, E=0 and yyz=0). From this assumption and Eq. 3.3 results e, - G(x) + zP(x) - 0 (3.4) and

Yyz-9 (x) +yP (x) +0 (x) -0 (3.5)

From equation 3.4 and 3.5 we obtain that G(x)-0
S(x) - 0 (3.6)
(x)--0(x)

The displacement field can now be represented in terms of seven functions of the x coordinate, instead of the ten functions used before.

u (x,y, z) - U(x) +yF(x) + z (x) +yza (x)
v(x, y, z) - V(x) - zO (x) w(x, y, z) - W(x) + yO (x)









61

Denoting the differentiation with respect to x by a prime, we represent the deformation field from Eq. 3.7 as


x U'+yF' #'+ya ex
Y F+V' +z a-O' - Yxy +ZK
Y \W) +y (a+) 0 a (3.8)


The deformation field can also be expressed as

E-E+yE (3.9) where


ex O U( F'

(E)- {Kx; ) I (3.10) Yx, a a0, From the constitutive relations the stresses developed in the kth layer are given by


()- [ e-xo)+z [1~C (3.11)


and


T - Q55 ( + W) + Y s (a +') (3.12) where


1 12 D16


.016 026 066









62

The terms Qij in Eq. 3.13 are stiffness terms and can be found in Whitney (1987).

By definition, the force (N) and moment (M) acting on the laminate are

h

(N, Mx) - f a,[ 1, z] dz
h
2
h


2
vX- f "r XYdz
h


2
x xzdz


and they are shown in Fig. 3.2.

The stiffness coefficient terms are defined as

h

(Ai, Bij,Dij) - f ij [1, z, z2] dz (3.15)
h
2

By using equations 3.9-13 and the definition of stiffness coefficients, the Eq. 3.14 can be written in a more simplified form as


F- C ( E+yE) (3.16) where

(F) T- (Nx, N (3.17) X, 'Y'MX M'' (317









63










Z








































Figure 3.2 Force and moment resultants









64

and the matrix C is formed with the following stiffness coefficient terms A1 Ag B,11 B,16 0 A16 A66 B16 B66 0 [C] - B11 B16 D11 D16 0 (3.18) B16 B66 D16 D66 0
0 0 0 0 kAss5 The factor k presented in term C,5 is the shear correction factor (Whitney, 1987), and it is considered equal to 5/6 in the present work.



3.3 New Set of Force and Moment Resultants

A new set of force and moment resultants is defined

from the integration of the column vector of forces along the width of the beam:

b

F- Fdy (3.19)
b
2


and

b

- fFydy (3.20)
b
2

In the two last equations, the vector of forces F is given by equation 3.17, and b is the width of the beam.

Substituting equation 3.16 in Eq. 3.19-3.20, and assuming that the stiffness C of the beam is constant, we obtain









65

b/2
F-_f C Edy-bCE
-b/2 (3.21) /2 CEy3dy( )CE

-b/2

The explicit beam contitutive relations are:

SX Azz Az6 B z B16 0
N A16 A66 B16 B66 0 V F (3.22)
[] -b " B16 I, 1D16 0 V
MxY B16 B66 D16 D66 0 ax, 0 0 0 k2As55 and


SAl A16 B11 B16 0 xA16 A66 B16 B,66 0 F
0 v1 B, B,, D D,, o 0(323)
B16 B66 D16 D6 0 0
0 0 0 0 k2A5,



3.4 Equilibrium Equations

The Principle of Minimum Potential Energy is applied in the derivation of the equilibrium equations. The principle says that of all displacement fields which satisfy the prescribed constraint conditions, the correct state is that which makes the total potential energy of the structure a minimum (Tauchert, 1970). The total potential energy of the structure (ir) is obtained from the sum of the strain energy of the beam (0) and the potential of the external force (X):









66

" -4 + X (3.24) According to the Principle of Minimum Potential Energy, 6I - 6# + 6x - 0 (3.25) where the symbol 6 is the variational operator.

Considering only transverse loading acting on the beam surface in the z-direction, the potential relative to the external loading is given by

b
L 2
x - -f f q(x,y) w(x,y) dydx (3.26)
0 b
2

where q(x,y) and w(x,y) are the transverse loading and displacement, respectively, and L is the beam length.

Similar to the procedure adopted to define a new set of force resultants in Eqs 3.19-20, the transverse loading is divided into two parts and defined by q(x,y) - q(x) + y 2(x) (3.27) where

b
2
(x) - q(x, y) dy
b
2 (3.28) I(x) - f q(x, y) ydy
2









67

Substituting the definition of transverse displacement (Eq. 3.7) and transverse loading (Eq. 3.28) in Eq. 3.26 gives the result

L
X - - f((x) W(x) + <(x) (x)) dx (3.29)
0

From Eq. 3.16, the strain energy per unit area (OA) of the beam can be derived as


4 - 1 ETCE
2 (3.30)
- (E+yE) TC(E+yE)
2

The strain energy per unit length (OL) is obtained by integrating pA along the width (b) of the beam

b
2
1L Af A dy
b
2 b (3.31)
1 f (E+yE) TC(E+y dy
b
2

Because

b
2

f Er'CEydy- 0
2 (3.32)
b


2









68

the Eq. 3.31 takes the form of

b b

-lf (fC) dy+ f ( TCr) y2 dy
b b (3.33)
2 2
1- (b!TC + b3 fTC
2 12

The total strain energy in the beam is expressed as

L
S- foLdx
o (3.34)
1 b
- 0 bETCE + --E9C dx 2Jk 12

The total potential of energy in the whole structure (w)

is obtained by substituting the Eq. 3.29 and 3.34 into

equation 3.24.


S bETC+ c-(x) w(x) (x)(x) dx (3.35)
2 12

Applying the Principle of Minimum Potential Energy we

obtain the following equilibrium equations


an dNx
-0
8U dx
6 s d Tx
-0
6V 8x
+ -0
6W dx
6 - d - (3.36)

6F "x dx0 8Tn d ( #x- Nxy)

-6 -o
SO dx 87t dMx
8o dx a c dx X Ry 8a Xr 7fxx








69

Next, from each equilibrium equation in Eq. 3.36 results one corresponding natural boundary condition, which is shown in the next equations as


d -N -0 Rx6U- 0
dx
d xo -0 x, 8 v- 0 dVx
dx
d x (3.37) d x - 1x 6 F- 0 d( -N- 60
d ( V- TXY)-o
dx --4 x-Mxy) 60 - 0 dx
dx-x M x Ixb# - 0 dx



By using Eqs. 3.22-23, the new set of force and moment resultants can be expressed in terms of seven unknown functions: U, V, W, F, 0, a and 0. Substituting for the resultants in terms of displacements in the equilibrium equations 3.37 results in a system of seven ordinary differential equations for the displacements.



3.5 Finite Element Solution

The system of ordinary differential equations obtained from the equilibrium equations 3.37 is solved in this work using the finite element method. A new finite element has been developed to model the beam problem. This finite element has three nodes and seven degrees of freedom (U, V, W, F, 0, a, 9) at each node. The middle node is statically condensed









70

when solving the problem for displacements, but it is considered when calculating the strain energy release rate, in order to obtain a more accurate solution. The nodal forces and

moments for the ith node of the structure are Fx, Fy, Fz, Mz, Mx W and T.

A quadratic variation of all seven displacements is assumed along the element length. Denoting by X any specific displacement, and by Xi this specific displacement at the ith node, the displacement and deformation within the element are defined as

X- a X+ a2X2+a3X3 (3.38) X'-bX1+b2X2 +bX3

In the last equation, the terms a1, a2 and a3 are interpolation functions in the form of a1 - 1-3x+2I
a2-x(2x-1) (3.39) a3 - 4x(1-2x)

and bl, b2 and b3 are their derivatives with respect to the x direction

b, - (4x-3)/L
b2- (4x-1)/L (3.40) b3-4 (1-2x)/L








71

where

x-x/L

In the last equation, L is the element length and x is the coordinate along the beam axis.

By using Eq. 3.38, and the terms a. and bi (i=1,2,3) defined in Eqs 3.39 and 3.40, the deformation field (Eq. 3.10) within each finite element as a function of the nodal displacements is expressed in the form of

()- [a] (q) (3.41)
(T) - [9] (q)

where q is the vector of the displacements

[ q] T' [U1, V, W, F, II, a1 1,
U2, V2 W2, F2" 2, "2, 2, (3.42) U3, V3, W3, F,31,a33,8]

and & and a are two strain-displacement matrices, both with dimension of (5 x 21).

The first strain-displacement matrix, a, is composed by the following terms:


'b1 0 0 0 0 0 0 b2 0 0 0 0 0 0 b3 0 0 0 0 0 0 0 b 0 a, 0 0 0 0 b2 0 a2 0 0 0 0 b3 0 a3 0 0 0 0 0 0 0 by 0 0 0 0 0 0 b2 0 0 0 0 0 0 b3 0 0 0 0 0 0 0 al -bi 0 0 0 0 0 a2 -b2 0 0 0 0 0 a3 -b3
0 0 bi 0 al 0 0 0 0 b2 0 a2 0 0 0 0 b3 0 a3 0 0 (3.43)









72

and the second one, A, of:

0 0 0 b0 0 0 0 0 0 b2 0 0 0 00 0b3 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
000 0 0 b 0 0 0 0 0 b2 0 000 0 0 b3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 00 0 0 al bi 0 0 0 0 0 a2 b2 0 0 0 0 0 a3 b



(3.44)

Similarly to Eq. 3.10, the total deformation within each element is given by the expression e-e+yd (3.45) or, after using the strain-displacement matrices, by e- (a+yg) q (3.46) With

a-;a+yd (3.47) the element stiffness matrix is defined as

b

K - f f aDa dy dx O b
b (3.48)
- ff (a+yd) TD (a+yd) dy dx
0 b
2









73

In Eq. 3.48, D is the elasticity matrix for composite materials. The terms in matrix D are the same terms shown in Eq. 3.18. Each term of the matrix D is defined as

b
(Aj" B,IDi) - i, 11z, z'2] dz (3.49)
a

where a and b are respectively the lowest and highest coordinates in the z direction in the considered element. By defining two element stiffness matrices in the form of

L
K-bfaTDa dx
o (3.50) Kb3f TD f dx
0

the equation 3.48 is reduced to K- K+ K (3.51) The elements of the stiffness matrices K and K are shown in the Appendix. The dimension of the stiffness matrix k are reduced from (21 x 21) to (14 x 14) by using static condensation and the stiffness matrix of each finite element is assembled in a global stiffness matrix, representative of the stiffness of the whole structure.

The final representation of the problem is

F - K q (3.52)









74

where K is the total stiffness matrix, and F and q are vectors with the nodal forces and displacements, respectively. In the present work, the equation 3.52 is solved using the Gauss Elimination Method (Bathe, 1982).



3.6 Strain Energy Release Rate

The strain energy release rate is calculated by using the expression for J-integral. A zero volume path is delineated surrounding the crack tip, as shown in Fig. 2.8, and the Jintegral expressed in terms of the strain energy per unit length (U) of each element around the crack, as in Eq. 2.55, j - U(1) + U(4) U(2) _ U(3) (3.53) The strain energy per unit length of the beam (Eq. 3.34) is given by


U 1 bETD + b_ D9~ (3.54)
2 12

The expression for the strain energy per unit length of the ith element is obtained by substituting Eq. 3.46 in the previous equation.


U) q biTD + _ DS) q (3.55) 12









75

In Eq. 3.55, the vector q contains the displacements of the nodes of the ith element, including the middle node. The reason for considering the middle node in the calculation of the strain energy is because the static condensation process assumes that no forces or couples are applied to the intermediate node. This assumption is not true for the case where bottom and top offset elements are put together to model the beam.



3.7 Rigid and Gap Elements

The possibility of interference between the crack surfaces in the delaminated region, when the beam is under loading, principally under torsion loading, requires that gap elements or rigid elements be placed in appropriate positions in order to monitor the contact.

Two types of rigid elements are anticipated to be used depending on which side, left or right, of the beam presents the interference phenomenon. Figure 3.3 illustrates the interference situation and the rigid element actions.

In matrix form, the equations of the rigid element placed on the left side can be represented as









76

0 0 0 0 1
0 0 0 0 b Wt: (Fz) t
2 () Tt
00 0 0 -1 Wb (3.56)

0 0 0 0 -- e.
b

b b sF, O 1 b- 1 0 0
2 2

where t and b represent top and bottom elements, respectively, and the force F is shown in Fig. 3.3.

Similarly, the equations of the rigid element placed to on the right side are

0 0 0 0 1
0 0 0 0 --_:t 'Wt TF 2 e Te
o 0 0 0 - Wb - (F,)b (3.57)

2 b Tb
1 - 0
2 2


The gap element is placed in the structure to avoid the interference in the delaminated region while the structure is under vertical loading. In matrix form, the constraint equations relative to the gap element are:


0 0 1 ' ( (Fz) t (3.58)

1 -1 0) F)(0 O









77

crack


M





(a)
bF

-F
F
zIp- b --.


y rigid element -
bF
F

F
Interference rigid element action
(b)


bF
F b


rigid element

yF
-F bF interference rigid element
(c) action


Figure 3.3 Rigid elements: (a) Beam under torsion

(b) Interference on left side (front view) (c) Interference on right side (front view)









78

The rigid element and gap element matrices in Eqs. 3.563.58 are assembled in the global stiffness matrix using the same procedure adopted to assemble the element stiffness matrices. In Eqs. 3.56-3.58 Fz is the contact force between the interfering elements, Tt or Tb is the couple due to the contact force and F is the force in the rigid or gap element in the z direction.














CHAPTER 4

NUMERICAL RESULTS AND DISCUSSION



4.1 Numerical Results

This chapter presents the numerical results obtained from the finite element described in the previous chapter. Three different tests were performed, and their numerical results compared with those obtained from the beam theory.

The first test of the beam finite element consists of a specially orthotropic cantilever beam subjected to several different end loading conditions. This test is divided into two parts: one in which the beam is modeled using two sublaminates (top and bottom) and the other, where the beam is formed by just the top sublaminate. The purpose of theses tests is to prove the validity of the finite element developed in modeling the beam problem, including the torsion problem.

In the second numerical test, a specially orthotropic delaminated cantilever beam is put under torsion loading. Initially, the delamination is supposed to be in the middle plane of the beam, and the result in terms of the strain energy release rate is compared with the closed-form solution. The delamination is then placed between different layers along the thickness, and a study is performed to understand the


79









80

effect of the delamination position on the strain energy release rate.

The last numerical test consists in the simulation of a quasi-static impact test in specially orthotropic and quasiisotropic simply supported beams. Some experimental observations for this specific tests are available in our Composite Laboratory (Kwon, 1991), and the present numerical study will help interpret the experimental results. Again, the focus was on the effect of the delamination position on the values of the strain energy release rate.



4.2 Material Properties and Stiffness Coefficients

In all the numerical tests performed in the present study, the beam is supposed to be made up of material with the following properties:

Longitudinal Elastic Modulus (El): 14.00 GPa Transverse Elastic Modulus (E2): 1.00 GPa Shear Modulus (G12): 0.53 GPa Poisson's Ratio (V12): 0.30 Poisson's Ratio (V23): 0.55

The above values are typical values for high performance graphite/epoxy unidirectional composites.

Table 4.1 lists the laminate stiffness coefficients for a specially orthotropic beam with thickness equal to 2 mm, and the reference X-Y plane considered to be in the middle









81

(MIDDLE), on the top (BOTTOM) or on the bottom (TOP) surfaces of the beam.



4.3 Cantilever Beam

The dimensions of the beam are 0.2 m x 0.02 m x 0.002 m (length x width x thickness), and the beam is subjected to a set of six different end loading conditions, as shown in Fig.

4.1 and Fig. 4.3.



4.3.1 Model with two sublaminates

Initially, the beam subjected to the forces shown in Fig. 4.1 is modeled by using two connected sublaminates (top and bottom), as shown in Fig. 4.2. Table 4.2 presents the numerical and beam theory solutions, in terms of the displacements at the tip of the beam for each of the loading conditions.

The expressions for beam theory results can be found in Timoshenko (1970), Reissman and Pawlik (1980) and Sankar (1991), or can be derived from the beam equations presented in Sankar (1991). In the following expressions the superscripts t and b stands for top and bottom, respectively.


4.3.1.1 Beam under Force 1 (F11

In this case, we have just the U displacement. All other displacements are zero.









82





Table 4.1 Stiffness coefficients




STIFFNESS POSITION OF THE X AXIS COEFFICIENTS
TOP BOTTOM MIDDLE A11 0.2818173E8 0.2818173E8 0.2818173E8


A16 0.0 0.0 0.0

A6 0.1066621E7 0.1066621E7 0.1066621E7

B11 0.2818229E5 -0.2818229E5 0.0 B16 0.0 0.0 0.0 B6 0.1066643E4 -0.1066643E4 0.0


D11 0.3757714E2 0.3757714E2 0.9394285E1

D16 0.0 0.0 0.0

D6 0.1422219E1 0.1422219E1 0.3555547


K2 A55 0.8888511E6 0.8888511E6 0.8888511E6








83




F
z

X FF











(a)



F F
36












(b)




Figure 4.1 Loading conditions at middle plane
(a)Perspective view (b)Front view








84





L
z





cantilever beam






z





top offset elements











bottom offset elements







Figure 4.2 Model with 42 elements






Table 4.2 Displacements for model with 42 elements


FORCES DISPLACEMENTS AT THE TIP OF THE BEAM U (m) V (m) W (m) F (rad) 1 (rad) a (rad) 0 (rad)


1 FEM 3.55E-7 0.0 0.0 0.0 0.0 0.0 0.0


BEAM 3.55E-7 0.0 0.0 0.0 0.0 0.0 0.0 2 FEM 0.0 1.51E-4 0.0 1.06E-3 0.0 0.0 0.0


BEAM 0.0 1.64E-4 0.0 1.06E-3 0.0 0.0 0.0 3 FEM 0.0 0.0 0.0142 0.0 0.1064 0.0 0.0


BEAM 0.0 0.0 0.0142 0.0 0.1064 0.0 0.0 4 FEM 0.0 -1.06E-3 0.0 0.0106 0.0 0.0 0.0


BEAM 0.0 -1.07E-3 0.0 0.0107 0.0 0.0 0.0 5 FEM 0.0 0.0 -0.1064 0.0 1.0645 0.0 0.0


BEAM 0.0 0.0 -0.1064 0.0 1.0645 0.0 0.0


6 FEM 0.0 0.0 0.0 0.0 0.0 -34.7343 6.6082


BEAM 0.0 0.0 0.0 0.0 0.0 -34.7343 7.0312


0o
(.'1









86

L (4.1) U-FL
1 b (A t+A b11)

4.3.1.2 Beam under Force 2 (F21

Under Force 2 the beam has the displacements V and F. The remaining five displacements are zero.


( 4L3 2L ) (4.2) Sb (A tl+A b11) b (A 5ss+A b


6L2 (4.3) F--F,
b3 (A t 1+A b11)



4.3.1.3 Beam under Force 3 (F31

From this case results the displacements W and 4.

L3 2L (4.4)
3b(D t11+Db11) b(A tss+A b55)


1 L2 (4.5)
-F
2 b(D t+D 1b)



4.3.1.4 Beam under Force 4 (F")

The displacements under Force 4 are V and F. The remaining displacements are zero.

6L2 (4.6) V--F43 (A
'b~ (A ezz+A zz,)









87

312L (4.7)
b3 (A tl+A b11)



4.3.1.5 Beam under Force 5 (F,1

Under this loading condition the displacements are W and 0, and all others are zero.

1 L(4.8) W--Fs2 b(D tll+Dbn)


(4.9)

-F5 L
b (D tz+D bz)




4.3.1.6 Beam under Force 6 (FJ)

The displacements under Force 6 are 0 and a. The displacement Ois given by the expression (Sankar, 1991b)


-FL 1+ D Dt+Db66 (4.10)
6 4b(D t66+Db66) D 55+Db55 where

(4.11)
D tss+D ss55 k (A t55+A bss)









88

and


k2_-5 (4.12)
6


while the displacement a is defined as 12 (4.13)
a-Cxcosh (x) +C2sinh (x) -C3 12 (4.13) b312 (D ezz+Db 11)


where


.12 D55 (4.14)
b3 2 (D t11+Db1) D*66


(4.15)
C2--C1l sinh(x) and


D***
C3
D * *66


In Eqs. 4.13-4.16,


XL 48 (D 66+Db 66) (D tss55+D b55) (4.17)
b2 (D tl+D b) D 66



L2. 12D**66 (4.18) b2 (D te1+Db11)









89


(A 5+1 (D t6+D b66) (4.19)
1 66



D**"- 12 b2K2 (A t5+A b55) 1 (D66 +D 66) (4.20) 6612 4




4.3.2 Model with top sublaminate

In a second test, the beam is modeled using top offset elements only (Fig. 4.3). The beam is subjected to the set of loading conditions displayed in Fig. 4.3. Table 4.3 shows the numerical and closed solutions for this test. The sources for the beam theory solutions presented in Table 4.3 are the same

used before: Timoshenko (1970), Reismann and Pawlik (1980) and Sankar (1991b). Some expressions are derived from the beam equations.



4.3.2.1 Beam under Force 1 (F1,

The tip displacements under F1 are U, W and 0. The remaining displacements are zero.

L
U - F At (4.21)



(h) L2
2 _ (4.22) b( 11)
4








90




z






Z




delaminatlon


(a)




F4








2 F
(b)



Figure 4.3 Loading conditions at the bottom plane
(a) Perspective view (b) Front view




Full Text
xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID EJ66KK3NG_Y4VXL3 INGEST_TIME 2015-04-13T19:13:29Z PACKAGE AA00029935_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES



PAGE 1

FINITE ELEMENTS FOR FREE EDGE DELAMINATIONS AND DELAMINATED ANISOTROPIC BEAMS BY MARCO ANTONIO SANTOS PINHEIRO A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1991

PAGE 2

In memory of my father-in-law Fernando Pinto still and forever a beacon to my family

PAGE 3

ACKNOWLEDGMENTS I would like to express my sincere gratitude to my advisor Dr. Bhavani V. Sankar for his continuous support, encouragement and constructive criticism during all my doctoral studies. His competence and understanding of his students always impressed me and are examples to be followed in my future life. Gratitude is extended to Professors U.H. Kurzweg, L.E. Malvern, M. Rao and C.T. Sun, members of my doctoral supervisory committee, for their assistance during my research. My appreciation goes to Dr. C.E. Taylor for his participation in the initial part of my work. Special thanks are owed to the Brazilian Army for giving me the opportunity to come to the United States of America to pursue doctoral studies. I would like to thank my friends at the Department of Aerospace Engineering, Mechanics and Engineering Science, whose friendship have been motivation for the completion of this work. Finally, I am grateful to my wife, Fernanda, my daughter Simone, my son Luis Marco, my parents Jose and Nely, and my mother-in-law Geny. They have been always a source of encouragement, support and love. iii

PAGE 4

TABLE OF CONTENTS Page ACKNOWLEDGMENTS iii LIST OF TABLES vi LIST OF FIGURES vii ABSTRACT ix CHAPTERS 1 INTRODUCTION 1 1.1 Motivation and Objectives 1 1.2 Literature Survey 4 1.3 Scope and Outline of this work 8 2 AN OFFSET BEAM FINITE ELEMENT FOR MODELING FREE EDGE DELAMINATIONS IN COMPOSITES 10 2 . 1 Introduction lo 2 . 2 Laminated Beam 16 2 . 3 Fabrication Process 16 2.4 Free Edge Problem 22 2.5 Finite Element Solution 26 2.6 Finite Element Models 28 2.7 Strain Energy Release Rate 30 2.8 Stresses Ahead of the Crack 39 2.9 Numerical Results 43 2.9.1 Strain energy release rate 44 2.9.2 Stresses ahead of the crack 45 2.9.3 Conclusions 49 3 MODELING DELAMINATION IN COMPOSITE BEAMS AND APPLICATION TO TORSION PROBLEM 56 3.1 Introduction 55 3 . 2 Beam Eguations * 53 3 . 3 New Set of Force and Moment Resultants 64 3.4 Equilibrium Equations 65 3 . 5 Finite Element Solutions .'.'.*.*.'.* 69 3 . 6 Strain Energy Release Rate .*.'.'.*.*.*.*.* 74 3 . 7 Rigid and Gap Elements .'.'.*.*.*.*. 75 iv

PAGE 5

4 NUMERICAL RESULTS AND DISCUSSION 79 4.1 Numerical Results 79 4.2 Material Properties and Stiffness Coefficients ... 80 4.3 Cantilever Beam 81 4.3.1 Model with two sublaminates 81 4.3.1.1 Beam under Force 1 (F^) 81 4.3.1.2 Beam under Force 2 (Fj) 86 4.3.1.3 Beam under Force 3 (Fj) 86 4.3.1.4 Beam under Force 4 (F^) 86 4.3.1.5 Beam under Force 5 (Fj) 87 4.3.1.6 Beam under Force 6 (F^) 87 4.3.2 Model with top sublaminate 89 4.3.2.1 Beam under Force 1 (F^) 89 4.3.2.2 Beam under Force 2 (Fj) 93 4.3.2.3 Beam under Force 3 (Fj) 95 4.3.2.4 Beam under Force 4 (F^) 95 4.3.2.5 Beam under Force 5 (F5) 97 4.3.2.6 Beam under Force 6 (F^) 98 4.3.3 Comments on the finite element solution .... 98 4.4 Delaminated Cantilever Beam under Torsion 102 4.4.1 Delamination at the middle-plane of the beam 105 4.4.2 Delamination position varying along the thickness 108 4.5 Simply Supported Beam under Transverse Force 112 5 SUMMARY AND FUTURE RESEARCH 116 5.1 Summary 116 5.2 Future Research 118 APPENDIX 119 REFERENCES 128 BIOGRAPHICAL SKETCH 130 V

PAGE 6

LIST OF TABLES Table Page 2 . 1 Comparison of energy release rate obtained from different methods 46 2.2 Relative difference (%) between the results from present method and other methods 47 4.1 Stiffness coefficients 82 4.2 Displacements for model with 42 elements 85 4.3 Displacements for model with 21 elements 92 4.4 Strain energy release rate (G) for cantilever beam. 110 4.5 Strain energy release rate (G) for simply supported beam 114 vi

PAGE 7

LIST OF FIGURES Figure Page 1.2 Types of delamination 3 (a) Free edge (b) Transverse delamination 2.1 Free edge delamination 11 2.2 Sublaminates 13 2 . 3 Finite element model for strain energy release rate 15 2.4 Finite element model for peel stresses 15 2.5 Values of h for top sublaminate 19 2.6 Force and moment resultants (top sublaminate) 23 2.7 Forces and degrees of freedom 29 (a) Nodal and average forces (b) Degrees of freedom 2.8 Zero-volume path for J-integral 32 2.9 Crack-closure method 37 2.10 Model for spring constant calculation 40 2.11 Stresses ahead of the crack tip for specimen A ... 50 2.12 Stresses ahead of the crack tip for specimen B ... 51 2.13 Stresses ahead of the crack tip for specimen C ... 52 2.14 Stresses ahead of the crack tip for specimen D ... 53 2.15 Influence of temperature variation 54 2.16 Influence of moisture level variation 55 vii

PAGE 8

3 . 1 Laminated composite beam 59 3.2 Force and moment resultants 63 3.3 Rigid elements: 77 (a) Beam under torsion (b) Interference on left side (front view) (c) Interference on right side (front view) 4.1 Loading conditions at middle plane 83 (a) Perspective view (b) Front view 4.2 Model with 42 elements 84 4 . 3 Loading conditions at the bottom plane 90 (a) Perspective view (b) Front view 4.4 Model with 21 top elements 91 4.5 Beam under loading F^ 96 4.6 Angle of twist along the length (L/h=100) 100 4.7 Angle of twist along the length (L/h= 10) 101 4.8 Delaminated beam 103 (a) Under torsion (b) Separated by parts 4.9 Finite element model 104 4.10 Delamination close to the top plane 109 4.11 Strain energy release rate for cantilever beam 111 4.12 Simply supported beam 113 4.13 Strain energy release rate for simply supported beam 115 A Stiffness matrix K 120 viii

PAGE 9

Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy FINITE ELEMENTS FOR FREE EDGE DELAMINATIONS AND DELAMINATED ANISOTROPIC BEAMS By Marco Antonio Santos Pinheiro December, 1991 Chairman: Dr. Bhavani V. Sankar Major Department: Aerospace Engineering, Mechanics, and Engineering Science The major objective of this study is to develop a technique to analyze delaminations in composite beams by using the finite element method. A method is proposed to obtain the strain energy release rate and the stresses ahead of the crack for the specific case of free edge delaminations. A shear deformable beam finite element with nodes offset to either top or bottom side has been developed. The thermal effects originating from the curing process are included in the formulation of the problem. Three different expressions to calculate the strain energy release rate are derived, one in terms of the force and moment resultants in the elements surrounding the crack-tip, and the two others in terms of the forces transmitted by a rigid element placed at the tip of the crack. The stresses ahead of the crack are obtained by using discrete spring ix

PAGE 10

elements in the uncracked portion of the beam. The results calculated are in good agreement with results obtained from other studies of the same problem. The advantage of using the present method is that it requires less computational effort, and is capable of modeling any type of geometry, loading, and boundary conditions. An efficient beam finite element has been developed based on a generalized beam theory for anisotropic beams, in which the equilibrium equations are assumed to be satisfied in an average sense over the width of the beam. This introduces a new set of force and moment resultants for the beam. The displacements in the beam are expressed in terms of three deflections, three rotations and one warping term. The effectiveness of the method is illustrated by solving problems where the beam theory solutions are known. The strain energy release rate is determined for the cases of delaminated cantilever beam under torsion and simply supported delaminated beam under transverse loading. Studies are performed to understand the influence of the position of the delamination on the strain energy release rate. X

PAGE 11

CHAPTER 1 INTRODUCTION 1.1 Motivation and Objectives Structural materials can bear defects proceeding from the fabrication processes or accumulate damage during their service life. Frequently, such defects and damage are not visually observable and have to be detected by using nondestructive evaluation techniques (NDE) . In a variety of applications, such as aerospace and aircraft structures, offshore platforms, bridges and tall buildings, the replacement of the imperfect parts is not a cost-effective solution for the problem. The costs involved in keeping the equipment out of order, the difficulty of access, and other reasons suggest that there is a need for research and development of methods to predict the occurrence and extent of the damages. Composite materials and structures, for example, are highly susceptible to damage from low-energy impacts. Incidents such as tool-drops, improper handling or hailstorms (aircraft structures) can cause internal damage and significant reduction in compressive strength of the material. Under subsequent transverse loadings, the damage will grow 1

PAGE 12

2 causing a redistribution of the stresses in the plies of a laminate and consequently reduction in the material stiffness and fatigue life. However, depending on the subsequent loading conditions and the damage characteristics, the structure often can still be used despite the existence of damage. Failure in a composite material often initiates in the form of a matrix crack or a delamination. The former refers to intralaminar failure or transverse crack, while the latter refers to interlaminar failure. Both types of failures are illustrated in Fig. l.l. Delamination, or separation of the layers, is a mode of failure unique to laminated materials, and it is the predominant failure mode in continuous fiber composites. Based on the location and direction of the growth, we can discern two types of delamination: the free edge delamination and the transverse delamination. In the edge delamination the growth initiates at the load-free edges of the laminate and propagates normal to the load direction. Under transverse delamination, also known as crack tip or local delamination, the growth starts from a matrix crack and propagates normal to the transverse crack from which it originates. The Fig. 1.2 shows both types of delamination. Since delamination is a very commonly observed failure mode in composite materials, its study is of great importance in the field of failure prediction for composite structures.

PAGE 13

Figure 1.1 Failures in composite materials (a) Matrix crack (b) Delamination delamination matrix crack (a) delamination (b) Figure 1.2 Types of delamination (a) Free edge (b) Transverse delamination

PAGE 14

4 1.2 Literature Survey In the past few years, several researchers have studied delamination occurrence and growth. Usually they predict delamination growth from the strain energy release rate (G) , and try to develop analytical and numerical methods capable of modeling this phenomenon. A brief description of some of these studies is important for accurate comprehension of the present work. Chan and Ochoa (1987) developed a two-dimensional finite element model to calculate interlaminar stresses and strain energy release rate in composites subjected to uniaxial tension, bending and torsion. The element is an isoparametric element with eight nodes and three degrees of freedom per node. The basic idea in their work is to represent the effect of longitudinal or twisting curvature as a "load effect." The study concentrates on establishing relationships of Gj , G,, and Gjjj as functions of crack length and stacking sequence. Pagano (1974) studied the case of fibrous composite laminates under a uniform axial extension and developed an approximate method to define the distribution of the interlaminar normal stress along the central plane of a symmetric laminate. A simple expression for the interlaminar normal stress was derived and the results obtained by using this expression were compared with those obtained by using a finite difference approximation of the elasticity equations. The two methods agreed well, and both results show that

PAGE 15

5 significant stresses exist only in narrow regions next to the free edges. The length of this region is approximately equal to the laminate thickness. Whitney (1986) proposed an approximation to the interlaminar normal stress distribution ahead of the crack and the influence of specimen geometry on Mode I strain energy release rate for edge delamination specimens. The specimens were subjected to tension and compression loadings and were of the classes [0^/90^]^ or [B/S^/e/so^]^. Whitney used a higher order plate theory, which considers transverse shear deformation and thickness-stretching. The effect of thermal residual stresses due to lamination fabrication was considered in the formulation of the problem. The numerical results suggested significant influence of the residual thermal stresses on both the stress distribution ahead of the crack and the Mode I strain energy release rate. In general, the thermal residual stresses increase the stresses at crack tip and the Mode I strain energy release rate for specimens under tension, and has the inverse effect for specimens under compression. Again the results confirmed that the effect of normal stress ahead of the crack dissipates within a distance of approximately twice the thickness from the crack front. Armanios and Badir (1989) used a shear deformation theory to calculate the interlaminar stresses and energy release rate associated with both Mode I and mixed-mode edge delaminations. The formulation used by Armanio and Badir does not consider

PAGE 16

the thickness strain effect. The analysis includes the effects caused by the residual thermal stresses and moisture stresses. The classes of laminates considered are the same ones used by Whitney (1986) . Results suggest that the residual thermal stresses tend to increase the interlaminar stresses and total energy release rate for both modes, while the moisture content tends to alleviate the thermal effects. The value of moisture content necessary to alleviate totally the interlaminar stresses and the strain energy release rate ranged between 0.70% and 0.80% in the examples considered, and it was found that this value has a weak dependency on both the stacking sequence and the mode considered. Usman and Gandhi (1989) studied the effect of moisture and temperature on the vibrational characteristics of composite materials. The nonlinear vibration of moderatelythick laminated composite plates under hygro-thermal conditions are considered by incorporating the effects of transverse shear and rotary inertia. The hygro-thermal changes along the thickness of the composite were explicitly incorporated in the form of the deformation field. The response of the composite under hygro-thermal environments was modeled by creating one model formed by two parts, one anisotropic solid (composite) and one ideal fluid (moisture) . The equations presented in this work are significantly different from the corresponding equations in classical laminated composite plate theory.

PAGE 17

Hu et al. (1989) developed a novel beam element to study the crack propagation in beam types of specimens. Finite elements are used to model the specimens, and the beam element has its nodes offset to one side, either bottom or top side. This approach allows modeling crack propagation without renumbering the nodes as the crack propagates, and is convenient in modeling the partial contact and slip at a delamination interface. The bottom and top un-delaminated portions of the beam were connected either by foundation springs or rigid elements. For some specimens the strain energy release rate was computed by using a zero-volume path integral and a crack closure method, and the results were found to be in good agreement. In the present study the beam theory is extended to analyze delaminated composite beams. Sankar (1991 b) developed a beam theory for laminated composite beams derived from the shear deformable plate theory. The equilibrium equations are assumed to be satisfied in an average sense over the width of the beam, which introduces a new set of force and moment resultants for the beam. The displacements in the beam are expressed in terms of three deflections, three rotations, and one warping term. The closed solutions for the cases of a cantilever beam subjected to end loads and a specially orthotropic laminate beam under torsional loading are presented. The solution considers the shear deformation and the edge effect. The results in terms of angle of twist agree well with available solutions.

PAGE 18

8 1.3 Scope and Outline of this Work Numerical and analytical methods have been used to express the complex stress state near free edges and to calculate the energy release rate during crack propagation. Among several numerical methods, the finite element method was selected to be used in the present work. Although requiring intensive computational effort, this method is very popular and provides accurate solutions. Also, this method is very convenient in modeling structural elements of variable cross section. In the present work, the typical element used to model the structure is the beam element. This element not only is very common in mechanical structures but also is the most appropriate to model some standard test specimens. The scope of this work can be described as to develop a method to analyze delaminations in anisotropic composite beams, using the finite element method to model the problem and selecting the strain energy release rate as a delamination parameter. The loading is considered independent of time. This dissertation consists of five chapters. Chapter 2 is focused on free edge delamination. A method is proposed for modeling this type of delamination using offset beam finite elements. The effect of temperature variation during the curing process (residual stresses) is considered in the problem. Chapters 3 and 4 are devoted to presenting a new and

PAGE 19

9 efficient method to analyze composite laminated beams, using a more complete and accurate beam finite element. The capability of this new finite element in modeling the problem was verified by solving the problem of a cantilever beam subjected to different end loading conditions, including torsion. Applications to delaminated beams under torsion loading and the simulation of a quasi-static impact test on a delaminated beam are also presented in these chapters. Finally, chapter 5 summarizes the principal conclusions of the present study, and presents suggestions for future studies in this field.

PAGE 20

CHAPTER 2 AN OFFSET BEAM FINITE ELEMENT FOR MODELING FREE EDGE DELAMINATIONS IN COMPOSITES 2 . 1 Introduction The separation of plies or delamination is a predominant failure mode in laminated composite structures. It occurs when interlaminar stresses exceed the interlaminar strength of the material, or as a result of the fabrication process, environmental degradation and improper handling of the material. Common examples of such failure regions during the normal use of composite structures are those near cut-outs, holes and free edges. In these regions, the abrupt changes in the geometry or material properties tend to intensify the stresses when the laminate is subjected to some external loads. Specifically, the occurrence of delamination in free edges has been receiving increasing attention of many investigators, in their efforts to understand and to prevent delaminations in composite structural designs. The reason is the presence of high interlaminar stresses, the "peeling" stress, in the neighborhood of a free boundary. The "free edge problem" usually is modeled by a plate subjected to a uniform axial strain ( e^ ) , as shown in Fig. 2.1. 10

PAGE 21

11 (applied strain) Figure 2.1 Free edge delamination

PAGE 22

12 During the usual fabrication process, the composite material is subjected to temperatures of several hundred degrees Fahrenheit, and then it is cooled down to the room temperature. The composite is assumed to be stress-free at the curing temperature. The change in temperature causes unavoidable stresses to the composite material layers when cooled to room temperature. Usually these stresses are called residual stresses or thermal stresses, and they are considered in the formulation of the present problem. A simple procedure is developed in this study to obtain an approximate solution for the peel stresses ahead of the crack tip and to calculate the strain energy release rate, for the case of free edge delaminations including thermal effects. The problem depicted in Figure 2.1 is modeled using a shear deformable beam finite element with offset nodes (Sankar, 1989) . The three-dimensional problem is modeled as a plane problem in the x-z plane assuming that the stresses and strains do not vary along the y-direction. This will be true for most of the length of the plate in the y-direction. The delamination plane divides the beam into two parts, the top sublaminate and the bottom sublaminate (Figure 2.2), with nodes offset to the bottom and top, respectively. The delamination plane is assumed to be the weakest one, therefore the crack propagates along this plane, without crack branchings .

PAGE 23

bottom sublaminate Figure 2.2 Sublaminates

PAGE 24

14 When calculating the energy release rate, the uncracked region of the beam is connected by common nodes, except the tip of the crack, where a rigid element is used to connect the top and bottom nodes (Fig. 2.3). This rigid element not only ensures the continuity of displacements and rotations at that section, but also brings in its formulation the generalized forces that are transmitted between the sublaminates . These forces are used in the calculation of the energy release rate. In evaluating the peel stresses ahead of the crack, the uncracked region of the beam is connected by spring elements (Fig. 2.4) , stiffness of which has to be judiciously selected. The strain energy release rate has been successfully adopted by researchers in characterizing delamination. In this study its value is obtained from three different methods (Sankar and Pinheiro, 1990) . The first one uses in its formulation the forces transmitted by the crack-tip rigid element. In the second method, the result is obtained by using a method similar to the crack-closure method used in twodimensional fracture problems. Finally, composites are very brittle materials, so that, the plastic zone in the vicinity of the crack tip is negligible, and hence the linear elastic fracture mechanics can be applied in which case the J-integral is equivalent to the energy release rate. Therefore, the third method consists in obtaining the J-integral from the beam model. The force and moment resultants in the region surrounding the crack-tip are used in this calculation.

PAGE 25

15 cracked region uncracked region Figure 2.3 Finite element model for strain energy release rate Figure 2.4 Finite element model for peel stresses

PAGE 26

16 The results from the three different methods are compared and found to be in good agreement. 2 . 2 Laminated Beam In this section a homogeneous orthotropic laminate subjected to a uniform initial strain = e^ applied over a length 2L is analyzed. The thickness of the laminate is 2h and the width is 2b. The thickness 2h is much smaller than 2L and 2L much smaller than 2b (Fig. 2.1). Due to symmetry conditions, just one half of the laminate is considered here. The problem is divided into two different phases. In the first phase no initial strain is applied to the specimen, which is subjected only to thermal stresses due to the curing process. In the second phase the cured specimen containing the edge delamination is subjected to the axial strain = in addition to the residual stresses. The equations derived in the next sections refer to a laminated beam situated above the delamination plane. The expressions for the laminated beam situated below the delamination plane differ only in the integration limits ( -h to 0 , instead of 0 to h ) . 2.3 Fabrication Process During the composite fabrication process there are heating and cooling phases. In the cooling phase, the temperature drops from the stress-free temperature to the room

PAGE 27

17 temperature. In this case, each layer influences the expansion or contraction of the other because of the difference in their coefficient of thermal expansion. A similar situation exists when the composite absorbs moisture during service. As a result of the temperature variation, each layer has residual stresses after the laminate fabrication. The effects of the residual stresses on the laminate as a whole are calculated and expressed in terms of thermal forces and moments. In the next phase, these thermal forces and moments are added to the external forces and moments. During the cooling process the in-plane and transverse displacements are assumed to be of the form: u{x,z)Uq{x) + j^ix) (2.1) v{y,z)v^iy) + z^ yiy) w(x, z) Wq ix) where Ug and Vp are the axial displacements of points on the x-y plane in the x and y directions respectively, and and are the rotations about the y and x axes respectively. It is assumed that the modification in the thickness of the laminate is negligible when compared with the modifications in the other two dimensions. Therefore, w is assumed to be constant through the thickness of the laminate (no thickness strain) .

PAGE 28

18 From Equations 2.1 the strains for the cooling process are derived as (2.2) A state of plane stress normal to the z direction is assumed during the curing process. Therefore, the thermal stresses developed within the k*^ layer (Agarwal and Broutman, 1980) are '^11^12^16 ^12^22^26 A6^26^66 XX » X e -C yy ' y -c xyj (2.3) Where Q.. are the stiffness coefficients and can be found in Whitney (1987), and Cx<^x^T (2.4) In the above equations a is the thermal expansion coefficient and AT is the change in temperature. By definition, the force (N) and moment (M) acting on the top laminate are .,M^) -Jo ^[1, z] dz 0 h {Ny,My) -fay[l,Z] dZ (2.5)

PAGE 29

19 In the present case, there are no net resultants acting on the laminate. Therefore, considering the definition of force and moment resultants, we have for the whole laminate dz (2.6) and MA i-s/C;), zdz (2.7) where n is the total number of layers on the top sublaminate and (h,^ h,^.^) is equal to the thickness of the k**' layer (Fig. 2.5). layer number T n n-1 plane of delamination Figure 2.5 Values of h for top sublaminate

PAGE 30

20 The stiffness coefficients of the top sublaminate are defined as (2.8) i/D ~ 1/2,. ..,6 In the above equation, a shear correction factor equal to 5/6 (Whitney, 1987) is introduced in the calculation of the terms ajj,b.. and d,.,., wherever the index i is greater then 3. Substituting Eq. 2.3 in Eq. 2.6 and 2.7, and using the definition of the stiffness coefficients in Eq. 2.8, we obtain the following expressions for the force and moment resultants xf (o°) •12 ^la V^12 C?i2j /_0 (2.9) ^12 ^12 ^22 ^22 ^22/ y (2.10) where N^^ and M^^ are the thermal force and moment in the x direction and so are N^^ and M^^ in the y direction.

PAGE 31

21 The thermal forces and moments are defined as [nJ,mJ]fiO^^C^^ Q^^Cy^QrsC^cy) n,z]dz (2.11) [<,m/]/(Cx2C, + 022Cy + 02sC^) dz (2.12) 0 Defining the matrices and [L][M]' ^11 ^11^ \,'^22 "^22^ Equations 2.9 and 2.10 take the form of and /0\ ) [L] + [M] ) m From Eq. 2.17 we obtain [N] (2.13) (2.14) (2.15) (2.16) (2.17) (2. 18) and, after substituting from Eq. 2.18 in Eq. 2.16, we obtain (2.19) T [p] r + [M] m -1 / 7^

PAGE 32

22 where [P][L] m m-"IM] Eq. 2.19 defines the thermal force and moment, and respectively, as [M] m-^l^'[P]{\A (2.20) From now on, the thermal force and moment resultants defined in Eq. 2.20 are considered as external forces and moments acting on the composite specimen (Agarwal and Broutman, 1980) . The finite element method is used to solve Eq. 2.20 and to obtain the solution for e°^(x) and K^ix) . For bottom element the derivation is the same, except for the limits of integration in Eq. 2.5, 2.8, 2.11 and 2.12. 2.4 Free Edge Problem In defining the in-plane and transverse displacements for the free edge problem (second phase), Eq. 2.1 remains the same, except for displacement in the y direction, which takes the form of v(y,z) e^y (2.21)

PAGE 33

23 Where e^, represents the applied strain in the y direction. Consequently, the strain e^^ in Eq. 2.2 is replaced by e^y (2.22) Because the dimension 2b in Figure 2 . 1 is considered much greater than the dimension 2L, it is assumed that a state of plane strain parallel to the x-z plane exists in the second phase of the free edge problem. The stress-strain relations for the k'^*' layer in this case are 0 0 0 (e, y XX 0 y xz, (2.23) The in-plane force and moment acting on the beam are shown in Fig. 2.6 Figure 2.6 Force and moment resultants (top sublaminate)

PAGE 34

24 P, M and V are defined as h [P,M,V][ [a^,ZO^,X^^] dz (2.24) 0 or, using the definitions given in Eq. 2.8 and 2.23, as (2.25) ^11 0 ^ XX \m 0 + 1° 0 y xz, [o) Where and are the force and moment induced by the applied strain e^. They are defined as h [P^^^^]1^(1,2) dz [ai2,23i2leo (2.26) 0 The thermal force and moment resultants defined in Eq. 2.20 are added to the problem in the form of external forces, and hence the problem takes the final form of (2.27) (p) M + 1 0 ; 0^ XX 0 1° 0 [o )

PAGE 35

or \m 311 ^11 ^11 ^11 0 0a 0 0 55/ ( [o ) [o j ) (2.28) The last equation can be expressed in a more concise form as F S e R (2.29) where S is the stiffness matrix, components of which are [5]b b. 11 " 11 c?ii 0 0 0a 557 (2.30) and R is the vector that represents the forces and moments induced by the applied strain and by the thermal effects {Ry pH [o ) lo j (2.31) The vector F refers to the external forces acting on the beam and its components are (F) M (2.32)

PAGE 36

26 while e is the vector of deformations in the form of (2.33) Finally, the Equation 2.27 can also be expressed as e c (F+R) (2.34) where c = s'^ and is the compliance matrix. 2.5 Finite Element Solution The Finite Element Method is used to solve Eq. 2.29, and a beam finite element with offset nodes is selected to model the problem. The derivation of this specific beam finite element was done by Sankar (1989). The finite element has two nodes and three degrees of freedom (\x,(f/ and w) at each node. The nodal forces and moment for the i*** node are denoted by p,., m,. and v,. . In the present study a slightly different method is used to derive the stiffness matrix. The average force and moment resultants (P, M and V) at the center of the element are given by (P\ ((P2 Pi) / 2\ M m^) / 2 W) [iv^ V,) / 2^ (2.35)

PAGE 37

The deformations at the center of the element are expressed in terras of the nodal displacements as (e)zx) (^(T^ + T,) /2+ IL] (2.36) where L is the element length. The beam finite eleraent forces and degrees of freedom are shown in Fig. 2. 7. By substituting from Eq. 2.33-2.34 into Eq. 2.29, the relations between nodal forces and displacements are obtained in the form of f k g R (2.37) where the vectors of nodal forces (f) and displacements (q) have the following components (g) ^[Ui,l|fi, W^, U2,^2'^2^ (2.38) The vector R and its components are defined by Eq. 2.20, 2.26 and 2.31. The elements of the stiffness matrix k are

PAGE 38

28 '11 '11 L 0 '11 '11 L 4 ^55 '11 *55 "55 L 0 *ii L ^1 *ii ,.4. '11 *55 '11 "^55 2 ^55 '11 ni ) *55 *55 *55 '11 L 0 *55 ) *55 *55 (2.39) 2 . 6 The Finite Element Models Fig. 2.3 and 2.4 show the models used in modeling the free edge delaminations . They consist of two sublaminates above and below the delamination plane. The assumptions are that the delamination plane is the weakest plane in the material and that there is no frictional resistance against sliding of the sublaminates. The first assumption means that the crack will propagate along the delamination plane, without crack branching. In the uncracked region, the top and bottom portions of the beam are connected by common nodes, when calculating the strain energy release rate (Fig. 2.3) , and by spring elements, when predicting the stress ahead of the crack (Fig. 2.4).

PAGE 39

29 Z (a) Figure 2.7 Force and degrees of freedom (a) Nodal and average forces (b) Degrees of freedom

PAGE 40

30 2.7 Strain Energy Release Rate In the calculation of the strain energy release rate (G) , one rigid element is placed in the tip of the crack connecting both sublaminates (Fig. 2.3). The rigid element has nine degrees of freedom. Apart from the three displacements at each node, the three generalized forces, F^, transmitted by the rigid element are also considered as unknown degrees of freedom. The constraint equations corresponding to the rigid element are 'O3 O3 -J3^ 03 03 J3 (2.40) where F, and are the nodal forces at nodes 1 and 2 and and are the nodal displacements. The symbol O3 represents the 3x3 null matrix and I3 the 3x3 identity matrix. Considering a zero volume path surrounding the crack tip as shown in Fig. 2.8, the expression for the J-integral relative to two-dimensional states can be found in Hellan (1984) as /du iwdx^-Pi^ds) (2.41)

PAGE 41

31 where Pi rijOji i-1,2 For the present case, and considering the path r very thin (zero-volume), the Eq.2.41 takes the form of J fwdz fn^o^u^ds-fn^x^^w,,ds (2.42) r r r where W is the strain energy density and is defined as -| (Oxx^xx + Oyyfiyy ^z;cYzx) (2.43) The path r has four segments, one in each sublaminate behind and ahead of the crack tip (Fig. 2.8). The values of n^^ (unit vector normal to the path) and ds for each element are: Elements 1 and 4: n^^ =-1; ds =-dz Elements 2 and 3 : n^^ = 1 ; ds = dz Considering the Element 1 in Fig. 2. 8 we have: h h h J^-fwdz^ ja^u,,dz^jx^^w_^dz (2.44) 0 o 0 where the subscript stands for the element number. The first term of the right hand side of Eq.2.44 is defined by using Eq.2.43 and the stress-strain relations given in Eq.2.2 3 as ti h (2.45) 0 0

PAGE 42

32 "1,4 2,3 L » element » 1 i — 1 T " elemi»rif:'4:>:' 1 elements ahead elements behind of the crack tip of the crack tip Figure 2.8 Zero-volume path for J-integral

PAGE 43

33 From the definition of the stiffness coefficients and Eq. 2.30 and 2.3 3 we obtain the following relation h e^Sef {Ai^^'^Qs^yz.') dz (2.46) 0 Therefore, the Eq.2.45 can be modified and expressed as h jwdz-^e^Se^ + a^^e^G^ ^ ^a,,e,' (2.47) The second term of the right hand side of Eq.2.44 is modified in a similar way to what has been done for the first term. h h {^xxU.x<^z-f (^iifixx^ + ^iaG^eJ dz (2.48) 0 0 or h h f^^^,:c<^z-f (^iie«'-^^i2G«.eo + ^55Yzx'-^55Yz/) dz (2.49) 0 0 Finally, h I'^xxU.xdz-e^^Se^ + a^^e.e^-a^^yJ (2.50)

PAGE 44

34 The third term of Eq. 2.42 is h f-^zxK.^Z-^x^i (2.51) 0 where and ip^ are the shear force resultant (Eq. 2.24) and the rotation at the crack tip, respectively, for the element number 1. Substituting Eq. 2.48, 2.50 and 2.51 into Eq. 2.44 we obtain j^-le^-'Se^ + (2.52) where U^. is a energy term defined by -| 322^02-355 Yz;c' (2.53) similar expressions can be obtained for the other elements in the form of -\e^^Se^ (2.54) For linear elastic fracture, the J-integral is equivalent to the energy release rate (Hellan, 1984) . This is valid for fractures in brittle composite materials, where the plastic zone in the vicinity of the crack is negligible. Hence,

PAGE 45

35 G J + J2 + J^i + (2.55) The Eq.2.55 can be calculated in a more convenient form, where the forces F. (1=1,2,3,4) are calculated from the finite element solution, and so are the compliance matrices c^ and c^^ from the input data. The elements 2 and 3 are in the uncracked region, and so, they are under the same deformations. We thus obtain the relations, cP-cF (2.57) Substituting the last equation, and also the equilibrium relations (F^ + F^ = Fj + F3) , in Eq. 2.56 we obtain G-^ (F/c,F, + fJc^F, F/c.F^ F.^C^F,) (2.56) G-1 {F,^c, F/c^) (F.-FJ (2.59) or, using again the equilibrium relations. (2.60)

PAGE 46

36 Actually, (F, Fj) is the force transmitted by the rigid element, F^, between bottom and top crack tip nodes. Therefore, the strain energy release rate (G) can be expressed in a more convenient form as G i.F/(c, + c^) (2.61) where F^ are considered unknown degrees of freedom and part of the finite element solution of the problem (Eq. 2.40). One can also evaluate the value of the strain energy release rate using a method analogous to the virtual crack closure method in two-dimensional fracture problems ( Hellan, 1984) . In this method, the energy necessary to "close" the crack is equal to the energy release rate. In Fig. 2.9, as the crack propagates, the crack tip Nodes 1 and 2 have displacements and Ag, respectively, and the crack tip moves to Node 5. If the length of the elements behind the crack tip is sufficiently small, AL, then the crack opening, (Aq^-Aqj) , will be approximately equal to the difference in displacements of Nodes 3 and 4,(Aq3-Aq^). Therefore, the problem is analyzed considering the displacements of Nodes 1 and 2 as the tip displacements of cantilever beams of length AL under applied forces F^ on each beam. For very short cantilever beams (thickness/ length «1)

PAGE 47

Figure 2.9 Crack-closure method

PAGE 48

38 the tip deflections (q^jp) can be calculated using the expression g^ip L [c] F (2.62) where L and F are the beam length and the forces applied at the tip of the beam, respectively, and c is the compliance matrix. In the present case, Ag,--ALc,F, <2.63) or Ag^-Ag^ -(c, + c^)F^ (2.64) AL Substituting the last equation in Eq. 2.61 gives ^lil^r^(Ag,-Ag,) (2.65) In this equation, the displacements in the tip nodes and the forces in the rigid element are part of the finite element solution. In the present study, the strain energy release rate was calculated by three different methods. These methods are represented by Eq.2.56, 2.61 and 2.64, and the results obtained from the three equations are in very good agreement.

PAGE 49

39 2.8 Stresses Ahead of the Crack In evaluating the stresses ahead of the crack tip, the uncracked region was connected by discrete spring elements (Fig. 2.4). The stress ahead of the crack in the i^^ element was obtained by where V and L are the vertical force in the spring element and the length of the element, respectively. The vertical forces are obtained from the Finite Element solution of the problem, and the element lengths are part of the input data. The determination of the spring constant for the spring elements is of significant importance in appropriately modeling the problem. The spring constants are derived from the change in thickness of a beam subjected to transverse forces. The exact change in thickness is obtained from the two dimensional elasticity theory (Timoshenko, 1970) . In the first step, just the top part of the beam is considered and it is represented in Fig. 2.10. From the elasticity theory we have, for this case "T^ (y y^-c^y-|c3) (2.67)

PAGE 50

Figure 2.10 Model for spring constant calculation

PAGE 51

41 where E* is defined as ^"33(2.68) and A^-^12V21-V23V32-V3lVl3-2V2iV32Vl3 ^^'^^^ ^1 ^2 ^3 The vertical displacement (v) is obtained by integrating the deformation in the y direction. The result is where is a constant. From Eq. 2.70 the displacements at the positions y=-c and y=0 are determined as v(-c)^ Cl (2.71) 16 v(0)CI The relative displacement between these two particular positions is V,W-O-WO). -M^ (2.72)

PAGE 52

42 For a beam with unit area we have that P = q. In the present case, c is equal to half the thickness of the beam (c = h/2). Therefore Eq. 2.72 is rewritten as where k is the spring constant. By comparing between Eqs. 2.73 and 2.74 we have the value for the spring constant per unit length considering just the top part of the element. k IMl (2.75) 13h ^ ' Considering both parts, top and bottom, we have two springs placed in a series arrangement (Fig. 2.10). The spring constant equivalent (k^^) to this kind of arrangement is given by P 13h (2.73) For a spring model we have PkV^ (2.74) k. eg (2.76)

PAGE 53

43 where the subscripts t and b stand for top and bottom respectively. Substituting Eq. 2.80 in the last equation we obtain the final expression for the spring constant used in the present study 2.9 Numerical Results Some results from the investigation of strain energy release rate and the interlaminar stresses in a free edge specimen can be found in Whitney (1987) , (Pagano, 1989) and Raju (1986) . These researchers computed results for the same composite specimen. The present study uses the same specimen to make the comparison easier. The free-edge specimen width to thickness ratio (b/h) is 25 (Fig 2.11-2.16), which is typical for conventional edge delamination test (EDT) specimens. Considering a coordinate system where and X2 are parallel and transverse to the fibers respectively, and Xj is through the thickness, the ply properties have the following relations and values: 32E* (2.77) 13 (ht + Afc) E1/E2 = 14 E2/E3 = 1 G^j/Ej = 0.533 G23/E2 = 0.323 = 0.3 = 0.55 a = -9 E-07/°C Otg = "3 = 23 E-06/°C yffj = yffj = 5560 E-06/% 0: = 0

PAGE 54

' s ' s 44 where E. is the elastic moduli in the i*^ direction, G.. is the shear modulus in the x,.-Xj. plane and v.. are the Poisson's ratios. Finally, the thermal expansion coefficients and the moisture swelling coefficient in the x,. direction are denoted by a. and 0.^ respectively. These properties and relations among them are typical of current high performance graphite/ epoxy unidirectional composites. The specimens considered in the present study have the following lay-ups: Specimen A: [903/03]^ Specimen B: [-eo/eOgZ-eo/Oj], Specimen C: [-45/452/-45/02]g Specimen D: [O2/-6O/6O2/-6O], A reasonable stress-free temperature for standard graphite/ epoxy laminates is 138°C (280°F) (Whitney, 1987) . Assuming 21°C (70°F) as the room temperature, the difference in temperature (AT) that causes residual stresses in the plies is equal to -117°C (210°F) . The critical axial strains (eg) experimentally determined in edge delamination tests are 0.3% for angle-ply laminates, and 1% for the special case of a bidirectional class of laminates (Whitney, 1986) . Therefore, these values of strains are used in the numerical examples. 2.9.1 Strain energy release rate The normalized strain energy release rate, G^, is given by

PAGE 55

In the present study the energy release rate was calculated using Eqs. 2.56, 2.61 and 2.65. The results obtained from these three methods were ^n good agreement and they are refered to as "Present" in Table 2.1. Comparison of the results from the present study with available results is also present in Table 2.1. The relative differences between the results from the offset-node method and each of the other three approaches are shown in Table 2.2. From Tables 2.1 and 2.2 it may seen that the present results compare very well with Raju (1986) , whereas the comparison with the other two methods is not good in a few cases. It should be mentioned that of the three methods we compared with, Raju (1986) used a detailed finite element analysis. 2.9.2 Stress ahead of the crack The interlaminar stresses ahead of the crack tip are evaluated using spring elements in the uncracked region of the specimen. The spring constant for these elements is calculated by Eq. 2.77 and found to be equal to 300x lo' N/m .

PAGE 56

• (0 TJ O Eh O o o o VO o O o P Z O CO o VO in in "* 0) u o ID o rH rH VO IVO G w O n o o 00 • ' • • • • • • P K H CM H o o o o o C n (fl < 1 1 • Q) cm rH fN o o o o rH (0 > rat fEY o O o o VO o o z O o o o ease WHIT • H CO • rH • • o • o • o CO • o rH a) )H >1 tp • ener At • rH • o H • o • H • o • rH H 1 rH rH 1 c 1 ai str n n n n o o • • • • • • • ((H CO o o o o rH rH o 0 c 0 (0 -H u t/l u fM o pa O VO fNJ o Com ATE o VO 1 1 r\j O » in T H !2 (NJ in O M O VO o VO O ID ble o VO 1 1 IM o CTV (-45, (0 ^ Eh

PAGE 57

47 T3 (0 T3 0 X! P u 0) 0) B 43 +J P 0 C 0) (0 04 u B 0) 0 Li 0 (0 1 1 \ rH T3 3 O U) J3 0) +J 0) 0) 4J +j C (U c (0 0) (U 0) ^ -p X o o H 0) II 0) 0 0 c c 0) 0) V4 0) • 0) TJ <»-• H O -H T! J3 T3 -P (U Q) 0) > e > •H -H P >^ -P (0 0) (C •H 43 -H 0) -P (U ct: O (i: CM
PAGE 58

48 Fig. 2.11-2.15 show the results for peel stresses for different laminate configurations, with and without thermal effects. From the results we find that the effect of residual thermal stresses is to cause increase in the maximum peel stress at the free edge. The effect is significant for specimen A ( [9O3/O3] ) and less significant for specimen C ( [-45/452/-45/0213 ). In Fig. 2.11 the stress distribution is compared with Whitney (1987) for A= 0°C and A= -117 °C for the specimen B ( [-eo/eOjZ-eo/Ojlg ). The agreement is excellent. For the bidirectional laminates, specimen A, the influence of temperature can be better appreciated in Fig. 2. 15, where the temperature variation goes from zero to -120 °C in increments of -30°C. Besides thermal stresses, some studies on effect of moisture absorption on the free edge stresses were also performed. In fact, the analysis procedure is similar except that a AT has to be replaced by )S AH, where yff is the swelling coefficient and AH represents the variation in the moisture concentration. Fig. 2.16 shows the influence of the moisture level variation on specimen A. The moisture level was allowed to vary from 0% to 1.2% of the laminate weight, which reflects feasible conditions. This result agrees with those obtained by Armanios and Badir (1989). Unlike the thermal stresses, the residual moisture tends to decrease the interlaminar stress close to the crack front.

PAGE 59

49 2.9.3 Conclusions The residual thermal stress tends to increase the total energy release rate and the interlaminar stress at the tip of the crack, while the residual moisture alleviates these effects. The results obtained from the offset nodes method are very close to those obtained by Whitney (1987), Pagano ( 1989) , Raju(1986) and Armanios and Badir (1989) using different methods. Therefore, the present offset-node method is a good method to predict the energy release rate and stresses ahead of the crack tip in the free-edge problems. The present finite element method has the facility to model any type of geometry, loading and boundary conditions, and avoids renumbering the nodes and elements as the crack propagates. The offset node method can be used in studying dynamic delamination, and can be extended to problems of sublaminate buckling and delaminations in composite plates.

PAGE 60

50

PAGE 61

51

PAGE 62

52

PAGE 63

53 n

PAGE 64

54

PAGE 65

55

PAGE 66

CHAPTER 3 MODELING DELAMINATION IN COMPOSITE BEAMS AND APPLICATION TO TORSION PROBLEM 3 . 1 Introduction The application of composite materials is increasing day by day. Not only traditional areas like aerospace and automobile structures, but also several other areas, such as robotics, marine industries and medical devices and prostheses, have shown that the application of composite materials is major factor in the improvement of their mechanical structures. The advances in the manufacturing processes of composite materials, the development of new matrices and fibers, and the progress in methods for study of the mechanical behavior of composites have been made it possible for the composite materials to replace traditional materials in several mechanical structures. Usually composite materials are considered as a possible and adequate solution whenever the weight savings and the capability of "tailoring" the material in according to specific needs are more critical than the costs. As the use of composite materials in structural applications increases, there is more need for structural 56

PAGE 67

57 analysis. Unlike an isotropic beam, stretching a laminated composite beam can produce bending, shear and twisting in addition to the extensional deformation. The coupling among the modes of deformation is the major obstacle to precisely modeling the mechanical behavior of beam like composite structures . In the present chapter, a new and efficient method to analyze delaminations in laminated composite beams is presented. The equations of equilibrium were derived by Sankar (1991) from the classical shear deformable laminated plate theory (Whitney, 1987). The equilibrium equations are assumed to be satisfied in an average sense over the width of the beam. From this assumption results a new set of force and moment resultants, which amplify the possibilities for modeling beams under different loading conditions, including torsion. The equilibrium equations are derived from the Minimum Potential Energy Principle. In the present study a new offset beam finite element is developed for modeling the problem. Initially, the capability of the new finite element in modeling the problem is verified by solving the problem of a cantilever beam under different end loading conditions. The results are compared with those from beam theory solutions found in Timoshenko (1970), Reismann and Pawlik (1980), Whitney (1987) and Sankar (1991). Next, the strain energy release rate for a delaminated specially orthotropic cantilever beam under torsion is

PAGE 68

58 calculated. The delamination is considered to be at the middle-plane of the beam. The results for G from the finite element method are compared with analytical solutions of the problem. The effect of position of delamination on the strain energy release rate values is presented. Finally, the problems of a specially orthotropic and a quasi-isotropic simply supported delaminated beam under a transverse force is considered. The strain energy release rate values for various positions of the delamination are computed by using the finite elements developed in the present study. The results will be useful in explaining delamination propagation in composite beams due to quasi-static impact. 3 . 2 Beam Equations Considering the laminated beam shown in Fig. 3.1, the Shear Deformation Theory (Whitney, 1987) for laminated plates assumes the displacement field as uix.y, z) -Uo(x,y) +zilr^(x, y) v{x,y,z)'Vo ^/y) +zty(x,y) (3.1) w(x,y, z) -Woi^.y) where Ug, Vq and Wg are the displacements in the middle plane of the plane in the x, y and z directions, respectively, and so, are tjj^ and t^^ with respect to the rotations about the X and Y axes. All five displacements and rotations in equation 3.1 can be expanded as a Taylor series, in the direction of the width

PAGE 69

59 Figure 3 . 1 Laminated composite beam

PAGE 70

of the plate. If just the first order terms of this series are retained, the displacements and rotations can be expressed in the form of Ug {x,y) -C7lx) +yF{x) Vo(x,y) +yG(x) w U.y) -W{x) +ye{x) (3.2) i|r^{x,y) -
PAGE 71

Denoting the differentiation with respect to x by a prime, we represent the deformation field from Eq. 3.7 as e XX 1 Y xy [y ZX) \ (j)' + ya' (^x°] + z yxy° + Z ^ xy J [ 0 j Y ° I 0 ) U'^yF' F+ V' The deformation field can also be expressed as E'-'E+yE where (e ) ( U' ] {f'\ y xyO V'+F 0 ! (E) ; a' ^xy a-e' 0 ^y ZX, U+e'j (3.8) (3.9) (3.10) From the constitutive relations the stresses developed in the k**' layer are given by ^"'\-[Q]\^''°Uz[Q] xy) V « xyoj (3.11) and ^zxQss (+^') + y c>55 (a+e') (3.12) where [Q] %1 Ql2 ^16^ ^12 Q22 Q26 Pl6 Q2S OsS, (3.13)

PAGE 72

62 The terms Q.. in Eq. 3.13 are stiffness terms and can be found in Whitney (1987) . By definition, the force (N) and moment (M) acting on the laminate are ' 2 2 {N^,M^)f X^[l,z]dz (3.14) 2 Ji 2 ' 2 and they are shown in Fig. 3.2. The stiffness coefficient terms are defined as {A^j,B^j,D^j)f Q,j[l,z,z^] dz (3.15) _ J3 2 By using equations 3.9-13 and the definition of stiffness coefficients, the Eq. 3.14 can be written in a more simplified form as F-C{E+y^) (3.16) where (F) ^{N^,N^,M^,M^,V^) (3.17)

PAGE 73

63 Figure 3.2 Force and moment resultants

PAGE 74

and the matrix C is formed with the following stiffness coefficient terms ^16 fill file 0 ^ A,s ^ee file fiee 0 fill file ^11 ^le 0 file fiee ^le ^ee 0 1° 0 0 0 The factor k presented in term C55 is the shear correction factor (Whitney, 1987), and it is considered equal to 5/6 in the present work. 3.3 New Set of Force and Moment Resultants A new set of force and moment resultants is defined from the integration of the column vector of forces along the width of the beam: 2 FI Fdy (3.19) 2 and b 2 (3.20) In the two last equations, the vector of forces F is given by equation 3.17, and b is the width of the beam. Substituting equation 3,16 in Eq. 3.19-3.20, and assuming that the stiffness C of the beam is constant, we obtain

PAGE 75

65 F' f C Edy-bCE -b/2 biz Ff C E y^dy -i/2 •(«) CE The explicit beam contitutive relations are: (3.21) N [F]xy M xy V \ xi -b ^11 ^16 -^11 -^16 ^16 -^6 ^16 ^66 •^ii -Bie -Oil ^16 ^16 ^66 -^16 -^66 0 0 0 0 0 icM, 55j a-Q' (3.22) and 0 0 V \ X/ 12 ^11 -^16 ^11 ^16 0 ^16 -^6 ^16 ^66 0 ^11 ^16 ^11 ^16 0 ^16 ^66 ^16 ^66 0 0 k^A, 55/ ' F' a' 0 0 (3.23) 3.4 Equilibrium Equations The Principle of Minimum Potential Energy is applied in the derivation of the equilibrium equations. The principle says that of all displacement fields which satisfy the prescribed constraint conditions, the correct state is that which makes the total potential energy of the structure a minimum (Tauchert, 1970) . The total potential energy of the structure (n) is obtained from the sum of the strain energy of the beam (0) and the potential of the external force (/) :

PAGE 76

66 71 (|) + X (3.24) According to the Principle of Minimum Potential Energy, where the symbol S is the variational operator. Considering only transverse loading acting on the beam surface in the z-direction, the potential relative to the external loading is given by where q(x,y) and w(x,y) are the transverse loading and displacement, respectively, and L is the beam length. Similar to the procedure adopted to define a new set of force resultants in Eqs 3.19-20, the transverse loading is divided into two parts and defined by Sti 6<|> + 6x 0 (3.25) b L 2 (3.26) 0 2 Qix,y) g{x) + y
PAGE 77

67 Substituting the definition of transverse displacement (Eq. 3.7) and transverse loading (Eq. 3.28) in Eq. 3.26 gives the result X -f{g(x)W(x) + (S(x)QU)) dx (3.29) From Eq. 3.16, the strain energy per unit area (0^) of the beam can be derived as <1)^ ^ E^CE T _ _ (3.30) ^ (E+yE) '^CiE+yE) The strain energy per unit length (0^) is obtained by integrating
PAGE 78

68 the Eq. 3.31 takes the form of 2 2 4( 2 (3.33) bE'^CE + E^C^I The total strain energy in the beam is expressed as <|> j^L^ 2Jl 12 ) (3.34) dx The total potential of energy in the whole structure (tt) is obtained by substituting the Eq. 3.29 and 3.34 into equation 3.24. bE'^CE+ — e'^CEq{x) W{x) ^(x) Q ix) X 2 dx (3.35) Applying the Principle of Minimum Potential Energy we obtain the following equilibrivim equations 6n dN, 6a 6£7 dx dN -0 ^-0 bx dV^ -+g-0 dx -N xy dN, dx bn bW bn bF bn 56 bn 54> bti dM^ ^ _ -V-M 0 dx ^ ^ diV^-M^) ^ ^ ^2— + 0 dx dM^ dx (3.36)

PAGE 79

69 Next, from each equilibrium equation in Eq. 3.36 results one corresponding natural boundary condition, which is shown in the next equations as dx -0 N^bU-0 bx -0 N^bV-0 dx — g V^bW-0 dx N^bF-0 dx (V,-M^) 66-0 dM, dx M^6(j)-0 dM, dx M^ba-0 (3.37) By using Eqs. 3.22-23, the new set of force and moment resultants can be expressed in terms of seven unknown functions: U, V, W, F,
PAGE 80

70 when solving the problem for displacements, but it is considered when calculating the strain energy release rate, in order to obtain a more accurate solution. The nodal forces and moments for the i*'' node of the structure are F^^, F^, F^, M^, M^, W and T. A quadratic variation of all seven displacements is assumed along the element length. Denoting by X any specific displacement, and by X. this specific displacement at the i*'' node, the displacement and deformation within the element are defined as x'-b^Xj^-^b^X^^b^X^ IJ.je; In the last equation, the terms a,, and aj are interpolation functions in the form of a^-l-3x+2}? a^-x{2x-l) (3.39) a^-4.x{l-2x) and b^, and b3 are their derivatives with respect to the x direction i^i(4x-3)/L jbj(4X-1)/L 2)3-4 (l-2x)/L (3.40)

PAGE 81

71 where In the last equation, L is the element length and x is the coordinate along the beam axis. By using Eq. 3.38, and the terms a, and b,. (i=l,2,3) defined in Eqs 3.39 and 3.40, the deformation field (Eq. 3.10) within each finite element as a function of the nodal displacements is expressed in the form of (e) [a] (g) ,3 v where q is the vector of the displacements U^.V^.W^,F^,^^,a^,Q^, , (3.42) and a and a are two strain-displacement matrices, both with dimension of (5 x 21) . The first strain-displacement matrix, a, is composed by the following terms: 00000 0 252 00000 0^)3 00000 0^ 0 2?, Oa^OO 0 0 0 32 00 0 Oi^jOajOO 0 0000 23iO 0 0000 2)2 0 0 0000 233 0 0 0 0 0 0 0 a, -2?, 0 0 0 0 0 a2 -2?2 0 0 0 0 0 53 -2)3 00 2?i0ai 0 0 00 2^2 0 32 0 0 00 2^3 0 33 0 0 (3.43)

PAGE 82

72 and the second one, a, of: 'OOOi?iOOO OOOijjOOO OOOJbjOOON 000000 0 000000 0 000000 0 00000 jb^O 00000 i^gO 00000 jb3 0 000000 0 000000 0 000000 0 0 0 0 0 0 i?^ 0 0 0 0 0 a2 0 0 0 0 0 2^3^ (3.44) similarly to Eq. 3.10, the total deformation within each element is given by the expression e-e+yS (3.45) or, after using the strain-displacement matrices, by e(a + yi) q (3.46) With a-a+yi (3.47) the element stiffness matrix is defined as b L 2 0 b 2 (3.48) L 2 2

PAGE 83

73 In Eq. 3.48, D is the elasticity matrix for composite materials. The terms in matrix D are the same terms shown in Eq. 3.18. Each term of the matrix D is defined as b {Aij,B^j,D^j) -Jq^j [1,Z,z2] dz (3.49) a where a and b are respectively the lowest and hiqhest coordinates in the z direction in the considered element. By defining two element stiffness matrices in the form of L K'bfa'^D a dx ° ^ (3.50) K-^ (S-'D a dx 12 J 0 the equation 3.48 is reduced to K -K + K (3.51) The elements of the stiffness matrices K and K are shown in the Appendix. The dimension of the stiffness matrix k are reduced from (21 x 21) to (14 x 14) by using static condensation and the stiffness matrix of each finite element is assembled in a global stiffness matrix, representative of the stiffness of the whole structure. The final representation of the problem is F K q (3.52)

PAGE 84

74 where K is the total stiffness matrix, and F and q are vectors with the nodal forces and displacements, respectively. In the present work, the equation 3.52 is solved using the Gauss Elimination Method (Bathe, 1982) . 3.6 Strain Energy Release Rate The strain energy release rate is calculated by using the expression for J-integral. A zero volume path is delineated surrounding the crack tip, as shown in Fig. 2.8, and the Jintegral expressed in terms of the strain energy per unit length (U) of each element around the crack, as in Eq. 2.55, The strain energy per unit length of the beam (Eq. 3.34) is given by The expression for the strain energy per unit length of the i*" element is obtained by substituting Eq. 3.46 in the previous equation. (3.53) (3.54) (3.55)

PAGE 85

In Eq. 3.55, the vector q contains the displacements of the nodes of the i^*^ element, including the middle node. The reason for considering the middle node in the calculation of the strain energy is because the static condensation process assumes that no forces or couples are applied to the intermediate node. This assumption is not true for the case where bottom and top offset elements are put together to model the beam. 3 . 7 Rigid and Gap Elements The possibility of interference between the crack surfaces in the delaminated region, when the beam is under loading, principally under torsion loading, requires that gap elements or rigid elements be placed in appropriate positions in order to monitor the contact. Two types of rigid elements are anticipated to be used depending on which side, left or right, of the beam presents the interference phenomenon. Figure 3.3 illustrates the interference situation and the rigid element actions. In matrix form, the equations of the rigid element placed on the left side can be represented as

PAGE 86

76 '0 0 0 0 1 ' 0 0 0 0 b 2 0 0 0 0 -1 0 0 0 0 _ b 2 1 i. b 2 -1 _ b 2 0 / (F ) . 0 , (3.56) where t and b represent top and bottom elements, respectively, and the force F is shown in Fig. 3.3. Similarly, the equations of the rigid element placed to on the right side are '0 0 0 0 1 0 0 0 0 _ b 2 0 0 0 0 -1 0 0 0 0 b 2 1 _ b 2 -1 b 2 0 > (F ) . 0 , (3.57) The gap element is placed in the structure to avoid the interference in the delaminated region while the structure is under vertical loading. In matrix form, the constraint equations relative to the gap element are: '0 0 1 ^ 0 0 -1 (F.)b .1 -1 0. . 0 , (3.58)

PAGE 87

77 interference crack t (a) rigid element t (b) rigid element — ^ ^ -F interference (c) bF 2 bF 2 rigid element action bF 2 rigid element action Figure 3.3 Rigid elements: (a) Beam under torsion (b) Interference on left side (front view) (c) Interference on right side (front view)

PAGE 88

78 The rigid element and gap element matrices in Eqs. 3.563.58 are assembled in the global stiffness matrix using the same procedure adopted to assemble the element stiffness matrices. In Eqs. 3.56-3.58 is the contact force between the interfering elements, or is the couple due to the contact force and F is the force in the rigid or gap element in the z direction.

PAGE 89

CHAPTER 4 NUMERICAL RESULTS AND DISCUSSION 4.1 Numerical Results This chapter presents the numerical results obtained from the finite element described in the previous chapter. Three different tests were performed, and their numerical results compared with those obtained from the beam theory. The first test of the beam finite element consists of a specially orthotropic cantilever beam subjected to several different end loading conditions. This test is divided into two parts: one in which the beam is modeled using two sublaminates (top and bottom) and the other, where the beam is formed by just the top sublaminate. The purpose of theses tests is to prove the validity of the finite element developed in modeling the beam problem, including the torsion problem. In the second numerical test, a specially orthotropic delaminated cantilever beam is put under torsion loading. Initially, the delamination is supposed to be in the middle plane of the beam, and the result in terms of the strain energy release rate is compared with the closed-form solution. The delamination is then placed between different layers along the thickness, and a study is performed to understand the 79

PAGE 90

80 effect of the delamination position on the strain energy release rate. The last numerical test consists in the simulation of a quasi-static impact test in specially orthotropic and quasiisotropic simply supported beams. Some experimental observations for this specific tests are available in our Composite Laboratory (Kwon, 1991) , and the present numerical study will help interpret the experimental results. Again, the focus was on the effect of the delamination position on the values of the strain energy release rate. 4.2 Material P roperties and Stiffness Coefficients In all the numerical tests performed in the present study, the beam is supposed to be made up of material with the following properties: Longitudinal Elastic Modulus (E^) : 14.00 GP^ Transverse Elastic Modulus (Ej) : 1.00 GP^ Shear Modulus (Gij) : 0.53 GP^ Poisson's Ratio (v^^) o.30 Poisson's Ratio (v^^) o.55 The above values are typical values for high performance graphite/ epoxy unidirectional composites. Table 4.1 lists the laminate stiffness coefficients for a specially orthotropic beam with thickness equal to 2 mm, and the reference X-Y plane considered to be in the middle

PAGE 91

81 (MIDDLE) , on the top (BOTTOM) or on the bottom (TOP) surfaces of the beam. 4.3 Cantilever Beam The dimensions of the beam are 0.2 m x 0.02 m x 0.002 m (length x width x thickness) , and the beam is subjected to a set of six different end loading conditions, as shown in Fig. 4.1 and Fig. 4.3. 4.3.1 Model with two sublaminates Initially, the beam subjected to the forces shown in Fig. 4.1 is modeled by using two connected sublaminates (top and bottom), as shown in Fig. 4.2. Table 4.2 presents the numerical and beam theory solutions, in terms of the displacements at the tip of the beam for each of the loading conditions. The expressions for beam theory results can be found in Timoshenko (1970), Reissman and Pawlik (1980) and Sankar (1991) , or can be derived from the beam equations presented in Sankar (1991) . In the following expressions the superscripts t and b stands for top and bottom, respectively. 4.3.1.1 Beam under Force 1 fF ^) In this case, we have just the U displacement. All other displacements are zero.

PAGE 92

82 Table 4.1 Stiffness coefficients 1 STIFFNESS COEFFICIENTS POSITION OF THE X AXIS TOP BOTTOM MIDDLE Ai1 0.2818173E8 0.2818173E8 0.2818173E8 ^6 0.0 0.0 0.0 0. 1066621E7 0. 1066621E7 0. 1066621E7 Bii 0.2818229E5 -0.2818229E5 0.0 16 0.0 0.0 0.0 0. 1066643E4 -0. 1066643E4 0.0 D„ 0.3757714E2 0.3757714E2 0,9394285E1 0.0 0.0 0.0 0. 1422219E1 0. 1422219E1 0.3555547 K2 A55 0. 8888511E6 0.8888511E6 0.8888511E6

PAGE 93

83 Figure 4.1 Loading conditions at middle plane (a) Perspective view (b) Front view

PAGE 94

84 Figure 4.2 Model with 42 elements

PAGE 95

85

PAGE 96

86 L (4.1) b {A \^^A 4.3.1,2 Beam under Force 2 (F ..) Under Force 2 the beam has the displacements V and F. The remaining five displacements are zero. bHA\^^A\^) b{A'=,^^A^,^) (4.3) F--F, b^{A \^*A \^) 4.3.1.3 Beam under Force 3 fF ,) From this case results the displacements W and
PAGE 97

87 12L (4.7) 4.3.1.5 Beam under Force 5 fF^ ) Under this loading condition the displacements are W and 0, and all others are zero. 1 7" 2 (4.8) 2 b(D\,*D\,) (4.9)
PAGE 98

and while the displacement a is defined as 12 o-C,cosh(x) +C,sinh(x) -C, where 12 55 C2 — Cj^A,sinh(;<) and In Eqs. 4.13-4.16,

PAGE 99

89 (4.19) (4.20) D •* 66 4.3.2 Model with top sublaminate In a second test, the beam is modeled using top offset elements only (Fig. 4.3). The beam is subjected to the set of loading conditions displayed in Fig. 4.3. Table 4.3 shows the numerical and closed solutions for this test. The sources for the beam theory solutions presented in Table 4 . 3 are the same used before: Timoshenko (1970), Reismann and Pawlik (1980) and Sankar (1991b) . Some expressions are derived from the beam equations. 4.3.2.1 Beam under Force 1 (F ^) The tip displacements under F, are U, W and 0. The remaining displacements are zero. L (4.21) 11 W ^ (4.22)

PAGE 100

Figure 4.3 Loading conditions at the bottom plane (a) Perspective view (b) Front view

PAGE 101

91 cantilever beam Z A top offset elements Figure 4.4 Model with 21 elements

PAGE 103

93 (|) L (4.23) 4.3.2.2 Beam under Force 2 (F .) In the case of a cantilever beam under Fj the tip displacements are V, F, ^ and a. The displacement V can be represented as VV^^ (4.24) where V, (F2) 4L3 ^ 2L \ 3 a c n 55 J (4.25) and (4.26) The displacement 0 is necessary to calculate and is defined as e F, 2 66 1 + 55 / where the term D^jj is defined as (4.27) 55 12 55 (4.28)

PAGE 104

The displacement F is given by 94 6L^ 3 a t (4.29) 11 Finally, the expression for the displacement a is a qcosh(x) +C2sinh(z) -C3 48 (4.30) where C'l-n 48 D' 55 (4.31) C2 -qA,sinh(^) (4.32) and 55 66 In Eqs. 4.30-4.33, (4.33) X ' L 48 D\eD\^ (4.34) 12£>' 66 11 / (4.35) 55 12 66 (4.36)

PAGE 105

95 D • * 66 66 (4.37) and the term D%, is defined as 55 (4.38) 4.3.2.3 Beam under Force 3 (F ,) Under force Fj, the beam presents three displacements: U, W and 0. They are defined as 4.3.2.4 Beam under Force 4 CF^ ) The displacements in this case are V, F, a and 6. The remaining displacements are zero. This case can be considered as the result of the addition of two different loadings, I and II, as shown in Figure 4.5. (4.39) and 2 (4.41)

PAGE 106

96 h/2b The displacements for the Loading I were calculated as explained in the Section 4.3.1.4. From Loading II results the displacements a and 8, which are calculated using the following expressions (Sankar, 1991b) : (4.42) a-C8sinh(A.L) where Cg is a constant defined as ^ 12/?„ 1 (4.43) ^ jb^n Xcosh(AL) 1

PAGE 107

97 where the constant A is calculated using Eq. 4.18. From the definition given in Eq. 3.20, we have that 2 M,f M,dy -b/2 -2^4 2 2b h (4.44) The rotation 6 is given by where the constant Cj is defined in Eq. 4.16. Considering the boundary condition e(0)=0, we have e(x) —c^aL Loading I does not produce a or e displacements. Hence, the a and 9 displacements calculated using the previous equations are the result for the total loading. The displacement F is given by Eq. 4.7 and the displacement V is obtained using the expression v-v, +6 — where V, is given by Eq. 4.6. 4.3.2.5 Beam under Force 5 (F) Under force Fg the displacements are U, W and 0. The remaining displacements are zero, where the displacement


PAGE 108

98 U -^<1) (4.48) -fs^^ (4.49) The displacement W is defined as ^--^.-H^ (4.50) 4.3.2.6 Beam under Force 6 (F^) In this case the displacements are V, F, a and d. 2 (4.51) Where the rotation 6 is calculated using the expression e b D 66 1 + 66 4 D 55 / (4.52) The rotation F is defined as F tan-1 ( — ) L (4.53) and the warping function a can be obtained using Eqs. 4.304.37. 4.3.3 Comments on the finite element so Tut inn Tables 4.2 and 4.3 show a very good agreement between the present finite element solutions and the analytical results derived from beam theory. This agreement can also be seen in

PAGE 109

99 Figures 4.6 and 4.7, in which is shown the angle of twist along the beam length for two cantilever beams with different lengths. The beams are under a unit torque, and in both cases they were modeled by ten finite elements. From these figures, we can see that in the region close to the support the agreement is not good, and probably it will require use of larger number of finite elements near the fixed support. The beam theory solution for this case (Sankar, 1991b) incorporates the edge effect, represented by the third term in brackets on the right hand side of Eq. 4.52, and is given by the following expression where e(x) and e(L) are the angle of twist at a distance x from the fixed support and at the tip of the beam, respectively. The tip rotation of a beam, e(L), is given by e ix) -F, X QiL) (4.54) 4 2? A 66 (4.55) In Eqs. 4.54 and 4.55 we have that (4.56) X -

PAGE 110

100

PAGE 111

101 M * « N TO O O O O O o o d d o (pej)e 'UOUB^OU

PAGE 112

102 and 55 12 55 (4.57) ^66 ^66 + 12 55 where is a shear correction factor, considered equal to 5/6 in the present work. 4.4 Delaminated Cantilever Beam under Torsion In the present application, the beam is assumed to have a delamination of length 20 mm. The dimensions of the beam are 0.60 m x 0.02 m x 0.002 m. The length of the beam and of the delamination were enlarged if compared with the beam used in the previous examples. The reason was to avoid the edge effect, so as to be able to model the problem using fewer finite elements. A unit torque is applied at the free end of the beam, and the strain energy released is computed. The numerical result is then compared with the closed-form solution obtained using the beam theory (Sankar, 1991b). Figure 4.8, (a) and (b) , shows the beam under torsion and the model used to obtain the closed-form solution, while Fig. 4.9 illustrates the finite element model used to solve the problem. In this model, middle-plane beam finite elements are used in the uncracked regions, behind and beyond the delamination, (regions 1 and 3) , while top and bottom offset elements are used in the cracked portion of the beam (region 2). Later, the

PAGE 113

103 Figure 4.8 Delaminated beam (a) Under torsion (b) Separated by parts

PAGE 114

104 Figure 4.9 Finite element model

PAGE 115

105 delamination will be assumed in different positions along the thickness of the beam, the strain energy release rate will be computed for each of these positions. 4.4.1 Delamination at the middle-plane of the beam The finite element model used to solve this case has forty-six nodes and forty-six elements, as shown in Fig. 4.9. There are twenty-two regular finite elements and twenty-four offset elements (twelve top elements and twelve bottom elements) . Referring to Fig. 4.8(b) and to the strain energy release rate (G) definition ( Hellan, 1984) we have ~ dA~ b dL^ where U is the strain energy of the beam and A is the crack area. In the present case, the area of the crack is formed by a rectangle with dimensions egual to b and Lj. The strain energy of the beam (U) is defined as U^Te(L) (4.59) where T is the torque and 5(L) is the rotation e at the tip of the beam. The rotation at the tip of the beam is calculated by adding the rotations in the three parts of the beam shown in Figure 4.7, and is defined by

PAGE 116

106 (4.60) Where the subscripts stand for the regions in the beam (Fig. 4.8(b) ) . By definition, the superior and inferior limits of the beam in the z direction. The parts l and 3 of the beam have the same integration limits, h/2 and -h/2, while for part 2 they are h/4 and -h/4. Considering these integration limits and solving Eq. 4.61, we obtain the following relations among the stiffness coefficients of the three parts of the beam: By using Eqs. 4.54, 4.55 and 4.60, and neglecting the edge effect, the rotation at the tip of the beam is now (4.61) where the integration limits h^ and h. stand respectively for (4.62) expressed as e(L) n + + 4b [d,, D, 1 (4.63) 4 A 55

PAGE 117

Referring to Fig. 4.8(b), we have that 107 L1+L3-L-L2 Substituting Eq. 4.64 in Eq. 4.63 gives r(L-L,) I (4. 64) e(L) 4i) 1 . 1 (4.65) Using Eq.4.59 we obtain the strain energy of the beam (U) as + { ^ 1^66 ^55 j b 8Z?. (4.66) 55/ Considering Eq. 4.58, we obtain the following expression for strain energy release rate (G) G rp2 1 D, 55/ 2^6 ^D^^i, (4.67) or, after simplifying the last equation, 3 r2 G 8 b^n (4.68) 66 Substituting in the last equation the values of r-1 Nm b-0. 02m Z?66 -0-3556 -^m N _3

PAGE 118

108 we obtain for the strain energy release rate the value of 2,636.0 J/m^, while the numerical result for the present case is 2,545.5 J/m^. The relative difference between numerical and closed-form solutions is 3.45%. 4.4.2 Delamination position varying along the thickness In this example, the delamination was assumed to be in different positions along the thickness of the beam. The procedure adopted here is similar to that used to calculate the strain energy release rate with the delamination placed in the middle plane, except for the integration limits when calculating the stiffness coefficients A.., B.^. and D... Figure 4.10 illustrates the adopted model for the case of delamination placed between layers close to the top. Table 4.4 presents the obtained results for different positions of the delamination along the thickness of the beam. These results are also depicted in Figure 4.11. From Figure 4.11, we can conclude that the possibility of crack propagation will be increased as the crack position approaches the middle plane of the beam. As the crack plane approaches the top or bottom faces of the beam, the strain energy release rate becomes negative. This could be due to errors in the energy calculations, and can be avoided by having shorter elements near the crack tip. This problem should be subjected to further investigations .

PAGE 119

109 B : bottom offset element Figure 4.10 Delamination close to the top plane

PAGE 120

110 Table 4.4 Strain energy release rate (G) for cantilever beam POSITION G (J /m } between layers 1 and 2 -17.5 between layers 2 and 3 173.2 between layers 3 and 4 653.8 between layers 4 and 5 1,409.9 between layers 5 and 6 2,195.1 between layers 6 and 7 2,545.5 between layers 7 and 8 2,195.5 between layers 8 and 9 1,409.9 between layers 9 and 10 653.8 between layers 10 and 11 173.2 between layers 11 and 12 -17.5

PAGE 121

4 Q. O o >s _ £ Q. o 1? v> 21. n ^ / to (M c « M 3 O ID d d (M E u I< u (0 < III _l 111 cc >o cc 111 z 111 < cc I»O 0> CO to CM TNoiivNii/Mviaa do noiivoot

PAGE 122

112 4.5 Simply Supported Beam under Transverse Force In this section we consider the problem of transverse loading on a simply supported beam (Figure 4.12). The dimensions of the beam, the length of the delamination, the material and the finite element model are the same as in the preceding example. A force equal to 100 N is applied to the center of the beam. The finite element model is shown in Figure 4.12. Two classes of laminates were analyzed, a specially orthotropic laminate and a [0/45/-45/90]s laminate. The delamination was assumed in different positions along the laminate thickness, and the results for each class of laminate are listed in Table 4.5 and shown in Figure 4.13. Similarly to the cantilever beam case, from the present results we can conclude that the possibility of crack propagation will be increased as the crack position approaches the middle-plane of the beam. This phenomenon is less pronounced in the specially orthotropic beam than in the other beam. Another conclusion is that for a given transverse force, the possibility of crack propagation is more in the [0/45/-45/90]s laminate. These results have an important application in understanding the propagation of delamination damage in composite laminates due to low-velocity impact as well as quasi-static indentation types of loading (Kwon, 1991) .

PAGE 123

100 N 1 1 elements t 1 1 elements R 4 4 R T T >— < T >— ( T — « B B B B \ / R 4) < R 12 elements R : regular element T : top offset element B : bottom offset element Figure 4.12 Simply supported beam

PAGE 124

114 Table 4.5 Strain energy release rate (G) for simply supported beam 1 POSITION G (J/m^) Specially orthotropic [0/45/-45/90]s between layers 1 and 2 406.9 1,945.0 between layers 2 and 3 1,071.0 5,858.0 between layers 3 and 4 1,971.0 10,070. 0 between layers 4 and 5 2,495.0 11,390.0 between layers 5 and 6 1,971.0 10,070.0 between layers 6 and 7 1,071.0 5,858.0 between layers 7 and 8 406.9 1,945.0

PAGE 125

115

PAGE 126

CHAPTER 5 SUMMARY AND FUTURE RESEARCH 5 » 1 Summary The present investigation was restricted to delaminations in laminated beams. Finite elements were developed to model specific cases of delaminations, and the numerical results were compared with results from other methods and approaches. A laminated shear deformable beam finite element with nodes offset to either top or bottom has been developed to analyze free edge delaminations. These offset elements were used to model the sublaminates above and below the delamination plane. A rigid element was used to connect the nodes at the crack tip, and the forces and moments transmitted by this rigid element are part of the finite element solution of the problem. The strain energy release rate was calculated using expressions derived from the J-integral, the crack tip force method and crack closure method. These expressions are in terms of displacements and forces in the elements surrounding the crack-tip, which are part of the finite element solution. The present method considers the effect caused by the thermal stresses induced in the composite layers during the curing process. 116

PAGE 127

117 There are some advantages in using the present method. The finite element method has great flexibility in modeling any kind of geometry, loading and boundary conditions. This makes it easier to model common specimens used in laboratory experiments and compare results. In modeling the cracked region of the beam by using offset node elements, one can avoid renumbering the nodes as the crack propagates. It tremendously saves computer storage and time involved in redefining nodes. The delamination in anisotropic beams was analyzed using a new beam finite element to model the problem. The formulation of this new beam finite element is based on a beam theory for laminated beams in which the equilibrium equations are assumed to be satisfied in an average sense over the width of the beam. A new set of force and moment resultants for the beam were introduced from this assumption. Two problems of practical interest were solved using this method. One of them was the problem of a specially orthotropic delaminated cantilever beam under torsion loading, and the other was the case of specially orthotropic and quasi-orthotropic simply supported beams under transverse loading. In both problems the delamination was assumed to be in different positions along the thickness of the beam, and the strain energy release rate for each position was computed. The strain energy release rate was found to have the maximum value in the middle plane of the beam, for both the cantilever beam under torsion and the

PAGE 128

118 simply supported beam under transverse force. Practically, this indicates that the probability of crack propagation becomes higher as the position of the delamination approaches the middle plane of the beam. 5.2 Future research The present method has almost all its results validated by comparison with other methods and approaches. The method needs to be applied in an extensive body of experimental tests. The extension of this method to the study of buckling in laminated composite beams will add an important possibility in composite structural mechanics. One suggestion that naturally arises from this study is to extend the presently developed method to three-dimensional fracture problems, to compute strain energy release rate in delaminated plates. The study of plates with multiple delaminations can be an important sequel to this research. The final suggestion is the extension of this method to dynamic problems, which will bring significant benefits to both research in impact dynamics and design of composite materials. Then, this work will have achieved an important goal in the process of simulation of delaminations and crack propagation.

PAGE 129

APPENDIX STIFFNESS MATRIX The subroutine STIFF is used to form the stiffness matrix K of the beam finite element used in modeling anisotropic laminated beams. The matrix K results from the addition of two matrices, K and K. K K + K Figure A illustrates the formation process of matrix K and it is important for understanding the subroutine steps. The final result is the matrix C, with dimension of 14 x 14. The input to the subroutine STIFF are: B: width of the beam SV: vector with the stiffness terms for top, bottom and regular elements. XLE: element length IXl: integer that defines the element type: l for top element, 2 for bottom element and 3 for regular element. The output of the subroutine STIFF are: CC: matrix (7,7) (see Fig. A.l) CA: matrix (7,14) (see Fig. A.l) C: matrix (14,14), representative of the stiffness of the beam finite element. 119

PAGE 130

SUBROUTINE STIFF Input data matrix k CHAT') 14 15 21 1 14 15 21 CB ACB CAB CCB 1 14 15 21 CH ACH CAH CCH 1 14 15 21 C AC CA CC 120 14 14 Figure A Stiffness matrix K

PAGE 131

SUBROUTINE STIFF ( B , SV , XLE , IXl , CC , CA , C ) IMPLICIT DOUBLE PRECISION (A-H, 0-Z) DIMENSION SV(30) ,CB (14,14) ,CH( 14,14) ,C (14, 14) DIMENSION CI (14, 7) ,C2 (14, 14) ,CC(7,7) , AC (14, 7) ,CA(7,14) DIMENSION ACB(14,7) ,ACH(14,7) , CAB (7, 14) ,CAH(7,14) DIMENSION CCB(7,7) ,CCH(7,7) 11=10* (IXl-1) A11=SV(II+1) A16=SV(II+2) A66=SV(II+3) Bll=SV(II+4) B16=SV(II+5) B66=SV(II+6) Dll=SV(II+7) D16=SV(II+8) D66=SV(II+9) A55=SV(II+10) DO 1 1=1,14 DO 1 J=l,14 C(I, J)=0.0 CONTINUE DO 5 1=1,14 DO 5 J=l,14 CB(I, J)=0.0 CH(I, J)=0.0 C2(I,J)=0.0 CONTINUE DO 10 1=1,7 DO 10 J=l,7 CC(I,J)=0.0 CCB(I,J)=0.0 CCH(I, J)=0.0 CONTINUE DO 15 1=1,14 DO 15 J=l,7 AC(I,J)=0.0 CA(J,I)=0.0 ACB(I, J)=0.0 CAB (J, I) =0.0 ACH(I, J)=0.0 CAH(J,I)=0.0 C1(I, J)=0.0 CONTINUE 'BAR' MATRICES CB(14,14) CB(1,1)=A11*(70./XLE) CB(1,2)=A16*(70./XLE) CB(1,4)=A16*(-15.) CB(1,5)=B11*(70./XLE) CB(1,6)=B16*(-15. ) CB(l,7)=B16*(-70./XLE) CB(1,8)=A11*10./XLE

PAGE 132

CB(1,9)=A16*10./XLE CB(1,11)=A16*5. CB(1,12)=B11*10. /XLE CB(1,13)=B16*5. CB(1,14)=B16*(-10./XLE) CB(2,2)=A66*(70./XLE) CB(2,4)=A66*(-15. ) CB(2,5)=-CB(1,7) CB(2,6)=B66*(-15) CB(2,7)=B66*(-70./XLE) CB(2,8)=CB(1,9) CB(2,9)=A66*10./XLE CB(2, 11)=A66*5. CB(2,12)=-CB(1,14) CB(2, 13)=B66*5. CB(2,14)=B66*(-10./XLE) CB(3,3)=A55*(70./XLE) CB(3,5)=A55*(-15.) CB ( 3 , 10) =A55*10 . /XLE CB(3, 12)=A55*5. CB(4,4)=A66*4.*XLE CB(4,5)=CB(1,6) CB(4,6)=B66*4.*XLE CB(4,7)=-CB(2,6) CB(4,8)=-CB(1,11) CB(4,9)=-CB(2,11) CB(4, 11)=A66*(-XLE) CB(4,12)=-CB(1,13) CB (4 , 13 ) =B66* (-XLE) CB(4,14)=CB(2,13) CB(5, 6)=D16*(-15) CB (5, 5) =D11* (70 . /XLE) +A55*4 . *XLE CB(5,7)=D16*(-70./XLE) CB(5,8)=CB(1,12) CB(5,9)=-CB(1,14) CB(5,10)=-CB(3,12) CB(5,11)=CB(1,13) CB(5,12)=D11*10./ XLE+A5 5 * ( -XLE ) CB(5,13)=D16*5. CB(5,14)=D16*(-10./XLE) CB(6,6)=D66*4.*XLE CB(6,7)=D66*15 CB(6,8)=-CB(1,13) CB(6,9)=-CB(2,13) CB(6,11)=CB(4,13) CB(6,12)=-CB(5,13) CB(6,13)=D66*(-XLE) CB(6,14)=D66*5. CB(7,7)=D66*70./XLE CB(7,8)=CB(1,14) CB(7,9)=CB(2,14) CB(7,11)=-CB(2,13)

PAGE 133

CB(7,12)=CB(5,14) CB(7,13)=-CB(6,14) CB (7 , 14) =D66*10 . /XLE CB(8,8)=A11*70. /XLE CB(8,9)=A16*70. /XLE CB(8, 11)=A16*15. CB (8 , 12 ) =B11*70 . /XLE CB(8,13)=B16*15. CB(8, 14)=CB(1,7) CB(9,9)=CB(2,2) CB(9,11)=-CB(2,4) CB(9,12)=-CB(1,7) CB(9,13)=-CB(2,6) CB(9,14)=CB(2,7) CB(10,10)=CB(3,3) CB(10,12)=-CB(3,5) CB(11,11)=CB(4,4) CB(11,12)=-CB(1,6) CB(11, 13)=CB(4,6) CB(11, 14)=CB(2,6) CB ( 12 , 12 ) =D11*70 . /XLE+A55*4 . *XLE CB(12, 13)=D16*15, CB{12,14)=CB(5,7) CB(13,13)=CB(6,6) CB(13,14)=-CB(6,7) CB(14,14)=CB(7,7) DO 20 1=1,14 DO 20 J=I,14 CB(J,I)=CB(I, J) CONTINUE DO 18 1=1,14 DO 18 J=l,14 CB(I, J)=CB(I, J) *B/30. CONTINUE MATRIX ACB(14,7) ACB(1 , 1 )=All*(-80. /XLE) ACB(1 r2 l=A16*(-80. /XLE) ACB(1 ,4 i=A16*(-20. ) ACB{1 r5 =Bll*(-80. /XLE) ACB(1 ,6 )=B16*(-20. ) ACB(1 =B16*(80./XLE) ACB(2, i; =ACB (1,2) ACB(2, 2 =A66*(-80. /XLE) ACB(2, 4, =A66*(-20. ) ACB(2, 5] =-ACB (1,7) ACB(2, 6] =B66*(-20. ) ACB ( 2 , 7] =B66*(80./XLE) ACB ( 3 , 3] =A55*(-80. /XLE) ACB ( 3 , 5) =A55*(-20. ) ACB (4, 1) =A16*20.

PAGE 134

ACB(4,2)=A66*20. ACB (4 , 4) =A66*2 . *XLE ACB(4,5)=B16*20. ACB(4 , 6) =B66*2 . *XLE ACB(4,7)=B66*(-20.) ACB{5,1)=ACB(1,5) ACB(5,2)=-ACB(1,7) ACB(5, 3)=A55*20. ACB(5,4)=-20.*B16 ACB (5,5) =D11* (-80 . /XLE) +A55*2 . *XLE ACB(5,6)=D16*(-20,) ACB(5,7)=D16*80. /XLE ACB(6, 1)=ACB(4,5) ACB(6,2)=-ACB(4,7) ACB(6,4)=ACB(4,6) ACB(6,5)=D16*20. ACB (6 , 6) =D66*2 . *XLE ACB(6,7)=D66*(-20. ) ACB(7,1)=ACB(1,7) ACB(7,2)=ACB(2,7) ACB(7,4)=-ACB(2,6) ACB(7,5)=ACB(5,7) ACB(7,6)=D66*(20. ) ACB (7 , 7 ) =D66* (-80 . /XLE) ACB(8,1)=ACB(1,1) ACB(8,2)=ACB(1,2) ACB(8,4)=ACB(4,1) ACB(8,5)=ACB(1,5) ACB(8,6)=ACB(4,5) ACB(8,7)=ACB(1,7) ACB(9,1)=ACB(1,2) ACB(9,2)=ACB(2,2) ACB(9,4)=ACB(4,2) ACB(9,5)=-ACB(1,7) ACB(9,6)=-ACB(4,7) ACB(9,7)=ACB(2,7) ACB(10,3)=ACB(3,3) ACB(10,5)=A55*20. ACB(11,1)=-ACB(4,1) ACB(11,2)=-ACB(4,2) ACB(11,4)=ACB(4,4) ACB(11,5)=-ACB(4,5) ACB(11,6)=ACB(4,6) ACB(11,7)=20.*B66 ACB(12,1)=ACB(1,5) ACB(12,2)=-ACB(1,7) ACB(12,3)=-ACB(5,3) ACB(12,4)=20.*B16 ACB ( 12 , 5) =D11* (-80 . /XLE) +A55*2 . *XLE ACB(12,6)=ACB(6,5) ACB(12,7)=ACB(5,7)

PAGE 135

ACB(13, 1)=-ACB(4,5) ACB(13,2)=ACB(4,7) ACB(13,4)=ACB(4,6) ACB(13,5)=-ACB(6,5) ACB (13,6) =D66*2 . *XLE ACB(13,7)=D66*(20.) ACB(14,1)=ACB(1,7) ACB(14,2)=ACB(2,7) ACB(14,4)=ACB(4,7) ACB(14,5)=ACB(5,7) ACB(14,6)=ACB(6,7) ACB(14,7)=ACB(7,7) DO 24 1=1,14 DO 24 J=l,7 ACB(I, J)=ACB(I, J) *B/30. CONTINUE .CAB (7, 14) DO 25 1=1,14 DO 25 J=l,7 CAB(J,I)=ACB(I,J) CONTINUE .CC(7,7) CCB(1, 1)=A11*160. /XLE CCB (1,2) =A16*160 . /XLE CCB(1,5)=B11*160./XLE CCB(1,7)=-B16*160./XLE CCB(2,2)=A66*160. /XLE CCB(2,5)=-CCB(1,7) CCB(2,7)=B66*(-160./XLE) CCB(3,3)=A55*160. /XLE CCB(4,4)=A66*16.*XLE CCB(4, 6)=B66*16. *XLE CCB (5 , 5) =011*160 . /XLE+A55*16 . *XLE CCB(5,7)=D16*(-160./XLE) CCB(6, 6)=D66*16.*XLE CCB (7 , 7) =066*160 . /XLE DO 28 1=1,7 DO 28 J=I,7 CCB(I, J)=CCB(I,J)*B/30. CONTINUE DO 26 1=1,7 DO 26 J=I,7 CCB(J,I)=CCB(I, J) CONTINUE 'HAT' MATRICES CH(14,14) BB=B**3/12. CH(4,4)=A11*70./XLE CH(4,6)=B11*70./XLE CH(4,11)=A11*10./XLE CH(4,13)=B11*10./XLE

PAGE 136

CH(6,6)=(D11*70./XLE)+(A55*4.*XLE) CH(6,7)=A55*(-15.) CH(6,11)=CH(4,13) CH(6,13)=(D11*10./XLE)+A55*(-XLE) CH(6,14)=A55*{-5.) CH(7,7)=A55*70./XLE CH(7,13)=-CH(6,14) CH(7,14)=A55*10./XLE CH(11,11)=CH(4,4) CH(11,13)=CH(4,6) CH(13,13)=CH(6,6) CH(13,14)=-CH(6,7) CH(14,14)=CH(7,7) DO 55 1=1,14 DO 55 J=I,14 CH(J,I)=CH(I,J) CONTINUE DO 60 1=1,14 DO 60 J=l,14 CH(I, J)=CH(I, J) *BB/30. CONTINUE .ACH(14,7) ACH(4,4)=BB*(All*(-80./XLE))/30. ACH(4,6)=BB*(Bll*(-80./XLE))/30. ACH(6,4)=BB*(Bll*(-80./XLE) ) /30. ACH(6 , 6) =BB* (Dll* (-80 . ) /XLE+A55*2 . *XLE) /30 ACH(6,7)=BB*(A55*20.)/30. ^^^/-^u. ACH(7,6)=BB*{A55*(-20.) ) /30. ACH(7,7)=BB*(A55*(-80./XLE))/30. ACH{ll,4)=BB*(All*(-80./XLE))/30. ACH(ll,6)=BB*(Bll*(-80./XLE))/30. ACH(13,4)=BB*(Bll*(-80./XLE))/30. ACH(13 , 6) =BB* (Dll* (-80 . ) /XLE+A55*2 . *XLE) /30 ACH(13,7)=-BB*(A55*20.)/30. ' ACH(14,6)=BB*(A55*20.) /30. ACH(14,7)=BB*(A55*(-80./XLE) ) 730. .CAH(7,14) ' DO 65 1=1,14 DO 65 J=l,7 CAH(J,I)=ACH(I, J) CONTINUE .CCH(7,7) CCH(4,4)=BB*(A11*160.'/XLE) /30. CCH(4,6)=BB*(B11*160./XLE) /30. CCH(6,6)=BB*( (011*160. /XLE)+(A55*16.*XLE) ) 730 CCH(7,7)=BB*(A55*160./XLE)/30. ^LE))/30,

PAGE 137

DO 70 1=1,7 127 DO 70 J=I,7 CCH(J,I)=CCH(I, J) 70 CONTINUE C ADDING THE 'BAR' AND 'HAT' MATRICES DO 30 1=1,14 DO 30 J=l,14 C(I, J)=CB(I, J)+CH(I, J) 3 0 CONTINUE DO 32 1=1,14 DO 32 J=l,7 AC(I, J)=ACB(I, J)+ACH(I, J) CA(J,I)=CAB(J,I)+CAH(J,I) 32 CONTINUE DO 34 1=1,7 DO 34 J=l,7 CC(I, J)=CCB(I, J)+CCH(I, J) 34 CONTINUE C STATIC CONDENSATION CALL INVDET(CC,7) CALL MATMUL(AC,CC,C1,14,7,7) CALL MATMUL(C1,CA,C2,14,7,14) DO 40 1=1,14 DO 40 J=l,14 C(I,J)=C(I,J)-C2(I,J) 40 CONTINUE CALL INVDET(CC,7) RETURN END -

PAGE 138

REFERENCES Agarwal, B. D. and Broutman, L. J., Analysis and Performance of Fiber Composites. John Wiley & Sons, New York, NY, 1980. Armanios, E. A. and Badir, A., Hvarothermal Influence on the Edge De lamiantion in Composites . Technomic Publishing Company, Lancaster, Pennsylvania, 1989 Bathe, K. J., Finite Element Procedures in Engineering Analysis, Prentice-Hall, Englewood Cliffs, New Jersey, 1982. Chan, W. S. and Ochoa, O. 0., "An Integrated Finite Element Model of Edge-Delamination Analysis for Laminates Due to Tension, Bending, and Torsion Loads," Proceedings of AIAA/ASME/ASCE/AHS/ASC 28^^ SDM conference, Monterey, California, 1987, pp. 27-35. Hellan, K., Introdu ction to Fracture Mechanics ^ McGraw-Hill, New York, New York, 1984. Hu, S., Sankar, B. V. and Sun, C. T. , "A Finite Element Study of Dynamic Delamination Growth in Composite Laminates," Proceedings of AIAA/ASME/ASCE/AHS/ASC 30**' SDM conference, Mobile, Alabama, 1989, pp. 1250-1255. Kwon, Y.S. , "Indentation-Flexure Damage in Graphite/ Epoxy Laminates," Master Thesis, University of Florida, Gainesville, Florida, 1991 Pagano, N. J., "On the Calculation of Interlaminar Normal Stress in Composite Laminate," J. Composite Materials. Vol. 8, 1974, p. 65. Pagano, N. J. , Interlaminar Res ponse of Composite Materials . Elsevier Science Publishing Company, New York. New York. 1 Q O O ' , I. S., Q3DG A Computer Program for Strain Energy Release Rate for Delamination Growth in Composite Laminates," NASA CR 178205 (1986), 103 p 128

PAGE 139

129 Reismann, H. and Pawlik, P. S., Elasticity Theory and Applications . John Wiley & Sons, New York, 1980. Sankar, B. V., "A Finite Element for Modeling Delaminations in Composite Beams," Computers & Structures, Vol. 38, N° 2, pp. 239-246, 1991a. Sankar, B. V., "A Beam Theory for Laminated Composites and Application to Torsion Problems," AeMES TR-91-1-01, 1991b. Sankar, B. V. and Pinheiro, M. A., "An Offset Beam Finite Element for Fracture Analysis of Delaminations," Proceedings of AIAA/ASME/ASCE/AHS/ASC 31^*" SDM conference, Long Beach, California, 1990, pp. 1227-1233. Tauchert, T. R. , Energy Principles in Structural Mechanics . McGraw-Hill, New York, 1974. Timoshenko, S. P. and Goodier, J. N. , Theory of Elasticity . Mc Graw-Hill, New York, 1970. Usman, M. and Gandhi, M. V., "The Effect of Moisture and temperature on the Vibrational Characteristcs of Composite Materials," Proceedings of the American Society for Composites, 4*'' conference, 1989, pp. 555-579. Whitney, J. M., "Stress Analysis of a Mode I Edge Delamination Specimen for Composite Materials," AIAA Journal, Vol. 24, N° 7, 1986, pp. 1163-1168. Whitney, J. M., Structural Analysis of Laminated Anisotropic Plates, Technomic Publishing Company, Lancaster, Pennsylyania, 1987.

PAGE 140

BIOGRAPHICAL SKETCH Marco Antonio Santos Pinheiro was born on August 3, 1947, in Rio de Janeiro, Brazil. He is married to Fernanda Elizabeth Pinto Pinheiro and has two children, Simone, 20, and Luis Marco, 17. He joined the Brazilian Army in 1966, entering the Agulhas Negras Military Academy as a cadet. He graduated in 1969 and was commissioned as lieutenant in the Army Ordnance Corps . In 1975 he entered the Military Institute of Engineering, located in Rio de Janeiro, and received the bachelor's degree in mechanical engineering in 1977. As a result of the grade points average achieved, he received the "Henry Ford Award," given by the Ford Brazil Company. This award consists of a specialization in automotive engineering, with visits to automotive plants and research centers during a period of six months . After being technical director and leader of a military plant in the industrial city of Sao Paulo for two years, working and participating in the process of development and fabrication of armored cars and other heavy equipments for the Army, he was invited to be an instructor at the Military 130

PAGE 141

131 Institute of Engineering, in the city of Rio de Janeiro. He joined the faculty at this Institute in 1980 and taught undergraduate level courses in the field of automotive engineering until 1989. In the same period, he taught undergraduate courses at a private engineering school, the University Gama Filho, also in the city of Rio de Janeiro. These courses were in the area of mechanical materials and designs. He started his master's degree program in 1985 at the same Institute, defending his thesis on "Dynamic Analysis of Armored Cars" in October of 1987. After receiving his master's degree in mechanical engineering, he was selected by the Brazilian Army to pursue doctoral studies at the University of Florida, in the Department of Aerospace Engineering, Mechanics and Engineering Science. . '