If.

1.114. $4...

5.. a
a.

 

‘1‘ SH". Ix
Hahn»? ,

:zwwwrwum

.‘c

33.7

.I ’
. £34..

llnln:
u {a

 

 

 

 

THEBB
2

2004-
méampoo LIBRARY

Michigan State
University

This is to certify that the
thesis entitled

SENSITIVITY STUDY ON THE ENERGY ABSORBED
DURING A QUASl-STATIC CRUSH OF A COMPOSITE
TUBE

presented by

Narasimhan Swaminathan

has been accepted towards fulfillment
of the requirements for the

M. S degree in Mechanical Ergineering

 

 

Mew

Major Professor’s Signature

0/? /03

Date

MSU is an Affirmative Action/Equal Opportunity Institution

PLACE IN RETURN BOX to remove this checkout from your record.
To AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE

DATE DUE

DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/01 c:/CIRCIDateDue.p65—p.15

 

 

SENSITIVITY STUDY ON THE ENERGY ABSORBED DURING A QUASI-STATIC

CRUSH OF A COMPOSITE TUBE.
By

Narasimhan Swaminathan

A THESIS
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
MASTER OF SCIENCE
MECHANICAL ENGINEERING

2003

ABSTRACT

SENSITIVITY STUDY ON THE ENERGY ABSORBED DURING A QUASI-STATIC

CRUSH OF A COMPOSITE TUBE.
By

Narasimhan Swaminathan

In this work we have captured the main energy absorbing features during a tube crush
experiment by simulating the same in ABAQUS/STANDARD. Mechanisms such as mate-
rial failure, delamination, friction between plies and friction between the initiator and the
tube are modeled and their effect on the total energy absorbed is analyzed by an ANOVA
(Analysis of Variance) study. Results of such analysis may be useful in reducing cost and

time in modeling, while numerically simulating a crash event for practical purposes.

Analogy between the behavior of concrete material model available in ABAQUS and
failure of composite materials is used to simulate material failure. Delamination is imple-
mented using the penalty based interface element developed as a user element subroutine.
This unique interface element technology enables modeling delamination and friction be-

tween plies and hence their effect on the total energy absorbed.

Sensitivity of the response to various factors is calculated using ANOVA by taking ratio
of sum of squares for an effect to the total sum of squares. Results of ANOVA concluded
that delamination is an important factor and needs careful consideration for accurate numer-
ical modeling of crash events. Sensitivity of the response to material strain energy release
rate, friction and delamination related parameters were studied and it was found that mate-

rial strain energy release rate and friction play a vital role in improving the response.

To my loving parents and sister.

iii

You are only to perform your duty without an eye on their fruits.
—Gita

iv

ACKNOWLEDGMENTS

I would like to express my sincere gratitude to my advisor Dr. Ronald C. Averill. He
has motivated and guided me to the successful completion of this thesis. I would like to
thank him for giving me an opportunity of being his Teaching Assistant for the graduate
and undergraduate level finite element course, because of which I was able to lay a strong
foundation on the fundamental concepts in the subject. I am indebted to him for the moral
support that he provided during difficult times encountered, and for the time he spent with
me during the course of my master’s thesis. As his student I was fortunate to have an oppor-
tunity of working for his company, Red Cedar Technology, as a CPT (Curricular Practical

Training) which gave me a great practical experience.

I would like to thank Dr. Dahsin Liu and Dr. Farhang Pourboghrat for being in my

thesis committee and for their valuable suggestions and guidance.

I am extremely thankful to Michigan State University for having funded my Master’s
program through Teaching Assistantships without which I would not have been able to

dedicate time and effort on this work.

I would like to thank my parents and my sister for their constant words of encourage-

ment and support of my ideas.

Last but not least, I am thankful to all my friends and colleagues who were always of

immense moral support.

TABLE OF CONTENTS

Page

List of Tables ...................................... viii

List of Figures ...................................... ix
Chapters:

1. Introduction .................................... 1

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

1.2 Aim of the Thesis .............................. 5

1.3 Mechanics Of Crushing .......................... 6

1.4 Literature Review .............................. 7

2. Fracture Mechanics Based Material Model .................... 17

2.1 MATERIAL MODEL ........................... 17

2.1.1 *ELASTIC ............................. 19

2.1.2 *CONCRETE DAMAGED PLASTICITY ............ 19

2.1.3 *CONCRETE TENSION STIFFENING ............. 21

2.1.4 Fracture Mechanics ........................ 23

2.1.5 *COMPRESSION HARDENING ................ 28

2.1.6 *CONCRETE TENSION DAMAGE ............... 28

2.1.7 *CONCRETE COMPRESSION DAMAGE ........... 30

3. Interface Element Technology .......................... 31

3.1 Modeling Delamination .......................... 31

3.2 Application of Penalty and Lagrange multiplier method ......... 33

3.3 Interface Element to model Delamination ................. 37

3.3.1 Specification of the Interface Element for Input File ....... 41

4. Analysis ...................................... 51

4.1 The Model ................................. 51

4.2 Concrete Material Properties ........................ 53

4.3 Interface Element Properties ........................ 53

4.4 Factors and Their Levels .......................... 54

4.4.1 Material Strain Ratio ........................ 55

4.4.2 Strength Of the Interface ...................... 56

vi

4.4.3 Relative Displacement Ratio (RDR) ................ 58

4.5 Calculations ................................ 61

4.5.1 Calculation of B .......................... 61

4.5.2 Calculation of Material Strain Energy Release Rates ....... 61

4.5.3 Calculation of Mode 11 Strain Energy Release Rates ....... 62

4.6 Preliminary Analysis ............................ 62

4.6.1 EXPERIMENT l ......................... 63

4.6.2 EXPERIMENT 2 ......................... 65

5. Analysis of Variance ............................... 67
5.1 ANOVA ................................... 67

5.2 Results and Discussion ........................... 70

5.2.1 Effect of Strain Ratio ....................... 76

5.3 Conclusion ................................. 83

5.4 Scope for future Research ......................... 84
Bibliography ...................................... 86

vii

LIST OF TABLES

Table Page
4.1 Interface Properties .............................. 54
4.2 Values used for Experiment 1 ......................... 63
4.3 Results for Experiment 1 ........................... 65
5.1 Factors and their levels ............................ 71
5.2 Response for Each Experiment ........................ 72
5.3 Contributions of each factor ......................... 73
5.4 Effect of Strain Ratio ............................. 78
5.5 Effect of Strain Ratio, (Detail) ........................ 82

viii

LIST OF FIGURES

Figure Page
1.1 Picture showing typical composite tube and initiators ............ 3
1.2 Picture showing set up for the crush experiment ............... 3
1.3 Picture showing crushed composite tubes .................. 4
1.4 Cross Section of an Axial Crushing of a Composite Tube, From [2] ..... 8
1.5 Load vs Displacement curve ......................... 8
1.6 Illustration of semi—apical angle ....................... 12
2.1 Uniaxial Tensile and Compressive behavior of a Composite ......... 18
2.2 Illustration of Dilatation angle ........................ 20
2.3 Tension Stiffening, (Type: Strain, option) .................. 23
2.4 Tension Stiffening, (TypezDisplacement) .................. 25
2.5 Figure Illustrating the size effect of fracture energy ............. 27
2.6 Simple Linear Hardening in Compression .................. 28
2.7 Illustration of Damage Material Model ................... 29
3.1 Application of Lagrange Multiplier method to a bar problem ........ 34
3.2 Bilinear Model ................................ 38
3.3 Updating tensile stresses ........................... 39
3.4 Updating shear stresses ............................ 40
3.5 Interface Element enables to model friction between plies ......... 41
3.6 Specification Of Interface element nodes ................... 45

ix

3.7 Illustration for derivation of y ......................... 47
4.1 The Model ................................... 52
4.2 Illustration of Strain Ratio Factor ....................... 56
4.3 Purpose of Strain Ratio Factor ........................ 57
4.4 Illustration Of Relative Displacement Ratio factor .............. 58
4.5 Purpose of Relative Displacement Ratio Factor ............... 59
4.6 Force vs Deflection with %inch radiussed initiator ............. 64
4.7 Force vs Deflection with {Binch radiussed initiator ............. 64
4.8 Illustrating delamination ........................... 66
5.1 Labeling the responses ............................ 69
5.2 Force vs Deflection for Experiment ((1)) ................... 74
5.3 Force vs Deflection for Experiment (c) .................... 74
5.4 Force vs Deflection for Experiment(ac), High friction coefiicient ...... 75
5.5 Force vs Deflection for Experiment (c), Low friction coefficient ....... 76
5.6 Effect of Strain Ratio ............................. 77
5.7 Deformed shape , Strain Ratio = 1.2 ..................... 79
5.8 Deformed shape, Strain Ratio = 2.4 ..................... 79
5.9 Deformed shape, Strain Ratio = 4.8 ..................... 80
5.10 Force vs Displacement, Strain Ratio = 1.2 .................. 80
5.11 Force vs Displacement, Strain Ratio = 2.4 .................. 81
5.12 Force vs Displacement, Strain Ratio = 4.8 .................. 81

5.13 Improvement in convergence through high Strain Ratio ........... 83

xi

CHAPTER 1

Introduction

1.1 Introduction

Aerospace, automobile industries and other researchers have developed various tech-
niques to prevent damage to the passenger cell during vehicular collision. This is directly
related to the energy absorption capability of the vehicle. Substantial efforts have been
made to improve crash resistance in a vehicle and this is termed as crashworthiness. Crash-

worthiness can be achieved in two ways

1. Minimum deformation of the exposed structural member

2. Maximum deformation up to a specified limit

In the former case the frame contents may be subjected to a high level of momentum and
force transfer which may be unacceptable, in the later case a large deformation in the frame
is allowed and the degree of force and momentum transfer to the passenger cell is main-
tained within acceptable limits. This can be achieved by using energy absorption devices at
strategically chosen regions of the vehicle. Since these devices undergo irreversible dam-

age they are discarded and replaced with new ones after a crash event.

One class of energy absorption devices are structures made of materials that undergo
substantial plastic deformation like steel or aluminum. During a vehicular collision, ki-
netic energy is used up in plastically deforming the structure and thereby reducing force

and momentum transfer to the passenger cell. Another class of structures include those that

undergo fracture and other complex failure mechanisms thus using up most of the kinetic

energy and maintaining acceptable limits of forces in the passenger cell.

Specially manufactured composite materials are excellent energy absorption devices.
They also tend to reduce the weight of the vehicles of which they are made due a high
strength to weight and stiffness to weight ratio. This has proved to be a boon in the
aerospace and the automobile industry, enabling manufacturers to manufacture lighter ve-
hicles and improve fuel economy. Composite tubes are placed at specific locations of the

vehicle to absorb the impact of the collision.

Experiments involving axial crushing of composite tubes have provided useful informa-
tion regarding different failure mechanisms that contribute to energy absorption capacity of
composite structures. These experiments typically involve crushing a composite tube with
a flat platen or an initiator. The energy absorbed is the area under the load-displacement

curve.

For a general pictorial description of the crush experiment, please refer to pictures 1.1,
1.2 and 1.3. Figure 1.1 shows typical composite tubes and initiators. Figure 1.2 shows the
general set up for the experiment while Figure 1.3 shows the formation of fronds and axial

crack between these fronds once the tube has failed.

Unlike materials that absorb energy by plastic deformation composite materials absorb
energy by a wide range of failure mechanisms and other factors like friction. Experiments
involving axial crushing of tubes have brought to light different factors that contribute to

energy absorption. Some of them are:

 

 

Figure 1.1: Picture showing typical composite tube and initiators

I s hit-9“,.

‘r

 

Figure 1.2: Picture showing set up for the crush experiment

(90/0/90). stitch mm. —- 6‘?
(90/0/90). stitch ma: (90/0/90). “m“ m“ T 6° (m a“

 

Figure 1.3: Picture showing crushed composite tubes

0 Friction between the initiator and the composite structure

0 Matrix deformation and failure

0 Fiber failure

0 Crack propagation due to Delamination

0 Internal friction between failure surfaces

Each of the above mentioned factors are known to contribute to the total energy ab-
sorbed and empirical formulae have been derived to explain each of them. Though based
on experimental results, it is necessary to identify the factors that contribute the most to
energy absorption. This would facilitate designers, analysts and manufacturers to create

better predictive models and design better crashworthy structures.

1.2 Aim of the Thesis

It is clear that a large volume of experimental results have provided researchers with
useful information regarding many features of composite failure. This trend however, im-
portant, one must not overlook the importance of numerical simulation of crash events. The
Finite Element Method and related software packages offer themselves as versatile tools to
model and simulate a crash event and has opened new horizon to the field of numerical
optimization and other sensitivity analysis. Thus numerical simulation has been able to
reduce experimental testing of structures over full range of design variables, to a more fun-

damental level.

The aim of this thesis is to represent quasi-static axial crushing of a composite tube,
and to perform a sensitivity analysis on a few parameters that have been thought to affect
the energy absorption of composite tubes during a crash event. This will enable researchers
to decide on certain attributes of the material while numerically simulating a crash event. It
will also provide better insight into those attributes that affect energy absorption the most,
" or in reducing time, effort and money spent in modeling attributes whose effects can be

neglected.

Analogy between the concrete material model available in ABAQUS and the failure of
composite structures is implemented to simulate the axial crushing of a composite tube.
This readily simulates the energy dissipated due to material failure. Delamination is imple-
mented with the interface element that was developed for Langley research Center (NASA)
by Dr. Antonio Pantano and Dr. Ronald C. Averill [1]. This unique and innovative technol-

ogy simulates delamination and also enables modeling the internal friction between plies

after delamination and hence its effect on energy absorption. The interface element incor-
porates delamination through fracture mechanics concepts and allows for the specification
for MODE I and MODE 11 strain energy release rates along with the corresponding strength
values. This enables a study of the effect of delamination itself as a factor for the sensitivity

study.

1.3 Mechanics Of Crushing

The axial crushing of the composite tube is most useful to study the mechanics of the

phenomenon. A brief description of the mechanism from [2] is given as follows.
Quasi-Static Crush of a circular tube

The experiment is performed by crushing the composite tube with a steel platen (Fig-
ure 1.4). Tubes fail due to the progressive collapse and formation of continuous fronds
that spread outward and inward causing a mushroom type failure. There is axial splitting
that is caused by the tensile failure along the circumference (hoop direction). Initially the
tube behaves elastically and the load raises to a peak value Pmax and then drops abruptly
(Figure1.5). At this stage a central intrawall crack of length LC is formed. In the post
crushing regime the load oscillates about a mean crushing load indicating progressive fail-
ure. The formation of the first peak and drop of the load is due to the formation of two equal
lamina bundles, bent inward and outwards, due to flexural damage. At this stage triangular
debris of pulverized material begins to form and its formation is due to the friction between
the bent bundles and the steel platen. The wedge is formed when the load starts oscillating

and subsequently its size remains a constant throughout the deformation process.

The behavior of the fibers in the tube depend on their orientation. Axially aligned fibers
bend inwards or outwards, with or without fracturing, according to their flexibility. Fibers
aligned in the hoop direction can only expand outwards by fracturing and inwards by either

fracture or buckling.

Delamination occurs due to shear and tensile separation between plies. The axial lami-
nae split into progressively thin layers, forming, therefore, cracks normal to the fiber direc-
tion, mainly due to fiber buckling, finally resulting either in fiber fracture or in intralarninar
shear cracking which splits the laminate into many thin layers. Cracks propagate usually
through the weak regions of the structure of the composite material through the resin rich
regions or boundaries between the hoop fibers, resulting in their debonding, or through the

interface between the hoop and axial fibers causing delamination.

1.4 Literature Review

The aim of this work is to determine the factors that contribute most to the energy ab-
sorbed during the crush of a composite tube through AN OVA (Analysis of Variance) study.

It is in order to mention a few of the many works that has already been done in this field.

As the composite tube is axially crushed many energy—absorbing mechanisms are no-
ticed and efforts have been made to predict their contribution to the total energy absorbed
through experiments and empirical relationships. Significant work in this field has been

done in [2, 3]. Three main failure modes were identified here.

P

l

 

 

PLATEN

 

 

PULVERIS I
WEDGE

‘.
|

FRON DS

 

CRACK LENGTH,L

Figure 1.4: Cross Section of an Axial Crushing of a Composite Tube, From [2]

Pmax_

Crushing Load

 

 

Axial Displacement

Figure 1.5: Load vs Displacement curve

1. Progressive crushing associated with large amounts of crush energy. This was called

as MODE 1.

2. Brittle failure of the component resulting in catastrophic failure causing little energy

absorption (MODE 2).

3. Folding and hinging similar to metal tubes with a medium energy absorption.

Experiments were performed on square tubes made of two different materials ( A and B)
and were crushed using a flat platen and Mode 1 and Mode 2 failure were studied. Material
A was a fiber glass composite with chopped strand mat orientation in the plane of the mat
and Material B was a commercial glass fiber and vinyl ester composite. Starting from the
exterior of the shell, the plies were laid up in the sequence [(90 /O/ 2Rc) / (2Rc / 0/ 90) /RC_75]

for material B.

Principal source of energy absorption in MODE 1 were intrawall crack propagation,
bending of fronds due to cracking, axial splitting between fronds, flexural damage of plies,
frictional resistance between laminates, frictional resistance due to the penetration of the

debris wedge, frictional resistance of the fronds sliding along the platen.

In MODE 2 not much of energy absorption was observed as the tube collapsed due to
buckling. A failure analysis was done and empirical formulae were developed to predict
the energy absorbed in various failure mechanisms.

It was concluded that,

1. Energy of friction between the annular wedge and the fronds and between the fronds
and the platen was estimated to be about 44% and 48% of the total, for materials A

and B.

2. Energy absorbed due to bending of fronds was 48% and 32% for A and B respectively

3. Energy consumed due to crack propagation was about 7% and 18% for A and B

respectively

4. Energy absorbed due to axial splitting at the four corners of square tube was 1% and

2% for A and B respectively

Haq and Newaz in [4] have studied the effect of failure modes on energy absorption in
unidirectional PMC tubes. The major failure modes considered were axial cracking of
tube wall, transverse cracking of fronds and delamination. Crushing standard pultruded
fiberglass/polyester tubes with flat platen and initiator of two different radii performed the
experiments. A hoop constraint was also used to see its effect on the peak load and hence
the energy absorption. It was shown that the initiator induces a splaying crushing mode
due to which a stable crushing mode is obtained. Longer the axial cracks the more compli-
ant do the fronds become thereby not absorbing sufficient energy as they fail in buckling.
The peak load was shown to be higher in case of using a flat platen than the one obtained
by an initiator. It was also shown that external hoop constraint reduced the peak load and
improved the sustained crushing load, avoided major transverse crack and caused delam-
ination. It was concluded that initiation of delamination in the early stages of crushing is
important to reduce the undesired high peak load value. It was also shown that the specific

energy absorption in tubes increases as the radius of the initiator decreases.

10

In [5], test of axial splitting crush were performed on composite tubes and energy ab-
sorbing mechanisms were identified as, propagation of splitting cracks, friction, bending
and delamination. Tests were performed to determine the fracture toughness, and other
properties related to tension and compression. In order to evaluate energy absorption in
a sharp bending due to radiussed initiator, a fixture was designed and empirical relation-
ships were derived to find the amount of energy absorbed due to friction as a function of

the bending angle. Also a formula was developed to find the load causing the axial splitting.

In their paper [6] Farley and Jones studied the effect of material and geometrical prop-
erties on the energy absorption. An equation similar to the buckling of a column was
formulated for the sustained crushing load. This equation had six geometrical and material
parameters associated with crushing of composite tubes, that could be useful in preliminary

design and understanding the crushing behavior of composites.

In [7] it was stated that a good correlation between the experiment and the built model
requires accurate representation of the structure, correct boundary conditions and the cor-
rect constitutive relationship. A fractional factorial study was done to see the effect of

material properties on the energy absorption.

It was concluded in [8] that dynamic testing dissipated less energy absorption than

quasi-static testing. Similar observations were recorded in [9] and it was specified that cir-

cular tubes had higher energy absorption than square tubes.

11

 

Composite Tube

(1, is the semi apical angle

Figure 1.6: Illustration of semi—apical angle

Experiments on square and circular frusta were performed in [10]. Effect of semi-apical
(Figure 1.6) angle and wall thickness was studied. It was reported that increase in the semi-
apical angle of the frusta tends to decrease the specific energy absorption while an increase
in wall thickness tends to increase the specific energy absorption. It was also reported that

specific energy absorption for circular frusta is higher that for square frusta.

In [11] crush experiments were performed on braided and CSM glass reinforced com-
posite tubes. Initiator plug were used to initiate progressive crushing and force vs deflec-
tion curves were obtained. Analytical model for the tube crush was developed based on
axisymmetric elastic deformation. The steady—state axial load was found using energy bal-

ance approach.

Composite tube crushing poses a high degree of non-linearity due to failure. Due to this,
the choice of the material for numerically modeled analysis becomes an important consid-
eration. The material model must be capable of predicting the material failure (matrix and

fiber) and delamination. Bi-phase rheological models are obtained by superimposing the

12

effects of an orthotropic material phase and of a unidirectional material phase, the fibers.

Each phase is assigned a different rheological law.

Composite materials present difficult challenges for crash modeling, this is because
their failure can be attributed to many mechanisms. One of the most important issues con-
cerning composite failure is delamination. As pointed out in [12] most common methods

to model delamination are

1. Using Spot welds to represent an interface

2. Cohesive Failure Model

3. Using fracture mechanics based approaches like VCCT (Virtual Crack Closure Tech-

nique)

As pointed out in [12] the spot weld method does not have a strong physical basis, the co-
hesive failure model is difficult to implement and the VCCT approach, though accurate re-
quires a high computational ability. These models are compared in [12] with MSC—Dytran
and it was concluded that in computationally using these models to predict delamination
the mesh size should be small and that the dynamic property data required to predict the

delamination growth are not easily obtained.

In [13] a material model capable of predicting damage accumulation and braider tow
scissoring was developed as a user material subroutine in ABAQUS/STANDARD. It was
found that the tri-axial braid [0/ i601 circular tubes exhibited crushing response that closely

resembles ideal force/deflection curve. A parametric study was done and it was concluded

l3

 

that friction played a major role and the radius of the initiator is also a crucial factor affect-

ing the energy absorbed.

The VCCT approach was also used in [14] to determine the strain energy release rates
in composite skin/stringer specimens for various combinations of uniaxial and biaxial load-
ing conditions. A method was developed to calculate the strain energy release rates as a

function of delamination lengths.

In the paper by Farley and Jones [15], the crack propagation was achieved by discon-
necting the coincident nodes of the finite element. It was pointed out that for successful
modeling of the energy-absorption capability, factors like crack initiation and growth, mate-
rial and geometric non-linearity, fracture of individual lamina or laminae, modeling around

bifurcation loads, element type and level of descretization, need to be considered.

Dr. Antonio Pantano and Dr Ronald. C Averill presented a novel interface element
technology [1] for modeling delamination in composite materials and will be described in

detail and also used in this work to simulate delamination.

The simulations involved in this work are implicit and quasi-static in nature and hence
inertial effects are neglected. This causes serious convergence problems during the analysis
and as mentioned in [16] it was required to use viscous damping that damps out the high

frequency responses in the specimen.

In [17] reasons that are given for such convergence problems. In continuum based soft-

ening models used in statics is that the partial differential governing equations change from

14

elliptic to hyperbolic in the softening regime. Due to this, the boundary/initial conditions
applicable to one class of differential equations are not appropriate to the other. As a result
it is not possible for any rate-independent softening to proceed within any non-zero volume
of the domain of solution [18]. Thus it is required to introduce what is known as ”localiza-
tion limiters”. These help to resolve the mathematical ill-posedness issues that occur due

to softening.

As this work involves the use of concrete material model to represent a composite struc-
ture, a reasonable understanding of numerical modeling of concrete yielding and failure
was required. Literature on this was found in a consolidation of many works, in a Report

on Finite Element Analysis of reinforced Concrete [19].

Constitutive relationships established to represent the behavior of concrete are mainly

of 3 forms.
0 Elasticity Based Models
0 Plasticity Based models
0 Plastic Fracturing Models

In the Elasticity based models the stress is a function of strain or vice-versa. This is

generally expressed as

Gij =Fij€k1 (1‘1)

15

 

In the above equation F is a tensorial material response function and is generally a
tangential or secant constitutive matrix. These models are suitable for monotonic and pro-

portional loading and fail to predict the response correctly under unloading conditions.

The plasticity-based models assume that concrete flows on yielding and use plasticity
theories like yield surfaces and flow-rules to represent the behavior. Incorporating con-
crete’s deviations from normal metal plasticity like pressure dependent yielding and in-

crease in volume under plastic straining (dilatancy) better represents these plastic models.

It was later investigated and concluded that plasticity models only partly represent the
concrete material behavior. The fact that, under low confining pressures concrete experi-
ences micro cracking and micro fracturing, which is accompanied by a reduction in stiff-
ness, lead to the development of the plastic fracturing model representation. The concrete
material model in ABAQUS is a Damaged Plasticity model and is based on the plastic

fracturing theory. Further description of this model is found in the next chapter.

16

CHAPTER 2

Fracture Mechanics Based Material Model

2.1 MATERIAL MODEL

Composite materials are known to absorb energy through a wide range of failure mech-
anisms. These can be broadly classified into delamination and material failure. In the
model considered for analysis both of these have been implemented. The material failure
is implemented using the concrete model in ABAQUS [20], and a novel interface element
technology that was developed in [1] is used to simulate delamination. In the following

paragraphs detailed explanation of the concrete material model is given.

The Concrete Damaged Plasticity Material Model

It is well known that composites have in general a brittle nature. Under tensile loading the
stress-strain relation is linear within the elastic limit and then material failure occurs which
is due to matrix or fiber failure, or both. Under compressive loading the material behaves
linearly within the elastic limit and post—elastic region in compression can be assumed to
be a hardening surface (Figure 2.1). This is true because there will be load transfer across
cracks formed due to compressive failure. The above-described model fits the concrete
material model available in ABAQUS to model reinforced and plain concrete structures

subjected to cyclic or dynamic loading under low confining pressures.

The model considers concrete as an isotropic material and hence the analogy between

concrete and composite only approximately determines the material and therefore this work

17

 

Stxess

Uniaxial Tension

 

Strain>

Uniaxial Compression

 

Figure 2.1: Uniaxial Tensile and Compressive behavior of a Composite

does not take into account the anisotropic nature of composite. This sacrifice does not in
any way affect the sensitivity analysis that is to be performed which is the main goal of this

work.

The Concrete material model used, as mentioned previously is a damaged plasticity
model. It is most convenient to describe and understand the model with the appropriate
keywords that is used to implement the model in ABAQUS (Please see [20]). Fundamental
concepts and theory are also explained with the keywords. The model is implemented using

the following keywords.

18

o *ELASTIC

o *CONCRETE DAMAGED PLASTICITY

o *CONCRETE TENSION STIFFENING

o *CONCRETE COMPRESSION HARDENING

o *CONCRETE TENSION DAMAGE

o *CONCRETE COMPRESSION DAMAGE

2.1.1 *ELASTIC

This keyword specifies the elastic properties of the isotropic material. The Young’s

modulus and poisson’s ratio form the data card to this keyword.

2.1.2 *CONCRETE DAMAGED PLASTICITY

The Concrete damaged plasticity model in ABAQUS uses concepts of isotropic dam-
aged elasticity along with isotropic tensile and compressive plasticity to represent the in-
elastic behavior of concrete. The keyword is used to define the flow potential, yield surface

and viscosity parameters. The various data card to be mentioned for the keyword are:

c Dilation angle, w: Dilation angle , in degrees in the p — q plane (p - First invariant of
the stress tensor, q- Von-Mises equivalent deviatoric stress, see Figure 2.2). Dilatancy
is a concept applicable to brittle materials that represents an increase in the volume
on undergoing inelastic strains. This is mainly due to the formation of microcracks.

The term dilatancy has in general not known to be applied for composite failure and

19

 

 

>P

 

Figure 2.2: Illustration of Dilatation angle

for all analyses considered it is given a value of 30 O.

0 Flow potential eccentricity: A flow potential describes the accumulation of plastic
strain as the body undergoes plastic deformation. This model uses the Drucker-
Prager flow potential function and requires a small positive number representing
the flow potential eccentricity at which the hyperbolic flow potential approaches an

asymptote. The default is 0.1.

0 Ratio of initial equi-axial compressive yield stress to the initial uniaxial compressive

yield stress: This data specifies the yield surface and has a default value of 1.16. This

20

value remains as a default for this analysis.

0 Kc: This value represents the ratio of the second stress invariant on the tensile merid-
ian to that on the compressive meridian at initial yield for any given value of the
pressure invariant such that the maximum principal stress is negative. The default

value of g is used.

c Viscosity parameter a: is used for the viscoplastic regularization of the concrete
constitutive equations. This is given a value of IE — 4 to facilitate numerical conver-

gence. Please see [17,18,20—23] for further details on the purpose of using this value.

2.1.3 *CONCRETE TENSION STIFFENING

When both concrete and the reinforcement are subjected to tensile stresses, large cracks
form in concrete. The bond failure between the reinforcement and concrete matrix also oc-
cur at the same time and relative movement between reinforcement and concrete takes
place. The shear forces at the contact feed tension stresses to the concrete between the
cracks. The concrete hangs on the structure and contributes to the overall stiffness of the

structure. This is referred to as tension stiffening.

The post failure region is modeled so that the stress decreases gradually from the ulti-
mate tensile stress to zero. This describes the cracks as cohesive cracks; this means that the
cracks are capable of transferring stress. This may be a valid assumption in case of com—

posite structures also, because there is definitely stress transfer by the interlocking fibers

21

as the matrix fails. The cracking model is a smeared cracking model. This means that at a
point when stress reaches the crack detection surface the crack is accounted by the way in
which it affects the stress distribution in the structure. Cracks are assumed to occur parallel
to each other in a finite element and this is numerically modeled by reducing the stiff-
ness in a direction perpendicular to the crack [18]. ABAQUS offers three ways to specify
the post failure region. These are available as options under the *CONCRETE TENSION

STIFFENING keyword.

0 TYPE=STRAIN: The first one is TYPE=STRAIN option (Figure 2.3). This gives the
post failure stress as a function of cracking strain across the crack. This approach is in

general both physically and computationally unacceptable due to following reasons

described in [18].
- The softening zone has zero volume and width.

— The inelastic strain and fracture work are zero.

— The solutions are un-objective. This means that the results do not converge to a
unique value as the mesh is refined. This is due the formation of narrower crack

bands and localization of strains.

Due to these disadvantages this method of specification was not used in the current

model. Please refer to [18] for more information.

o TYPE=DISPLACEMENT: This option overcomes the disadvantages of the previous
option. Here the stress is a function of the crack opening displacement rather that the
cracking strain (Figure 2.4). This method of defining the post failure region belong

to the crack band models and are acceptable both physically and computationally.

22

G>q

 

 

éfl----—--__---------------

Figure 2.3: Tension Stiffening, (Type: Strain, option)

The option is based on the model developed by Hillerborg in, [24]. It is based on the
assumption that the stress across the crack does not go to zero as long as the crack
is narrow. The data card that needs to be given for this option is the crack opening
displacement or the width of the crack at which the stress linearly decreases to zero.
The model can be easily calibrated for the current model by knowing the value of the
strain energy release rate for the material. Before describing the Hillerborg’s model

a few failure theories from fracture mechanics are in order.

2.1.4 Fracture Mechanics

Many theories have been developed that are used to predict the formation and prop-
agation of cracks and to simulate them using Finite Element Method. Fracture me-

chanics generally predicts material failure and crack propagation using energy based

23

approach.

Stress intensity factor approach.

In this theory the stresses in the crack tip theoretically approach infinity according to
the expression 0 = K / m where r is the distance from the crack tip and K is the
stress intensity factor that depends on load and crack geometry. Once K reaches a
critical value of Kc called the critical stress intensity factor, the crack grows. Using
this method along with FEM to predict crack propagation requires finer mesh at the

crack tip. This limits its application to complicated problems.

Energy Balance Approach

This was the model originally developed by Griffith. He concluded that the energy
Gc required to create free surfaces in a body was in general a material constant this

is given as the energy required to produce a unit area of crack in a material.

When work is done to deform a material a part of the energy is stored as elastic strain
energy and the remaining is available for forming free surfaces. If W is the external
work done and U is the strain energy stored, then W — U is available to form free
surface. If W — U is equal to CC * (free area produced due to crack) then the

crack propagates quasi-statically (refer to [18]).

24

>0

 

 

u+u

Figure 2.4: Tension Stiffening, (Type=Displacement)

Hillerborg’s Model

In the Hillerborg’s model the stress does not go to zero as long as the crack is narrow,
the limit being given by the crack opening displacement “0- Figure 2.4 illustrates the

model.

So the amount of energy required to open a unit area of this crack is given by

l
50'qu (2.1)

This energy can be equated to the strain energy release rate Gc of the material to

obtain the crack opening displacement.

25

 

The distribution of the stress as a function of the crack opening displacement can be

assumed to be linear for brittle materials like concrete and composites

For a linear relationship we therefore have

1
30'qu = Gc (22)
G

o TYPEzGFI: This option is the same as the TYPE=DISPLACEMENT option, except
that one can directly specify the strain energy release rate for the ultimate tensile
strength. This option has been used for all analysis in this work . The specification of
the strain energy release rate is independent on size of the finite element and hence
overcomes the disadvantages of strain localization and mesh dependency of the re-

sults.

Consider a square finite element of dimension 1 by l as illustrated in Figure 2.5. The
stress in the element reaches maximum value at a strain value of (if. The stress then
decreases to zero at a strain of 108f. After failure the finite element is in two pieces
separated by a crack of dimension 108 fl = “0- Thus we have converted the strain
value to a crack opening displacement assuming that the inelastic strain in the two
pieces in negligible. The strain energy release rate can be calculated according to
Hillerborg, [24] as
1

50-14110 = CC (2.4)

Specifying the strain energy release rate in this way resolves the illposedness is-

26

 

 

 

 

 

 

CRACK

 

 

 

 

 

I

Figure 2.5: Figure Illustrating the size effect of fracture energy

sues of the governing partial differential equations to a certain extent. This is true
because the crack no longer localizes on to a point but to a characteristic dimen-
sion in the specimen. This is therefore a kind of ”localization limiter”. Please refer
to [17,18, 20] for more detailed mathematical explanation. The specification of the
Stress-Crack opening Displacement or the Stress—Strain Energy release rate is mesh
objective. This however modifies the plastic modulus of the stress-strain diagram as

the mesh descretization is changed.

For the current work the value ” 10” will be reffered to as Strain Ratio (SR) . It is the
ratio between the strain at complete loss of strength to the strain at failure. This ratio
will be considered to be a factor that affects the response.The strain energy release

rate is dependent on the strain ratio through the expression

27

 

 

I...

Figure 2.6: Simple Linear Hardening in Compression

CC = %0u(l)(SR)ef (2.5)

2.1.5 *COMPRESSION HARDENING

This keyword is used to specify the hardening behavior of concrete in compression.
It requires the specification of maximum stress in compression as a function of inelastic
strain. For the sake of simplicity a linear hardening is assumed after maximum stress in

compression is reached. Please see Figure 2.6.

2.1.6 *CONCRETE TENSION DAMAGE

This keyword specifies the decrease in elastic stiffness as the stress decreases to zero.

The degradation of stiffness increases with increase in plastic strain. This specification

28

 

 

 

 

 

 

I”
E(I-—dt) 7’ 1’:
1’ I I
I ’ I
I’ |
I; ” I
’1 I : 8
81 T 8
A p e1
8ck

 

Figure 2.7: Illustration of Damage Material Model

adds plasticity to the model and is governed by scalar damaged elasticity as shown in Fig-
ure 2.7. The undamaged inelastic strain is decomposed into plastic strain and cracking
strain. The data is mentioned by specifying the scalar damage variable (0 S d g 1) as a
function cracking strain (Eckl- This is not used in the current research, due to difficulty in
obtaining the softening data. This keyword also has the option to specify the compression

stiffness recovery (0 S we 3 1) that may be required if loading is cyclic.

29

2.1.7 *CONCRETE COMPRESSION DAMAGE

Similar to the tension damage, the model also allows modeling crushing in concrete
by stiffness degradation. Option for tension stiffness recovery is also provided to enable
modeling cyclic loading. This is not used in this model due to difficulty in specification
of softening data. The specification of this may lead to convergence difficulties due to de-

crease in stiffness during compression.

30

CHAPTER 3

Interface Element Technology

3.1 Modeling Delamination

The interface element developed to join independently modeled finite element meshes
and simulate delamination in composite structures has been used to incorporate delamina-
tion, [1]. The interface element is a penalty-based model, this means that the constraints
are enforced through a penalty parameter. A brief description of the interface element and

its method of specification for the model are in order.

The motivation to develop the interface element based on penalty method comes from

the Lagrange multiplier method, which is described below.

The following examples from [25] illustrates the two methods.
Lagrange Multiplier Method

Consider the minimization of a function f (x, y) = 4r2 — 3y2 + ny + 6x — 3y + 5 subject to
the constraint G(x, y) 2 2x + 3y = 0. Using the Lagrange multiplier method the equation is
F L(x, y) = f (x, y) + M2): + 3y), by taking the partial derivative of this equation with respect
to x,y, and A and equating to zero, three equations are obtained that can be solved for x, y

and A.

61?
—a—L:8x+2y+6+27t=0 (3.1)
X

31

aFL

___—_—6y+2x—3+3}\.=0 (3.2)
3y
BFL
__ = 2x 3 = 0 3.
a). + y ( 3)
This givesx = —3, y = 2 and A = 7.
Penalty method

Consider the minimization of the same function, using the penalty method we consider
the equation in the following manner, Fp(x, y) = f (x,y) + %(2x+ 3y)2. By differentiating

this equation with x, y and equating to zero we get x and y in terms of the penalty parameter .

BF
31-)— - 8x+2y+6+2y(2.x+ 3y) (3.4)
8F
-—B— —6y+2x—- 3 +3y(2x+3y)= 0 (3-5)

This gives x = :125—6—3—1637 andy =% and clearly as y —> co the value of x and y

reach the value given by Lagrange multiplier method and this also happens to be the exact
solution to the problem. Though the accuracy of the problem depends on the value cho-
sen for y when used in finite element method or the potential energy method this greatly

reduces the size of the stiffness matrix.

A high value of y predicts results more accurately. The reduction in the size of the ma-
trix results because the function is not differentiated with respect to the penalty parameter
and hence there is lesser number of primary degrees of freedom to solve for, unlike the La-
grange multiplier method. Lagrange multiplier method does in fact impose the constraints

accurately and hence serves as a limit to the accuracy of finite element solutions obtained

32

by the penalty method.

3.2 Application of Penalty and Lagrange multiplier method

Consider a problem of an axial load applied to a bar. To solve this problem the Total
potential energy (T PE) of the bar is developed as the sum of the strain energy and the work
done by external forces. We have for a bar of length I acted upon by external forces, the

potential energy is given by:

1 I du 2
111.. E[015.4(5) dx—fl-ui (3.6)

, where E is the Young’s modulus, A the area of cross section and u the displacement as a
function of x, f,- the applied external forces on the nodes and the corresponding displace-
ments. Now consider two elements of the bar, each of length I joined together at the node

as in Figure 3.1

The TPE of the bar will alone be

1 I du 2 l 2’ du 2 4
n1 —§/OEA(E;) dx+§l1 1511(5) car—Elm,- (3.7)

In ordinary problems this functional is minimized after substituting approximate func-
tions, for the displacement field, u. Now, if we have the interface element, which has its
own displacement field V, connecting the two bar elements, the potential energy of the bar
is subjected to two constraints V — u2 = O and V — u3 = 0 where uz and u3 are the global

primary degrees of freedom for the second node of the first bar element and first node of the

33

21

3‘ . .
l.

 

 

“2W3 u.

,nr. , .rP

Interface Element

 

 

 

 

Figure 3.1: Application of Lagrange Multiplier method to a bar problem

second bar element. Thus in using the Lagrange multiplier method to minimize the TPE

subject to the constraints specified we have the TPE of the system to be

1 I du 2 1 21 du 2 4
[1:5/0EA(E) dx+§/l 511(3) dx—iglfiui+}.1(V—u2)+}\.2(V—u3) (3.8)

Now to solve this problem we first replace a with its approximation functions. Using

Lagrange interpolation functions we have

2
u: 2 ujllij (3.9)
i=1
or (3.10)
u1(l 3;) +1426?) (3.11)

34

Thus substituting for u we get

_E_1__A1 (111:)de E__2A2 21 dry d__Wj
[’(2 :“id—fll >3“ ., _)+d. 2 /, (:1uid—x—‘H2u ,- .7...)

(3.12)
—Pu4 +3.1 (V — u2) +A2(V — N3)

Now taking partial derivative with respect to u1,u2,u3,u4,V,Al,A2 and equating to

zero, we will get a FE model that has 7 equations in seven variables looking like

 

 

 

 

 

 

{914—1 4,4 0 0 00 MM) (f1)
—§lfh £11111 0 0 0—10 u2 0
o 0 55214232142004 1.3 0
o 0 £2142 E2{£200 o u4 = P (3.13)
0 0 0 o 011 v o
o —1 0 0 10 0 A1 0
(0 0 —1 0 10 OJUZ) (0)

Clearly the above system is symmetric, banded and not positive definite. Solving we

get,

f1 =P (3.14)
A1 =P (3.15)
A2 = —P (3.16)
Pl
u =u =V=-—— (3.17)
2 3 EIAI
Pl Pl

 

= 3.18
“4 E141+52Az ( )

Now we do the same thing using two penalty parameters to impose the constraints

35

_EA dry dwj EA uidw dw
l_l_2fol(i2 ldx121__x)d +22____2/21(i§ 1(211

(3.19)
"P”4+ le-(V”“2)2+ YfW -u3)2

Partially differentiating the above equation with respect to u1,u2,u3,u4 and V we get

5 equations and the FE model looks like

(’51? £19 0 0 0\l0)lf1\
-§l{11§1,‘f‘—1+yl 0 0 —y1 u2 o

 

 

 

 

 

 

O O EzléeryZ -§2142 —yz u3 = 0 (3-20)
0 0 £2142 Ezf‘l 0 u4 P
\ 0 -Y1 —Y2 0 Y1+Y2l KVl k“)

The above matrix is symmetric, banded, sparse and also positive definite after applying

the boundary conditions. Solving, we get

f1 =—P (3.21)
Pl
V — <— 1 + I )P (3 23)
Y1 E141 °
1 1 1
u = ——+—+— P (3.24)
3 (Er/‘1 Y1 Y2)
l 1 l 1
= __ ._ p .2
“4 (5141+Y1) +(152/‘2+\l2)P (3 5)

And as 71 and 72 tend to infinity the values u3, V and u4 take values obtained using the
Lagrange multiplier method. Thus the Lagrange multiplier method sets the limit of accu-
racy obtained using the penalty method. Dividing the solutions obtained from the Lagrange

multiplier and the penalty method for 144 we get

36

enalty _1_ _1_
“ii 11 + Y2

ulagrange

 

l (3.26)

l +
E131 E242

Now if y is made dependent on the model using the following relations, Y1 = B (—11-1E A )
penalty

and 72 = B (@2142), we get, 512—875}? 2 1+ 11;. As we see, for a high value of B the solu-
4

tion approaches that given by the Lagrange multiplier method.

3.3 Interface Element to model Delamination

The above-described Interface element can be used to model delamination. One feature
of the interface element that facilitates delamination is the ability to release desired portion
of the interface element smaller than the finite element itself. The general model based
on which it is developed is the bilinear-softening model (Figure 3.2). Here the stresses
increase to the strength of the material, so that it is proportional to the penalty parameter.
Then the stiffness of the interface is reduced by decreasing the penalty parameter in the

form, V = (1 — D)y , here D is a damage parameter and has a value between 0 and 1.

In Model delamination, the stress increases to strength of the material, until the relative
displacement between the top and bottom finite element mesh is 80 . Then the interface
accumulates damage and when the relative displacement between the meshes is 5f, the in—
terface is completely debonded. As in the figure, when the relative displacements between
the two FE meshes are between and 50 and 8f, if there is unloading, it unloads to the ori-

gin. This is obtained by using , D, the damage parameter in the following manner

0 2 (1 -— D)y5 (3.27)

37

   
   

     

.4

 

 

Sf’g

om».---
o \

Figure 3.2: Bilinear Model

The value of D is computed from geometry as

= 5f(5 - 50)
8(5f —80)

Thus if one knows any two of the four interfacial properties, Gc,ot,5f,50 the interface

0(8) (3.28)

element can be described completely.

8 0
CC = 12—! (3.29)
0t
5 = —— (3.30
0 Y >

Where, Gc is the strain energy release rate, or is the ultimate strength of the interface,
50 is the relative displacement at the onset of failure and 5 f is the relative displacement
when there is complete debonding of the surfaces. This bilinear interface model is applied

to a sub-region of the interface element. This enables simulating delamination by releasing

38

>Q

Glc

 

b--------------

 

20 820 81f 82f

Figure 3.3: Updating tensile stresses

the desired portions of the interface element.

In general during an axial crush of a tube both Mode I and Mode II delamination are

seen to take place. In order to simulate this , a quadratic failure criterion is used as follows.

ems—321
2
(:72) + (52:6)2 = 1 (3.32)

Where, 8Z0 : le,8x0 Z %,O’Z = YzOz,sz = YxOx
Based on this criterion updated values of T and S are calculated. Thus the interlayer
strengths and the relative displacements at which they occur are updated for the mixed

mode failure. Figures 3.3 and 3.4 illustrate this.

39

)0

 

 

S_
S'— , E
E Glc
' >8
6:10 8xo Std 81d :1
Figure 3.4: Updating shear stresses
Thus,
I I
5' = “5ch (3.34)

Similarly the quadratic failure criterion,

“)2 (3)2
_ + _ :1 (3.35)
(Glc 02c

is assumed for the strain energy release rates and the values of 61 c and 62c is updated for

the model.

As mentioned in the various literature reviews it is also important to model friction be-

tween the interfaces. This is taken care of by ABAQUS itself. The model is constructed

40

Nodes with same co-ordinates

  
   
 
 
  

Interface Element

Axisymmetric Elements

Free surfaces

Figure 3.5: Interface Element enables to model friction between plies

in such a way that the layers of elements are joined by the interface elements. At the in-
terface, the elements of one mesh are in reality not joined to the other, though the nodal
co-ordinates are the same, they are nodes of different elements. In ABAQUS this allows us
to model the two layers as contacting surfaces and hence include friction. This is illustrated

in Figure 3.5.

3.3.1 Specification of the Interface Element for Input File

ABAQUS allows the specification of user element using the following keywords

o *USER ELEMENT

41

o *ELEMENT
o *UEL PROPERTY

The data card for the keywords their options and specification for the model is described

below

*USER ELEMENT, NODE=12, TYPE=U1, PROPERTIES=10,
IPROPERTIES=6,COORDINATES=2, VARIABLES=4

1, 2
This keyword specifies the use of a user element that is a FORTRAN program written for

ABAQUS/STANDARD.

Options

0 NODE=12: Specifies the total number of nodes that form a single interface element.
This is actually the summation of N, M, NPS values of the UEL PROPERTY key-
word. This may vary from one interface element to another but must be unique for a

particular TYPE of interface element.

a TYPE=U1: Specifies the type of the interface element. Each type may have a unique
property definition. Subsequent types may be given names as U2, U3 and so on. If
all the interface elements on the interface have the same property, as is the case with

the analysis performed for this work, they all may be of TYPE=U1.

42

PROPERTIES=10: This specifies total number of properties that are required for the

interface element. The default of 10 must be used.

IPROPERTIES=6: This indicates the total number of integer properties required.

This is a default value for the program and must not be changed.

COORDINATES=2: Total number of coordinates for a node of an interface element.
Since this is a 2D interface element it has a value 2. This is a default and cannot be

changed.

VARIABLES=4: Total number of solution dependent state variables that needs to be

stored by the program. This is also a default value and must not be changed.

DATA CARD: The data card 1,2 indicates that two degrees of freedom are active for

the nodes of the interface element.

*ELEMENT, TYPE=U1, ELSET=INTFACE1

3000,423,424,425,426,15001,]5002, 15003, 15004,634,635,636,63 7

3001,426,42 7,428,429, I 5004, I 5005, I 5006,] 5007,63 7, 638,63 9, 640

3 002,429,430, 43 I ....................................................

This keyword specifies the interface element number and the node numbers (already

defined) which form the interface, (basically the connectivity for the interface elements).

43

All the elements with a unique property may be assigned to an element set INTFACEI

(ELSET=INTFACE1 )

0 DATA CARD In the example illustrated 3000,3001 are interface element numbers.
423 through 426 are nodes on the left mesh from top to bottom of interface element
3000. 15001 through 15004 are nodes that form the interface element 3000. They
must be at least four in number and correspond to the value of NPS in the UEL
PROPERTY definition. 634 through 637 are nodes that form the nodes of the right
hand side mesh, from top to bottom. The specification for interface element 3001 is
now straightforward from the explanation and Figure 3.6. It must be noted that the
interface element nodes must coincide with finite element nodes at least at the ends,
the intermediate nodes may not. However for the sake of simplicity in this case the
nodes of the finite elements and interface elements coincide as indicated by Figure

3.6.

*UEL PROPERTY, ELSET=INTFACE1

**BETA, THKl, THK2, E1, E2, Srl, GIC, Sr2
0.6670, 0. 191112505, 0. 191112505,25E+9,25E+9,2OE+6, 134. ,23E+6
**GIIC, FRICTION, N, M, ND, NEP, NPS, IFDAMG

154.1, 0.0, 4,4, 9, 5, 4, 1

The keyword specifies that the property described in the data card for this keyword is ap-

plicable for the element set INTFACEl

—-—> Interface Element

 

634

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

3000 635 ___, Axrsymmetnc Element
636
637
3001 638
639
640
1 j
r
* Slave Surface
* Master Surface

Figure 3.6: Specification Of Interface element nodes

45

DATA CARD: The specification of values for this datacard is most crucial for the anal-

ysis and is explained in sufficient detail in the following paragraphs.

0 BETA, (B)
Based on the explanation given on delamination implementation and the formulation
of the interface element, it should be clear that the value of Beta (B) is related to the
penalty parameter. The value of B specified here will be used in the code to enforce

the constraint.

The derivation to specify B value for the analysis is as follows. The Lagrange multi-

plier can be related to the penalty parameter in the following manner.

A =yC(u1,u2) (3.36)

,where C (uluz) is the constraint.
The Lagrange multiplier can be physically understood to be the force that is required

to enforce the constraint.

Hence we have, (Figure 3.7)

d1
F zy/O (V—u)dl (3.37)

V is the primary degree of freedom for the interface element and u for the finite ele-
ment nodes. Now, if we assume that the crack propagates a length of dl for a crack

opening displacement of 60 at the maximum stress 0’“ , then

46

 

 

Figure 3.7: Illustration for derivation of y

F = §50d1 (3.38)

, dividing both sides by the area, we get,

Y = (2)(0u)(thk) (3.39)
50

 

The penalty parameter, 7, is related to beta (B) through the geometry of the structure

as follows.

 

r = B ((E)ghk)) (3.40)

E ,The Young’s modulus
d ,The length of a finite element along the interface
thk, Thickness of finite element mesh.

Equating the two values of y, we get

_ 20nd

B_ E50

 

(3.41)

Hence with appropriate values of ou,d,E,50 (material/model dependent) one can

specify the value of B required for the analysis. Since delamination is modeled for

47

both Model and ModeII any one of the values may be used to define the value of B.
In using axisymmetric elements the thickness value may be specified as the circum-

ference of a circle with radius at the interface element.

Ambiguity in y definition
From the Figures 3.2, 3.3 and 3.4 the values of y is given by:

y = a = Y1 (3.42)

Essentially the slope of the stress/relative displacement graph. But in this section we

have derived the relation,

 

Y: (2)(0u(08r)7)(’hk) =12 (3.43)
0
Clearly,
72 = (2)(thk)y1 = (constant)y1 (3.44)

Also, y is just a large number that enforces the constraint, so the specification of y is

consistent as long as it’s definition is made clear. 80 we can assume that the slope of

the Stress/displacement is w—dYsI'W = y].

TH K 1 and THK 2
These are the thickness of the left and right hand finite element meshes and have a

value 21tr, where r is the radius to the location of the interface element in the model.

E 1 and E2
These are the Young’s modulus of the left and right hand side meshes. For this anal-

ysis these values are the same.

48

0 SI'I, GIC
These are the ultimate tensile stresses and critical strain energy release rates for Mode

I delamination.

o Sr2, GIIC
These are the ultimate shear stresses and critical strain energy release rates for Mode

II delamination.

0 Friction
This specifies the coefficient of friction between the two surfaces. Since the code
is successful in modeling contact between the two surfaces only for cases of small
sliding, this is not used in the current analysis and a value of 0.0 is given. Instead the

two meshes are modeled separately and contact is established by ABAQUS directly.

0N

Number of nodes on the left hand mesh

0M

Number of nodes on the right hand mesh

0 NEP

Number of evaluation points, Must be five and is a default

0 NPS

Total number of nodes forming an interface element. This must have a minimum

49

value of four.

0 ND
Specifies the number of subdivisions in an interface element. This must have a value
greater than one. Greater the number of subdivision more accurately is the delami-
nation predicted but more expensive is the solution. This enables releasing portions

of the interface smaller than the finite element itself.

0 IFDAMG
Giving a value of l enforces the constraints and includes friction , a value of 0 simply

prevents overlapping of the two meshes.

50

CHAPTER 4

Analysis

In this chapter, details of the Finite Element Analysis are described. Results of a few
experiments that help to intuitively understand the effect of various parameters that affect

the energy absorbed during the crush experiment are also discussed.

4.1 The Model

The Finite Element Model was based on simple axisymmetric (CAX4) elements and
the initiator was modeled as an analytical rigid surface. Figure 4.1 illustrates the model.
All nodes at the top of the tube were encastered (fixed in all degrees of freedom) while
the initiator was given a displacement in the upward (or) 2 direction. Tube of length of
4 inches, outer diameter and a thickness of 2.5 inches and 0.075 inches respectively was

chosen for the study.
As specified in section 3.3 the interface elements naturally allow us to model friction

between the plies once delamination has taken place. Due to severe contact convergence

problems only two plies were used here to simulate delamination.

51

AXIS

  

   

Axisymmetric Elements

    

     

 

 

 

 

 

 

   

 

 

 

 

Initiator K,

Figure 4.1: The Model

 

 

L
.1

 

 

52

4.2 Concrete Material Properties

The following properties were chosen for simulating material failure in the concrete

material model.

B, Y ' M d 1 2E N
oungs ouus +10%:
u, Poisson’s ratio 0.3

Strain ratio, see section 2.1.4, Page 27 1.2 or 2.4

Ultimate Tensile stress, on 155E + 6A]:
m

Ultimate Compressive stress, am: 1555 + 6%
m

A visco—plastic parameter of IE — 4 was used to facilitate numerical convergence, avoid

mesh-dependency and prevent strain localization. Please see [17, 18, 21—23] for details.

4.3 Interface Element Properties

The interface element required to simulate delamination was given the following prop-

erties as in Table 4.1. (Section 3.3.1, Page 41)

It can be noted that two values are given for ModeII interface strength and four val-
ues of strain energy release rates are specified. These values are used to see the effect of
delamination on the energy absorbed. This method of specification of the factors enables
automatic ANOVA studies using softwares. Detailed explanation of these factors are given

in the following section.

53

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Beta, B 0.967619047
THK1,THK2 0.193506399m
E1,E2 2.0E+1031Nz
Srl 80E +63]?
GIC 352g
Sr2 20E +63% (or) 20E+7I—n"lz
GIIc 50% (or) 100% (or) 500% (or) 10000?"
FRICTION 0.0
N 4
M 4
ND 6
NEP 5
NPS 4
IFDAMG 1

 

 

 

 

Table 4.1: Interface Properties

4.4 Factors and Their Levels

For preliminary analysis the following factors were considered.
1. Radius of the initiator (% inch and 156 inch)
2. Friction between the initiator and the tube (0.0 and 0.3)
3. Material Strain Ratio (1.2 and 2.4)
4. Strength of the interface (20E + 6 N and 20E + 7 N )
r? m—2
5. Relative Displacement Ratio (5 and 10)
6. Friction between the plies (0.0 and 0.2)

Each factor has 2 levels and this amounts to performing 26 or 64 experiments to per-

form the analysis of variance. The contribution of a certain factor, in percentage is obtained

54

by taking the ratio of sum of squares for that factor to the total sum of squares [26]. The
next chapter deals with the details of the procedure to perform the ANOVA and obtain the

contributions.

The first two factors can be considered as model related factors, the third is material de-
pendent factor while the last three are related to delamination. The model and the friction
factors are very direct and are available in the ABAQUS input file to be changed directly.

The specification of other factors are explained below.

4.4.1 Material Strain Ratio

This factor is considered to study the sensitivity of the response to the material strain
energy release rates. Please see Figure 4.2. From section 2.1.4, page 27 we see that the

strain energy release rate depends on the strain ratio through the expression
G — lO (l)(SR)e (41)
c — 2 u f .

Therefore for a certain 8 f , Cu and I, each value of SR we will have a different strain energy

release rate value. Please see section 2.1.3, page 21 to implement this factor.)

In general it is observed that higher values of negative slope in the post failure region
of material models cause serious convergence problems. If the ANOVA study indicates
that the response is significantly sensitive to the strain energy release rate (Area under the

stress vs crack opening displacement diagram) then it would be possible to maintain the

55

 

 

 

>8
sf 1.213f 2.413f

Figure 4.2: Illustration of Strain Ratio Factor

strain energy release rate while lowering the value of the negative slope to obtain quicker

numerical results. Figure 4.3 illustrates this.

4.4.2 Strength Of the Interface

This factor has two levels with each level depicting the ModeII strength. The lower
strength causes delamination while the higher one does not. By including this as a factor,
the effect of delamination on the total energy absorbed can be studied. This method of

implementing this factor enables automation of the ANOVA study using softwares.
Experiments where delamination is not required to take place may also be run by modi-

fying the input file to not include the interface element. This will save sufficient computation

time and is used in the current work.

56

 

 

 

 

 

 

 

 

Gui Strength of the Material
0 - ...................... , ,
“1 uOi Crack Opening Displacement
Gc Strain Energy Release Rate
OIu2' u02 combination gives quicker solutions
0112 -__ .........

 

 

“or “02

I _ 1 _
3oulu01 —— Guzuoz -Gc

Figure 4.3: Purpose of Strain Ratio Factor

57

 

 

60 580 1080 ’5

Figure 4.4: Illustration Of Relative Displacement Ratio factor

4.4.3 Relative Displacement Ratio (RDR)

The purpose of this factor is to study the effect of ModeII strain energy release rate on
the total energy absorbed. The values of the Relative Displacement Ratio refers to the
ratio of the displacement between the layers at failure to the displacement at complete loss
of interface strength (26) (Figure 4.4). In ModeII these refer to the sliding displacement.

The following expressions together with Figure 4.5 illustrate this factor

1
GIIc = %(ou)(RDR)50 (4.3)

It is expected that the conclusions from ANOVA for this effect will be helpful in the

same sense as the Strain Ratio factor. The RDR factor is also implemented by changing the

58

 

 

 

 

 

 

 

 

oui Strength of the Interface
0 ... ..................... , , .
“1 Ofi Maximum Relative Displacement
GIIc Mode 11 Strain Energy Release Rate
Ou2 ,5t‘2 combination gives quicker solutions
5112 --- .........

 

P-—--——- —-—--——-—-—-—-— —-

 

8“ 5f2

lama“ =_10u,tsf2 =oiic
2 2

Figure 4.5: Purpose of Relative Displacement Ratio Factor

59

value of the Strain Energy release rate in the *UEL PROPERTY keyword. The strain en-

ergy release rate is dependent on the Model] strength and the RDR value. The 50 term can

be eliminated using a constant value for the penalty parameter. The following calculations

clarify this and are useful when using a software to automate the ANOVA study done for

this work.

We know from section 3.3.1, equation 3.39

_ (2)(THK)(0u)
Y" 8
0

 

Cu _ Y
S ‘ <2><THI<>

and from equation 3.41 ,

_ (2)(0u)(d) _ (0(4)

’3‘ (5)186) ‘" (E)(THK) :CONSTANT

 

Therefore B will be a constant for all the analysis.
Now from equation 4.3

0116 = %(ou)(RDR)50

(or)

(o£)(RDR)(THK)
Y

GIIC I

 

(4.4)

(4.5)

(4.6)

(4.7)

(4.8)

In this way, while automating ANOVA studies, the GIIC values can be made as a variable,

dependent on RDR and 0;, factors.

60

4.5 Calculations

In this section all the required calculations are shown to facilitate rerunning of analysis

for future research.

4.5.1 Calculation of B

From the previous sections we know B is a constant for all the analysis that will be

performed for the AN OVA study.
d = L
2 (Y)( mg) (4.9)
(E )(TI-I K )
Where,

d, is the length of a finite element along the interface
L, is the Length of the tube in meters (4” = 0.1016m)
Ne, is the number of elements along the length of the tube (210 for model used)

(2)(THK)(ou) __ 2 =1: 0.193506399 :1: 2E +7
80 " 1E - 6

 

 

= 7.74025596E + 12 (4.10)

7.74025596E + 12 . ‘21‘0'0'1016
_ 2E+ 1040193506399

 

= 0.967619047 (4.11)

4.5.2 Calculation of Material Strain Energy Release Rates

We know that,

 

1
GC,SR = iouiefXSRlal (4.12)
2
0u(SR)(l) _ 2.4025E+16*(1.2)*4.808E—4 _
25 - 2*25 +10 —346.5366 (4.13)
2
.4 , . _
Gu(SR)(l) = 2 0255+ 16* (2 4) :14 808E 4 : 693.0732 (414)

 

2E 2*2E-l-10

61

4.5.3 Calculation of Mode 11 Strain Energy Release Rates

From section 4.4.3

 

 

 

2
o RDR THK
G”c.Gu,RDR=( "X Y“ ) (4-15)
4E +14 * 5 *O.193506399
011 = = 50 4.16
c,2E+7,5 7.74025596E+12 ( )
4E +14 * 1040193506399
II = =1 4.
G C425+7110 7.74025596E+12 00 ( 17)

The following calculations are required only if analysis requiring no delamination are
performed with a high value for the interface strength (usefitl while automating ANOVA
through softwares). Analysis where no delamination is required can be easily and effi-

ciently run by modifiring the input file to not include the interface element.

4E+16*5*0.193506399 _

 

 

on = _ 5000 4.18
c.2E+8.5 7.74025596E + 12 ( )
4E +16=1= 10 a: 0.193506399
0,1635%,” — 7.74025596E+ 12 "‘ 10000 (4'19)

4.6 Preliminary Analysis

As indicated previously, due to a high degree of non-linearities, convergence was a se-
rious problem with all the analysis in this work. All the analyses did not run to completion
due to these non-linearities. Several analysis were performed to determine the robustness
of the finite element software to various values of these factors. Two such experiments that
enabled to determine final values of one factor and also helped in deciding the the total final

initiator displacement are discussed below.

62

4.6.1 EXPERIIVIENT 1

In this experiment two tubes of same dimensions were crushed using two initiators of
different radii (I36ind’ and %inch). The initiators were forced to cover a distance of 1.5

inches. Other properties used for this experiment are tabulated in Table 4.2

No delamination element was used for this experiment. The Energy absorbed for higher
radiussed and lower radiussed initiator was 319.38] and 192.32J respectively. This indi-
cates that the lower radiussed initiator increases the energy absorbed.The peak load is also
higher for the lower radiussed initiator. Please see Figure 4.6 and Figure 4.7 and Table 4.3

From the above results it can also be observed that the contribution of friction is approxi-
mately the same for both the experiments. This concludes that the interaction between the

radius effects and friction effects are negligible.

The time required for the analysis using the 1% inch radius initiator was as long as 2
days and around 20, 000 increments was needed. It hence was agreed that a 156 inch radius

initiator would be used for all the analysis in the ANOVA study.

 

E,Young’s Modulus 2E + 10 £2

u,Poisson’s ratio 0.3
Strain ratio 1.2
N

Ultimate failure stress,ou 155E + 6;?
Friction coefficient 0.3

 

 

 

Table 4.2: Values used for Experiment 1

63

 

 

 

 

 

 

 

 

L rote. v0 pathetic-i 1min 1116 meta
1x102)
‘0. I I l I l 1 I I l T INT I T I I
0.00 — —
A
z 0.00 — _
V
o «—<
l 4.00
2.00 —
.l
o.mIIIIIILIIII_IIIIII
0.00 0.00 10.00 24.00 32.00 [x101
Dlsplaoement (m)

 

 

 

Figure 4.6: Force vs Deflection with 1—36inch radiussed initiator

_
C

 

 

rote. v1 tit-placate: 1 sue inch xdiuIJ

[1110']
a.” I T T I I I I l’ i

 

5.00 -— #

4.00 r— A

son —

Force (N)

2.” —

1.00 a

 

 

 

o_w r I I 1 1 I; 1 L 14 L i I r
0.00 10.00 20.00 30.00 [x10‘]
Displacement (m)

 

 

 

Figure 4.7: Force vs Deflection with Tisinch radiussed initiator

Usual features of atypical force deflection diagram can be noticed in Figure 4.7. It can
be divided into three main zones.
Zone I: In this region there is a steady increase in the load as the tube is deformed. This
indicates that the tube is deforming elastically. The first drop in the load around 4KN cor-
responds to the beginning of material failure.
Zone II: This zone indicates a steady increase in the load along with failure taking place in
the tube. The zone directly leads into the sustained crushing zone (ZONE III).
Zone 111: In this zone sustained crushing is seen to be taking place, where the load oscil-

lates about a mean crushing load indicating progressive failure of the tube.

4.6.2 EXPERIMENT 2

In one of the experiments the 1% inch radius initiator was used and delamination was
allowed to take place with an interface strength of 2E + 7%, a strain ratio of 2.4 was used
in this case. No friction was used, between the plies and between the tube and initiator.
This analysis died with the displacement of the initiator being only 0.69375 inches. This
was the analysis with the lowest value for the final initiator displacement. After careful
inspection of the deformed configuration (Figure 4.8) and the ABAQUS message file it

was concluded that contact between the plies caused severe convergence problems. Since

 

 

 

 

Radius Total Energy Energy Absorbed due Percentage contribu-
Absorbed to Frictional Effects tion of fiiction

%inch 319.381 1451 45.4

%1nch 192.321 92.21 47.7

 

 

 

 

 

Table 4.3: Results for Experiment 1

65

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

:3"an 5. All futon cates: B at. A! low IWIII (bl
I C 1 m.rh-:t).odb Wt Landau! 6.3-1 Bod Oct 00 11-02 I? m 200)
l

I! I It '1
Inge”? USO- Itop 1'1.- - 0.5419

Dolor-d V02: 0 Dolor-um knlo Patton 01.009900.

 

 

 

Figure 4.8: Illustrating delamination

the delamination code used in this work does not have the capability to write required data
into a restart file, analysis involving delamination could not be restarted if they died due
to problems in numerical convergence. It is however true that responses used for ANOVA
must be for a certain final initiator displacement. It was therefore agreed that response up

to 0.69375 inches of initiator displacement would be used for the ANOVA study.

66

CHAPTER 5

Analysis of Variance

5.1 ANOVA

Analysis of variance (ANOVA) is a well known method in the field of design and anal-
ysis of experiments. The use of ANOVA here, is to know which factors contribute sig-
nificantly to the variability in the response. This is called Factor Screening. For sake of
continuity a brief description on the procedure to perform ANOVA is described in this chap-
ter. For details on ANOVA please see [26]. Finally results of ANOVA are tabulated and
are described in detail. A section describing future extension of this work for automated
analysis of variance is given with a motivation for incorporating ”localization limiters” in

material models for composite materials involving softening.

Let a response say Energy absorbed, during crushing of the tube be assumed to be
dependent on five factors such as, A, B, C, D, E. These factors can be qualitative (Pres-
ence/absence of delamination) or quantitative (numerical values). Let each of these factors
have only two levels each. This means that we want to study the response for 2 values of
each of these factors. This will give us a total of 2 >1: 2 =1: 2 =1: 2 =1: 2 = 25 = 32 experiments to
be performed and are referred to as a full factorial study. Apart from each of the factors
individually affecting the response (Main Eflect), there may be interactions (Interaction
effect) among these factors. If it is known in advance that interaction of certain factors do

not affect the response, then a subset of the 25 = 32 experiments could be performed and

67

this is called fractional factorial study.

If the experiments were performed in a laboratory with testing equipments, each of
these experiments would be performed more than once in order to negate the effect of
random errors that creep in due to several uncontrollable factors. However, performing
experiments on a computer give rise to deterministic responses. Hence each of the 32 ex-

periments need to be performed only once.

The basic idea of Analysis of variance stems from dividing the total variability into
its component parts. The total corrected sum of squares (SST) is taken as a measure of
variability. It is the sum of squares of the differences between individual treatment values
and the grand averages. The total sum of squares is taken as the sum of squares of the
differences between treatment averages and the grand averages (SSTreatment) plus the
sum of the squares of the differences of observations within treatments from the treatment
average (SS E). Clearly SSTreatment is a measure of differences between treatment means,
while SSE can be only due to random error. Since computer experiments are deterministic,

this term vanishes. So we have:

SST : SSTreatment (5'1)

The SST = SSTreatment are equal to the sum of squares due to each ejfect (including

Main eflect and interaction effects). This gives:

ssT = SSTreatment = SSA +553 +ssC +ssD +555 +ssAB ........... +ssABCDE (5.2)

68

 

Bl B2

 

C1 C2 C1 C2

 

D1 D2 D1 D2 D1 D2 D1 D2
El E2 E1 E2 E1 E2 E1 E2 E1 E2 E1 E2 E1 E2 E1 E2

 

 

Al ( 1 ) e d de c ce cd cde b be bd bde bc bce bcd bcde

 

A2 0 ae ad ade ac ace acd acde ab abe abd abde abc abce abcd abcde

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1, 2 indicate low and higher levels of the factor

Figure 5.1: Labeling the responses

Here, SS A,SS B,SSC,SS D,SS E are the sum of squares of the main effects while SS A 3....
etc, are the sum of squares of the interaction effects. Thus the contribution of each effect to
the total variability is simply found by taking the ratio of the sum of squares for that effect

to the total sum of squares.

ss
ABCE (5.3)

ContributionABCE 2 SS
T

Calculation of the sum of squares for each effect is done using contrasts for each effect

in the following manner.

Each response is labeled as a, b, c, d, e, ab, ...... abde, ..... abcde. Please see Figure 5.1.

The presence of an alphabet indicates a higher level for that factor and absence of that

69

alphabet indicates a lower level of that factor. For example, abce indicates the response
that was obtained for the experiment that used high levels of factor A, B, C, E and a low

level of factor D. A simple method to calculate the contrast is given by,
ContrastABCD = (a —1)(b — 1)(c— 1)(d— 1)(e+1) (5.4)

= abcde+cde+bde+ade+bce+ace+abe+e+abcd+cd+bd+ad+bc+ac+ab+ (l)

—a—b—c-abc—d——-abd—acd—bcd—ae—be—ce—abce—de—abde-—acde—bcde
(5.5)

Note: (I ), indicates response, when all factors are at a low level

Once the contrasts are estimated, the effects and the sum of squares are estimated as fol-

lows:

I
SABCD = 2—4(C0ntrastABCD) (5.6)

1
ssABCD = 53(c0n1ras1ABCD)2 (5.7)

5.2 Results and Discussion

The various factors considered to study the sensitivity of the response are tabulated in
Table 5.1. ANOVA was performed using method described in the previous sections.

Low values, of factor E is assumed, as higher values did not give satisfactory responses
interms of final initiator displacement. The analysis failed to converge for higher values
for this factor and restart was not possible due to reasons specified in section 4.6.2, Page
65. Even with the current values, each analysis completed to different extents. It was
hence decided to use responses for only 0.69375 inches of the final initiator displacement

(Experiment 2, previous chapter) for the ANOVA. Due to this reason the force deflection

70

 

 

 

 

 

 

 

Factor Label Factor Description Levels Values
A Friction between in1- 2 0.0 and 0.3
tiator and tube
B Strain Ratio 2 1.2 and 2.4
C Delamination 2 Presence and absence
D Relative Displacement 2 5 and 10
Ratio
E Friction between plies 2 0.0 and 0.2

 

 

 

 

 

Table 5.1: Factors and their levels

diagrams shown henceforth may not show sustained crushing or (ZONE III). It can also
be noted that factors D and E have relevance only if delamination takes place, hence only

20 experiments need to be performed.

The responses for each experiment and the contribution to the variability from each

factor is given in Tables 5.2 and 5.3.

It is found that the contribution from the interaction effects to the total variability of the
response are negligible. The Main Effects of factors A, B, C have significant contribution

while the other factors have no contribution at all.

From Table 5.3 it is evident that about 61.97% of the variability in the response can
be attributed to friction between the tube and the initiator, 10.39% to the strain ratio factor
and 25.55% to delamination. It is observed that the presence of delamination decreases
the total energy absorbed and and has a significant contribution to the variability in the
response. It can be hence concluded that numerical modeling of delamination is crucial for

accurate represention of the failure of composite tubes. In this work it was observed that

71

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Experiment Response(Joules)
(l) 25.67
a 48.79
b 35.54
ab 62.5
c 40.76
ac 74.21
be 51.55
abc 88.2
d 26.5
ad 49.9
bd 36.2
abd 62.8
cd 40.76
acd 74.21
bcd 51.55
abcd 88.2
c 27.2
ae 50.1
be 39.7
abe 64.3
ce 40.76
ace 74.21
bce 51.55
abce 88.2
de 27.1
ade 51.1
bde 37.2
abde 65.8
cde 40.76
acde 74.21
bcde 51.55
abcde 88.2

 

 

Table 5.2: Response for Each Experiment

delamination takes place during the early stages of the crush thus decreasing the overall

stiffness of the tube causing low peak forces and hence low energy absorption. Figures 5.2

72

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Factors Contribution (%)
A 61 .97
B 10.39
C 25.55
D 0
E 0.06
AB 0. 18
AC 1.73
AD 0.01
AB 0
BC 0
BD 0
BE 0.01
CD 0
CE 0.06
DE 0
ABC 0
ACD 0.01
ADE 0.01
AEB 0
ABC 0
BCD 0
BDE 0
BBC 0
CDE 0
BDA 0
ABCD 0
ACDE 0.01
BCDE 0
CABE 0
DBAE 0
ABCDE O

 

 

 

 

 

Table 5.3: Contributions of each factor

and 5.3 illustrates this.

73

 

 

rozco " Deflection to: tarpon-until” I

 

 

[1110'
20$ I I l I T I I I l

 

 

 

 

I. -.
1.50 -— —
A j. -(
Z
V
g 1 00 — —
a ’- .I
0.50 _
0.” I I m L 1 I I I 1 I
0.00 4.00 0.00 1200 10.00 20.00 [1110‘]
Displacement (m)

 

 

Figure 5.2: Force vs Deflection for Experiment ((1))

 

 

force VI Deflection tor Dwarf-eta]

 

 

 

1.50 >— —

Force (N)

LW 1- —l

0.50 -

 

 

0.“) 1 I 1 I l 1 I 1
0.00 5.00 10.00 15.00 20.00 [1110‘]

Displacement (m)

 

 

 

Figure 5.3: Force vs Deflection for Experiment (c)

74

The friction factor (A) plays a vital role and positively increases the energy absorbed.
This can be observed from Figures 5.4 and 5.5. Fiction between the plies factor (E) does
not affect the response significantly and may be neglected in this case. However, it is im-
portant to increase the number of plies to know the effect of this factor in greater detail.

Material strain energy release rate ( factor B) describes about 10.39% of the variability

—

 

 

 

 

 

 

 

[ Pore. Va Dolloctlon (or Wort-nut]
[1110']
_ 'T ' l—' 1* 1 r 1 W
5.00 —- —
4.00 — —
A
Z r -
V
g 3.00 i- —
o - _
“- 2.00 —
1.00 _
0.” I I L I l I PLJ 1 1 I 1 1 I
0.00 0.00 10.00 24.00 32.00 [1110‘]
Displacement (m)

 

 

 

Figure 5.4: Force vs Deflection for Experiment(ac), High friction coefficient

in the response and plays a vital role. This suggest that a higher material strain energy
release rates increases the energy absorbed by the tube during failure. The sensitivity of
the response to this factor enables to improve the convergence of the analysis, as specified
in section 4.4.1, Page 55. This phenomenon is demonstrated in the next section. Delami-
nation strain energy release rate (D) does not affect the response at all and hence suggests

that it can be neglected. However, it may be important to consider more plies for a more

75

—
f

 

 

Porto V. Deflection to: Input-mun I

010']
3.20

 

fiI‘I‘ITT'I'I‘
2.”
2.40
2.00
1.”

TWFITIIIIII

Force (N)

1 .20
0.”
0.40

llllllLllLJllJl

 

 

a.” 1 I l l l l 1 l l l 1 J 4 L 1
0.00 0.00 1000 24.00 32.00 [x10‘]
Displacement (m)

 

 

 

 

Figure 5.5: Force vs Deflection for Experiment (c), Low friction coefi‘icient

convincing conclusion on this factor.

5.2.1 Effect of Strain Ratio

It was mentioned in 4.4.1, page 55 that decreasing the negative slope, in post failure
region of the stress—strain diagram can cause improved convergence, if the ANOVA results
indicate high sensitivity to the Strain Ratio factor. From Table 5.3 it can be seen that there
is 10.39% of contribution through the strain ratio factor to the variability in the energy ab-

sorbed.

To verify the effect of strain ratio, three experiments were conducted. For all three
experiments strain energy release rate of 3465366 was considered. The strain at complete

loss of strength was considered at three different strain ratios (1.2, 2.4 and 4.8). Please see

76

Figure 5.6.

 

 

 

o
A
01
I
I
I
I
I
I
O2 — ‘ i
I I
I I
I I
I I
O’ I I
3— ! l‘~
I I
I I I
I l
I I I
I I I
I I I
I I I
I I I
I I I
i I i ‘\ ~e
8122 831.281 2.4:2 4.81::3

1 1.251 102.41 104.851 0
302 1 =32 I:’2.=?3 3: c

Figure 5.6: Effect of Strain Ratio

For a constant strain energy release rate the modified strength was found using the fol-

lowing derivation.

l
_ 1 “f

 

_ /ZGCE
O'f — [(SR) (5.10)

For the three different strain ratios, the strengths, energy absorbed and number of incre-
ments required to complete the analysis for 1.25 inches of final initiator displacement are

tabulated in Table 5.4. It is observed that maintaining the strain energy release rate while

77

lowering the negative slope in the post failure region improves convergence. An increase
in the Strain Ratio by a factor of 2, improves convergence by nearly 50% (in terms of the
number of increments) while the decrease in the energy absorbed is 13%. This indicates
that a careful choice of the strain ratio, will help improve convergence without compromis-

ing on the results.

 

 

 

 

SR Strength, in Pa Response, (J) No. of Increments Time to complete
1.2 155E + 6 84 5833 7.7 Hrs

2.4 10960113 + 6 73 2837 2.5 Hrs

4.8 77.5E+6 51 1998 1Hr40min

 

 

 

 

 

 

 

Table 5.4: Effect of Strain Ratio

Figures 5.7 through 5.12 show the deformed shapes and the force deflection diagram
of all the three analyses. Large values of the strain ratio, cause strain localization as can
be observed in Figure 5.9, at regularly spaced intervals. This strain localization manifests
itself as huge drops in the load deflection diagram (Please see Figure 5.12). Such a behav-
ior does not represent the response curve correctly and hence places a limit on the choice
of the Strain Ratio value. The cause of this strain localization is attributed to the reduced
viscoplastic effect. In other words, the introduction of rate dependency by viscoplastic
regularization of the constitutive equations acts as ”localization limiters”. Quicker conver—
gence due to high Strain Ratio negates the viscoplastic effects and causes strain localization.

Please see [21—23] for further details.

78

geld-1t 5. A“ laser.
- mrthM-I.

8! It -1
nights? 50)): I!" “H.- u 1.
outer—d Va: 0 Detox-:10. mu

 

Figure 5.7: Deformed shape , Strain Ratio = 1.2

 

 

 

 

i .-
I .
I.
on
.-
== ,.,
II o 0 Q
E: ..:.:“‘
II 9 “
" ’0‘“.
3: 9"“
u g'.‘
I
E. '6“-
:: '-'I:'
.1: 3":
II -:.I
:1 «3'
I: g":
:0 '03!
.- :Q.l
I. . .. C
33: I“:
on 0..
=III 0...:
'3: :0..-
’:--I ...0 .
u“. 0.0 0..
‘5‘». 0%....
\“ O 0 0 0 O
C. O .. . O O
‘0‘:o.e. 0’ 0.0.0.
‘9’:’o"- e“’¢"¢‘0’
’~'~?0-'a':3=3l‘ 3’9
.~:.::.:===“

All lactate Io- end In col-inacl
“06"

I (M t '1 1:),
I ooi: nu_:l2 cab n/sl-ndud G.) I 'No law
I

S: s at 1
mgr-u? 20)? It" ‘21.- I 1,000
Dolor-d vu: U Deter—uln- Scnlc furor: CLOOOooO

8
Or-

ueod thh Ie e coat I
1.46:0: 2:1 :00:

 

Figure 5.8: Deformed shape, Strain Ratio 2 2.4

79

 

 

 

 

 

 

s
’ II
V ' ~' ' . . O ‘..t
0 I “l!_” w W's/It“. ' 3 'fl .‘ 1 :00)
L—I

'1 \ it -1
biz—n? 1990. In. T!- I 1.000

Dolor-d vu- I: Dolor-ue- ken "not: 0).. .000

 

 

Figure 5.9: Deformed shape, Strain Ratio = 4.8

 

 

 

000']

w t l T T f] 1 I I l I I r
+— -l
2” t _.
2.40 - j
A 2.00 — —
E, - .
g 1.00 — —
u. 1.20 - _
0.00 i
0.40 —

 

 

o-m 11L141111111111
0.00 0.00 10.00 24.00 3200 [x10‘1

Dleplaoement (m)

 

 

 

Figure 5.10: Force vs Displacement, Strain Ratio = 1.2

80

—-—- luau: Milo 2.0

000’]
2.”._Irrlrlrrr [[1

 

2.40

2.00

I r I I j
l

1.00

1.20

Force (N)
. Tfi i
1

0M —

0.40 _

 

 

0.” 1 I 1 l 1 l 4 J 1 l 1 l 1 l 1
0.00 aoo 1000 24.00 32.00 0110‘]
Displacement (m)

 

 

 

 

Figure 5.11: Force vs Displacement, Strain Ratio = 2.4

_

 

(x10‘
“1’ i'iFT‘i‘i'l"
2.00~

I.” —

1.20— -

Force (N)

OM— #

0.40 .2

 

 

°.m L 1 1 i 1 l 1 l 1 l L i 1 l 1
0.00 0.00 10.00 24.00 32.00 [x10fl
Displacement (m)

 

 

 

 

Figure 5.12: Force vs Displacement, Strain Ratio = 4.8

81

In Table 5.5, error in the Energy absorbed (compared to Energy Absorbed for a Strain

Ratio of 1.2), and the improvement in number of increments required to complete the anal-

ysis is tabulated for several values of Strain Ratio. Figure 5.13 is a graphical representation

of the same data. For a Strain Ratio 1.3, the error in the Energy absorbed is 8.43% while

there is a 35% reduction in the number of increments required to complete the analysis and

a 22% reduction in time required to complete the analysis.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

SR Response, No. of In- Hrs. to Strength, Error % Im- % Im-
(J) crements complete MPa in Re- prove- prove-
sponse, ment in ment in
% N0. of time
Inc.

1.2 84 5833 7.7 155 0 0 0

1.3 76.92 3846 6.00 148.919 8.43 34.06 22.08

1.4 75.19 3318 5.5 143.502 10.49 43.12 28.57

1.5 76.31 3189 4.25 138.6362 9.16 45.33 44.81

1.6 74.22 2793 3.17 134.233 11.64 52.12 58.87

2.0 76.86 3012 4.5 120.062 8.5 48.36 41.56

2.4 73 2837 2.5 109.655 13.10 51.36 67.53

2.8 70 2521 3.42 101.47 16.67 56.78 55.63

3.2 68.33 2567 3.42 94.9177 18.65 55.99 55.63

3.6 62.48 2130 3.25 89.489 25.62 63.48 57.79

4.0 62.16 2128 3.00 84.896 26.0 63.52 61.04

4.4 63.94 1376 2.0 80.946 23.88 76.41 74.03

4.8 51 1998 1.66 77.5 39.29 65.75 78.44

 

Table 5.5: Effect of Strain Ratio, (Detail)

It is true that number of increments is not the best representation for the convergence

rate nor is it directly proportional to the time required to complete the analysis. This is

because each increment may have different number of iterations depending on whether the

guessed solution is within the convergence radius. The trend that has been established in

82

 

 

 

 

 

 

80 T I I I I I I
.k
’ \
’ \
70'- ’ \ -4
/ \
I w
+ — — u-K
/
60 u— / u
/
/+- — — *+/
/
4‘ xf/ /
so~ / ‘ \ , , —
I \4"
a 4-
x
§4o~ ’ a» a
9 /
:2 I l/
1-
I
30" ' / _
i /
i /Q- -- - “k .. \ I
w
' /
/
20" I — ‘1' -4
l I A“ _.
1' ,r I
I
it "' w
" 4— Error in Energy Absorbed
" --+— Improvement in convergence
L l g L 1 I I

 

0 i
1 1 .5 2 2.5 3 3.5 4 4.5 5
Strain Ratio

Figure 5.13: Improvement in convergence through high Strain Ratio

this study only indicates that it is possible to obtain quicker converged solution with a lower

sacrifice in the accuracy of the response (Energy Absorbed), with a careful choice of the

Strain Ratio.

5.3 Conclusion

As a result of the work done, it is evident that obtaining responses for quasi-static crush
problems is very challenging due to high degree of non-linearity. However , ”localization
limiters” like visco plastic algorithm come to our rescue and help to obtain meaningful so-
lutions with some effort. Convergence due to contact may be overcome by having a finer

mesh for the slave surface than the master surface (please see Figure 3.3.1 .

83

 

 

ANOVA results enable us to conclude that delamination is a significant factor affecting
the response. It decreases the overall stiffness of the tube and decreases the energy that can
be absorbed by the tube. Numerical modeling of delamination is crucial to accurately rep-
resent failure of composite tubes during crushing. Factors like friction contribute to nearly
50% of the total energy absorbed and this is consistent with the literature [2] and [3]. Other
factors like friction between plies and the RDR factor do not affect the response in this

case. Use of several plies may lead to a different conclusions for these factors.

5.4 Scope for future Research

As described, ANOVA was done using Ms Excel and all the analysis needed to be car-
ried out separately from the command prompt. It is also possible to automate the process
using softwares. This was attempted in this work using a software called HEEDS (Heirar-
chial Evolutionary Engineering Design System) developed at Red Cedar Technology, East
Lansing, Michigan, for advanced optimization problems. HEEDS is also capable of per-
forming AN OVA once the variables (Factors) and their levels are specified. HEEDS allows
the user to tag the design variables and the response. The user can also specify a tool or
many tools that will be used to evaluate the design and output the response. Once all such
evaluations are carried out HEEDS writes out in an output file the contribution of each ef-
fect. All factors except delamination (presence and absence) are quantitative and can be
directly tagged, to set up the problem in HEEDS. The delamination factor can be imple-
mented using the methods described in sections 4.4.2 and 4.5.3, by varying the strength
of the interface strength to effect delamination. This work may be extended to perform

advanced ANOVA studies using HEEDS.

84

Though convergence and mesh sensitivity are crucial problems in softening material
models in numerical analysis , methods to overcome these can be implemented in mate-
rial models used to simulate failure in composite structures. Some methods available for
this purpose are suggested and explained in [18, 21—23]. As a motivation for furthering
improvement in building material models involving softening it should be interesting and
useful to incorporate concepts of ”localization limiters” and objective methods of specify-

ing the softening behavior.

85

BIBLIOGRAPHY

[1] A penalty based interface technology for connecting independently modeled substruc-
tures for simulation growth of delamination in composite structures, by Dr. Antonio
Pantano and Dr. Ronald C. Averill.

[2] Mamalis, A. G., Manolakos, D. B., Demosthenous, G. A., Ioannidis, M. B., Crash—
worthiness of composite Thin-Walled structural Components. Technomic Publishing
Co., INC. Lancaster Basel.

[3] Mamalis, A. G., Manolakos, D. E., Demosthenous, G. A., Ioannidis, M. B., The
Static and Dynamic Axial Crumbling of thin walled fiberglass composite square
tubes. Composites Part B, 28B(1997), 439—45 1.

[4] Abdel-Haq, M and Newaz, G. M., Role of Failure Modes on Energy absorption in
Unidirectional PMC Tubes. Journal of Composite Materials, Vol. 35, No.11/2001.

[5] Yuan, Y. B., Vlegelahn, G. L., Modeling of Crushing Behavior of Fiber-
glass/Vinylester Tubes. Advanced Composite Materials: New Developments and
Applications Conference Proceedings, Detroit, Michigan USA Sept. 30 - Oct. 3, 1991.

[6] Gary L.Farley and Robert M. Jones, Analogy for the effect of Material and Geo-
metrical Variables on Energy-Absorption Capability of Composite tubes. Journal of
Composite Materials, Vol. 26, No. 1/ 1992.

[7] Scott Burr, Daryl Hertema, Jeff Sweeney, Taking a hit:- Modeling the parameters of
thermoplastics in car crashes. Plastics Engineering, September ’97.

[8] Mamalis, A. G., Monolakos, D. E., Demosthenous, G. A., Ioannidis, M. B., Axial
Collapse of thin-walled fiberglass tubular composite components at elevated strain
rates. Composites Engineering, 3, 1994.

[9] Mamalis, A. G., Yuan, Y. B., Viegelahn, G. I., Collapse of thin walled composite sec-
tions subjected to high speed axial loading. International Journal of Vehicle Design,
1992, 13, 564.

86

[10] Hannapel, F. W., Yuan, Y. B., Vlegelahn, G. L., Martin. C. D., Crashworthy char-
acteristics of fiberglass/polyester frusta under quasi-static axial loading. Advanced

Composite Materials: New Developments and Applications Conference Proceedings,
Detroit, Michigan USA, Sept. 30-Oct. 3, 1991.

[11] Shu Ching Quek, Anthony M.Waas, Jennifer Hoffman, Venkatesh Agaram, The
Crushing Response of braided and CSM glass reinforced composite tubes. Composite
Structures, Vol. 52, 2001, 103-112.

[12] David C. Flemming, Modeling delamination growth in composites using MSC
Dytran.

[13] Shawn J. Beard, Energy absorption of braided composite tubes, A Dissertation sub-
mitted to the department of aeronautics and astronautics and the committee on grad-
uate studies of Stanford University, April 2001.

[14] Ronald Krueger , Pierre J. Minuguet and T. Kevin O’Brien, A method for calculating
Strain Energy Release Rates in preliminary design of composite skin/stringer

debonding under multi-axial loading. Composite Structures: Theory and Practice,
ASTM STP 1383, 2000, pp. 105-128.

[15] Gary L.Farley and Robert M. Jones, Prediction of energy absorption capability of
composite tubes. Journal of Composite Materials, Vol, 26, No. 3/1992.

[16] E.Haug & A. De Rouvray, Crash Response of Composite Structures. Engineering
Systems International SA, 94578 Rungis-Cedex, France

[17] Don Simons, Constitutive Model V & V for softening Materials.

[18] Zdenek P. Basant and Jaime Planas, Fracture and Size effect in concrete and other
Quasi-brittle Materials. CRC Press

[19] State of the Art Report on Finite element Analysis of Reinforced Concrete. Published
by American Society of Civil engineers.

[20] Hibbit, Karlsson, & Sorenson, ABAQUS V 6.3 Manuals.

[21] Etse, G, William, K., Failure Analysis of Elasto—Viscoplastic Material Models. ASCE
Journal of Engineering Mechanics, Vol 125, 1999, 60—69.

87

[22] Milan Jirasek, Zdenek P. Bazant., Inelastic Analysis of Structures. John Wiley & Sons;
(November 12, 2001) .

[23] Needleman, A., Material rate dependence and mesh sensitivity in localization prob-
lems. Comp. Meth. Appl. Mech. and Engineering, Vol 63, 1988, 69—85.

[24] Hillerborg. A., Modeer. M., Peterson, P. E., Analysis of crack formation and crack
growth in concrete by means of fracture mechanics and Finite Elements. Cement and
Concrete Research, Vol. 6, 1976, pp. 773-782.

[25] J. N. Reddy, An introduction to the Finite Element Method-2nd ed. McGraw-Hill
series in Mechanical Engineering.

[26] Douglas C. Montgomery, Design and Analysis of experiments. John Wiley & Sons,
INC.

[27] T. Kevin 0’ Brein. Characterization, Analysis and Prediction of Delamination in
Composites using fracture Mechanics.0ral Reference number: ICF1009420R

[28] J. N. Reddy, Energy and Variational methods in Applied Mechanics. A Wiley-
Interscience Publication, John VVrley & Sons.

[29] Dr. Robert M. Jones, Mechanics of Composite Materials. Taylor & Francis

[30] Alastair F. Jhonson and Anthony K. Pickett, Impact and Crash modeling of Compos-
ite Structures: A challenge for Damage Mechanics.

[31] P. Ladeveze and E Le Dantec, Damage Modeling of the Elementary ply for laminated
composites. Composites Science and Technology, Vol. 43, No. 3/1992, pp 257-267.

[32] Kupfer, H., Hilsdorf, H.K, Rusch, K., Behavior of concrete under biaxial stresses.
Journal of the American Concrete Institute, Vol. 66, No.8, August 1969, pp. 656-666

[33] Chen, W. F., Plasticity in Reinforced Concrete. McGraw Hill, (January 1982).

[34] Sultan T. Mesmar, Roger S. Crouch, Fracture Regularization by extension of the
explicit microplane model for concrete to include Duvaut-Lions visco-plasticity by,
European Congress on Computational Methods in Applied Seiences and Engineering,

88

ECOMAS 2000.

[35] Computational Inelasticity by Simo and Hughes. Springer Verlag, (August 1998).

[36] Anderson, T. L., Fracture Mechanics, Fundamentals and Applications, CRC Press.

[37] Leslie M. Moore and Bonnie K. Ray, Statistical methods for sensitivity and perfor-
mance analysis in computer experiments. Proceedings of the 1999 Winter Simulation
Conference.

[38] Jerome Sacks ,Susannah B. Schiller and VVIlliam J. Welch, Design of Computer Ex-
periments. Technometrics, Feb. 1989, Vol. 31, No. l.

[39] “William J. Welch, Robert J. Buck, Jerome Sacks, Henry P. Wynn, Toby J. Mitchell and

Max D. Morris, Screening, predicting and Computer Experiments. Technometrics,
February 1992, Vol. 34, No. 1.

[40] Lubliner, J., J. Oliver, S. Oller and Ofiate, A Plastic- Damage model for concrete.
International Journal of Solids and Structures ,Vol. 25, No. 3, pp 229-326, 1989.

[41] Lee, J., and G. L. Fenves, Plastic Damage Model for Cyclic loading of Concrete
Structures. Journal Of Engineering Mechanics, Vol. 124, No. 8, pp, 892-900, 1998

[42] Tim Rzenitzek, Dr. Heiner Miillersch'o'n, Dr. Frank C. Giinther and Michal Wozniak,
Two-stage Stochastic and Deterministic Optimization.

[43] MI. Y, Crisfield, M. A., Davies, G. A. O., and Hellweg, H. B., Progressive Delamina-
tion Using Interface Elements. Journal of Composite Materials, Vol. 32, No. 14/1998.

[44] Xianqiang Lu and Dahsin Liu, Finite Element Analysis of Strain Energy Release Rate
at Delamination Front. Journal of Reinforced Plastics and Composites, Vol. 10, May
1991.

[45] Bruno. D, Greco. F, Lonetti. P, Modeling of delamination and contact in composite
laminates.

[46] Reedy. Jr, E. D., Mello F. J., Guess. T. R., Modeling the initiation and growth of
delaminations in composite structures. Journal of composite materials, Vol 31, No.
8, 1997, pp. 812 - 831.

89

[47] Hull, D, A unified approach to progressive crushing of fiber-reinforced composite
tubes. Composite Science and technology, Vol. 40, (1991), 377-421.

[48] Vreede, P. T., Van Veldhuizen, Kragtwijk, S. P., Carleer, B. B, Finite Element anal-
ysis of tubular energy absorber in automobiles. University of Twente, Enschede, The
Netherlands.

[49] Ribeaux, M, Warrior, N. A., Effect of Impact Damage on the Specific Energy Absorp-
tion of Glass/Polyster Composites. University of Nottingham, NG7 2RD UK.

90

 

S

H11111111111111W