MSU

LIBRARIES
m

 

 

RETURNING MATERIALS:
Place in book drop to
remove this checkout from
your record. FINES wiII
be charged if book is
returned after the date
stamped beIow.

 

 

 

 

 

 

STUDIES ON THE BENDING
OF ELASTIC PLATES

BY

Thomas C. Assiff

A DISSERTATION

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

DOCTOR OF PHILOSOPHY

Department of Mathematics

1982

 

ABSTRACT

STUDIES ON THE BENDING
OF ELASTIC PLATES

BY
Thomas C. Assiff

The bending of thin elastic plates with clamped
edges is studied under two related theories. The
classical theory. resulting in the familiar biharmonic
boundary value problems, is well-known. The so-called
improved theory of Timoshenko and Reissner differs from
the classical theory by; taking into account the effect
of shear deformations. A set of three linear elliptic
partial differential equations of the second order

results.

In addition to showing the existence of solutions
to the problem formulated in the improved theory, a
detailed analysis is made to establish the relationship
between the two theories, in terms of a single small
Parameter e. Standard functional-analytic methods,
along with a perturbation technique of Babuska. result
in integral estimates comparing solutions for the clamped

plate problem under the two reapective theories.

Thomas C. Assiff

The feasibility of numerical approximations to
solutions using the finite element method in the improved
theory is also examined. It is shown that accuracy

is adversely affected by small values of a.

Finite elements that are customarily used to obtain
approximate solutions in the classical theory are limited
by the requirement that they have continuous first
derivatives. The improved theory may be regarded as
one of the so-called penalty function methods, by whixfla
the smoothness requirements on the finite elements may
be relaxed. The behavior of the finite element solutions
in the improved theory for small values of e is
studied, as is their usefulness in approximating solaxtions

of the problem in the classical theory.

Additional asymptotic analyses are made, which
serve to illustrate the above theoretical results.
Numerical examples are presented, using piecewise linear

elements for problems for clamped beams and plates.

 

ACKNOWLEDGEMENT

The author wishes to thank Dr. David Yen, of
the Michigan State University Department of Mathematics
for the expert guidance and inspiration he provided
throughout the process leading to this thesis. His

confidence has been essential.

Appreciation goes also to Mrs. Cindy Lou Smith
for her superhuman effort, against all odds, in typing

thiS‘work.

A special thanks goes to the author's parents,
Mr. and Mrs. John Assiff for their unwavering support

and infinite patience.

ii

TABLE OF CONTENTS

AC MOWLEGMM. O O O O O O I O O O O O O 0

LIST OF TABLES O O O O O O O O O O O O O O

L I S T 0F F IGIJRES O O O O O O O O O O O O O 0

1.

INTRODUCTION . . . . . . . . . . . . .

l-l. Statement of Problems and
Background . . . . . . . . . .

1-2. Organization of the Dissertation
1-3. Notation and Function Spaces . .

1-4. Equations of Equilibrium for
a Clamped Plate in the Improved
Theory 0 O O O O O 0 O O O O O

EXISTENCE, CONVERGENCE AND FINITE
ELEMENT APPROXIMATIONS FOR SOLUTIONS OF
PROBLEM ( I ) o o o o o o o o o o o o o o

2-1. Existence of Solutions to Problem

2-2. Convergence of Solutions of
Problem (I) to Solutions of
PrOblem (C) o o o o o o o o o o

2-3. Convergence of Solutions of
Problem (I) to those of Problem
for the Beam and Circular Plate
with Axisymmetric Loading . . .

2-4. A Finite Element Theory for
Problem (I) . . . . . . . . . .

iii

Page
ii

vi

19

29

29

45

48

56

TABLE OF CONTENTS (cont'd.)

Page
3. BENDING OF A CLAMPED TIMOSHENKO BEAM . . . . 70
3-1. Asymptotic Expansion for the
Solution to the Clamped Beam
Problem in Timoshenko's Theory . . . 70
3-2. Principal Error for the Clamped
Beam in Timoshenko's Theory . . . . . 80

3-3. Construction of the Element
Stiffness Matrix for the Clamped
Timoshenko Beam . . . . . . . . . . . 103

3-4. Numerical Results for the
Clamped Beam . . . . . . . . . . . . 112

4. BENDING OF CLAMPED PLATES IN THE
IMPROWD MORY O C C O O O O I O O O O O O O l 2 5

4-1. Asymptotic Analysis of the

Solution Us . . . . . . . . . . . . 125
4-2. The Clamped Circular Plate . . . . . . 130
4-3. Principal Error . . . . . . . . . . . . 142

4-4. Construction of the Stiffness
Matrix for the Clamped Plate in
the Improved Theory . . . . . . . . . 154

4-5. Numerical Results for the
Clamped Square Plate . . . . . . . . 168

5. CONCLUSIONS AND FURTHER REMARKS . . . . . . . 186

BIBLIWMPW O O O I O O O O O O O O O O O O O O 1 92

iv

LIST OF TABLES

Table Page

1. Errors for beam problem (I),
pm = Sinh 2x 0 O O O I O O O O O O O O 118

2. Errors for beam problem (I), extrapolated,
pm = Sinh 2x 0 I I O O I O O I O O O O 119

3. Errors for beam problem (C),
pm = Sinh 2x 0 O O I O O O O O O O O O 120

. Errors for beam problem (I). p/D = 1 . . . . 121
. Errors for beam problem (C), p/D = l . . . . 122
6(X). o o 123

4

5

6. Errors for beam problem (I), p/D

7. Errors for beam problem (C), p/D = 6(x). . . 124
8. Square plate deflections, polynomial load . . 172
9. Square plate deflections, uniform load . . . 173

10. Square plate deflections, point load . . . . 174

LIST OF-FIGURES

Figure Page
1. Piecewise linear "roof" function . . . . . . 104
2 Nodes for local difference Operators . . . . 143
3. Type 1 triangular element . . . . . . . . . . 154
4. Type 2 triangular element . . . . . . . . . . 155
5. Support for "tent" function . . . . . . . . . 156
6. Center deflections extrapolated to e = O 185

vi

CHAPTER 1 - INTRODUCTION

l-l. Statement of Problems and Background

This thesis concerns the bending of thin elastic
plates, studied under two related theories. The
classical theory of plate bending is usually attributed
to G. Kirchhoff and A.E.H. Love [12,14]. The "improved"
theory, develOped by E. Reissner and R.D. Mindlin [19.16].
is derived from Timoshenko's beam theory [23]. It takes
into account the effect of shear deformation, which is
neglected in the classical plate theory. The classical
plate theory can then be regarded as the limit of the
improved plate theory as the shear rigidity becomes
infinite.

Throughout this work, a plate with the clamped
edge conditions is assumed. The governing boundary value
problems are compared, and a study is made to determine
how solutions to the problem in the improved theory tend
to those in the classical theory, in the limit mentioned

above.

The feasibility of using the finite element method
to obtain approximate solutions in the improved theory

is then studied.

Approximating solutions by finite elements in the
classical theory is hampered by differentiability
requirements not imposed in the improved theory} The
question then addressed is whether finite element
approximations constructed using simple elements such
as piecewise linear ones, can be used in the improved
theory to approximate the solutions to the problem in
the classical plate theory. In this way, the improved
theory is used as a penalty function method approach
to approximate the solution in the classical plate theory,
thus avoiding the construction of complicated finite
elements needed to satisfy the differentiability

requirements.

The two boundary value prOblems are cast in their
strong and weak forms, and functional-analytic methods
are applied to show existence of solutions and to derive
error estimates in terms of integral norms. Convergence
is demonstrated for solutions in the improved plate
theory to solutions in the classical theory by estimating
the difference, in terms of a small parameter c, which
measures the reciprocal of the shear rigidity. The
estimates hold for domains with smooth boundaries and their
sharpness is studied.. The problems of the clamped
beam and circular plate are studied as special cases,
for which improved estimates are presented, along with

asymptotic analysis and explicit solutions.

Estimates based on the approximating prOperties
of finite element spaces are used to guarantee convergence
of the finite element solutions to solutions of the
boundary value problem in the improved theory for fixed
6 > O. The dependence upon negative powers of e
observed in these estimates poses a practical computational
difficulty. This dependence is studied asymptotically

for the clamped beam and plate.

Using piecewise linear finite elements, a numerical
study of the clamped beam in the Timoshenko theory is
carried out. Solutions are then used to illustrate the
adverse effect of small 6, and a method is developed

for recovering reliable results.

A further numerical study of the clamped square
plate is offered as an example where the error estimates
may not strictly hold due to the presence of corners.
Again, piecewise linear elements are constructed within
the improved theory, and the problem is solved for
constant, variable, and point loads. These results are
analyzed and conjectures are offered toward explaining

some of the difficulties that they present.

The well known governing equation for the plate
in the classical theory is an inhomogeneous biharmonic

equation

and the clamped edge requirement provides homogeneous
Dirichlet boundary conditions. The variational form of
the biharmonic equation, from which the finite element
approximations derive, involves energy integrals of
squares of second derivatives. For the energy to be
finite,the finite element approximations must be
continuously differentiable. This is a cumbersome

restriction.

The governing equations for the same clamped plate
in the improved theory of Timoshenko and Reissner are
expressed as an elliptic system of three coupled second
order partial differential equations, along with
homogeneous boundary conditions. The boundary value
problem will be stated precisely in l-4c The system
depends on the small parameter c, in such a way
that as e» O, the biharmonic governing equation is

recovered.

The variational form of the clamped plate problem
in the improved theory, as derived in 2-1, imposes the
less restrictive condition that a finite element
approximation need only be continuous. This reduced
restriction on the finite elements is at the cost of
accuracy, which is adversely affected when the shear

rigidity is large (e small).

One of the objectives of this thesis is to

attempt to use simply-constructed finite elements in

5
the improved theory with small parameter c, to
approximate solutions in the classical theory thereby
avoiding the restrictive requirements usually imposed
by the latter. This indirect approach for applying "non-
conforming" elements to the biharmonic problem can be
viewed as an alternative to satisfying the "patch test"
(see [21]). The cost, again, is that small 6 requires
excessively small mesh size h in the finite element

construction.

1-2. Organization of the Dissertation

The remainder of this chapter is organized as
follows. Section 1-3 includes notations and definitions
of function Spaces, along with an informal discussion of
the finite element method. 1-4 introduces the equations
of equilibrium for the clamped plate in the improved
theory. Comparison is made to the classical theory,
and a physical interpretation is offered for the
parameter c. Problem (C) and problem (I) are formally

stated.

In section 2-1, weak forms of the boundary value
problems are derived and the existence of solutions
is shown for the improved theory, making use of the
Lax-Milgram theorem, a standard Hilbert Space result.

The variational property that a solution minimizes its

6

potential energy is stated and easily proved, although
this is also a standard result of the calculus of

variations [15].

Section 2-2 quantifies the relationship between
the classical and improved theories, by demonstrating
convergence, as 3 tends to zero, of solutions to the
clamped plate problem in the improved theory to those
in the classical theory. Direct integral estimates are
derived, verifying the energy estimate claimed by
‘Westbrook [25]. The sharpness and implications of

this and related estimates are discussed.

Section 2-3 contains improvements on the estimates
of 2-2 for the special cases of the clamped beam and
the axisymmetrically loaded clamped circular plate.
These estimates are derived by a method introduced by
Babuska [:3]. These results are contrasted with those
of 2-2, and compared to asymptotic behavior discussed

in Chapter 3 and 4.

In section 2-4, estimates are provided, following
again the procedure of Babuska, to show convergence of
finite element approximations to solutions of the clamped
plate bending problem in the improved theory. These
estimates show the presence of e-1/2 multiplying the
mesh size h (for piecewise linear elements) in the

error bounds, implying that accuracy for small e may

require an excessively fine mesh. Finally estimates

7

are also provided to show the convergence of finite
element approximations, formed in the improved theory

for small e. to solutions in the classical theory.
Again, error bounds contain terms proportional to e-l/Qh
(for piecewise linear elements) as well as terms with
positive powers of e, suggesting that convergence

occurs only when 6 and h both tend to zero, in a

somewhat restricted way.

Chapter 3 is a study of the clamped beam in the
Timoshenko improved theory. In section 3-1, a regular
perturbation expansion is deve10ped for the solution
in powers of e. For symmetric loads, the expansion
actualLy truncates, showing that the solution is linear

in 6.

Section 3-2 studies the discretization error due
to the use of piecewise linear finite elements on the
Timoshenko beam problem. Consistency between the finite
element system and the boundary value problem are shown
to depend adversely on the small parameter 3. An
observation of the structure of this adverse e

dependence leads to its removal.

Section 3-3 presents the construction of the
finite element stiffness matrix associated with piece-
‘wise linear elements, in terms of the element stiffness
matrix which is constructed explicitly from the potential

energy functional.

8
Section 3-4 contains numerical results for the
clamped beam under constant, variable and point loads.
Effectiveness in approximating exact solutions in the

classical and Timoshenko theories is analyzed.

Chapter 4 is a study of the clamped plate in the
improved theory. The difficulty of deriving an
asymptotic series in e for the general case is

discussed in section 4-1.

Section 4-2 studies the clamped circular plate.
Explicit solutions are constructed using a constant
load and a non-axisymmetric load, '§'= cos 9. The
latter serves to demonstrate the sharpness of the main

estimates for general plates, exhibitng boundary layer

behavior in its twisting moment.

Section 4-3 is analogous to 3-2, containing analysis
of discretization error. Again, consistency with the
continuous problem is shown to depend on 6, through

powers of its reciprocal.

In section 4-4, the element stiffness matrix is
constructed as in 3-3. The programming involved in
assembling the global stiffness matrix is briefly

outlined.

Section 4-5 contains numerical results for the square
clamped plate under constant, variable and point loads.
The possible effects of dOmains with corners are

discussed. Kondrat'ev's fundamental work [13] predicts

9
.singular behavior in the shear forces at the corners of
a square clamped plate in the classical theory. Possible

implications for the improved theory are considered.

Chapter 5 contains a summary of the results and a
discussion of their implications, as well as conjectures

and suggestions for further work.

1-3. Notation and Function Spaces

Several different notations are used regarding
partial differentiation. Let u1 and u2 be
sufficiently differentiable functions. The partial
derivative of ul 'with respect to x is denoted by

any of the following equivalent forms

bu
1 = = (1.0)
EX "u1,x _ D u1 .
Similarly,
Bu
1 _ s (0.1)
y = ul'y — D ul 0

Likewise, higher derivatives are expressed in the

same fashion:

 

azu
l (2,0)
an ED u
ax2 1,xx l
3
a u
1 E - D‘l'2)u

 

u =
BXBYZ 1.XYY l

10

The general form involving the multi-index is

(C1 Ia ) ala|u
Dau = D 1 2 u =--——¥L- ,
l 1 01 02
ax By
where o = (al,c2) and |a| = c14-a2 .

Vector notation is also used when convenient.

u
Vu = l,x
l u
l,y
u
1 _
v u - ul'X-i-uz'y
2
vzu = u 4-u
l 1,xx l,yy
v ul = ul,xxxx+2ul,xxyy+ul,yyyy .

Vector functions are generally denoted by capital letters,

 

 

eogo
r-u j
T 1
U = [ul,u2,u3] = u2
[“3]
A

N denotes the outward normal to the boundary where
this can be defined and a/BN denotes the normal

derivative at the boundary.

11

In the discussions of the beam, similar notation
‘will be used although only one independent variable is

involved

a” El?
n
as
u
c
E

Q.

ll
£3

xi

In fact vector notation may be used to maintain the

analogy to the more general case

E? “3.. E.e‘lé?‘

V1.1

4
c:

II
«JG

and capital letters will denote vector functions of

two components

u
U=[uou]T= l o
l 2 u
2
It will also be understood that whenever sufficient

smoothness is lacking, the governing equations will

be interpreted in a weak sense.

The standard order of magnitude symbols will be

used, and interpreted in the following sense
A = onon)

if lA/bn| remains bounded as b 4 O ,

12

with all other parameters in expression A
held constant. Since more than one parameter will be
involved at times, the expression O(bo) will be used
instead of the ambiguous 0(1).

Consider the bounded open connected region Q in
R2 and its boundary BO. 6 represents the closed
region occupied by the undeformed plate, hence it is the
domain for the dependent variablesi transverse displacement
and angular displacements. Assume the boundary is such
that O is Lipschitzian, as defined in [23]. Roughly,
this means that, locally, the domain 0 is all on one

side of an, and that the region is free of cusps.

Let C°(5) be the space of all real infinitely
differentiable functions on 0 where all their derivatives
can be extended continuously to the boundary 60. Let
c;(fi) be the subspace of C°(5) consisting of all
functions with compact support in 0. Let HO(O) be
the usual Hilbert space L2(O), i.e. the space of

square integrable functions on O, with norm given by

\M@=§§¥m
Q

and associated inner product

(u,v)O = If uvdA .
0

Let 1.2 1 be an integer.

l3

Define H‘(O) to be the closure of C"(§) under

the norm

nun: = 03%|“ “Danni,

Likewise, define Hé(n) to be the closure of C3(5)

under the norm H-Hi.

Define the semi-norms

|u|i = 2) ”Dan“: , for 1.2 1 integer
0|=£ .

Product spaces and their norms are defined in the usual

way:

3
(H‘) t t L

‘with norm defined by

Iluui= nulni+nu2ui+uu3n§. where U: u e (363 .'

 

Analogous definitions will hold for
1. 3 1. 2 1. 2

A convenient space related to the potential energy for

the clamped plate in the improved theory is defined for

each fixed 6 > O. For

l4
_ T 1 3
U — [111.u2ou3] 6 (HO) I

define

2

 

Hun: = IUliM'l [a] («13+ (3:) a. .

HUHe is clearly a norm and the following theorem shows
3
it to be equivalent to HUHl on (H1)

0
3
that He and (Hg) consist of the same set of

. It follows

functions.

Theorem: The norms H-Hl and H-He are equivalent.

In fact, for domain 0 ‘with largest dimension unity,

éuvui g Mi 3 (1+ 2e‘1>nvui .

Proof: Lemma (c) from section (2.3) provides that

for u 6 H3 ,

nun; 3%; umxug
Hung 3% umyng

SO

Hun:g%<11u.x\\§+uu.yu§) .

Let V = [v v v ]T 6 (H1)3
1' 2' 3 O '

Then 15

2 2 2
“le1 ano+ m1

2 2 2 2
||v1l|o+ uvzno+ uv3no+ Ivll

and by Lemma (c).

1 2 2 l 2 2
.<. z<l|v1,xllo+llvl,yl|o)+z(llv2,x||o+llv2,yllo)

1
+ quv )+ M:

2 2
3mm, +11%,“O
5 2 5 2
g; Iv:l g; “vue .

Also,

2 2 -
”V“e IV|1*'€

1 V1 2
'vv34'(v2) ”O

 

V1

V2

4.

 

II [V1
3 O v

2 O

 

 

 

2 -1 2
s |V|1+e (llvv3||o+2||vv

 

 

 

5)

2 -1 2 2
S |V|1+e (|V|1+2|V|1HVHO+IIVHO)

2 -l 2

|V|1+e <|Vl1+HVHO>
2 -l

s|V|1+Ze <Iv|§+uvu§>

g (1+ 26‘1)uvuf .

 

16

In discussing the finite element method, it will
be necessary to work with certain finite dimensional
subSpaces of Hk, ‘which are generated by bases made up
of finite element functions. These subspaces are defined
in terms of their ability to approximate functions from
infinite dimensional spaces. In practice these spaces
are generally made up of piecewise polynomials with
certain degrees of continuity. The domain in question
is subdivided into simple regions (say triangles or
rectangles) called finite elements, and a space of
functions is defined as all functions which are polynomial
of a prescribed degree over each element, and which
satisfy certain prescribed continuity conditions across
the boundaries between elements. The space is then
shown to be generated by a finite basis consisting of
functions (also refered to as finite elements) which
are themselves piecewise polynomial, and which have the
addition prOperty that they are zero except over a

certain few adjacent elements of the domain.

Although finite element spaces are occasionally
non-polynomial in nature, these are by far the exception,
and then are usually an attempt to deal with difficulties
(such as singularities) in specific problems. The
definition, which will follow, is taken from Chapter 4
of [ 3]. It and related ones omitted here, provided

the common thread which ties all finite element

1?
generated subspaces: namely, a formal expression of a
finite dimensional space's ability to approximate

functions belonging to an infinite dimensional space.

In the following definition, the parameter k
represents the degree of continuity (in a mean square
sense) present in the functions belonging to the finite
dimensional spaces being defined. In the context of

piecewise polynomial spaces, t. indicates the degree

of the polynomials present over each element (t 2

means linear, t = 3 means quadratic, etc.)

The parameter h is the mesh size. That is, each
element in the subdivided domain has a largest dimension
or diameter. h is the maximum of these diameters over
all subdivisions. Finite element subspaces are usually
constructed so that the mesh size can be varied. Each
different mesh size leads to different subspaces. The
finer the mesh size (smaller h) the more subdivisions
are necessary, and the higher the dimension of the
corresponding subspace becomes. As the dimension of the
subspace increases and the mesh size decreases, the
ability to approximate a particular function improves.
As h tends to zero, the system of subspaces should
provide a sequence of functions which tend to the desired
one. The definition which follows expresses an expectation
of this approximation prOperty; an expectation which is met

by standard types of finite element spaces, such as

18

piecewise polynomials defined over regular subdivisions
such as triangles or rectangles.

As defined in Chapter 4 of [ 3], let Sfi'k(0)

denote any linear system of functions with the following
prOperties: for t > k.2 O,

i) sfi'km) c Hkm)

there exists ¢ 6 Sfi'k(0) such that

u
ug-cvus s eh ngn, .
where U = min(t-s,z-s) and
C is independent of g and h .

Such a system is called a (t,k)-system. In the present
1

context, a (t,k)-system will be a subspace of H0 or

H]' spanned by finite element base functions, constructed
over an apprOpriately subdivided domain. Primary
importance in this thesis is given to the Spaces generated
by piecewise linear base functions ("roof" functions)
constructed over a domain subdivided by right triangles
formed from squares, all bisected by the diagonal with

a common direction. For these piecewise linear base

functions, t = 2 above.

19

1-4. Equations of Equilibrium for a Clamped Plate in

the Improved Theory

Consider a thin homogeneous isotropic plate of
arbitrary shape, with its largest lateral dimension
unity. The complete list of assumptions which comprise
the classical theory, sometimes known as "small
deflection, thin plate" theory, can be found in [22],
and others. The assumptions used to derive the improved
theory are the same, except for the deletion of the

assumption:

"Lines originally perpendicular to the
center surface of the plate before
deformation remain perpendicular to

the deformed surface after“.

The deletion of this requirement allows the inclusion
of the shear deformations, which characterize the

improved theory of plates.

The equilibrium equations for plate bending

under the improved plate assumptions are:
2 2 B I OW _
2[(1-u)v ¢x+<1+m§x1 -G H(¢x+—ax) - o
D 2 69 I OW
LAJ. -— l- + 1+ -GH +-— :0

G'H(v2w4-¢) = -p

where

 

a

20

vertical displacement
angular displacements (rotations) due

to bending

shear modulus

a constant introduced by R.D. Mindlin [16]
to bring about agreement between wave speeds
of short transverse waves obtained from
two- and three- dimensional analysis of
infinite plates.*

GH3

-€TI:ET = bending modulus

Poisson's ratio, 0 g_u g 0.5
plate thickness
transverse load over surface of plate

function of p, varying roughly linearly
2

from n2 = .76 at u = o, to n = .91 for u= 0.5.

2.1

 

2
Defining a = H 2 , equation (1) becomes:
6(l-u)x
l 2 -l
(l-4.2a) -2-[(1-p.)v ¢x+ (1+ 10%] -e (wx'Fg—X) = O
1 2 an _ -1 a_w =
(1-4.2b) 2[(1-11)v(:Y+(1+11)ay] e (¢y+ay) o
(1-4.ZC) €—l(V2W+¢) = -g
where
a a
¢._._‘”_x+_‘"z,
BX By

Equations (l-4.2) for the "Improved" plate theory
may be compared with the equation for classical plate

bending. The classical equation is:

(1-4.3) v4w=§.

In order to compare, let the load p be at least
twice differentiable. Differentiating (1-4.2a) with
respect to x, and (1-4.2b) with respect to y, and

adding, gives
v2¢- e-l(¢+ VZW) = 0

with (1-4.2c) this becomes

(1—4.4) v24) = -§ .

22

Now taking the Laplacian of (l—4.2c), combining with

(1-4.4)o

(1-4.5) v4w=§-e v22 .

With the angular displacements eliminated, (l-4.5)
represents an equation to determine the transverse
displacement under the improved plate theory assumptions.
The striking feature of (l-4.5) is its similarity to the
governing equation for the plate under the classical theory,
equation (l-4.3), which neglected shear affects.
Equation (l-4.5) merely has the additional term -eV2 5.
due to inclusion of shear effects, which vanishes as
6 tends to zero. Since the classical governing
equation is recovered as the limit as 6 tends to zero,
the inference is that e is a measure of the effect of
including shear. The solution ‘w governed by (l-4.5)
is still coupled to *x and wy through boundary
conditions, but assuming the boundary conditions are
chosen to be compatable between classical and improved
problems, it is reasonable to expect the displacement
w produced by the improved theory to approach the
displacement produced by the classical theory, as 6
tends to zero. This convergence will be examined for

the clamped plate boundary conditions in both an integral

sense, and asymptotically in a pointwise sense.

23

A physical interpretation of this convergence
is worth noting. The process of convergence as 5 tends
to zero may be viewed within the context of a fixed
material (meaning u, x, and G are fixed). Then
6 = H2 2 tending to zero means that the plate
6(1-u)n
thickness tends to zero. Consider the transverse load
p being scaled also, in such a way that ‘5' remained
constant, independent of plate thickness H. That is,
let p be prOportional to H3. It must be remembered
that in both models, classical and improved, all
in-plane forces are neglected so the limiting process
above should not be viewed as leading physically to a
description of the "membrane" problem. Rather, as the
thickness decreases, intuition suggests that the
transverse load is carried less and less by the shear
stresses and more and more by the bending moments. In
the limit there is no shear, and the improved plate model

becomes identical with the classical plate theory of

bending.

It will be useful to express equations (l-4.2)

in operator form:

(1-4.6) LU

ll
'11

where

 

 

 

L = LB+ 6: Ls
' 2 2
1 2 _a__ 1 A a 1
'2’[(1-H)V + (1+H.)ax2] 2(1+H)BXOY O
B2 B2
I.B = "(Hmaxay 2[(1-u)v2+(1+u);—5y2] 0
L o o o -
1
" .31.
-1 O '-Bx
= _ __§_
LS O 1 By
1 _3_ v2
L Bx By -
r" r-
wfl 01
= p F = O c
U WY
2
w -
_ A _ D.J

 

 

 

 

 

 

 

Define, also the following related energy expressions:

T
for U - [u1,u2,u3]

Bw3
PS,(UW)=II[(—+ul?£-)( +w

O

+

T
and w = [wl.w2.w3]

l)

Bu Bw
3 3
Bjy + w:)]dA

By 4-112 )(—

(—

25

1 Bu1 Bu2
PB(U,W) = :- QX[(1+H)(W + W

Bw Bw
)(—1+-3)

Bu Bu Bw' Bw
+<1-u)<—3—,% - ——?-)(-—-1- - —-2—)

Bu Bu aw Bw
1 2 1 2
+ (l—p)( BY 4- 3X

T
PL(F.W) = -jflj‘wr dA= £f§w3 dA

_ -1

Ps(U!V) o

The quadratic forms derived from the bilinear energy
expressions have some physical interpretation.
e-lPS(U,U) represents the potential energy due to

shear deformation present in the solution U a [u1.u2.u3]T

where u3 gives the displacement and u and u

1 2
give angular rotations due to bending. Note that the
Bu Bu
quantities 7§?'+ u1 and .35; + u2 are respectively

the angular rotations due to shear in the x- and y-
directions. PB(U,U) represents the potential energy
due to bending in the solution U. PL(F,U) represents
the potential energy due to the applied transverse

load p. B€(U,U) represents the total energy due to
both bending moments and shear. If U represents

the solution to the classical plate bending problem,

 

 

 

 

 

26

since rotations due to

then U=[ 5!- Bw W]T.

'Bx' ‘3'
shear are zero. Consequently, PS(U,U) vanishes and

_ 2 2 2 2
PB(U,U) - H’ w'xx+ 2w,xy+w,yy+ 2p(w,xxW.yy-W.xy)dA
= ‘H‘ (w +w )2-2(l- )(w w -w2 )dA
'xx 'yy “L 'xx 'yy 'xy '

This is the strain energy (up to a normalizing constant)
due to bending for the classical plate [22]. Under
the clamped plate assumption, the boundary conditions

w = B_w_ = 0 cause the term H‘ w,

BN w, -w? dA to

xx YY xY

vanish. Then

' 2
PB(U,U) = U (w,xx+w,yy) dA .

The boundary value problems for the clamped
plate under load f = p/D can now be stated precisely.

The problem in the classical theory is

 

BW’ Bw
- .. __9. ___9 T
Problem (C). find UO - [ - BX , By , wO]
satisfying
v w = f on 0

w = O on BO

0 on B0 .

27

The problem in the improved theory is

 

 

. _ T
. Problem (I): find Ue — [¢x,wy,w]
satisfying, on O.
2 2

at 61
.1 2 x y 1 fig;
21(1-11)v (x+(1+u)(ax2 + axay” e (ixJ'ax 0

52¢ . a2
1 2 _X __x - -1 A! _
2[(l-u)v wy+(1+u)(axay+ ay2)] e: ((ty-I-ay O

Bw Bv
-1 2 .1 _x .._
e (vw+Bx+By)-f

and

Problem (C) and Problem (I) will be interpreted
in their weak forms when sufficient smoothness of
solutions is lacking. The weak forms, derived in 2-1,

are as follows. For F = [O,O,-f]T, f =-E .

 

D
Bw Bw
. _ o o T
Problem (C). find UO - [-—-ax , -—3y , wo]
2

where ‘w e H and

O O

PB(UO,V) = PL(F,V)

for all
— fl _

where

28

v 6 H

T

.. _ 13
Problem (I). find Ue — [ix'Iy’W] 6 (HO) where

B€(U€.V) PL(F.V)
for all

3
v 6 (H3)

CHAPTER 2 — EXISTENCE, CONVERGENCE AND FINITE
ELEMENT APPROXIMATIONS FOR SOLUTIONS OF PROBLEM (I)

Chapter 2 studies the existence of solutions to
problem (I) and their convergence to solutions of problem
(C), as the shear rigidity tends to infinity. The Special
cases of a clamped beam and a circular plate with
axisymmetric loading are studied as examples where
improvement in the estimates for convergence can be
achieved. Finite element approximations to solutions of
problem (I) are developed, and their behavior investigated

as the shear rigidity tends to infinity.

2-1. Existence of Solution to Problemg(I)

Let wb

plate problem

be the solution to the classical clamped

(2-1.1)

WO=W=O on an.

The existence of a unique solution w0 6 H3 is guaranteed

(see [ 8]), at least for :8 6 Ho and smooth domain (2,

problem (C) being an elliptic boundary value problem with

29

3O

homogeneous Dirichlet boundary conditions. For the

weak form of (2.1.1), let v e H: then.

0 = [I (v4wo-%)v dA
0
using Green's Identity

2 2 B Ba
{11‘ (av B-Bv a)dA =a£ (01 fi-B 'a‘fims

(2-1.2) O = H‘ (vzwovzv-g v)dA .
O

i In terms of the energy definitions of Chapter 1, (2-1.2)
2

states that for every v 6 HO, wo satisfies
where U = [-w -w ‘w ]T
O O,x' O,y' O
_ T
V - [‘Voxo‘vaylv] 0

(2-l.3) can be derived also as the Euler equation for

problem of finding the minimum over H2

O of the quadratic

functional

(2-1.4) JCLASS = PB(U.U) - 2PL(F,U)

Let Ue = [¢x,¢y,w]T satisfy equations (1-2.2) along

with the homogeneous boundary conditions. That is,

31

U6 = [0,0,0]T on B0 .
Then for any V = [v v v ]T 6 (H1)3
1' 2' 3 O '

(2-1.5) H (VTLUe-VTF)dA
o

= $05. {%[(1_p)v2¢x+ (1+u)§:'1 -e-l(¢x+%-¥)}Vl

1 _ 2 m _ -1 a.
+ {21(1 u)v ¢Y+(1+u)ay] e (¢Y+ay)}v2
-l 2 _
+ {e (v W+¢)o+§]v3 dA — o .
Integration by parts yields
1 1 ix av1
If -‘2'(1-u)V\)x ° vv1--2-(1+u)v ° (F) 33"-
' O Y
-1 Bw
- c (¢x+-5_§)Vl
l l wx 3V2
-7(l-u)v~)y - sz -§(1+u)v - (1y) 3;-
-l Bw
- (WY+?§)V2
-l -l wx
-€ W'VVB-e (v)-vv3+Dv;dA=o
Y

.1 1 ¢x\ vl
In] -2(1-u)[v¢x - vv1+va - vv2] --2-(1+u)v - (WY) v - (v2)

«I(<::)+w)-<(::)-wa>w>a=o

32
BUX BV

 

__l_ _1_ _X_
(2-1.6) XQJ‘ {'2(1'P~) TX OX + by by + 3x BX
Bt sz
BY Y

aw BN Bv Bv

.; X By x1+ 2
’2‘1+“)(—ax+ HM +ay)}dA

-1 _

New consider

th Bv2 wa sz w Bv1 t Bv1

(2'1'7) if ( Bx By *' By Bx By x Bx y

by integration by parts

B¢ BV Bt B¢ Bv B t
= If (_ x v2 ___£L.V _ __XZ__l._ ___X.V )dA
Bx By ByBx 2 By' Bx BxBy 1

by integration by parts

B¢ Bv BW Bv B
If <-3%a%+tfa?2'?%7+#7m

Inserting this into (2-1-6), it becomes

-1 _
-PB(Ue ,V) -€ PS (U6 ,V) + PL(F,V) - O

-1 _
(2-l.8) PB(U€,V)+S PS(U€,V) - PL(F,V)

1 3
v V 6 (Ho) .

33

(2-1.8) can be written
_ 1 3
(2-l.9) Be (Us ,V) - PL(F,V) V V 6 (HO) .

(2-1.8) or (2-l.9) can be derived as the Euler equation

3
for the problem of finding the minimum over (H3) of

the quadratic functional

_ -1
(24.10) JIMPRDWE _ PB(U,U) +6 PS(U,U) -2PL(F,U)
or
(2-1.11) Jmpaovan = B€(U,U) -2PL(F.U) .

T .
As before, let Uo - [‘wb,x"wb,y'wb] be the solution

of problem (C). That is, wb solves

v4w = E on 0

O
(2—1.12)
WC = 0
BW on B0 .
_0- = O
BN

Another weak form will be needed for U0. Consider,
3

_ T 1

34

. _ l - *
(2-1.13) PB(UO.v) — H 2(1+“)("”o,xx w0,yy)(vl,x+v2,y)
+--i(1- )(-w 4vw )(V -V )
2 H 0.xx O,yy 1.x 2,y
+ l'U-‘H-H‘w —w )(V +V )dA
2 0.xy O,yx l.y 2.x
v
= H v(v2w0) - (v1)dA .
2
by integration by parts .
4 Bi.
Now v wO-D - 0, so

4 E ..
If (v wO-D)v3dA — O 0
Integrating by parts here also.

If (-v(v2w0) ° vv3 -§ v3)dA = O

and
(2‘1 14) If -v(v2w ) ° (W +(v1)) -2 v dA
° 0 3 v2 D 3
v
= ~fj' v(v2wo) - Vila) dA .
Using (2-1.l3). (2-l.l4) becomes

V
Pswo'v) = U v(v2w0) ' (VV3+(V:))dA+H ‘5 V3 dA

or

35

(2-1.15) PB(UO,V)

V

_ 2 o 1
- PL(F,V)4-ff v(v WC) .vv34- v2 dA

1 3
v V 6 (H0) .

The existence of solutions for problem (I) is now
considered. The Lax-Milgram Theorem is used below to snow
tnat (2-l.9) has a unique solution, Ue’ for each e > O.
In fact, since B€ is symmetric, it represents an inner
product on He' In this context the Lax-Milgram theorem

reduces to the Riesz representation theorem.

Theorem 2 (Lax-Milgram): Given a bilinear form
B(U,V) defined on S)<S, where S is a Hilbert space

'with norm H°H. if:

i) There exists M > 0 independent of U and

v such that |B(U,V)| g MHUHHVH v U,v e 5 XS, and

ii) There exists y > 0 independent of V such that

B(V,V) 2 «AMP. v v e s,

then the equation
B(U,V) = F(V) , V V 6 S

has a unique solution U E S for every continuous linear

functional F defined on S.

The theorem is proved in Gilbarg and Trudinger [ 9].

36
A number of lemmas are required before Theorem 2
can be used to show the existence of a unique solution

to problem (1).

Lemma (A):

 

Hul'x+u2'y“0 S fl “”1
\\u1,y+u2,xuo g J5 lUll

Hul,x"u2,y"O‘S-J§ IUIl ’

my
nu”: u2.yllé s <uu1.xuo+ ““2.y“0)2
s 2‘““1.x”c23+ ““2.y“g)
s ZlUIi -
Similarly.
Hul'y+u2'xllg g 2|U|i

2 2
Hula: ‘uz'yno S 2|U‘1 '

Lemma (B):

HUHEHVHe 2 lU|1|V|1+ 6'1 JPS(U.U) A/PS(V.V)

37

Proof: For any real numbers a,b,c,d,

22 22

03 (ad-bc)2 = a d +bc -2 abcd

that is ,

22 2

ad +bc222abcd.

Adding 32c2+b2d2 gives

2 2 22

22 22 22 +bd+2abcd,

ac +bd+ad+bc22a2c

that is,
(a2+b2)(c2+d2) 2 (ac+bd)2 .

Setting a = 1n|1, c = |v|1, b = e'l/z ./pS(U,US,

d = e‘1/2 ./PS(V,V)
(|U|1|v|1+ 6-1/2 ./PS(U,US ./PS(V,V))2

g muff e'lps(u,um|v‘§+ e‘lps(v,V)) = nun: "v“:

1
O 0

Lemma Q2): For u E H
2 1 2
llullo s 7H!» xHo

Hung g %uu.yu§

38

Proof: Denote

D
II

{(XOY) E O, for x constant]

0

y {(Xay) E O. for y constant} .

Let (x.y) 6 0. Let (XOIY) 6 BO, where xO < x and

v g 9 xo<g<x. (§.Y)€O. then

X

u(x.y) = y -§§ (§.y)d§
x
O
_ x au
|u(X.y)| - If 33 (§.y)dg|
x0

x
g,f lgg-(g.y>|d§
x
O
and by Schwartz inequality,

1/2 1/2

.g (jx LEE (§.y )Izdg)
x0 ag o

|x-x0|

Now,
Hung = If IU(x.y)l2dy dx
0

X
s x“: Ix Iméfllz aux-xotdy dx
0

g_ff I |%%12d§|X-xo|dx dy
Q

39

 

Similarly
2 1 2
Hung sinuwuc,
. l l 1
Lemma (D): For u3 6 HO, and ul 6 H , u2 E H
and for all 0 < p < l.
2 2 2
If (u3'x4-u1) dA-Z (1-p) If u3.x W‘S I u1 1.xdA
O O O
2 2 2
Inf (“3,y+“2) 5A2 (1‘9) 1‘0.) u3,37‘11‘3 a.) n2 2,ydA

that is,

2 2 2 2
PS(U,U) 2 (1-p)|u3|l-E If (“1,x+u2,y)dA .

Proof:

U (“3,x+“1’2d“

-2 If (“i.xi'zu'u3.x)dA

If (u; x--2u3u1 x)dA (by integration by parts)

= (1-10) If “3.x GNP N “3.x “'2 U “3“Lx GA

2 <1-p> H nix dA+% n u§ da-z m u3u1,x «m

(by Lemma (C))

4O

2 4ui x
= (1‘9) U ‘13.di +2 EU (“3 %“3u1,x+' 2 ”A
2 2
“E.” u1.x dA
=(1-p)j‘fu§.x dA+523 ff (u3—% “1.x’2‘1A

2 2
"E if u1.x dA
2(1-P> II°§.di--'Hui -

Similarly

H <u3.y+u2>2dA2 <1-p) 11%.), CIA-€- U “iy dA-

Theorem 3. Equation (2-l.9), which is problem (I),
has a unique solution U6 6 36' For each e > 0 and for

0 3
each F E (H )

Proof: Clearly PL(F,V) is a continuous linear
1 3
0)
form B€(U,V) = PB(U,V)4-e

. It remains to show that the bilinear
-1

functional on (H
PS(U,V) is continous and
coercive (i.e. satisfies conditions i) and ii) of the
Lax-Milgram Theorem) on H€)<H€. The result then follows

from Theorem 1.

_ T _ T
Let U - [u1,u2,u3] and V — [vl,v2,v3] 6 He'

i) To show B€(U,V) is continuous on HexHe ,

41
lB€(U.V)|

l

= [PB(U,V) + e‘ PS(U,V)|

-1
g |PB(U,V) | + e |PS(U,V)|

l
"5 If (1+ p)(u1'x+ u2,y) (Vl,x+ v2.17)

+(1-p)(u u )(v1,x-V2 )

LX- .17

2,y
+-(1-u)(u1'y4-uz'x)(vl'y4-v2'x)dA|
+ 8-11“. (u3,x+ u1)(v3'x+ V1)

+-(u 4-v2)(u3'yi-v2)dA|

3,y
and by Schwartz inequality,

1
S 3(1 + MUM”; + uz'yllollvl'x + Vz'yuo

£11

2 "“13” ‘12 ,ynou V1.x - v

+ 2,yHO

+ ‘3”- " H)”u1.y+ “2.xu0“v1.Y+ Vz'xno
+ 3‘1Hu3 ,x+ u1|loUu3,x+ Vluo

+ e‘lnu3,y+u2uouv3,y+vzn 0
and by lemma (A)
—1
S (3-”)|U‘1‘v|1+€ “u3'x+u1"0\‘v3'x+vll‘o

+ 61111.13 ,3, + uznouv3 ,y+ vzuo

g (3 -p)[|U|1|v|1+ e'1./1=S(U.U) ./1>S(v.v)]

41
|B€(U,V)|

l

= IPB(U.V) + e” PS(U,V)|

-1
g [PB(U,V) | + e [PS(U,V)|

l
= b- ‘U (1+ u)(u1’x+ u2,y)(v1,x+ v2,y)
+(l"m(‘11,x-‘12,y)(vl,x"V2.y)

+ (1 — u) (ul'y+ uz'x) (vl'y+ V2.38“)

+ e-ll.” (u3,x+ u1)(V3,x+ V1)

+( +_v2)(u3'y+v2)dA|

113 'y

and by Schwartz inequality.

1
S 7(1 + H)‘|u1'x+ uz'yuonvl'x+ Vz'yuo

V

1..“
+ 2 Hul z'yuo

,x - u2 .ynollvl.x '-
+ -§-<1 - u>11u1,y+ u2,xl\oilv1,y+ VLXHO
+ e-1||“3 ,x+ u1"o”“'3 ,x+ Vlno

.+e‘1nu

3 .y+ u2|IOHV3 .y+ v2H o

and by lemma (A)

g <3 - w IU|1|V|1+ e'luu3,,.+ u11\onv3,x+ vluO

+e‘1u

u3,y+ 1511011va vzno

g (3 -u) [|U|1|v|1+ e‘l./Ps(U,U) ./1>S(v,V)]

42
and by Lemma (B) .

g <3 -u)||UHe||V|I€
g 3||U||€HVlle

Note: |B€(U.V)I g MHUHellvue where M is independent

of e and u aswell as U and V.

ii) to show B€(V,V) is coercive, that is

B€(V.V) 2 CHV": , for all v 6 (H1)3

0 where C is

independent of V and e.

B€(V,V)

1

PB(V.V) + e' PS(V.V)

)2

=%N (1+p.)(vl'x'I-vz'y)2+(1-(.J.)(v1'x-v2'y

+(1-u)(v1'y+v )2 dA

2.x

)ZdA

-1 2
+ e H. (v3'x+v1) + (v3'y+v2

-1 2 2
+6 I! (v3'x+v1) + (v3’y+v2) dA

l 2 2 2 2
2(1 '11) H. v1,x+ v2,y+ v1,y+V2,x

-2 +2v dA

v1,xv2,y l,yv2,x

-l 2 2
+ e If (v3'x+v1) + (v3'y+v2) dA

43
and by integration by parts,

_ 1 2' 2 2 2
- 2(l-H) If Vl'x+V2'y+V1'Y+V2'x dA

-1 2 2
-+e If (v3'x4-v1) 4-(v3'y4-v2) dA
__ _1_ _ . 2 2' 2 2
_. 2(1 p) j! vl.x+v2,y+vl,y+v2,x dA
-1 1 2 2
+ (e ‘18) H (V3.x+vl) +(V3.y+v2) dA
1 2 2
+18 U (v3'x+v1) +(v3'y+v2) (SA
and by Lemma (D), with p ='%
1 2 2 2 2
2-2-(1-11) If vl,x+V2,y+v1,y+v2,x dA
+ (3‘1.“%%) If (v3 x+vl)2+(v3 4-:2)2 dA
1
+ 36 If v3 x4—v§ -9 ff v1 x+ 2y dA
1 2 2 2 2 2
2 (7241-11) -3) if vl,x+v2,y+v1.y+v2ox dA
1 2 2
6 If V3,x*'v3,y’dA
-1 1 2 2
+ (e ’18).” (V3.x+vl) +(v3.y+v2) dA
35 . 1
andfor O<eS-2—, and51nce O<u<§
l 2 2 2 2 2 2
B€(V,V) 2 36 U vl'x+v2'y+v1'y+v2 x+v3 x+v3'y dA
+..€_-1 (v +v)2+(v +v)2 dA

-1

_._ 3% (M1. . pSwm) = 3% “V“:

44
Note that the constant of coercivity C = g%- for the
range of O < e < %§-. and since 3 is considered a
small parameter, this is more than sufficient. The
constant can be improved slightly by more judicious choice
of certain constants in the estimates. as well as allowing
C to vary with u. However these improvements are
slight. Any improvement will also be valid for a smaller
range of 6. Likewise. if necessary, coercivity may be

shown for large a, at the cost of reducing the constant

C.

Theorem 4: U6 minimizes the functional

J(U) = B€(U,U) -2PL(F,U) over He .

Proof: Since Ue satisfies B€(U€.V) =

PL(F.V) V V 6 He.

J(U) B€(U,U) -2Be (Us ,U)

B€(U-Ue,U-U€) -B€(U€,U€) .

Clearly J(U) is minimized by U = Ue' since B€(V,V)
is positive definite by virtue of property ii) (coercivity)

of Theorem 3.

45

2-2. Convergence of Solutions of Problem (I) to

Solutions of Problem (C)

From (2-1.15) and (2-l.9), the solution vectors
U0 and Us' for problems (C) and (I) respectively,

. l 3
satisfy for all V e (HO) .

(2-2.1) - PB(UO,V)

= PL(F,V) +j‘0j‘ v(v2wo) - vv3+ vi) dA

(2-2.2) Be(U€.V) = PL(F,V)

Since PS(UO,V) = O for all V 6 (H3)3 , (2-2.l) can

be rewritten as

(2-2.3) B€(UO,V)

v
_ 2 . 1
- PL(F.V)4-fj v(v WC) vv34- v2 dA .
Q
Subtracting gives
B(U—UV)- (vzw)° v+vldA
e e 0' _ If V O V 3 v2 '

for all v 6 (H3)3 .
Then

|B€(U€-UO,V)|

V
$.(ff vazwo)lzduzul/2 If Vv34, v: 2dA 1/2

46

wh ich is

(2-2.4) (Raw€ -U0.v)| g Hv(v2wo)\|o(Ps(V.V))l/2 .

1

From Be(V,V) = PB(V,V)4-e- PS(V.V) it follows

PS(V,V).g.e B€(V,V) so (2-2.4) becomes
(2-2.5) |13€(U€ -UO.V)I g 61/2IIV(V2wo)HO(Be(V.V))1/2 .

Setting V = U€--UO in (2-2.5) results in

(2-2.6) (136m€ 410.11e 'Uo”1/2 g el/ZHVWZWOHIO

or
(2-2 7) B (U -U U -u ) g eHv(v2w )u2
° 6 e 0' e O 0 0

which is Westbrook's energy estimate.

From (2-2.7). using coercivity and norm equivalence,

1/2

(2-2.8) HUS -UO”€ S C 6 "V(V2wo)no

(2-2.9) ”Us -UO||1 g c el/vawzwonlo .

where the constant C appearing above will be used in a
generic sense, being different in different contexts, but

being alway independent of e and the functions involved.

The estimate (2-2.8) is sharp, but the sharpness of

(2-2.9) is an Open question. It can certainly have no

power greater than CB/Z. As shown in 4-2, the solution

47
to the problem of the circular plate in the improved
theory with load p/D = cos 9 exhibits a term of order
C(Cl/Q) over a boundary layer region. This causes
HUquOH1 = 0(33/4), ruling out a dependence of US on
e which is analytic at e = 0. The 31/2 dependence

appears in the radial derivative of the rotation component

Bw aw
in the e-direction, that is, -—§-. -——& corresponds to

Br Br
the twisting moment. The presence of this dependence
prohibits improvement of 31/2 to e in (2-2.9), and

may be responsible for the lack of an easily constructed
asymptotic expansion for U6 for small 6. as discussed
in Chapter 4. It also may contribute to difficulties in
using numerical approximations to U€ in order to
approximate U0.

While (2-2.9) cannot be improved from 31/2 to e in
general, this can be done in two cases of special interest.
First, for the clamped beam, which may be considered a

limiting case of the clamped rectangular plate: and second,

for the circular plate with axisymmetric loading.

Estimates (2-2.8) and (2-2.9) are useful whenever
the boundary B0 is smooth. In that case, partial
differential equations estimates [21] show that for

w satisfying

0

48

(2—2.10) HwOH4 g cnfino .

Then

uv<v2wo>uo g cuwon3 g cuwon4 g cnguo < . .

If the boundary is not smooth,(2-2.10) is not guar-
anteed. and the presence of singular behavior of solutions at
corners is well known. For example, in the important
case of a boundary with right angles, it is precisely
the square integrability of the third derivatives, needed
for (2-2.8) and (2-2.9). which is jeopardized. The
corner difficulty is discussed further in Chapter 4, and
the possibility of extending (2-2.8) and (2-2.9) in
the presence of a domain with corners, is investigated

numerically in Chapter 4.

2-3. Convergence of Solutions of Problem (I) to those
of Problem (C) for the Beam and Circular Plate with

Axisymmetric Loading

The description of the clamped beam can be derived
from the equations and energy expressions for the clamped
plate by simply deleting all dependence on y, and making
the obvious corresponding modifications to the definitions
and proofs regarding function spaces, norms, etc. The
purpose of this section is to derive estimates analogous

to those for the general plate in 2-2, but with some slight

49

but important improvement. This improvement over the
general clamped plate estimates is also shared by the
axisymmetrically loaded circular clamped plate, and hence

it is included here also.

In what follows, the notation used in 2-2 for the
general plate will also be used here. For the beam,
it will be understood that v and v- mean -é§-, etc.

Thus (2-l.9) and (2-1.15) are to be interpreted as

B€(U€.V) = PL(F,V)

and
1/2 I! I
PB(UO,V) = PL(F,V)+_J'1/2wo (v2+v1)dx ,
where

aw
U = [ :] , ‘wm’= 51, ‘w = w” = O at x = $145!.

* V1 1 2
n.1,]. v1.1a -

For the circular plate with axisymmetric loading, the
notation used in 2-2 is valid in its general form. However,
the description of this case would ordinarily be presented
in polar coordinates, as is done in Chapter 4. For the
sake of the following estimates, this is not done, in

order first, to maintain similarity to the beam formulation.
and second, to demonstrate why this improvement fails

for the general plate, and to comment on its connection

with the lack of a satisfactory asymptotic eXpansion.

50

Consider the weak forms satisfied by Us and U0,

respectively, equations (2-l.9) and (2-1.15)

v
1
(2-3.1) B€(U€,V) = PL(F,V) . V V = V2 6 He ,
v3
where
_ -l
B€(U€.V) - PB(U€.V)+-e PS(U€.V)

(2-3.2) PB(UO.V)

v1

_ 2 v1

- PL(F.V) +‘H' v(v wo) . vv3+ (v2) dA

V V = v2 (EH
v

Following a method similar to that used by Babuska [3 ]

set
g1 n1
Ue = Uo-e§4-n , 'where g = g2 and n = n2 0
g3 n3

g o T] 6 H6
Then (2-3.2) becomes
B€(Uo-eg+n,V) = PL(F,V)
or, Since PS(UO.V) = O for all V 6 He .

PB(UO,V) -e Be(§,V)+B€(n,V) = PL(F,V)

51

Subtracting (2-3.2) gives

V
-e B€(§.V) +B€(n.V) = ~fo v(v2wo) - (vv3+ (v:))dA

V
Be(n.V) = e 1’]_,,(§..V)+P3(§.V)-fl~ v(v2wo) - (W3 9(v:))dA

O

§
(2-3.3) B€(n.V) = e PB(§.V)+fnj' ((v§3+(;:) -v(v2wo))
. v1
° vv dA .

(2-3 4) v 5 + g1 = V(V2W)
. 3 :2 0 .

Now choose g 6 He' so that

The existence of g (or lack of it for the general plate)

will be shown below.

(2—3.3) then becomes
B€(n,V) = e PB(§,V), for all V 6 He .
Now choose V = n. Then

B€(T]0T'|) = € PB(gon)
g e .FT‘TPB a; .F——TPB(n.n
S e Hglll damn)

so dividing by JB€(n,n) gives

52

(2-3.5) JB€(n:n) S 6 “g” -

It then follows that

(2.3.6) unne so a Ilélll

(2-3.7) HnHl g,c e ”EH1 . . where g 6 He
satisfying (2-3.4) is independent
of e. depending only upon ‘w

0 0

Now U€-U can be estimated using (2-3.6) and (2-3.7).

o
'with norm equivalence. Since Ue-UO = -e§4-n .
(2-3.3) HUS-U0“ = \l-egmlle

g 6H§||€+ HnHe g c 61/2 M1
(2-3.9) ||Ue - 0H1 = !!-e§+nH1 g eHEHl+HnHl g C e llglll

where g depends only on the classical

solution wo.

Although (2-3.8) has the same 91/2 dependence as
does (2-2.8). (2-3.9) shows the improvement over (2-2.9),

1/2

from s to c.

This improvement was uncovered by splitting the error

term and explicitly determinine one part. then using it to

53

estimate the other part more finely. This is analogous

to assuming an asymptotic series form for a solution, then '
explicitly solving for the consecutive terms, using each
one to derive the next. The success of this indirect
method of estimating U€-U hinges on the existence of

0
g 6 He' a solution of (2-3.4).

For the beam, the analog of (2-3.4) is

dgz dawb
(2-3.10) —+ g =——— .
dx ]. dx3

Along with the homogeneous boundary conditions :1 = g2 = 0

at x = :t%- assigned to all members of He , a solution

can be formed by defining

3
d‘wO dgz

g =—-—
1 dx3 dx

where 52 is the solution, for any function g(x), of

(14% ( )
dx4

§2=0 at x= 1%

3
dgz d.w
dx

ll

m

n

x

I

H-
Mhd

By the boundary conditions assigned to :2, it is clear
that §1 = O at x = :t% . One possibility is to choose

g(x) to be identically zero. However, if g(x) is

54
fl
chosen to be ‘%7' then § agrees with the first order
term of the formal asymptotic expansion of U6 given in

Chapter 3.

For the axisymmetrically loaded plate, (2-3.4)

can be solved analogously, by defining

5

(2-3.11)

where g3 satisfies

v :3 = 9

:3 O on B0
B; B( 2w )
3 V o

755' -—Tfif_— on B0

Here it is less clear, but still true that :1 and 52

satisfy the necessary boundary conditions :1 = g2 = O.
g1
g2

tangential and normal to the boundary, both of which must

Considered as a two dimensional vector, has components

vanish. By virtue of the boundary conditions assigned to
E3. the normal component of (2-3.11) vanishes:

g2 ' an V 0 an g3 ‘ °

As for the tangential component, since g3 vanishes

on the boundary, its gradient is perpendicular to the

boundary. And since ‘wb is the classical solution to the

55

axisymmetrically loaded circular plate problem, it is
a function of radius only, as is vzw , hence v(v2wb)

is also perpendicular to the boundary. Thus the tangential
components vanish for each term on the right side of

(2-3.ll). Since the vector has zero normal and tangential

components at the boundary,

The difficulty in applying this indirect method
for the general plate problem is inherent in (2-3.4).
Because :3 must vanish on the boundary, its gradient will
always be normal to the boundary. Because g1 = g2 = O,
the entire left side of (2-3.4) will be normal to B0.

whereas, in general v(v2

wb) 'will have both tangential
and normal components. Except for the two Special cases
above, (2-3.4) is generally inconsistent with the required
boundary conditions on g. It is of interest that (2-3.4)
is exactly the condition which prevents the evolution

of even the first order term of an asymptotic expansion
of Us in powers of 3. Such inconsistencies are
usually an indication of boundary layer phenomena. In
one dimension, boundary layers can be handled by matching
expansions valid near the boundary to those valid away
from the boundary. In two dimensions, even for simple

geometries this is not very hOpeful, since the boundary

layer thickness will tend to vary.

56

The presence of difficulties near the boundary is
not altogether unexpected. Even solutions of problem (C)
exhibit shear singularities at corners of the boundary,
when they are present. This is discussed in Chapter 4.
The boundary phenomena have implications which may affect

approximate solutions even well away from the boundary.

2-4. A Finite Element Theory for Problem (I)

Before using a finite element method to approximate
solutions of problem (I), it must be determined how the
quality of these approximations depends on the shear
rigidity, measured by e-l. This is illustrated by

error estimates in terms of e and the finite element

mesh size h.

In addition, the value of using finite element
approximations to solutions of problem (I), in order to
approximate solutions of problem (C), is investigated
using methods similar to Babuska and Aziz [3 ]. When
successful, this indirect method avoids the use of
complicated elements normally required to approximate

solutions of problem (C) directly.
Consider the weak form of problem (I). For each

_ T . . .
e > 0, U6 - [¢x,¢y,w] satisfies (2-l.9), i.e.,

B€(U€,V) = PL(F,V), for all V 6 He .

57
Let Sh be an N-dimensional subspace of He with
basis {¢1,¢2,...,¢N} as yet unspecified. Let

N
Uh = .2 qiq’i °
i=1

In order that Uh best approximate US (in the
.sense of the energy functional B€(-,°)) the qi are
chosen to satisfy the algebraic system of linear

equations

lpooo'NO

(2-4.1) B€(Uh,¢i) = PL(F,¢i) , i
that is

N o
(2‘402) i§1 qu€(¢j'¢i) = PL(FI¢i) o 1 = lpooo'N 0

Equations (2-4.2) can also be derived by minimizing

the functional
(2-4.3) J = B€(U,U)-2PL(F,U)

over all U e Sh' i.e. by chosing q1,q2,...,qN to

minimize
N N N
(2-4.4) J = B€(j2:31 qj¢j. 1:31 qi<Pi) -2PL(F, 1:31 qi¢i)
N N N
= jéi £E>1qjque(¢j'¢i)'-2 i2: quL(F'¢i) -
Setting ‘ggi = 0 produces (2—4.2), noting that Be is

a symmetric,positive definite bilinear form. The existence

58

of a solution Uh e Sh to (2-4.l) and its property of

minimizing J follow from these corollaries to the

theorems in 2-2:

Corollary to Theorem 3.

Let S be any closed subspace of He' There exists a

unique solution Us satisfying

B (U ,V) = P (F,V) for all V e S .
e S L

Corollary to Theorem 4.

US minimizes
J(U) = B€(U,U)-2PL(F,U) over S .

The proofs are identical to those of the theorems

themselves.

The following Theorem, with its elementary proof,

define the sense in which Uh e Sh is best approximation

to U .
3

Theorem 5: Let S be a subspace of He' If

Us 6 S satisfies
B€(USIV) = PL(FIV) a V V E S i

then

. 59
B€(U€-US,U€—US)gseme-vme-V) v ves .

 

m
E
I
<:
c
I
5
ll

B€(U€,U€)-ZB€(U€,V)4-BS(V,V)

B€(U€,U€)-2PL(F,V)4-B€(V,V)

The right hand side is minimized provided
-2PL(F,V)4-Be(V,V) is minimized. By the Corollary to

Theorem 4, Us minimizes J(V) = Be(V,V)-2PL(F,V).

The system in (2-4.2) is often written in matrix

form:

ll
'11!

(2-4.5) K0

where '
Kij = BS (¢j a¢i)

is refered to as the stiffness matrix. Note that K

is symmetric and positive definite.

T
Q = [q10q200000qN]

~

T
F = [PL(FI¢i)] I

the discrete load vector.

Similarly, (2-5.4) can be written

(2-4.6) J = QTKQ-zf‘o

60

and (2-4.5) derived from it by -§§L = O.
i
Equation (2-4.5) is used to determine the finite
element approximation to Ue' The construction of K

and is‘" will be described in Chapters 3 and 4.

The basic elements used here will be vector functions
_ T _ T
of the form ¢i — [mi,0,0] or ¢i -[O,¢1,O] or
¢i = [O,O,oi]T. The functions oi, used to construct
all three types, will be piecewise linear in the numerical
examples presented later.
Consider a (t,k)-system Sfi'k as described in
Chapter 1. Let S = St'kxst'kxst'k, and assume for
h h h h 3

each fixed h, that Sh c:H€. Equivalently, Sh c (33) .
By the approximation properties of Sfi'k, there exists

W’e Sh such that

p _ .
”Us --WHS g C h “Uqu , where p. - min(t-s,q-s) .
Taking s = l,
(2-4.7) HUa -wu1 g c hulerllq .

Now applying Theorem 5, with the coercive and continuity

properties of Be ,

61
“Us ”’11“: 5- C Be(Ue ' h'Ue ‘01:)
g c same --w,Ue -W)
s cuue mu:

_<. c e‘luue-wni

-1 2p 2
go 6: h ||U€||q .

That is,

(2-4.8) Hue -UhlIe g C 12"”2 hm”0qu
and therefore,

(2-4.9) nue - hnl _<. c .-1/2 w‘uueuq .

In the important case, t = 2 (piecewise linear elements)

-1 2
(2-4.10) ||U€ -Uh|l‘3 g C e / hl|U€H2
(2-4 11) ”U -U H g c e‘l/2 h||U 1|
' e h 1 e 2 '
The presence of e-l/z on the right side of these

estimates indicates that difficulty can be expected in
the accuracy of finite element approximations to solutions

of problem (I), when the shear rigidity is large.

62

Now consider the approximation of U0 by Uh’

Let U = ,wO]T solve problem (C), i.e.

o [‘wo,x"'wo,y

v4w =-2 on O

O D
3W0 .
wO=—S-N-=O on 80 .

Let Sh be defined as before. By the approximation

prOperties of Sfi'k, there exists W1 6 Sh such that

(2-4.12) HUG-WINS s C h“\onl|q .

where p = min(t-s,q-s) .
Using (2-2.8) and norm equivalence,

”U -W = “U -U +U -w
6 €

lHe 0 0 Inc

2 Nu, -vone+ IIUo-Wlue

g c(e1/2Hv(v2wo)|lo+ e-l/ZHUO “W1“1) '

Now setting 3 = l in (2-4.12), (2-4.l3) becomes

1/2 -1/2

(2.4.14) Hue -wl||€ g C(e ||v(v2wo)H0+ e hullU

Dug) ,

p. = min(t-l,q-l) .

Of primary interest here is the case t = 2 (piecewise

linear finite elements).

Then (2-4.l4) gives

63

1/2 -1/2

(2-4.15) ler-Wlllegme llv<v2wo)||o+e huv

ouz)
Since Hv(v2wo)“0 g CHUOHZ' (2-4.15) can be written more
concisely as

(2-4.16) MUe -W1H€ g cal/2+ 61/2 h)HUOI|2 .

If piecewise linear finite elements are abandoned
in favor of piecewise quadratics, the power on h in
(2-4.14) can be improved to 2, but only in the presence

of enough smoothness in U to make HUOH3 < e. The

0
question of smoothness for the solution of problem (C)
‘will be discussed later for the case of domains with
non—smooth boundaries. If the boundary is infinitely

differentiable, the result mentioned in 2-2 holds, namely
llwo||4 3 Cano -

In terms of the vector U0, and with f ='§'6 Ho

(2.4.17) ”no“, 2 cuguo .

The load p may be smoother, allowing potentially better
approximation, but the main limitation is the degree of

piecewise polynomials used in the finite elements.

A compromise can be reached between piecewise linear
and piecewise quadratic elements by using the quadratics

to generate the third component (displacement) and linear

64

elements for the first two. In this case a power h2

in (2-4.l4) is retained as if quadratics had been used for

all three components.

Using (2-4.16), which is based on piecewise linear
finite elements, a comparison;can be derived between Us

and U the finite element solution of (2-4.l). By

h!
Theorem 5, coercivity and continuity of Be' and (2-4.16) ,

Hue 41h“: S c sews - h,U€ -Uh)
g c sews -W1,U€ -w1)
g cuue -w1u§
g Mel/2+ 61/2 h)2||00||§ .
That is
(2-4.18) Hue -Uh||€ g cal/2+ 6-1/2 muuou 2

and by norm equivalence

l 2 -1 2
(2-4.19) HUe-Uhlllgme / +6 / mucouz .
Finally, by (2-4.18) and (2-2.8).

(2-4.20) HUG-Uh“, s HUG-USN,“ "Us “Uh“e
g.C(el/2+'€-1/2'h)”Uo“2 '

Similarly, by (2-4.19) and (2-2.9),

65

(24.21) HUG-.1111”1 g Mel/2+ 61/2 h)HUo||2 .

Estimates (2-4.20) and (2-4.21) show first that no
convergence is expected either for h fixed and 3

going to zero, or for 3 fixed with h going to zero.
Rather both must go to zero, and it is even important that
h go to zero faster than 31/2. Both terms on the

right can be made to converge at the same rate by choosing

e proportional to h. Then

1/2

(2-4.22) HUG-U g c h ”u

h“ 1 0H 2 -

While (2-4.22) is appealingly simple, and setting a
proportional to h seems a natural way of disposing of

the freedom in picking e, in practice, highly unreliable
numerical output will result. The mesh size h is, in
practice, finite. In fact, it is severly limited from
becoming too small, due to the rapidly increasing dimension
of the stiffness matrix. Since the coefficient of h]‘/2
in (2-4.22) may be very large, (2-4.22) guarantees little
if any accuracy. A method of extrapolation as well as

some inferences from the asymptotic and numerical results

will prove more reliable.

(2-4.21) also reinforces the expectations expressed
in [44,18] that letting 6 become extremely small "spoils"
the solution. This is seen through the negative power of

e multiplying h. It appears h must be made small

66

enough to overcome this adverse 6 factor, or a way

must be found to make inferences about U from finite

0
element solutions Uh generated by e values which are
not "too small". Since h is severly limited the second

approach must be followed.

It should be pointed out that the requirement

3
S c (Hg) is somewhat restrictive. In practice the

h
finite elements usually have their support over a poly-
gonal domain. Generally, nodes of the subdivision of the
original domain 0 are placed on the true boundary. If

0 ‘were a smooth non-convex region, being subdivided

by triangles, for example, there would be points outside

0 where a finite element trial solution would be non-zero.
Since functions in (H35)3 must satisfy homogeneous boundary
conditions, the containment Sh c.(H(]5)3 'would fail, as
would the estimates derived. In this case it is necessary
to carefully analyze the so called "skin" region near the
boundary. The containment Sh c (H35)3 can be guaranteed
with polygonal subdivisions of 0, provided 0 is

convex, or'a finite union of convex sets. Of course,

boundary regularity, as discussed earlier, is a concern

in addition to that mentioned here.

Estimates (2-4.18) - (2-4.21) have counterparts
showing slight improvement for the special cases of
section 2-3, namely the clamped beam and the clamped

circular plate with axisymmetric loading.

67

Recalling the method of 2-3, Ue = UO-e§+-n, where

g is a vector function depending only on U0. ‘With

the same definition of Sh as before, there exists ‘W1

and W: in Sh such that

(2-4.23) "U() —W

1H1 S C hu||Uo‘,'q )1 (t- 1.q- 1)

(2.4.24) Hé-qul g c 11““qu u = (t- l,q—1).

Now using (2-3.7), (2-4.23). (2-4.24)

+eW

(2.4.25) Hue—w 2H,

1

=|Wb-e§+n-Wr+ewflu

1.<. HUG-W1“, + en: -w2ne+ M16

-1/2

gc<e 1\vo-wllll+smug-5111+snail)

.<. «61/2 h“HUo||q+ 61/2 huugnq+ enéllp
Taking q = 2, t = 2.

(2-4.26) HUS-W +6 W

l Zue

-1/2 1/2

sue human,” huzuz+eugul>

Using Theorem 5, and coercivity and continuity of Be

and (2-4.26)

68

2
(2-4.27) \er -UhH€ g c B€(U€ -Uh,U€ -Uh)

g c same-w +3 W2,U€-W +e W2)

1 1

2
1+ 3 que

g cuUe -w

s curl/2 hnuonz+ 61/2 huznz+ auguy2 .

(2-4.28) Hue — Uh“e

-1/2 1/2

gen: hnuouzn huguznngul)

and (2-4.28) with norm equivalence gives

(2.4.29) Hue..uhu1

g cw” hnuouz+ .1/2 hngnz+ engug .
N...
(2-4.30) lon-Uhll€ _<. HUG-USN; \er -UhH€

combined with (2-4.28) and (2-3.8) gives

(2-4.31) HU0--Uh||e

-1/2

g cm hnuo\|2+ 61/2 hlléH2+ eugul+ el/Zugnl) .

However, using (2-4.29) and (2-3.9),

69
(24.32) HuO - hnl

g cue-V2 hHUOHZ+ e1” hugnz+ eH§l|1> .

(2-4.32) is improved over (2-4.21) by having 6 rather
than el/Z in the term independent of h. The additional
term 31/2 huguz is of higher order than the term before

it, for small 6 and h.

CHAPTER 3 - BENDING OF A CLAMPED
TIMOSHENKO BEAM

The clamped Timoshenko beam is the one-dimensional
analog of problem (I). Its formulation can be derived
from that of problem (I) by deleting all dependence on
y. Similarly, the clamped beam in the classical theory
can be described. Both reduce to two point boundary
value problems which can be compared and analyzed in
more detail than the two dimensional boundary value

problems posed by problems (C) and (I).

An asymptotic analysis of the two beam problems
is carried out in 3-1. 3-2 contains an asymptotic
pointwise analysis of discretization error generated by
approximating the Timoshenko beam solution by piecewise
linear finite elements. In 3-3 the element stiffness

matrix is constructed, and 3-4 contains numerical results.

3-1. Asymptotic Expansion for the Solution to the

Clamped Beam Problem in Timoshenko's Theory

Define the Operator L by

__d_2_ o 1 1

dx2 -1 d"

(3-1.1) L = + e 2
o o __d_ __<_1_

dx dxz

71

and

(3‘102) U = p F = o

The differential system

(3-l.3a) LU = F

represents the Timoshenko beam equations when f1 = 0

and f2 = g. Consider also the clamped boundary condition
0 1

(3-1.3b) U = at X = i=3 .
0

In what follows, f1 and f2 ‘will both be allowed to
be non-zero, and will be considered infinitely differentiable.
This more general right hand side will need to be

considered in the discussion of principal error in 3—2.

Due to the homogeneous boundary conditions, a

solution U to the system LU = :: can be Written as
(3-1.4) U = U(1)+U(2) '
where
(3-1.5) LU(1) = f1.
0
(3-1.5) mm = f0

72

and
O
(3-1.7) U = UH") = U(2) = at x = t-;- .
0
Let
a B
(3-1.8) u”) = . U”) = .
u w

Expanding each in powers of the small parameter c. and
solving for the coefficients produces a uniform asymptotic

(1) (2)

expansion for each solution U and U .

First, let

(3-1.9) 0(1) = 2 e1 0;“

that is,

(3-1.10a) a = Z 6:1 a. , a. = O at x = *% '
i=0 1 1

(3-1-10b) u = 2 61 u. , u = O at x = *%
i=0 1 1

Substituting equations (3-1.10) into (3-l.5) and equating

powers of 6 yields

I _

(3-1.lla) cod-no - O
I II _

(3-l.llb) cod-uo - O

I II _
l _ v
I I _
(3-1.13b) az-I-u2 - 0

etc., in general

I _ N
(3—1.14a) civl-ui - ci_1

I a _ '
(3-1.14b) .cxi-i-ui — 0 for

in: 2,394,... o

Differentiating (3-l.11b), substituting into (3-1.12a),

differentiating (3-1.12a) and subtracting (3-1.12b) yields

(3-1.15a) u" = f

Since u = O at x = :h%' (3-1.lla) implies

O

(3-1.15b) u’ = -a = O at x

and u is uniquely determined. (3-l.11a) then gives

0

(3-1.15c) a = -u’

In a similar manner, 111 is seen to satisfy

(3-1.15a) u”" = O

74

with

_ I _ m _ l
(3-1.15b) u1 - O and u1 - f1--uo at x - ‘12 .
Then

_ - III _ I
And u2 satisfies
(3-l.l7a) u"” = 0

(3-1.17b) u =0 and u’ = -u”’ at x= :% .

2 2 1
Then
__ __ III__ I
(3-l.l7c) a2 - ul u2 .

In general, for i = 2.3.4,...

(3-1.18a) u’.'"= O

_ I _ _ III = _1_
(3-1.18b) ui — O and u. — u. 1 at x :hz .
Then

__ III _ I
(3-1.18c) oi — ui-l ui .

(2)

In exactly the same manner U is determined.

Let

75

(3-1.19) U(2) = 2') e1 1152’
. 1.
i=0
0
where U12) = at x= ti o
1. 2
O
that is,
(3-1.20a) 5= 2 31:31. 5. =0 at x= :%
i=0 1
(3-l.20b) w=chwI w.=0 at x=¢-;'-.
i=0 1 1

Substituting equations (3-1.20) into (3-1.6) and equating
powers of e as before produces the following conditions

which determine the qi and wi:

(3-1.21a) wig: f2

(3-l.21b) wO = wc’) = O at x = 1%

(3-l.21c) BO = -w6

(3-l.22a) wi”: —f£

(3-1.2213) wl = o, wi = -wg at x = 1%
(3—1.22C) Bl = dwg-rwi

(3-l.23a) WE”: O

(3-1.23b) w2 = o, wé = —fé-wi’ at x = I;
(3-1.23c) 52 = -fé-w’i’-wé

76

(3-l.24a) WE”: O

(33-1-2413) W3 = O, wé = -w’é’ at x = ¢%
(3-1.24c) 53 = .wléuwé

and for i = 3.4.5,...

(3-1.25a) WE”: O

(3-1.25b) wi = 0. w; = w§:_1 at x— 1%-
(3-l.25c) Bi = aw:_l-wf

Some useful observations can be made from these expansions.
If f2 is an even function (f2(-x) = f2(x)). then the
series in (3-1.lO) truncate after only two terms for the
following reasons: First. if f2 is even. then each

wi is even. Next, the function defined by

II!
0

__ I _ m - - ‘ I = _ II_ III
A(x) - -f2 wl is constant, Since A (x) f2 w1

by (3-l.22a). However by (3-l.23b) wé(%) = 13%) =

I

2 2
that -wé %-) = wé( -%) . So they both vanish. Thus from

A( -%) = wé( --]2-'). Since w is even, w is odd, so

(3-l.23a) and (3-1.23b), w satisfies

2
(3-1.26a) “1;”: o
(3-l.26b) w2 = wé = O at x = gr;-
implying that
(3-1.27) w a O .

77

The argument above also implies that A(x) vanishes,

and from (3-1.23c) it follows that

It then follows by the recursive nature of relations

(3-1.25) that

(3-1.29) Bi E‘Wl a O for 1 = 3.4....
To summarize, if :62 is an even function, then the
. (1)_ ‘1
solution U - to the system
‘w
° 1
LU: I B=W=O at x=¢—
f2 2
is given simply by
(3-1.30a) w = wO-I- e: wl
__I_ III_ I
(3-1.30b) B - ‘wb e(wo ‘w1)

'where w' and ‘w are determined in (3-l.21) and

O l
(3-l.22).
In particular. if f2 = g . p being the transverse
load applied along the beam, then 5 represents the
‘w
solution to the clamped beam problem in Timoshenko's

W

theory, and wb

beam problem in the classical theory. Similarly each "i

represents the solution to the clamped

78

represents a solution to a particular clamped beam

NIH

problem in the classical theory with slope at x = :h
determined by the previous term. New if p is a
symmetric load, then ‘w and B differ from ‘wo and awé,
respectively, by a term simply linear in e.

No such conclusion can be reached for p being an
antisymmetric loading. However, it is useful to notice

that since any load p can be written as a symmetric

plus an antisymmetric part, i.e.

(3-l.31) p = ps+pA

and under the homogeneous boundary conditions the

3
solution can be likewise decomposed
‘w
B B B
(3-1.32) = s + A
‘w ws “3
where
r
s O
(3-l.33a) L =
Ws ps _.
fiA 0 1
(3-l.33b) L = .
WA PA I

 

It can then be seen that all nonlinear dependence on e
is due solely to the antisymmetric part of the load.

Moreover, if the quantity of interest is the displacement

79

at the beams midpoint x = 0, then 'w(O) = ws(0)4va(0)
is simply linear in 3 since “3(0) is linear, and
wh(0) vanishes.

A similar truncation takes place for where

fl is an odd function.

A second observation will prove useful in dealing
with principle error in 3-2. If

(3-1.34) £2 = -£' ,

then by comparing (3-l.21) and (3-l.15), it can be

seen that

(3-l.35) w = —u

0 0 and BO = -a

o O

This means that the leading term of U = U(l)4-U(2)

vanishes. Thus the solution to the system
f1

I
-fl

(3-1.36) LU

is of first order in the small parameter c, i.e.

U = 0(6) .

In 3-2, it will be observed that terms representing
error from discretizing the problem by finite elements
‘will solve systems with the special form of (3-l.36).
This observation will lead to great improvement to the

form of that error.

80

To summarize, the solution to the clamped beam
problem in Timoshenko's theory can be uniformly expanded
as powers of e. The leading term represents the solution
(slape and displacement) for the clamped plate problem

in the classical theory. The difference

I -w’
- O is of order e. in fact is prOportional
‘w ‘wo

to e for symmetric loading. Other results useful in

controlling discretization error are also uncovered here.

3-2. Principal Error for the Clamped Beam in

Timoshenko's Theory

In order to examine the error due to discretizing
the boundary value problem by use of finite elements,

that is, by replacing

ll
H-
Mud

(3-2.1) LU=F, U=0 at x

‘with a system of linear algebraic equations,

(3-2.2) KQ = f'

it is useful to keep in mind that equations (3-2.2) are
simply a type of finite difference equations which are

selected to satisfy a variational principle rather than
being chosen in the more conventional way by replacing

derivatives in the operator L by certain standard

difference quotients.

_ 81
Viewing (3-2.2) as a system of finite difference
equations, an analysis of the local discretization error
can be carried out by determining the error as a power
series of the discretization parameter, in the case of

finite elements, the mesh size h.

To compare solutions to the discrete problem and
the continuous one at a node xi of the finite element
grid. it is necessary to examine two consecutive rows of
the stiffness matrix, namely. the (2i-1)th and 2ith,
where the unknowns qu-l and q2i represents the
approximations to ¢(xi) and ‘w(xi) respectively. This
.can be done either by superimposing element stiffness
matrices constructed in 3-3, as they would be assembled
in the g1dba1 stiffness matrix, or by computing the

elements Kij directly from the energy functional, i.e.

The former is simpler here, with rows given in (3-3.20).
NOte that it is necessary to divide equation (3-2.2) by

h in order to observe a "difference" form.

After dividing by h, designate these two

rows of equation (3—2.2) by

(3-2.4) LhUh = Fh

where

 

r l e-1
(3-2.5) “fi‘1 6 1
-1
_£L_.
2h
-1
(-—22-+ 263 )
h h
L =
o
-1
(-15+-€ )
h
-1
6:
2h
h 2
(3-2.7) f = S'f

82

-1/2

The difference form of Lh

h

L by the vector function

 

 

q21-3

q21-2

q2i-l

qfifil

q21-4-2

is made apparent by multiplying

 

83

w(xi-h)

Mxi)

C}!
II

(3-2.8)
w(xi)

¢(xi+ h)

 

 

w(xi+ h) .4

b

which corresponds to Uh. In this way Lh can be
treated as a difference Operator applied to vector

functions U = [¢,w]T.

For convenience, define the following difference

Operators at a node xi:

g(xi-i-h) + 4g(xi) + g(xi -h)

 

 

0°th = 6
g(x.+h)-g(x.-h)
Dllg] = 1 2h 1
2 9(Xi+h) - 29(xi) + g(xi-h)
D [g] = 2 .

 

h

These can be seen to arise naturally by multiplying

Lhfi and rearranging

h 432w] + e‘1<n°[¢1 +D1Iw1)
(3-2.9) L f:

-e‘1(D1w] +Dzlw1)

84

To simplify notation, let Lh represent both an

Operator on a vector function U = , and the matrix
‘w
¢
defined in (3-2.7). That is, for U = define
w ,

the operator Lh by

(3-2.10) LhU = L U

where LhU is the matrix product given eXplicitly in

(3-2.9). Thus (3-2.10) written explicitly is

-nzm + e'1(n°m +D1tw1)

(3—2.11) LhU
-e'1(D1[¢1+D2[w])

Since

-¢” + e-IW-i-W')
(3-2.13) LU =
-6-1(¢l+wfl)

the similarity between Lh and L is Obvious.

Note that for all local analysis which follows,
function evaluation is understood to take place at node Jot

unless otherwise indicated.

The variational method, in generating the finite
difference equations by use of piecewise linear finite

elements, "prefers" tO replace all the derivatives in

85
L by central differences. but elects to replace the
function evaluation of W by a "Simpson's average",

((x+h)4-4¢(x)4-w(x-h)
6

Error due to replacing a continuous system by a
discrete one is examined by substituting the solution
to the continuous problem into the discrete Operator.
and comparing to the result of applying the discrete
Operator to the discrete solution. That is, if LU = F

is discretized by
then the quantity

Lh(Uh-U) = Fh-LhU

is a measure of consistency between the two systems.

Using Taylor series about a grid point, (Fh-LhU) is
expressed as a series in h (the mesh size). The leading
term being, say 0(h), indicates consistency. Further
refinement can be done by defining this leading term to

be hF

1, where Fl independent of h, and defining

e1 to satisfy the continuous problem

(with apprOpriate boundary conditions) then computing

86
Lh(Uh-U-he1) = Fh--LhU--the1 which in turn is

expected to be, say 0(h2), and so on. In this way

an error expansion

=Uh- Z} e. lhi

i=1
is obtained.

The following theorem gives the computation needed

to produce each Lhei, given Lei' The notation f(B)

indicates the fith_derivative of f.

f1 0
Theorem 6. If LU = F = and U = at x==
f2 0

MHA

then for u, f1 and f2

following asymptotic form holds:

infinitely differentiable. the

 

(3-2.20)
f
LhU = 1
f2
.1
(2c) (20-l)
h2a 2 f1 -2a f2
+ :‘11(2—a+2'TT (2(1)
2 f2 J
r- q
2a -—(a 2-1)<f‘2°‘12)+ $213))
+ e-l' Z)-—lL-T-
_ (2a+2 :
0-1 20 (f(2a-1)+f(2a-2)) J

 

 

Proof:

(3-2.21)

then

(3-2.22)

87

Using Taylor series eXpansions

 

 

 

 

 

c (a)
1Hx+h)= h‘z.(x)
c=0 '
a (a)
u(x-h) = h ua' (x)
c=O '
2a
0 1 h u
D [‘1]: “*3 “El (2a):
20 (2c+l)
l u
D [u] = “2:30 (2a+1)l
(20-2) (2a)
2 _ h u
D [u] - 2 c=1 (2c)!

Using (3-2.22) in (3-2.11). taking U = [v.wIT.

(-

LhU=

 

rearranging

-2

 

 

 

 

 

 

h(2c-2)¢(2c)-
(1:1 (20):
O 4
F 1 ‘ h2a¢(2c) h2qw(2c+l)
1 “-30:51 (2“): +<::=o (n+1):
e-
h2cw(2d+l) h(2c-2)w_(2c)
a=o (2a+1): + 2 0:1 (2a):

 

 

J

88

 

 

 

(3-2.23)
_¢II+ e-l(¢+wl)
'LhU =
—1 , u
-e (w 4“W )
2c _2W(2a+2)
+ 23 h
_ (20+2)!
a-l 0
2c FL20+2)3(26+]J ¢(2c) + (2a+2)w(20+1)
+ 3'1 23—13-—
a=1(2c+2):
(2Q+2)¢(2c+1) + 2W(2c1+2) J
1) £1
Now since L = , the first term on the right
‘w f2
f'1
of (3-2.23) is . Moreover, by differentiating and
f
2

combining the equations in LU = F the following identities

are found:

(3-2.24a) w'” = -f]’_-f2
(3-2.24b) w’” = -efg+ fi+ f2
(3-2.24c) (+w’ = e(¢”+£l)
(3—2.24d) ¢’+w” = -ef2

Repeated differentiation of these give

89

(3-2.25a) W“ = 4:54) ”593’” . for a 2 3
(3-2.25b) w”) = -€ 1393"” +fia-3)+féB-4). for a 2 4
(3-2.25c) ((BHw‘B”)

= -3 féB-l), for B 2 l

Inserting these into (3-2.23). then

(3-2.26)

 

 

 

 

 

 

 

 

 

 

 

 

r- ‘1 F' .1
III (204-2)
h f1 h2 21 r12“ ‘21
L U = + --,- + Z
_ f 4. _2 W
2 o ‘1" O
L .4 b d
_1 2 F4¢'+4w”’
+ e P
4' ‘IM
-4 ”-2w
L. W .J
(2a+2¥f2a+l) ¢(2a)+-(2d+2)w(2a+1)
-1 h“
+ —
‘1 G232 (2a+2):
L _(2a+2)¢(20+1).-2w(20+2) J
f‘ II I1
f1 hz 2“31‘212
= f + E I
2 L_ 2f2 a
r. H
hm 2:9“) - 2cfé2a-l)
+ Z ———.
L 2 J

 

 

 

 

 

O
+ 6-1h2
4! I
2f1+2f2
r 1
(2c+2)(2-2c) (2c-2) (2c-3)
3 (f1 +-f2 )
_1 hZG
+ e 0232 (2a+2):
' (2a-1) (20-2)
L 2c(f1 4' f2 ) .4 .
That is,
f 2 2f?” - 2cfé2a'1)
h 1 h ‘1
L U = + Z W
f c=l ° (2c)
2 2f2

P 2 u
2 __4(c3-l)(f§2a-2)+_fé2c-3))
O

+ -1 23 h
6 (1:1 __T(2a+2 :

 

 

L 20(f£2c-1)+_fé2c-2)) J

which is equation (3-2.20).

To compute the principal error, the expression,

LhUh must be compared to LhU. LhUh = Fh can be

expressed as a Taylor series by expanding the integrals

defining components of Pb.

The two components Of the load term on the right

side of (3-2.4) are (see (3-3.3) for definition of cpi)

h _
h 1 1/2 1 xi+h
(3-2.27b) F21 - 3111/2 f 1 dx - 31‘ f cpl dx

91

Define

X
(3-2.28) f('1)(x) = j f(s)ds
x.
l

X

(3-2.29) f“2)(x) = j f(-l)(s)ds .
X.

1

Integrating by parts

(3-2.30)

xi+h
I h f(x)¢i(x)dx

i
xi xi+h
= f(x)¢.(x)dx4-
I . h l J‘x.
i i

f(x)¢i(x)dx

x. x.+h _
f1-1)(X)cp}_(x)dx-J‘ 1 f('l)(x)cp:]-L'(X)dx
X.

i i
x.+h

= -lj‘ 1 f(-l)(x)dx+% f 1 f(-1)(x)dx

XO-h x.
J. .1

£“” (xi+h) - 2f(‘2) (xi) + f“” (xi-h)
‘ h

 

f(-2)

EXpanding as a Taylor series

- j
(3-2.31) f(-2)(xith) = Z £(3‘2)(xi) (—1.h.—)—

j=o 3'

Then (3-2.30) becomes

92

xi+h1
(3-2.32) f f(x)¢i(x)dx

x.-h
i

i=2
j even

2c-l
_ (20-2) h
— 2 2 f (xi) W

c=l

2a+l
_ (2a) b
‘ 2 ago 1 (xi) (2c+2)!

20+1
hf<xi)+2 23 f12a)(xi)(—hm

c=l

and (3-2.27b) becomes

(3-2.33) F

h 1
21 if

h2c

W '

f(xi) + 2 231 f‘za) (xi)
(3:

Inserting these expressions for Fh into (3-2.34)

 

 

 

)- O 1
(3-2.34) LhUh =
h20f(20)
f1’2 (2c+2)'
L c=l '_ -
. l 2 . .
By Theorem 6, if U = 6 (HO) satisfies
'w

0
LU

ll
'11
II

, then
f

(3-2.35)

93

 

 

 

 

O
f
_ 2 1
_26f(2o-l)._4(c -l) 6-1 f(201-3)
h2a 3
+ £1 (2_a+' '2"):
2f(2a)4-2a e-I f(2a-2) L
0 W
2a
h (2c)
“2 Z ‘('2'a_+2'): f
L a=1 _
- 2 T
_2af(2c-l)-_4(c -l) 6-1 f(2c-3)
2c 3
+ Z) h
afl_QmQH
" 2 -1 f_(2<::—2) J
L a e

 

Subtracting (3-2.35) from (3-2.34)

(3-2.36)

 

Lh(Uh-U)
2a
b
= -031 (2a+_2)':
L 2a 6

 

v— 2 1
_2af(2c-l)-_4(c -l) e-l f(2c-3)

-l f(2O-2)

The first thing to be noticed in (3-2.36) is that

Lh(Uh-U) = 0(h2) .

 

94

meaning that the discrete and continuous problems are

consistent, with the error being quadratic in h.

The second is that the factor 6-1 is present in
the leading and higher order terms. Since e is a
small parameter, this is undesirable. However, by

considering an "error" expression
-1 2
(3-2.37) e= (1+E—f21‘—)Uh—U

in place of Uh-U, considerable improvement is attained.

 

 

 

 

 

 

 

 

 

 

 

 

-1 2
1.11((1+€12h )Uh-U)
P o 1 — o 1
—1 2
= (l-I-elzh ) _
h2cf(2c) h2cf(2a)
f+2 2 . f+2 ,
L c=l ( c+2).d L “=1 (20+2).J
F' 2 fl
_2af(2a-l) _4(a3-1) e:-1 fag-3)
2a
_ 23 _Tl__
d=l (2c+2)!
2a 3.1 f(2a-2)
L .—
O
-1 2 O -1 2c
=€12h + 6 h2 '
f c=2 ( o).

95

 

 

 

-2cf(2a'1)— 4(c2- H) -lf(2a-3)
2a 3
- z; __D___
c=l (2c+2):
2a 6-1 f120-2)
,
Mf<2a-1)++3(02—1)e-1f(2a-3)
2 f 2c
=h— + 2 _ll___
2 O c-2 (2c+2)!
- (a-l)(2a-l) -l f(2c-2)
3 6 J
L
The leading term is now free of 3.1. This gan now be
h e
used to generate the principal error term 121. in e,

by letting e1 satisfy the equation
f’ o 1
(3-2.39) Lel = , el = at x = i3
0 O

The quadratic term of the error e having
the expression

2

-l 2 h e
e h Uh 1
+ 12 ) ’U"

(3-2.40) (1 12

 

0(h4). To recover the 0(h4)

 

been recovered,

is expected to be term,
the expression

M26
(3-2.41) Lh ((1+— —1fi—)uh U-—— 121
is considered. Note first that, from (3-2.39), e1
has e-dependence through positive powers only. By

equation (3-2.20)

96

r- 2 .1
2f(2a+1)_4(a3-1) e-lf(20-l)
2

h C1
+ Z)-—————yr

 

 

2a 6-1 f(2a) J

r 2 ’
2f(2a+1)_ 4((1 3—1) e-1f(2c-l)

 

2c e‘1 f(2a) J

 

(n+1) (20-1-1) f(2a-1)
2 f’ h2c 3
(2O+2):
O

P

-_20(a-2)(d+1)(20+1) f(2c-3;
9

(2o+2 :
“=2 (c-l)(a+l)(2c+1) f(2a-2)
3

h

 

 

Now (3-2.43) and (3-2.38) give

97

 

 

(3-2.44)
-1 2 2
h ___e h 1:.
L “1+ 12 )iV-U-lzel)
a r- '1
F-f”’+4e'1 f’ _(a-1)3(2a-1) f(zen-1)
='—T + Z) -———-2--5-T
6- -1 0:3 (201+ .
-4 6 fl! 0 d
c J b .

 

 

 

2(a+l)(2a3-3a2

h2c 9
+ 9-1 :8
a=3 725127?

1
+4a-6) f(2c-3)

__Lg:l)(2a2+afgl f(2c—2)
L 3 d

 

 

Let e2 be the solution Of

 

+fm-4e-1f' O
(3-2.45) Le = , e = at x= *1:
+4 6 f”
then the expression
e-1h2 Uh hzel h4
(3-2.46) (1+T) -U- 12 + 6—: e2

is expected to be 0(h6). Consider the e dependence

in e2. Ordinarily the presence of 6-1 on the right
side of (3-2.45) would indicate that e2 is 0(e-l).
-1 fl
However, 3 is multiplied by a vector f where
2
f = -f’, a condition which causes the adverse e

2

98

dependence to vanish in the solution e This follows

2.
from the asymptotic analysis Of the system
f
LU = 1 from 3—1. Hence e2 has dependence only on
f
2

positive powers of e, and again the principal error

term remains free of adverse e dependence.

This special relationship present in (3-2.45) in
the e-l-term provides another advantage, this time in the

determination of e3. To compute Le3, the expression
4

gT-Lhe2 must be added to equation (3-2.44). Having
e71 in Le would be expected to lead to 6-2

2

dependence in Lhez, but once again the form Of the 6-1

dependence in Le2 and in (3-2.20) combine to leave only
fill

 

e-l dependence in Lhe2 (coming from the term in
0
Le2).
4 h III _ I
h L e2 - h4 f + 48 1h4 —f
"‘ET'” “'5? sf I
O f
2c 2f(2a-1)-%-c-1(a+3)(a-1)f(2
h
+ 2 (2a 2Yl6l
“=3 2(d+2)e'1 f(2a-2)

with (3-2.44)

o-3)

 

99

 

 

 

 

 

 

 

 

 

 

 

(3-2.48)
—l 2 2 4
Lh((1+2.fill—)uh-u-§.2. e1+2—: e2)
1
( (c-l)(2c-l) + 2 )f(2c-l)
_ 2(2a+2)? (2d-2):61
= h2c
c=3
L 0 J
_2(d+ll(2G3-3a2+4G-6) f(20-3;
2c 9
-1 h
+ 6 Ag; (2a+2): 2
_ (c-l)(20 +d+2) f(2O—2)
L- 3 .4
F_4(d+3)(d-1) fad—3f
2d 3
+ '1 Z) h
e a_3 (2d—2):6:
’ 2 (2d-2)
J. (n+2) f J .

 

The leading term which will be used to generate e3 is

 

 

P u
g_ (s) -l (3)
6 3 f -8 e f
(3-2.49) ‘5?
L _
again exhibiting the Special form in its 6-1 term.
This special prOperty of the 6-1 term can be

shown to propagate over and over as the error terms are
uncovered. Although tedious, an expression can be

develoPed for determining the general principal error

100

function ek base on extending patterns of coefficients

which hold for the first several terms.

Define ek e (Hé)2 by

 

 

f(2k-l) f(2k-3)
-l
(3-2.50) Le = — 2k 6
k 0 _f(2k-2)
for k = 2.3.4,... and
f!
O
The lending error term is
20 thLe
(3-2 52) k k
° (2k+2):
Ck satisfies
2 k-l 2k+4

j=l 2j+2 3

the first few terms of Oj are:

01 = 1

G2 = ‘“%

03 ='%

C4 = ‘“%

05 = 5

o = .végl .

101

Expressing the principal error itself,

 

 

 

 

 

 

 

23'
-l 2 m 20.h e.
.§__}.‘_ h - 1 J 0 2m+2
(3-2.54) U-(l+- 12 )U — jg; (2j+2): + 0(e h )
where ej is of order C(60) for each j.
For m.= 3,
- 2
(3-2.55) U-(1+-“-lTlh—)uh
hze h4e 4h6e
= 1 - ———2-+ -———§-+ 0(e0h8)
12 6: 3-8! °
-lh2
To see the significance of the factor (1+612 ), rewrite
(3-2.55) as
(3-2.56) 1 —U—Uh
e-lh2
(1+ 12 )
hze h4e 4h6e
= 1 <——1--—2-+—3+o<e°h8))
-l 2 12 6: 3-8:
(1+_€__£.)
12
1 e-lh2
expanding _1 2 for small values of -jEf—-,
e h
1+» 12
-1 2 i
(3-2.57) U( 23 (-m) )-U‘h
. 12
i=0
- 2 4 6
-1 2 1 h e h e 4h e
= _€ 11 l _ 2 3 O 8
(Z ( 12) )(12 6: +3.8: +O(€h))

i=0

102

2
h _ -1 h
( ’2U+ '1e +1 e )-1-1-4—+ 0( '3h6)
’ 3 e 1 s 2 12: 3 °

(3-2.58) is the usual error eXpression which would have

- -l 2
been Obtained without involving the factor (1+HEIEE-).

As expected it is "contaminated" by negative powers of

e, in fact to larger and larger negative powers in the
higher order terms in h. (3-2.58) even seems to violate
the theoretical estimates in Chapter 2, which require
U--Uh to be at worst 0(6-1/2h) (although these are
integral estimates and (3-2.58) is pointwise at a node,

the e-dependence shou%d be comparable). However, using
-1
e h
1 . _ I .

e h e h
]u+-—I§—' Ink—i§-—

e-lh2
12
-1'2
6 h,
1+H-jjr—

(3-2.59) U- = U

2 4 6
l h e1 h e2 4h e3

0 a
““':I'2 ( 12 " 6: *‘4L-8:"°(€ h 1"

e h
Ink-jE§—-

In this closed form the e dependence is clearly C(60)
for small 9. However, the leading error term on the
right of (3-2.59), though technically 0(h2) for fixed

6, becomes nearly independent of h for very small 6.

103

This appears numerically by values of Uh being driven
to zero as 3 becomes small.

-1 2

It appears that the factor (lure 12 ) "collects"

 

the error terms of all orders which contain the adverse

_ -1 2
dependence on e 1. allowing the quantity (1+___.€12h

)Uh

to approximate U much more reliably than would Uh
itself for small e. Numerical results confirm this

effectiveness.

3-3. Construction Of the Element Stiffness Matrix
for the Clamped Timoshenko Beam

Consider a beam clamped at x = :h%- and subdivided

by nodes Xi= ih--% in n = 3 equal segments.

As in 2-4, the discrete form of the improved plate

equations can be derived by minimizing
(3-3.1) J(U) = B€(U,U) -2PL(F,U)

over a finite dimensional space generated by a basis

of finite element functions. Taking
(3-3.2) U= Z qi (bi .

where, for i = l,3,5,...,

Define

(x-xi_1)/h for xi-l < x < xi

(3-303) mi

(x.

1+1-x)/h for xi < x < x

i+1

”i are the sO-called "roof” functions of Figure l.

J
J

1—1 xi 1+1

Figure l. Piecewise linear "roof" function.

With (3-1.2), (3-1.l) becomes

2n-2 2n-2
(3-3.4) J(U) = B€( j§1 qi¢i' j§1 qjcpj)
2n—2
- 29L(F. j§1 qj¢j)
2n-2 2n-2
= £31 jE], qiqj B€(¢i.<|>j)
2n-2

- 2 $21 quL(F,¢j) .

105

Or, in matrix form,

(3-3.5) J(U) = QTKo-ZQT§

where

Kij B€(¢i.¢j)

Q = [ql""'q2n-2]

W2
ll

j PL(F.¢j)

J(U) is minimized over the finite element space by

choosing Q to satisfy
(3-3.6) KQ = fi .

While K and fi are given explicitly in terms
of the base functions ¢i' it is usually more convenient
to construct them indirectly. This is done by constructing
energy expression J(U) as a sum of the energies
computed over the individual elements of the domain. The
fact that U is polynomial when restricted to a single
element can be exploited for computing the integrals

‘which make up K and E.

J(U) =‘Z)J(U(e))

I

the sum over all elements of

106

J(U(e)) = B€(U(e).U(e)) -2PL<F.U(e))

(e)

where U = U restricted to’a particular element.

(e)

First, each U is expressed in polynomial form for

the purpose of computing J(U(e)). Then U‘e) is
expressed in terms of the finite element base functions
having their support over that element of the domain.

By transforming the polynomial coefficients into the

coefficients of the base functions,

(3-3.7) J(U‘e)) = Q(e)TK(e)Q(e) _ 20(e)T§(e)

where 0(a) represents those components of Q which
multiply the base functions having support over that

element.

Once these elemental energies are computed they
need only be assembled into the total energy expression
J(U). The stiffness matrix K is the resulting

assembly of the K(e). and § is formed from g(e).

Consider the finite element space consisting of

piecewise linear elements. To compute K‘e) and g(e)

over a specific element xi S-x-S~Xi+h' let
¢(e) a1x+a3
(3-3.8) u‘e) = ) = .
(e
w a2x+a4d

 

Also

107

(e)

Fq

(3—3.9) (e)

 

The relation between coefficients is found as follows.

A

At x = Xi' c"i = 1' Qp1+1 =
alxi-i-a3
a2xi+a4

At x = xi+h, (pi = O, @14-

a1(xi+h)+a3
32(Xi+h)+a4
thus
r
x. O l
1
O x. O
i
xi+h C) 1
(D x.+h O
i

 

Inverted, this is

1 “’1 + q3 q"1+1

_qée)w

(e) 1

(e)

 

i + ‘34 ¢i+ld

 

 

 

 

 

 

0, so
F" '1
qie)
(e)
qu
1 = 1, so
‘fi Fqge)
- (e)
g
.J L' 4
1 r- 1 '-
0 a1
1 a2 =
0 a3
1 a4
d L .4 "

 

 

 

q

 

108

 

 

 

 

 

 

 

 

(3-3.1o)
F _, r' ’ (er
a1 -1 O l O W ql
a O -l O l g(e)
2 1 2
= h (a)
a3 x +h 0 -xi 0 q3
(e)
-34.} b O xi+h 0 -xi.. :14 J
or
(3-3.11) A = PQ(e) .
where
(3-3.12) ' 1
—1 O 1 O
0 -l O 1
x +h O -x. O
J.
L O x +h O -xiJ
F’ “ (er
a1 ql
(e)
a2 (e) _| q2
A: ' Q _ ( ) '
e
a3 q3
(e)
a q‘
L. 4.. 4 .4

 

 

 

Next, J(U(e)) is computed from the polynomial

109

(3-3.13) J(U(e))

= pBw‘e) 41(6)) + e'lpSm‘e) 41(3)) - 29L(F.U‘e))

x.+h 2 x.+h

_ J. -l J. 2
—j‘ aldx+e I (a1x+a3+a2) dx
x. . x.
1 i
x.+h
1 .2
-2 Ix D (a2x+ay)dx
i
‘ 3211+ -lh{(x2+hx -o-l"—2-)a2+a2+a2
' 1 e i i 3 1 2 3

+ (3132+ a1a3) (2xi+h) ‘l" 23233}

X-+h x.+h
- 1 ;E _ 1 ‘2
2a2 J. Dxdx 234 j‘ D dx .
' X.
1 J.
In matrix form. (3-3.13) is
(3-3-14) J(U(e)) = ATN(e)A-2A1E(e)
where
(3-3.15) N(e) =
P1+ ‘1(x2+hx +124 ‘1(x +11) 63-1”: +21) 0‘
e i i 3 e i 2 i 2
e'1(x +321) 6-1 6-1 0
h
“1 h. -1 -1
_ 0 0 O 0.1

 

 

110

(e) xi

wl
u

(3-3.16)

 

 

Now (3-3.ll) can be used to eliminate A from

(3-3.14).

(3_3.17) J(U(e)) = Q(e)TPTN(e)PQ(e)._20(e)TPT§(e) .
By comparing (3-3.17) to (3-3.7), one has

(3-3.18) K(e) = PTN(e)P

(3-3.19) E‘e) = 19?“) .

In more complicated settings, the matrix multi-
plication in (3-3.18) and (3-3.19) are left to be done
by computer (by choosing local coordinates, P and
N(e) each can be expressed independent of the element
location). This is done in Chapter 4 for the clamped

(e)

plate problem. Here K can be computed directly.

111

(3-3.20) K(e) = PTN(e)P
Pl/h o -1/h 07
o o o o
- -1/h o l/h o
b o o o o_

 

 

q

P’h/B -1/2 h/B 1/2
_1 -1/2 1/h -1/2 -1/h
h/EV -1/2 h/B 1/2
_ 1/2 -1/h 1/2 l/h d

 

 

Once computed, the K(e)

can be superimposed to
give the global stiffness matrix K. This superposition
is guided by the necessity of matching the local

ée) with their prOper place in the global

coefficients g
Q vector. For the one-dimensional problem, the K‘e)
are superimposed by placing each successive 4)<4 matrix
over the previous one by moving it two entries down

the diagonal; that is, by adding the (1,1) entry of the
second matrix to the (3.3) entry of the first, and so
on. This is generally always performed by computer.

NOte that K is a symmetric band matrix, having a
bandwidth of seven. regardless of the dimension of K.

The dimension of K, of course, depends on the mesh size.

For elements of the domain which border on the boundary,

112

(e) can be used, but with two

the same formula for K
rows and corresponding columns deleted. For example,
the right-most element of the beam would have

q;8)= qée) = 00

by virtue of the clamped (essential)
boundary conditions. This corresponds to the fact that
there is no base function nonzero at x = %*, hence
no corresponding coefficient in the global vector 0.

K(e) for this right-most element contributes only

Hence,
its upper left corner 2)(2 submatrix into the global

matrix K.

3-4. Numerical Results for the Clamped Beam

Numerical results for the clamped beam verify the

expectations expressed in Chapters 2 and 3.

In each of Tables 1 through 7, the first two columns
give values of e and h, respectively, used in the
computations. Whenever more than one value is given
for e or h leading to a single computation, a
Richardson extrapolation has been performed using the
values listed in accordance with error forms given
in sections 3-1 and 3-2. The number of values listed
indicates the number of extrapolations carried out,

e.g. two values of h, say, 1/32 and 1/64, ‘would
be used for one extrapolation; five values, 1/4, 1/8,

1/16, 1/32, 1/64 ‘would indicate four extrapolations to

113
produce the final computed results in the right hand

columns.

Quantities indicated by $8 h and 'w are

30h
computed from we h and ‘we h by multiplying by
I l

(1+-e-1h2/12), as discussed in section 3-2. That is,

before extrapolation,

- 2
_ Q.
we,h- (1+ 12 )w

60h
e,h 12 e,h °

Functions are all evaluated at the beams midpoint.

Tables 1, 2, 4, and 6 give relative errors comparing

numerical solutions and exact solutions to problem (I).

As expected, errors relative to *e for both we,h
and $e,h in Table 1 appear to be quadratic in h.
After one extrapolation, these errors become quartic,
agreeing with the form predicted in 3-2. These same
patterns appear in the errors for we.h and wc'h for
the two symmetric loads p/D = l and p/D = 6(x). This

6, but fails for e g_2'15, indicating

holds for e > 2-
that the small values of e are causing non-quadratic
terms in the principal error to interact with the

leading term. we,h and we,h then lose all approximating
value (in fact they tend to zero as e tends to zero).

However, we h and we h maintain excellent accuracy
0 O

114
for all small values of e, and remain quadratic
in h, except for the most crude mesh sizes (see also

Tables 4 and 6).

Tables 2, 4, and 6 show that great improvement in
accuracy can be gained by using extrapolations on b.
In fact extrapolations using all available mesh sizes
extends the range of small a for which we'h and
we,h give reliable accuracy. For example, an error
of less than one percent can be achieved by we,h for

11

e as small as about 2' with mesh sizes to h = 1/64

‘with four extrapolations.

Tables 3, 5, and 7 give relative errors comparing
numerical solutions to problem (I) with exact solutions

to problem (C).

For large values of a, Table 3 shows that neither
we'h nor $3,h bears any resemblance to ‘w’, even
after extrapolating on 3, indicating an interaction
of the higher order terms with the linear term in the
asymptotic series described in 3-1. For small e,
$e,h approximates ‘wé well, with or without extra-
polations. However, we'h again loses its ability to
approximate W6 as 3 decreases. For 6 = 2-% we,h
is too large. Then it decreases past the exact solution
as e decreases. Hence there is some 3 which could

be considered optimal. This is of no use, however, without

some method of determining the optimal e value.

115
Similar results occur for the solutions
represented in Tables 5 and 7, except that for large
values of e, we,h and we'h can be used to approximate

w accurately simply by extrapolating once to remove

0

the linear dependence upon a. As shown in 3-1, when

p/D is symmetric, the dependence of W“ upon 6 is

S

linear.

As indicated by the error estimates of Chapter 2
and by the asymptotic formulas of Chapter 3, the results

show that U6 is a useful approximation to U6 and

,h

U but caution must be exercised when e is small.

0 VI
While Richardson-type extrapolations improve accuracy
and extend the range of small e for which results
remain reliable, limitations are still encountered.
-lh2

Because the factor (1+-e /12) multiplied by

the solution U removes the adverse effect of small

e,h
e, the quantity fie,h approximates both U6 and U0
much more reliably than does Ue,h over the entire
range of small 9. In the case of the clamped beam,
the limitations described above are overcome by this

algebraic adjustment to the numerical solutions.

The exact solutions to problems (C) and (I) for
the clamped beam under the three different loads
considered in 3-4 are given below, along with the

quantities (center deflections or slopes) used for

116

comparison in Tables 1-7. In each case, the domain is
l 1
-§SXS7

For p/D = 4 sinh 2x,

Gosh 1-3 sinh l) 3

WC = % sinh 2x+ 4 x- (cosh l-sinh l)x

 

w _e sinh 2x+ 123 [(-cosh 1+ Zél+4e)sa.nh l) x

i:
u

 

 

e - O l+lZe
+ (3 cosh 1-3-4 sinh l)x3]
_ I 126: (3 cosh 1-4 sinh l)
we ’ “’0“ 1+12¢l 4

+ (-3 cosh 1+4 sinh 1.)le

, -1 1 .3. - -
wo(0) —2+4 cosh 1-4 811111 1 :3 .43692635 E 2

—_ ’ _.3_€_ _ -
¢€(O)— WO(O)+1+125 (3 cosh l 4 sinh l) .

For p/D = l

wO = 5%.: (1-4x2)2
we = w0+% (l-4x2)
we = -w0
wow) = 3,1371- ~ .26041667 3-2

_ E
w€(0) — wo(0)+8 .

117

For p/D = 6(x),

W0:
£3_}_{3.+——1 O X '1;
12 6 192 ' S $2
e(-'2'+i') . -%gx< 0
we—wo=
1
_ l
we — 4"O
.._1 _ _
wO(O) .. 192 - .52083333 E 2

.9

w€(0) = w0(0) +4

118
Table 1.

Errors for beam problem (I),

p/D = 4 sinh 2x

 

6: h (Vim-“V“ Welh-‘bevfie
25 1/4 -.2531 E0 -.2530 E0
1/8 -.6342 E-l -.6344 E-l
1/16 -.1588 E-l -.1587 E-l
1/32 —.3971 E-2 -.3968 E-2
1/84 -.9928 E-3 -.9922 E-3
20 1/4 -.2s37 E0 -.2500 E0
1/8 -.6365 E-l -.6243 E-l
1/16 -.1593 E-l -.1561 E-l
1/32 -.3983 E-2 -.3902 E-2
1/64 -.9958 E-3 —.9754 E-3
2‘5 1/4 -.2996 E0 -.1829 E0
1/8 -.7972 E-l -.4137 E-l
1/16 -.2028 E-l -.1007 E-l
1/32 -.5091 E-2 -.2500 E—2
1/84 -.1274 E-2 -.624O E-3
2-15 1/4 -.9945 EC -.4906 E-l
1/8 —.9771 EC .7152 E-3
1/16 -.9142 E0 .9968 E-3
1/32 -.7212 E0 .3004 E-3
1/84 -.4000 E0 .7831 E-4
2‘60 1/4 -1.0- -.4868 E-l
1/8 —1.0 .8349 E-3
1/16 -1.0 .1028 E-2
1/32 -1.0 .3084 3-3
1/84 -1.0 .8031 E-4

Errors for beam problem (I), extrapolated, p/D = 4 sinh 2x

119

Table 2.

 

e h
25 1/84 -.9928
1/4. 1/8 -.2641
1/32, 1/64 -.6493
20 1/84 -.9958
1/u. 1/8 -.3094
1/32. 1/84 -.7574
2’5 1/84 -.1272
1/4, 1/8 -.6430
1/32. 1/84 -.1885
2‘15 1/84 -.4000
1/4, 1/8 -.9713
1/32, 1/64 -.2909
2‘60 1/64 -1.0
1/4. 1/8 -1.0
Laz,1flm -LO
2'5 1/84 -.1274
1/4. 1/8, 1/16.
1/32, 1/64 -.1145
2'6 1/84 -.1744
1/4, 1/8, 1/16.
1/32, 1/64 -.2688
2‘7 1/84 -.2863
1/4. 1/8. 1/16.
1/32. 1/64 -.6089
2‘8 1/84 -.5307
1/4' 1/8' 1/16'
1/32. 1/84 -.1153

E-3
E-3
E-7

E-3
E-3
E-7

E-2

E-9

E-2

E-8

E-2

E-7

E-Z

E-5

(fie h"V€)/W€ . (Eelh"V€)/We

 

-.9922

.2051
.5288

.5784
.1428

.7831
.1731
.4271

.8031
.1734
.4279

-.4439

-.1055

-.1265

E-3
E-3
E-7

E-3
E-4
E-8

E-
E-
E-

Ule.)

E-4
E-l
E-S

E-4
E-l
E-S
E-3

E-11

E-3

E—ll

E-3

E—lO

E-3

E-10

120
Table 3 .

Errors for beam problem (C), p/D = 4 sinh 2x

 

I I I I
5 h ”6.11 “wom‘b ”gm'wovwo
25 1/64 4.08 4.08
25.24 1/64 4.06 4.06
2"5 1/64 1.11 1.12
2"5,2"6 1/64 .173 .176
2"15 1/64 -.40 .1577 E-2
2"15.2"16 1/64 -.74 .8058 E-4
2"60 1/64 -1.0 .8031 E-4
2‘60,2‘61 1/64 -1.0 .8031 E-4
2'5,2‘6,2‘7,2‘8 1/64 -.9934 E—2 .7562 E-3
1/4,1/8,1/16.

1/32,1/64 .6729 E-3 .6762 E-3

121
Table 4.
p/D = 1

Errors for beam problem (I),

h-w€)/w€ (weth-w€)fi'€

 

h (w
6 0
1/54 -.6358
1/4, 1/8 -.6622
1/82, 1/54 -.9657
1/64 -.2034
1/4. 1/8 -.6738
1/82, 1/54 -.1663
1/54 -.6506
1/4, 1/8 -.5714
1/32, 1/54 -.1690
1/4. 1/8 -.9714
1/82, 1/54 -.2909
1/54 -1.0
1/4, 1/8 -1.0
1/32, 1/54 -1.0

-.6341
-.1349
-.8038

-.5708
-.2088

-.2387

-.1346
“01901

-.2107
-.1513
-.2761

E-ll
E-12
E-ll

E—ll
E-12
E-ll

E-11
E-12
E-10

E-8
E-ll
E-8

E-8
E-ll
E-8

122
Table 5.

ll
H

Errors for beam problem (C), p/D

 

e h (w£,h-w0)/WO (we,h-WO)NO
5
2 1/64 .1536 E4 .1536 E4
25.24 1/64 -.9785 E-3 .2627 E-8
1/4, 1/8 -.3055 E-4 -.1455 E-10
2"5 1/64 .1498 El .1500 El
2‘5.2‘6 1/64 ~ -.2925 E-2 —.1911 E-10
1/4, 1/8 -.5302 E-l -.4370 E-12
-15
2 1/64 -.3791 E0 .1465 E-2
2'15.2‘16 1/64 -.7431 E0 -.1056 E-8
1/4, 1/8 -.9997 EC -.8740 E-12
2‘60 1/64 -1.0 -.2107 E-8
2‘60,2‘61 1/64 -1.0 -.9481 E-9
1/4, 1/8 -1.0 p.2132 E-12

123

 

Table 6.
Errors for beam problem (I), p/D = 5(x)
h (we.h-w§)Ne (weLLh-wenws—
25 1/64 -.6358 E-6 -.5105 E-ll
1/4, 1/8 -.6621 E-8 -.1065 E-12
1/32, 1/64 -.8088 E-ll -.6476 E-ll
20 1/64 —.2034 E-4 -.4608 E-ll
1/4, 1/8 -.6738 E-S -.1531 E-12
1/32, 1/64 -.1661 E-8 -.5819 E-ll
2‘5 1/64 -.6506 E-3 -.7572 E-ll
1/4, 1/8 -.5741 E-2 -.2089 E-12
1/32, 1/64 -.1690 E-S -.9533 E-ll
2‘15 1/64 -.4000 E0 -.1380 E-8
1/32, 1/64 —.2909 E0 -.1811 E-8
2"60 1/64 -1.0 -.1997 E-B
1/4. 1/8 -1.0 -41428 E-ll
1/32. 1/64 -1.0 -.2617 E-8

124
Table 7.

5(x)

_ Errors for beam problem (C), p/D

 

6 h (w€,h -wO)NO (we '1'! -wO)NO
25 1/64 .1536 E4 .1536 84
25.24 1/64 -.9785 E-3 .2179 E-8
2'5 1/64 .1498 E1 .1500 81
2‘5,2‘6 1/64 -.2925 E-2 -.1803 E-lO

1/4, 1/8 -.5302 E-l -.3944 E-12

-15
2 1/64 -.3991 E0 .1465 E-2
2‘15,2‘16 1/64 -.7431 E0 -.1023 E-8

1/4, 1/8 -.9997 80 -.7461 E-12
2‘60 1/64 —1.0 -.1997 E-8
2‘60,2‘61 1/64 -1.0 -.8660 8-9

CHAPTER 4 - BENDING 0F CLAMPED PLATES
IN THE IMPROVED THEORY

This chapter studies the problems of the bending
of clamped circular and square plates. The solutions
Us of problem (I) for a circular plate under axi-
symmetric and non-axisymmetric loading are studied
analytically as 8 tends to zero. Numerical solutions
are obtained for U6 for a square plate under several
types of loadings. The behaviors of such solutions
as 8 tends to zero are noted and discussed, as they
serve to illustrate complexity encountered when

dealing with plates with non-smooth boundaries.

4-1. Asymptotic Analysis of the Solution Us

 

Because the main estimates of Chapter 2 guarantee
that solutions U€ of problem (I) converge to solutions
Uo of problem (C) (in the H1 sense), it is natural
to seek to examine this convergence asymptotically.

A glance at the governing equations reveals the

appearance of a singular perturbation problem (small

parameter multiplying the highest order derivatives).

125

126

Singular perturbation problems usually have solutions
exhibiting boundary layer behavior. Such a solution
‘would be quite regular throughout the domain, but would
change drastically in a narrow region near the boundary,
in order to satisfy the boundary conditions. Schemes
are available for handling one-dimensional boundary
regions. They usually involve two asymptotic series

for the solution, one valid near the boundary, one
valid far away from the boundary, and a matching process

to connect them.

Such schemes tend to be complicated even for
problems in one dimension. Little can then be expected
of such approaches for a two dimension problem, where
even the very boundary region thickness depends in some
unknown way upon the geometry of the boundary as well

as on the boundary value problem itself.

On a more hopeful note, the one-dimensional analog
of the clamped plate problem, namely the clamped beam
problem, exhibits no such boundary layer phenomena at
all. In an attempt to imitate the asymptotic analysis

of the clamped beam, consider a formal power series

. _ '1'
expansion of U6 - [¢x,¢yiwl .
U8 = Z E:iUi = Z eiwxi'w i'wi]T
i=0 i=0 y
where U = [W W w ]T 6 (H1)3
1 xi' yi' i 0 °

127

Substituting these into (1-4.6), assuming

sufficient smoothness,

or

-1 _
LBU€+ e LSUE: - F .

Collecting terms with like powers of 8 yields

(4-1.1) LSUb = 0
(4-1.2) LSUl = F-LBUO
(4-1.3) L302 = -LBU1
(4-1.4) LSU3 = -LBU2

and so on. Now

 

 

 

 

 

'- awo a F" 3
-w .. .——— 0
x0 0X
0w
L U - —¢ - -—JQ = 0
S 0 yo 0y
Mxo 5+ My + vzw 0
L- ax BY 0.; L d

The first two components say

128

(4‘15) U0= ‘87:” By '

and the third is identically satisfied as a consequence

of the first two.

Using (4-1.5). (4-1.2) is

 

 

P aw]. 1 " _§_ 2
”*x1 " 0x 3x V wOT
_ 3‘11 = .4. v2...
-¢yl BY by 0
a 5
L 0" BY 1.4 L d

 

 

Differentiating the first component by x, the second

by y, and adding to the third results in

Since Uo = [0,0,0] on an. U0 is in fact the solution

to problem (C) as expected.

Repeating the same process with (4-l.3),

 

84 31
2 x1 yl _
V ( 3x + by ) - O

returning now to take the Laplacian of the third

component in (4-1.2) yields

129

However U1 is over-determined on the boundary by

(4-l.2). Since U must satisfy le = w 1 ='w1 = O

1 y
on 00. the first two components of LSUl' as a vector,

can only be normal to the boundary. This is because

aw’ aw
T- T _._1 __1. ..
[wxl'wyl] - [0.0] , and 0x.' 0y -vwi has no

tangential component at the boundary since ‘wl = 0 there.

2

On the other hand, [a/Bx v‘w , a/By vszJT = v(v2wb)

0
in general will have both tangential and normal components

at the boundary. Thus U cannot meet the boundary

1
conditions, and the regular perturbation procedure

breaks down.

As pointed out in Chapter 2, the overly restrictive
condition of (4-l.2) is exactly that which could not
be met when attempting to apply Babuska's method for
improving the energy estimates in the general case.
The presence of a boundary layer is even more strongly

suggested.

The only two cases where (4-l.2) can be expected
to be satisfied are the clamped circular plate with
axisymmetric load, and the clamped beam. In these cases
the expansion can proceed. This was carried out for

the beam in Chapter 3.

The presence of a boundary layer can be verified
in one instance, namely that of a circular plate
with non-axisymmetric loading. This case with

g = cos 9 is carried out explicitly in 4-2.

130

4-2. The Clamped Circular Plate

The important special case of a clamped circular
plate is treated here in the improved theory (problem (I)).
Even for most simple geometries exact solutions of
problem (I) are rare, as are exact solutions to problem (C).
By enabling the construction of exact solutions, the
circular plate provides insights into the behavior of
solutions to both problems, which otherwise would be

missed.

In addition, the clamped circular plate with a
non-axisymmetric load allows first hand examination
of behavior which prevents improvement in the main energy
estimates of 2-2, and which also prevents a simple

asymptotic series expansion from being formed.

For solving the circular plate problems it is
necessary to transform the boundary value problems into
polar coordinates r and 6. Assume plate radius

unity.

Problem (C) becomes

U = -20 ‘1 .I‘awo w ]T
0 ar" r 06 ' O
. 4 - 2
‘7 w0 " D
(4-2.1) wo = 0
at r = 1
310 _. 0
ar 7

131

2 2
where v2=-L2-+-:ES%+-J-'2--a-2- .
0r r 06
Problem (I) is
T

 

 

 

 

 

2 2
5 W 5‘1! 3 ‘4!
83—; 15 :--— .
ar Zr 69 r
2
__2_3_‘”_e+ .13 *e
21:2 8 2r 81:86
2 2
3 1) 01! 3 V
+“lbw—2 2r+'_12—_ee'+'21?arag)
Zr 89 Zr
.. e‘l((r+g:) - 0
2 2
la ‘98 1 5"’8 1 a *8 1
(4'2'21’) 2 2 *5 r *7" 2 " 2%
at r 06 Zr
2
+ _3_.a_“'a 1 _1_ a ‘*'r
21.2 9 2r arae
2
3 B
+ H( _%_¢%_._l_.__w_e+_l_ ‘1’
0r 2r r 2r2 9

 

-1fi11n1183m312
(4-2.2c) e (arZ r 5r r2 09 r
5*
.1. l__9 _ E
+’r‘hr+r 9)- 7D
with Uh a [O,O,O]T at r = 1.

The first special case is to consider an axisymmetric
load. Then we a O, and all derivatives with respect

to- 9 can also be dr0pped from (4-2.l) and (4-2.2).
In (4-2.l)

2
(4-2.3) V2 =a—2'+
0r

Hue
9.19

In (4-2.2), (4-2.2b) is identically satisfied. The

remaining equations are

 

2
0 w aw
»r l. r' l ._ -1 BW’ _
(4-2'43) 01:2 + r ar - r2 1r e Wr+ar) - 0
2 av
4.9.3:. .16_W _1_=. .1. _ -2
(4-2.4b) e (arz + r r + :r + r (r) - D .

Note that the dependence on H vanishes.

Let 'E»= l. The solution to problem (C) is given

D
by
2
(4-2.5a) wo ='é% (1-r2)
aw
(4-2.5b) —-2 = £- (1-r2) .

ar 1

0‘

133

The solution to problem (I) with '§'= l is

given by
-_1_ 2 2 8 2 _ .6. 2
3w
_ _r_ 2 - _0

As in the case of the clamped beam with symmetric
loading, the dependence on e is simply linear. Also

the estimates of 2-3 agree with these solutions.
NOw consider the non-axisymmetric load '§'= cos a.
'The solution can be obtained to problem (C) in a

straightforward manner by assuming a form.'wo = R(r)cos a.

This separation of variables leads to solution

(4-2.7a) wO = €%(l-—r)2(2r+-l)cos 6
0w

(4-2.7b) 3:?- = j'—,(1)(1-r)(8r2 — r- 1)cos e
aw

(4-2.7c) 5—60 = -19%(l-r)2(2r+l)sin e

In order to solve problem (I) for load ‘§’= cos a,

it is helpful to recall from Chapter 1 that

(4-2.8) v4w = g-e v ‘3 .

134

(4-2.8) is valid in polar as well as rectangular
coordinates. ’Although a second boundary condition is
needed at r = l, in addition to w = O, a solution can

be found in terms of an undetermined constant a:

(4-2.9) w = r(r-l)[-a(r+ 1) +316(r-l)(2r+ 1) +-§r]cos e .
For convenience define the quantity

(4-2.10) C(r) = (r-l)[-a(r+1)+§%(r-l)(2r+1)+-§- r] .
That is,

(4-2.11) w = rC(r)cos 6

Also define dependent variables g, n by

(4-2.12a) t; cos a = ¢r+g—¥
(4-2.12b) 1. sin e = ¢e+%%¥ .

Substituting (4-2.11) and (4-2.12) into (4-2.2) results in

(44-139) r29” +r5’ + < -%+%- (1:2);
+ '1'? rn'+E—;3 n

r2(C+rC’)”+r(C+rC’)’-2rC’

135

(4-2.I3b) 1—59- r2 HEB

2 2
= -r2C”-3rC’
(4-2.13c) r§’+-g+-n = -er .
From (4-2.13c)
(4-2.l4) n = -r§’-—§-—er
(4-2.13a) reduces to
(4-2.15) r2§”4-3r§’-2%;;E3-§
= fiIr2(C+rC’)”+r(C+rC’)' -2rC’-(1-u)er].
Setting
(4-2.16) P = r;
and
(4-2.17) p = or, where a -——4£Z—-

=f1-TJE

(4—2.15) becomes the inhomogeneous Bessel Equation

136

(4-2.18) pZP”+-pP'-(1+-p2)P =

= epz[l:3-'-E e p2+8(-a-§6+%) [l—E'LJ" J— P- (1")4)€]-

A general solution to (4-2.18) is given as a particular

solution plus a solution to the homogeneous equation.

(4-2.19) P = PHOM+ PPART .

The homogeneous solution is of form

(4-2.20) PHOM = A11(p)4-BK1(p)

where I1 and K1 are modified Bessel functions whose
prOperties are well-documented (see [ 1]). A particular
solution is

(4-2.21) P = —2 8‘2 e pz55--8(-a--;-+-—)<1"'l

PART ‘3' 30 3 99'

From (4-2.16) through (2-2.21) g(r) is reconstructed,

and boundary conditions can be applied.

-1
(4-2.22) g(r) - r[A11(cr)+BK1(cr)-I-PPART] .
Since -% K1(or) is singular at the origin, it

does not belong to H1 (or even HO). Therefore, set

B = 0. Next, the edge condition at r = 1 is

(r(l) = 0 .

137

From (4-2.12a), it follows
5(1) = —2a+§ .

From (4-2.22),

 

-8ae+%:+§e 2 -2a

Now g(r) can be substituted back into (4-2.l4) to
determine n(r), where the edge condition 89(1) = 0
determines the constant a. The solutions g(r) and
n(r) can be substituted into (4—2.13b) to verify
that all three equations (4-2.l3) are satisfied. Edge
condition 86(1) = 0 implies n(l) =

I’ 1(0)

11 8 e2 __1_ _
S)GW+8€(-a-3O+3

_ £”_
0 - 3 (- 2a- 8ae4-lse

Using the Bessel function identity
I’(x) = I (x)-—l-I (x)
1 O x l

and solving for a,

10(0)
4 -1e__1_6 -162 532)

-3'0 3 a e+Il(c) (156+3e
I O(0) _1

(24-8€)'T—TET'- 20 '1-160 8

'where a = -—JZ——— .
Jl-u J?

138

Now define a by
(4-2.25) a = ca .
To summarize,
(4-2.26a) w = {go-(1-x)2(2r+ 1)
+er<1-r)(§(r+1)-i§-Hcos e

I (or)
_ 11—30 a 1 2r 4
(4'2'26b) ¢r+a 3r {a (_ 15 r TIE?)— - 3 + 15)

2 ~ I (or)
88 -3a+1 l _
+ ——3 (__r —Il(a) 1)}cos 9
(ar)
1.6'w _ fJ130a-11IO

(or)
30a- 11 I1 _ §r~4
" 6("_“151: —)'1 (c 15 )

:2 Io(cr)

m I1(a)

I 1(or)
+ —(1-3a)€2(m-+1)}Sin 9 .

“341-35) e «E

While (4-2.26) gives explicit solutions, the
dependence upon 8 is hidden in the quantities 5 and
a. It is therefore desirable to express (4-2.26)
asymptotically. Considering 8 a small parameter, then
a is large. Expansions of Bessel functions for large

argument satisfy

139

_livii 181-188 -9>
8x

 

               

J2vx 2'(8x)2
from this
I (x)
O l 3 l
(4-2.28) ~ 1+—+ —+ O(—)
Il(x) 2x 8x2 x3

using (4-2.28) in (4-2.24) and expanding, (4-2.25)

gives

(4-2.29) a =—15(11-faC—_H -e(8e+£(-g—‘Eln+me3/2>

Observing that the first term of w in (4-2.26a)
is

(4-2.30) . w = -r)2(2r+-1)cos 9

r
0 90‘1

(4-2.26a) can be rewritten asymptotically, with the

help of (4-2.29), to express comparison to *wo

(4-2.31a) ‘w = 60(l-r)(r+-1l)4-O(e 3/2) .

w0+ 3

Likewise, (4-2.31a) can be used to provide analgous

comparisons derived from (4-2.26b) and (4-2.26c)

0w0

(4-2.3lb) (r = "5—;- - fiu-rzmos 8

63/2 ”JIZ-lL-)cos 9
“1 Eiflf + 0(63/2)

15.5 r I

 

 

140

I (_JZ_;£_)
1 awO 8 O 1- 8 -
(4-2.31C) $9 = -1_" T- — .1.—5' 5111 e
Il(—._ACE._-_)
Jl-u J3

+ 5%(3 -r2)sin 9+ 0(63/2)

While the leading error terms in each of (4-2.3l)

are O(e) this is not the case for the radial derivative

of (4-2.31c). In fact, using identity 16(x) = Il(x),

3w

.3. .1.—2

(4-2.32) ar<we+r 06
I (L)

1
15Jl-p. I ( J2 )
1 Jl-HJ?

The first term on the right is only O(JE). This can

be seen using (4-2.27) to expand the ratio of Bessel
functions, for r bounded away from zero,

I (or) -a(1-r)
1 ~ e 3(l—r) 3 1-r)(11r+5) ,,,,
(4-2.33) _W" .E (1.. .. _L )

Bar 12802r2

 

indicating that

This J2" dependence, being confined to a boundary layer

region, means that

141

= O(eB/a)

(4—2.34) ”Us -U0||1

so estimate (2-2.9) cannot be improved from 51/2 to c.
It also indicates the presence of a boundary layer

phenomena, since near the boundary

:2 r _ J2 (l-rz
(4-2.35) 111(JIT'I) .. e NJ?
11(4L) J;
J'1-—fe'

Exactly at the boundary, this ratio of Bessel functions
is one. Away from the boundary, its value diminishes

1/2)

exponentially,as does the O(e error term it

multiplies. A boundary layer could be defined as that

region where the term of 0(61/2)

is significant, say
when the ratio in (4-2.25) is larger than 1/10. A
boundary layer so defined would have a thickness

proportional to .JE, based on (4-2.35).

5*
The quantity '35? represents the twisting moment.
For all geometries of clamped plates, and for all loads,
some twisting moment is expected, hence some boundary
layer phenomena may be present in the improved theory.
The only clamped plate which is an exception is the
circular plate with axisymmetric load. In 2-3, it was

shown that such phenomena cannot occur. The absence of

a twisting moment appears to be the reason.

142

4-3. Principal Error

The following is the analog of section 3-2 for
the clamped plate with piecewise linear finite elements
in the improved theory. The eXpansions for principal
error are considered only as far as the first term,

i.e. LhUh-LhU is considered. Consistency is shown,
1

depending on e- , The form of the leading term leaves
e-lh2
little hope that a factor such as (1+”_737—) ‘will

emerge. In fact, because this eXpression reduces to
that for the beam when all dependence upon the

-l 2
y-coordinate is removed, the factor (l4w3I52—0 is

the only conceivable one which could be present.

 

- 2
Numerical results suggest that computing (1+-€1:h )Uh-U

as was done for the beam does in fact seem to reduce a
large share of the error, but no claim can be made, as
was done in the beam case, that all adverse dependence

upon 3 is removed.

The differential operator, L, from (1-4.6)
is replaced by the finite difference Operator Lh,
generated variationally using linear finite elements. A
Taylor series expansion is then carried out to check
consistency and observe the form of the discretization

error, as was done in section 3-2.

143

The finite difference Operator is defined locally
in terms of the values of ul, uz, u3 at a central
node (subscripted O) and its neighboring nodes
(subscripted n, ne, e, s, sw, w, respectively. See

Figure 2 below.).

 

Figure 2. Nodes for local difference Operators.

 

By observing the three consecutive rows of the
global stiffness matrix K associated with a particular
node (the central node described above), a difference.
form can be found analogous to that determined for the
beam in 3-2. As With the beam, the difference form,

Lh, can be seen as a replacement for the differential
Operator, L. The following difference quotients arise

naturally in Lh as replacements for their differential

counterparts in L: define

-—- + + + + +
12[usw u‘w us 6uO un ue+une

D(1'O)[u] = -]'—[-u

6h ]

- + - + +
sw Zuw “3 un 2ue une

144

(4-3.1)

(0.1) _.JL
D [u] - 6h[-usw+ uW - 2us + 2un - ue + une]

D(2'O)[u] = —1-2-[uw-2u +ue]

h 0
(0.2) _ i
D [u] - h2[us--2uo+un]
(1.1) _ l _ _
D [u] - —2h2[usw uw us+2uO-un-ue+une] .

The Operator LhU = [L?U,L2U,L§U]T is defined by

the following replacements:

_ 1- .112 1
LlU - -ul,xx"( 2 )ul,yy ( 2 )uz'yx4-e (u14-u3'x)

is replaced by

(4-3.2a) LII‘U = -D(2'°)[ul] l1—3“-‘>D(°'2)[u11

+ e-l(D(O’O)[u1]4-D(1'O)[u3])

=-L:E _
L2U ( 2 )u2,xx u2.yy

— (igfl9u )

+u-l(u2+u

1oYX 30y

is replaced by

145

(4-3.2b) LgU = —(i§B)D(2'°)[u2] -D(0'2)[u2]

_ 1+3 (1,1)
( 2 )D [ul]

-1 D(O'O)[u2]4-D

[113])

__12
L3U — e {-V u3-ul'x-u2'y)

is replaced by

(4-3.2c) 1.13%; e-1(_D(2.0)[u3]_D(0.2)

.[u3]

_ D(100) [ul] -D(O'1)[u2]) .

To analyze the discretization error, it is first
necessary to eXpand the difference quotients of Lh
as formal Taylor series about the central node, in

powers of h. The replacements are (function evaluation

is at the central node)

(4-3.3)

D(o.o>[u]

=u+-6]= Z 11-—[(D(1'O)+D(O'1))k+D(k’O)‘t-D(O'k)]u

146

(4-3.3) (continued)

= D(1'O)u+.l'. Z L[(D(1'O)+D(o'l))k+l
3 k=2 (k+1):
k even
+ 2D(k.+l,0)-_D(O,k+l)]u
D(o,1)[u]
_D(0.1)u+_1_ 23 h" [(D(1,O)+D(O,l))k+l
— 3 k=2 (k+l):
k even
+ 2D(O,k+l)”Ema-1.0)]u
D(2,0)[u]
k
_ (2,0) h (k+2,0)
-D u+2 1:22 ———(k+2):D u
k even
-D(O'2)u+2 23 ———hk D(O'k+2)u
- k=2 (k+2):
k even
_D(1.1) +2 23 hk [(D(l,0)+D(O,l))k+2
- u k=2 (k+2):
k even

_ D(k+2.o) _D(O,k+2)]u .

147

For the purpose of generating principal error terms

as was done for the beam problem, the system

where

T
F — [f1,f2,f3]

is considered, instead of setting fl = f2 = O, which
occurs in the original problem (I). The more general
right hand side is needed for computing terms of the
error beyond the first, although they will not be

considered here.

Certain differentiation identities are useful in

simplifying the expression LhU. They are

(4-3.5)

1,x 2,y 3 1.x 2.37
4 _ 2
‘7 u3 " fa'ev f3+f1.x+fz.y
_ g - 2
ul+u3'x - efl+2[(l LUV u1+ (1+“)(ul,x+u2,y)’x]
_. .6. _ 2
u2+u3,y - ef2+2[(1 NV u2+(l+p)(u1'x+u2'y).y]
v2u +11 -+u = -6f
3 1,X 2,y 3 '

(4-3.

for

(4-3.

148

The following formal series result from using
5) and definitions (4-3.3) in the expressions

LhU given by (4-3.2)

6a)

Lyn - Llu

k=2

k 2,0
k -[(l-p.)(D( + )+D(0,k+2)) + (1+“)D(k+2,0)]u1

h—
k: (k+2)(k+l)

 

even

+

(1,0) +1)(0,1))k+2 (k+2,0) +1)(0,k+2))uz]

(1w) [-(D u2+ (D

 

2(k+2) (k-I-l)

 

6

(2.0)

(l-p)(D1'O)+D(O'1))k(D +D(O'2))u1

 

12

(k,0)+D(O,k) (2'0)+D(O'2))ul

(1-(1) (D ) (D

 

12

(2,0)u +D(1,1)u

(l,O)+D(0,l))k 2)
' 1

(1+u)(D (D

12

(1+“, (D(k.0) + D<o.k) ) (D<2.0)u + Dun)“l

12

1 2)

c-1(D(1'O)4-D(O'1))k[-(k-1)D(1'O)4-2D(O'1)]u3

 

6(k+l)

(k+1,0)

e‘1[-(k-3)D — (k+1)D(l'k) -2D(O'k+l)]u3

 

6(k+l)

149

(4-3.6b)

L211 — LZU

ID‘T

 

= Z
k=2
even

(1+-p.) [-(D

:3‘

---[(l “)(D(k+2,0)‘|_D(0,k«I-2))+(1"_H)13(0,}<-1-2)]L12
(k+2)(k+l)

(1,0)+D(0,l))k+2 (k+2,0)+D(0,k+2))u1]

ul-I- (D

 

+ 2(k+2) (k+1)

(1,0) +1)(o.1) k+D(k,0) +D(°'k)]f2

[(D )

 

 

 

 

 

 

6
(l-u)(D1°)D'(°'(1))kD(2'°) +1.,(0.2)m2
+
12
12
(1+“)(D(1,O) D.(o.k1)) (D(1,1)ul+D(o,2)u2)
+
12
(l+p)(D(k'O)+D(O'k))(D(1'1)u1+D(0'2)u2)
+ 12
-1.(D(1.0) D.(o.k1)) (2D(1,0) _ (k_1)D(o,1))113
+ 6(k+1)
+

 

-1[_(k_3)D(o'k+l) - (k+1)D(k'1) _ 2d(k+1'0)]U3
6(k+l)

150

 

(4-3.6c)
h
L3U-L3U
_ Z hk (D(l,0)+D(0,l))k+l(f +f )
"' _ 3.(k+1): 1 2
k-2
even

+ 1‘L£(D(1'O)+ D(0.1),1<+1(D(2.0) +D(o.2))(111+ u

2 2)

.1.-5E(D(1'0) + D(O'1) )k+1(D(200) +D(1'l))u1

4.
1+ (1.0) (0.1)k+1 (1.1) (0.2)
+'TE(D +D ) (D +D )u2
+ 2(BUGLO) +D(O.k+1))(fl+f2)_3(D(0.k+1)fl+D(o.k+1)f2)

+ (1-“) (D(k+l'0) +D(O'k+1) ) (D(2'o) +D‘0'2) ) (111+ u2)

- -3-(1-u) (D

2 (2'0)+D(0'2))(D(0'k+1)u

(k+l,0)
1+D u2)

+ (1+“) (D(k+1.0) +D(o.k+1) ) (D(2.0) +D(1.1) )“1

(k+1,0) +D(O'k+l) ) (D(1,1)

+ (1+u)(D +D(O'2))u2

(o.k+1) (D(2.0)u + D(1.1)u ,

3
' 2‘1“”) 1 2

(k+1,0) (D(1,1)u + D(0,2)u

' %(1+p)D 1 2)

_ e-1(D(1.0) + D(0.1) )k+2u3

151

(4-3.6c) (continued)

_ 22-1 (11:43.) (13(k+2,0)_{_D(0,k+2))113

+ E:

'1(D(1'k+1)+D(k+1'1))u3}
(4-3.6) is analogous to (3-2.20), the result of
Theorem 6. It can be observed that (4-3.6) shares with
1U and LgU-LZU
have leading terms which are free of e‘l, i.e. they

(3-2.20) the prOperty that L?U-L

are O(eohz). The 3.1 dependence in the leading term

is limited to the third component LgU-L3U.

An important difference is that while the right
side of (3-2.20) is expressed solely in terms of f1

and f the right side of (4-3.6) contains direct

20
dependence on U = [u1,u2,u3]T, the solution of LU = F,
in addition to f1, f2, and f3. ul, u2 and u3
cannot be removed entirely by identities, as was done

in deriving (3—2.20). (4-3.6) is eXpressed in a way

1)

which separates the terms of order O(e' from those

of order C(60).
In order to apply (4—3.6) to compute
LhUh-LhU = §-LhU, it is necessary to derive a formal

expansion for the components of B. By tedious integration,

it can be shown that

E = [O,O,f]T

152 _

where
~-i
(4-3.7) f - 2 If fm dA
h
k k .
h k+2 3 (k-j,j)
=f+2 23 .23 [(- )+(-1) ]D f
k=2 (k+3). j=0 3+1

k even

where T is the finite element "tent" function, with
value 1 at the central node, 0 at all others and

linear in-between.

h

An expression for L Uh-LhU can be written from

(4-3.6) and (4-3.7), in order to define the leading

principal error term, which is the 0(h2) term of

§ -LhU

(4-3.8) LhUh - LhU

As observed, the leading (0(h2)) term of each of the
first two components is free of 6.1. Since this

dependence is of concern, consider only this dependence

from the third component for the term proportional to h2:

153

(4-3.9)

2 -1

h e .1
6 (2 u3,xxxx+'u3,xxxy+2u3,xxyy+ 3,xyyy

2 -1
_ -h e .1 4

3,xxxy+ 3,xxyy+u3,xyyy

_..-_hie_:_1.

6

1

(ff-+1.1 )+O(eo) .

3.xxxy+u3.xxyy+u3.xyyy

At this point in the one dimensional problem the factor

(1+-e'lh2/12) could be employed to remove the 3'1 term

-hze-1f/12. In this case the additional mixed derivatives
remain. From this analysis, it is unclear whether the

use of this special factor has any advantage. It appears
that the adverse effect of small '3 on finite element
solutions is unavoidable. The apparent message here is
that very small a values should not be used in numerical
computation, as their presence may cause error which
cannot be controlled. Information about solutions of
problem (I) for small 6 values will have to come by

‘way of inferences from solving the problem for more

moderate values of the parameter.

154

4-4. Construction of the Element Stiffness Matrix

for the Clamped Plate in the Improved Theory

The construction Of the stiffness matrix for the
clamped plate is similar to that for the clamped beam
of Chapter 3, except for being much more complex, both
in the computation of the element stiffness matrix

and in the assembly of these into the global matrix.

The domain is subdivided into isosceles right
triangles of two types, based upon their orientation.
A local coordinate system will be assigned to each
triangular element of the domain having its origin at
the triangle's centroid. The legs of each triangle
will be h, the mesh size. The coordinates of the
vertices are indexed (locally) (X1,Y1), (X2,Y2),
(x3,Y3) starting at the right angle, and proceding
counterclockwise. The two types of elements are

given by Figures 3 and 4, respectively.

 

 

(x2,Y2) where (x131) = (331. -%)
(X2,Y2) = (313.33%)
P (3:333) = (-359. -%>
(X3,Y3) (X131)

Figure 3. Type 1 triangular element.

155

(Xl'Yl) (x3 .23)

 

 

' _ 1:. _h_
where (X1,Y1) — {-3 . 3)
h .Zh
(X oY ) = (-_I -—')
(X2,Y2) 2 2 3 3
_. 2.11 1

Figure 4. Type 2 triangular element.

Note that by replacing h by -h in the coordinates
for a type 1 element, the coordinates for a type 2
element are produced. In this way formulas deveIOped
for type 1 elements become valid for type 2 elements

by simply negating h.

The approximating finite element solution U to
the solution of problem (I) is expressed as a sum of

piecewise linear finite elements

where N is the number of parameters needed to represent

the solution over the entire domain, and

O O @i

156
to be determined by the indexing scheme. The function
”i is the finite element "tent" function generating
the piecewise linear solutions. It equals 1 at a
certain node (depending on the indexing), 0 at all
other nodes, and is otherwise linear. Its support is
the six elements adjacent to its central node (see

Figure 5).

 

 

 

 

 

 

Figure 5. Support for "tent" function.

As with the beam prOblem of 3-3, the finite element
solution is determined by finding the qi which minimize

the energy functional

(4-4.1) J(U)

B€(U,U) - 2PL(F,U)

QTKQ - 201']?

giving

ll
'11!

(4-4.2) KQ

157

(in the same notation as in 3-3). K and F are again

determined by summing the energy over the individual

elements
(4-4.3) J(U) = Z J(U(e))

= :3Q(e)TK(e)Qe._ZQ(e)T§(e) .
K(e) and §(e) are most conveniently found by

(e)

expressing U as a linear polynomial in the local
coordinates over a single element, then computing the
energy in terms of the coefficients, and finally

transforming to an expression in terms of the qi.

 

 

 

 

 

 

Let
F- 1 P P H
(e) 7
1.11 31 a4 a7 X
(4-4 4) U(e) = 11““ = a a a Y
' 2 2 5 8
(e)
u a a l
L 3 3 L 3 a6 9.. L _4 °

Let ”1' $2: $3 be the "tent" functions which are
non-zero at (X1,Yl), (X2,Y2), (X3,Y3), respectively.
Let q1,q2,...,q9 be coefficients, indexed locally,

described as follows

(4-4.5)

Equating (4-4.4) and (4-4.5) at the three nodes of

U(e)

the element,

(4-4.6)

 

 

 

 

Inverting gives

.ge)

uge>

 

 

L .3

 

 

P‘ (MD‘LMS‘

158

 

H‘ “4gztskr

h‘ (MD‘

 

2.5
3

 

 

 

 

 

 

 

 

 

159

 

 

 

 

O Dmd 5“.

(Dwatnhacnha

J

 

(4.4.7)
"a a a“ ' T ’P l
1 4 7 q1 ‘14 ‘17 h
a2 as as = q2 q5 q8 O
a a a q q q —-1-
3 6 9 3 6 9 h
L. d (_ .1 L
(4-4.7) can be eXpressed as
(4-4.8) A = 120(9)
where

 

 

T
A — [a1,a2,a3,a4,a5,a6,a7,a8,a9]
T
0(6) = [qloq20q3vq4oqsoq60q7oq80q9]
P1 . I 1 H
'1? I | O 1 -F I
-..._....l ...... I _____
- 1 ' l: '
P - h I : h I = O
c---—-l------——|-—---
1 . 1 1
_3 I I 3 I I 3 Id
1 0 O1
I = O l O

O O 1 .
L -

 

 

(e)

To determine K from

(4-4.9) B€(U(e)'U(e)) = Q<e)TK(e)Q(e)

160

first define a matrix N by computing

(4-4.1o) B€(U(e),U(e)) = ATNA .
Then
ATNA = 0(e)TK(e)0(e)
and
A = PQ(e)
imply

Q(e)TK(¢:‘.)Q(e) = Q(e)TP'I’m,Q(e) ’
that is
(4-4.11) K(e) = P NP .

This matrix multiplication of two nine by nine matrices
is best left to the computer. Importantly though,

it need be computed only once, and then used repeatedly
for the different elements of type 1 (for type 2,
negate h).

(e) 'U(e))

Computing B€(U involves integrals,

over the single elemental triangle, of form

Prs = If XrstXdY .

161

Formulas for Prs are available for triangular regions
with local coordinates (see Holland and Bell [10])

saving some tedious work here as well as in computing F.

Only the following are required for B€(U(e),U(e)):

h2
Poo=-2_' P01-1310=0
(4-4.12) 4 4
P = P =-h— P -'h—
20 oz 36 ' 11 ' 72 '

Inserting the polynomial forms for ul, u2, u3 from

(4-4.4),
(4-4.13)

B€(U(e) 'U(e))

_ 1+“ 2 1-u _ 2
“ 2 If (“1,x4'“2,y) dA“ 2 If (“1.x “2,y) dA
(e) (e)
.+.l:E ff (u d-u )ZdAi-e-lff(u +u )2+(u +u )2dA
2 lay 20x ‘ 31x 1 I 3.Y 2
(e) (e)
- -h—2-{a2+a2+2 a a +1211; a2+fl a2+ (l- )a a}
" 2 1 5 “ 1 s 2 4 2 2 “ 2 4
-1 2
+ €—-7h—{a§+ a; + 2a3a7 + a§+ ag+ 2a6a8]
e-lh4 2

{ai+a +a a +a2+a2+a2a5} .

36 4 1 4 2 5

162

(4-4.13) can be arranged in matrix form ATNA 'where

The result is

N (is symmetric.

NY)A

-1h4
36

 

 

where

LO

 

 

163

j

1 o o %- o o o o o
o 1 o o % o o o o
O o o o o o o o 0
-§- 0 o 1 o o o o 0
NY: 0 -21— o o 1 o o o o
o o o o o o o o o
o o o o o o o o o
o o o o o o o o o

0 O 0 O O O O O O
L 3

 

 

with N defined in (4-4.l4) and P defined in (4-4.8),

K(e) can be computed by

(e) T

K = P NP .

§(e) can be handled

The element load vector
in a similar manner although it will generally depend

on the position of the element in the domain

(e)T (e) _ (e) _ (e)
Q 17‘ -PL(F,U )—-H‘ E113 dA

 

 

(e)
x
_ .2
— [a3 a6 a9] (I D Y dA
(e) 1
’1 .1. 1“ ' “
'5' 'h 3 x
_1. .1. B
(e)
l o l 1
Lh 34 L _J

 

 

164

If the centroid of the elemental triangle is
(XC,YC) relative to the global coordinates (x,y),
then x=X+Xc and y=Y+Yc so p/D can be

expressed in local coordinates for the sake of

integration.
E.-.E. _.E
D - lD(x,y) — D(X+Xc, Y+Yc)
so
"(c-2) _ T
F - [Oiolfllololf2'0'0.f3]
where
(- H P l l 1‘1 P a
f1 h "h 3 X
- 1 .1. 2
f2 - O h 3 If D Y dA
1 l (e)
f3 5-? O -3- 1
L _J “ L- _J

 

 

 

 

 

 

is a function of (Xc,Yc).

If ‘5' happens to be polynomial, then the formulas
of Holland and Bell can be applied to the integration.
If numerical integration is called for, account should

be made of any error introduced.

The numerical studies of section 4-5 are carried
out for a square clamped plate under constant load, point
load, and a particular fourth degree polynomial load.

The plate covers the unit square with center at the

165

origin. Because all three loads considered are
symmetric in x and y, the problem is solved over
the first quadrant only.

The clamped edge condition imposes essential

_ __1_
3 - 0 at x - 2 and

boundary conditions 111 = u2 = u
at y ='%. Therefore no unknowns qi are associated
with nodes at the boundary. At each interior node,
three unknowns are present, namely qi. qi+1. qil+2
associated with ul, uz, u3 respectively. However,
the use of symmetry to restrict the problem to one
quadrant introduce constraints along the "edges" x = 0
and y = 0. IBecause 111 must be an odd function of x
and an even function of y, an essential "boundary

0. Similarly

condition" is imposed: u1 = 0 at x

u = 0 at y = 0. Thus along these two edges, each

2
node will possess only two unknowns instead of the
customary three. Special care must be taken to prOperly
index the variables (qi for i = 0 mod 3 no longer

represents just displacements u for example).

3'
An algorithm can be outlined for the indexing

procedure to be used.

1. Place nodes in a rectangular grid over
the quarter plate positioning three consecutive
indices at each node including those on the

boundaries of the quarter square. The elements

166

will be right triangles with diagonals
parallel to that of the full square. Work
horizontally to the right along each row to
reduce bandwidth in the stiffness matrix.
This original indexing information will be
retained so that each qi for i = 0 mod 3,
for example, will be readily associated with
a local index which is 0 mod 3 for the
purpose of assembling the element stiffness

matrices and element load vectors.

Systematically delete indices associated with
homogeneous essential boundary conditions,
storing the original index corresponding to
each new one. Call these "new" global

indices.

For each element, associate the nine local
indices to the corresponding original global
ones. Indices which have been deleted in

step 2 indicate rows and columns of the element
stiffness matrix which will not be entered

into the global matrix. Indices which have
survived the deletions of step 2 indicate
entries of the element stiffness matrix to be
entered into the global array. The positions

in the global matrix to receive the entries

167

from the local matrix are indicated by the
"new" global indices associated with the
surviving original indices. In this way,
elements adjacent to the boundaries of the
quarter square (and there are many) need
not be handled as Special cases (of which

there also would be many).

As the stiffness matrix is being assembled, the load
vector may be done simultaneously, in a similar manner.
For each element the centroid (Xc,Yc) is located and

the local loads computed.

The number of indices grows rapidly as the mesh
size shrinks. If h ='%% (n being the number of
elements along a side of the quarter square) the number

of original indices is 3(n+-l)2.

Even after deleting the boundary indices, the
number of "new" indices is still 0(n2) and so the
stiffness matrix has 0(n4) entries, many of which
are zero due to the bandwidth being O(n). Also the
matrix is symmetric, and while the fact alone can be
used to reduce the number of operations in assembling
the matrix it does not cut down on the size of the

array required.

To conserve on Operations in the matrix solving

routine (a Cholesky, or "square root" routine was used

168
for the symmetric positive definite banded stiffness
matrix) it is possible to make use of apriori knowledge
that certain entries are zero. This has been done

in the computations of 4-5.

More economical storage is also gained by arranging
the upper diagonals of the stiffness matrix
"rectangularly" in an array. In this way, storage
requirements are 0(n3) a substantial savings, but
at the cost of computer time and programming complexity.
This has also been done for the work in 4-5. The array
size required is (3n2-2n) by (3n4-2). For h = l/18,
this is 225 by 29. Naturally a step 4 in the process
is required to associate a pair of "new" indices with

a position in the rectangular array.

4-5. Numerical Results

Numerical solutions to problem (I), the clamped plate
problem, are computed below for a square plate under
three different load distributions. The computations
are carried out using piecewise linear finite elements,
as described in section 4-4. The plate has sides of
(dimensionless) length one. Mesh sizes of h = 1/4, l/B,
l/16, and 1/32 have been used for computing solutions for

each 6.

169
Deflections at the plate's center have been displayed

in Tables 8, 9, and 10, for 6 values ranging from 2"1

to 2.8. A Richardson extrapolation has been performed
using consecutive mesh sizes in accordance with

section 4-3. These appear at the right in the tables.

It should be noted that exact solutions of problem (I)
are not generally available for test prOblems generated
‘with simple load functions. Thus inferences as to the
reliability of solutions must come from the numerical
results themselves, rather than by computing errors
relative to exact solutions. As in the case of the
clamped beam of section 3-4, it is expected that numerical
solutions will fail to maintain accuracy when 6 becomes
very small (while the use of the correcting factor
(1+-e’lh2/12) seems to increase the range of e for
which the results are reliable, there is no solid
theoretical basis for its use here, in contrast with

the case of the beam).

The estimates of section 2-4 can be interpreted as
saying that in order to maintain a given accuracy in
the numerical approximation of solutions to problem (I),
as 6 becomes smaller, a finer and finer mesh size,
h, must be used. h is restricted from becoming too
small by the rapid growth of the system matrix. h = 1/32
produces 736 system variables distributed over the

quarter-square, and a stiffness matrix stored

170
"rectangularly" in an array, 736 by 50. In general,
for h = l/2n, an arraYo (3n2-2n) by (3n4-2) is
required. Thus, there exists a practical limitation on

the size of e for which accurate results can be obtained.

The following rather arbitrary criterion will be
applied to the center deflections in Tables 8, 9, and
10, to determine the reliability of the numerical

outputs:

The computed values for the center
deflection will be considered reliable
for all values of e for which there
is less than one percent difference ‘
between the two values computed using

a given mesh sizes and a mesh size half

as large, respectively.

In order to investigate the behavior of solutions
for smaller 6, a less restrictive criterion may be
allowed. Consider the center deflections computed with
three consecutive mesh sizes, say, h = 1/8, 1/16, and
1/32. As in the above criterion, a one percent
difference will be acceptable between the result of a
Richardson extrapolation using h = 1/8 and l/16,

and one computed using h = 1/16 and 1/32.

In all the computations in Tables 8, 9, and 10,
the value p = 0.3 is taken for Poisson's ratio. The

problem of the clamped plate is solved for a plate

171
covering the region -l/2 g x g 1/2 and -l/2 g y g l/2.
For the three loads tested, values can be given for

center deflections in the classical theory.

For the polynomial load

p/D = 24(x4+ 12x2y2+y4) - 36(x2+y2) + 5

the solution to problem (C) is

2 2

w = 2-8(4x -1)2(4y2-1)

O
The center deflection is then given by

wo(0,0) = 2‘8 = .390625 3-2 .

For a uniform load, p/D = 1000, a value for the
center deflection has been computed to three significant

digits (see [24]) to be
wo(0,0) = 1.26

For a point load, p/D = 1000 5(x) a center

deflection to three significant digits (see [24]) is

wo(0,0) = 5.60 .

172

Table 8.

Square plate deflections, polynomial load

 

h ‘w€(0,0) h . we(o.0)
1/6 .113116 E0

1/8 .113371 E0

l/16 .112286 E0 1/6. 1/16 .111925 E0
1/32 .111707 E0 1/16, 1/32 .111514 E0
1/& .578558 E-l

1/8 .585893 E-l

1/16 .581983 E-l 1/6. 1/16 .580680 3-1
1/32 .579396 E-l 1/16. 1/32 .578534 E-l
1/4 .301841 3-1

1/8 .311751 3-1

1/16 .311419 3-1 1/6, 1/16 .311308 E-l
1/32 .310471 3-1 1/16, 1/32 .310155 E-l
1/4 .162721 E-l

l/B .174236 E-l

1/16 .175911 E-l 1/8, 1/16 .176469 E-l
1/32 .175849 E-l 1/16, 1/32 .175829 E-l
1/4 .918538 E-2

1/8 .104783 E-l

1/16 .107770 E-l l/B, 1/16 .103799 3-1
1/32 .108283 E-l 1/16. 1/32 .108454 E-l
1/4 .544297 E-2

l/B .685689 3-2

l/16 .730850 E-2 1/6. 1/16 .745903 3-2
1/32 .741368 E-2 1/16, 1/32 .744874 E-2
1/4 .331845 E-2

1/8 .484342 E-2

1/16 .548202 E-2 l/B, 1/16 .569489 3-2
1/32 .566086 3-2 1/16. 1/32 .572047 3-2
1/2 .200441 E-2

1/8 .354541 3-2

1/16 .443066 E-2 1/8. 1/16 .472574 E-2
1/32 .472872 3-2 l/16, 1/32 .482807 3-2

173

Table 9.

Square plate deflections, uniform load

 

h w€(0,0) h we(0,0)
2'1 1/6 .400761 32

l/B .389808 32 '

l/16 .384852 32 1/6. 1/16 .383200 32

1/32 .383036 32 l/16, 1/32 .382431 32
2"2 1/4 .205338 32

l/B .201426 32

1/16 .199252 32 1/8, 1/16 .198527 32

1/32 .198396 32 1/16, 1/32 .198111 32
2‘3 1/6 .107468 32

1/8 .107151 32

1/16 .106406 32 1/8, 1/16 .106157 32

1/32 .106042 32 1/16. 1/32 .105920 32
2‘4 1/3 .582425 31

l/B .598567 31

1/16 .599000 31 1/8, 1/16 .599145 31

1/32 .598029 31 1/16, 1/32. .597706 31
2'5 l/h .331305 31

l/B .359290 31

1/16 .365063 31 1/8, 1/16 .366597 31

1/32 .365854 31 1/16, 1/32 .365727 31
2"6 1/4 .198141 31

1/8 .234978 31

1/16 .245902 31 1/6. 1/16 .249543 31

1/32 .248390 31 1/16. 1/32 .249220 31
2'7 1/4 .121882 31

l/B .165642 31

l/16 .183124 31 1/8, 1/16 .188951 31

_ 1/32 .187941 31 1/16, 1/32 .189600 31

2'8 1/4 .741327 30

l/B .121003 31

1/16 .147093 31 1/6. 1/16 .155791 31

1/32 .155821 31 1/16, 1/32 .158731 31

Square

174

Table 10.

plate deflections, point load

 

e h w§(0,0) h w§(0,0)

2'1 1/4 .191177 33

1/8 .249895 33

1/16 .305900 33 1/8. 1/16 .324568 33

1/32 .361238 33 l/16. 1/32 .379683 33
2'2 1/4 .973868 32

1/8 .127631 33

1/16 .155866 33 1/8, 1/16 .165278 33

1/32 .183578 33 1/16. 1/32 .192816 33
2‘3 1/4 .504340 32

1/B .664676 32

1/16 .808343 32 1/8, 1/16 .856232 32

1/32 .947386 32 1/16. 1/32 .993734 32
2"4 1/6 .268515 32

1/8 .358268 32

1/16 .432908 32 1/8, 1/16 .457788 32

1/32 .503007 32 1/16, 1/32 .526373 32
2‘5 1/4 .148785 32

1/8 .203995 32

1/16 .244711 32 1/8, 1/16 .258283 32

1/32 .230522 32 1/16, 1/32 .292459 32
2‘6 1/4 .861552 31

1/8 .125035 32

1/16 .149833 32 1/8, 1/16 .153099 32

1/32 .168855 32 1/16, 1/32 .175196 32
2"7 1/4 .513293 31

1/B .826778 31

1/16 .101134 32 1/6. 1/16 .107353 32

1/32 .112471 32 1/16, 1/32 .116233 32
2"8 1/4 .304324 31

1/8 .574003 31

1/16 .749797 31 1/8, 1/16 .808394 31

1/32 .835674 31 1/16, 1/32 .864300 31

175

Tables 8 and 9 indicate that for both a uniform
load and the polynomial load, changes in the center
deflections for problem (I) occurring between mesh sizes
h = 1/16 and h = 1/32 are less than one percent for
6.2 2"5 ‘while changes of 1.5% and 1.01% respectively
occur at e = 2'6. indicating a need for finer mesh
sizes in order to regain accuracy. Under the first

‘5 and above would be considered a

criterion, 3 = 2
range of e for which numerical results are expected to

be reliable for the two loads mentioned.

While one percent may be considered a rough estimate
of the accuracy of the numerical results relative to the
exact solutions of prOblem (I), improved accuracy can be
expected when these "reliable" values are subjected to a

Richardson extrapolation.

By using the second criterion to compare previously
extrapolated values, the range of e for which reliable

results can be expected is extended to 6‘2 2-7.

Table 10 gives the center deflections under a
point load for problem (I). No such indication of
convergence is seen as the mesh is refined, for any range
of e. The center deflections steadily increase as the
mesh size shrinks. This, however, is not due to the
adverse effect of small 3. It indicates the failure
of the finite element method to approximate the singular
nature of the exact solution of problem (I) under a point

load.

176

‘As shown in Chapter 1, the displacement component,
w, of the solution to problem (I) can be partially
decoupled from the lepe components (x and 3y,

leading to the differential equation
4 2
(4-5.1) V‘w = p/D-ev p/D .

Viewed as a distribution, with p/D prOportional to
5(xgyh ‘w is expected to have a term proportional to 6.

‘with a logarithmic singularity at (0,0).

Physically, an infinite shear occurs under a point
load, and an infinite displacement results from this
infinite shear. Table 10 reflects the failure of the
finite element method with polynomial base functions to
accurately describe this singular behavior under the

point load.

It should also be mentioned that, in contrast, the

governing equation for problem (C),

leads to a singularity in W’ of form r2]J1r, 'when

0
p/D represents a point load. The classical theory results

in continuous displacements in response to a point load.

The question of approximating solutions of problem
(C) by numerical solutions of problem (I) poses several

obstacles.

177
While the estimates of Chapter 2 indicate that

solutions of problem (I) converge to solutions of problem
(C) as 3 tends to zero, this convergence cannot be
directly exploited numerically because of the limitations
just discussed. The numerical solutions of problem (I)
needed for this convergence are for small values of 3.
But this is exactly the range of 3 for which these

numerical solutions are unreliable.

The question then becomes: can approximations

to solutions of problem (C) be inferred from numerical
solutions of problem (I), using moderate sizes of e?

To address this question, it is necessary to postulate

the form through which solutions of problem (I) depend

on e. for small 6. Correctly postulated, this form

can be used to extrapolate the numerical data produced

for moderate values of e to produce a value corresponding

to e=0.

Without the help of an asymptotic expansion,
different assumptions may be made leading to quite
different conclusions. The only guideline is that
any assumed form must satisfy the estimates of section 2-2,
namely, that the dependence of U€-U on 3 cannot

0
be "more singular" than 31/2.

In fact, even this may be violated on a domain
with corners. While estimate (2-2.9) holds for domains

with non-smooth boundaries, the right hand side contains

178
the factor HV(V2WO)“O 'which may be infinite. In
the case of the particular polynomial load being consider
numerically in this section, WC is also a polynomial,

so that estimate (2-2.9) is binding on the solution,

Us' to problem (I).

However solutions to problem (C) in general exhibit
singular behavior at the corners of the square, even

for smooth loads (including uniform loads).

Kondrat'ev [13] studied the effect of corners on
solutions of elliptic boundary value problems, finding
that, while solutions are analytic up to regions of
the boundary which are smooth, singularities are introduced
into solutions at angular points of the boundary, the

severity of singularity depending upon the angle.

For the biharmonic equation with homogeneous
boundary conditions (such as for prOblem (C)) Kondrat'ev
deveIOps an equation to determine the singularity at a
corner, as a power of r, the radial distance from the
corner. The severity of the singular depends only on
the angle of the corner, not on the forcing term of the
differential equation. For a right angle corner, the

singularity has the form

2+1

where z is a complex solution to the equation

179

sin-Egg = 22 .

Besides the integer solutions 2 = 0 and z = 11

it can be shown that solutions lie between

3
Z< Re(z) < 1
and

-l < Re(z) < "% .

Since the integers represent analytic solutions

the solutions which may be singular have exponents either

%< Re(z+1) < 2
or

0 < Re(z+-l) < %-.

The latter is rejected because such a singularity will

prevent the solution from being in H3. The former

however provides a singularity of form

7/4-P5

w~r where 0<5<;]1'-.

Roughly Speaking, then

180

BW’ 3/44-5
-—— ~ r
Br

3.231 r-l/4 + 6

BIZ

63w

BIB

~

r-5/4-+6

‘while such a solution has square integrable second
derivatives (required for prOblem (C)), the third
derivative has infinite energy. The third derivative

corresponds to the shear force.

It may be conjectured that while this shear
singularity is confined to the corner in the restrictive
setting of problem (C), this same shear effect is
released in the less restrictive problem (I). No longer
a singularity confined to a corner, its effect spreads

throughout the boundary as a boundary layer phenomena.

Being no longer confined to isolated points, its
effect is also more readily propagated to the region of
the domain away from the boundary. This may account
somewhat for the difficulty in using numerical results
for problem (I) in order to approximate solutions to

problem.(C).

Even in the cases where Hv(v2wo)Ho is finite,
some irregular dependence upon 3 can be expected. In
the case of the circular plate with p/D = cos a, from

section 4-2, it was illustrated that a boundary layer

181
behavior occurs, though this would not show up in the
center deflection. Without the symmetry of a circular
domain, it is unclear whether or not such behavior would

be confined to the boundary region.

The simplest approach to the question of approximating
solutions of problem (C) is to assume a simple power
series in e:
i

U -U = Z a. e o
e 0 i=1 1

While this cannot be expected to describe the solutions
over the entire domain it may’be valid at least at the
center of the square, or perhaps only for the displacement

component.

Based on this simple form (see ‘w on Figure 6).

6A
Richardson extrapolations yield the following results:
For the polynomial load, extrapolation using

-1 -2 -3 -4 -5

e = 2 , 2 , 2 , 2 , 2 produce a center deflection

Since the exact solution to problem (C) is

wk is in error by 4.6%.

The values Of 3 above were chosen to satisfy

the first criterion for reliability. A slightly more

182
accurate result is obtained by using smaller 3 values
(satisfying the second criterion), 3 = 2'4, 2’5, 2 , 2 .
This produces a center deflection

WA = .396349 E-2 ,

which is in error by 1.5%. Clearly, caution must be
exercised in interpreting results based on such ad hoc

assumptions.

A second approach is to make inferences about the

e dependence from the numerical data itself.

First, it may be Observed that, while problem (I)
does not in general admit a regular series expansion
for small 3. it does so easily for large 3. The form
is

-i
Us — e V-1+Vo+i§1 6 vi .

The functions Vi are independent of e. and have no
direct relationship to the terms (if they exist) in an
expansion of Us for small 3. However this form may
serve as a framework to help clarify the behavior of Us

for the moderate values of 3 being used to approximate

U0. Let w: be the third component of eV_1+VO.

This "outer" expansion indicates an almost linear
relationship between U3 and e for e sufficiently
large. This is observed numerically in the center

deflections for e as small as about 2'3 (see Figure 6).

183
Since numerical values for a large are quite
reliable, the linear portion of the graph of Figure 6 can
be determined accurately. The non-linear region near

3 = 0 can be viewed as a departure from the linear form.

It is noted that in the nonlinear region for which

3 is still large enough to produce reliable results
(2-7 g_e g_2-3, roughly), there is a very good
correlation between the values ‘wéo)(0,0)-rwe(0,0) and
an exponential of the form. C1 e-CZJE. These constants
can be determined by taking logarithms and applying a
linear regression. The value of C1 will be the
extrapolated value wB (see Figure 6), which approximates
wb.

It is observed that the constants C1 and C2 are
very sensitive to the data used to produce them.
Moreover, while the numerical values used to construct
the linear portion of the graph are generally very
accurate relative to their size, their precision may
.520)

not be adequate to produce values, , ‘which are

equally accurate for small 6.

This sensitivity makes it difficult to attain
consistently reliable results, and again caution must
be exercised in interpreting the results. Generally,
values for center deflections extrapolated to e = 0
are within 2-3% of the center deflection for the classical

theory (generally being on the low side).

184
Qualitatively, the sensitivity of this method may
be due to the large (in fact infinite) slope at e = 0
for the curve, 'weB, following the exponential form.
Reflecting an O(fe?) dependence, the curve will tend to
be displaced disprOportionately at e = 0, when slight
changes are introduced at 6 values in the moderate

range.

Whether or not this assumed exponential form is
correct, such sensitive behavior is allowed by the main

estimate (2-2.9).

185

 

 

 

Figure 6. Center deflections extrapolated to 3 = 0.

CHAPTER 5 - CONCLUSIONS AND FURTHER REMARKS

The relationship between the classical theory of
plates and the improved theory of Timoshenko and Reissner
has been investigated for the clamped plate, in terms
of the small parameter 3. a measure of the shear
effects included in the improved theory. As 3 tends
to zero, the solutions in the improved theory tend to

solutions in the classical theory.

While this continuity is demonstrated by the
estimate (2-2.9), strong evidence is offered to show that
the 3 dependence is not analytic in 3. nor is it
uniform over the surface of the plate. Rather, a boundary
layer occurs, possessing different 3 dependence than
is present throughout the plate. There is reason to
believe that this is related to the shear distribution
about the boundary in the solutions in the classical
theory. For plates with corners it may also be tied
to the shear concentrations which are known to occur
at the corners in the classical theory. An explanation
of these relationships might shed light on the role of

shear forces in both theories.

186

187
While solutions to problems in the improved theory
are rarely attainable in closed form for even very
simple loads and for domains with very simple geometry,
solutions developed by other means, such as Fourier
series, may be used to clarify the possible variations in
the 3 dependence over the plate surface, as well as

to determine the sharpness of (2-2.9).

It would also be useful to develop an estimate
analogous to (2-2.9), but estimating the expression

HU€-U The example of section 4-2 shows that

o‘lo-

HUe-U = 0(33/u) because of a term in a derivative

0“ 1
Of U3' If» "Ue"Uo"0 = 0(3). then the numerical
difficulties are lessened, in approximating center
deflections for the classical theory by those for moderate
values of 3 in the improved theory. It still must

be remembered that of primary interest to engineers

are the bending moments, which are estimated in (2-2.9).

By formulating the clamped plate problem in the
improved theory, it is possible to approximate solutions
using piecewise linear finite elements instead of the
more complicated elements generally required to solve
plate bending problems in the classical theory. However,
the accuracy of the results are adversely affected by
the small values usually chosen for the shear parameter
3. While this "spoiling effect" Of small 3 cannot be
avoided, in general, it may be possible to lessen its‘

effect.

188
Through the term 3-1PS(U,U) in the energy expression
which is minimized by the solution Ue to the clamped
plate problem in the improved theory, it is clear that
each candidate function U is "penalized" according

m (:1) 2

dA. The
2

numerical approximation U3 h to U3 minimizes the
I

 

to the value of PS(U,U)

 

same energy expression over a finite subspace, hence is
similarly "penalized". Since all three components of
U3,h are constructed with piecewise linear elements,
vu3 is constant within any single element. For a

fixed mesh size, h, as 3 tends to zero,
2
u

-l l

3 (( vu3+ (“2)
n1
unbounded unless (u
2

over the element. Then continuity requirements at the

dA over each element must become

 

 

) tends to a constant distribution

boundaries of the elements, tend to force

u
(“2) to the same values for all elements. ‘With

homogeneous boundary conditions at the plates edge,

u
(111) is forced to zero, and likewise, u3 must follow
2

suit.

It is tempting to blame this "mismatch" of polynomial

u1

“2
of the numerical results as 3 tends to zero. The

degrees in the eXpression vu34-( ) for the behavior

values are, in fact, "driven to zero" as 3 4 0.

189
Unfortunately, this phenomenon is inherent in the
.penalty function method form, as can be seen from a

simple algebraic example.

Consider the problem: Minimize

J = (x--l)2-I-3-1 y2

over 32 . Clearly the minimum occurs at the point (1,0)

for all 3 values. Consider the same minimization
problem restricted to the subspace y = hx. Of course
if h =_0, the minimum over the entire space R2
'would belong to the subspace, y = 0. Otherwise, a
"best" approximation to (1,0) would be found by
minimizing

-1 2 2

(x-1)2+3 hx .

The minimum along y = hx is (1/(1+ 3-lh2), h/(l+ (1112)).
For a fixed h, as 3 4 0 this point is forced toward
the origin, not to a location close to (1,0). The
measurement of closeness "penalizes" the y coordinate
much more than it does the distance from (1,0).

The similarity of the minimizing point above and

3'1h2

the factor (l+-—I§—fl, which proved useful in controlling
error for the clamped beam problem cannot be overlooked.

If the coordinates of the minimizing point are multiplied

190

by (1+-3-lh2) the resulting point is (l,h), a very
good approximation to (1,0) for small h, and this

approximating point is not adversely affected by small 3 .

Two questions thus arise related to the clamped
plate problem. First, can reliable results be obtained
for smaller values of 3 by approximating the third
component of U6 by piecewise quadratic elements instead
of linear ones? The estimates of section 2-4 suggest
that improvement will take place. Correcting the
"mismatch" mentioned above may improve accuracy even
more so. The second question is whether by approximating
by quadratic or higher degree finite elements, a
factor such as (1+-3-1h2/12) can be found to alleviate
the adverse 3 dependence altogether as was done in

the case of the beam.

Since the improved plate theory has been chosen
as a penalty function method approach to solving the
clamped plate problem in the classical theory, in order
to both avoid restrictive smoothness requirements and
bypass the "patch test" usually required Of solutions
failing these requirements, a more general question is:
can the patch test itself be relaxed by virtue of a
penalty function method? If such an approach can be
made rigorous, it would Open the door to the use of many
types of finite elements applied to higher order

equations.

191

The problem of extending this work to the cases
of simply-supported plates and plates with free edges
is hindered by the failure of the bilinear energy form
B€(U,U) to be coercive. The problem can perhaps be
approached by using a device of Cornwell and Yen [5],
namely, by attaching torsional springs to the plate
surface. As the spring constant is reduced to zero, the
simply-supported (or free-edged) plate is recovered. In
this limit, it must be shown that estimates similar to

(2-2.8) and (2-2.9) are recovered.

In order to successfully apply the penalty
function method to solve the problem of bending of a
clamped plate in the classical theory by using numerical
solutions in the improved theory, either a better
understanding must be obtained for the behavior of
solutions for small 3, or a better numerical process
must be developed which is not so adversely affected

by small values of 3.

BIBLIOGRAPHY

BIBLIOGRAPHY

Abramowitz. M., and I.A. Stegun, eds., "Handbook
of Mathematical Functions", Dover, (1965).

Aziz, A.K., ed., "The Mathematical Foundations of
the Finite Element Method with Applications to
Partial Differential Equations", Academic Press,
(1972).

Babuska, I., and A.K. Aziz, "Survey Lectures on
the Mathematical Foundations of the Finite
Element Method", in "The Mathematical Foundations
of the Finite Element Method with Applications
to Partial Differential Equations", Academic
Press, (1972).

Balakrishnan, A.V.. "A Computational Approach to
the Maximum Principle", J. Computer and
System Sciences (1971), Vol. 5, pp. 163-191.

Cornwell, P.E., and David H.Y. Yen, "Boundary
Value Problems in the Improved Theory of
Elastic Plates. I: Existence of Eigenvibrations
for Plates of Arbitrary Shape", SIAM J. Appl.
Math. (1976), Vol. 30, NO. 3. pp. 469—482.

Courant and Hilbert, "Methods of Mathematical
Physics", Vol. 1 and 2, Interscience.

Forsythe, G., and C.B. Moler, "Computer Solution
of Linear Algebraic Systems", Prentice—Hall,
(1967).

Garding, L., "Dirichlet's Problem for Linear Elliptic
Partial Differential Equations", Math.
Scandinavica (1953), Vol. 1. pp. 55-72.

Gilbarg, D., and N.S. Trudinger, "Elliptic Partial

Differential Equations of Second Order",
Springer-Verlag, (1977).

192

10.

11.

12.

13.

14.

15.

16.

17.

18.

19.

20.

21.

22.

193

Holland, I., and K. Bell, eds., "Finite Element
Methods in Stress Analysis", Tapir, (1969).

Joyce, D.C., "Survey of Extrapolation Processes
in Numerical Analysis", SIAM Review (1971),
Vol. 13, pp. 435-490.

Kirchhoff, G., "VOrlesungen fiber mathematische
Physik", Vol. l, 8.6. Teubner, (1876).

Kondrat'ev, V.A.. "Boundary Problems for Elliptic
Equations with Conical or Angular Points",
Trans. Moscow Math. Soc. (1967), Vol. 16,
pp. 227-313.

Love, A.E.H., "A Treatice on the Mathematical
Theory of Elasticity", Dover Publications, (1944).

Mikhlin, S.G.. "The Problem of the Minimum of a
Quadratic Functional", Holden-Day, (1965).

Mindlin, R.D., "Influence of Rotatory Inertia and
Shear on Flexural Motions of IsotrOpic Elastic
Plates", J. Appl. Mech., Vol. 18, (1951),
pp. 31-38.

Nayfeh, A., "Perturbation Methods", John Wiley and
Sons, (1973).

Reddy, J.N., "0n the Accuracy and Existence of
Solutions to Primative Variable Models of
Viscons Incompressible Fluids", Int. J.
Engng. Sci. (1978). Vol. 16. pp. 921-929.

Reissner, E., "The Effect of Transverse Shear
Deformation on the Bending of Elastic Plates",
J. Appl. Mech., Trans. ASME, Vol. 67, (1945),
p. A-69.

Stakgold, I., "Boundary Value Problems of Mathematical
Physics", Vol. 1 and 2, MacMillan.

Strang, G., and G.J. Fix, "An Analysis of the Finite
Element Method", Prentice-Hall, (1973).

Szilard, R., "Theory and Analysis of Plates:
Classical and Numerical Methods", Prentice-Hall,
(1974).

194

23. Timoshenko, 8., "0n the Correction for Shear of
the Differential Equation for Transverse
Vibrations of Prismatic Bars", Phil. Mag.,
Series 6, Vol. 41, (1921), pp. 744-746.

24. Timoshenko, 8., "Theory of Plates and Shells",
lst ed., McGraw-Hill, (1940).

25. Westbrook, D.R., "A Variational Principle with
Applications in Finite Elements", J. Inst.
Maths. Applics. (1974), Vol. 14, p. 79-82.