vivlvav V v u .

J1_—1u'3 3311131 '33'3

    
   

 
   
  

3.3 33.3.3“

 

[TI 3.3“ 1 u 1

|l -,43' 1'1. 1.331313
‘1' 11311333331 113121511." 111* 311' 3"
1? 13.1111311111111111IL 331" 1 '1f
111.13% "; 11;?112‘111113112133111‘1 1111111

' 4 1

331

wfihflmflian
33 1.3 13'333'II333331‘13u‘3;

 

.94.
o
4.

a4
\p._a
.336 1"

.v

1.13 .
.Fwwwv3

 

.311.‘

   

     

    

      
    
    
   

     

   

WI.”

1.‘7 111%

331111.313 1111333333333” 1333'3‘3’3" '3333 3333333333'33333'13
I!”

M

11’
1'1

‘w‘3
33Ifi

11111.11. 1: 1 111..

"(Mug 33" 3.3". 1.33 331333333 .113 ”3333' 33333333

     
   
     
         
  
 
  
 
   
  
   
 
 

313%} 3.331% 11"3
11.13:: gm” ’33

3'“ 33' 19.33333.

H'IIIIE

3 333113"
1.3)" |11II {IIIII 3‘3

“333:

1‘; . {1.1M

313.331.. MY 3 311111.. i 3.3313 3' 3:: '3'33 "'

'33 333333 #3313“ ‘i "’ 333333'33'33 -' .’ 1 3333'3i " 3' 3-3“ 3 33 ‘3}:

14" _; :‘3 [431:3 " 3333' “3333‘ 3"3333 II, “1" 3:4 33?:“33:

‘nhuwm‘*w33” 1: " a}

I)! 334; IIIIII'II 3"" 2|. .141 3331. 3‘333 P333]: 8;“

:3 33.33.31. 33.3; 3.31,: ‘1'“ 3 3 "1 . iiz!‘ “33;? 33333353 3333 £313.13“, :3“ .

132 “33331133‘13'33'331'3‘333' 1.33 .‘ 12' 3313 I o 315;“; 31:33: :33. ‘ gi-

111111.133 3""1’ . 1“" ' ‘1" 11:3: 1‘ " 3111.31“ . .1131

3513:3333 1'3‘13 ' “3&3“ 1. 31131;“?53811.1“11‘133‘11IR.133'5'I21. 333433.?

".3333 3"333 i 1 ‘11:” '1‘." ‘21:. inf . 3: 333333 1:13 Pkg LI £151“; 123"

“35W 3 1 : ,L,&m Iflfigth
, “l. :I’E,r;1v2;' '3: ‘ZI 3

 
 
  
  
  

'2‘!” «=43

 

33.3.3333

., 1. .
3‘13\I333‘.33333

2.!
1 22°: :‘1' :.

u
.7

VS}

 

3... g: \\1. vi“ w- -—-1

   
 

11 1'1

3 3- A3
1 41.111.143.113
1' '1' C -

 

.1.

  
           
  
  
  
       
      
     
  
       
      

.-. “A ’1'."
:w no-

' rv».
. . .4. ._:v1
1’

.fi

11..

       
  
      
 
     
  

| I‘? V I
13313:}!
\ v3 ‘ -I.‘
33:3 33 3 ;. I7”

‘: t

w
".0

—-

  

A... 1-

‘ 3. 311.13%
.1' .; .1.” .
:31". ...,‘...II
"13'.v‘1';.1".
I V . .
113191~=v

'3‘ll .‘fl .g‘N‘nf‘;

‘0

 
   

..._.
. .

.31'3‘3

I”. J
I ,1 .
.

    
     

':f

1' ’ I??? .
"‘11 .113 111’ L
'1'- 33.4.3 33 h 3 3 F3353
W ,L I ngy~
I3a$mbE.
. 2:19 R‘fi .I‘
"3133M."
:31"?
131-

 

 
  

 

-:
$75.: :1
:719‘3:
a: y
.::r:".':‘.‘ .-

11F

. . A
l-V {:1

.1 1...“

        

  

    

  

   

  

‘1“ '
.3 ”a ‘1'.) '3'.
.934‘9

“0

 
  
  
 
 
   
 

 

{1.03,
:3 .I‘n'ud I
2 '93}: 1|
. 3'! $32.01" .
3,3333. . 3 '. 331k“ ’
.l V VI} 4: “.
{1° ‘21,:

   
 
   
 
    
 
   

 
 
 

 
   
 
  

  
    
  
  

   
        

       
       

   
     
   

  
    
  

  

     
 
   

    

  

-‘. é" -;
1.2.11 .. -~ 3.33:9; I 333331331333?
'L3t 5.1 .- 3"3'3333’3 "333. .1 . ‘I‘t‘flvjéjry
"1 1133111 ' 3 333-33333?
1'3'33'13i3’i'3‘2.1;3333'33.3 1: 13“:'3T‘3",!, 1.5-. 1131‘ 33533333.”?
1311 1.1331 {3111111.ng11 1.121111111111111“.
“Mm%hflfl1 3 “@3J1
3333333 :33 "‘ H3151 . {L312
.93.”? 33} {1393333MZ-2133141I‘H3‘ ‘éI ,
31"? . 1 .1 @133
1.1 111‘“ "19' ‘

33.131111...“

.1 . 1:21.113
33 33 13313333131311;

  

THtSl3

STUNATE

MMMMMMMM MMMMMMMMMMMMMM MMMMMMM

31293 01682 6269

 

 

 

 

 

This is to certify that the
dissertation entitled
A C0 Zig-Zag Finite Element for Analysis of
Damaged Laminated and Woven Composites

presented by

Venkateshwar Rao Aitharaju

has been accepted towards fulfillment
of the requirements for

Ph . D . degree in Mechanics

 

 

Mam

Major professor

Date ////3/?8

MSU is an Affirmative Action/Equal Opportunity Institution 0.12771

 

 

LIBRARY

Michigan State

Unlverslty

 

 

PLACE IN RETURN BOX
to remove this checkout from your record.
To AVOID FINES return on or before date due.

 

DATE DUE

DATE DUE

DATE DUE

 

D—Ec5é f 21309

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1M CICIMMM

 

 

A (:0 ZIG-ZAG FINITE ELEMENT FOR ANALYSIS
OF DAMAGED LAMINATED AND WOVEN COMPOSITES
By

Venkateshwar Rao Aitharaju

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Materials Science and Mechanics

l 997

ABSTRACT

A C0 ZIG-ZAG FINITE ELEMENT FOR ANALYSIS
OF DAMAGED LAMINATED AND WOVEN COMPOSITES

By

Venkateshwar Rao Aitharaju

This dissertation describes the research work carried out towards
developing analytical and numerical models for accurate analysis of composite laminates
and woven composites.

Three different higher order theories (cubic, quadratic, and partly cubic)
including zig-zag effects and transverse shear continuity are presented for analysis of
laminated composites. These theories are cast in such a way that the unknowns in the
model will involve only displacements and rotations at the top and bottom surfaces of the
laminate. It is advantageous to formulate the theories in this way as it provides the
opportunity to refine the associated finite element model by stacking sublaminates in the
thickness direction. In addition, a new transverse normal strain is derived by assuming
the transverse normal stress to be constant through the thickness and using Reissner's
mixed variational principle. After a comparison of results obtained from the above
theories, a quadratic zig-zag theory is recommended for general lamination schemes

based on accuracy and computational effort in the finite element formulation.

A new CO finite element for analysis of laminated beams based on the
quadratic zig-zag theory is developed. This element satisfies all the requirements of a
computationally efficient finite element model for analysis of multilayer laminates,
namely 1) the number of degrees of freedom is independent of the number of layers; 2)

the degrees of freedom are displacements and rotations; 3) no shear correction factor or

penalty numbers are present in the formulation. Consistent shear fields are developed for
the present finite element formulation to eliminate shear locking and thus extend the
applicability of the present element to extremely thin applications. It is observed that the
new formulation developed for transverse normal strain eliminates Poisons’s ratio
locking. Numerical results demonstrate that this new element is robust, accurate and
computationally efficient for analysis of multilayer laminates.

For analysis of composites with damage a new formulation, referred here
as the Variable Zig-Zag Kinematic Displacement Model (VZZKD) is developed. In the
present model the zig-zag fields are varied in an element. Numerical results Show that
present model can predict the response of composites with delamination very accurately.
But in the case of ply-damage, the results are very much dependent on the damaged layer
thicknesses and their placement in the laminate. Overall the present element holds great
promise for analysis of laminated composites with damage.

In the case of woven composites, a new analytical/numerical model which
is simple in terms of modeling and efficient in computational effort is presented to
estimate the effective stiffness properties. The model assumes the fiber geometry to be
sinusoidal and accounts for the wavineSs in both the directions. The problem of
estimating three dimensional stiffness properties is divided into three subproblems and
the superposition method is used. The results for the stiffness properties demonstrate that
the present model, which is both simple and computationally efficient, can give very
accurate results compared to a complex three dimensional finite element models.

The elements that are developed for the analysis of laminated composites
and woven composites are implemented into the MARC commercial finite element
program through a user defined subroutine (U SEELM). This demonstrates that the
present special elements developed for analysis of laminated composites and woven

composites are favorable for implementation in general purpose finite element programs.

Dedicated to

my beloved Parents

ACKNOWLEDGEMENTS

The present work is carried out under the guidance of Professor R.C.Averill. I shall be
forever indebted to him for his technical and professional guidance during the course of
this work. I express my sincere gratitude for the financial support provided during the
present research work.

Special thanks are also given to Professor D. Liu, N. J. Altiero and P. M. Duxbury for
serving on my Ph.D committee and for providing technical guidance during the course of
research. I would like to thank my CMRG members Y.B. Cho, Y.C. Yip, L.M. Smith, D.
Eby for the valuable discussions I had with them. I would also like to thank Bharat and
their family and all my friends in Spartan Village.

Finally, my deepest appreciation goes to my parents and my wife, for their support and

encouragement during this work.

TABLE OF CONTENTS

LIST OF TABLES ............................................................................................................. ix

LIST OF FIGURES ............................................................................................................. x
CHAPTER 1

INTRODUCTION ........................................................................................... l

1.1 Introduction ................................................................................................. 1

1.2 Analysis of Laminated Composite Panels .................................................. 3

1.2.1 Literature Review on Analysis of Thick Laminated Composite Plates

.......................................................................................................... 3

1.2.2 Goals and Objectives ...................................................................... 9

1.3 Analysis of Woven Composites .................................................................. 9

1.3.1 Literature Review on Analysis of Woven Composites .................. 10

1.3.2 Goals and Objectives .................................................................... 13

1.4 Layout of the Thesis .................................................................................. 13
CHAPTER 2

HIGHER-ORDER ZIG-ZAG THEORIES FOR ANALYSIS OF

LAMINATED COMPOSITE PANELS .............................................. 16

2.1 Introduction ............................................................................................... 16

2.2 Multilayer Beam Geometry and Notation ................................................. 17

2.3 Higher Order Zig-Zag Theories ................................................................ 17

2.4 Assumed Displacement Functions ............................................................ 23

2.5 Approximation for the Normal Stress, ..................................................... 24

2.6 Results and Discussion ............................................................................. 26

2.6.1 Eleven Layer Sandwich Panel ....................................................... 27

2.6.2 Composite Armored Vehicle (CAV) panel ................................... 29

2.7 Preliminary Conclusions ........................................................................... 30
CHAPTER 3

FINITE ELEMENT ANALYSIS OF LAMINATED COMPOSITE BEAMS...

............................................................................................................ 5 1

3.1 Introduction ............................................................................................... 5 1

3.2 Formulation of the Theory ........................................................................ 52

3.2.1 Multilayer Beam Geometry .......................................................... 52

3.2.2 Displacement Field ....................................................................... 52
3.2.3 Strain-displacement and Constitutive Relations ........................... 55
3.2.4 Approximation for the Normal Stress, ......................................... 56
3.3 Finite Element Formulation ...................................................................... 58
3.3.1 Finite Element Model ................................................................... 58
3.3.2 Consistent Shear Strain Fields ...................................................... 59
3.3.3 Poisson’s Ratio Stiffening .............................................................. 61
3.4 Results and Discussion ............................................................................. 61
3.4.1 Five Layer Unsymmetrical Laminate ........................................... 62
3.4.2 Composite Armored Vehicle (CAV) Panel ................................... 66
3.5 Preliminary Conclusions ........................................................................... 67
CHAPTER 4
ANALYSIS OF LAMINATED COMPOSITES WITH LOCAL DAMAGE...
............................................................................................................ 94
4. 1 Introduction ............................................................................................... 94
4.2 Formulation of the VZZKD model ........................................................... 97
4.3 Finite Element Formulation ...................................................................... 99
4.4 Results and discussion ............................................................................ 102
4.4.1 SimplySupportedBeanKO/90/0)SubjectedtoSinusoidallyVaryingLoad
103
4.4.2 Cantilever Beam Subjected to Point Load (Cross-Ply Laminate) 105
4.5 Preliminary Conclusions ......................................................................... 108
CHAPTER 5
THREE DIMENSIONAL ANALYSIS OF WOVEN COMPOSITES ....... 138
5.1 Introduction ............................................................................................. 138
5.2 Fabric Geometry - Configuration in Unit Cells ...................................... 139
5.3 Effective Moduli Theory ......................................................................... 143
5.4 Finite Element Analysis for Effective Stiffness Properties ..................... 146
5.5 Results and Discussion ........................................................................... 149
5.5.1 Fiber Volume Fraction ................................................................ 150
5.5.2 Effective Moduli and Poisson’s Ratio Estimation - Effect of Wavi-
ness Ratio ..................................................................................... 150
5.5.3 Effective Moduli and Poisson’s Ratio Estimation - Effect of Width
to Thickness Ratio (alt) .............................................................. 152
5.6 Preliminary Conclusions ......................................................................... 152
CHAPTER 6
CONCLUSIONS AND RECOMMENDATIONS ...................................... 169
6.1 Conclusions ............................................................................................. 169
6.2 Recommendations for Future Work ........................................................ 173
APPENDICES

vii

_H§J-‘

APPENDIX A .................................................................................................................. 176
APPENDIX B .................................................................................................................. 180

BIBLIOGRAPHY ....................................................................................... 18 1

viii

Table 2.1:

Table 2.2:

Table 2.3:

Table 2.4:

Table 2.5:

Table 2.6:

Table 3.1:

Table 3.2:

Table 3.3:

Table 3.4:

Table 4.1:

Table 4.2:

Table 5.1:

Table 5.2:

LIST OF TABLES

Material Properties of Eleven Layer Sandwich Panel .................................... 32

Lay-Up Sequence of Eleven Layer Sandwich Panel ...................................... 32
Normalized Center Deflections of Simply Supported Beam under Sinusoidal

Load (eleven layer sandwich panel) ............................................................... 33

Material Properties of TACOM Composite Armored Vehicle (CAV) Panel ..33

Lay-up Sequence of CAV Panel ..................................................................... 34
Normalized Predicted Center Deflections of CAV Panel by Various

Theories ......................................................................................................... 34
Material Properties for Five Layer Unsymmetrical Laminate ....................... 69
Lay-Up Sequence for Five Layer Unsymmetrical Laminate .......................... 69
Material Properties of TACOM Composite Armored Vehicle (CAV) Panel ..70

Lay-up Sequence of CAV Panel ..................................................................... 70
Material Properties of (0/90/0) Laminate ..................................................... 110
Nondimensionalized Center Deflections of Simply Supported Beam

Subjected to E Sinusoidal Load with Various Damaged 90° Layer
Thicknesses h and Damage Ratios, A. ......................................................... 110

Material Properties of Tow and Resin [Whitcomb (1989)] .......................... 154

Effect of Waviness Ratio on Moduli and Poisson’s Ratio (Moduli and
Poisson’s ratios are normalized with respect to a 0/90/90/0) laminate) ....... 155

Figure 2.1.

Figure 2.2(a).

Figure 2.2(b).

Figure 2.3(a).

Figure 2.3(b).

Figure 2.4(a).

Figure 2.4(b).

Figure 2.5(a).

Figure 2.5(b).

Figure 2.6(a).

Figure 2.6(b).

LIST OF FIGURES

Multilayer beam geometry and coordinate system. ................................... 35

Through-thickness variation of inplane displacement at the end of a
simply-supported beam (eleven layer sandwich panel, L/h=4). ................ 36

Through-thickness variation of inplane displacement at the end of a
simply-supported beam (eleven layer sandwich panel, Uhfl, constant
transverse normal stress assumption) ......................................................... 37

Through-thickness variation of inplane normal stress at the center of a
simply—supported beam (eleven layer sandwich panel, L/h=4). ................ 38

Through-thickness variation of inplane normal stress at the center of a
simply-supported beam (eleven layer sandwich panel, IJh=4, constant
transverse normal stress assumption) ......................................................... 39

Through-thickness variation of transverse shear stress at the end of a
simply-supported beam (eleven layer sandwich panel, Uh=4). ................ 40

Through-thickness variation of transverse shear stress at the end of a
simply-supported beam (eleven layer sandwich panel, L/h=4, constant
transverse normal stress assumption) ......................................................... 41

Through-thickness variation of transverse normal stress at the center of a
simply-supported beam (eleven layer sandwich panel, Llh=4). ................ 42

Through-thickness variation of transverse normal stress at the center of a
simply-supported beam (eleven layer sandwich panel, L/h=4, constant
transverse normal stress assumption) ......................................................... 43

Through-thickness variation of transverse normal strain at the center of a
simply-supported beam (eleven layer sandwich panel, Llh=4). ................ 44

Through—thickness variation of transverse normal strain at the center of a
simply-supported beam (eleven layer sandwich panel, IJh=4, constant
transverse normal stress assumption) ......................................................... 45

Figure 2.7.

Figure 2.8.

Figure 2.9.

Figure 2.10.

Figure 2.11.

Figure 3.1.
Figure 3.2.

Figure 3.3

Figure 3.4.

Figure 3.5.

Figure 3.6.

Figure 3.7.

Through-thickness variation of inplane displacement at the end of a
simply-supported beam (CAV panel, L/h=4, constant transverse normal
stress assumption). ..................................................................................... 46

Through-thickness variation of inplane normal stress at the center of a
simply-supported beam (CAV panel, L/h=4, constant transverse normal
stress assumption). ..................................................................................... 47

Through-thickness variation of transverse shear stress at the end of a
simply-supported beam (CAV panel, Llh=4, constant transverse normal
stress assumption). ..................................................................................... 48

Through-thickness variation of transverse normal stress at the center of a
simply-supported beam (CAV panel, L/hfi, constant transverse normal
stress assumption). ..................................................................................... 49

Through—thickness variation of transverse normal strain at the center of a
simply-supported beam (CAV panel, L/h=4, constant transverse normal

stress assumption). ..................................................................................... 50
Multilayer beam geometry and coordinate system. ................................... 71
Topology of finite element. ........................................................................ 72

Sublarninate stacking schemes used for five layer unsymmetrical
laminate (LAMl). ...................................................................................... 73

Normalized transverse deflections at the center of simply-supported beam
(five layer unsymmetrical laminate). ......................................................... 74

Through-thickness variation of inplane displacement at the end of simply-
supported beam (five layer unsymmetrical laminate, L/h=4). ................... 75

Through-thickness variation of transverse displacement at the center of
simply-supported beam (five layer unsymmetrical laminate, L/h=4). ....... 76

Through-thickness variation of inplane normal stress at the center of
simply-supported beam (five layer unsymmetrical laminate, L/h=4). ....... 77

Figure 3.8.

Figure 3.9.

Figure 3.10.
Figure 3.11.

Figure 3.12.

Figure 3.13.

Figure 3.14.

Figure 3.15.
Figure 3.16.

Figure 3.17.

Figure 3.18.
Figure 3.19.

Figure 3.20.

Through-thickness variation of transverse Shear stress at the end of
simply-supported beam (five layer unsymmetrical laminate, L/h=4). ....... 78

Through-thickness variation of transverse normal stress at the center of
simply-supported beam (five layer unsymmetrical laminate, L/h=4). ....... 79

Through-thickness variation of transverse normal strain at the center of
simply-supported beam (five layer unsymmetrical laminate, L/h=4). ....... 80

Alternate sublarninate stacking scheme used for five layer unsymmetrical
laminate (LAM2) ....................................................................................... 81

Through-thickness variation of inplane displacement at the end of simply-
supported beam with various stacking schemes (five layer unsymmetrical
laminate, Uh=4) ......................................................................................... 82

Through-thickness variation of transverse displacement at the center of
simply-supported beam with various stacking schemes (five layer
unsymmetrical laminate, Llh=4). ............................................................... 83

Through-thickness variation of transverse shear stress at the end of
simply-supported beam with various stacking schemes (five layer
unsymmetrical laminate, Uh=4) ................................................................ 84

Through-thickness variation of inplane displacement at the end of simply-
supported beam (five layer unsymmetrical laminate, Llh=10). ................. 85

Through-thickness variation of inplane stress at the center of simply-
supported beam (five layer unsymmetrical laminate, L/h=10). ................. 86

Sublaminate stacking schemes used for CAV panel. ................................ 87

Through-thickness variation of inplane displacement at the end of simply-
supported beam (CAV panel, L/h=4). ........................................................ 88

Through-thickness variation of transverse displacement at the center of
simply-supported beam (CAV panel, L/h=4). ............................................ 89

Through-thickness variation of inplane normal stress at the center of
simply-supported beam (CAV panel, Uh=4). ............................................ 90

xii

Figure 3.21.

Figure 3.22.

Figure 3.23.

Figure 4.1
Figure 4.2
Figure 4.3

Figure 4.4

Figure 4.5

Figure 4.6

Figure 4.7

Figure 4.8

Figure 4.9

Figure 4.10

Figure 4.1]

Figure 4.12

Through-thickness variation of transverse shear stress at the end of
simply-supported beam (CAV panel, L/h=4). ............................................ 91

Through-thickness variation of transverse normal stress at the center of
simply-supported beam (CAV panel, Uh=4). ............................................ 92

Through-thickness variation of transverse normal strain at the center of

simply-supported beam (CAV panel, Uh=4). ............................................ 93
Topology of Finite Element ..................................................................... 111
Damaged and Undamaged Regions of a Laminate .................................. 112

Present Modeling Approach for Damaged and Undamaged Regions ..... 112

Tip-loaded Cantilever Beam (0/90/0) Beam with Damage in the Off-axis
Plies Near the Clamped End .................................................................... 113

Through-thickness Variation of Inplane Displacement in the Damaged
Region of Clamped Beam (0/90/0 Laminate) .......................................... 114

Through-thickness Variation of Inplane Displacement at the
Interface of Damaged and Undamaged Regions of Clamped Beam
(0/90/0 Laminate) ..................................................................................... l 15

Through-thickness Variation of Inplane Displacement in the Undamaged
Region of Clamped Beam (0/90/0 Laminate) .......................................... 116

Through-thickness Variation of Inplane Normal Stress near the Fixed End
in the Damaged Region of Clamped Beam (0/90/0 Laminate) ................ 117

Through-thickness Variation of Inplane Normal Stress in the Middle of
Damaged Region of Clamped Beam (0/90/0 Laminate) .......................... 118

Through-thickness Variation of Inplane Normal Stress in the Undamaged
Region of Clamped Beam (0/90/0 Laminate) .......................................... 119

Through-thickness Variation of Transverse Shear Stress near the Fixed
End in the Damaged Region of Clamped Beam (0/90/0 Laminate) ........ 120

Through-thickness Variation of Transverse Shear Stress in the Middle of
Damaged Region of Clamped Beam (0/90/0 Laminate) .......................... 121

xiii

Figure 4.13

Figure 4.14

Figure 4.15

Figure 4.16

Figure 4.17

Figure 4.18

Figure 4.19

Figure 4.20

Figure 4.21

Figure 4.22

Figure 4.23

Figure 4.24

Figure 4.25

Through-thickness Variation of Transverse Shear Stress in the Undamaged
Region of Clamped Beam (0/90/0 Laminate) .......................................... 122

Through-thickness Variation of Inplane Displacement in the Damaged
Region of Clamped Beam ((0/90)4s Laminate) ........................................ 123

Through-thickness Variation of Inplane Displacement at the Interface
of Damaged and Undamamged Regions of Clamped Beam ((0/90)4S
Laminate) ................................................................................................. 124

Through-thickness Variation of Inplane Displacement in the Undamaged
Region of Clamped Beam ((0/9O)4s Laminate) ........................................ 125

Through-thickness Variation of Inplane Normal Stress near the Fixed End
in the Damaged Region of Clamped Beam ((0/90)4s Laminate) ............. 126

Through-thickness Variation of Inplane Normal Stress in the Middle of
Damaged Region of Clamped Beam ((0/90)4s Laminate) ....................... 127

Through-thickness Variation of Inplane Normal Stress in the Undamaged
Region of Clamped Beam ((0/90)‘;s Laminate) ........................................ 128

Through-thickness Variation of Transverse Shear Stress near the Fixed
End in the Damaged Region of Clamped Beam ((0/90)4s Laminate) ...... 129

Through-thickness Variation of Transverse Shear Stress in the Middle of
Damaged Region of Clamped Beam ((0/90)4s Laminate) ....................... 130

Through-thickness Variation of Transverse Shear Stress in the Undamaged
Region of Clamped Beam ((0/90)4S Laminate)......................... ............... 131

Tip-loaded cantilever beam (0/90)4s beam with dalamination in the Center
of Beam near the clamped end ................................................................. 132

Through-thickness Variation of Inplane Displacement in the Delaminated

Region of Clamped Beam ((0/90)4s Laminate) ........................................ 133
Through-thickness Variation of Inplane Displacement at the Interface
of Damaged and Undamamged Regions of Clamped Beam ( (0/

90),,s Laminate) ....................................................................................... 134

xiv

Figure 4.26

Figure 4.27

Figure 4.28

Figure 5.1
Figure 5.2
Figure 5.3

Figure 5.4

Figure 5.5
Figure 5.6
Figure 5.7

Figure 5.8

Figure 5.9
Figure 5.10

Figure 5.11

Figure 5.12

Figure 5.13

Through-thickness Variation of Inplane Displacement in the Undamaged
Region of Clamped Beam ((0/9O)4s Laminate) ........................................ 135

Through-thickness Variation of Inplane Normal Stress near the Fixed End
in the Delaminated Region of Clamped Beam ((0/90)4s Laminate) ........ 136

Through-thickness Variation of Inplane Normal Stress in the Undamaged
Region of Clamped Beam ((0/90)4s Laminate) ........................................ 137

Various Fiber Architectures of Woven Composites ................................. 156

Divisions of Unit Cells of Plain Weave and 5-Hamess Satin Weave ....... 157

Fiber Geometry in Region A, Band C .................................................... 158
Fiber Geometry in an Element in Region-B ............................................ 159
Fiber Geometry in an Element in Region-C ............................................ 160
An Element in Multilayer Laminate ........................................................ 161
Unit Cell of Symmetrically Stacked Plain-Weave ................................... 162

Variation of Volume Fraction of Fiber in the Plain-Weave Composite with
Waviness Ratio ......................................................................................... 163

Various Finite Element Meshes Used for Analysis of Unit Cell ............. 164
Variation of Inplane Modulus (E11) of Plain-Weave with alt Ratio ....... 165

Variation of Inplane Shear Modulus (G12) of Plain-Weave with

alt Ratio ................................................................................................... 166
Variation of Poisson’s Ratio ((912) with alt Ratio .................................. 167
Variation of Poisson’s Ratio (V 13) with alt Ratio ................................... 168

XV

CHAPTER 1
INTRODUCTION

1.1 Introduction

Fiber reinforced composite laminates have excellent mechanical properties such as
high specific strength and specific stiffness in the fiber directions. Initially, the use of lam-
inated composites was limited to the applications where inplane properties are important.
But as the advancement of manufacturing technologies and materials continues, compos-
ite laminates are finding applications where three dimensional properties are important. In
the third direction, i.e., thickness direction, laminated composites have poor mechanical
properties and are prone to interlaminar delamination. The compressive strength of lami—
nated composites with delarninations is many times less than the original strength. Besides
the effect of delamination, the response of thick composites is dependent upon transverse
effects. For efficient design of these laminated composites structures, a rigorous theory is
required to predict the delamination and transverse effects accurately.

One way of reducing the delamination type of failure is to provide reinforcement in the
thickness direction. Towards this, woven composites have been developed. The compos-
ites thus formed have good properties in mutually orthogonal directions as well as more

balanced properties and better impact resistance than unidirectional laminates. The ability

of these composites to drape and conform to irregular shapes makes them especially
appealing. The stiffness and strength behavior of woven fabric composites depends on the
fabric architecture. A number of parameters are involved in determining fabric architec-
ture such as type of weave, density of yarns, characteristic of fibers and matrix, factors
introduced during weaving such as crimp angle, etc. Many weave architectures are possi-
ble and some understanding of the behavior of composites as a function of weave architec-
ture is helpful in selecting an efficient weave for a specific application. Hence analytical
models are necessary to study the effect of various parameters on the behavior of woven
composites and to select an efficient fabric architecture.

The fiber reinforced composite laminates can sustain much higher loads after the onset
of local damage because of material property directionality. The direction of the principal
stress may not coincide with the direction of maximum strength. Also the highest stress
may not be the stress governing the failure of the material. An estimate of ultimate
strength of composite laminates is essential for efficient design of these structures.

From a theoretical point of view, the prediction of interlaminar shear stresses critical
for delamination failure, the selection of an efficient architecture of a woven composite for
a given strength and stiffness and the prediction of progressive failure of composites, can
be managed using three-dimensional analysis where each layer is modeled using conven-
tional three dimensional brick elements. But such an approach quickly becomes computa-
tionally expensive for structures having a large number of layers. Therefore it is highly
desirable to have a plate/shell model which can capture the local variations of displace-

ment and stresses through the thickness of laminate.

The goal of the present research is to develop analytical/numerical tools to analyze
laminated composites and woven composites. First, analysis of laminated composites
using a new layer-wise theory is attempted. A new analytical/numerical method is devel-
oped to estimate the effective stiffness properties of woven composites. For both cases,
new finite elements are developed and incorporated in the MARC finite element program
[MARC (1994)].

A brief review of literature on the analysis of laminated composites and woven com-
posites is presented below. Also, the differences among the present approach for analysis

of laminated composites and woven composites with the existing methods are highlighted.

1.2 Analysis of Laminated Composite Panels

A brief literature review for analysis of laminated composite panels and goals of the

present research are given below.

1.2.1 Literature Review on Analysis of Thick Laminated Composite Plates

Theoretical development in analysis of thick composites has attracted the attention of
many researchers. The theories developed for analysis of laminated composites are
divided into two categories; 1) equivalent single layer theories; 2) discrete layer-wise the-
ories. In equivalent single layer theories, the heterogeneous laminate is replaced with a
statically equivalent, Single, heterogeneous layer. In discrete layer-wise theories, a unique
displacement field is assumed in each layer and suitable continuity conditions for dis-

placements and stresses are enforced at the ply interfaces.

The theories belonging to equivalent single layer approach are developed by expand-
ing the displacement field in a polynomial series through the thickness. First, second and
third order theories have been proposed [see Yang, et a1. (1966), Whitney and Sun ( 1973)
and Reddy (1984)]. The first order shear deformation theory (FSDT) by Yang, et al.
(1966) is an extension of first order shear deformation theory for isotropic plates [Mindlin
(1951)] to composite laminates. FSDT allows a constant transverse shear strain through
the thickness and needs a shear correction factor to adjust transverse shear energy. This
shear correction factor is dependent on the lamination scheme and layer properties for
accurate results [Whitney (1973)]. Reddy (1984) developed a simple third order shear
deformation theory satisfying zero shear traction on the top and bottom surfaces of the
laminate.

The desired properties of a computationally efficient finite element for the analysis of
multilayer laminates are: 1) the number of degrees of freedom in the model must be inde-
pendent of the number of layers; 2) the assumed displacement model must be of Co conti-
nuity and simple to use; 3) no arbitrary shear correction factors must be present in the
formulation; 4) the element must pass the patch test; and 5) nodal degrees of freedom must
be standard translations and rotations (i.e., no higher-order or generalized degrees of free-
dom that do not have clear physical meaning). The finite'element formulation based on
first order shear deformation theory (Pryor and Barker, 1971) satisfies all the requirements
of the desired model except for the determination of the shear correction factor for a given
lamination scheme. The finite element model based on Reddy’s third order theory [Phan
and Reddy (1985)] does not usually need a shear correction factor, but C1 continuity is

required for the transverse displacement. The finite elements having displacement conti-

nuity more than C0 are not ideal for general purpose finite element programs. These equiv-
alent layer theories predict sufficiently accurate global structural response, i.e., transverse
deflections, vibration frequencies and buckling loads. But when the layer properties vary
drastically through the thickness, the results from these theories are often inaccurate due
to the increased presence of transverse shear effects. The reason for this discrepancy in
results is due to the C1 continuity of inplane displacements through the thickness assumed
in the model, leading to double-valued interlaminar shear stresses at the lamina interfaces.
To overcome the drawback of discontinuous shear stresses at the lamina interfaces, dis-
crete layer- wise theories have been proposed. These theories are divided into two catego-
ries depending on the number of unknowns in the model. The theories wherein the number
of unknowns in the model is dependent on the number of layers [c.g. Srinivas (1973);
Robbins and Reddy (1993); Toledano and Murakarni (1987a); Li and Liu (1992)] are often
prohibitively expensive for multilayer laminates.

Robbins and Reddy (1993) proposed a generalized laminate plate theory which could
account for any degree of approximation of the distribution of the inplane and transverse
displacement through a proper selection of variables and functions. Based on this theory,
he showed that finite elements with Lagrange and Hermite interpolation functions through
the thickness can be developed. The developed finite element model satisfies all the
requirements of a desired model except that the number of unknowns in the model is not
independent of the number of layers. Even though these discrete layerwise theories are the
most accurate for modeling discrete layer effects, they are expensive to use for a real struc-

tural analysis where computational efficiency is an important factor.

Discrete layer theories in which the number of unknowns in the model does not
depend on the number of layers in the laminate are described in [e. g. Murakarni (1986); Di
Sciuva (1986); Liu and Li (1992)].

Liu and Li (1996) showed that the zig-zag form of inplane displacement distribution
through the thickness cannot be represented by simply increasing the order of the polyno-
mial series as is usually done in equivalent single layer theories. Hence a new class of dis-
crete layer-wise theories has been developed wherein a piece-wise linear displacement
function is superimposed over a linear displacement field [Di Sciuva (1986); Murakarni
(1986)] or a cubic displacement field [Di Sciuva ( 1987); Cho and Parmerter (1993);
Toledano and Murakarni (1987b)] through the thickness of the laminate. The piece-wise
function in these theories is determined by imposing transverse shear continuity at the
lamina interfaces. These models are often called zig-zag theories. Toledano and Murakarni
(1987b) additionally assumed transverse shear stresses to be quadratic functions across
each layer and employed Reissner’s mixed variational principle. Di Sciuva (1986) devel-
oped a four node 32 degrees of freedom rectangular finite element for the analysis of mul-
tilayered composite laminates based on a linear zig-zag theory. The transverse
displacement is made C1 continuous across the element edges. Averill (1994), based on a
first order zig-zag theory, developed a C0 finite element and circumvented the Cl continu-
ity requirement by incorporating the concepts of interdependent interpolations [Tessler
(1981)] and penalty functions [Reddy (1980)]. Carrera (1996) developed a C0 Reissner-
Mindlin element which satisfies all the properties of a desired model for analysis of multi-
layer laminates based on a first order zig-zag theory. In this formulation, two different

fields along the laminate thickness direction are assumed for displacements and transverse

shear stresses, respectively. Di Sciuva (1993b) formulated a four node 40 degree of free-
dom rectangular plate based on the cubic zig-zag theory. Herrnitian functions were used to
approximate the transverse displacement. Rao and Meyer-Peining (1990) developed a
mixed bending element based on the Toledano and Murakarni model (1987b). Averill and
Yip (1996) developed a CO finite element based on cubic zig-zag theory using interdepen-
dent interpolations for transverse displacement and rotations and penalty function con-
cepts.

In developing finite element models, the displacement shape functions are usually
taken in the form of polynomials satisfying well-known conditions to ensure convergence,
completeness and continuity conditions [Reddy (1993)]. Unfortunately, elements derived
rigorously based on these basic paradigms can behave in unreasonably erratic ways in
many important situations. These difficulties are mostly pronounced in the low order finite
elements. The problems encountered are Shear locking, Poisson locking, etc.,. Many rem-
edies have been tried and some of them resulted in acceptably accurate elements. These
remedies include reduced integration (Zienkiwiecz et al. 1971; Hughes et al. 1978), mixed
formulation (Lee and Pian 1978), addition of nonconforming modes (Bathe 1976), energy
balancing (Fried 1974), etc. Some of these violate the well-known norms for finite ele-
ment formulation, but are accepted because the elements thus formulated are more accu-
rate than rigorously formulated elements. Averill and Reddy (1992) presented an
analytical technique for identifying the exact form of shear constraints that are imposed on
an element when its side-to-thickness ratio is very large and developed a robust four node
plate element based on Reddy’s third order theory. The cause of poor bending behavior of

shear deformable elements in thin plate or beam applications was identified by Prathap

(1994) through field consistency arguments. He introduced consistency paradigms for the
displacement functions and developed robust locking-free finite elements.

From the previous literature review it can be seen that:

1. The discrete layerwise theories are the most accurate compared to equivalent single
layerwise theories and the discrete layer theories wherein the number of degrees of
freedom is independent of number of layers (zig—zag theories) are the likely choice for
analysis of multilayer laminates.

2. Many of the zigzag layer-wise theories require C1 continuity requirement for finite ele-
ment analysis.

3. Only a few models are available having a CO continuity of displacements. Cl continuity
requirement is often circumvented using interdependent interpolations and penalty
function concepts.

4. Shear locking and Poisson’s ratio stiffening are the problems encountered in developing
displacement based finite elements. Except for the field consistency concepts of Prat-
hap (1994), all other tricks to circumvent locking violate well-known norms for finite
element formulation.

5. The finite element models developed based on the higher order theories are not compat-
ible with conventional plate/shelllbeam elements available in commercial packages.

6. Conventional plate/shell elements formulated from these higher order theories are not
convenient in stacking in the thickness direction, when refinement in the thickness
direction is desired. Multipoint constraints are often used which requires a lot of effort

in constraint formulation and implementation.

Keeping these aspects in mind, the goals and objectives of the present research with

respect to analysis of composites are given below.

1.22 Goals and Objectives

1. Development of a zig-zag discrete layer-wise beam theory and associated finite element
model wherein the number of degrees of freedom is independent of number of layers in
the laminate.

2. The developed zig-zag theory should give rise to a CO finite element formulation.

3. No shear correction factors or penalty numbers must be present in the formulation.

4. The element must be free from shear locking and Poisson’s ratio stiffening.

5. The element must be compatible with plate/shell/beam elements available in commer-
cial finite element programs.

6. The element must be able to capture three dimensional effects during the progressive

failure analysis of composites.

1.3 Analysis of Woven Composites

The effective properties of woven composite materials can be determined by analyzing
a unit cell. This type of model is popular because of its ability to capture the effects of
complicated fiber architectures in the unit cell. These models are based on an assumption
that a composite structural element can be formed by assembling the unit cells in all the

three directions. Hence, these unit cell models are valid for thick composites.

10

1.3.1 Literature Review on Analysis of Woven Composites

Ishikawa (1981) and Ishikawa and Chou (l982a and l982b, 1983) have proposed three
types of models to analyze woven fabric composites. These are the mosaic model, fiber
undulation model (crimp model) and the bridging model. These models are called classi-
cal models Since the basic assumption involved is that classical lamination theory is valid
for every infinitesimal piece of repeating unit cell of the woven fabric region.

In the mosaic model [Ishikawa (1981)] a fabric composite is idealized as an assem-
blage of pieces of asymmetrical cross ply laminates. The two dimensional extent of the
laminate is simplified by considering two one-dimensional models, i.e., a parallel model
and series model. In the parallel model, a constant strain state (iso—strain) is assumed in the
laminate and elastic stiffness matrices are derived which give upper bounds on in-plane
stiffness. In the series model or iso—stress model, the assumption of constant stress is made
and elastic stiffness matrices are derived. The assumption of constant stress leads to lower
bounds for in-plane stiffness of the fabric. In the above Studies the undulation of the fiber
is not considered.

In the fiber undulation model [Ishikawa and Chou (1983, l982a and l982b)], the
undulation of the fabric is considered in deriving the stiffness matrices. It is an extension
of the series model considered above [Ishikawa (1981)]. The undulation of the fibers in the
warp direction is not taken into account. This model is suitable for weaves with lower
order repeats in the unit cell. The in-plane stiffness constants obtained using this model are
lower than those obtained by the mosaic model (parallel model).

The bridging model [Ishikawa and Chou (1982b)] is an extension of the fiber undula-

tion model applied to general satin weave composites wherein, there are straight thread

11

regions surrounding an interlaced region and interlaced regions are thus separated from
one another. This model is a combination of series and parallel models. It is found that the
results obtained using the above model agree well with the experimental results for satin
weave composites.

Raju and Wang (1994) using classical laminate theory (CLT), obtained stiffness coeffi-
cients and thermal expansion coefficients for plain weave, 5 and 8 -harness satin weave
composites. The fiber geometry was assumed to be sinusoidal. They derived closed form
expressions for CLT stiffness matrices from which they evaluated elastic moduli, Pois-
son’s ratio and coefficients of thermal expansion.

The analytical models mentioned above are all based on classical laminate theory,
hence no prediction can be made for average transverse moduli, E33, 613 and average
Poisson’s ratio’s \713 and \723.

Naik (1994) developed a micromechanics model based on the iso-strain assumption
that discretely models the yarn architecture in the unit cell of a textile composite to predict
the overall mechanical properties. The analytical technique was implemented in a code
called TEXCAD. The calculated overall stiffnesses compared well with available 3-D
finite element results and test data for woven and braided composites.

Whiteomb (1989) used a three dimensional finite element model to analyze plain
weave composites. He found that the inplane moduli decreased linearly with increasing
tow waviness. Also, waviness was found to cause large normal and shear strain concentra-
tions in the composites when subjected to uniaxial load. He found some inconsistencies
between the finite element and experimental results for Poisson’s ratios. It is mentioned in

his paper that to reduce the number of elements in the model, compatibility in the center

l2

portion where matrix resin is present is not maintained. The accurate three dimensional
finite element modeling of woven composites requires enormous modeling effort and it
leads to a large number of elements when tow and resin constituents are discretely mod-
eled.

Foye (1989) used a fabric analysis method to predict the stiffness properties of a vari-
ety of simple and complex weaves. In this analysis, the unit cell of a woven composite is
divided into smaller rectangular subcells in which the reinforcing configurations are easier
to visualize, define and analyze. Each subcell is modeled as an inhomogeneous hexahedral
finite element. The effective properties are discontinuous in the subcell and numerical
integration is used to evaluate the stiffness matrices. After assembling the subcells, the
boundary conditions on all the faces of the unit cell are prescribed in terms of displace-
ments corresponding to six independent unit strain cases of 3-D elasticity. The stiffness
results thus obtained are found to be in close agreement with the experimental results.

Dasgupta, et al. (1990) analyzed plain weave composites for stiffness and thermal
properties using the homogenization technique. A two-scale asymptotic expansion was
used to relate the microscale behavior with macroscale behavior of the composite. The
microscale boundary value problem, defined over a periodic unit cell of the composite,
was solved using the finite element method.

From the above literature review it can be seen that,

1. The analytical methods take advantage of approximating the tow undulation by a sim-
ple geometry.
2. The analytical models which are used for analysis of woven composites are based on

either iso—strain or iso-stress assumptions.

l3

3. The analytical models used for analysis of woven composites are based on classical
lamination theory, hence no predictions can be made for transverse modulus and Pois-
son’s ratios.

4. The finite element models for woven composites are very complex and a rigorous mod-
cling effort is required for discrete modeling of tow and resin constituents.

5. The number of finite elements required for analysis are large as the fiber and resin are
discretely modeled. The computational effort involved is very huge.

6. Finite element models can not take advantage of tow geometry by simple continuous
functions.

Based on the literature review, the objectives of the present research for analysis of

woven composites are given below.

1.3.2 Goals and Objectives

1. Development of a model which combines the best features available in analytical and
numerical models for stiffness analysis of woven composites.

2. The model should be computationally efficient and should predict results with reason-
able accuracy.

3. The model should be easily extended for other fiber architectures.

1.4 Layout of the Thesis

The main goals of the present research are to develop analytical and numerical tools

for the analysis of laminated and woven composites. In the case of analysis of laminated

.Mw

l4

composites, three new zig-zag theories are developed, namely fully cubic, quadratic and
partly cubic, and are presented in Chapter 2. These new theories are assessed for their suit-
ability for finite element analysis of laminated composites. A new formulation has been
developed for the transverse normal strain by assuming the transverse normal stress to be
constant in the laminate to improve transverse normal strain and average transverse nor-
mal stress predictions. The accuracy and computational efficiency of the above theories
are evaluated by analyzing complex laminated cases.

A new finite element based on the quadratic zig-zag theory is developed and presented

in Chapter 3. In developing the finite element formulation, consistent shear fields are used
to eliminate shear locking and to make the element robust.
For damage analysis of composite laminates, a new scheme is proposed and implemented.
In the new scheme, the zig-zag displacement form is varied in an element and it is here
called a Variable Zig-Zag Kinematic Displacement (VZZKD) model. The formulation of
the model is presented in Chapter 4. The present model is tested for the analysis of lami-
nated beams with ply failure and delamination.

For analysis of woven composites, a new analytical/numerical model is developed
which has the best features available in analytical and numerical models developed by pre-
vious investigators. The new proposed model is discussed in Chapter 5. In the present
model, the fiber waviness and cross section are assumed to be sinusoidal. The waviness in
two directions is considered. The model is cast in a simple finite element model and the
results obtained for stiffness of the woven composite is compared with the results from the
complex 3-D finite element models. The developed model has lot of potential in progres-

sive failure analysis of woven composites.

15

The final chapter presents the conclusions of the present research. Recommendations
for future study are also presented in Chapter 6.
The tables and figures used to present the results of each chapter are presented imme-

diately after the corresponding chapter.

CHAPTER 2
HIGHER-ORDER ZIG-ZAG THEORIES FOR ANALYSIS
OF LAMINATED COMPOSITE PANELS

2.1 Introduction

The aim of the current study is to present and assess three new higher order shear
deformation theories (fully cubic, quadratic and partly cubic) including zig-zag effects.
Also, their suitability for finite element formulation and implementation is investigated.
These theories are based on a polynomial expansion of displacements through the thick-
ness of the laminate with a piecewise linear function superposed upon it. The maximum
order of the polynomial is limited to three in the inplane displacement field and one in the
transverse displacement to make these theories computationally viable. The fully cubic
theory includes terms in the series expansion of inplane displacements up to the cubic
term. The quadratic theory contains higher order terms up to the quadratic term in the
inplane displacement, and the partly cubic theory contains terms up to cubic with the qua-
dratic term absent. These models are designated as HZZT3, HZZT2, and HZZT3, respec-
tively. AS the transverse displacement effects through the thickness are small, the
transverse displacement is assumed to vary linearly across the thickness in all the three

models. In order to improve the predictions of transverse normal strains and their associ-

l6

l7

ated energies, the transverse normal stress is assumed to be constant through-the-thickness
of the laminate in a manner consistent with Reissner’s mixed variational principle [Reiss-
ner (1984)]. These models will be assessed by analyzing a simply supported beam sub-
jected to Sinusoidal load. The results obtained from these theories are compared with exact

elasticity solutions.

2.2 Multilayer Beam Geometry and Notation

The geometry, notation and coordinate system of the laminate consisting of N layers is
shown in Figure 2.1. The integer k denotes the layer number, with layer one located at the
bottom of the laminate. The origin of the coordinate system is taken at the bottom of the
laminate. The thickness of the beam is assumed to be h. 2“ denotes the distance from the
reference surface to the bottom of the kth layer. Each lamina (ply) is assumed to be ortho-

tropic and is perfectly bonded to adjacent plies. The length of the beam is assumed to be L.

2.3 Higher Order Zig-Zag Theories

The higher order zig-zag theories considered in the present study are based on the fol-
lowing assumed displacement fields:

Fully cubic zig-zag theory (HZZT3)

k-l
uirx. z) = ubrx) + wax) + z’er) + z3erx) + 2 (z — 20:10)
i= 1 (2.1)

u:(x, z) = (l — §)wb(x) + (%)w,(x)

18

Quadratic zig-zag theory (HZZT2)

k-l
ufitx, z) = ago) + zwb(X) + 22¢(x) + 2 (z - z.)&.(x>

i: 1 (1.2)
u:(x, z) = (I — %)wb(x) + Lama)
Partly cubic zig-zag theory (HZZT3)
k-l
u:(x, z) = abut) + zwb(x) + z36(x) + 2 (z — z,)§,.(x)
i: 1 (1.3)

uf(x, z) = (1 — 79’9””) + (3W1!)

where ukx‘ u“z are the displacements in x and 2 directions, respectively, in the layer k. ub,

111,, are the displacements and rotations at the bottom of the laminate. wb and wt are the
transverse displacements at the bottom and top of the laminate. o, 9 are the higher order
terms in the expansion of the inplane displacement field. The zig-zag functions 5,,- are
included only in the inplane displacement field.

The above higher order shear deformation theories are to be formulated in terms of
displacements and rotations at the bottom and top surfaces of the laminate. It is advanta-
geous to formulate the theories in this way as it provides an opportunity to refine the
model by stacking the sublaminates in the thickness direction. This advantage can be con-
veniently exploited in finite element analysis when refinement in the thickness direction is

desired for investigating thickness effects and also when the layerwise material properties

vary drastically [see Averill and Yip (1996) and Cho and Averill (1997)].

19

The zig-zag functions, 3;- , for each theory are determined in terms of other variables

by enforcing shear stress continuity at each layer interface:
1: = t atz = zk (1.4)

In case of HZZT2 and HZZT3 theories the shear strain expression is approximated by
assuming the contribution for shear strain from the transverse displacement is constant in

all the layers. That is, the shear strain is approximated as:

k Buk Buk

— ._x __Z
Y“ - Oz + 3x
k (1.5)
k _ BuJr 1 Bwb 8w,
7“ ’ 35 I 2(a_x E?)
Using Equation (2.4), the zig-zag functions E,- , are given by
HZZT3 theory:
A A A A aWb A 8W,
HZZT2 theory:
_ A ". .. 1 Bwb aw,
gi — ai+b,¢+ai§(E—+-5;) (17)
HZZT3theory:
_ A A A l awb 8W!
a, — ai+Ci9+ai§($+-a—x-) (18)

 

 

 

 

20
where
G,‘ i-l
a. =[ . —1][1 + 6k]
‘ 1+1
G k=l
Gi 1-1
8; = . -1 21+ 8k
1+1 ‘
G k=l
Gi i-l
a, =[ ”14134.3 2,] (1.9)
G k=l
(f Z- i—l
d, =( i+l-l]{(l—Z)+ 2d,.)
G k=l

 

Gi is the transverse shear modulus of the ith layer. From the Equations (2.6) to (2.9) it can
be seen that the zig-zag functions {3,- depend on the ratio of the shear modulus of adjacent
layers, and the zig-zag effects will be small if this ratio is small.

Substituting the zig-zag function i; for each theory into the inplane displacement
fields (Equations (2.1) to (2.3)), we can write the inplane displacement field as,

HZZT3 theory:

k-l k-l k-l

(u,")(x. z) = u.+v.[z+ 2 (z—z.)a.-]+¢[z’+ 2(2-2.)13.]+9[z3+ 2 (2-292,)
i=1 i=1 i=1

(1.10)

+ (86%,)[kit (z _ 2321,] + (35:31ng (2 - 29:3,]

i=1

21

HZZT2 theory:

k-l k-l
uxk(x,z) = ub+wb[z+ 2 (z-zi)&i]+¢[zz+ 2 (z—zi)l3i]+

i=1 i=1

$6311“. —‘¥M2<z- w]

i=1

(1.11)

HZZT3 theory:

k—l k-l
ukaCrZ) = ub+wb[z+ 2 (Z—Zi)ai]+e(z3+ Z (z—zi)2‘i)+

i=1 i=1

Mai”. M—‘E‘M’iu a)

i-l

(1.12)

The inplane displacement field for HZZT 3 involves O and 0. These functions can be

eliminated in favor of the displacement and rotation at the top surface of the laminate:

N
N dux
u = u ,

The inplane displacement field for HZZT3 involves derivatives of the transverse dis-

 

= v, (1.13)
z=h

 

placement and thus leads to a Cl theory.

The inplane displacement fields for HZZT2 and HZZT3 involve it) or 0 and
1 819,, 3w,
iii; t a:
(2.13). Thus the inplane displacement field for HZZT2 and HZZT3 does not involve any

). These terms can be eliminated by applying the conditions in Equation

derivatives of the transverse displacement and leads to a CO theory. Admittedly, the elimi-

3w 3w,
nation of the term 56—" + —) in favor of a new variable may be viewed as inconsis-

Bx 3x

22

tent. However, as the form of the assumed displacement field is somewhat arbitrary, no
variational crimes have been committed. Further, it will be shown that this approach leads
to a theory that is both computationally efficient and accurate.

The kinematic displacement field for all the theories can now be written in terms of
u b’ “1’ Vb, \Vt, wb’ Wt’ Wb, x and w,, xof the laminate. AS mentioned before, these
degrees of freedom facilitate through-thickness stacking of the sublaminates in finite ele-
ment analysis. The final form of the displacement fields for the above theories can be writ-

ten as:

u: = Ear IZY'IWBY 1,57 a =1’2 13:1:2 "Y =1,2 (114)
k__ k '

where the index a is used to denote the degree of freedom associated with inplane dis-
placement field (u, w), B is the index used to denote the transverse displacement field
(w, w,x) and 'y is the index used to denote the top and bottom faces of the laminate. For
example, 1111 is used to denote ub, 522 is used to denote w, and W22 is used to denote
(35—3). Summation on repeated indices is implied. 12”, 1’67 and K ’57 are the shape func-
tions associated with the thickness direction. They are given in Appendix A for all the
models. For HZZT2 and HZZT3 models the thickness shape functions, 11;? are zero for
all B and 'y.

The strains in the layer k can be calculated as

au" au" at." air"
k k k
exx = a—xx 'sz = 'a—zx'P—af 822 = a—Zz (1.15)

23

According to the principle of minimum total potential energy, the potential energy is

stationary at equilibrium i.e.,

OH = 5(U-l-V) = O (1.16)

Where U is the strain energy stored in the system and V is the work done by external

forces. The strain energy of the system can be written as,

N lit
1 k k k k k k
U = 25 M (a... e..+o.. saw... 7..) dx dz 0-17)
k=l ZK_l

The stress in the lamina k is given by

{o}" = iQi’iei" (1.18)

The coefficients in matrix Q are given in Appendix B. The potential energy of the applied

loads in the transverse direction can be written as

1,2 (1.19)

V = Jtlflwlfl dfl B
Q

where t” , t12 are the distributed transverse loads on the bottom and top surfaces of the

laminate, respectively.

2.4 Assumed Displacement Functions

To obtain Navier’s solution for the case of a simply supported beam under sinusoidal

load, the following functions are assumed for flow, Wm.

24

.. _ ~ it):

“at — ucW cos—
(1.20)

W - ii) sinzt—x

ly — l‘Y L

Substituting these functions into the potential energy and minimizing with respect to

it

a? and filly, we get a system of linear equations. These equations are solved for the

unknowns “at and w”.

2.5 Approximation for the Normal Stress, 022

In all the above formulations, the transverse normal strain is assumed to be constant.
For laminates wherein the layer properties vary drastically or in the case of a sandwich
panel with a thick, compliant core the transverse normal stress evaluated by the present
models is not predicted satisfactorily. To correct the distribution of normal stress and its
related energy in the laminate, a modification is proposed by assuming the transverse nor-
mal stress to be constant throughout the thickness of the laminate. The value of assumed
constant normal stress, 6'22 is calculated by using Reissner’s mixed variational principle

(Reissner, 1984), which yields the following condition:

N lit
2 j (atz—EZMZ = 0 (1.21)
k=1 zK_l

where 8:2 is the strain derived from the displacement field and is given by

k _ (wt—W17)

ezZ — h (1.22)

25

k . . . . . .
82: is the transverse strain in the layer k evaluated from the constitutive equation assum-

ing a constant variation of transverse normal stress, 6 and is given by

22’

.. k
kc Gzz Q12 k

ea = -—k-——kexx (1.23)
Q22 Q22
Substituting Equations (2.22) and (2.23) into Equation (2.21) and solving for ("ya , we get:
_ 1 817 BW
622 = §[Wr - Wb + (737)5117 +( a xB-J)Tfi7] (1.24)
where
N Zr 1
R = 2 J k dz
k = 1 Zr lQ22
(1.25)
N 2!: N Zr
Q1211: Q__i_2 Jk
=2 17" =2 1751‘!
k=1zK_ lQ22 k=lzx_ lQ22
using thiso zz , 8:: in each layer can be calculated from Equation (2. 23). The new normal

strain expression is given as:

 

k 1 (Wt-Wb) 1 (Barry) (5%?—
. = __ . ._ on, )
” Q’s. R ail 3‘ “I

+ leaf )(T-M-Q’le’évll

From the Equation (2.26), we can see that a; is no longer constant through the lami-

(1.26)

 

nate thickness but varies depending on the inplane and transverse stiffnesses and the dis-

26

placement field. The terms owing to the inplane displacement field in the transverse
normal strain are significant when the laminate is composed of layers with different trans-
verse modulus. This normal strain expression is used in the strain energy expression. The
procedure mentioned above leads to improved transverse normal strain energy and it is
observed to produce consistent results for normal strains and average transverse normal
stress. From a computational point of view, the above approach of assuming 0‘; to be
constant through thickness of laminate eliminates Poisson’s ratio locking in the associated

finite element model.

2.6 Results and Discussion

This section describes the results obtained from the three higher order theories and
comparison with the exact solutions (Pagano, 1969). Pagano’s solution has been modified
slightly to make it appropriate for beams as opposed to plates in cylindrical bending.
FSDT results are also obtained using Whitney’s shear correction factor (Whitney, 1973)
for the purpose of comparison. All the stresses are evaluated using the constitutive rela-
tions.

AS mentioned before, zig-zag effects are dominant when there is a difference in shear
modulus of adjacent layers. To study these effects, an eleven layer sandwich panel and a
panel from the US. Army composite armored vehicle (CAV) are used in the present study.
These laminates will be used to study the importance of higher order terms in the inplane
displacement field. A sinusoidally varying load is applied over the top surface of the lami-

nate:

1101

all

D01

27

. 1tx
t -_ — 1.27
12 Sin ( )

2.6.1 Eleven Layer Sandwich Panel

A simply supported eleven layer sandwich panel is considered with a length-to-thick-
ness ratio of four. The material properties and layup scheme are given in Table 2.1 and
Table 2.2, respectively. The material properties and layup sequence of the face sheet have
been studied previously (Murakarni, 1986), and provide a good test for the above models
because the face sheet consists of layers of distinctly different stiffness properties. Mate-
rial 1 is compliant in tension/compression and shear. Material 2 is stiff in tension/compres-
sion and shear. Material 3 is stiff in tension/compression and compliant in Shear.

Figure 2.2(a) shows the variation of inplane displacement field through the thickness
at the simply supported edge for the three forms of higher order theories using a constant
normal strain assumption. The elasticity solution and FSDT solution are also plotted for
the sake of comparison. It can be seen that the HZZT3 theory predicts results very close to
the elasticity solution. As the laminate is symmetric, the quadratic term in the inplane dis-
placement field vanishes, so the results from the HZZT3 are expected to be the same as the
HZZT3 results. The difference in results between HZZT3 and HZZT3 can be attributed to
the approximation made in the shear strain field for the HZZT3 theory. Figure 2.2(b)
shows the inplane displacement field for a constant normal stress assumption. Table 2.3
Shows the normalized center deflections of the beam for all the theories with constant nor-
mal strain and constant normal stress assumptions. The transverse displacements are nor-

malized with respect to the exact elasticity solution. It can be seen that the constant normal

28

stress assumption yields results closer to the exact solution in comparison to the constant
transverse normal strain approach. The constant normal stress assumption has improved
the transverse displacement at the top of the laminate by 4%. It can also be seen that the
difference in results between quadratic and partly cubic theories is not significant. FSDT
with Whitney’s shear correction factor overestimates the transverse displacement by 6.8%.
Figures 2.3(a) and 2.3(b) Show the through-thickness variation of inplane stresses in the
face Sheet at the center of the beam. Figures 2.4(a) and 2.4(b) show the variation of shear
stresses at the simply supported edge with constant strain and constant stress assumptions.
FSDT results are not shown as they are in large error compared to the exact solution. The
through-thickness variation of transverse normal stress with constant normal strain and
with constant normal stress assumption are shown in Figures 2.5(a) and 2.5(b), respec-
tively. Figures 2.6(a) and 2.6(b) Show the normal strain distribution through the thickness
with constant normal strain and constant normal stress assumptions, respectively. The con-
stant normal stress model predicts the normal strain in the core to be the average of the
exact solution.

From a comparison of the results with constant nOimal strain and constant normal
stress assumption, we can see that the constant normal stress assumption improves the pre-
diction of displacements, strains and stresses. The transverse normal strains are especially
improved, providing better understanding of mode-l delamination. For brevity, only the

results with normal stress assumption are presented for the subsequent case.

29

2.6.2 Composite Armored Vehicle (CAV) panel

In this example, a panel from the US Army composite armored vehicle (CAV) is stud-
ied. The laminate is divided into three sections and it is composed of fifty-five layers. The
details about the layup sequences are given below.

1. Inner Shell: Four plies Of S-2 Glass/Phenolic Fabric with stacking sequence [902/02]
with thickness of 0.01 in. for each ply. Thirty seven layers of SZ-Glass/8553-40 epoxy
tows with stacking sequence {[0/90][45/45/0/90]4[45/0/45],[90/0l-45/45]4} with

thickness of 0.021 in. for each ply.

2. Armor core: One layer EPDM rubber with thickness of 0.0625 in. and one layer of
ceramic tile of thickness 0.7 in. are used. This region is used to absorb impact energy.

3. Outer shell: Twelve plies of 8-2 Glass/855340 epoxy fabric with stacking sequence [0/
90/45l-45l0/90]s with thickness of each ply 0.01 in.

To facilitate comparison of the present models with exact solutions, the layers oriented
at 450 are changed to 900 and the layers oriented at -450 are changed to 00. The material
properties and laminate sequences are given in Table 2.4 and Table 2.5, respectively. The
length-to thickness ratio is again chosen ”to be four. Figure 2.7 and 2.8 Show the variation
of inplane displacement and inplane stress through the thickness of the laminate. The nor-
malized transverse displacements at the center of the beam for various models are given in
Table 2.6. In this case the results obtained from HZZT2 are closer to the exact solution
than HZZT3 results. This shows that for highly unsymmetrical laminates the quadratic
term is dominant compared to the cubic term. The transverse displacement predicted by

the fully cubic model is in error with the exact solution by 12%. The quadratic model and

p.117:
resp
ill:
the)
hill

mm:

Side:
std

ills

2.7

30

partly cubic model predict deflections with an error of approximately 15% and 16%
respectively. For this case FSDT predicts deflections which are 215% in error with the
exact solution. The transverse shear Stress in the laminate at the simply supported edge is
shown in Figure 2.9. In this case, the predictions of Shear stresses by all three models are
rather poor, though the through-thickness average of the shear stress is adequately approx-
imated. This suggests that the current zig-zag models are not sufficient for predicting
transverse shear stresses in such a complicated laminate, and motivates the usage of a sub-
larninate approach to capture the shear stress variation [Averill and Yip (1996), Cho and
Averill (1997)]. Figure 2.10 shows the variation of transverse normal stress through the
thickness of the laminate at the center ofthe beam. The transverse normal strain through
the thickness at the same location is shown in Figure 2.1 1. The transverse normal strain
predicted by the current constant transverse normal stress model is very reasonable con-
sidering the simplicity of the model. The through-thickness variation of transverse normal

strain is predicted accurately, even though the value of the strain predicted in the rubber

layer is in error by nearly 100%

2.7 Preliminary Conclusions

Three different higher order theories (cubic, quadratic and partly cubic) including zig-
zag effects and transverse shear continuity are presented and studied. The theories are
developed in a such way that degrees of freedom associated with the problem will involve
engineering degrees of freedom (displacement and rotations) at the top and bottom of the
laminate. This type of formulation is advantageous as it provides an opportunity refine the

model in the thickness direction using multiple sublaminates. This advantage can be con-

31

veniently exploited in finite element analyses. Unlike the full cubic theory, the quadratic

and partly cubic theories are C0 and lead to convenient finite element formulation. Based

on the numerical results, the following conclusions can be drawn.

1. All higher order theories proposed herein predict the displacements and local stresses
very accurately compared to FSDT.

2. For symmetric laminates, as the quadratic term in the inplane displacement field van-
ishes, the partly cubic theory predicts better results than quadratic theory. But for
highly unsymmetric laminates, it is observed that the quadratic theory predicts better
results than the partly cubic theory. A quadratic zig-zag theory is preferred over a partly
cubic zig-zag theory for analysis of structures with general lamination schemes. Also,
as the quadratic theory is complete up to quadratic term in the inplane displacement
field, it is expected to converge more readily to the elasticity solution when refinement
using multiple sublaminates in the thickness direction is performed.

3. The assumption of a constant transverse normal stress leads to improved prediction
capabilities, especially in the cases where the material properties vary drastically

through the thickness.

32

Table 2.1Materia1 Properties of Eleven Layer Sandwich Panel

 

 

 

 

 

 

I Property Material-1 Material-2 Material-3 Core
Ex 1 .0e06 33 .EO6 25.0e06 5.0e04
El 1.0e06 l .0e06 l .0e06 5.0e04
ze 0.2e06 8.0e06 0.5e06 2.17e04
v n 0.25 0.25 0.25 0.15

 

 

 

 

 

 

 

Table 2.2 Lay-Up Sequence of Eleven Layer Sandwich Panel

 

 

 

 

 

 

 

 

 

 

 

 

 

 

N11121:; Material I'D “1:211:88
1 I 1 0.010
2 2 0.025
3 3 0.015
4 1 0.020
5 3 0.030
6 Core 0.800
7 3 0.030
8 1 0.020
9 3 0.015
10 2 0.025
1 l 1 0.010

 

 

 

 

 

 

 

33

Table 2.3 Normalized Center Deflections of Simply Supported Beam under
Sinusoidal Load (eleven layer sandwich panel)

Table 2.4 Material Properties of TACOM Composite Armored Vehicle (CAV) Panel

 

 

 

 

 

 

 

 

 

 

 

 

Constant Normal Constant Normal
Strain Stress

Theory

Wb Wt Wb Wt
HZZT3 1.0291 0.9483 0.9859 0.9876
HZZT2 1.0270 0.9465 0.9839 0.9853
HZZT3 1.0279 0.9473 0.9842 0.9856
FSDT 1.0680

 

 

 

Material

 

SZ-Glass
Phenolic
Fabric

S-2 Glass
8553-40
Tows

S-2 Glass

8553-40
Fabric

 

 

EPDM
Rubber

Ceramic

   
 
 

 

 

 

 

 

 

 

 

 

Ex 3 .0E06 6.2E06 3 .0E06 3 .0E03 5 .0e06

El 1.2806 1.0e06 1.1e06 3.0e03 1.25e05
ze 4.6605 0.3e06 3.9e05 l .0e03 8.5e03
VXZ 0.18 0.29 0.18 0.45 0.15

 

 

 

34

Table 2.5 Lay-up Sequence of CAV Panel

 

 

 

 

 

 

 

I Region Material Type “2311:3658 J
S-2 Glass/Phenolic 0.01 per ply
(90/90/0/0)

Inner She" S-2 Glass/8553.40 Epoxy Tows 0.021 per ply
((0/90) (90/0/0/90)4 (90/0/90) (90/
O/O/9O)4)
EPDM Rubber 0.0625

Arm“ core Ceramic 0.70

Outer Shell S-2 Glass/855340 Epoxy Fabric 0.01 per ply
(0/90/0/90/0/0/90)2

 

 

 

 

 

 

Table 2.6 Normalized Predicted Center Deflections of CAV Panel by Various

 

 

 

 

 

 

 

Theories
Theory “’1; wt I
HZZT3 0.8924 0.8793
HZZT2 0.8940 0.8544
HZZT3 0.8750 0.8380
FSDT 3. 15 1

 

 

 

 

 

35

 

 

Figure 2.1. Multilayer beam geometry and coordinate system.

36

 

 

 

 

 

1.0 . e r , - r r
\\ .
\‘ ,g"
0.3 _ Elasticity \ ,2! -
-----~»~ HZZT3 \ x,-
- - - - HZZT2 \ ,r’f
--- HZZT3 ~\ ,1 /
- - - F DT . ,
0'6 I S \ I’// F
\ [I /
1 \x/
x):
’3’, \
0.4 ’4’.- \
I}, \
I}! ‘
If; \‘
0.2 ~ ,9" \ -
. / \
\
00 . 1 . 1 . ‘\ l .
4.06-06 -2.0e-06 0.09+00 2.06-06 4.06-06

Li

X

Figure 2.2(a). Through-thickness variation of inplane displacement at the end of
simply-supported beam (eleven layer sandwich panel, L/h=4)

0.l

z'n

SM

37

 

 

 

 

 

1.0 ‘ - ‘ ' T , . i ' i
\\‘
. . \ ,7
Elastrcrty ‘\ I/ _
0.8 - ........... HZZT3 \ 4,,
---- szrz . ,I,
-—- szra \ ,r/
—-— FSDT \‘ x/
0.6 i— \ I”/ 3;" ..
\ / ’-
\ ’ '
Z/h ‘ II// 1
[AK 3'...
,I / \
0.4 F I / ,3 -
II/ '2 \
I / ,." \
7’; 3"", \\\
0.2 '- I; ’3... \ -
\
F \
0.0 l 1 . l . \ 1 .
4.06-06 -2.0e-06 0.06-1-00 2.06-06 4.06-06

U

Figure 2.2(b).Through-thickness variation of inplane displacement at the end of
simply-supported beam(eleven layer sandwich panel,L/h=4,constant transverse normal stress)

0.1

1M)
0]

38

 

    
   
  

 

 

 

 

0.10 . r r . . r
. /
l
i
l / — Elasticity
. ------------ HZZT3
!____ ---- HZZT2
— — - HZZT 3
— - -— FSDT
2’“ 0.05 ~ .
I...
l
l
l
I I '
l
0.00 r l 4 711 a l .
-100.0 -50.0 0.0 50.0 100.0

0'

XX

Figure 2.3(a). Through-thickness variation of inplane normal stress at the center of simply-
supported beam(eleven layer sandwich panel, L/h=4)

39

 

   
     
 

 

 

 

 

 

0.10 . . s . -
I /
l
l
1
l — Elasticity ‘
'— ‘ — ‘ ----------- HZZT3
- - - - HZZT2
— —— - HZZT3
— - — FSDT
2’“ 0.05 ~ _
i—
1 l
l
1
/ i .
1
0.00 L 1 r 11 1 1
-100.0 -50.0 0.0 50.0 100.0

are:

Figure 2.3(b). Through-thickness variation of inplane normal stress at the center of simpl -
supported beam(eleven layer sandwich panel, Uh=4,constant transverse normal stress

40

 

 

 

 

 

 

Figure 2.4(a). Through-thickness variation of transverse shear stress at the end of simply-

1.0 r M ._ . , ..‘ ,
fl“ '
k 1
0.8 ~
— Elasticity
szra
--—- HZZT2
0.6 L —-'- HZZT 3
2m .
0.4 r
0.2 e
b
.'a
l
1
0.0 1 . lr—— r . 1 . 1 . 11 1
-1.0 -0.5 0.0 0.5 1.0 1.5
0'

supported beam(eleven layer sandwich panel, Uh=4)

 

41

 

 

 

 

 

 

 

 

1.0 . r - LL. . r . , \ . , .
“RD
0 1|:
'8 —— Elasticity ll
------------ HZZT 3
- - - - HZZT 2
— — — HZZT 3
0.6 — d
0.4 — ~
0.2 _ r
ii: i
0.0 i 1 . 1F— 1": 1 1 1 9 3‘; l 1
-1.0 -0.5 0.0 0.5 1.0 1.5 2.0
c

X!

Figure 2.4 b).Through-thickness variation of transverse shear stress at the end of simply-
supporte beam(eleven layer sandwich panel, Uh=4,constant transverse normal stress)

42

 

 

1.0 / z'/J
/'/—"/
0'8 ’ —— Elasticity *
------------ HZZT3
- - - — szrz
— - - HZZT3
0.6 - -

 

 

 

 

 

0.4 M _
0.2 _ M
/ //
0.0 . . . a/ . .
-1.0 0.0 1.0 2.0 3.0
0'

22

Figure 2.5(a). Through-thickness variation of transverse normal stress at the center
of simply-supported beam(eleven layer sandwich panel, Uh=4)

43

 

 

 

 

 

 

1.0 . l ' i p i ' i
Elasticity
0'8 “ ---------- HZZT3 “
---- HZZT2
——- HZZT3
0.6 ~ ‘ a
z/h
0.4 r -
0.2 - _
0.0 1 1 a 1 l J l r
0.0 0.2 0.4 0.6 0.8 1.0

0’

22

Figure 2.5(b). Throu h-thickness variation of transverse normal stress at the center of
simply-supported beam(e even layer sandwich panel, Uh=4, constant transverse normal stress)

44

 

 

   

 

 

 

 

1,0 . . . fi\ . . . , a . . ,
0.8 » .
— Elasticity
HZZT3
0.6 r szrz ~
——- HZZT3
2’“ M
0.4 r —
0.2 r u
0.0 4 m r 1; r r L 1 1 1 1 l L 1 a
-1.0e-05 0.0e+00 1.0e-05 2.0c-05

822

Figure 2.6(a). Through-thickness variation of transverse normal strain at the center of simply-
supported beam(eleven layer sandwich panel, Uh=4)

45

 

 

   
 

 

 

 

 

 

1_o - . . . e . , . . . ,
_ “‘1
0.8 — _
— Elasticity
------------ HZZT 3
- - - - HZZT 2
0'6 — — - HZZT3
z/h i
0.4 r —
0.2 ~ —
i l
0.0 1 r L 1 \ l 1 1 i r r r 1 r
-1.0e-05 0.0e+00 1 .0e-05 2.0e-05

8

22

Figure 2.6(b). Throu h-thickness variation of transverse normal strain at the center of simply-
suppcrted beam(e even layer sandwich panel, Uh=4, constant transverse normal stress )

46

 

   
   

 

 

 

   

 

1.0 . l r l r
— Elasticity
0.8 r HZZT 3 7
- - - - HZZT2
- — - HZZT 3
— - - FSDT
0.6 r _
0.4 c r
0.2 ~ a
0.0 r l . 1 .1 l r \_ . l
-2.0e-05 -1 .06-05 0.0e+00 1 .06-05 2.06-05 3.06-05

U

Figure 2.7. Throu h-thickness variation of inplane displacement at the end of simply-
suported beam(fi five layer CAV panel, Uh=4, constant transverse normal stress)

47

 

 
  
 
  

 

 

 

 

 

 

 

 

1.0 .
— Elasticity
........... l-lZZT3
0.8 ---— HZZT2 ‘
— - - HZZT3
- - - FSDT
0.6 a
J
I
0.4 _
0.2 -
0.0 r r 1 1
400.0 -50.0 0.0 50.0 100.0

Figure 2.8. Throu h-thickness variation of inplane normal stress at the center of simp -
supported beam fifity five layer CAV panel, Uh=4, constant transverse normal stress

48

 

1.0 . \ , . , . 1
\
\\|l
Ii
1 . .
0.8 — l 32:23” —
i - - - - HZZT2
E — —- HZZT3
0.6 . j -

 

0.4 ~

 

 

 

0.0 . . . . .
-1.0 0.0 1.0 2.0 3.0 4.0

 

Figure 2.9. Throu h-thickness variation of transverse shear stress at the end of simpl -
supported beam fifity five layer CAV panel, Uh=4, constant transverse normal stress

49

 

 

 

 

 

 

1.0 I F ‘l I I
—-- Elasticity
0.8 . ........... HZZT3 —
---- HZZT2
—-- HZZT3
0.6 ~ _
Z/h M
0.4 ~ ~
0.2 L -
0.0 r l 1 1 r 1 1 1 1
0.0 0.2 0.4 0.6 0.8 1.0

0'

22

Figure 2.10. Through-thickness variation of transverse normal stress at the center
of simply supported beam(fifty five layer CAV panel, Uh=4,constant transverse normal stress)

50

 

 

 

 

 

 

 

 

 

1.0 r r
— Elasticity
-------~ HZZT 3
0.8 — — - - - HZZT 2 ‘
— - - HZZT 3
0.6 - —
Z’h / /
0.4 _ r
0.2 — _
M
0.0 . r 1 1 1 L
-1 .0e-04 0.0e+00 1 .0e-04 209-04

8

H

. Figure 2.11 Through-thickness variation of transverse normal strain at the center of
Simply-supported beam(fifty five layer CAV panel, Uh=4, constant transverse normal stress)

CHAPTER 3
FINITE ELEMENT ANALYSIS OF LAMINATED
COMPOSITE BEAMS

3.1 Introduction

A new C0 finite element for accurate analysis of laminated composite beam structures
is developed in this chapter. The element formulation is based on a quadratic zig-zag
layer-wise theory developed previously in Chapter 1. The theory assumes a zig-zag distri-
bution of the inplane displacement field through the thickness and satisfies the interlami-
nar shear stress continuity across the layer interfaces. In developing the finite element
formulation, the shear strain fields are made field consistent and thus the shear locking
phenomenon is eliminated. A new transverse normal strain is derived by assuming the
transverse normal stress to be constant through the thickness of the laminate. This assump-
tion is shown to remove Poisson’s ratio stiffening. The results obtained from the present
finite element are found to be in good agreement with exact elasticity solutions available
for simply-supported beams. A multi sublarninate approach which is simple to implement
With the present element is Shown to improve the predictions of the present model for

Complex laminated structures.

51

52

3.2 Formulation of the Theory

3.2.1 Multilayer Beam Geometry

The geometry and coordinate system of a laminate consisting of N layers is shown in
Figure 3.1. The integer k denotes the layer number, with layer one. located at the bottom of
the laminate. The origin of the coordinate system is taken at the bottom of the laminate.
The length and thickness of the beam are assumed to be L and h, respectively. 2“ denotes
the distance from the reference surface to the bottom of the kth layer. Each lamina (ply) is

assumed to be orthotropic and is perfectly bonded to adjacent plies.

3.2.2 Displacement Field

The initial form of the displacement field for the quadratic zig-zag theory developed in

[Aitharaju and Averill (1997a)] is given below.

k-l
aid. 2) = ub(x) + zvbm + 22¢(X) + 2 (z — dis-(x)
i= 1 (3.1)

uf(x, z) = (1 — 19W”) + (and)

Where ukx, ukz are the displacements in x and 2 directions, respectively, in the layer k. ub,
‘15 are the inplane displacement and rotation at the bottom of the laminate. wb and wt are
the transverse displacements at the bottom and top of the laminate. 11> is a higher order
tfirm in the expansion of the inplane displacement field. The zig-zag functions £1. are

i“Cluded only in the inplane displacement field. The zig-zag functions £1. are determined

53

in terms of other variables by enforcing the continuity of transverse shear stresses at each

layer interface:

1.’ = 1' atz = zk (3.2)

In determining the shear strain, the contribution from the transverse displacement is

assumed to be constant in all layers. That is, the shear strain is approximated as:

k _ Bu: Bu: _ Bu: 1 dwb 3w, 33
sz‘a—z+5;"a7+§(a—x+§;) H
The zig-zag functions i,- are given by:
_ . ._ A 1 dwb 3w,
gr - ain+bz¢+ai§('-a-;+§;) (3.4)

where

 

. '_1

. G‘ ' .
WI 2..)

k=l

i i—l
~ G
b, = [Gi+l—1][2z,+ 2 8k]
k=l

Gi is the shear modulus of the ith layer. It can be seen from Equation (3.5) that the zig-zag

(3.5)

 

fllllCtious depend on the ratio of shear moduli of adjacent layers, and zig-zag effects will

be Small if this ratio is small. The inplane displacement field can now be written as:

54

k-l k-l
uxk(x,Z) -_- ub+wb[2+ 2 (Z‘Zi)&i]+¢[zz+ Elk-2,013,)

i=1 i=1

 

 

(3.6)

1 8WI; aw: k-l ,.

+ §( 8x + 8x) 2:1“ -Z,°)a,-

1 =
. . 1 aWb aw; . . .
In the above displacement field, the functions (1) and §(-a—x + 3;) are eliminated in
favor of the displacement and rotation at the top surface of the laminate:
N
du ]
N x
x Z = h [ dz Z = h
. . . . 1 awb aw: . .

Admittedly, the elimination of the term §(-a—x— + 3;) in favor of a new variable may

be viewed as inconsistent. However, as the form of the assumed displacement field is

Somewhat arbitrary, no variational crimes have been committed. Further, it will be shown
that this approach leads to a theory that is both computationally, efficient and accurate.

The kinematic displacement field can now be written in terms of

u b, at, \V b’ Wt, w b’ wt of the sublarninate using C0 functions. Here, as before, the sub-

Scripts b and t are used to denote bottom and top, respectively, of the sublarninate. These

degrees of freedom facilitate stacking of the sublaminates when refinement in the thick-

I“e88 direction is desired. The final form of the displacement fields can be written as:

a I"

my a7 a=1,27=l,2

(3.8)

_ k
WYKI7

Where the index a is used to denote the degree of freedom associated with inplane dis-

Placement field (u, w) and 'y is the index used to denote the bottom and top faces of the

55

sublarninate. For example, an is used to denote ub, 321 is used to denote W1; and W2 is
. . . . . . k k

used to denote wt. Summation on repeated indices IS implied. 1a.! and K I? are the shape

functions associated with the thickness direction. They are given in Equation (A2) of

Appendix A.

3.2.3 Strain-displacement and Constitutive Relations

The strains in the layer k can be calculated as

au" au" au" 3.."
k _ x k = x z k = 1
8H ' a'x' 7“ 32 + a”? 822 732 (3-9)

The stress in the lamina k is given by

{o}" = [Q]"{e}" (3.10)

The coefficients in matrix Q are given in Appendix B.
While evaluating the transverse shear strains in the energy expression, Equation (3.3)

is used. But during the post-processor, the transverse shear strains are calculated with the

1 awb aw: . . . . . .
rm i( + —) in Equation (3.3) replaced With other displacements and rotation van-

E’ 3x

ables using Equation (3.7). This procedure is observed to produce continuous transverse

shear stresses across the ply interfaces.

56

3.2.4 Approximation for the Normal Stress, 022

In the above displacement formulation, the transverse normal strain is assumed to be
constant through the thickness. For laminates wherein the layer properties vary drastically
or in the case of a sandwich panel with a thick, compliant core, the transverse normal
stress evaluated by the present model is not predicted satisfactorily. To improve the distri-
bution of normal stress and its related energy in the laminate, a modification is pr0posed
by assuming the transverse normal stress to be constant throughout the thickness of the

laminate. The value of assumed constant normal stress, 6 is calculated by using Reiss-

ZZ’

ner’s mixed variational principle [Reissner (1984)], which yields the following condition:

N ZK
k k
2 j (ea—523cc = (3.11)
k=l z,“I

In Equation (3.11), 8:2 is the strain derived from the displacement field and is given by

k = (Wt—W17)

822 h (3.12)

kc

an is the transverse strain in the layer k evaluated from the constitutive equation assum-

ing a constant variation of transverse normal stress, 6 and is given by

22’

k 5 Q12 k
8,: = é—Texx (3.13)
922 922

Substituting Equations (3.12) and (3.13) into Equation (3.11) and solving for 6 we get:

22."

_ 1 817
“u = §[w,—w,,+( 3:3)507] ‘3'”)
where
N Zr 1
R = z I —k-dZ
k: 121:4sz (315)
N ZK Qk .
12 I:
Say - Z l T 1,, d
k: 1 z‘Msz
using this fizz , 8:: in each layer can be calculated from Equation (3.13). The new normal
strain expression is given as:
-k 1 (Wt-wb) 1 Barry Say k k
822 = 7—7— + T[(-5?) (T—Q121a7)] (3.16)
Q22 922

From Equation (3.16), we can"see that E; is no longer constant through the laminate
thickness but varies depending on the inplane and transverse stiffnesses and the displace-
ment field. The terms owing to the inplane displacement field in the transverse normal
strain are significant when the laminate is composed of layers with different transverse
modulus. This normal strain expression (Equation (3.16)) is used in the strain energy
expression. The procedure mentioned above leads to improved transverse normal strain
energy and it is observed to produce consistent results for normal strains and average
transverse normal stress. From a computational point of view, the above approach of
assuming 0; to be constant through the thickness of the laminate is shown to eliminate

Poison’s ratio locking in the associated finite element model.

58

3.3 Finite Element Formulation
3.3.1 Finite Element Model

According to the principle of minimum potential energy, the potential energy is sta-

tionary at equilibrium, i.e.:
511 = 5(U+V) = O (3.17)

where U is the strain energy stored in the system and V is the potential energy of the exter-

nal forces. The strain energy stored in the system can be written as

121% {I f (cxxexx + ozzezz + O'xz'sz) dx dz (3.18)
= Q, lit 1

The potential energy of the applied loads in the transverse direction can be written as

V = ItBWB do 13= 1,2 (3.19)
(2

where tl , t2 are the distributed transverse loads on the bottom and top surfaces of the lam-
inate, respectively.
The topology of the finite element developed is shown in Figure 3.2. Within each ele-

ment, the displacements and rotations are approximated as follows:

59

‘2“ (3.20)
W7=2Nifi2w a=1,2 7:1,2

i= 1
where N,- are the linear Lagrangian shape functions. Substituting Equation (3.20) into the

expression of the first variation of total potential energy, I'I(L‘¢ Way) , and collecting the

(l‘Y’

coefficients of variations of 17a i and ii) we obtain the stiffness matrix and force vector

7 Yi’

(for details, see Reddy 1993). A two point integration rule (full integration) was used
along the beam length direction. All the stresses are evaluated from the constitutive equa-
tions.

The present element has been incorporated into the MARC commercial finite element
code through a user defined subroutine USELEM (MARC Analysis Research Corpora-
tion, 1994). Mesh generation, imposition of boundary conditions and loads is carried out
using the host program. The host program calls the element subroutine for each element in
turn and assembles the stiffness matrices. The host program optimizes the calling
sequence to the element to reduce the front width for frontal solution or the band width for
band solvers. The host program treats the present element as a 12 degree of freedom four

node element.

3.3.2 Consistent Shear Strain Fields

The finite element model incorporates the field consistency concept of Prathap ( 1994)

for shear strain fields. All the variables in the model i.e., ub, ut, \Vb, Wt, wb’ w are

t

approximated using linear Lagrangian shape functions. With this approximation, the trans-

6O

verse shear strain field in the element is not consistent, which leads to shear locking
behavior. To see it clearly, the transverse shear strain fields need to be converted to the nat-
ural (elemental) coordinate system. Prathap and Somashekar (1988) termed the transverse
shear strain fields in the natural coordinate system as covariant shear strains. Examining
the covariant shear strains in the kth layer, 722 can be written as:

k 1

k _ 1 _ .. I
Ygz - u§.z+§(wb.§+wt.§) — “av anzfi‘whwwli) (3.21)

a=1,2 7:1,2

where E is the local coordinate in the x- direction( —1 .<_ E S 1 ).

As the variables 17a and W7 are approximated using linear Lagrangian functions, we

7

can write the shear strain fields as:

 

 

k l—§ - k 1+ - k 1 _ _
'Yéz = ( )“avllay,z+( zgjua721a7,2+Z(W72-W71)
(3.22)

2
_ é - _ 1k 1 - _ 1k 1 _ _

‘ 5 (“(172 T “(171) 017,: + imam] + “0172) ay,z + 2(W72 " Wyl)
The coefficients of the % term in the above transverse shear strain field are clearly incon-
sistent since they contain inplane displacements and rotations alone and will lead to spuri-
ous constraints when thin beams are modeled. Based on the consistency concept, these
inconsistent terms are eliminated and only consistent shear strain fields are retained. The

derived consistent shear field can be written as

k 1 _ _ k 1
752 = 5(uayl + “a72)lay,z + 2(W72 - W71) (3-23)

61

This consistent shear strain field is similar to eliminating inconsistent terms through selec-
tive/reduced integration. But the consistency concept removes only relevant inconsistent
terms and therefore will not produce spurious zero energy modes. The shear expression

given in equation (3.23) is used in the finite element formulation.

3.3.3 Poisson’s Ratio Stiffening

The phenomenon of Poisson’s ratio stiffening is observed when low-order continuum-
type elements are used to model fiexure using only one layer of elements through the
depth of the beam (Prathap 1994). When thin beam flexure is modeled, the transverse nor-
mal stress in the beam should be vanishingly small, which leads to 822 = —vexx. As the
bending strain, an varies linearly through the thickness and reverses its sigh as it passes
through the neutral axis, cu must also vary linearly through the depth of the beam and
have opposite sign in the top and bottom half of the beam. If all is constant in 2, then the
above constraint can not be satisfied for bending behavior. This results an error in deflec-
tion prediction of the order of v2 (usually about 10%). In the present formulation, the
variation of an derived from the constant normal stress assumption is of the same order as

an , thus the Poisson’s ratio stiffening is eliminated.

3.4 Results and Discussion

The accuracy of the present zig-zag finite element is investigated by analyzing a sim-
ply supported beam subjected to a sinusoidally varying load and comparing the results

with an exact elasticity solution. A five layer unsymmetrical laminate and a panel from the

62

US. Army Composite Armored Vehicle (CAV) are used. Pagano’s elasticity solution
(Pagano, 1969) has been modified slightly to make it appropriate for beams as opposed to
plates in cylindrical bending, and is used to assess the accuracy of the present numerical
solutions. FSDT results with Whitney’s shear correction factor (Whitney, 1973) are also
presented for comparison. The finite element results are obtained with a mesh of twenty
elements along the length of beam. All the stresses are evaluated from the constitutive

equations at the Gauss points nearest the point of interest.

3.4.1 Five Layer Unsymmetrical Laminate

A simply supported five layer unsymmetrical laminated beam is considered with sinu-
soidal loading on the top surface of the laminate. The lay-up and material properties are
given in Table 3.1 and Table 3.2, respectively. The material properties and layup sequence
of the laminate have been studied previously by (Murakarni 1986), and provide a good test
for the above models because the laminate consists of layers of distinctly different stiff-
ness properties. Material 1 is compliant in tension/compression and shear. Material 2 is
stiff in tension/compression and shear. Material 3 is stiff in tension/compression and com-
pliant in shear.

Various sublarninate stacking schemes used in the present analysis are shown in Figure
3.3. In the case of a l-sublarninate model, one four node element is used through the entire
thickness of the laminate. The 2-sublaminate model has two elements thrOugh the thick-
ness. Three elements through the thickness are used in 3-sublaminate model. The interface
between the sublaminates is made to coincide with the layer interfaces, a scheme denoted

as LAM]. At the interface the layers are perfectly bonded, so the inplane and transverse

63

displacements are made continuous. But the rotation components are left free to allow the
zigzag form of inplane displacement. To achieve this a special assembly procedure is
employed.

Figure 3.4 shows the variation of normalized mid-span transverse displacement at the
top of the laminate with aspect ratio, Uh, for various thickness refinements. The trans-
verse displacements are normalized with respect to the exact elasticity solution. It can be
seen that the error in the solution of transverse displacement with one sublarninate through
the thickness is 6% for Uh=4. Two sublaminates through the thickness reduces the error to
1.2%. Three sublaminates through the thickness gives results very close to exact elasticity
solution with a difference of only 0.7%. At aspect ratio, Uh=1000, the present solution
matches with the exact solution, even for one layer of elements through the thickness,
showing no shear or Poisson’s ratio stiffening. Figure 3.5 shows the variation of inplane
displacement at the simply supported edge with various thickness refinements for the case
Uh=4. The elasticity and FSDT results are also plotted for the sake of comparison. It can
be seen that with two sublaminates, the present element predicts results very close to the
elasticity solution. The through-thickness variation of transverse displacement at the cen-
ter of the simply supported beam is shown in Figure 3.6. Figure 3.7 shows the variation of
inplane stress at the center of the simply supported beam.

The through-thickness variation of shear stresses at the end of simply supported beam
is shown in Figure 3.8 for Uh=4. It can be seen that the l-sublaminate model gives only
the average value of the shear stress. Two and three sublarninate models yield improved
shear stress predictions. As the present element does not have transverse shear stress as

degree of freedom at the nodes (See Averill and Yip, 1996), the shear stress evaluated at

64

the interfaces of sublaminates is discontinuous. Also the present element does not satisfy
exact zero shear traction at the top and bottom of the laminate. The transverse shear stress
results with five sublaminates through the thickness (one element per ply) matches with
the exact elasticity solution very well.

Figure 3.9 shows the variation of transverse normal stress at the center of the simply
supported beam for several thickness refinements. The transverse normal stress is constant
in a sublarninate (element), and accurately predicts the average stress within the sublarni-
nate. The transverse normal strain in the laminate at the center of the simply supported
beam is shown in Figure 3.10. The transverse normal strain results predicted by the
present model, even with one element through the thickness, agree with the overall trends
of the exact elasticity solution. From these results, it can be concluded that besides elimi-
nating Poisson’s ratio stiffening, the constant transverse normal stress assumption has
improved the transverse normal strain predictions considerably compared to a constant
value through the thickness. Also, except for the shear stress distributions, all other com-
ponents of stress and strain are predicted very well with only two sublaminates through
the thickness.

To avoid the special assembly procedure for multiple sublaminates, an alternative sub-
laminate stacking scheme (LAM2) is also tried. Figure 3.11 shows the alternative 2-sub—
laminate scheme. In this scheme, the sublarninate interfaces are kept at the middle of
layer-3. As the rotation components are continuous in a layer, the usual assembly proce-
dure is used for all the displacements and rotations. Figure 3.12 shows the through-thick-
ness variation of inplane displacement for 2-sublaminates with sublaminate scheme

LAM] and LAM2. The results obtained in Figure 3 for sublarninate stacking scheme,

65

LAM], are again plotted here. It can be seen that there is no difference in results for
inplane displacement from the two stacking schemes. The through thickness variation of
transverse displacement obtained from LAM] and LAM2 stacking schemes is shown in
Figure 3.13. Figure 3.14 shows the variation of transverse shear stress obtained from two
different laminate stacking schemes. From these comparisons, it can be concluded that by
keeping the sublarninate interfaces inside the layer (LAM2 stacking scheme), no special
assembly procedure is needed and the results obtained are closer to the elasticity solution
than those of the LAM] scheme (for the discretization chosen here). However, when it is
required to capture drastic variations between plies, the LAM] stacking scheme is pre-
ferred over the LAM2 stacking scheme.

The above results show that the present element predicts the displacements and
stresses very close to the exact solution at an extreme aspect ratio of 4, whereas FSDT
even with appropriate shear correction factor does not predicts the results with reasonable
accuracy. In order to show the comparison of results between the present element and
FSDT at less severe aspect ratios, through-thickness distribution of inplane displacement
and inplane stress for an aspect ratio of 10 are presented in Figure 3.15 and 3.16, respec-
tively. Even at this aspect ratio the distribution of inplane stress predicted by FSDT is in
error with exact solution. The present element with l-sublaminate predicts the maximum
stress with an error of 0.9% compared to exact solution, whereas FSDT predicts the stress

with an error of 16%.

66

3.4.2 Composite Armored Vehicle (CAV) Panel

In this example, a panel from the US Army composite armored vehicle (CAV) is studied.

The laminate is divided into three sections and it is composed of fifty-five layers. The

details about the layup sequences are given below.

1. Inner Shell: Four plies of S—2 Glass/Phenolic Fabric with stacking sequence [902/02]
with thickness of 0.0] in. for each ply. Thirty seven layers of SZ-Glass/8553-4O epoxy
tows with stacking sequence {[0/90][45/45/0/90]4[45/0/45],[90/0l-45/45]4} with

thickness of 0.021 in. for each ply.

2. Armor core: One layer of EPDM rubber with thickness of 0.0625 in. and one layer of

ceramic tile of thickness 0.7 in. are used. This region is used to absorb impact energy.

3. Outer shell: Twelve plies of S—2 Glass/855340 epoxy fabric with stacking sequence [0/

90/45l-45/0/90]s with thickness of each ply 0.0] in.

To facilitate comparison of the present models with exact solutions, the layers oriented at
450 are changed to 900 and the layers oriented at -450 are changed to 00. The material
properties and laminate sequences are given in Table 3.3 and Table 3.4 respectively. The
length to thickness ratio is chosen as four.

In this problem, as the thin rubber layer influences the results considerably, the rubber
layer is discretely modeled using one element for the cases when multiple sublaminates
are used (LAMl stacking scheme). Figure 3.17 shows the sublaminate schemes used for
the problem. As mentioned before, a special assembly scheme is used to assemble the dis-

placements and rotations at the interface of sublaminates.

67

Figure 3.18 shows the through-thickness variation of inplane displacement with vari-
ous thickness refinements. One element through the thickness is able to capture the behav-
ior of inplane displacement very well. A 3-sublaminate model predicts results very close
to the exact elasticity solution. The through-thickness variation of transverse displacement
at the center of the beam is shown in Figure 3.19. FSDT theory is unable to provide even
reasonable deflection results for this case, even when an appropriate shear correction fac-
tor is used. The inplane normal stress variation at the center of the simply supported beam
is shown in Figure 3.20. The through-thickness variation of transverse shear stress at the
end of the simply supported beam is shown in Figure 3.21. A 7-sublaminate model is
required to get acceptable results for shear stresses. Figure 3.22 shows the variation of
transverse normal stress at the center of the simply supported beam. The through-thick-
ness variation of normal strain is shown in Figure 3.23. Three sublaminates through the
thickness is able to capture the normal strain variation very well. Again, It can be seen
that, except for transverse shear stress distribution, all components of stress and strain are

obtained very accurately with only. a few sublaminates through the thickness.

3.5 Preliminary Conclusions

An efficient C0 finite element was developed based on a quadratic zig-zag theory. The
element formulation includes zig-zag effects in the inplane displacement field and trans-
verse shear continuity at the ply interfaces. The element satisfies all the requirements of a
computationally efficient finite element model for analysis of multilayer laminates,
namely 1) the number of degrees of freedom is independent of the number of layers; 2) the

degrees of freedom involved are only the displacements and rotations; 3) no shear correc-

68

tion or penalty factors are present in the formulation. The element is developed in such a
way that it can be stacked in the thickness direction conveniently, when the refinement in
the thickness direction is desired. In the finite element formulation, all the strain fields are
made field consistent and thus shear locking is eliminated. A new transverse normal strain
field is derived by assuming the transverse normal stress to be constant through the thick-
ness of the laminate. This approximation is shown to remove Poisson’s ratio stiffening.
Also it is shown that constant normal stress assumption improves the predictions of nor-
mal strain. Overall, the element is both accurate and efficient, and provides excellent capa-

bility for analyzing thick and/or complex laminated structures.

69

Table 3.] Material Properties for Five Layer Unsymmetrical Laminate

 

 

 

 

 

 

 

Property Material- 1 Material-2 Material-3
Ex 1.0e06 33.0e06 ' 25.0e06
Ez 1.0e06 21 .0e06 1 .0e06

ze 0.2e06 8.0e06 0.5e06
v H 0.25 0.25 0.25

 

 

 

 

Table 3.2 Lay-Up Sequence for Five Layer Unsymmetrical Laminate

 

 

 

 

 

 

 

 

N11112:: Material I.D “£3588
1 1 0.10
2 2 0.25
3 3 0.15
4 1 0.20
5 3 . 0.30

 

 

 

 

70

Table 3.3 Material Properties of TACOM Composite Armored Vehicle (CAV) Panel

 

 

 

 

 

 

 

 

 

 

 

SZ-Glass S-2 Glass S-2 Glass EPDM
Material Phenolic 8553-40 8553.40 Ceramic
. . Rubber
Fabr1c Tows Fabr1c
Ex 3.0E06 6.2EO6 3.0E06 3.0EO3 5.0e06
1 .2e06 l .0e06 1 . 1e06 3.0e03 1.25e05
ze 4.6e05 0.3e06 3.9e05 1 .OeO3 8.5e03
VXZ 0.18 0.29 0.18 0.45 0.15

 

 

 

 

 

 

Table 3.4 Lay-up Sequence of CAV Panel

 

 

 

 

 

 

 

 

Region Material Type Thickness
(1n.)

Inner Shell S-2 Glass/Phenolic (90/ 0.01 per ply
90/0/0)
S-2 Glass/855340 Epoxy 0.021 per
Tows ((0/90) (9O/0/0/9O)4 ply
(90/0/90) (90/0/0/9O)4)

Armor Core EPDM Rubber 0.0625
Ceramic 0.70

Outer Shell S-2 Glass/855340 0.01 per ply
Expoxy Fabric

 

(O/90/0/9O/0/0/90)2

 

 

 

 

71

 

 

Figure 3.1. Multilayer beam geometry and coordinate System

72

C C
(“It’wlta‘l’lo (“Zt’WZv‘I’ZJ

 

_———-——__—__——d

P——————_———_-—-

P——————_——————-

 

 

 

C
(111b,W1b,\V1b)e (“2b1w2b,‘|’2b)

Figure 3.2. Topology of finite element

73

 

 

 

 

(“van’WlOl (“20W2W21)1
a
Layer 5
Layer 4
, ___________ Layer 3
Layer 2
g ---------- ‘. Layer 1
(“ltvwrtr‘VnQl (“Zb’WZb’W2b)l

1- Sublarninate scheme

 

2 :
(“ivwlv‘llrd 0121""20‘1’20
C C

-———-—-—-—_d

 

(“1b9w1b9‘l’lb)2 (“zb’w2bfil'2b32
finbd’u)’. luszL’th‘PnY
l

t ---------- 3

1 1
(“lb’wlb’W1b) (“2b9w2b’W2b)

 

 

 

at the sublarninate interface

1 2 1 2
W1t¢wlbandw2t¢w2b

2- Sublarninate scheme

Figure 3.3. Sublarninate stacking schemes used for
five layer unsymmetrical larninate (LAMl)

74

 

 

 

 

 

 

1.20 . - . 1' ' . * .
-- 1-sublaminate
1.15 - -—- 2-sublaminate ‘
\ — 3-sublaminate 4
\\ — ' _ FSDT
1.10 — \\ 1
\\
\
w/Wm1 .05 . ‘\ ‘
‘\
\
1.00 — 571152.13,“ _____________________ _
0.95 — .
0.90 . 1 . 1 1 . . 1 . ‘ ‘
0.0 10 100 1000
L/h

Figure 3.4. Normalized transverse deflections at the center of simply-supported
beam (five layer unsymmetrical laminate)

75

 

 
  
 
    
     

 

   

 

 

1.0 I .
-. — Elasticity
------------ 1 -sublaminate
0.8 - ~\ ---— 2-sublaminates
.\ — - - 3-sublaminates
— — -— FSDT
0.6 -
z/h .
0.4 ~
0.2 r
1
0.0 ‘ ‘ ‘
-0.0000010 0.0000000 0.0000010

U

Figure 9.5. Through-thickness variation of inplane displacement at the end of
Simply-supported beam (five layer unsymmetrical laminate, Uh=4)

1.0

0.8

0.6

0.4

0.2

0.0
5.5

76

 

l

l

  

— Elasticity

------ 1 -sublaminate
- - - - 2-sublaminates
— — - 3-sublaminates
— - — FSDT

l

4

l

.—-_-—--—-—-_---—-—-—-—--—-_--_
l A l I

 

 

 

06-06 6.006-06 6.506-06 7.006-06 7.506-06
U

2

Figure 3.5. Through-thickness variation of transverse displacement at the center of
Simply-supported beam (five layer unsymmetrical laminate, Uh=4)

77

 

 

1.0 . T ’ '
—-— Elasticity
---------- - 1-sublaminate
0.8 _ ---- 2-sublaminates _
——- 3-sublaminates
L —-— FSDT
0.6 F 4

 

 

0.2 F

 

 

.i

0.0 . l . l I I] . l L

-30.0 -20.0 -10.0 0.0 10.0 20.0
0'

X!

 

 

 

 

Figure 3.7. Through-thickness variation of inplane normal stress at the center of
simply-supported beam (five layer unsymmetrical larninate, Uh=4)

78

 

   

 

 

 

 

 

1.0 ' l ' \T r v 1 l T j 1
l
L ! —-— Elasticity
! ------------ 1-sublaminate
0-3 ” . | . ---—- 3-sublaminate ‘
' ' ——- 5-sublaminate
l
. -- h —-— FSDT
0.6 r -
0.4 ~ -
E 1
l
0.2 '- ! _(
I
_____________________ .3
0.0 r L I L I r _I A
-2.0 -1.0 3.0 4.0 5.0 6.0

 

Figure 3.8. Through-thickness variation of transverse shear stress at the end of
simply-supported beam (five layer unsymmetrical laminate, Uh=4)

79

 

 

   

 

 

 

 

 

 

 

1.0 fl 1
i
— Elasticity E
------------ 1-sublaminate I
0.8 - - - - - 2-sublaminates : -
— - - 3-sublaminates I
— - — 5-sublaminates 4. —J
i
0.6 — , _
l
I
r ' — ' - ' - 7 I
i 1 l
0.4 ~ ' l -
._ l _L.__.__'
0.2 ‘
0.0 ‘ J
0.0 0.5 1.0
022

Figure 3.9. Through-thickness variation of transverse normal stress at the center
of simply-supported beam (five layer unsymmetrical laminate, Uh=4)

80

 

1.0 . ' ' ' '
— Elasticity
0 8 ----------- 1-sublaminate “
- ’ - — - - 2-sublaminates
— — - 3-sublaminates
— - — 5-sublaminates
0.6 ~ A

 

 

 

 

 

 

 

0.4 *- r
0.2 r _
0.0 1 . . 1 1 . . I \ 1 . J "a L a
-5.08-07 0.06+00 5.06-07

8

21

Figure 3.10. Through-thickness variation of transverse normal strain at the center
of simply-supported beam (five layer unsymmetrical laminate, Uh=4)

81

 

 

 

 

 

 

 

 

(“It’wlv‘l’loz (02¢.W2t.\ll2()2
Layer 5
t ________________
Layer 4 2
(“1b1w1b9‘l’1b)2 (“21:“szsz
Layer 3
____(“iv‘i"_1t;‘.l’1_t)l_ _ _ _ _ 1119121111213: -
Layer 2 1
Layer 1 ———————————————— i
. 1 1
(ulb’wle’lb) (“zb’Wzle'zfl

at the sublarninate interface ‘1’: t = \fiband W2t = ng

Figure 3.11. Alternate sublarninate stacking scheme used for five
layer unsymmetrical laminate (LAM2)

82

 

  

 

 

   

 

1.0 I '
0. 8 _ -— Elasticity -
---------- 2-sublamiantes (LAM1)
— — - 2-sublaminates (LAM2)
0.6 ~ I
z/h
0.4 ~ -
0.2 L -
0.0 ‘ ‘ I 4
.0.0000010 0.0000000 0.0000010

Figure 3.12. Through-thickness variation of inplane displacement at the end of
simply-supported beam with various stacking schemes (five layer unsymmetrical laminate, Uh=4)

83

 

 

 

 

 

1.0 . 1
0.8 — _
0.6 - -
2’“ / — Elasticity 1
----------- 2-sublaminates (LAM1)
0,4 - .. —-- 2-sublaminates (LAM2) _
/
.3: I
0.2 ~ ,l 4
/
/
/
:1" /
0.0 :' I 1 . 1 1 1 1
0.00000600 0.00000620 0.00000640 0.00000660 0.00000680

U

2

Figure 3.13. Through-thickness variation of transverse displacement at the center
of simply-supported beam with various stacking schemes (five layer unsymmetrical Iamiante, Uh=4)

84

 

 

 

 

1 .0 ‘ “._. ' T ' I
\ .. — Elasticity
\\ \‘ .......--... 2-sublamiantes (LAM1 scheme)
\ '\“ — — - 2-sublaminates (LAM2 scheme)
0'8 b \\."'-... 2
\
\
\
\
\
\‘-.
0.6 * \'-._ _
\
z/h \
K..—
0.4 F ‘
’1....,,_.::,,,
02 ’- ///".’ ........ ..
. r”.
l
0.0 1 1 1 1 .
0.0 1.0 2.0 3.0
on

Figure 3.14. Through-thickness variation of transverse shear stress at the end of
simply-supported beam with various stacking schemes (five layer unsymmetrical laminate, Uh=4)

85

 

 
  
   
    

 

 

 

1.0 1 1

—- Elasticity
0.3 - ------------ 1-sublaminate .

—--- 2-sublaminate

——- 3-sublaminate

—-- FSDT
0.6 ~ —
0.4 ~ ~
0.2 - 1
0.0 1 1 1 I 1 \ l 1
-0.000020 -0.000010 0.000000 0.000010 0.000020

U

X

Figure 3.15 Through-thickness variation of inplane displacement at the end of
simply-supported beam (five layer unsymmetrical laminate, Uh=10)

86

 

1.0

0.6 F

 

 

—— Elasticity

------------ 1 -sublaminate
- - - — 2-sublaminate
_. _ — 3-sublaminate

 

 

 

 

 

 

0.4 — —-— FSDT -
0.2 — -
0.0 l 1 l 1 l P

-100.0 -50.0 0.0 50.0

0'

Figure 3.16 Through-thickness variation of inplane stress at the center of
simply-supported beam (five layer unsymmetrical laminate, Uh=10)

100.0

87

 

 

 

 

 

 

 

 

 

 

 

 

 

Outer Shell Outer Shel]
— — — - _ — — — Ee-f-armc— Cerarmc
. _ TEES”. _ _ :REbb‘ng: Rubber
Inner Shell Inner Shell
1- Sublarninate 3- Sublarninate

Figure 3.17. Sublarninate stacking schemes used for CAV panel

88

 

1.0 . r

— Elasticity
............ 1 -sublaminates
——- 3-sublaminates
\ =. —-— FSDT

0.8 ~

 

 

 

 

] 1‘ 1 1 1 l

 

 

0.0 * ‘ *
-2.0e-05 -1 .06-05 0.06+00 1 .0e-05 2.06-05 3.06-05
U

X

Figure 3.18. Through-thickness variation of inplane displacement at the end of
simply supported beam (CAV panel, Uh=4)

89

 

1 e o ' l .‘i a l

0.8 ~

0.6 -

 

 

 

Y — Elasticity I

 

 

 

 

1’
0.4 — :
II ------------ 1-sublaminate
l — — - 3-sublaminate
l
0.2 ~ l -
l
l
l
l
0.0 v' A 1 l 1 1
0.000110 0.000120 0.000130

U

Figure 3.19 Through-thickness variation of transverse displacement at the center
of simply-supported beam (CAV panel, Uh=4)

1.0

0.8

0.6

0.4

0.2

0.0
-100.0

90

 

 

 

 

 

‘1:—'—"7.

 

 

—— Elasticity
----------- 1 -sublaminates

— - - 3-sublaminates
— - - FSDT

 

 

 

 

 

50.0

100.0

Figure 3.20. Through-thickness variation of inplane normal stress at the center of
simply-supported beam (CAV panel, Uh=4)

91

 

 

 

 

 

 

1.0 t... 1 Y I l r
l
\ —— Elasticity
-------- ' 1 -sublaminate
0.8 - ——- 3-sublaminates “
--— 7-sublaminates
0.6 — / ‘
2’" 'wl— _________________________
-!.\~\‘ '
0.4 - ‘\-\ ..... |1 '
-, \
\\‘ l
\_ l
l 1
I l
0.2 — / ______ I 7
.. /’ I
—/.’/ |
‘ l
'1. 'z' I
0.0 4'” 1 1 ‘ i ‘
0.0 2.0 4.0 6.0 8.0
a

Figure 3.21. Through-thickness variation of transverse shear stress at the end of
simply-supported beam (CAV panel, Uh=4)

92

 

1.0

— Elasticity
0.8 ~ ............ 1-sublaminate
— —- 3-sublaminates

———————J

0.4 r

0.2 L

 

 

 

 

0.0 4 ‘
0.0 0.5 1.0

0'

Figure 3.22. Through-thickness variation of transverse normal stress at the center
of simply-supported beam (CAV panel, Uh=4)

1.0 1 ‘ '
—— Elasticity
------------ 1 -sublaminate
0.8 ~ ——- 3-sublaminates ‘
0.6 ~ ‘
l
0.4 ~ "
0.2 ~ ‘
1
0.0 1 1 1 l 1 1 4
-0.00010 0.00000 0.00010 0.00020

93

 

 

 

 

 

 

 

 

 

Eu

Figure 3.23. Through-thickness variation of transverse normal strain at the center
of simply-supported beam (CAV panel, Uh=4)

CHAPTER 4
ANALYSIS OF LAMINATED COMPOSITES
WITH LOCAL DAMAGE

4.1 Introduction

Unlike metals, composite materials suffer permanent loss of structural integrity under
the influence of external loads, giving rise to discontinuties in the form of micro cracks.
This process of permanent loss of integrity is called ‘damage’. There are three types of
damage that occur in composites, namely, (1) transverse cracking (2) delamination and (3)
fibre breakage. The transverse cracking and delamination are the most observed failure
modes. To estimate the ultimate loads of composites, it is essential to predict the damage
initiation and evolution.

There are well-established methods for predicting the onset of damage or first ply fail-
ure in a composite laminate. It is more difficult to predict subsequent failures after the ini-
tial damage has occurred, since the stress analysis of a composite laminate with a large
number of small cracks becomes intractable. In progressive failure analysis the problem of
studying the progressive failure is tackled in a indirect way. The stiffness properties of
damaged layers in a laminate are reduced according to a particular stiffness reduction

method. This stiffness reduction method depends on a particular type of failure predicted

94

95

in the composite. Thus the problem of stress analysis of a damaged composite is converted
to the problem of stress analysis of a composite with discontinuous material properties in
the plane of the composite and as well as in the thickness direction.

The delaminations in a composite can be modeled by embedding a compliant layer
between the delaminated layers. This approach is called the compliant layer approach or
the embedded layer approach. Such an approach has been shown to be valid for modeling
delaminations, provided the delaminations remain closed (i.e. no mode -1 deformation;
see Lu and Liu 1992). In this approach the material properties of the compliant layer are
reduced to a small value compared to the original properties. Thus stress analysis of com-
posites with delaminations is also converted to a stress analysis of laminates with discon-
tinuous material properties both in the plane and thickness direction of the composite.

An accurate prediction of stress and deformation states in composite laminates con-
taining material discontinuties (drastic variation in material properties) both in the plane
and through the thickness direction is very important to analyze the composites with
delaminations or to investigate progressive failure. The distribution of displacements and
stresses near the damaged region, where the material properties are discontinuous, is very
complicated and the analytical models are unpractical and perhaps unable to model them
explicitly. The most suitable tool to analyze composites with damage is the finite element
method. The key issue involved in the analysis of composites with damage is that it con-
tains three-dimensional effects.

From the theoretical point of view, the prediction of displacements and stresses in
composites with damage can be managed using three-dimensional finite element analysis,

where each layer is modeled using conventional brick elements. But such an approach

96

quickly becomes computationally expensive for .structures having a large number of lay-
ers. Therefore, it is highly desirable to have a shell type model that can capture local vari-
ations of displacements and stresses accurately with one element through the thickness
and at the same time is reasonably accurate compared to a full three-dimensional model.

Structural finite elements available in commercial finite programs do not require a
large number of elements for analysis of composites with damage, but are largely inaccu-
rate due to the following two reasons; First, the commonly assumed linear variation of
inplane displacement is not adequate in the regions where damage is present. Second,
these elements often do not adequately represent the effects of transverse shear and normal
stress. Hence higher order effects in the inplane displacement and a good representation
for shear and normal stresses are a must for improved structural analysis of damaged com-
posites.

The aim of the present chapter is to develop a reasonably accurate and computationally

efficient approach for modeling laminates with damage. Recently the authors developed a

C0 4-node zig-zag beam element (Aitharaju and Averill, 1997b) for analysis of composite
beams based on the quadratic zig-zag theory. The results obtained from this element are
found to be in good agreement with reference solutions. The topology of the beam ele-
ment is shown in Figure 4.1. The nodes are located on the top and bottom surface of the
element. The degrees of freedom at these nodes facilitate stacking of sublaminates when
refinement in the thickness direction is desired. This type of sublaminate scheme was used
efficiently by (Averill and Yip 1996, Aitharaju and Averill 1997b) to obtain results for
complex laminate cases. In this present chapter, the analysis of composites with damage

will be attempted.

97

The displacement form in the present model contains material parameters (shear rigid-
ities) of layers in the laminate. This material dependent displacement form is due to the
satisfaction of transverse shear stress continuity at the laminate interfaces. When damage
is present in the laminate, we encounter a material discontinuity (discontinuity in shear
modulus) across damaged and undamaged element boundaries, leading to a displacement
discontinuity in the model. In displacement finite element models the displacement form
must at least be continuous and hence a new scheme is developed to make the displace-
ments continuous across the elements. This forces the zig-zag displacement field to vary in
an element.The scheme developed here is called a Variable Zig-Zag Kinematic Displace-
ment model (VZZKD). The results with the new scheme will be compared with exact elas-
ticity solutions wherever available or a reference solution obtained by stacking multiple
sublaminates. A brief review of the displacement field for quadratic zig-zag theory and

finite element formulation for VZZKD model is presented below.

4.2 Formulation of the VZZKD model

The initial form of the displacement field for a quadratic zig-zag theory as developed

in Chapter 3 is given below

k k k k
“(211 +“112'l'll’b13 +‘l’114

(4.1)

N??- 317:"

k k
Klwb-l-sz,

where uvub are the inplane displacements at the top and bottom surfaces of the beam. wb

and wt are the transverse displacements at the bottom and top surfaces, respectively. II) are

98

the shape functions associated with thickness direction and contain zig-zag functions.

These zig—zag functions depend upon shear moduli of layers in the laminate. Kp are the

through-thickness shape functions associated with the transverse displacement. The above
equation is a result of satisfying transverse shear stress continuity at the ply interfaces for
a quadratic zig-zag theory. For details see Chapter 3.

It can be seen from Equation (4.1) that the displacement form is depend on the mate-
rial properties (shear moduli) of the layers in the laminate. When damage is present, there
is a discontinuity of inplane displacement at interfaces between damaged and undamaged
regions (see Figure 4.2). To make displacements continuous at the interface, the material
properties at the interface are replaced with an average property (see Figure 4.3). Keeping
the computational efficiency in mind, the material properties are averaged by one of the

following two averaging schemes.

Scheme -1
2 2
k G]: +612
Gm,1 — T7 (4.2)
G1 +62
Scheme -2
k k
20 G
k 1 2
av2 = “-7—; (4.3)
01-1-02

where G" is the shear modulus of the layer k at the interface i (i=1, left and i=2, right).

Gavmk is the average shear modulus at the interface and is used in deriving the kinematic

99

displacement field. It can be seen that the two schemes yield undamaged shear properties

when damage is not present.
The average shear modulus obtained from the two schemes is used in obtaining the
kinematic displacement field given in Equation (4.1) at the damage interface. Elsewhere,

the through-thickness shape functions are obtained using the actual shear properties.

Hence the functions ka are now functions of the spacial coordinate, x.

4.3 Finite Element Formulation

Within each element, the displacements, rotations and the shape functions I: are

approximated as follows:

2
2N1. Wri 7:1,2

i=1

Q
-<
11

2 (4.4)
WY: 2N1. W71 7:1,2

i=1

2
k k
1p: 2N, 1p, p=l,4

i=1

where N,- are the linear Lagrangian shape functions.

The strain-displacement relations in a layer k can be written as

100

1 Bu: 011,, an, 3111,, aw”,

k
in - s; - X’rta’2+$’3+‘5;’4
01’; 01’; 01’; 31f,
+uba—x+ut§x-+Wb.a—x+wta_x

(4.5)

In determining the shear strains, the contribution from the transverse displacement is
assumed to be constant in all the layers. The shear strains are approximated as:
,, au" au" 31" 1

Ya = _a_zx+§x_l = a—zx+§(wb,x+wt,x)

k k k k 1
“bll, z "' “112.1 1‘ ‘l’bI3, z 1’ W114, z '1' 5W1), x + W1, 1:) (4'6)

k k k k 1
“1111+ uth + \VbJ3 + 01,14 + 50%,): + w“)

In the above displacement formulation, the transverse normal strain is assumed to be
constant. But it is observed that when material properties vary drastically through the
thickness in a laminate, the constant normal strain assumption leads to incorrect transverse
normal strain energies. To improve the transverse normal strain and its related energy in
the laminate, a modification in the transverse normal strain is proposed by assuming the
transverse normal stress to be constant through out the thickness of the laminate. The
assumed constant normal stress is obtained using Reissner’s mixed variation principle

(Reissner, 1984).

N ZK

k k
2 l (fizz-Ezbdz = 0 (4.7)
k=l 211-1

where 8; is the strain derived from displacement field and is given by

(4.8)

 

(if: is the transverse normal strain in the layer k as evaluated from the constitutive equa-

tions assuming a constant variation of transverse normal stress, 6'2z , and is given by

._ k
k 0 k
82: = .2 _ .128 (4.9)
k k xx
Q22 Q22

Substituting Equations (4.8) and (4.9) in Equation (4.7) and solving for 6a , we get

 

w,- Wb Bu,
azz=( )R+ (i):- :)[511N1+512N2]+( 73)[521N1+522N2]

h
a a
(13’)[S31N1+ 53,111,] +(— 3"’——‘)[s,,,1vl + 54,111,] +
4.10
S {ill/—l S 8N2 S 3N1 3N2 ( )
+“b[11"a_ch + 123x——':]+“[ 213x 1’ 22 37]
BN 3N EN EN
“WP 31 '37] +532 3 x2]+‘l’:[541_ 8x '—1+S42— 372]
where,
h l
S = J— dz
oC’éz
]
R = = .
S (411)
" C"
1 12
5,. — {1 1,..—- dz]
J S 0 JCI;2

102

Using the 622. in Equation (4.10), the transverse normal strain ékzz in the layer k can be

evaluated.

A 4-node fine element based on the VZZKD model is developed and implemented in
the MARC finite element program through User Subroutines. The present element will be
used to analyze laminated beams with damage. The results from one element through the
thickness will be compared with reference solutions. The reference solution will be
obtained either from the exact elasticity solution, whenever available, or by modeling with

multiple four node elements in each layer through the thickness of the laminate.

4.4 Results and discussion

To test the accuracy of the present zig-zag finite element in predicting the response of
composite laminates with damage, two example problems are considered and results from
the present model with one element through the thickness will be compared with reference
solutions. In the first example, a simply supported beam with lamination scheme (0/90/0)
subjected to sinusoidally varying load is considered. The results from the present element
with one element through the thickness and FSDT will be compared with exact solutions

for various thickness ratios of the laminate and various stiffness degradation coefficients
of the 900 layer. In the second example, a cantilever beam with (0/90/0) and a (0/90)4s
lamination scheme subjected to tip load will be considered. The off-axis plies are assumed
to be damaged near the fixed end over a specified distance. The displacements and stress

obtained from the present element with one element through the thickness will be com-

pared with reference solutions.

103

4.4.1 Simply Supported Beam (0/90/0) Subjected to Sinusoidally Varying Load

A simply supported composite beam with (0/90/0) lamination scheme subjected to
sinusoidal load is considered. The material properties of 0 and 90 layers are given in Table

4.1. The thickness of 0 and 90 layers are assumed to ho and hgo, respectively. The thick-

ness ratio 71 is defined as,
_ h
h = 1’ (4.12)

The value of 1.1 is varied from 1.0 to 0.01., keeping the total thickness of the laminate
1.0. The aspect ratio of the beam considered for the analysis is 20. Present results are
obtained with one element through the thickness of the laminate. To simulate damage and
delamination in the composite beam, the stiffness properties of the 900 layer are degraded

by a damage index value. The damage index, A. is defined as:

= 1:90 = ffl (413)
E9011 G901!

where Ego and 690 are undamaged inplane modulus and shear modulus, respectively. E90d
and 690d are the inplane and shear modulus of damaged layers. The damage ratio is varied
between 1.0 and 105.

Due to the symmetry, only one half of the beam is considered for analysis. TWenty ele-

ments are used along the length of the beam. For this case, the available exact solutions

104

(Pagano, 1969) will be used to compare the results obtained from present the element and
FSDT.

The normalized center deflections obtained from the present element and FSDT are
presented for various damage ratios in Table 4.2. FSDT results are obtained using Whit-

ney’s shear correction factor (Whitney, 1973). From the Table 4.2, it can be seen that the
present element is able to predict deflections accurately up to a damage ratio of 104. At the
damage ratio of 105 the present model over predicts the deflections by 21% at 1.1 = 1.0. As

the thickness of 900 layer is reduced compared to 00 layer, the center deflection results
predicted by the present model are in close agreement with exact solution. The error at
high damage ratios and significant damage layer thickness is due to the fact that the
present element is not able to represent the shear fields exactly in the damaged layer. Also
at these high damage ratios and significant damage layer thicknesses, the shear energies in
the damaged layer are significant compared to the total bending energy of the laminate.

FSDT results show a similar trend but the error in center deflection prediction is signif-
icant even at small damage ratios. For example, for 7. = 100 and 72 = 1.0 , the error is

112%. This error increases dramatically for higher values of 7t and hence the results are
not given in the table. The error is only 1.8% at 7. = 100 and 71 = 0.0] . At the damage

ratio of 104 and I; = 0.01 , FSDT results are still 600% in error with respect to the exact
solution. These results show that FSDT results are not reliable for damage ratio’s higher
than 100 and for moderate damage layer thicknesses.

The above example problem gives an idea about the performance of the present ele-

ment and the element available in commercial packages for analysis of composites with

 

105

ply damage and delamination. The analysis of composites with small damaged layer
thickness with degraded properties represent a delamination problem. The progressive
failure analysis of composites with ply-failure is represented by keeping the thickness of
damaged layer significant compared to undamaged layers. It can be seen clearly that the
present element predicts the response of composites with delamination very accurately. In

the case of ply-damage, the global results obtained with the present element are very accu-

rate except for damage ratios higher than 104 and when the damaged layer thickness is sig-

nificant.

4.4.2 Cantilever Beam Subjected to Point Load (Cross-Ply Laminate)

Ply-Failure
In this example, a cantilever beam with (0/90/0) lay-up subjected to a unit tip load is

considered (see Figure 4.4). The aspect ratio of the beam is taken as 20. The off axis plies
(900) are assumed to be damaged up to a length 1* from the clamped end. In the present
case 1* was taken as 0.25L. The material properties of 0 and 90 plies are the same as given

in Table 4.1. The stiffness of damaged plies are reduced by 103 times the original stiffness.

As there is no available elasticity solution for the above problem, the reference solu-
tions are obtained by using 3 four node elements as developed in Chapter 2 in each layer
of the laminate. Twenty elements along the length of the beam are used. Displacements
and stresses are obtained using one VZZKD element through the thickness of the laminate

and are compared with the reference solutions.

106

Two types of averaging schemes as given in Equations (4.2) and (4.3) are tried. The
averaging scheme 1 yields the average shear modulus close to undamaged value and aver-
aging scheme 2 yields the average shear modulus close to the damaged shear modulus.

Figure 4.5 shows the inplane displacement in the damaged region at a distance of 2h
from the interface. Reference solutions are also presented. It can be seen that the present
solution with one element through the thickness yields results very close to the reference
solution for the two averaging schemes. Figure 4.6 shows the variation of inplane dis-
placement through the thickness at the interface between damaged and undamaged regions
using two averaging schemes. Figure 4.7 shows the variation of inplane displacement
through the thickness at a distance of 2h from the interface in the undamaged region. In all
the cases the inplane displacement obtained by the present element with one element
through the thickness matches very close with reference solutions for both the averaging
schemes. As both the averaging schemes provide solution which are very close, the
remaining set of results will be presented using averaging scheme 1 only.

Figures 4.8 to 4.10 show the through-thickness variation of inplane normal stress in
the damaged region near the fixed end, at 2h from the damaged/undamaged interface, and
in the undamaged region at 2h from the interface, respectively. It can be seen that present
results with one element are very close to the reference solutions.

Figures 4.1] to 4.13 show the through-thickness variation of shear stress in the dam-
aged region near the fixed end, at 2h from the interface, and in the undamaged region,
respectively. It can be seen that the shear stress results from the present element attempt to
capture the average of the reference solution, but not very accurately. Also, the present ele-

ment does not satisfy zero shear stress at the top and bottom of the laminate.

107

From the above results it can be seen that the present VZZKD model with one element
through the thickness gives all displacements and stresses (except shear stresses) very
close to the reference solutions for composites with damage.

Similar sets of results are presented for a (0/90)4S laminate. Figures 4.14 to 4.16 show

the inplane displacements from the present element and reference solutions. The inplane
displacement obtained from the present model at different sections does not match well
with the reference solutions. The comparison of inplane stresses from the present model
and the reference solutions are given in Figures 4.17 to 4.19. Figures 4.20 to 4.22 show the
transverse shear stress comparison. The reason for inaccurate predictions of inplane dis-
placement and stresses for this laminate is because the through-thickness variation of
transverse shear stress in the damaged region is too complicated to be represented by the

present model. In this case, multiple elements through-the-thickness should be used.

Delamination
In the above example problem, we have seen that the results obtained for laminate

sequence (0/90)4s using VZZKD model for ply-failure are not very close to the reference
solution. This is due to the reason that the laminate sequence (0/90)4S is complicated to
predict the response. In the present example problem the same (0/9O)4s clamped beam is
considered for delamination study. A thin delamination of thickness 0.0] is introduced at

the center of the beam. The properties of delaminated layers are reduced by 103. The

Delamination is assumed to extend up to a distance of 0.25 L from the clamped end (See

Figure 4.23).

108

Figures 4.24 to 4.26 show the comparison of through-thickness variation of inplane dis-
placement in the delaminated region at different sections along the length of the beam
using present VZZKD model and the reference solutions. From the figures we can see that
the results from the present VZZKD model match very well with the reference solution.
Figures 4.27 to 4.29 show the inplane stress comparison. The inplane stress obtained by
the present model with one element through the thickness match very close with reference
solution.

The present delamination study shows that present VZZKD can be used to model the

laminates with delamination accurately.

4.5 Preliminary Conclusions

In the present chapter, a reasonably accurate and computationally efficient model for
analysis of composites with damage is developed. This model facilitates the analysis of
composites during progressive failure due to ply damage and with delaminations. Such an
analysis can also be carried out by modeling each layer with a three dimensional brick ele-
ment. But this approach is computationally very expensive when a large number of layers
is present. The aim of the present model is to analyze damaged composites with just one
element through the thickness of the laminate. Towards this a new model is developed and
it is referred to here as the VZZKD mode]. Two example problems are studied and from
the results the following conclusions can be drawn.

1. When the damaged layer thickness is small (1% of thickness), the present model yields
very good results even for large damage ratios (up to A = 105 ). Hence the present

model can be used to accurately analyze composites with local delaminations.

 

109

2. When the damaged layer thickness is significant compared to undamaged layers, the
results from the present model are in close agreement with reference solutions up to a
damage ratio of 104. At the damage ratio of 105 the present model may contain signifi-
cant error compared to the exact elasticity solution. Hence at damage ratios lower than
104 the present model with one element through the thickness can be used very reliably

for analysis of composites with ply-damage and simple lamination schemes.

3. The present results show that the accuracy of the present model depends upon the
placement of damaged layers in the laminate and the percentage thickness of each dam-
aged layer. In case of laminates with multiple damaged layers of significant thickness,
the local stresses and inplane displacements predicted, by the present model might be in

error.

4. The commonly available FSDT elements available in commercial packages may not be
reliable after a damage ratio of 20 when the damaged layer thicknesses are significant
compared to undamaged thickness. When the damaged layer thickness is small the
results can be reliable up to a damage ratio of 100. This study shows that one has to be
very cautious while using the FSDT elements available in commercial packages for

progressive damage or delamination analysis of composites.

110

Table 4.1: Material Properties of (0/90/0) Laminate

 

 

 

 

 

 

 

 

 

Property 00 900
EX 25606 1 .0606
E2 1 .0606 l .0606
ze 0.5606 0.2606
vxz 0.25 0.25

 

 

Table 4.2: N ondimensionalized Center Deflections of Simply Supported Beam
Subjected to a Sinusoidal Load with Various Damaged 90o Layer Thicknesses 71 and

Damage Ratios, 7t .

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

71 Model A
1 10 50 100 103. 104 105

Present 0.9914 0.9960 0.9971 0.9978 0.9961 0.9991 1.2120
1'0 FSDT 0.9996 1.0708 1.5294 2.1234

Present 0.9880 0.9937 0.9969 0.9971 0.9964 1.001 1.2211
0'5 FSDT 0.9995 1.0695 1.5752 2.2454

Present 0.9848 0.9896 0.9954' 0.9962 0.9954 0.9978 1.1198
0'2 FSDT 0.9990 1.0332 1.1006 1.8089

Present 0.9945 0.9872 0.9935 0.9949 0.9954 0.9958 1.0514
0'1 FSDT 0.9988 1.0130 1.1777 1.4395

Present 0.9821 0.9828 0.9850 0.9869 0.9946 0.9948 0.9958
0'01 FSDT 0.9986 0.9991 1.005 1.018 1.541

 

 

111

e C
(ulvwltt‘l’lt) (“ZttWZtt‘VZO

 

_-————————————.

-————————-————‘

 

 

 

e
(We‘ll/16111196 (“2b1W2b1‘l’2b)

Figure 4.1. Topology of finite element

 

 

112

Damaged Region Undamaged Region

 

 

 

 

 

GIl \ \ 621
G}. / (,2.
G15 A (325

 

 

Kinematic displacement field

Fig.4.2 Damage and Undamaged Regions of a Laminate

 

 

 

 

 

Damaged Region Undamaged Region
G11 Gll G21 Gzl
G 1 2 G 1 2 G22 G22
G G
/ A 2 2
614 G14 624 624
G 15 G 1 5 G25 G25

 

 

Fig.4.3 Present Modeling Approach for Damaged and Undamaged
Regions

 

113

21

 

 

 

 

 

 

 

4 a»
\ Damaged Layer

Figure 4.4. Tip-loaded cantilever beam (0/90/0) beam with damage in the
off-axis plies near the clamped end

114

 

   
   

 

 

 

 

 

1.0 l f f ' r
— Reference Solution
- - - - 1- Sublarninate (scheme 1)
— — - 1- Sublarninate (scheme 2)
0.8 r -
l
l
0.6 e -
z/h 1
0.4 — -
0.2 ~ _
0.0 1 L 1 1 1 l 1
-0.000020 -0.00001 0 0.000000 0.00001 0 0.000020

L]

X

Figure 4.5 Through-thickness variation of inplane displacement in the damaged
region of clampled beam (0/90/0 laminate)

1.0

0.8

0.6

115

 

 

f Y Y I T V I I I V V I

Reference Solution

 

i - - - - 1- Sublarninate (scheme 1)

— — - 1- Sublarninate (scheme 2)

 

 

 

0.4 - ~
0.2 1 ~
1 1
0.0 1 1 1 1 1 1 1 4 1 1 1 1 1
4.06-05 -2.06-05 0.06+00 2.06-05
U

4.0e-05

Figure 4.6 Through-thickness variation of inplane displacement at the interface
of damaged and undamaged regions of clampled beam (0/90/0 laminate)

116

 

 

 

 

 

1 .0 ' l V l ' l

— Reference Solution

- - - — 1- Sublaminate(scheme 1)
0 8 _ — — -— 1- Sublamiante(scheme 2) _
0.6 ~ _
0.4 ~ —
0.2 — 4

1

0.0 1 l 1 1 1 1 1
-0.000040 -0.000020 0.000000 0.000020 0.000040

U

X

Figure 4.7 Through-thickness variation of inplane dispalcement in the undamaged
region of clamped beam (0/90/0 laminate)

117

 

   
 

1.0 1 1

— Reference Solution

- - - - 1- sublarninate
0.8 ~ 7
0.6 e _

0.4 ~ 4

 

0.2 F

 

 

 

I

0.0 J 1 1 1 1
-200.0 -100.0 0.0 100.0 200.0

 

0'

XX

Figure 4.8 Through-thickness variation of inplane normal stress near the fixed
end in the damaged region of clamped beam (0/90/0 laminate)

1.0

0.8

0.6

0.4

0.2

118

 

L - - - - 1-sublaminate

 

— Reference solution

 

 

l 1 J

 

0.0 1 l 1 1
-200.0 -100.0 0.0 100.0 200.0

U

Figure 4.9 Through-thickness variation of inplane normal stress at the middle of
damaged region of clamped beam (0/90/0 laminate)

119

 

 

 

 

1.0 1 1

Reference Solution
0'8 ‘ - - - - 1- Sublarninate ‘
0.6 — ..

0.4 _ *

 

 

 

 

 

0.2 F a
0.0 1 l 1 J_ 1 1
-100.0 -50.0 0.0 50.0 100.0
0'

Figure 4.10 Through-thickness variation of inplane normal stress in the undamaged
region of clamped beam (0/90/0 laminate)

120

 

 

 

1.0

0.8

W'_-—T‘———‘I' —

 

0.6 F

0.4 -

Reference Solution ~
- - - - 1 - Sublarninate

 

 

0.2

——Y——'—T————V—

 

 

 

 

 

0.0
0.0

0.5 1.0 1.5 2.0
0

Figure 4.11 Through-thickness variation of transverse shear stress near the
fixed end of clamped beam (0/90/0 laminate)

121

 

 

 

 

 

 

 

 

1.0
0.8 F 7
J
Reference Solution
0.6 F ---- 1-Sublaminate ‘
0.4 F -
:1
1
l
l
0.2 _ : -
l
l
l
l
l
l
0.0 : ‘ 1 r ‘ ‘
-1.0 0.0 1.0 2.0 3.0

0'

X2

Figure 4.12 Through-thickness variation of transverse shear stress in the middle
of damaged region of clamped beam (0/90/0 larninate)

122

 

1.0

0.8 F

0.6 F

0.4 F

 

0.0

—— Reference Solution
- - - - 1- Sublarninate

_--—-—————----—-——-——---—-------q

l l

 

 

-1.0

0.0 1.0 2.0
0'

Figure 4.13 Through-thickness variation of transverse shear stress in the
undamaged region of clamped beam (0/90/0 laminate)

3.0

123

 

 

 

 

0.00004

 

1 .0 l ‘ \ l l
‘1 Reference Solution
0.8 F l - - - - 1- Sublarninate .
\
0.6 I -
1'
l
l
l
l
0.4 1 _
\ \ ,
I
0.2 I ..
I
,\
0.0 - 1 1 1 ‘ 1 1
0.00004 -0.00002 0.00000 0.00002
u

X

Figure 4.14 Through-thickness variation of inplane displacement in the damaged

region of clamped beam ( (0/90)48 laminate)

124

 

  

 

  
  

1.0 1 r r
\
\ Reference Solution
0-8 7 ‘ - - - - 1-sublaminate ‘
0.6 F F

 

 

 

0.4 F F

0.2 F _

O 0 I l \\

-0.00004 -0.00002 0.00000 0.00002 0.00004
. ' u

Figure 4.15 Through-thickness variation of inplane dis Iacement at the interface
of damaged and undamaged interfaces ( (0 90),, larninate)

125

 

   

1.0 F 1 1 1

0.8 F — Reference Solution ‘
- - — - 1 - Sublarninate

0.6 F 1

 

 

 

0.4 F F

0.2 F +

0.0 1 4 1 m l 1 1 1 L g 1 1 1 \ l 1 4 1

-0.00010 -0.00005 0.00000 0.00005 0.00010
U

X

Figure 4.16 Throu h-thickness variation of inplane displacement in the
undamag region of clamped beam ( (0/90)88 laminate)

126

 

 

 

 

 

1.0 ,
0 8 _ — Reference Solution

' - - — - 1- Sublarninate
0.6 _ .

 

 

 

 

  

 

 

 

 

0.4 F ‘
0.2 r d
0.0 1 1 1 1 ‘1‘ 1

400.0 F200.0 0.0 200.0 400.0

Figure 4.17 Through-thickness variation of inplane normal stress at the fixed
end of clamped beam ( (0/90)& larninate)

1.0

0.8

0.6

0.4

0.2

127

 

 

 

 

 

—— Reference Solution
------------ 1- Sublarninate

 

 

 

 

 

 

 

-400.0 -200.0

0.0 200.0 400.0

0

Figure 4.18 Through-thickness variation of inplane normal stress at the center
of damaged region of clamped beam ( (0/90)8, laminate)

128

 

 

 

 

 

 

 

 

 

 

1.0 1 r - 1 1
\
\
\
" ' ' ' l
0.8 F 5
Y\
— Reference Solution
‘1 - - - - - - - 1- Sublarninate
0.6 — 1 ~
0.4 - 1 _ _ -
- “\
\
‘ -l
0.2 — ' —
— - - “ .l
\
\
0.0 l 1 l 1 11
-200.0 -100.0 0.0 100.0
a

XX

Figure 4.19 Through thickness variation of inplane normal stress in the
undamaged region of clamped beam ( (0/90/0)83 laminate)

200.0

1.0

129

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

I 7 fl j
I
I
l
0.8 fl
: l
l :1
0.6 ' d
Reference Solution
- - - - 1- Sublarninate 1
0.4 4
, :i
0.2 d
' \
i
I
0.0 ‘ ‘ ‘ l I
0.0 1.0 . 2.0 3.0
o t

Figure 4.20 Through-thickness variation of transverse shear stress at the fixed
end of clamped beam ( (0/90)“ laminate)

130

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1.0 I I I a r
l 1
l
l
I

0.8 ~ : d
I
l :L
I

0.6 ~ ' 4
: Reference Solution
: --—- 1- Sublarninate +
l
I

0.4 F 1 -
I
l I
I
l

0.2 ~ : ~
I
l l
l J
l

0.0 H 1 I . 1 . I

—1.0 0.0 1.0 2.0 3.0
all

Figure 4.21 Through-thickness variation of transverse shear stress in the
damaged region of clamped beam ( W90)“ larninate)

4.0

131

 

 

 

 

 

 

 

 

 

 

 

1.0 i I , T
I >—
I
I
i
I

0.8 — g ‘
i a
I
I
I
l

0.6 — : -
I Reference Solution
: ---- 1- Sublarninate
:
I

0.4 ~ : q
' 2
I
I
I
I

0.2 r g r
I
' >-
I
I
I

0.0 L I l 4 1 4

-1.0 0.0 1.0 2.0
OX2

Figure 4.22 Through-thickness variation of transverse shear stress in the
undamged region of clamped beam ( (0/90)« laminate)

3.0

132

1.0

 

gl I 0.25L I I
L

< i»

 

Figure 4.23. Tip-loaded cantilever beam (II/90L;s beam with dalamination
in the center of beam near the clamped end

133

 

  

 

 

 

 

 

1.0 a . v \ l f s i
0.8 — Reference Solution _
-—- 1-sublaminate
0.6 ~ _
0.4 - I
0.2 ~ _
\
0.0 ‘ L 1 * ‘ ‘ . l A X L L 1
-4.0e-05 -1 .5e-05 1 .Oe-05 3.5e-05

U

Figure 4.24 Through-thickness variation of inplane displacement in the
delaminated region of clamped beam ( (0/90)“ laminate)

134

 

 

 

 

 

 

1.0 \ I r i
\
\ Reference Solution

\ —-—- 1- Sublarninate
0.8 r a
0.6 — ~
0.4 — ~
0.2 ~ _
0.0 . 1 I 1 A 1 \
-4.0e-O5 -2.0e-05 0.0e+00 2.06-05 4.0e-05

U

Figure 4.25 Through-thickness variation of inplane displacement at the interface
of damaged and undamaged regions of clamped beam ( (0/90)“ laminate)

135

 

1.0 \ Y , , , , , f

Reference Solution
\ — — - 1- Sublarninate

 

0.4 ~

0.2 ~

 

 

 

I

o. L I r A l
-4.00e-05 -2.00e-05 0.00e+OO 2.006-05 4.006-05
' U

 

0

Figure 4.26 Through-thickness variation of inplane dispalcement in the undamaged
region of clamped beam ( (0/90)“ larninate)

136

 

 

 

 

 

 

 

 

1 .0 \ I I I
\
\
\
\
\

0'8 h _ Reference Solution d

\ \ — - - 1- Sublarninate J
0.6 ~ ‘
0.4 l _

\
_\_
0.2 L d
\
\
\
\
\ \
0.0 A 1 A 1 r L \
-200.0 -100.0 0.0 100.0 200.0
0

XX

Figure 4.27 Through-thickness variation of inplane normal stress near the fixed
end in the delaminated region of clamped beam ( (0/90)“ laminate)

137

 

 

 

  
  

 

 

 

 

1.0 , \ \ a l v .
\
\
_ \
0.8 - _
\
—— Reference Solution
——— 1- Sublarninate
0.6 ~ ~
0.4 P —
\—A
0.2 - _
\
\
\
\
0.0 * 1 A 1 4 \
-200.0 -100.0 0.0 100.0

XX

Figure 4.28 Throu h-thickness variation of inplane normal stress in the
undamage region of clamped beam ( (0/90)43 laminate)

200.0

 

CHAPTER 5
THREE DIMENSIONAL ANALYSIS OF WOVEN
COMPOSITES

5.1 Introduction

In the present chapter, a new analytical/numerical model is developed to analyze
woven composites. The aim of the present study is to obtain reasonably accurate results
with minimum modeling and computational effort. The unit cell of the woven composite is
considered to evaluate the effective stiffness properties. The unit cell of the woven com-
posite is discretized with eight node brick elements with one element through the thick-
ness of the cell. Based on assumed fiber geometry, the fiber volume fraction and average
fiber inclination in an element can be determined as a function of its spacial location. The
average stiffness properties of each element in the unit cell are determined from the fiber
volume fraction and constitutive properties of layers in the element using effective moduli
theory [Chou, et al. (1972)]. These effective stiffness properties are given as input to the
finite element model. A superposition method is used to determine the overall effective
properties of the woven composites. A brief description of assumed fiber geometry and

effective moduli theory are given below.

138

139

5.2 Fabric Geometry - Configuration in Unit Cells

All woven fabrics consist of two sets of interlaced threads known as warp and fill
threads. The type of fabric can be identified by the pattern of repeat of interlaced regions.
Two basic geometric parameters are required to characterize a fabric. nfg denotes that a
warp thread is interlaced with every nfgth fill thread and nwg denotes that a fill thread is
interlaced with every nwgth warp thread. For non hybrid fabrics, the geometric parameters
nfg and nwg are equal (ng=nfg=nwg). Fabrics with 11g 2 4 are called satin weaves. Figure
5 .1 shows fiber architectures of plain weave (ng=2), S-hamess satin weave (ng=5) and 8 -
harness satin weave (ng=8). The unit cells of the woven composites are shown in dark
boxes.

The unit cell of a woven composite can be broadly divided into regions characterized
by no waviness (cross-ply region; denoted region A), waviness in one direction (either
warp or fill has undulation; denoted region B) and waviness in two directions (both warp
and fill threads undulate; denoted region C). Obviously there are subcategories within
each of these three primary categories that distinguish whether the warp or fill tows are on
top or are undulating, etc. For brevity, only representative cases are described in detail
here (see Figure 5.2).

In the present study, the geometry of tow cross sections and undulations are assumed
to be sinusoidal. Figure 5.1, which is a plan form of the woven composite, does not show
the details of tow cross sections and undulations. The division of a unit cell of plain weave

and 5- harness satin weave with regions A, B and C is shown in Figure 5.2. By assembling

140

the regions A, B and C suitably, keeping the subclassifications in mind, one can generate
unit cells of any woven composite.

The geometry of the fiber in the regions A, B and C is shown in Figure 5.3. The origin of
the coordinate system is fixed at the lower left comer of the cell and the x3 axis is assumed
to lie along the thickness direction of the composite. The thickness of the woven compos-
ite, t is assumed to be 2h. It is necessary to mention that only one of the possible fiber
geometries in regions A, B and C is shown in Figure 5.3. The other fiber geometries
belonging to each group can be easily visualized by rotating the fiber geometries about the

x3 axis.

Region - A (Cell dimensions a0/2 x a0/2 x 2h)
The geometric parameters, h1(x1,x2), h2(x1,x2) and h3(x],x2) of warp and fill threads

are given below (see Figure 5.3).

“o “o
for 0<x2<§ and O<x1<3

h3(x1, x2) = 0
Based on the above geometric parameters, the thickness of warp and fill fibers can be cal-

culated as

tf(x1,x2) = h .

Region -B (Cell Dimensions an x a0/2 x 2h)

141
In this region, the geometric parameters h1(x],x2), h2(x1,x2) and h3(x],x2) and
h4(x1,x2) of warp and fill threads are given below.

“0
for O<x2<—

h1(x x) 3+cos 1t_ ,1 O<x <2
1' 2 = au 2 1 2

_ “x1 h (5.3)
h2(xl,x2) — [3-COS[a—u ]§ O<x1<au

h3(xl,x2) = h2(xl,x2)-h O<x1<au
nx__1 h a“
h4(x1,x2) = 1 + cos “u 2 -2—<x<au

As geometric parameters depend only on the x] coordinate, the thickness of warp and

fill fibers can be written as,

“0
for O<x2<—
2
tw(x1) = h O<xl <au (5.4)
a

The local angle between the warp tow and the xl-axis of the global coordinate system
is given by

1 dh2(xl, x2))

_. a
9w(x1) = tan ( dx 0<x1<au and O<x2<—0
l

2

—1[ . (7‘11)th (5.5)
= tan Sin — —
au 20“

In this region, the fill fiber is assumed to be oriented along the x2 axis.

142

If we divide the region B into small rectangular solids with m elements in the
xldirection and n elements in the xz-direction and one division through the thickness, the
geometric parameters of warp and fill threads in an element P in the domain 0<xl<au/2

and O<x2<a02 can be written as (see Figure 5.4):

t: = h on all edges
t}; = tf(x'1") on edge 1 and edge 4 (5-6)
= tf(x'1" + l) on edge 2 and edge 3

Assuming the element P is sufficiently small, the average thickness of warp, fill fibers and

the resin is taken to be

25 = h

-P 1 1

If = -2-[tf(x',")+tf(x'1"+ )] (5.7)
if = 2h—l5—l;

The average warp fiber orientation in the element can then be written as,

65 = $[9,(x’,">+ew(x',"“)1 (5.8)

Region- C (Cell Dimensions an x an x 2h)
If we divide the region C into small rectangular solids with m elements in the x1 -
direction and n elements in the x2- direction and one division through the thickness, the
thickness of warp and fill fibers in an element P in the region 0<x1<aul2 and O<xz<aul2

can be written as (see Figure 5.5):

143

t: = tf(x'2') on edge 1 and edge 2
t: = tf(x; + l) on edge 3 and edge 4

and (5.9)
t; = tf(x'1") on edge 1 and edge 4

P m+l

If = tf(x1 ) on edge2andedge3

Assuming once again that the element P is sufficiently small, the average thickness of

warp, fill fibers and thickness of resin in the element P is given by:

n+1

iw = %[tf(x;)+tf(x2 )]

-P l m I

If = §[tf(x'1n)+tf(r1 + )] (5:10)
-P .P .P

t, = 2h—tw—tf

The average fiber orientations of warp and fill fibers in the element P can be written as:

9,, = %[Ow(x'1") + 9W(x'l" + 1)] about the x2 axis
(5.11)

n+1

9] = %[(9”,(x;)-+-9W(x2 )] about the x1 axis

5.3 Effective Moduli Theory

The basic assumption of the effective moduli theory proposed by Chou, et al. (Chou, et al.
1972) is a combination of Voigt’s (constant strain) and Reuss’s (constant stress) hypothe-
sis. In particular, it is assumed that: l) the normal strains parallel to the layers and shear
strains in the plane of layers are uniform and equal for all the layers, and the correspond-
ing stresses are average values; 2) the normal stress perpendicular to layers and the shear

stresses in the planes perpendicular to the layers are uniform and equal in each material,

144

and the corresponding strains are average values. With these assumptions both the equilib-
rium at the interfaces and compatibility of the material are satisfied.

Figure 5.6 shows an element in a multilayer laminate with N layers. The layer number
in the element is denoted with index k. The plane of the layering is assumed to be the x1-
x2 plane and the laminate is stacked in the x3 direction. According to the effective moduli
theory, the assumptions regarding strain and stress quantities are given below. The overbar

quantities represent the average quantities in the laminate thickness.

_ _ k _ _ _ k . _ _ k

E11 ‘ 311’ 822 ’ 822’ 812 ‘ 812

_ _ k . _ _ k , _ _ k

c’13 " C513' 623 ‘ “23' 633 ‘ 033

and (5.12)

_ k k _ k k _ k k

“11 = V011; °22=V°2zi “12: V012

_ k k - k k .. k

813 = V 813* ‘523. = V 823’ 833 = Vk533

k k ..
where 0' i J and 81.]. are the stress and strain components in the layer k. 0' i J and éij are the
average stress and strain in the laminate. V1‘ is the volume fraction of the layer It. In the
above equation, indicial notation is used, wherein summation on repeated indices is

implied(1SkSN).

Hooke’s law applied to each layer can be written as
o = C e (5.13)

The stress and strain components in a layer can be written as,

k k k k k k k T
9 = {011,522,033,523,513’012} (5-14)

and

145

k k k k k k k T
§ = {311,522,539823,513,512} (5°15)

The relations between average stresses and average strains using average constituent

material constants can be written as:

(:3 = E g (5.16)
Equations (5.12), (5.13) and 16 represent 12N+12 equations in 12N+12 variables.
The average constitutive constants, Cij can be obtained as [Chou, et al. (1972)]:
_ N k k , .
C, = 2 V Cu 1,) = 1,2, 3,6 (5.17)
k = 1
N k k
6,,- = [ 2 V (CS/A ):|/A i,,j = 4,5 (5.18)
k =
where
N N N 2
k k k k k k
A = [ 2 V (CL/A M 2 V (Cigs/A )]—[ 2 V (CZS/A )] (5.19)
k = l k = 1 k = l
and
k _ k
A " Ck44clgs‘c4sclg4 (520)

If the fibers in the medium are inclined, the constitutive properties in the fiber directions
are transformed to the global coordinate system (x1,x2,x3), depending on the fiber orienta-

tion.

 

146

5.4 Finite Element Analysis for Effective Stiffness Properties

The unit cell of the woven composite is first identified and it is divided into regions A,
B and C. If there is a symmetry in the unit cell, it is exploited in the analysis. The regions
A, B and C are discretized with three dimensional brick finite elements. Based on assumed
geometry in the regions A, B and C, the average thickness and orientations of warp and fill
fibers are determined in the elements depending on their spacial location. From the aver-
age thickness of warp and fill fibers, the average thickness of resin in an element can be
found. Knowing the average fiber orientations, the material properties of the warp and fill
fibers in the elements are transformed into global x1x2x3 coordinates. Using the volume
fractions of layers and the constitutive matrix in the global directions, the effective proper-
ties of the elements are determined from Equations (5.17) and (5.18). These elemental
properties are given as input to evaluate the stiffness matrix of the elements for the finite
element formulation. Both the extension and shear loadings are considered for the woven
composite to determine the extensional and shear stiffnesses. The dimensions of the unit
cell are assumed to be 2a1x2a2x2a3, respectively. The origin of the coordinate system is
fixed at the center of the unit cell.

For woven composites, when the material properties along two inplane directions are
the same, only six of the nine elastic constants (E11, E33, (712, \712, (713, 1713) have to be
computed.

To determine the moduli and Poisson’s ratios, three subproblems with the following
boundary conditions are considered.

Subproblem 1:

147
U1(a1,x2,x3) = U U2(x1,a2,x3) = O U3(x1,x2,a3) = O (5.21)

Subproblem 2:

U1(a1,x2,x3) = 0 U2(x1,a2,x3) = V U3(x1,x2,a3) O (5.22)

Subproblem 3:
U1(al,x2,x3) = O U2(x1,a2,x3) = O U3(x1,x2,a3) = W (5.23)

with

U1(—al,x2,x3) = O U2(x1,—a2,x3) = O U3(x1,x2,—a3) = O (5.24)
for all subproblems. U1, U2 and U3 are the displacements in x1, x2 and x3 directions,
respectively. For the three subproblems, the sum of the normal reaction forces on bound-
ary surfaces, x1=-a1, x2=-a2 and x3=-a3 are calculated and denoted as Fl(n)’ F200, F300 for
each case. The subscript n represents here the subproblem number.
The effective modulus in the x1 direction is determined as follows.

From the nine constraint forces, the coefficients on, B and the resultant load P1 are calcu-

lated from the following equations:

F1(1)+a F1(2)+B 171(3) .-. p1
F2(1)+°‘ F2(2)+B F2(3) = 0 (5.25)
F3(1)+a F3(2)+B F3(3) = o

The average normal strains and the average Poisson’s ratios can be determined as:

g = _ g = 1‘: g = 122
11 2a1 22 2a2 33 2a3
- - (5.26)
_ _ ‘822 _ _ ”833
V12 ' — v13 ‘ ‘—

811

148

The effective Young’s modulus in the x1 direction (F311) is calculated from the following

energy balance equation:

1 1.. _ 2
fl” = 21511 (811) (2012a2203) (5.27)

where 2a] 2a2 2a3 represents the volume of the unit cell.
The remaining Poisson’s ratios and Young’s moduli in other directions can be calculated
in a similar manner.

To estimate 612, 623 and 613, the following boundary conditions are used.

Fora”:
(5.28)
U3 = O on x3 = —a3,a3
Fora-23
U = Kx U = Kx on x = -az,a and x3 = —a3,a3
2 3 3 2 2 2 (5.29)
Fora-13
U = Kx U = Kx on x1 = —al,a and x = —a3,a3
l 3 3 l l 3 (5.30)

U2=0 0n x2=-az,az

where K is a scalar that can be chosen arbitrarily.

The reaction forces for the above boundary value problems are summed on the boundaries
and average stresses are calculated. The shear modulus is determined from the calculated

average stress and the strain applied.

149

5.5 Results and Discussion

In the present paper, effective properties of plain-weave composites will be predicted
and the results from the present model will be compared with those of Whiteomb (Whit-
comb 1989). Whiteomb analyzed symmetrically stacked plain weave composites using 3-
D finite elements. The unit cell used for analysis is shown in Figure 5.7, where resin
regions have been removed to increase clarity of the tow structure.

It should be noted that the effort needed to accurately model this simple woven archi-
tecture using a traditional three-dimensional finite element model is large and requires the
use of sophisticated state-of-the-art 3-D modeling packages. When the fiber and resin
regions are discretely modeled, the number of elements in the model becomes large (Das-
gupta, 1990). Above all, the effort in-prepan'ng the input data for material properties of the
elements can become overwhelming, as the elemental local coordinate systems vary from
one element to the other due to tow undulation. The present model has the ability to cap-
ture the effects of different layers in a single element through the thickness, so one needs
to discretize the unit cell region in the xl-xz plane only. Thus only a small number of ele-
ments is required to represent the unit cell. Also, based on the assumed geometry in the
regions of the unit cell, one can calculate analytically the volume fractions and average
tow inclinations of each layer in the element as a function of spacial location. From the
volume fraction and constitutive properties of each layer, the average stiffness properties
can be calculated and given as input in computing the elemental stiffness matrix. The sub-
problems for estimating the stiffness properties only differ by load terms, so the finite ele-

ment method can be effectively used to solve the subproblems efficiently. The complete

150

procedure can be automated, and the computational effort involved is much less than that
of a full three dimensional model [Whiteomb (1989); Dasgupta (1990)].
In the present work, plain weave composites with waviness ratio’s of 0.167, 0.25 and

0.5 are considered. The waviness ratio, A , of a plain-weave is defined as

A - ——a“ 531
- (ao+au) (° )

where au is the fiber undulation length and 30 is the length of straight portion. The width of

the plain-weave is defined as a=a0+au.

5.5.1 Fiber Volume Fraction

The portion of the plain-weave used for analysis is divided into regions A, B and C. As
the tow geometries in these regions are defined, the volume of tows can be calculated. In
Figure 5.8, a plot of fiber volume fraction versus waviness ratio, A , is given. It can be seen
from the figure that a waviness ratio of 0.3 can have a maximum fiber volume fraction of
0.88 when voids are not present. Also from the figure, it can be seen that fiber volume

fraction decreases linearly with waviness ratio.

5.5.2 Effective Moduli and Poisson’s Ratio Estimation - Effect of Waviness Ratio

Figure 5.9 shows the discretization of the quarter model of a unit cell with 8- node
brick elements with various levels of mesh refinement. One element through the thickness
is used. Three types of finite element meshes are used. The coarse and medium mesh con-

tain 64 elements (162 nodes) and 100 elements (242 nodes), respectively. The fine mesh

151

contains 196 elements (450 nodes) in the model. The finite element mesh used by Whit-
comb (Whiteomb 1989) has 96 20-node brick elements (595 nodes) in the model.

While estimating the shear modulus, antisymmetric boundary conditions must be used
on the symmetry edges. These antisymmetric boundary conditions can be imposed
through multipoint constraints. But in the present study, to reduce the computational effort
involved in multipoint constraints, approximate boundary conditions given in Equations
(5.28) to (5.30) are used on the quarter model.

The material properties used for tow and resin regions are given in Table 5.1. Table 5.2
compares the effective stiffness properties obtained by the present model with Whiteomb’s
results as a function of waviness ratio for various mesh refinements. The moduli and Pois-
son’s ratio of the plain weave are normalized with respect to the equivalent properties of a
(0/90/90/0) laminate with the same dimensions. It can be seen from the table that the
present model with 64 elements is able to give reasonably accurate results compared to
Whiteomb’s results. We can also see that, in the case of E1 1, E33, 512 , the results
obtained by the present model agree very closely with Whiteomb’s results for all waviness
ratios. The maximum difference between the results is around 5%. In the case of
912, 1732, V13 there is considerable difference between the present results and Whit-
comb’s results. The Poisson’s ratio values, \712 given by Whiteomb decrease with wavi-
ness, which is unexpected. But the present results show a slow increasing trend. Whiteomb
pointed out this discrepancy in his results (Whiteomb, 1989) and mentioned that test val-
ues show an increase of \712 with increasing waviness, which is predicted by the current
model. The present transverse shear moduli results (313 ) decrease with waviness in con-

trast with Whiteomb’s results. The present model is not able to capture an increase in

152

transverse shear modulus with waviness ratio. Except for the transverse shear modulus
613, all other stiffness pr0perties predicted by the present theory agree with the predic-
tions of a complex three dimensional finite element model. The results demonstrate that
the present model which is both simple and computationally efficient can give very good

stiffness predictions compared to complex three dimensional finite element models.

5.5.3 Effective Moduli and Poisson’s Ratio Estimation - Effect of Width to
Thickness Ratio (alt)

Figures 5.10 to 5.13 show the variation of E1 1, 612, \712 and\723 versus width to
thickness ratio (alt) for the plane weave composite for various waviness ratios. The results
are obtained using the coarse mesh. The material properties given in Table 5.1 are used in
this case. It can be seen that effective stiffness increases with increasing alt ratio. In case of
E33 , the present model predicts no change with alt ratio. It is observed that the Poisson’s

ratio, \713 is almost independent of waviness ratio.

5.6 Preliminary Conclusions

The effective properties of woven composites are estimated using a new analyticall
numerical model. The fiber undulation and cross section are assumed to be sinusoidal. The
unit cell of the plain weave is identified and it is discretized with three dimensional brick
elements. Based on the assumed tow geometry, the volume fraction of the fiber and its
average orientation in an element are determined spacially. Using the effective moduli the-

ory, the effective properties of the element in the unit cell are determined. These effective

153

properties are given as input to the finite element model. The problems of determining the

elastic stiffnesses and Poisson’s ratios are divided into three subproblems. These subprob-

lems differ only by the load cases, hence the finite element method can be effectively used

to solve these subproblems efficiently. From the results, the following conclusions can be

drawn.

1. Except for the transverse shear modulus 613, all other stiffness properties predicted by
the present theory agree with the predictions of a complex three dimensional finite ele-

ment model.

2. The present model is simple in terms of finite element modeling compared to a full 3-D
finite element model in which tows are represented discretely. As the tow and resin
regions are not modeled explicitly in the current model, the number of elements
required is small, and hence the analysis can be carried out with less computational

effort compared to a full three dimensional model.

3. The present analysis procedure can be easily extended for other fiber architectures (5 -

harness satin weave, 8- harness satin weave, etc.).

Table 5.1: Material Properties of Tow and Resin [Whiteomb (1989)]

154

 

 

 

 

 

 

 

 

 

 

 

 

 

 

E E E G G G
Material 11 22 33 V12 V23 V13 12 23 13
mo10 1.1010 no10 no10 mo10 xlo10
Tow 13.4 1.02 1.02 0.30 0.49 0.30 0.552 0.343 0.552
Resin 0.345 0.345 0.345 0.35 0.35 0.35 0.128 0.128 0.128

 

 

155

Table 5.2: Effect of Waviness Ratio on Moduli and Poisson’s Ratio (Moduli and

Poisson’s ratios are normalized with respect to a 0/90/90/0) laminate)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Waviness Approach E11 E33 012 G 13 V12 V13
Present ”suns 0.9380 0.9341 0.9665 0.9527 1.0343 1.0050
(meshl, 64 elements)
Present results
0.9364 0.9361 0.9668 0.9537 1.0247 1.0048
(mesh2, 100 elements)
0.167
Present results
0.9340 0.9367 0.9745 0.9546 1.0198 1.0029
(mesh3, 196 elements)
(Whiteomb 1989)
0.92 0.95 0.96 1.10 0.87 1.100
(20-node 96 elements)
Present results 0.9069 0.9035 0.9529 0.9289 1.0399 1.0076
(meshl, 64 elements)
Present results
0.9064 0.9054 0.9535 0.9307 1.0384 1.0074
(mesh2, 100 elements)
0.25
Present results
0.8978 0.9061 0.9617 0.9318 1.0249 1.0033
(mesh3, 196 elements)
(Whiteomb 1989)
' 0.88 0.93 0.94 1.14 0.81 1.14
(20-node 96 elements)
Present resmts 0.8170 0.8146 0.9157 0.8578 1.0876 1.0160
(mesh 1 , 64 elements)
Present results
0.8169 0.8152 0.9170 0.8612 1.8 1.14
(mesh2, 100 elements) 0 6O 0 5
0.5
Present results
0.7746 0.8157 0.9242 0.8636 1.0242 1.0007
(mesh3, 196 elements)
("lntcomb1989’ 0.75 0.84 0.87 1.22 0.60 1.28

 

(20-node 96 elements)

 

 

 

 

 

 

 

 

 

 

IIIIIIII Fill Fiber

 

G Unit Cell

Resin

1|

 

 

 

 
 
 
 

 

 

 

 

 

 

 

 

 

""‘LLUJ:

”E

 

 

___——

__——

—___——
IIIIIIIIIIIIIIIIIIIIII

  

 

156

 

   

 

 

 

 

 

pun-W:

WEE

I: -III
Ellm l-Illlh WI II-

 

 

0.)

Q

H].U.U...lllllI é)
IIIIII |_ 1::
E '1:

E (U

U)

. - m

.. :
:llllé ob
3

   

JIIIIL.IIIIIJIIIILllfllJIIIIUIIIIIJIIIIIl IIIIII

 

 

 

 

 

 

 

 
  
 

""""||||

 

' '-:::::::: """" _ . ' . ' 11... '''' 11:112.. .'.1u.-11 .

b) 5- Harness Satin Weave

Figure 5.1 Various Fiber Architectures of Woven Composites

157

283 .53 “scene as 283.52.: e 8:8 =5 .3 9.38:5 dm 953.:

933-533 mmofiam -m 3

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

om mum m m
m<m<m<m<m<m
um m mumum
m<m<m<m<m<m
m m mum m
m<m<m<m<m<m
m m m m m
m<m<m<m<m<m
m mom m m
m<m<m<m<m<m
m m m m m

 

 

 

 

 

2.2625
95. E 30:38,? nonr—
5585

0:0 E 32363 Hopi
=o_wom 39380

0:603:35 Am

 

 

 

 

 

 

 

 

 

 

 

 

158

 

o as m .< sees. a tease .85..— .nm 8:3..—

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

a new»: a .585 .3 a tease 12E em 2&3

AOaH—NN «fin—vafixv AOaQNN «8,—XV

s.....€:.....\ a

. \1
.\

mx \\ .
: swam
m owum llY , ENE? .ESO
:4 owem

ENMENR Mingus ENJENX 183%

 

 

 

 

 

 

159
>2”
A

  

 

 

 

 

160

o :3»:— 5 seesaw .3 a tease age..— .m.m 953m

AOM—NK «TTEfiNV ACE—Nun «afixv

 

 

 

 

 

 

 

 

 

 

 

 

 

 

161

33:88— uemaEE—Z a 5 Eons—m .3. Rm 953%

 

 

 

 

 

 

 

 

 

 

 

 

162

885%.; Seam 25:35.3 3 :8 =5 .3 2:3...

 

 

owem— bHoEEIAm
\

   

ownm \mboEENAm

 

 

33% 805.53 5t» «ions—co 959$ 52m 2: E aux—E we :euoaum «Ea—oar ..o .853..an .w.m «Earn

.« .ozmm mmo£>m>>

o: . md . so . so . No

-

0.
go
o

 

1
1
C‘!
o

163
I I
L I
<0. 9.
O O

JBQL-J IO uogroerd etunloA

1
L
m.
C

1

 

 

 

L
Q
F

 

 

 

 

 

 

 

 

 

 

a) Coarse Mesh (64 Elements)

b) Medium Mesh (100 Elements)

164

 

 

0) Fine Mesh (196 Elements)

 

 

 

Figure 5.9. Various Finite Element Meshes Used for Analysis of Unit Cell

 

165

 

as: E as, 283 sea as 3% 8.3.62 2335 a6 8:33» 53 2:2..—

 

     

 

 

 

 

 

 

 

w a...
o? . ow . om . m4 . ow . 1 o and
m
Wu:
1 - . w
md 0 $6553: ow o e
. MN.O H wwGE>w>>Dlm W
SE n wwmcsmaolo m
- - omd m
MS
.3
A... I 8.0 v...
m.
c 0
F b p p b Okoo

166

33. an a? 283 :3.— 2 39 9.3.82 .85 28.9.: e 8:35» .:.m 2:5

Dow

 

 

8 . 8

omd 0 «355346114
mud n wmo£>m>>n+lm
mm to u mwo£>m>>®l®

 

 

L

 

I

J

l

 

oood

Ohod

omod

o
o:
O.
0

Q ‘smnpow Jesus euelduI

00790.
x

”Ol-

0:.0

167

 

33. S 5? AS 3 3:2 8.883.: :6 gears, .28 28:.

 

 

 

 

{a
co _. om om 0? ON 0
- and u mmo£>m>>oll® -
. mud n mwocsmamlm
E to u amasaaolo 1
An 0 -

 

    

 

 

omod

mmod

ovod

mvod

omod

‘oueu s,uossgod

ZLA

168

33. i as, A: 3 2.5: 6.586.— ? Seats» fie 2:9..—

00..

om

 

     

omd 0 383.2,:
mud n mmocsmamlm
word 0 mmocSmBQImV

l

 

 

 

 

omd

omd

OIIBH S,UOSS!Od

93 - _

CIA

CHAPTER 6
CONCLUSIONS AND RECOMMENDATIONS

6.1 Conclusions

The main goals of the present research to develop analytical and numerical tools for
analysis of laminated and woven composites have been achieved. The accomplishments of
the present research for analysis of laminated composites are:

1. Three new zig-zag higher order theories (fully cubic, quadratic and partly cubic)
for analysis of laminated composites and sandwich structures have been devel-
oped. Interlaminar shear continuity is enforced to make the number of unknowns
in the theories independent of number of layers in the model. These theories are

assessed for their suitability for finite element analysis.

2. Full cubic zig-zag theory requires C1 continuity for finite element analysis which
is not favorable for application in general purpose finite element packages. A spe-
cial formulation has been developed for quadratic and partly cubic theories to yield
C0 finite element formulations which are well suited for general purpose finite ele-

ment codes.

169

 

170

A new transverse normal strain has been derived by assuming the transverse nor-
mal stress to be constant through the thickness of the laminate using Reissner’s

mixed variational theory.

Closed form solutions (Navier Solutions) have been presented for the above theo-
ries and their accuracy has been assessed by comparing the results with exact elas-
ticity solutions in the case of a simply supported beam subjected to sinusoidally
varying load. Based on the computational efficiency and consideration of general
lamination schemes, a quadratic zig-zag theory is recommended over the other the-

ories.

The assumption of constant transverse stress leads to improved prediction capabil-
ities (transverse normal strain and the components of deflections and stresses),
especially in the cases where the material properties vary drastically through the

thickness of laminate.

A finite element based on a quadratic zig-zag beam theory has been developed and
the element satisfies all the requirements of a computationally efficient finite ele-
ment model for analysis of composite laminates, namely 1) the number of degrees
of freedom is independent of the number of layers; 2) the degrees of freedom
involved are only displacements and rotations; 3) no shear correction factors or

penalty numbers are present.

171

The Shear locking phenomenon has been eliminated by assuming consistent shear
strain fields in the formulation. The assumption of constant normal stress elimi-

nated Poisson’s ratio stiffening.

The results obtained from the present finite element were found to be in close
agreement with exact elasticity solutions for complex laminated cases. Overall the
element is both accurate and efficient, and provides excellent capabilities for ana-

lyzing thick and/or complex laminated structures.

The resulting element is compatible with other beam elements available in com-
mercial packages. The present element is incorporated into the MARC commercial
finite element code through user subroutine USEELM. The host program treats the

present element as a 12 degree of freedom element.

The following conclusions are drawn from the damage analysis of laminated compos-

ite structures using a damage model developed in the present study:

1.

A reasonably accurate and computationally efficient model for analysis of compos-
ites with damage has been developed. This model facilitates the analysis of com-
posites during progressive failure analysis and composites with delaminations. The
aim of the present model is to analyze damaged composites with just one element
through the thickness of the laminate. Towards this a new model was developed

and referred to here as the VZZKD model.

When the damaged layer thickness is small (1% of thickness), the present model

yields very good results even for large damage ratios (up to ll. = 105 ). Hence the

 

172

present model can be used to analyze composites with local delaminations accu-

rately.

3. When the damage layer thickness is significant compared to undamaged layers, the
results from the present model are in close agreement up to a damage ratio of 104.
At the damage ratio of 105 the present center deflections results are in 20% error
with exact elasticity solution for a simply-supported (0/90/0) laminate subjected to
a sinusoidally varying load. Hence at damage ratios lower than 104 the present
model with one element through the thickness can be used very reliable for analy-

sis of composites with simple lamination schemes containing ply-damage.

4. The present results show that the accuracy of the present model depends upon the
placement of damaged layers in the laminate and the percentage thickness of the
damaged layer. In the the case of laminates with multiple damaged layers of signif-
icant thickness, the local stresses and inplane displacements predicted by the

present model may contain significant error.

5. The commonly available FSDT elements in commercial software packages can not
be reliable after a damage ratio of 20 when the damaged layer thicknesses are sig-
nificant compared to undamaged layer thicknesses. When the damaged layer thick-
ness is small the results are reliable up to a damage ratio of 100. This study shows
that one has to be very cautious in using the FSDT elements available in commer-

cial packages for progressive damage or delamination analysis of composites.

From the analysis of woven composites following conclusions are drawn:

173

1. A new analytical/numerical model has been developed for predicting effective
properties of woven composites. The present model is simple in terms of finite ele-

ment modeling and in terms of computational effort.

2. Except for the transverse shear modulus 613, all other stiffness properties pre-
dicted by the present model agree with the predictions of a complex three dimen-

sional finite element model for plain-weave composites.

3. The present model can be easily extended for other woven architectures (5 -har-

ness satin weave, 8- harness satin weave).

6.2 Recommendations for Future Work

Zig-Zag theories are the most attractive theories for analysis of laminated composites
when the factors of accuracy and computational efficiency are considered. The main char-
acteristics of these models are that the number of unknowns in the model are independent
of the number of layers in the laminate. Unfortunately, the currently available zig-zag the-
ories give rise to finite element models that are either too computationally expensive or
poorly suited for general purpose analysis. The present quadratic zig-zag theory provides a
good solution to the above problem. The higher order zig-zag theories that are developed
in the present research are for analysis of laminated composite beams. For a detailed study
in the case of three dimensional structures like plates and shells, the present theories
should be extended by one more dimension. This can be achieved in a straight forward
manner.

The elements developed in the present research work are very accurate for analysis of

 

174

composites with multiple layers and complex lamination schemes. But these elements are
very expensive when used in the regions wherein significant transverse shear and
transverse normal stress effects are not expected. Hence for efficient modeling one
should model only the critical regions with the present element and noncritical regions
with conventional Mindlin elements (FSDT) available in the commercial programs. The
intermediate region between critical and noncritical regions are to be modeled with
transition elements that maintain the continuity requirements. Thus, there is a need for
development of a transition element for the present element for efficient analysis of large
scale structures.

For damage analysis of composites, there is a need to accurately model them after
first ply failure. The present VZZKD element with only medium computational effort can
not be expected to yield very accurate results in the case of complex damages that occur
in composites. The shear strain fields used in the present model are not sufficient enough
to capture the shear strain that exists in the damaged regions of the composite. Some
more research should be carried out to improve the shear strain fields of the present
element in the damaged regions with out increasing much computational effort.

For woven composites, the present model can be used for any fiber architecture to
obtain effective stiffness properties. But for the textile composites, where the fiber
architectures are very complex, new models are required for efficient and reliable
analysis. These models should substitute brute three dimensional analysis, which are the
only option at this time.

The present VZZKD element developed for damage analysis can also be used to

analyze woven composites. After discretizing the unit cell of the woven composite with

 

175

the present VZZKD element, the inclined tow regions in an element can be approximated
with a layer of average thickness lying parallel to the global coordinates (the global
coordinate axes are parallel to the length and thickness directions). The properties of tows
can be transformed into global coordinates using the average inclination in an element.
Thus the analysis of a woven composite can be approximated as an analysis of a
composite laminate with material properties discontinuous both in the plane and in the
thickness direction. The VZZKD element which is specially developed for analysis of
such laminates can be used to analyze woven composites with less computational effort
compared to a full three dimensional finite element model.

The new finite element models developed here for analysis of laminated composites
and woven composites have a lot of potential to be used in explicit analysis, implicit
analysis, structural optimization etc. The accuracy and efficiency of the present model for
the above analyses of composites can be a major area of future research as only very few

elements are available for such advanced analysis at this time.

‘ r M-_.'

T __“._-_‘.__.‘

APPENDICES

APPENDIX A

176

APPENDIX A
The thickness shape functions 1;, lgyand KEY for the HZZT 3 theory are given
below
k ~ ~ k-l ~ ~
111 = l-bzz-dz3+ 2(z-z,)(—bI3,-—de,.)
i=1
If, = bz2+dz3+ 2 (z-z,)(— b8,+d&i)
i= 1
k-l
k

21 = z + (— if: -113): + (— 57—82;)? + 2 (z -z,-)(a.--i>.-(af +138) —e.(2'~i+2§)
i=1
k—l
[:2 = 3122+Ez3+ Z (z-zi)(&l3,-+E&i)

i=1

Jll = 0
le = 0
k-l

1’51 = (_a;_}3;,)z’+(_a;—;921)z3+ 2 (z—z,)(—3?—Bp)13,-+(—E;—£ii2)8,.+&,-
i=1 ~
J’éz = (-Zz§—b?1)zz+ (-E§—dc~1)z3 + 2 (z—z,)(-a§—b?1)h.-+(—E§-d€~1)5.-+ 3,-

i=1

k — _§
I(ll —(1 h)

k Z
[(12-11
[(31-0
KI22=0

where

N
__ 2124.213,- 31124428,.

[5 b] = i=1 1:]
C a N A N
122+ 2(h-z,-)b,- 123+ 201—298,.

i=1 i=1 -

 

 

 

177

k—l
g 12+ 2 (h—z,)a,

i=1
NL-l

P = 2 (ll-29d,“
i= (A.l)

u—nu—c

k-
9: 201-2091
i=1

-1

E

3.-

‘i:

II
{M

to:
II
2
b.
1
«5

~.
ll
p—o

 

178

The thickness shape functions 1:”, 1E7 and K '67 for the HZZT 2 theory are given below.

k-l
11;] = l—hz2+ 2(z—z,)(—hl3,-—36i)
i=1
k-l
[If2 = hz2+ 2 (z—zi)(— 135,436,.)
i=1
~ k—l ~ ~ ~ ..
"£1 = z+(— ia—éb)z’+ 2 (z—z.)<a.-13.-(af+b§)-2.(Ef+d§):
i=1
k-l
1152 = fizz-1- 2 (z—z,)(&I3,-+E&i)

i=1

k
111: O
k
112 = 0
k
J —0
:1 (A2)
k z
K11=(1—}-1)
k 2
K12 = 72
k
K21: 0
k
K22 = 0
where
— N N --1
2h+ 213,- a

 

 

N
h + X (ix-298,- 2 (h—z,)a,

‘ i=1 i=1 -

179

The thickness shape functions 1:”, J [[37 and K EY for the HZZT3 theory are given below.

k—l
1’;] = 1—13z3+ 2(z—z,)(-1313,-Zze,)
i=1
k—l
172 = 523+ Z (z-zi)(— 1313,+Zi&,)
i=1
~ k—l ~ ~ ~ ~
J21 = Z+(—&}—b§)23+ 2 (z—zi)(&i-E,-(sz+b§)—£‘,-(Z'f+d§):
i=1
k-l
[$2 = Zzz3+ 2 (z—zi)(&hg+EE,-)

1:

(A3)

where

I——"'"I
O! 1::
RI an
I___l

ll
L»

It
2 u

 

 

N
h + 2 (h—z,)13, 2 (h-—z,.)e,.

" 1:1 i=1 -

APPENDIX B

180

APPENDIX B

Assuming the state of plane stress in x-z plane, the relation between the stress and strain in

a lamina can be written as

xx
ZZ

XZ

where the coefficients Qij are given below.

 

k k k k
Q11 = Ex/(l—vxzvxz)

k k
Q21 = 912

k

Q11
k

921
0

 

k _
k
Q12 0 EXX
Q22 0 822 (B-1)
0 Q; . sz
k k k k k
Q12 = Ex ver/(l'-V.rzvz1c)
k k k k
Q22 = 52/“ -szVu) (B.2)

k k
933:6

 

BIBLIOGRAPHY

 

181

BIBLIOGRAPHY

Aitharaju, V. R. and Averill, R. C. (1997a) “An Assessment of zig-zag Kinematic Dis-
placement Models for the Analysis of Laminated Composites”, Submitted for publica-
tion in Mechanics of Composite Materials and Structures.

Aitharaju, V. R. and Averill, R. C. (1997b) “A C0 Zig-zag Finite Element for the Analysis
of Laminated Composites Beams”, Accepted for publication in Journal of Engineer-
ing Mechanics, AS CE.

_~.-.“
I.

Aitharaju, V. R. and Averill, R. C. (1997c) “Three Domensional Analysis of Woven Com-
posites” Submitted for publication in Composite Science and Technology. 1'

Averill, R. C. (1994) “Static and Dynamic Response of Moderately Thick Laminated
Beams with Damage”, Composites Engineering, 4, 381-395.

Averill, R. C. and Reddy, J. N. ( 1992) “An Assessment of Four-Noded Plate FInite Ele-
ments Based on Generalized Third-Order Theory”, International Journal for Numeri-
cal Methods in Engineering, 33, 1553-1572.

Averill, R. C. and Yip, Y. C. (1996) “An Efficient Thick Beam Theory and Finite Element
Model with Zig-Zag Sublarninate Approximations”, AIAA Journal, 34, 1626-1632.

Bathe, K. J. and Wilson, E. L. (1976) Numerical Methods in Finite Element Analysis,
Prentice-Hall, Engle wood Cliffs, NJ.

Carrera, E. (1996) “C0 Reissner-Mindlin Multilyared Plate Elements Including Zig-Zag
and Interlaminar Stress Continuity”, International Journal for Numerical Methods in
Engineering, 39 1797-1820.

Cho, M. and Parmerter, R. R. (1993) “Efficient Higher Order Composite Plate Theory for
General Laminations Configuration’, AIAA Journal, 31, 1299-1306.

 

Cho, Y.B. and Averill, RC. (1997) “An Improved Theory and Finite Element Model for
Laminated and Sandwich Beams Using First Order Zig—zag Sublaminale Approxima-
tions”, Composite Structures, 37, 281-298.

Chou, P. C., Carleone, J. and Hsu. C.M., (1972), “Elastic Constants of Layered Media”,
Journal of Composite Materials, 6, 80-93.

182

Dasgupta, A. Bhandarkar, S. Pecht, M. and Barker, D., (1990), “Thermoelastic Properties
of Woven-Fabric Composites Using Homogenization Techniques”,. Proceedings of
American Society for Composites, Fifth Tech. Conference, Lansing, MI, June 12-14,
1001-1010

Di Sciuva, M. (1986) “Bending, Vibration and Buckling of Simply-Supported Thick Mul-
tilayered Orthotropic Plates: An Evaluation of a New Displacement Model”, Journal
of Sound and Vibration, 105, 425-442.

Di Sciuva, M. (1987) “An Improved Shear-Deformation Theory for Moderately Thick
Multilayred Anisotropic Shells and Plates”, Journal of Applied Mechanics, ASME, 54,
589-594.

Di Sciuva, M. (1993a) “Development of an Anisotropic, Multilayered, Shear-Deformable
Rectangular Plate Element”, Computers and Structures, 21, 789-796.

Di Sciuva, M. (1993b) “A General Quadrilateral Multilayered Plate Element with Contin-
uous Interlaminar Stresses”, Computers and Strcutures, 47, 91-105.

Foye, R. L. (1989), “Mechanics of Fiber Reinforces Composites”, Fiber-Tex 1988, NASA
Conference, Publication 3038, 237-247

Fried, I. (1974) “Residual Energy Balancing Technique in the Generation of Plate Bending
Finite Elements”, Computers and Strcutures, 23, 409-431.

Hughes, T.J.R., Cohen, M. and Haroun, M., (1978) “Reduced and Selective Integration
Techniques in the Finite Element Analysis of Plates”, Nuclear Engineering and
Design, 46, 203-222.

Ishikawa, T. (1981) “Antisymmetric Elastic Properties of Composite Plates of Satin
Weave Cloth”, Fiber Science and Technology, 15, 127-145.

Ihikawa, T., and Chou, T.W., (1982a), “Elastic Behavior of Woven Hybrid Composites”,
Journal of Composite Materials, 16, 2-19

Ishikawa, T., and Chou, T.W., (1982b), “Stiffness and Strength Behavior of Fabric Com-
posites”, Journal of Material Science, 17, 3211-3220.

Ishikawa, T., and Chou, T.W., (1983), “One Dimensional Micromechanical Analysis of
Woven Hybrid Composites”, AIAA Journal, 21, 1714-1721.

Lee, C. Y. and Liu, D., (1992) “ Layer Reduction Technique for Composite Laminate
Analysis”, Computers and Structures, 44, 1305-1315.

Lee, S. W. and Pian, T.H.H., “Improvements of Plate and Shell Finite Elements by Mixed
Formulations”, AIAA Journal, 16, 29-34.

. =iq

183

Li, X. and Liu, D. (1992) “An Interlaminar Shear Stress Continuity Theory for Both Thin
and Thick Composite Laminates”, Journal of Applied Mechanics, ASME, 59, 502-
509.

Liu, D. (1995) “Zig-Zag Theory for Composite Laminates”, AIAA Journal, 33, 1163-1165

Liu, D. and Li, X. (1996) “An Overall View of Laminate Theories Based on Displacement
Hypothesis”, Journal of Composite Materials, 30, 1539-1560.

Lu, X. and Liu, D. (1992) “Interlaminar shear slip theory for cross-ply laminates with non-
rigid interfaces”, AIAA Journal, 30, 1063-1073.

MARC Analysis Research Corporation (1994) Volume D, User Routines and Special Rou-
tines.

Mindlin, R. D. (1951) “Influence of Rotary Inertia and Shear on Flexural Motion of Isotro-
pic Elastic Plates”, Journal of Applied Mechanics, ASME, 18, 31-38.

Murakami, H. (1986) “Laminated Composite Plate Theory with Improved In-Plane
Responses”, Journal of Applied Mechanics, ASME, 53, 661-666.

Naik, R.A., 1994, “Analysis of Woven and Braided Fabric Reinforced Composites”, NASA
CR-194930.

Pagano, N. J. (1969) “Exact Solutions for Composite Laminates in Cylindrical Bending”,
Composite Materials, 3, 398-411.

Phan, N. D. and Reddy, J. N. (1985) “Analysis of Laminated Composite Plates Using a
Higher Order Shear Deformation Theory”, International Journal for Numerical Meth-
ods in Engineering, 21, 2201-2219.

Pratap, G. (1994) “The Displacement-Type Finite Element Approach- From Art to Sci-
ence, Progress in Aerospace Sciences, 10, 295-405.

Pratap, G. and Somashekar, B. R. (1988) “Field-and Edge Consistency Synthesis of a 4-
Noded Quadrilateral Plate Bending Element, International Journal for Numerical
Methods in Engineering, 26, 1693-1708.

Pryor, Jr. C. W., and Barker, R. M. (1971) “A Finite Element Analysis Including Trans-
verse Shear Effects for Application to Laminated Plates”, AIAA Journal, 9, 912-917.

Raju, LS. and Wang, J .T. (1994), “Classical Laminate Theory Models for Woven Fabric
Composites”, Journal of Composite Technology and Research, 16, No. 4, 289-303.

Rao, K. M. and Meyer-Peining, HR. (1990) “Analysis of Thick Laminated Anisotr0pic
Composite Plates by The Finite Element Method”, Composite Strcutures, 15, 185-213.

184

Reddy, J. N. (1980) “A Penalty Plate Bending Element for the Analysis of Laminated
Anisotropic Composite Plates”, International Journal for Numerical Methods in Engi-
neering, 5, 1187-1206.

Reddy, J. N. (1984) “A Simple Higher-Order Theory for Laminated Composite Plates”,
Journal of Applied Mechanics, ASME, 51, 745-752.

Reddy, J. N. (1993) An Introduction to the Finite Element Method, Second Edition
McGraw Hill, New York.

Reissner, E. (1984) “On a Mixed Variational Theorem and Shear Deformable Plate The-
ory”, International Journal for Numerical Methods in Engineering, 23, 193-198.

Robbins, DH. and Reddy, J .N. ( 1993) “Modeling of Thick Composites Using a Layerwise
Laminate Theory”, International Journal for Numerical Methods in Engineering, 36,
665-677.

Simo, J. C. and Hughes, T. J. R. (1986) “On the Variational Foundations of Assumed
Strain Methods”, Journal of Applied Mechanics, ASME, 53, 51-94.

Srinivas, S. ( 197 3) “A Refined Analysis of Composite Laminates”, Journal of Sound and
Vibration, 30, 495-507.

Tessler, A and Dong, SB. (1981) “On Hierarchy of Conforming Timoshenko Beam Ele-
ments”, Computers and Strcutures, 14, 335-344.

Toledano, A. and Murakarni, H. (1987a) “A Composite Plate Theory for Arbitrary Lami-
nate Configurations”, Journal of Applied Mechanics, ASME, 54, 181-189.

Toledano, A. and Murakarni, H. (1987b) “A Higher-Order Laminate Plate Theory with
Improved In-plane Response”, International Journal of Solids and Structures, 23, ll 1-
131

Whiteomb, J .D., (1989), “Three Dimensional Stress Analysis of Pain Weave Composites”,
Paper presented at the 3rd Symposium on Composite Materials: Fatigue and Fracture,
Orlando, FL, 417-438.

Whitney, M. J. (1973) “Shear Correction Factors for Orthotropic Laminates Under Static
Loa ”, Journal of Applied Mechanics, ASME, 40, 302-304.

Whitney, J. M. and Sun, C. T. (1973) “A Higher-Order Theory for Extensional Motion of
Laminated Composites”, Journal of Sound and Vibration, 30, 85-97.

Yang, P. C., Norris, C. H., and Stavsky, Y. (1966) “Elastic Wave Propagation in Heteroge-
neous Plates”, lntemational Journal of Solids and Structures, 2, 665-684.

 

_ J1.

 

 

185

Zienkiwiecz, O. C., Taylor, R. L. and T00, J. M. (1971) “Reduced Integration Technique in
General Analysis of Plates and Shells”, lntemational Journal for Numerical Methods
in Engineering, 3, 275-290.

F_——IZ.__

 

4'"

"‘1111111111111