• TABLE OF CONTENTS
HIDE
 Title Page
 Acknowledgement
 Table of Contents
 Abstract
 Introduction
 General theories of nonlinear...
 Three-dimensional link element
 Linear shell element
 Nonlinear shell element
 Nonlinear shell element perfor...
 Conclusions and recommendation...
 Appendix
 Reference
 Supplemental bibliography
 Biographical sketch
 Copyright














Title: Nonlinear gap and Mindlin shell elements for the analysis of concrete structures
CITATION THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00090188/00001
 Material Information
Title: Nonlinear gap and Mindlin shell elements for the analysis of concrete structures
Series Title: Nonlinear gap and Mindlin shell elements for the analysis of concrete structures
Physical Description: Book
Creator: Ahn, Kookjoon,
 Record Information
Bibliographic ID: UF00090188
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: alephbibnum - 001583924
oclc - 23011869

Table of Contents
    Title Page
        Page i
    Acknowledgement
        Page ii
    Table of Contents
        Page iii
        Page iv
    Abstract
        Page v
        Page vi
    Introduction
        Page 1
        Page 2
        Page 3
        Page 4
        Page 5
        Page 6
        Page 7
        Page 8
        Page 9
        Page 10
        Page 11
        Page 12
    General theories of nonlinear analysis
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
    Three-dimensional link element
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
    Linear shell element
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
        Page 83
        Page 84
        Page 85
        Page 86
        Page 87
        Page 88
        Page 89
        Page 90
        Page 91
    Nonlinear shell element
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
        Page 102
        Page 103
        Page 104
        Page 105
        Page 106
        Page 107
        Page 108
        Page 109
        Page 110
        Page 111
        Page 112
        Page 113
        Page 114
        Page 115
        Page 116
        Page 117
        Page 118
        Page 119
        Page 120
        Page 121
        Page 122
        Page 123
        Page 124
        Page 125
    Nonlinear shell element performance
        Page 126
        Page 127
        Page 128
        Page 129
        Page 130
        Page 131
        Page 132
        Page 133
        Page 134
        Page 135
        Page 136
        Page 137
        Page 138
        Page 139
        Page 140
        Page 141
        Page 142
    Conclusions and recommendations
        Page 143
        Page 144
        Page 145
    Appendix
        Page 146
        Page 147
        Page 148
        Page 149
        Page 150
        Page 151
        Page 152
        Page 153
        Page 154
        Page 155
        Page 156
        Page 157
        Page 158
        Page 159
        Page 160
        Page 161
        Page 162
        Page 163
        Page 164
        Page 165
        Page 166
        Page 167
        Page 168
        Page 169
        Page 170
        Page 171
        Page 172
        Page 173
        Page 174
        Page 175
        Page 176
        Page 177
        Page 178
        Page 179
        Page 180
        Page 181
        Page 182
        Page 183
        Page 184
        Page 185
        Page 186
        Page 187
        Page 188
        Page 189
        Page 190
        Page 191
        Page 192
        Page 193
        Page 194
        Page 195
        Page 196
        Page 197
        Page 198
        Page 199
        Page 200
        Page 201
        Page 202
        Page 203
        Page 204
        Page 205
        Page 206
        Page 207
        Page 208
        Page 209
        Page 210
        Page 211
        Page 212
        Page 213
        Page 214
        Page 215
        Page 216
        Page 217
        Page 218
        Page 219
        Page 220
        Page 221
        Page 222
        Page 223
        Page 224
        Page 225
        Page 226
        Page 227
        Page 228
        Page 229
    Reference
        Page 230
        Page 231
        Page 232
        Page 233
        Page 234
        Page 235
        Page 236
        Page 237
    Supplemental bibliography
        Page 238
        Page 239
    Biographical sketch
        Page 240
        Page 241
        Page 242
        Page 243
    Copyright
        Copyright
Full Text












NONLINEAR GAP AND MINDLIN SHELL ELEMENTS
FOR THE ANALYSIS OF CONCRETE STRUCTURES


















By

KOOKJOON AHN


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


UNIVERSITY OF FLORIDA


U::i: L "At; 7 OF FLORIN LIfRPMR17


1990














ACKNOWLEDGEMENTS


I would like to express my deep gratitude to professor

Marc I. Hoit for his invaluable guidance and support. I also

thank professors Clifford 0. Hays, Jr., Fernando E. Fagundo,

John M. Lybas, and Paul W. Chun for being on my committee. I

also express my gratitude to professor Duane S. Ellifritt

for his help as my academic advisor at the start of my Ph.D.

program.

Thanks are also due to all the other professors not

mentioned above and my fellow graduate students, Alfredo,

Jose, Mohan, Prasan, Prasit, Shiv, Tom, Vinax, and Yuh-yi.

Finally, I am thankful to every member of my family,

especially my wife and son, for their patience and support

in one way or another.

The work presented in this dissertation was partially

sponsored by the Florida Department of Transportation.















TABLE OF CONTENTS

page

ACKNOWLEDGEMENTS ...................................... ii

ABSTRACT ...................................... ........ v

CHAPTERS

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

1.1 General Remarks ............................ 1
1.2 Link Element ............................... 2
1.3 Shell Element .............................. 5
1.4 Literature Review .......................... 5

2 GENERAL THEORIES OF NONLINEAR ANALYSIS .......... 13

2.1 Introduction ............................... 13
2.2 Motion of a Continuum ...................... 14
2.3 Principle of Virtual Work .................. 16
2.4 Updated Lagrangian Formulation ............. 18
2.5 Total Lagrangian Formulation ............... 22
2.6 Linearization of Equilibrium Equation ...... 26
2.7 Strain-Displacement relationship
Using von Karman Assumptions ............... 28

3 THREE-DIMENSIONAL LINK ELEMENT .................. 34

3.1 Element Description ...................... 34
3.2 Formation of Element Stiffness ............. 43
3.3 Solution Strategy .......................... 51
3.4 Element Verification ....................... 52

4 LINEAR SHELL ELEMENT ............................ 59

4.1 Introduction ............................... 59
4.2 Formulation of Shape Functions ............. 59
4.3 The Inverse of Jacobian Matrix ............. 64
4.4 Membrane Element ........................... 66
4.5 Plate Bending Element .................... 73

5 NONLINEAR SHELL ELEMENT ........................ 92

5.1 Introduction ............. ................. 92
5.2 Element Formulation ........................ 93


iii











5.3 Finite Element Discretization .............. 100
5.4 Derivation of Element Stiffness ............ 113
5.5 Calculation of Element Stiffness Matrix .... 115
5.6 Element Stress Recovery .................... 119
5.7 Internal Resisting Force Recovery ......... 122

6 NONLINEAR SHELL ELEMENT PERFORMANCE ............ 126

6.1 Introduction ............................... 126
6.2 Large Rotation of a Cantilever ............ 126
6.3 Square Plates .............................. 133

7 CONCLUSIONS AND RECOMMENDATIONS.................. 143

APPENDICES

A IMPLEMENTATION OF LINK ELEMENT .................. 146

B IMPLEMENTATION OF LINEAR SHELL ELEMENT ......... 170

C IMPLEMENTATION OF NONLINEAR SHELL ELEMENT ....... 219

REFERENCES ............................................ 230

SUPPLEMENTAL BIBLIOGRAPHY ............................ 238

BIOGRAPHICAL SKETCH ................................... 240














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

NONLINEAR GAP AND MINDLIN SHELL ELEMENTS
FOR THE ANALYSIS OF CONCRETE STRUCTURES

By

KOOKJOON AHN

August, 1990

Chairman: Marc I. Hoit
Major Department: Civil Engineering

Segmental post-tensioned concrete box girders with

shear keys have been used for medium to long span bridge

structures due to ease of fabrication and shorter duration

construction.

Current design methods are predominantly based on

linear elastic analysis with empirical constitutive laws

which do not properly quantify the nonlinear effects, and

are likely to provide a distorted view of the factor of

safety.

Two finite elements have been developed that render a

rational analysis of a structural system. The link element

is a two-dimensional friction gap element. It allows opening

and closing between the faces of the element, controlled by

the normal forces. The Mindlin flat shell element is a

combination of membrane element and Mindlin plate element.

This element considers the shear responses along the element









thickness direction. The shell element is used to model the

segment itself.

The link element is used to model dry joints and has

shown realistic element behavior. It opens under tension and

closes under compression. The link element has shown some

convergence problems and exhibited a cyclic behavior.

The linear Mindlin shell element to model the concrete

section of the hollow girder showed an excellent response

within its small displacement assumption.

The nonlinear Mindlin flat shell element has been

developed from the linear element to predict large

displacement and initial stress (geometric) nonlinearities.

The total Lagrangian formulation was used for the

description of motion. The incremental-iterative solution

strategy was used. It showed satisfactory results within the

limitation of moderate rotation.

Three areas of further studies are recommended. The

first is the special treatment of finite rotation which is

not a tensorial quantity. The second is the displacement

dependent loadings commonly used for shell elements. The

third is the material nonlinearity of concrete which is

essential to provide realistic structural response for safe

and cost effective designs.














CHAPTER 1
INTRODUCTION



1.1 General Remarks



In the past few decades segmental post-tensioned

concrete box girders have been used for medium to long span

bridge structures. Highway aesthetics through long spans,

economy due to ease of fabrication, shorter construction

duration are some of the many advantages of precast segment

bridge construction.

The segments are hollow box sections, match cast with

shear keys in a casting yard, then assembled in place,

leaving the joints entirely dry. The shear keys are meant to

transfer service level shears and to help in alignment

during erection.

Current design methods are heavily based on linear

elastic analysis with empirically derived constitutive laws

assuming homogeneous, isotropic materials. The behavior

under load of the bridge system is very complex. Analyses

which do not properly quantify the nonlinear effects

including the opening of joints in flexure, are likely to

provide a distorted view of the factor of safety existing in

a structural system between service loads and failure. The










potential sliding and separation at the joints due to shear,

and by deformations generated by temperature gradients over

the depth and width of the box further complicate the

problem [1].

Two finite elements have been developed that render a

rational analysis of the system. The link element is a two-

dimensional friction gap element. It allows sliding between

the faces of the element, controlled by a friction

coefficient and the normal forces. It also accounts for zero

stiffness in tension and a very high stiffness under

compression. This link element was borrowed from rock

mechanics and newly applied to this problem to model the dry

joint between the segments. The Mindlin flat shell element

is a combination of membrane element and Mindlin plate

element. This element considers the shear responses along

the element thickness direction. The shell element was used

to model the segment itself. This element can handle large

displacement and geometric nonlinearities.



1.2 Link Element



A link element is a nonlinear friction gap element used

to model discontinuous behavior in solid mechanics. Some

examples are interfaces between dissimilar materials and

joints, fractures in the material, and planes of weakness.

These have been modeled using constraint equations, discrete










springs and a quasi-continuum of small thickness [2]. The

following characteristics of prototype joints were

considered.


1. Joints can be represented as flat planes.

2. They offer high resistance to compression in the
normal direction but may deform somewhat modeling
compressible filling material or crushable
irregularities.

3. They have essentially no resistance to a net tension
force in the normal direction.

4. The shear strength of joints is frictional. Small
shear displacements probably occur as shear stress
builds up below the yield shear stress.


A model for the mechanics of jointed rocks was

developed by Goodman [3]. The finite element approximation

was done as a decomposition of the total potential energy of

a body into the sum of potential energies of all component

bodies. Therefore, element stiffness is derived in terms of

energy.

The Goodman element was tested for several modeled

cases.


1. Sliding of a joint with a tooth.

2. Intersection of joints.

3. Tunnel in a system of staggered blocks.


A problem with the Goodman's two dimensional model is

that adjacent elements can penetrate into each other.

Zienkiewicz et al. [4] advocate the use of continuous










isoparametric elements with a simple nonlinear material

property for shear and normal stresses, assuming uniform

strain in the thickness direction. Numerical difficulties

may arise from ill conditioning of the stiffness matrix due

to very large off-diagonal terms or very small diagonal

terms which are generated by these elements in certain

cases. A discrete finite element for joints was introduced

which avoids such theoretical difficulties and yet is able

to represent a wide range of joint properties, including

positive and negative dilatency (expansion and compaction

accompanying shear) [3].

The element uses relative displacements as the

independent degrees of freedom. The displacement degrees of

freedom of one side of the slip surface are transformed into

the relative displacements between the two sides of the slip

surface. This element has been incorporated into a general

finite element computer program [5]. The use of relative

displacement as an independent degree of freedom to avoid

numerical sensitivity is discussed in detail [6]. An

isoparametric formulation is given by Beer [2]. A four-node,

two-dimensional link element and a eight-node plate bending

element were used to model the dry jointed concrete box

girder bridge with shear keys [7].










1.3 Shell Element



The shell element is formulated through the combination

of two different elements, the membrane element and the

Mindlin plate bending element.

The Mindlin plate element is different from the

Kirchhoff plate element in that the former allows transverse

shear deformation while the latter does not.

The nonlinearities included in the formulation of the

flat shell element is for large displacement and geometric

nonlinearity due to initial stress effects. The large

displacement effects are caused by finite transverse

displacements. These effects are taken into account by

coupling transverse displacement and membrane displacements.

The initial stress effects are caused by the actual stresses

at the start of each iteration. These stresses change the

element stiffness for the subsequent iteration. These

effects are evaluated directly from the stresses at the

start of each iteration and included in the element

stiffness.



1.4 Literature Review



The purpose of nonlinear analysis is to develop the

capability for determining the nonlinear load-deflection

behavior of the structures up to failure so that a proper










evaluation of structural safety can be assured. There are

two general approaches for nonlinear analysis. The first

approach is a linearized incremental formulation by reducing

the analysis to a sequence of linear solutions. The second

approach is mathematical iterative techniques applied to the

governing nonlinear equations [8].

The advantage of the incremental approach results from

the simplicity and generality of the incremental equations

written in matrix form. Such equations are readily

programmed in general form for computer solutions [9].

A generalized incremental equilibrium equation for

nonlinear analysis can be found in [10, 11, 12]. The

formulation is valid for both geometrical and material

nonlinearities, large displacements and rotations,

conservative and displacement dependent (nonconservative)

loads.

There are two frames for the description of motion. The

difference lies in the coordinate systems in which the

motion is described. These are the total Lagrangian

formulation which refers to the initial configuration [10,

11] and the updated Lagrangian formulation which refers to

the deformed configuration [12]. There have evolved two

types of notations in the description of motion. A

correlation is given these two notations, the B-notations

and the N-notations, currently used in the Lagrangian

formulation of geometrically nonlinear analysis [13]. A










short history of early theoretical development of nonlinear

analysis can be found in [9, 14].

One form of updated Lagrangian formulation is the

corotational stretch theory [15].

Shell elements are often derived from governing

equations based on a classical shell theory. Starting from

the field equations of the three-dimensional theory, various

assumptions lead to a shell theory. This reduction from

three to two dimensions is combined with an analytical

integration over the thickness and is in many cases

performed on arbitrary geometry. Static and kinematic

resultants are used. These are referred to as classical

shell elements. Alternatively, one can obtain shell elements

by modifying a continuum element to comply with shell

assumptions without resorting to a shell theory. These are

known as degenerated shell elements. This approach was

originally introduced by Ahmad, Irons, and Zienciewicz [16,

17]. Other applications can be found in [8, 18-25].

In large rotation analysis, the major problems arise

from the verification of the kinematic assumptions. The

displacement representation contains the unknown rotations

of the normal in the arguments of trigonometric functions.

Thus additional nonlinearity occurs. Further difficulties

enter through the incremental procedure. Rotations are not

tensorial variables, therefore, they cannot be summed up in

an arbitrary manner [17]. One of the special treatment of










finite rotation is that the rotation of the coordinate

system is assumed to be accomplished by two successive

rotations, an out-of-plane rotation followed by an in-plane

rotation using updated Lagrangian formulation [26, 27].

Usually the loadings are assumed to be conservative,

i.e., they are assumed not to change as the structure

deforms. One of the well known exceptions is pressure

loading which can be classified as conservative loading or a

nonconservative loading [28]. Another is the concentrated

loading that follows the deformed structure. For example, a

tip loading on a cantilever beam will change its direction

as the deformation gets larger. As loading is a vector

quantity, the change in direction means that the loading is

not conservative. Sometimes this is called a follower

loading.

The governing equation for large strain analysis can be

used for small increments of strain and large increments of

rotations [29]. This can be regarded as a generalization of

nonlinearity of small strain with large displacement. If

large strain nonlinearity is employed, an important question

is which constitutive equation should be used [9].

The degree of continuity of finite element refers to

the order of partial differential of displacements with

respect to its coordinate system. Order zero means

displacement itself must be continuous over the connected

elements. Order one means that the first order differential










of displacement must be continuous. Thus the higher order

the continuity requirement, the higher the order of assumed

displacement (shape, interpolation) function.

Mindlin-Reissner elements require only Co continuity,

so that much lower order shape functions can be used,

whereas in Kirchhoff-Love type elements, high order shape

functions must be used to satisfy the C1 continuity.

Furthermore, since Mindlin-Reissner elements account for

transverse shear, these elements can be used for a much

larger range of shell thickness. The relaxed continuity

requirements which permit the use of isoparametric mapping

techniques gives good computational efficiency if formulated

in the form of resultant stresses [30]. Unlike compressible

continuum elements, which are quite insensitive to the order

of the quadrature rule, curved Co shell elements require

very precisely designed integration scheme. Too many

integration points result in locking phenomena, while using

an insufficient number of quadrature points results in rank

deficiency or spurious modes [30]. While Gauss point stress

results are very accurate for shallow and deep, regular and

distorted meshes, the nodal stresses of the quadratic

isoparametric Mindlin shell element are in great error

because of the reduced integration scheme which is necessary

to avoid locking [31].

The degenerate solid shell element based on the

conventional assumed displacement method suffers from the










locking effect as shell thickness becomes small due to the

condition of zero inplane strain and zero transverse shear

strain. Element free of locking for linear shell analysis

using the formulation based on the Hellinger-Reissner

principle with independent strain as variables in addition

to displacement is presented in [32].

Shear locking is the locking phenomenon associated with

the development of spurious transverse shear strain.

Membrane locking is the locking phenomenon associated with

the development of non-zero membrane strain under a state of

constant curvature. Machine locking is the locking

phenomenon associated with the different order of dependence

of the flexural and real transverse shear strain energies on

the element thickness ratio, and it is therefore strictly

related to the machine finite word length [33]. Some of the

solutions are as follows:


1. Assumed strain stabilization procedure using the Hu-
Washizu or Hellinger-Reissner variational principles
[33].

2. The assumed strain or mixed interpolation approach [34,
35].

3. Suppressing shear with assumed stress/strain field in a
hybrid/mixed formulation [30]. Suppression of zero
energy deformation mode using assumed stress finite
element [36].

4. Coupled use of reduced integration and nonconforming
modes in quadratic Mindlin plate element [37].

5. Higher order shallow shell element, with 17 to 25 nodes
[38, 39].









6. Global spurious mode filtering [40].

7. Artificial stiffening of element to eliminating zero
energy mode, special stabilizing element [41].

In the faceted elements, due to the faceted

approximation of the shell surface, coupling between the

membrane and the flexural actions is excluded within each

individual element, the coupling is, however, achieved in

the global model through the local to global coordinate

transformation for the elements [39].

In geometrically nonlinear analysis with flat plate

elements, it is common to use the von Karman assumptions

when evaluating the strain-displacement relations. The

assumption invoked is that the derivatives of the inplane

displacements can be considered to be small and hence their

quadratic variations neglected. However, this simplification

of the nonlinear strain-displacement relationship of the

plate, when used in conjunction with the total Lagrangian

approach, implies that the resulting formulation is valid

only when the rotation of the element from its initial

configuration is moderate. Thus for the total Lagrangian

approach to handle large rotations, simplifications of the

kinematic relationship using the von Karman assumptions is

not permitted [39].

Some of the special solution strategies to pass the

limit point are given in references [25, 42-48]. A limit

point is characterized by the magnitude of tangential









stiffness. It is zero or infinite at a limit point. Thus

conventional solution strategies fail at the limit point.

Arc length method was introduced in reference [42], and

applied in the case of cracking of concrete [43]. This was

improved with line search and accelerations in references

[44, 45]. Line search means the calculation of an optimum

scalar step length parameter which scales the standard

iterative vector. This can be applied to load and

displacement control and arc length methods [44].

The traditional solution strategies are iterative

solutions, for example, Newton-Raphson, constant stiffness,

initial stiffness, constant displacement iteration, load

increment [46] along with Cholesky algorithm with shifts for

the eigensolution of symmetric matrices [47] for element

testing for spurious displacement mode.

The vector iteration method without forming tangent

stiffness for the postbuckling analysis of spatial

structures is also noted [48].

The linearized incremental formulation in total

Lagrangian description has been used for this study of large

displacement nonlinearity including initial stress effects.

The special treatment of finite rotation is not included in

the current study. Material nonlinearity is also excluded.













CHAPTER 2
GENERAL THEORIES OF NONLINEAR ANALYSIS



2.1 Introduction



The incremental formulations of motion in this chapter

closely follow the paper by Bathe, Ramm, and Wilson [11].

Other references are also available [9, 10, 12, 14, 15, 49,

50, 51].

Using the principle of virtual work, the incremental

finite element formulations for nonlinear analysis can be

derived. Time steps are used as load steps for static

nonlinear analysis. The general formulations include large

displacements, large strains and material nonlinearities.

Basically, two different approaches have been pursued

in incremental nonlinear finite element analysis. In the

first, Updated Lagrangian Formulation, static and kinematic

variables, i.e., forces, stresses, displacements, and

strains, are referred to an updated deformed configuration

in each load step. In the second, Total Lagrangian

Formulation, static and kinematic variables are referred to

the initial undeformed configuration.

It is noted that using either of two formulations

should give the same results because they are based on the










same continuum mechanics principles including all nonlinear

effects. Therefore, the question of which formulation should

be used merely depends on the relative numerical

effectiveness of the methods.



2.2 Motion of a Continuum



Consider the motion of a body in a Cartesian coordinate

system as shown in Fig. 2-1. The body assumes the

equilibrium positions at the discrete time points 0, dt,

2dt, ..., where dt is an increment in time. Assume that the

solution for the static and kinematic variables for all time

steps from time 0 to time t, inclusive, have been solved,

and that the solution for time t+dt is required next.

The superscript on left hand side of a variable shows

the time at which the variable is measured, while the

subscript on left hand side of a variable indicates the

reference configuration to which the variable is measured.

Thus the coordinates describing the configuration of the

body using index notation are


At time 0 = x

At time t = txi

At time t+dt = t+dtx














P t+dt
P( x )


t
P( Xi)


0
P( Xi)


Fig. 2-1 Motion of a Body









The total displacements of the body are

At time 0 = ui

At time t = tui

At time t+dt = t+dtu

The configurations are denoted as

At time 0 = C

At time t = tc

At time t+dt = t+dtc

Thus, the updated coordinates at time t and time t+dt are
txi = Oxi + tui

t+dtx = Oxi + t+dtu

The unknown incremental displacements from time t to

time t+dt are denoted as (Note that there is no superscript

at left hand side.)

u = t+dtui tui (2.1)



2.3 Principle of Virtual Work


Since the solution for the configuration at time t+dt

is required, the principle of virtual work is applied to the

equilibrium configuration at time t+dt. This means all the

variables are those at time t+dt and are measured in the

configuration at time t+dt and all the integration are

performed over the area or volume in the configuration at

time t+dt. Then the internal virtual work (IVW) by the









corresponding virtual strain due to virtual displacement in
t+dtC is


S t+dt t+dt t+dt
t+dt ij t+dt eij (+dt dV) (2.2)


where,
t+dt rj = Stresses at time t+dt measured in the
t+dt j configuration at time t+dt.

= Cauchy stresses.

= True stresses.
t+dt
t+dt eij = Cauchy's infinitesimal(linear) strain tensor
referred to the configuration at time t+dt.

= Virtual strain tensor.

6 = Delta operator for variation.

and the external virtual work (EVW) by surface tractions and

body forces is


EVW = t+dt t ] 6 [ t+dt u (t+dt dA)
t+dt k t+dt k (t+dt



St+dt t+dt b 6 t+dt u (t+dt dV)
St+dt t+dt k t+dt uk

(2.3)
where,
t+dt tk = Surface traction at time t+dt measured in
t+dt the configuration at time t+dt.

t+dt Uk = Total displacement at time t+dt measured in
t+dt the configuration at time t+dt.









6 t+dt Uk = Variation in total displacement at time
t+dt t+dt measured in configuration at time t+dt

= Virtual displacement.
t+dt p = Mass density per unit volume.
t+dt

t+dt bk = Body force per unit mass.
t+dt k


and all the integration is performed over the area and the

volume at time t+dt.


2.4 Updated Lagrangian Formulation


In this formulation all the variables in Eqs. (2.2) and

(2.3) are referred to the updated configuration of the body,

i.e, the configuration at time t. The equilibrium position

at time t+dt is sought for the unknown incremental

displacements from time t to t+dt.

The internal virtual work, the volume integral in Eq.

(2.2) measured in the configuration at time t+dt can be

transformed to the volume integral measured in the

configuration at time t in a similar manner that is given in

reference [52]


IVW t+dt t+dt e (t+dt dV)
t+dt ij t+dt i (t dV)



t Sij 6 Et ij (t dV) = EVW (2.4)
I ~t t~d










where,
t+dt Sj = Second Piola-Kirchhoff (PK-II) stresses
t measured in the configuration at time t.

6 t+dt ej = Variations in Green-Lagrange (GL) strain
t tensor measured in the configuration at
time t.


The PK-II stress tensor at time t+dt, measured in the

configuration at time t can be decomposed as

t+dt t t
Sij =t Sij + t Sij j + t Sij (2.5)



because the second PK-II stress at time t measured in the

configuration at time t is the Cauchy stress.

From Eq. (2.1), the total displacements at time t+dt

measured in the configuration at time t is

t+dt t
t Ui = t Ui + Ui = t ui (2.6)


This is true because the displacement at time t measured in

the configuration at time t is zero. In other words, the

displacement at time t+dt with respect to the configuration

at time t is the incremental displacement itself.

And the GL strain is defined in terms of displacement as


Eij = 1 (Ui,j + Uji + Uk,iUk,j) (2.7)

E and U are used in the places of e and u to avoid confusion

between general strain and incremental strain, and between









general displacement and incremental displacement used in
this formulation. It is noted that these finite strain

components involve only linear and quadratic terms in the

components of the displacement gradient. This is the

complete finite strain tensor and not a second order

approximation to it. Thus this is completely general for any

three-dimensional continuum [52].

Then the GL strain tensor at time t+dt measured at time

t can be calculated as


+dt = i [(ui+ i)+ u) + ( j + tuj),i

+{(uk + tuk),i)( tuk + uk),j)]


= tuij + tuj,i+ tuki tuk,j]


= teij + tij

= tij (2.8)

where,

ttij = tei + t ij

= Incremental GL strain in tC.

tej = (ui,j + tuji )

= Linear portion of incremental GL strain in t.

= This is linear in terms of unknown incremental
displacement.
= Linearized incremental GL strain in tc.










ttij = t (uk,i tuk,j)

= Nonlinear portion of incremental GL strain.

The variations in Green-Lagrange strain tensor at time
t+dt measured in the configuration at time t can be shown as
using Eq. (2.8).

t+dt t
6 t ij = 6 ( tij + Ei ) = 6 t (2.9)


6teij = 0 because teij is known. There is no variation in
known quantity.

Then using the Eqs. (2.5), (2.8) and (2.9), the

integrand of Eq. (2.4) becomes

t+dt t+dt t
t Sij 6 qt =ij ( t i + tSij ) 6 tCij

= trij + tSij)(6 teij + 6tij)

=tij(6eij + 6t 1ij) + t~ij 6 te + rij 6 tij
t t
=tij 6 ti ij + ti 6eJ 6 (2.10)
tj t t ttij 6 t tii j` t


The constitutive relation between incremental PK-II

stresses and GL strains are


tSij = tCijkj t'kl


(2.11)








Finally the equilibrium Eq. (2.4) from the principle of
virtual work using Eqs. (2.10) and (2.11) is


I tCijkl tkl 6 tcij tdv

+ tri 6 t Vij tdV


= EVW J ij 6 eij tdV
t 1t
(2.12)
where, the external virtual work must be transformed from
t+dtc to tC. This is not applicable to conservative loading,

i.e., loading that is not changed during deformation.


EVW= t+dt tk [ t+dt uk (tdA)


+ t+dt t+dt b t+dt td) (2 )



and this is the general nonlinear incremental equilibrium
equation of updated Lagrangian formulation.


2.5 Total Lacrangian Formulation


Total Lagrangian formulation is almost identical with
the updated Lagrangian formulation. All the static and
kinematic variables in Eqs. (2.2) and (2.3) are referred to









the initial undeformed configuration of the body, i.e, the

configuration at time 0. The terms in the linearized strain

are also slightly different from those of updated Lagrangian

formulation.

The volume integral in Eq. (2.2) measured in the

configuration at time t+dt can be transformed to the volume

integral measured in the configuration at time 0 as [52]


t+dt 'i 6 t+dt eij t+dt
t+dt ]i t+dt


o t Sij t+dt j (dV) (2.14)



where,
t+dtS = Second Piola-Kirchhoff stress tensor
0 measured in the configuration at time 0.

6 t+dtij = Variations in Green-Lagrange (GL) strain
0 tensor measured in the configuration at
time 0.

The PK-II stress tensor at time t+dt, measured in the

configuration at time 0 can be decomposed as

t+dtsij= tsij+ *ij (2.15)
o1 ) o) 1 o01j (2.15)

From Eq. (2.1), the total displacements at time t+dt

measured in the configuration at time 0 is










t+dtu = u + Ui (2.16)
o 0i o0 o


Then the GL strain tensor at time t+dt measured at time
0 can be calculated as
t+dt = [(tui + oui)j + (tuj + ou i



=tt+
+(( e*k+ oUk),i)(( Uk+ Uk) j)]

= Eij + e + i

= ij + Eij (2.17)

where,
tt tutu tu tu)
oij = I (t u + + )
0 o o rJ o k,i k,j

= GL strain at time t in oC.

ij = oeij + oij

= Incremental GL strain in oC.

oeij = (oi, + oUj1 + uk,i ouk,kj+ k,j ouk,i)
= Linear portion of incremental GL strain in oC.
=This is linear in terms of unknown incremental
displacement.
= Linearized incremental GL strain in oC.

oij = i (oUk, ouk,j)

= Nonlinear portion of incremental GL strain.

The variations in Green-Lagrange strain tensor at time
t+dt measured in the configuration at time t can be shown as









using Eq. (2.17).

Stdtij = 6 ( ij+ ij ) = 6 (2.18)
o i= oec olu olt


6 te, = 0 because ej is known. There is no variation in
o ij o-J
known quantity.

Then using the Eqs. (2.15), (2.17) and (2.18), the

integrand of Eq. (2.14) becomes

t+dt t+dt ts
t+dt*i 6 t+dti ( tSi + S ) 6 e
o o o o0J o'J

= (ij oSij) oeij + 6 oiij)

= Sij(6 eij + 6 ij) + oSij 6 e 6 +

= Si 6 eij + tS i 6 ei + ij 6 o7ij (2.19)


The constitutive relation between incremental PK-II

stresses and GL strains are


Sij = oCijkj ockl (2.20)


Finally the equilibrium Eq. (2.14) from the principle


of virtual work using Eqs. (2.19) and (2.20) is











I oCijkl oekl 6 j dV


+ tsi 6 oij OdV


= EVW tij 6 e dV
o 1 oe ij
(2.21)

where, the external virtual work must be transformed from
t+dtc to oC. This is not applicable to conservative loading,

that is, loading that is not changed during deformation.


EVW = t+dt t t+dt uk (dA)



+ t+dt p t+dt b t+dt uk (odV) (2.22)



and this is the general nonlinear incremental equilibrium

equation of total Lagrangian formulation.


2.6 Linearization of Equilibrium Equation


The incremental strain from time t to t+dt is assumed

to be linear, i.e., ekl = ekl in Eqs. (2.11), (2.12),

(2.20), and (2,21).









For the updated Lagrangian formulation,


tSij = Cijkj tekl (2.23)

and,


ItCijkl tekl teij d


tI t
+ tij 6 tij tdV


= EVW rj 6 e tdV
(2.24)

For the total Lagrangian formulation,


oSij = oCijkj oekl (2.25)

and,


IoCijkl oekl 6 oeij odV
ddV


+ [sij ij dV


= EVW Sj e odV
o 0J o ij
(2.26)
It should be noted that the surface tractions and the
body forces in the calculation of external virtual work may
be treated configuration dependent when the structure
undergoes large displacements or large strains. If this is









the case, the external forces must be transformed to the

current configuration at each iteration [10, 11, 12].



2.7 Strain-Displacement Relationship
Using the von Karman Assumptions


The nonlinear strain terms can be simplified for the

plate or shell type structures using von Karman assumption

of large rotation.

In the mechanics of continuum the measure of

deformation is represented by the strain tensor Eij [52] and

is given by using index notation.


2Eij = ( ui,j + uj,i + k,iuk,j ) (2.27)

where,

ui = Displacement in i-direction.

uij = aui / axj

xi = Rectangular Cartesian coordinate axes, i=1,2,3.

uk,iuk,j = ul,iul,j + u2,iu2,j + u3,iu3,j

The von Karman theory of plate is a nonlinear theory

that allows for comparatively large rotations of line

elements originally normal to the middle surface of plate.

This plate theory assumes that the strains and rotations are

both small compared to unity, so that we can ignore the

changes in geometry in the definition of stress components

and in the limits of integration needed for work and energy









considerations [53]. It is also assumed that the order of

the strains is much less than the order of rotations.

If the linear strain eij and the linear rotation rij

are defined as

2eij = uij + uji (2.28)

2rij = u ui (2.29)

Then the sum of Eqs. (2.28) and (2.29) gives

2(eij + rij) = 2uij (2.30)

and the subtraction of Eq. (2.29) from Eq. (2.28) gives

2(eij rij) = 2uj (2.31)

From Eqs. (2.30) and (2.31), it is concluded that


uk,j = ekj + rkj (2.32)

uk,i = eik rik (2.33)

Eq. (2.33) can be rewritten as


uk,i = eki + rki (2.34)

since eik = eki from the symmetry of linear strain terms and

rik = -rki from the skew symmetry of the linear rotation
terms.

The strain-displacement Eq. (2.27) now becomes


2Eij = 2eij + (eki + rki)(ekj + rkj)


(2.35)









by substituting Eqs. (2.30) through (2.34) into Eq. (2.27).

Thus the nonlinear strain terms have been decomposed into

linear strain terms and linear rotation terms.

From the assumption on the order of strains and

rotations


eki << rki and ekj << rkj (2.36)

Thus Eq. (2.35) can be simplified as by ignoring eki and

ekj.

2Eij = 2eij + rkirkj (2.37)


The straight line remains normal to the middle surface

and unextended in the Kirchhoff assumption, but it is not

necessarily normal to the middle surface for the Mindlin

assumption. For both assumptions the generic displacements

u,v,w can be expressed by the displacements at middle

surface.

For the Kirchhoff plate [20],


u(x,y,z) = uo(x,y) z[Wo(x,y),x]

v(x,y,z) = vo(x,y) z[wo(x,y),y] (2.38)

w(x,y,z) = Wo(x,y)
where,

uo, Vo, wo = Displacements of the middle surface
in the direction of x, y, z.

u, v, w = Displacements of an arbitrary point
in the direction of x, y, z.









Now the linear strain components eij and the linear rotation
components rij can be calculated using Eqs. (2.28) and
(2.29).

ell = (ul,1 + Ul,1) = ul,1 = ,x
el2 = (l,2 + u2,1) = 1(uly + Vx)

el3 = (u1,3 + U3,1) = (-Wox + Wx)
22 = i(u2,2 + u2,2) = u2,2 = Vy (2.39)
e23 = (u2,3 + u3,2) = (-woy + Wy)

e33 = (u3,3 + u3,3) = u3,3 =

The rotation terms r12, r13, r23 are the rotation quantities
about the axes 3(z), 2(y) and l(x), respectively. For the
plate located in the xy plane, the rotation about z axis rl2
is much smaller than rotation about x axis r23 and y axis
r13 and therefore rl2 is assumed to be zero here. And it is
noted further that wo(x,y) is the same as w(x,y) and is a
function of only x and y so that w,3 = w,z = 0.


lrl21 << lr231 or 1r13 (2.40)

rl = (ul,1 ul,1) = 0

r12 = 1(ul,2 u2,1) = i(u, Vx) = 0

r13 = I(ul,3 u3,1) = (-Wox Wx) = -Wx
r22 = 1(u2,2 u2,2) = 0 (2.41)

23 = 1(u2,3 u3,2) = (-Wo'y W'y) = -Wy
r33 = i(u3,3 3,3) = 0









The linear strain component eij is symmetric and the linear

rotation component rij is antisymmetric.


eij = eji
rij = -rji (2.42)


The strain components from Eq. (2.37) can be rewritten using

Eqs. (2.39) and (2.41).


Ex = ell + (r112 +

yy = e22 + (r122 +
Ezz = e3 + 1(r132 +

Exy = el2 + 2(rllrl2
Exz = e3 + (rllrl3

Eyz = e23 + U(rl2r13


r21 + r31 ) = ell + r31

r22 + r322) = e22 + r32

r232 + r332) = 1(r132 + r232) = 0

+ r21r22 + r31r32) = el2 + r31r32
+ r21r23 + r31r33)

+ r22r23 + r32r33) (2.43)


Egz term is assumed to be zero because it does not have the

linear term. Exz and Eyz terms are transverse shear terms

which can be ignored for thin plate. Then Eq. (2.43) can be

rearranged as follows using Eqs. (2.41) if all the zero

terms are removed.


Exx = ell + r312 = ell + (W,x)2
Eyy = e22 + r322 = 22 + )2

Ey = el2 + Ir31r32 = e2 + (W'x) (,y)

Exz = e13

Eyz = e23


(2.44)







33


Thus the decomposition of exact strain components has been

done using the Kirchhoff plate assumptions (2.38) and the

von Karman assumption (2.40) on the magnitude of rotation.

It is noted that all the inplane displacement gradients in

nonlinear strain terms are ignored through von Karman

assumptions [20]. This fact will be applied in chapter 5.













CHAPTER 3
THREE DIMENSIONAL LINK ELEMENT



3.1 Element Description



The link element used here is based on the two

dimensional element developed by Cleary [54]. The link

element is based on the following assumptions.

Any normal compressive force is transferred to the

other side of the link without any loss. To facilitate this,

a very limited amount of loss through displacement should be

allowed. Currently, this limited displacement is defaulted

to .001 units, while it is a input parameter. The link

separates in response to any net tension, losing its normal

stiffness.

To discuss the shear force transfer, some definitions

for friction are needed. The force to start one body sliding

along the other body is called the static friction force.

The force to keep it moving is the kinetic friction force.

There are two corresponding coefficients of friction, static

friction coefficient and dynamic friction coefficient, where

the static friction coefficient will generally be greater

than the dynamic friction coefficient.









Two laws of friction were used in the link element. The

first law is that the frictional force is proportional to

the normal force, with the constant of proportionality being

the friction coefficient. The second law is that friction

does not depend on the apparent area of the connecting

solids, i.e., it is independent of the size of the bodies.

The shear force is transferred through friction. The

uncertainty in friction is the factor which limits the

overall accuracy of the calculation. Therefore, it is

assumed that the static friction coefficient is proportional

to the dynamic friction coefficient. For nonmetallic

materials, the ratio of dynamic coefficient to static

coefficient is about 0.75.

The link element is composed of two surfaces. If the

shear force is less than or equal to the static friction

force, i.e., coefficient of friction times the normal force,

the shear force is balanced by the friction force and the

total force is transferred. This is shown in Fig. 3-1. But

if the shear force is greater than the static friction

force, one surface of the link element will move along the

other surface. In this case there will be a dynamic friction

force which is less than the shear force. This dynamic

friction force can only resist a portion of the shear and

the system is not in static equilibrium. Therefore, if the

shear force is greater than the static friction force, the

link element will lose its shear stiffness. This can also be










modeled with a body on roller and spring as shown in Fig. 3-

2. The spring model of the link element is shown in Fig. 3-

3.

The link element here has four nodes and each node has

three translational degrees of freedom in local u-, n-, and

w-directions. The total number of element degrees of freedom

is 12. The element degrees of freedom are shown in Fig. 3-4.

The equivalent "strain" for the link element is defined

as the average deformation at the center of the element. The

average deformation corresponding to the translational

degrees of freedom, i.e., uo, vo, and wo, can be directly

calculated from the joint displacements by averaging the

difference in nodal displacements at the ends of element in

local u-, n-, and w-direction in turn. The relative rotation

at the center of the element, ro, can be found using nodal

displacements in local n-direction and the element length as

shown in Fig. 3-5. This angle is not an "average" value but

the "relative" change in angle of the center line due to

rotation.

The two joint parameters must be introduced. These are

kn, the unit stiffness normal to the joint, and ks, the unit

stiffness along the joint.

The off-diagonal term kns to account for dilatation

during shearing is ignored because this joint element will

model the dry joint between concrete box girder bridge

segments. No significant dilation is expected in this case.










Some values of kn and ks were reported in geotechnique

area [3]. As the values are those for natural joints, they

do not directly apply to this case.

From the test results [7], it can be seen that the

shear stiffness of dry joint ranges from 70,000 to 286,000

psi per inch at ultimate. In case of frictional strength,

this can be interpreted as linear behavior between the

origin and the ultimate point.

It seems reasonable that the normal stiffness of the

element, kn, is assumed to be stiffer than the connected

material by the order of 103 to transfer the normal force

without any significant loss. The forces are either totally

transferred in compression or totally lost in tension. The

latter case has no problem related to the value of kn.

The shear stiffness parameter is more difficult to

define. The data available is so limited that even a

statistical treatment cannot be done. But in the analysis of

structural behavior up to the ultimate, these properties do

not have great influence because the forces are transferred

through friction.

The shear stiffness becomes zero upon sliding. But

there is some 'residual' shear force. This 'residual' force

is equal to friction force. Therefore, if shear displacement

is more than the displacement just before the sliding the

shear stiffness is set to be zero.

























Ff


P, N = External forces.

F = Friction force.
f


m = Friction coefficient






1) P < or = mN then P = Friction Force. In Equilibrium.

2) P > mN then the body moves but the frictional force

mN is acting against the other body.


Fig. 3-1 Friction Force







FRICTIONAL SPRING WITH SHEAR STIFFNESS
Force in spring = mN
N


-%AAN-


BEFORE


SLIDING


SFRICTIONAL SPRING WITH ZERO STIFFNESS


N


A F

mN


mN < F

AFTER SLIDING


Fig. 3-2 Spring Model of Friction Force


__ __






















Kn = Zero under tension.
= Very large under
compression.


SPRING MODEL FOR NORMAL FORCES


SPRING MODEL FOR SHEAR FORCES


Fig. 3-3 Spring Model of Link Element











11


K


Fig. 3-4 Element Degrees of Freedom of Link Element


L

12

I

34






CENTER OF ELEMENT
K


UK
NoLNo


uo = [(UK + UL) (UI + UJ)] / 2

ro= p a

S[(VK-VL) (VJ-VI)] / L


VK


O------



VI
o o-


Fig. 3-5 Element "Strain"


uo-_ZL









3.2 Formation of Element Stiffness


There are four nodes per element. Each node has three

degrees of freedom corresponding the translational

displacements in u-, n-, and w-direction resulting in 12

element degrees of freedom as shown in Fig. 3-4. The element

stiffness is derived directly from the physical behavior of

the element described in section 3.1.

The mathematical symbol {} is used for a column vector

and [] for a matrix. The nodal displacement column vector

{q)(12) is composed of 12 translational nodal displacements

corresponding to the 12 element degrees of freedom.


(q} = { ui vi wi Uj Vj Wj uk vk wk ul vl wl )T

The "strain" is defined as the average deformation at

the center of the element as shown in Fig. 3-5. The "strain"

column vector {e)(4) is


{e} = { uo vo wo ro )T


where,

uo = ( uk + ul ) / 2 ( u + j ) / 2

o = vk + 1 ) / 2 ( + Vj ) / 2

w = (wk + 1 ) / 2 ( i + j ) / 2

ro = ( vk v ) / L ( vj vi ) / L

where,

L = The length of the element.











uo, vo, wo = Average nodal displacements in local
u-, n-, w-directions.

ro = The relative angle change about local z axis.


Therefore the relationship between "strain" and nodal

displacements is


(e}(4) = [B](4,12) (q)(12)


The [B](4,12) matrix which gives strains due to unit

values of nodal displacements is



-0.5 0.0 0.0 -0.5 0.0 0.0 0.5 0.0 0.0 0.5 0.0 0.0
0.0 -0.5 0.0 0.0 -0.5 0.0 0.0 0.5 0.0 0.0 0.5 0.0
0.0 0.0 -0.5 0.0 0.0 -0.5 0.0 0.0 0.5 0.0 0.0 0.5
0.0 1/L 0.0 0.0 -1/L 0.0 0.0 1/L 0.0 0.0 -1/L 0.0



The "stress" is defined as the normal and shear stress

per unit of area. {s} is the average stress on the surface

due to the two nodal forces exerted in the plane of the

surface. This stress is in equilibrium with the stress on

the other surface of the element as shown in Fig. 3-6. m is

the moment of the nodal forces on one surface in local n-

direction about the center of the element. This moment is

also balanced by the moment of the nodal forces on the other

surface of the element. This moment is used to define the

distribution of the normal stress of the element as shown in

Fig. 3-7.










Pnl

fl1


P...'


uk


-~- -~-
-~- -~- -~-
-~- -~- -~-
-~- -~-


uj -nj



Pj J ~pnj


-W- Stresses

-- Nodal Forces


Local Coordinate System


Fig. 3-6 Nodal Forces and Stresses of Link Element


/


ni


-W.- -.- -w

-0--o 0--0 p











FORCE TRANSFER THROUGH ONE EDGE
OF THE LINK ELEMENT






P



L K
Vo ro

2 Vo


CENTER OF ELEMENT


p J


Fig. 3-7 Element "Strain", m









The "stress" column vector (s}(4) is


(s) = { Sx, sn, Sz, m )

The "stress-strain" relationship is


(s)(4) = [E](4,4) {e)(4)

where,



kx 0 0 0
[E] = 0 kn 0 0
0 0 kz 0
0 0 0 km



where km can be related to kn using the definition of

the moment m, i.e.,

m = (Sn)(L)(t)(0.5)(L)

= kmro = (km)(vo/(0.5)(L))



Thus, km = (knVo)(L)(t)(0.5)(L) / [Vo/(0.5)(L)]

= (0.25)(t)(kn)(L3)


where, t = Element thickness.

This assumes that there is no coupling between the

shear stress and normal stress.

The element nodal force column matrix (P)(12) is

composed of the 12 nodal forces shown in Fig. 3-6.










(P) Pui Pni Pwi Puj Pnj Pwj Puk Pnk Pwk

Pul nl Pwl T

Stress can then be related to nodal forces using the

definition of stress and force equilibrium between the two

surfaces of the element.

By the definition of stress,


sn = (1/Lt)( Pnk + Pnl ) (3.1)

sx = (1/Lt)( Puk + Pul ) (3.2)

sz = (1/Lt)( Pwk + Pwl ) (3.3)
m = Pnk(0.5)(L) Pnl(0.5)(L) (3.4)
where, Lt = (L)(t)

By force equilibrium of the two surfaces,

Pi = -P and Pj = Pk (3.5)

To express the element nodal forces in terms of the
stress, we use Eqs. (3.1) through (3.5) to find the force

recovery matrix [FR]. [FR] gives the nodal forces in

equilibrium with the element stresses.

From (Eq. (3.1) + Eq. (3.4)),

2Pnk = (L)(t) (sn) + 2(m)/L

Pnk = 0.5(L) (t) (sn) + (l/L)(m)









From Eq. (3.5),

Pj = -Pk
nj = -Pnk
= -0.5(L)(t)(sn) (1/L)(m)

From Eq. (3.1),


Pn1 = (L)(t)(sn) Pnk
= (L)(t)(sn) ( 0.5(L)(t)(sn) + (1/L)(m))
= 0.5(L) (t) (n) (1/L) (m)

From Eq. (3.5),

Pni = Pn1
= -0.5(L)(t)(sn) + (1/L)(m)

From the assumption that Puk = Pul and Eq. (3.2),


Puk = (L)(t)(sx)/2

Pul = (L)(t)(sx)/2

From Eq. (3.5),


Pui = ul = -(L)(t)(sx)/2

Puj = uk = -(L)(t)(sx)/2

From the assumption that Pwk = Pwl and Eq. (3.3),


Pwk = (L)(t)(sz)/2
Pwl = (L)(t)(sz)/2










From eqn 5,


Pwi = wl = -(L)(t)(sz)/2

Pwj = Pwk = -(L)(t) (s)/2


Therefore, the force-stress relationship is


{P)(12) = [FR](12,4) (s)(4)


where the force recovery matrix [FR](12,4) is


[FR] =


-Lt/2
0
0
-Lt/2
0
0
Lt/2
0
0
Lt/2
0
0


0
-Lt/2
0
0
-Lt/2
0
0
Lt/2
0
0
Lt/2
0


0
0
-Lt/2
0
0
-Lt/2
0
0
Lt/2
0
0
Lt/2


0
1/L
0
0
-1/L
0
0
1/L
0
0
-1/L
0


And this relationship is further expanded using the

stress-strain relationship and the strain-nodal displacement

relationship as follows.


(P}(12) = [FR](12,4)

= [FR](12,4)

= [Bt](12,4)


[E](4,4) {e)(4)

[E](4,4) [B](4,12)

[E](4,4) [B](4,12)


Then finally this can be symbolized as equilibrium

equation.


(q)(12)

{q)(12)










{P)(12)= [Ke](12,12) (q}(12)


where [Ke] = [Bt][E][B]

Here it is noted that [FR] = [Bt] and [Ke] = [Bt][E][B] just

as in the case of common finite element method.

The final element stiffness matrix [Ke] is


(L)t

4


kx
0
0
kx
0
0
-kx
0
0
-kx
0
0


0
2kn
0
0
0
0
0
0
0
0
-2kn
0


0
0
kz
0
0
kz
0
0
-kz
0
0
-kz


kx
0
0
kx
0
0
-kx
0
o
-kx
0
0


0
0
0
0
2kn
0
0
-2kn
0
0
0
0


0
0
kz
0
0
kz
0
0
-kz
0
0
-kz


-kx
0
0
-kx
0
0
kx
0
0
kx
0
0


0
0
0
0
-2kn
0
0
2kn
0
0
0
0


0
0
-kz
0
0
-kz
0
0
kz
0
0
kz


-kx
0
0
-kx
0
0
kx
0
0
kx
0
0


0
-2kn
0
0
0
0
0
0
0
0
2kn
0


0
0
-kz
0
0
-kz
0
0
kz
0
0
kz


This matrix can be

standard rotation.


rotated to any direction using the


3.3 Solution Strategy


The structural stiffness changes because of the slip

and debonding of the link. Therefore, the process of the

resistance of the total structure physically becomes

nonlinear. Correspondingly, special solution techniques for

nonlinear behavior are needed.

This can be done using the iterative solution technique

with initial stiffness or tangent stiffness. The latter can










be formed by assembling the structural stiffness at the

beginning of each iteration and this converges faster than

the initial stiffness.

A third solution strategy for this case is event-to-

event technique which is usually employed for the linear

stiffnesses between any two "events," which are defined as

the intersection point between two linear segments. This

also provides means of controlling the equilibrium error.

Any significant event occurring within any element

determines a substep. The tangent stiffness is modified in

each substep, and hence, the solution closely follows the

exact response.

3.4 Element Verification



3.4.1 SIMPAL

The finite element analysis program SIMPAL [55], is

used to implement and verify the element formulation. SIMPAL

was chosen for the initial implementation because that was

the original implementation done by Cleary [54]. This way,

the 3-D aspects could be implemented and verified using

Cleary's original program. A table of the element

verification is shown in Fig. 3-8 and Fig. 3-9.









LOADING


10




10
Y
1 2


NODE THEORY SIMPAL ERROR

DISP 2 -.1333 -.1333 .000

DISP 4 -.1333 -.1333 .000

STRESS N/A -80 -80 .000


NODE THEORY SIMPAL ERROR

DISP 2 -.1333 -.1333 .000


DISP 4 -.1333 -.1333 .000

STRESS N/A -80 -80 .000



NODE THEORY SIMPAL ERROR


DISP 2 -.1333 -.1337 .003

DISP 4 -.1333 -.1337 .003

STRESS N/A -80 -80 .000


* NODE 2 Y DISP =-0.1017-04 Z DISP = -0.8684-05
NODE 4 Z DISP =-0.8684-05 Y DISP =-0.1017-04
2 2
SQRT((.1017) + (.08684) )=0.1337


Fig. 3-8 Link Element Test Using SIMPAL


3
3


THICKNESS = .25
Ks = 3E6
Kn = 6E6


RESULTS








LOADING


4 6
Aft Am


.P -......


3 5
1 _~~


RESULTS


8
Aft


Ii
MP, -


Adh -


10


1 I


1 3 5 7 Y
1 3 5 7 Y


Z
2


4 6


3 5


Fig. 3-9 Combined Test Model for SIMPAL


Z
2


-~araooa a;a8oo~ uru -


ia--m i-- i-o -n


----------


HtS~ ^~ k


10
8


10










3.4.2 ANSR



The test examples used are the same as those used in

the initial element verification using SIMPAL. The results

from ANSR [56] are exactly the same as those from SIMPAL.

The link element was tested further using a modeled membrane

element composed of 22 truss elements as a membrane element

was not available at the time of element verification in

ANSR. The results are shown in Table 3-1 and the structures

used are shown in Fig. 3-10 and Fig. 3-11.












Table 3-1 Displacements of
for ANSR


Truss Model


Node Truss Truss Diff.
No. only w/ LINK (%)

10-x -.1027e-4 -.1049e-4 2.2

10-y -.1990e-5 -.2010e-5 0.9

ll-x -.9017e-5 -.9211e-5 2.2

11-y -.4906e-6 -.4973e-6 1.4

12-x -.9915e-5 -.1049e-4 2.2

12-y +.9742e-6 +.9703e-6 0.4











LINK ELEMENT


20




.2

20







.2


20


Fig. 3-10 Combined Test Model for ANSR

















































Fig. 3-11 Truss Model for ANSR














CHAPTER 4
LINEAR SHELL ELEMENT


4.1 Element Description



The shell element is formulated through a combination

of two different elements, the membrane element and the

Mindlin plate bending element [57].

The Mindlin plate element is different from the

Kirchhoff plate element in that the former allows transverse

shear deformation while the latter does not.

The common portions of the formulation of two elements

are


1. Formation of the shape functions.

2. Formation of the inverse of Jacobian matrix.


These processes can be done at the same time. The four-

to nine-node shape functions and their derivatives in rs-

space can be formed and then transformed into xy-space

through the inverse of Jacobian matrix.



4.2 Formulation of Shape Functions



The formulation of shape functions starts with three

basic sets of shape functions shown in Fig. 4-1.










1. The bilinear shape functions for four-node element.

2. The linear-quadratic shape functions for nodes
five to eight of the eight-node element.

3. The bubble shape function for node nine of
nine-node element.


These shape functions can be formulated directly from

the local coordinates of the element nodes through the

multiplication of the equations of the lines which have zero

values in the assumed displacement shapes and the scale

factor to force the shape function value to one at the node

for which the shape function is formed. The derivative of

each shape function with respect to r and s is then

evaluated from the shape function expressed in terms of r

and s.

If node nine exists, the value at node nine of shape

functions one to eight must be set to zero. The value of the

bilinear shape functions for a four-node element at node

nine is one fourth and the value of the linear-quadratic

shape functions for the five- to eight-node element at the

node nine is one half. This can be forced to zero using the

bubble shape function of the nine-node element because this

shape function has the value of one at node nine and zero at

all other nodes. Therefore the modification is the

subtraction of one fourth of the value the bubble shape

function has at node nine from the each shape function for

the corner nodes and the subtraction of one half of the










value of the bubble shape function of node nine for the

nodes five to eight, whichever exists.

If any of the center nodes on the edge of the element

(any one of nodes five to eight) exists, the bilinear shape

functions of four-node element must be modified further

because the value at center of the edge is one half in those

bilinear shape functions. This can be done by subtracting

one half of the linear-quadratic shape function for the

newly defined center node on the edge of the element from

the bilinear shape functions of the two adjacent corner

nodes. The value of any five node shape functions at the

corner node is zero. Therefore, no further consideration is

needed except for the shifting of the shape functions in the

computer implementation. These processes are shown in Fig.

4-2.

If any of the linear-quadratic shape functions of nodes

five to eight is missing, all the linear-quadratic shape

functions thereafter and the bubble shape function must be

shifted to the proper shape function number. For example, if

linear-quadratic shape function five is missing, then the

shape functions six to eight must be shifted to five through

seven and the bubble shape function must be shifted to the

node eight because all of the linear-quadratic shape

functions have been defined and numbered as shape functions

for the nodes five through eight and the bubble shape

function for the node nine.


















Four Node Element



Shape Function for Corner Node


Five Node Element


Shape Function for Edge Center Node








Nine Node Element




Shape Function for Element Center Node


Fig. 4-1 Three Basic Shape Functions























SF 1


SF 2


SF 3


SF 3


SF 4 = (SF 1) (1/4) (SF 3)


SF 5 = (SF 2) (1/2) (SF 3)


Fig. 4-2 Formation of Shape Functions









4.3 The Inverse of Jacobian Matrix



While the generic displacements are expressed in terms

of rs-coordinate, the partial differential with respect to

the xy-coordinate is needed for the calculation of strain

components. Thus the inverse of the Jacobian matrix must be

calculated. This can directly be found from the chain rule

using the notation (a,b) defined as the partial differential

of function a with respect to the variable b for simplicity.


f,x = (f,r)(r,x) + (f,s)(s,x)

f,y = (f,r)(r,y) + (f,s)(s,y)


In matrix form,



f,x r,x s,x f,r Jll-1 J12-1 f,r

f,y r,y s,y f,s J21-1 J22-1 f,s



The inverse of Jacobian matrix


But the terms in the inverse of the Jacobian matrix are not

readily available because the rs-coordinate cannot be solved

explicitly in terms of xy-coordinate. On the other hand, for

the isoparametric formulation, the geometry is interpolated

using the nodal coordinate values(constants) and the

displacement shape functions in terms of r and s. Thus the

generic coordinate x and y can be expressed in r and s









easily and explicit partial differentials of x and y with

respect to r and s can be performed. Therefore the Jacobian

matrix is computed and then inverted.

The Jacobian matrix is derived by the chain rule.


f,r = (f,x)( x,r) + (f,y)( y,r)

f,s = (f,x)( x,s) + (f,y)( y,s)


In matrix form,



f,r x,r y,r f,x Jll J12 f,x

f,s x,s y,s f,y J21 J22 f,



Jacobian matrix


nn
Let s be .
i=1

where nn = number of nodes (4 to 9).

From geometric interpolation equations,


x = Z fi*xi

y = Z fi*Yi


The terms in the Jacobian matrix are


J11 = x,r = (z fi*xi),r = Z ((fi,r) xi)

J12 = y,r = (Z fi*yi),r = Z ((fi,r) yi)

J21 = x,s = (Z fi*xi),s = Z ((fi,s) xi)

J22 = y,s = (Z fi*yi),s = Z ((fi,s) yi)











xi, yi are coordinate values of the element and are

constants and therefore can be taken out of the partial

differentiation.

The inverse of two-by-two Jacobian matrix can be found

-i
Jl-1 = r,x = J22 / det(J)
-i
J12 = s,x = -J12 / det(J)
-I
J21- = r,y = -J21 / det(J)
-i
J22- = s,y = J11 / det(J)

where det(J) = J11J22 12J21



4.4 Membrane Element



The formulation of the membrane element used for the

implementation follows the procedure shown on pages 115

through 118 in reference [57]. The ( ) symbol will be used

for the column vectors.

Nodal displacements are the nodal values of two in-

plane translations and denoted as {ui vi}T. The generic

displacements are defined as two translational displacements

at a point and denoted as { u v )T. By the word generic it

is meant that the displacement is measured at an arbitrary

point within an element. The generic displacements u and v

can be calculated using shape functions. The shape function

is a continuous, smooth function defined over the closed










element domain and is differentiable over the open domain of

the element. The shape function is also the contribution of

displacement of a node for which the shape function has been

defined to the generic displacement. Thus the generic

displacement at an arbitrary point can be found by summing

up all the contributions of all the nodes of the element.

The displacement interpolation equations are


u = Z fi Ui

v = Z fi Vi


In the isoparametric formulation the geometry is

interpolated using the same shape function assumed for the

displacement interpolation.

Therefore, the geometry interpolation is


x = Z fi xi

y = Z fi Yi


where,

fi = Shape function for node i.

xi, Yi = Coordinates of node i.

ui, vi = Displacements at the node i.

u, v = Displacements at an arbitrary point within

an element.









The three in-plane strain components for a membrane

element are


{ E ) = ( ex Cy 7xy )

These strain components can be found through the

partial differentials of the generic displacements with

respect to xy-coordinates.

ex = ux

Ey = v,y

7xy = uy + VX

Using the inverse of the Jacobian matrix, the strain

components can be evaluated.


eX

= u,x

= (u,r)(r,x) + (u,s)(s,x)

= (u,r)(J11-1) + (u,s)(J12-1)

= ((Efiui),r)(J11-1) + ((EfiUi),s)(J12-1)

= Z[(fi,r)(r,x) + (fi,s)(s,x)] ui



y
= v,y
= (v,r)(r,y) + (v,s)(s,y)

= (v,r)(J21-1) + (v,s)(J22-1

= ((.fivi),r)(J21-1) + ((Zfivi),s)(J22-1)

= [(fi,r)(r,y) + (fi,s)(s,y)] vi











7xy
= u,y + v,x

= [(u,r)(r,y) + (u,s)(s,y)] + [(v,r)(r,x) + (v,s)(s,x)]

= [((Zfiui),r)(J21-1) + ((Zfiui),s)(J22-1)]

+ [((Zfivi),r)(J11-) + ((Zfivi),s)(J12-1)]
= Z[(fi,r)(r,y) + (fi,s)(s,y)] ui

+ Z[(fi,r)(r,x) + (fi,s)(s,x)] vi


New notations are introduced here to simplify the

equations. These are ai and bi and defined as follows:


ai = (r,x)(fi,r) + (s,x)(fi,s) = fi,x

bi = (r,y)(fi,r) + (s,y)(fi,s) = fiY

Then the strain terms above become


ex = Zaiui = Zf,x ui

Ey = Zbivi = Zfiy vi

7xy = Zbiui + Zaivi = Zfi,y ui + Zfi,x vi
In matrix form,



ex ai 0 ui
yE = 0 bi vi
7xy bi ai


In symbolic form,


[E] = E[Bi]Cqi]









where,



ai 0 fix 0
[Bi] = 0 bi 0 fi',
bi ai fi'y fix


and,


ui
[q] = Vi
vi


Therefore the strain at an arbitrary point within an

element is


[e] = [Bl][q1] + [B2][q2] + ... + [B9][q99


= [ B1 B2 B3 B4 B5 B6 B7 B8 B9 ]


ql
q2
q3
q4
95
q6
q7
q8
q9


The size of the vectors and matrix are


[e(3,1)] = [B(3,18)][q(18,1)]


In the actual calculation, this can be done by summing

up the [Bij[qi] over all the nodes for the given coordinates

of the point under consideration, i.e., the coordinates of

one of the integration points.










The stresses corresponding to the strains are

{ a } = { ax ay rxy )T

The stress-strain relationship of an isotropic material

is



Ell E12 0
[E] = E21 E22 0
0 0 E33


where,

E11 = E22 = E / ( 1 2 )

E12 = E21 = pE / ( 1 p2)

E33 = G

where,

E = Young's modulus

p = Poisson's ratio

G = shear mQdulus = E / ( 2*(1+p))

The integration over volume must be introduced for the

calculation of the element stiffness and equivalent nodal

loads of distributed loads or temperature effects. As the

thickness of the element is constant, the integration over

volume can be changed to the integration over area.

The element stiffness related to the degrees of freedom

of the node i can be calculated through the volume


integration of Bi(2,3)E(3,3)Bi(3,2).










[Ki] = BT E Bi dV


As [B] and [E] are constant about z, the integration

through the thickness from -t/2 to t/2 can be performed on z

only and yields


T
[Ki] = [ Bi(t)EBi] dA
A B

T -
= [ Bi E Bi] dA
A


where,


E = tE


The size of membrane element stiffness is 18 by 18.


[K] =


[ E ] [ B1 B2 B3 ... B8 B9 ] dV
(3,3)
(3,18)


B8
B9

(18,3)


Equivalent nodal loads due to body forces on the

membrane element are calculated as


Pb Jv fTbdV =11 1-1


fTbIJI dr ds


in which {b} = ( 0 0 bz T or { 0 b 0 }T or { bx 0 0 } in










accordance with the direction of gravity in the coordinate

system used. The nonzero quantities bx, by, or bz represent

the body force per unit area in the direction of

application.

Equivalent loads caused by initial(temperature) strains

are



PO = BTEEO dV
JV
1 1
BTEcIJI dr ds
S-1 J-i

where,

(C0)= exxO yyO 0 0 0 )T
= {aAT aAT 0 0 0 )



4.5 Plate Bending Element



The formulation of the plate bending element used for

the implementation has followed the procedures shown on

pages 217 through 221 in reference [57]. The { } symbol will

be used for the column vectors.

Many plate bending elements have been proposed. The

most commonly used are Kirchhoff plate elements and Mindlin

plate elements.










Kirchhoff theory is applicable to thin plates, in which

transverse shear deformation is neglected. The assumptions

made on the displacement field are


1. All the points on the midplane(z = 0) deform only
in the thickness direction as the plate deforms in
bending. Thus there is no stretching of midplane.

2. A material line that is straight and normal to the
midplane before loading is to remain straight and
normal to the midplane after loading. Thus there is
no transverse shear deformation (change in angle
from the normal angle).

3. All the points not on the midplane have displacement
components u and v only in the x and y direction,
respectively. Thus there is no thickness change
through the deformation.


Strain energy in the Kirchhoff plate is determined entirely

by in-plane strains ex, Cy, and 7xy which can be determined

by the displacement field w(x,y) in the thickness direction.

The interelement continuity of boundary-normal slopes is not

preserved through any form of constraint.

Mindlin theory considers bending deformation and

transverse shear deformation. Therefore, this theory can be

used to analyze thick plates as well as thin plates. When

this theory is used for thin plates, however, they may be

less accurate than Kirchhoff theory because of transverse

shear deformation. The assumptions made on the displacement

field are


1. A material line that is straight and normal to the
midplane before loading is to remain straight but
not necessarily normal to the midplane after









loading. Thus transverse shear deformation (change
in angle from normal angle) is allowed.

2. The motion of a point on the midplane is not
governed by the slopes (w,x) and (w,y) as in
Kirchhoff theory. Rather its motion depends on
rotations Ox and 0 of the lines that were normal to
the midplane of the undeformed plate. Thus 0, and 0
are independent of the lateral displacement w, i.e.,
they are not equal to (w,x) or (w,y).


It is noted that if the thin plate limit is approached, -xz

= 7yz = 0 because there is no transverse shear deformation.
In this case the angles 0x and 0y can be equated to the

(w,x) and (w,y) numerically but the second assumption still

holds.

The stiffness matrix of a Mindlin plate element is

composed of a bending stiffness [kb] and a transverse shear

stiffness [ks]. [kb] is associated with in-plane strains ex,

cy, and 7xy. [ks] is associated with transverse shear

strains 7xz and 7yz. As these two groups of strains are

uncoupled, i.e., one group of the strains do not produce the

other group of strains, the element stiffness can be shown

as [82]



[k] = ( BbEBb ) dA + (BsEBs) dA


because BbEBs = BsEBb = 0 from uncoupling (corresponding E =

0). Each integration point used for the calculation of [kg]

places two constraints to a Mindlin plate element,

associated with two transverse shear strains lyz and yzx. If










too many integration points are used, there will be too many

constraints in transverse shear terms, resulting in locking.

Therefore, a reduced or selective integration can prevent

shear locking. Or, the transverse shear deformation can be

redefined to avoid such locking.

For example, a bilinear Mindlin plate element responds

properly to pure bending with either reduced or selective

integration. But full two-by-two integration is used for

pure bending, shear strains appear at the Gauss points as

shown in Fig 4-3. As the element becomes thin, its stiffness

is due almost entirely to transverse shear. Thus, if fully

integrated, a bilinear Mindlin plate element exhibits almost

no bending deformations, i.e., the mesh "locks" against

bending deformations.

Nodal displacements for the plate bending consist of

one out-of-plane translation and two out-of-plane rotations

and are denoted as { wi 0xi 8yi )T. The rotations are chosen

independently of the transverse displacement and are not

related to it by differentiation. Thus the transverse shear

strains jxz and 7yz are considered in the formulation

resulting in five strain components. The generic

displacements are defined as three translational

displacements and denoted as { u v w }T. By the word generic

it is meant that the displacement is measured at an

arbitrary point within an element. These generic

displacements are different quantities from the nodal










displacements and therefore must be related to the nodal

displacements.

The generic displacements u and v can be calculated as

functions of the generic out-of-plane rotations using the

small strain(rotation) assumption. The relationship between

generic displacements and rotation is shown in Fig 4-4.


u = zBy

v = -ZBx


The generic displacements Ox and By can be found using the

assumed displacement shape functions and the corresponding

nodal displacements Oxi and 0yi.

The generic displacement w does not need any conversion

because it corresponds to the nodal displacement wi.

In the isoparametric formulation the geometry is

interpolated using the same shape function assumed for the

displacement interpolation.

The displacement interpolation is


8x = Z fi 0xi

6y = Z fi 0yi

w = Z fi Wi


Similarly, the geometric interpolation is


x = Z fi xi

y = Z fi Yi


















Zero Shear
Strain


One Point Gauss Integration


I Non-zero
Shear Strain


Two Point Gauss Integration


Fig. 4-3 Shear Strains at Gauss Point(s)


i~Amlh









Z

+ u



y

X



Positive small rotational angle about y-axis gives

positive generic displacement in x-direction ( u ).

Shown is xz-plane.


Positive small rotational angle about x-axis gives
negative generic displacement in y-direction ( v
Shown is yz-plane.


Fig. 4-4 Displacements due to Rotations










where,

fi = shape function for node i

xi, Yi = coordinates of node i



Therefore,


u = Zoy = z z fi 6yi

v = -zOx = -z Z fi Oxi

W = E fi Wi


The five strain components for plate bending element

are { ex ey 7xy 7xz Iyz )T. These strain components can be

found through the partial differentials of the generic

displacements with respect to xy-coordinates.


e = uI,

Ey = v,y

7xy = uy + V,X

7XZ = u,z + W,X

7yz = vz + w,y

Using the inverse of the Jacobian matrix found, the

strain components can be evaluated.


Ex

= u,X = (ZBy),X

= (u,r)(r,x) + (u,s)(s,x)

= (u,r)(J11-1) + (u,s)(J12-1)









= ((ZzfiOyi),r)(J~l-1) + ((zzfiyi) ,) (J12-1)
= z Z[(fi,r)(r,x) + (fi,s)(s,x)] Oyi




= v,y = (-zOx),y
= (v,r)(r,y) + (v,s)(s,y)
= (v,r)(J21-1) + (v,s)(J22-1)
= ((-zzfioxi),r)(J21-1) + ((-zZfi0xi),s)(J22-1)
= -z Z[(fi,r)(r,y) + (fi,s)(s,y)] Oxi


7xy
= u,y + v,x = (zoy),y + (-zox),X
= [(u,r)(r,y) + (u,s)(s,y)] + [(v,r)(r,x) + (v,s)(s,x)]
= [((zzfioyi),r)(J21-1) + ((zZfiyi),s)(J22-1)]
+ [((-zZfisxi),r)(Jll-1) + ((-zzfixi),s) (J12-1)]
= z Z[(fi,r)(r,y) + (fi,s)(s,y)] Oyi
-z Z[(fi,r)(r,x) + (fi,s)(s,x)] exi


7xz = (u,z) + (w,x)
= (ZOy),Z + w,x = 9y + w,x

= Zfisyi + [(w,r)(r,x) + (w,s)(s,x)]
= Zfioyi + [((Zfiwi),r)(J11-1) + ((Zfiwi),s) (J12-1)
= Zfieyi + [((Zfi,r)wi) (J11-1) + ((Zfi,s)wi)(J12-1)]
= ZfiOyi + z[(fi,r)(r,x) + (fi,s)(s,x)] wi










lyz = (v,z) + (w,y)
= (-z0x),z + w,y = (-OX) + w,y

= (-Zfioxi) + [(w,r)(r,y) + (w,s)(s,y)]

= (-Zfioxi) + [((Zfiwi),r)(J21-1) +((Zfiwi),s)(J22-1)]

= (-Zfioxi) + [((Zfi,r)wi)(J21-1) +((Zfi,s)wi)(J22-1)]

= (-Zfi0xi) + Z[(fi,r)(r,y) +(fi,s)(s,y)] Wi


New notations are introduced here to simplify the

equations. These are ai and bi and defined as follows:


ai = (r,x)(fi,r) + (s,x)(fi,s) =

bi = (r,y)(fi,r) + (s,y)(fi,s) =

Then the strain terms above become


f.l
i1,x


ex = z aioyi

Cy = -z Zbixi

7xy = z Zbiyi z Zaixi

7xz = fi0yi + Zaiwi

Tyz = fioxi + Zbiwi


In matrix form,



eX
fy
7xy =
7XZ
Lyz

In symbolic form,


[e] = E[Bi][qi]


0
-zbi
-zai
0
-fi
i


zai
0
zbi
fi
0


wi
exi
Oyi










where,


[Bi] =




or,




[Bi] =


0
- zbi
- zai
0
fi
i


0
0
0
fi,x
fi'y


zai
0
zbi
fi
0


0
- zfi,y
- zfi,x
0
fi


zfi,x
0
zfi,Y
fi
0


and,


[qi]


wi
ixi
Oyi


Therefore the strain at an arbitrary point within an

element is


[e] = [B][ql] ] + ... + [B2[ + + 9][q9]


= [ B1 B2 B3 B4 B5 B6 B7 Bg B9 ]


q91
92
93
94
95
96
97
98
99










The size of the vectors and matrix are


[E(5,1)] = [B(5,27)][q(27,1)]


In the actual calculation, this can be done by summing

up the [Bi][qi] over all the nodes for the given coordinates

of the point under consideration, i.e., the coordinates of

one integration point.

The stresses corresponding to the strains are


{ o ) = ( ax ay rxy rxz yz )TT

The stress-strain relationship of an isotropic material

is



E11 E12 0 0 0
E21 E22 0 0 0
[E] = 0 E33 0 0
0 0 0 E44 0
0 0 0 0 E55


where,

E11 = E22 = E / ( 1 p2 )

E12 = E21 = pE / ( 1 p2)

E33 = G

E44 = E55 = G / 1.2

where,

p = Poisson's ratio

G = shear modulus = E / ( 2*(1+p))









The form factor 1.2 for the E44 and E55 terms is
provided to account for the parabolic distribution of the
transverse shear stress rzx over a rectangular section. This
form factor 1.2 can be shown from the difference in
deflections of a cantilever beam at its free end [58].
Let a beam have a rectangular cross section of
dimensions b by t with a length of L. If P is the transverse
shear force, then the parabolic distribution of the
transverse shear stress rz, is

r'X = (3P/2bt3)(t2 4z2)
where z = 0 at the neutral axis.
Then the transverse shear strain energy from the parabolic
distribution can be calculated by


Us = (1/2)V (zx 2/G) dV


= (1/2) [((3P/2bt3)(t2 4z2))2 / G] dV


= [(1/2)(3P/2bt3)2]/G (t2 4z22 dAdz


= (area)[(1(/)(3P/2bt3)2]/G (t2 4z2)2 dz


= bL[(1/2)(3P/2bt3)2]/G (t2 4z22 dz

= 1.2(P2L/btG)/2









While the transverse shear strain energy from the constant

distribution is


Us = (1/2) I(rx2/G) dV


= (1/2) ((P/bt)2 /G) dV


= [(1/2)(P/bt)2/G] (btL)

= (p2L/btG)/2

This result suggests the view that a uniform stress

P/bt acts over a modified area bt/1.2, so that the same Us

results. Therefore the deflection at the free end of a

cantilever beam with parabolically distributed transverse

shear stress will be 1.2 times that with constantly

distributed transverse shear stress.

The Mindlin plate is generalized from the Mindlin beam.

Thus the higher transverse shear stiffness from the

assumption of constant transverse shear stress has been

reduced by dividing the corresponding elastic constants by

the factor 1.2 for flat plate element. The reduced stiffness

will produce the more flexible response in shear that is

expected from the actual parabolic distribution.

The integration over volume must be introduced for the

calculation of the element stiffness and equivalent nodal

loads of distributed loads or temperature effects. As the

thickness of the element is constant, the integration over









volume can be changed to the integration over area. This can

be accomplished through the partition of the strain-nodal

displacement matrix [Bi] as follows:



0 0 zfj,x
0 zfi,y 0
[Bi] = 0 zfi,x zfi,y

fi,x 0 fi
fi,Y fi 0


The submatrices are named as follows:



BiA ZBiA
[Bi] =
BiB BiB


The element stiffness can be calculated through the
T
volume integration of Bi (3,5)E(5,5)Bi(5,3). Thus the [E]

matrix is to be partitioned as follows:



Ell E12 0 0 0
E21 E22 0 0 0
[E] = 0 0 E33 0 0

0 0 0 E44 0
0 0 0 0 E55



The submatrices are named as follows:


E EA 0
E =
0 EB










Then the stiffness of the element is



FT B -T T [EA 0 ziA 1
[Ki] =I B E Bi dV = zBiA [ A B dV
0 EB BiB

2 -T T
= [ BA EA BA] + [BB EB BB]


As [B] and [E] are constant about z, the integration

through the thickness from -t/2 to t/2 can be performed on z

only and yields


= T 3
[Ki] = [ BA(t3/12)EABA + BB(t)EBBB] dA


I -T- -
= [ BAEABA + BBEBBB] dA


where,
EA = (t3/12)EA and EB = tEB

Then [Ki] can be rewritten as matrix equation as

follows:



-T T EA 0 BiA
[Ki] = [ BiA BiB A dA
0 E B B iB d


-T- -
= BiE Bi dA









The size of plate element stiffness will be 27 by 27.


[K] =


[ E ]
(5,5)


[ B1 B2 B3 ... Bg Bg ] dV

(5,27)


B8
B9J

(27,5)


The strain-nodal displacement matrix from which

the constant thickness is taken out is defined as [Bi].


[Bi] =


0
0
0

fi,x
fi,Y


0
- fi,y
- fi,x

0
fi


fi,x
0
fiy

fi
0


BiA]
BiB


Equivalent nodal loads due to body forces on the plate

element are calculated as


Pb = fTb dV = 1
V -1 -i


fTblJI dr ds


in which {b} = { 0 0 bz )T or ( 0 by 0 }T or { bx 0 0
T in accordance with the direction of the gravity in the

coordinate system used. The nonzero quantities bx, by, or bz

represent the body force per unit area in the direction of

application.









Equivalent loads caused by initial strains are



PO = BTEeo dV
JV
1 1
= BTEO#IJI dr ds


where,

{40)= ( xx0 Oyy0 xy 0 0 )T
= { aAT/2 aAT/2 0 0 0 T


The stresses can be calculated from the equation


[a] = [E][E]


The corresponding generalized stresses, if desired, may

be computed from



M = ( Mxy My y QM Q y T = E ( B q 40 )

It is noted that the generalized stresses are actually

moment and shear forces applied per unit length of the edge

of the plate element. Therefore these can also be turned

into common stresses using the formulation for the bending

stress calculation. The moment of inertia for the unit

length of the plate is t3 / 12. Then the in-plane stress at

a point along the thickness can be calculated as


a = Mz / I = M(t/2) / (t3/12) = 6M / t2







91


The transverse shear stresses can be found as


r =Q/ t

But this may be multiplied by a factor of 1.5 to get the

maximum shear stress at a point on a neutral surface because

the transverse shear stresses show parabolic distribution

while the calculated stresses are average stresses coming

from the assumption of a constant transverse strain along

the element z axis.













CHAPTER 5
NONLINEAR SHELL ELEMENT


5.1 Introduction



The nonlinearities included in the formulation of the

Mindlin flat shell element are those due to large

displacements and those due to initial stress

effects(geometric nonlinearity). The large displacement

effects are caused by finite transverse displacements. These

effects are taken into account by coupling transverse

displacement and membrane displacements. The initial stress

effects are caused by the stresses at the start of each

iteration. These stresses change the element stiffness for

the current iteration. These effects are evaluated directly

from the stresses at the start of each iteration and are

included in the element stiffness formation.

The total Lagrangian formulation is used. If the

updated Lagrangian formulation is used, the element

coordinate system cannot be easily formed for the next

iteration because the deformed shell is not usually planar

[26]. The symbol {} is used for a column matrix (a vector)

and the symbol [] is used for a matrix of multiple columns

and rows throughout the chapter.









5.2 Element Formulation


The generic displacements of Mindlin type shell element

are translational displacements {u v w)T and denoted as {U}.

The displacements and rotations at a point on the midplane

are (uo vo wo 9x 0y)T and denoted as (Uo}. The generic

displacements can be expressed in terms of the midplane

displacements and z as

u = Uo(x,y) + zOy(x,y)

v = Vo(x,y) zOx(x,y) (5.1)

w = Wo(x,y)

The linearized incremental strain from Eq. (2.17) is

e = ( ui,j + uj,i + uk,i uk,j + k,j k,i)
(5.2)

This equation can be written out for the strain terms to be

used for shell element using the generic displacements {u,

v, w)T


exx =
eyy=

exy =


u,

U,y
i (
2'


ez = i (



eyz = (


+ tu u, + tv, v, + tw,x ,

+ t,y U,y + t,y V,y + tw,y Wy

U,y + V,x + tU, U,y + t,x V,y + tx w,y

+ u, tuy + V,x tvy + W,x ty )

,z + w,x + tU u,z + tVx ,z + tW,x w,Z

+ u,x tu,z + V,x tv, + ,x tw,z )

V,z + w,y + tu,y u,z + t,y v,z + tWy w,z

+ u,y tU, + V,y tV, + w,y tW, )

(5.3)










The derivatives of inplane displacements u and v with

respect to x, y, and z are assumed to be small and thus the

second order terms of these quantities can be ignored

through von Karman assumption from Eqs. (2.44) [20, 21].

Furthermore the transverse displacement w is independent of

z for the shell element which means that w,z is zero. Then

Eqs. (5.3) can be reduced to

xx = u, + tw,x w,x

eyy = Uy + ty Wy

exy = ( Uy + Vx + tw,x Wy + ,x wy )

exz = ( u,z + wx )

eyz = 1 ( v,z + W,y )

(5.4)

The incremental Green's strains, sometimes called

engineering strains, can then be shown as

x = exx = Ux + tw,x w,x

4e = eyy = Uy + tw,y Wy

Ixy = 2exy = Uy + v,x + tw,x Wy + Wx ty

7xz = 2exz = u,z + wx

yz e = e = + Wy (5.5)
It is noted that the linearized nonlinear strains are left

only for inplane strain terms.

By substituting Eqs. (5.1) into Eqs. (5.5), the Green's

strain can be expressed in terms of midplane displacements.




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs