explore the methods of finite strip, finite prism and finite layer for
the analysis of isotropic and laminated composite plate-type structures.
A new finite strip method

is developed for the analysis of thin

plates , 2-D elastic problems and folded-plate structures. Unlike other
finite strips, the displacement functions of the present finite strip
method are expressed as the product of the usual 1-D shape functions in
the

transverse direction and

longitudinal direction.These

a

set

of

computed

static modes

can

be

static modes
easily

in

the

obtained

by-

means of a continuous-beam computer program. The resulting finite strip
possesses a few superior features when compared with the semi-analytical
finite strip and the spline finite strip. These features include, just
to

mention

a

few,

the

capability

of

dealing

with

geometric

discontinuities, multi-point supports, and the ease of implementation.
The original spline finite strip is restricted to the analysis of
thin plates or thin-walled structures.

In this thesis, a third-order

spline finite strip is formulated using a third-order plate theory. This

Introduction
Formulation
Numerical Examples
Pre』iwinery assessment : Static anaiysis of beams
A square plate with span a and subjected to a
point load
A square plate with span a and subjected to UDL
A 3-bay flat-plate wode』
A S-bay flat-slab model
A 5-bay flat-slab model
A continuous-slab bridge
rnnn 1 ncii on

structure in which one opposite pair of the sides of such a subdomain
are in coincidence with the boundaries of the structure. The width of a
strip remains constant along the longitudinal direction.One dimensional
finite element shape functions are used in one direction while a series
of trial functions is adopted in the orthogonal direction.
Finite prism (Fig.1.2) : It is a three-dimensional subdomain of a
structure in which one opposite pair of the faces of such a subdomain
are

in coincidence with the boundaries of

the structure.

The cross-

section of a prism remains constant along the longitudinal direction.
Two-dimensional , isoparametric, finite-element shape functions are used
within the cross-section of a prism while a series of trial functions
is adopted in the orthogonal direction.
Finite layer

(Fig.1.3):It

is a three-dimensional subdomain of a

structure in which two opposite pairs of the faces of such a subdomain
are in coincidence with the boundaries of the structure. The span and
width of a layer remain constant across the thickness. One-dimensional
finite-element shape functions are used across the thickness of a layer
while

a

series

of

trial

functions

Is

adopted

in

each

of

the

two

orthogonal direction.
Throughout this thesis, readers are assumed to be familiar with the
basic terminologies, concepts and procedures of the displacement-based
finite

fibrous compositesy which consist of fibers of one material in a matrix
material

of another;

particles

of

one

particulate

material

in

a

composites
matrix

of

another;

and

laminated

composites, which are made of layers of different materials. A laminated
composite plate is made of layers of the same orthotropic material that
are bonded together with the material coordinates of each layer oriented
differently with respect to the plate coordinates.
1.3 General review of finite strip,finite prism and finite layer
1.3.1 Finite Strip
The finite strip method pioneered by Y.K.Cheung[3], is an efficient
analysis tool for plate-type structures with regular geometry and simple
boundary conditions. The method is one of the best choices available to
an

-continuity for plate-bending problem and (6) ease of implementation.
In addition, the basic concept of the method is simple and involves no
complicated mathematics.
Since the first paper of finite strip[4] appeared in 1968, intense
interest

has

been

analysis. Several

focused

different

on

the

approaches

development
have

been

of

finite

developed,

strip
and

an

enormous literature now exists. The state of the art has been recently

s h o u l d be n o t e d t h a t a X - s p l i n e f i n i t e s t r i p was a l s o d e v e l o p e d b y Chong
et.aă&#x20AC;?.[18,19] for irregular shaped plates.
In addition to the beam vibration functions and spline functions,
two other types of longitudinal trial functions have been developed :(i)
Trigonometric series are superimposed on beam shape functions in such a
way that plates with general boundary conditions can be dealt with [20].
Owing to orthogonality relationships of the trigonometric functions used
in the formulation, partial uncoupling of the equations can be achieved,
and

this

can

be

taken advantage

Gaussian elimination process.

of by using

condensation prior

to

(ii) Orthogonal polynomials generated by

the Gram-Schmidt process [21]. The method was applied to free vibration
of plates with general shapes and boundary conditions. The performance
in terms of versatility, efficiency and accuracy is, at most, comparable
to the semi-analytical finite strip method.
(2)

Various plate theories for finite strips
The semi-analytical finite strips were originally formulated within

the context of Kirchhoff's thin-plate theory. They were later extended
to the analysis of moderately thick plates using the first-order shear
deformation

theory,

i.e.

Mindlinâ&#x20AC;&#x2122;s

plate

theory.

In

Hindiin-type finite strip has been successfully applied

recent
to

years,

the shape

optimization of variable thickness prismatic folded plates and curved
shells[22] and to the study of boundary layer effects in the vicinity of
edges and vertices of rectangular and sectorial plates[129].

elastic structure is often extremely tedious and complicated to find, if
not

impossible.

In

addition,

the

solution

is

often

restricted

to

structures with very simple geometry and boundary conditions. Although
the

standard

finite

three-dimensional
computationally
three-dimensional

element

method

solution

to

the

expensive.

To

reduce

finite

element

can

generate

problem,
the

analysis,

this

an
method

computational
the

approximate

finite

is
cost

prism

usually
of

a

method

[31-34] was devised for cost-effective analysis of prismatic structures
with two simply supported ends. For prismatic structures with free ends,
a special finite prism method was proposed by Fu and Yang[35]. In this
method, the end conditions are satisfied in an iterative manner using an
equilibrium approach. Additionally, the usual trigonometric series are

complicated formulation and implementation and therefore has not gained
popularity.
In addition

to conventional 2-D finite element

shape functions,

spline functions of various order can be used within the cross-section
of a finite prism[37].
Finite prism
system

or

can be

a Cartesian

formulated

coordinate

within

system.

a

An

cylindrical
example

of

coordinate
cylindrical

finite prism is due to Yoseph et.al.[38] who proposed a hybrid finite
prism for the analysis of singular stress fields in the vicinity of a
circular cut-out in a laminated plate.

In this method, equilibrium is

only satisfied in an integral sense via a modified He11inger-Reissner
variational principle.
1.3.3 Finite Layer
The

finite

layer

method

is

an

efficient

three-dimensional rectangular structures of

analysis

layered

method uses a two-dimensional analytical solution

tool

for

construction. The
in

the plane of a

three-dimensional structure and a one-dimensional finite element shape
function through the thickness. The method has been used to analyze the
structural

and/or

thermal

behaviour

of

rectangular

sandwich

and

isotropic plates with simply-supported edges [39,40]， and vibration of
rectangular plates with various homogeneous boundary conditions [3]. In
addition, the method was employed to study semi-infinite elastic bodies
such

as

the

static

analysis

and

consolidation

of

layered

soil

[41,42], interaction between an elastic slab and its foundation [43,44].
In addition

to rectangular

structures,

the method

can be formulated

within a cylindrical co-ordinate system for studying the vibration of
thick-walled circular cylindrical tubes with various end supports [3]，
and recently, the analysis of pile foundation using an infinite layer
method

[45]•

When formulation

is done within a spherical

system, the method can be applied
most

of

these finite

layer

to 3-D spherical

studies,

structures[46]. In

either harmonic series or beam

vibration functions are adopted as the
structures with supports only around

transverse shear strain across the thickness coordinate z which
contradiction with

the usual

parabolic distribution

assumption

of

is in

parabolic distribution

(

is a commonly adopted assumption for thin and

moderately thick , laminated or isotropic plates). To account for this, a
shear correction factor is required.This factor, which is crucial to an
accurate analysis， is sensitive to the plate gemoetry, configuration,
boundary conditions and

loading.

However,

rigorous

and

theoretically

sound procedures for the calculation of this factor is not currently
available. Also cross-sectional warping is neglected as "plane sections
remain

d i s p l a c e m e n t f u n c t i o n ， and therefore, its accuracy is questionable for
analyzing plates with highly irregular geometry.
In the LT, the displacement variation through the thickness of each
layer

f u n c t i o n s , the COMSFUN can be expressed i n the form
Y (y) = ) [N(y)]{ a }
m
k=l
with

Y (y ) = 1
m

Y

m

m

(y ) = 0
n

m = 1,2,

(2.2)
, r ďź&#x203A; r ^ p+1

(2.3)

17

where

[N(y) ] are the u s u a l

displacement v e c t o r o f

beam shape f u n c t i o n s and { a

the j o i n t s a t

the

} i s the
mk
beam element k

ends o f

corresponding t o the l o n g i t u d i n a l mode Ym â&#x20AC;˘ The summation is carried out
with due regard

to the joint arrangement.

In the following chapters,

unless otherwise stated, we assume r=p+l.
In the preceding discussions, no rotational restraint is applied to
either end of the beam. Put

in other words, the boundary condition at

each end of the beam is either simply supported or free when computing
the corresponding shape functions for
respectively.

As

a

result,

no

the internal nodes or end node

restraint

is

imposed

on

the

first

derivative of the resulting COMSFUN. For 2-D plane stress problems, the
first derivative of the in-plane displacement variables is not involved
in the kinematic boundary conditions, and therefore the afore-mentioned
COMSFUN can be adopted directly as the corresponding longitudinal trial
functions. For thin-plate problems with a transverse clamped or sliding
clamped

edge

however,

the

longitudinal

coordinate)

vanishes

the

at

of

first

derivative

the

out-of-plane

corresponding

end

of

the

(with

respect

displacement

strip.

to

the

variable

Accordingly,

the

rotation at the corresponding end of the beam must be constrained

to

zero. For thin-plate problems with a simply-supported or free edge, no
rotational restraint
Table 2.1 summarizes

where t = plate thickness and varies along y-axis
All integrations are carried out over the entire surface area of
the strip. The explicit form of integrations in both directions has been
discussed in Chapter 2. Note that [Kb3

is a 4x4 matrix.

Load vector
The load vector of a strip with lateral load f(x,y) is given by

26

{F} = { {F}
where

{F} =
m
In

1

A

particular,

transverse

2

f(xďź&#x152;y)
J

,

ďź&#x152;{F}

X

r

}T

(3.6)

[N (x) ] Y (y) dx dy
m

assuming

direction

{F}

m

that

and

the

stepwise

(3.7)

loading

constant

is

constant

longitudinally

in

the

(

the

loading, denoted by f , is a constant within element k of the fictitious
k
beam), then the following load vector can be derived,
{F} = {d/2 d h z
m
where

d/2

-dhz)

1 2 /12
k

{ f f (1 /2
L
k k
k=1

-1 2 /12) { a } }
k
m k
(3.8)

1 /2
k

d = strip width

and {a } = joint displacements of element k of the fictitious beam for
m k

the mth COMSFUN.
Transit ion Strip
As shown

in Fig. 3.3,

transverse direction.

The

there

is a change

longitudinal

shape

in

rigidity across

functions

of

the

the panel

strip are different from those of the column strip. A transition strip
is therefore introduced

to satisfy the compatibility requirement along

the nodal lines between adjacent strips. The displacement function of a
transition strip is given by
r

w(x,y) = V
L
m=l

[( N

1

N o )(Y )
2

1 m

( No N
3

4

) (Y ) 3 { r >

2 m

m

(3.9)

where

(Y ) and (Y ) denote respectively the COMSFUN of nodal lines 1
1 m
2 m
and 2. Derivation of all element matrices follows the same track as for
the original strip.
Nodal point springs
For nodal point springs corresponding to w and dw/dx, its stiffness
can be added directly to their corresponding location along the diagonal
in

the

global

stiffness

matrix.

However,

for

a

point

spring

corresponding to a rotation about X-axis (du/dy), the spring stiffness
is

superimposed

onto

the

stiffness

of

the

strip using

the compound

finite strip method.
Suppose such a rotational point spring is attached to nodal line 1

corresponding strip stiffness matrix and then assembled into the global
stiffness matrix.
Global stiffness matrix equation
The procedure of assembling of the strip stiffness matrices and
load vectors follows that of Cheung[3]. The final matrix equation is
stored

between the two p o i n t loads are subjected t o constant bending moments
( i . e . i n the constant s t r a i n s t a t e o f bending )â&#x20AC;˘ Although they can be
analyzed by

the usual beam elements, a single strip, with different

number of computed shape functions,is used to generate the results given
in Tables

3.1 & 3.2.

As

the

beams are

statically determinate,

the

variation in thickness does not affect the moment distribution. Using
shape functions 2 & 3,which take into account the presence of the two
point loads, accurate displacements, rotations and bending moments can
be obtained for the prismatic and

the non-prismatic beams. Using the

first shape function alone, however, the mid-span moments in both cases
differ from the exact values by more than 30 percent. As expected, very
accurate mid-span deflection can still be obtained by using the first
shape function alone, even though the latter does not take into account
the presence of the point loads.
modes
y=l

load ( q=l ) on the top surface. Because of symmetry, half of the beam
is

discretized

into

eight

finite

strips

through

the

depth

with

increasing number of nodes along the span. The boundary conditions at
the clamped ends ( u, v=0 at y=0) and mid-span ( v=0 at y=l) can be easily
imposed in the usual finite element manner. Results of deflections are
summarized

in

Table

5,1

and

solutions using 8 strips and
both methods are

compared

with

the

spline

finite

strip

11 sections[13]. Results obtained using

in very good

agreement. Convergence of

stresses are

given in Fig.5.4 and compared with the spline finite strip solutions and
the finite element

solutions.

At

mid-span of

the

beam, 5 nodes are

sufficient to give very accurate stresses. However,at the top and bottom
surfaces

of

the

clamped

end

where

the

stresses

are

theoretically

infinite, 9 nodes are required to pick up the sharp changes in stresses.

Sliding base : u=0 along x=2.5
Shear wall supported on portal frame [Fig.5. 12)
The present finite strip is also applied to the analysis of a shear
wall with a portal frame connected

to

its base. The portal frame is

modelled with the usual beam elements. To ensure full compatibility of
rotations between beam elements of the frame and finite strips of the
shear wall, the latter is divided into finite strips with the following
modified displacement functions,

78

where

u(x,y)

[M(x)] Ym(y) { a }

(5.13a〕

v(x，y)

[N(x)] Ym(y) {

〔5.13b)

/
3 }

[N(x)] denotes the beam shape functions. Accordingly, the nodal

displacements are modified as
{oc}

= { u., , u n }

and {p}

= { v

dv/dx

m
1
2
m
m
1
1
2
2 m
where dv/dx. represents the nodal rotation of vertical fibres of the
shear wall. Because of the addition of in-plane rotation at each node,
the stiffness matrix of each of the beam elements of the portal frame
can be

directly

manner,

onto

attaches.

superimposed,

the

(In

in

the

stiffness matrix

the

usual

of

semi-analytical

the

finite

standard

finite

finite

strip

strip

element

to which

method,

the

it

beam

stiffness matrices are first transformed and subsequently added onto the
corresponding

finite

strip

stiffness

matrices).

Formulation

of

the

stiffness matrices follows the same track as before.
The

accuracy

of

assessed by comparing
analytical

finite

Experimental

the

finite

strip

in

this

example

is

its solutions with those obtained by the semi-

strip

results

present
method

are

also

and

the

adopted

for

finite

element

comparison.

In

method.
the

semi-

analytical finite strip method, the COMSFUN in (5.1) is replaced by the
analytical

series

of

vibration

of

a

vertical

cantilever

with

an

elastically-supported base. Furthermore, the finite element results are
computed using a 4-node beam-type membrane element with nodal in-plane
rotations. Details of the semi-analytcial finite strip and the beam-type
finite element can be found in [132].
Vertical deflections

and

stresses

are

depicted

in

Fig.5.13.

As

sharp variation of stresses are expected in the proximity of the wallframe interface, comparison is made at a small distance away from the
centre-line
using 5

of

the

transfer beam.

The present finite

longitudinal nodes with non-uniform spacing

compared with

the five-terms solutions of

strip results,

(Fig.5.12c), are

the semi-analytical finite

strips. Together with the finite element and experimental results, four
sets

A deep cantilever with step-change in thickness at
y=10. Thickness t = 0.1 for y^lO and t = 0.2 for
y>:10. (a) FS1 : A mesh of 10 strips with 9 nodes
per nodal line is used, (b) FS2 : A mesh of 20
strips with 5 nodes per nodal line.

The free vibration and buckling under uniform compressive stresses
of continuous beams are considered. Although the beams can be analyzed
by the usual beam elements, they are treated as folded plates in this
example. The cross-section of

the beam, which was modelled with five

finite strips, is an I-section made up of two equal flanges of width 0.5
and web of 1.0 and a constant thickness of 0.1 throughout. In contrast
to

the

cantilever

folded-plates

in

previous

examples,

significant

membrane action is expected in the component flat plates of the beams.
The boundary conditions at

the end

supports and at

the

intermediate

support are modelled as rigid in-plane diaphragm. These conditions can
be easily imposed, in the usual finite element manner, by constraining
the in-plane degrees of freedom at the corresponding locations in the
global

stiffness

conditions
splines
priori

require

while
for

matrix

(

tediuos

In

modification

continuous-beam

the

spline

vibration

semi-analytical

finite

finite
of

strip

the

method,

these

corresponding

local

functions
strip

are

determined

method.)

Results

a
are

compared with beam theory[76] in Table 6.3. Agreement between the two
sets of results is found

to be satisfactory, as only a maximum of 2

percent difference is found.
It is worth mentioning that in cases where the internal diaphragms
at the supports are flexible, the in-plane stiffness can be taken into
account by dividing the diaphragms into the usual 2-D finite elements or
boundary elements. After applying condensation, the stiffness matrix of
the

diaphragms,

corresponding

to

the

degrees

of

freedom

of

the

connecting nodes between the diaphragms and the folded-plate structure,
can be directly superimposed onto the global stiffness matrix. Unlike
semi-analytical

semi-analytical finite strip (with complete decoupling of the stiffness
matrix corresponding to each term of the trigonometric series),

it is

used to examine the accuracy of the present finite strip for computing
higher

frequencies

of

prismatic

thin-walled

structures.

The

present

results of the lowest twenty frequencies are given in Table 6.5 together
with

the

' exact'

finite

strip

solutions

due

to

Wittrick

and

Williams[80]. 3-D view of the first seven mode shapes are also given in
Fig.6.7 from which one can see the coupling of the transverse bending
and longitudinal bending. Despite the fact that a relatively coarse mesh
with twelve strips and 6 nodes is intentionally used

to predict the

deformation patterns of the complicated mode shapes, reasonable results
can still be obtained. Except for modes 17,19 and 2 0 ， t h e differences
between the two sets of results are less than 6 percent. A maximum of 10
percent difference is found at mode 20.

identical and correspond to a typical high-modulus fibre composite with

E

1

/

= 30

2

G

/ E,= 0.6
12
2

u

12

0.25

They are simply-supported at both ends with rigid in-plane diaphragms.
Each of these tubes has a wall breadth B ( measured between the centre
lines of a pair

of opposite walls), wall

thickness 0.01B,

and

span

length 10B for the single-cell tube and 5B for the twin-cell tube. As
the wall thickness is relatively small with respect to the span of the
tubes, significant out-of-plane bending of the wall is anticipated. In
the finite strip analysisďź&#x152;each tube wall

is modelled by

two finite

strips. In Table 6.6 the finite strip results, using 8 strips x 9 nodes
for the single-cell tube and 14 strips x 9 nodes for the twin-cell-tube,
are compared with the finite element results using 80 four-node shell
elements for

the single-cell tube and

70 elements for

the

twin-cell

tube. For the single-cell tube, the lowest eight frequencies obtained by
the finite strip method, with 55 percent of the number of unknowns of
the finite element analysis, are within 1 percent difference from the
finite

provides a good test case for evaluating the efficiency and accuracy of
the present finite strip.
In

the

finite element analysis, a very fine mesh of four-node

assumed-stress hybrid element was used to provide a benchmark solution.
The panel is divided into two elements along the depth of the stiffener,
four elements between stiffeners,
1296 elements,

36 elements along

length,

total

unknowns.

In the spline finite strip analysis, six spline sections is

longitudinally with a

method,

the stiffened

panel

total

1369 nodes and

the

making a
used

of

and

of

approximately 8000

504 unknowns.

is divided

into

For

18 strips

the

present

( two strips

between stiffeners and one strip for each stiffener) and seven nodes
longitudinally. There are a total of 436 unknowns after imposing the
boundary conditions. A comparison of the three sets of results is given
in Table 6.7. The results obtained by the present finite strip model,
with only 5.5% of the number of unknowns of

CHAPTER 7
A THIRD-ORDER SPLINE FINITE STRIP
AND
NONLINEAR FORCED VIBRATION
I n t h i s c h a p t e r the o r i g i n a l s p l i n e f i n i t e s t r i p , which i s based
upon the t h i n - p l a t e t h e o r y , i s m o d i f i e d i n such a way t h a t the present
spline

finite

strip

is

formulated w i t h i n

third-order plate theory.
finite strip

is first

The performance of

the

nonlinear

analysis

procedures. For reducing
analysis, a special
number

of

equations

context

of

Reddyďź&#x152;s

this third-order spline

investigated with particular attention on thin

plates as well as shear-deformable plates.
to

the

of

forced

It is subsequently extended

vibration

using

the computational cost of

reduction method
is

solved

is employed

within

each

time-stepping

the
so

time

time-stepping

that a smaller

step

and

hence

a

reduction of computing time can be achieved.
7.1 Introduction
The

transverse

shear

moduli

of

usually very

low compared with the

result

transverse

that

importance compared
the usual assumption

shear

modern

composite

materials

Inplane tensile moduli, with the

deformation

can

be

of

to homogeneous isotropic materials.
of

are

considerable
Additionally,

"plane section remains plane"

is no

longer

appropriate for thick isotropic plates and moderately thick laminated
composite plates because of the significant cross-sectional warping. The
classical Kirchhoff plate theory is only capable of dealing with plates
with very large span-to-thickness ratio. Although the shear deformation
theory developed by Mindlin and Reissner has been commonly used for the
finite

element

warping
gives

analysis

of

moderately

thick

plates,

cross-sectional

is not taken into account. Furthermore, Mindlin plate theory
constant

transverse

shear

strain

distributions

through

the

thickness and hence a shear correction factor is required to account for
the actual parabolic distribution.( For moderately thick laminated or
isotropic platesďź&#x152; transverse shear strains are commonly assumed as a
parabolic distribution through thickness). This factor depends on the
geometry

of plates

conditions and

as well

loading.

as

its

laminate configuration,

It plays a paramount role

boundary

in the analysis of

laminated composite plates. However, accurate and reliable analysis of

By selecting different terms from these expressions, a very large
number of third-order plate theories have been developed(see e.g.[84]
and references therein). To reduce the number of independent variables,
one can

impose the the

traction boundary conditions on the top and

bottom surfaces of a plate and neglect transverse normal deformation. A
third-order plate theory, which contains the same number of dependent
variables as the Mindlin plate theory, has been developed by Reddy[85],
Details

of Reddy's

third-order

plate

theory

and

its

implementation

will be discussed in the following section.
7.2 Third-order plate theories
By neglecting transverse normal deformation, the displacement field
of a plate starts off with
U(x, y, z) = u(x,y)+ 2 0

X

V(x,y, z) = v(x,y) + z ^

y

W(x’y,z) = w(x,y)

+ z2co + z 3 5
X
X
+ z2o) + z 3 6
y

y

(7,2a)
(7.2b)
(7.2c)

where u,v and w denote the displacements of a point (x,y) on the midplane, and

b
l

x

and

y

iIj

are the rotations of normals to mid-plane about

the Y- and X-axes, respectively. Cross-sectional warping is represented
by a)
x

to)

y

,6
x

,5
y

.

First consider the transverse shear strains which are given by
(7.3a)

115

C

xz

= U

z

+ W

(7.3b)

,x

where p a r t i a l d i f f e r e n t i a t i o n i s denoted by a comma such as ()

ďź&#x152;

=~ 2

The t r a c t i o n boundary c o n d i t i o n s on the top and bottom surfaces o f
a p l a t e w i t h t h i c k n e s s t can be expressed as
cr

material, the four traction-free conditions of (7.4) can be applied to
eliminate four out of the nine variables in (7.2). Hence the following
displacement field of the Reddy's third-order plate theory (TOP) can be
derived :

U(x,y,z)

ip - 4/3 ( z/t ) 2 (

V(x,y,z) = v +

4/3 ( z/t ) 2 ( 0 + w

)]

(7.5)

W(x,y,z) = w(x,y)

It is clear that (7.5) contains the same number of independent variables
as the Mindlin plate theory.
By substituting (7.5) into (7.3), it can be easily shown that the
transverse shear strains at mid-plane are actually represented by

corresponding unknown parameters.
It has been shown by many researchers that the transverse shearing
strains

must

be

interpolated

consistently

in

order

to

avoid

shear

locking for finite element analysis of thin plates. In the spline finite
strip model based on TOP,it is obvious that, by substituting (7.10) into
(7.6), the contributions to the transverse shear strains at mid-plane
are

of different

order polynomials.

This causes a problem

for

thin

plates where both G
and e
must tend towards zero. A consistent
xz
yz
interpolation scheme would be one in which identical polynomial terms
must be used for both the rotations
the

vertical

interpolation

displacement
scheme

is

field

possible,

xjj (or xb )and the derivative of
x
y
w
( or w
). While such an
.x

it

,y

would

not

be

convenient

to

implement. However, with the modified displacement field (MODTOP) this
problem does not exist because the transverse shear strains e

xz

or

and

and c

yz

(

can be interpolated consistently by replacing 〜 a n d 〜

of equation (7.10) by A

and A

respectively.

7.4 Equilibrium Equation
Taking

into

consideration

the

von-Karman nonlinear

strains

and

initial imperfection, the strain-displacement expressions can be written
as

118

w

V w
,y 0,)

(7.11)

where w。 denotes the out-of-plane initial imperfection.
In a condensed form, the above relationships can be written as the
summation of linear strains, nonlinear strains and initial strains ； that

convergence of centre bending moments, obtained by the MODTOP model, are
given in Table 7.1. Again, very rapid convergence can be observed.
Results

of

analyses

of

moderately

thick

isotropic

plates

and

laminated composite plates with different span-to-thickness ratio are
presented in Tables 7.2 and 7.3.

In all cases, the plates are subjected

to uniformly distributed loads. A quadrant of each isotropic plate is
divided into strips with 4 spline sections. Results given in Table 7.2
indicate the good agreement between the spline finite strip results and
the

thick-plate solutions by Kant

[89] for

isotropic plates. Unlike

thin-plate problems, both models show rapid convergence
For

composite

laminates,

Reddy[85] are used

analytical

for comparison.

solutions
It

of

in this case.

TOP

is apparent

that

developed
there

by

is no

significant difference among the three sets of results given in Table
7.3.

It is important to see that there

is no major difference between

the results of the two spline finite strip models in analyzing sheardeformable plates.

Two isotropic square plates with different thickness ( a/h=5 and
100) are subjected to uniform distributed loads. All edges of the plates
are either simply-supported or clamped. Only a quadrant of the plate is
analysed with four finite strips and

four spline sections.Information

for the thin plate with a/t=100 is provided in Fig. 7.3. For presentation
of results in Fig.7.4 and Fig.7.5 the following

normalised quantities

are used ：
Q = q a 4 / E t4

and

From the results it

W = w(0,0) / t

is clear that the present finite strip results are

in reasonable agreement with the Levy* s solution for thin plates [91]
and

During the initial stage of edge compression, the loading is mainly
carried b y the bending resistance of the plate. As the MODTOP model has
been

shown

to

perform

better

in

especially with clamped edges， one
that

a

finer

mesh

is

required

for

bending

problems

of

thin

plates,

can easily observe from Fig.7.7(b)
the

TOP

model

to

achieve

similar

accuracy as the MODTOP model during

the initial stage of loading. When

the

the

edge

compression

resistance

becomes

transverse

shears

increases,

significant.

affects

only

importance

As
the

consistent

bending

of

the

membrane

interpolation

behaviour

of

the

of

strip

models, the difference between results of the two models reduces as the
edge compression increases.
o n MODTOP

are

in

In general,

reasonable

agreement

the finite strip results based
with

those

provided

b y Yamaki

[93].

Preliminary

Remarks

The third-order spline finite strip is applied to a few linear and
nonlinear static problems. For the analysis of shear-deformable plates,
there is n o major difference between the performance of the finite strip
models based o n TOP and MODTOP. However, the MODTOP model
in

the efficiency of the algorithm for transforming the complete system to
the reduced system. In this work, the R i t z vectors are adopted because
they have been shown to b e accurate, efficient and relatively easy to
generate when compared with vibration mode shapes for
vibration

Fig.7.16 Response of
the
square plate subjected
to a
exponential load of maximum amplitude SOq^. t=0.001 sec per time step.
(a)
vertical
displacement
response
response•Legend same as Fig.7.9

(b)

inplane

displacement

•Integration
reQion
(b)

Fig.7.17

(a) A t y p i c a l Ba-spline f unc t i o ns (b) a t y p i c a l
Integration
region
and
the
splines
that
intersect i t .

‘‘

CNJ
-M

c

0)

E0)

-f-*

c
①

•

TT"

CO

C
①

c
①

m

CO

£
①
l O

E0) E0)

C(2)-continuous
Fig.7.18

Discretization of a multi-connected domain into
spline finite strips with modified spline
function. Each finite strip can be treated as a
finite
element
for
assembly.
Displacement
unknowns at node k of element i are
(u) ，
rI
cI

{ q } = nodal displacements of the element k computed from
k j
the j-th term of the global function Q(x, 2).
ne = number of elements

(y) represents the i-th term of the spanwise Bs-spline function with
,

and

.

being

the

corresponding

unknown

displacement

parameters. Kinematic boundary conditions at y=〇 and y=a are taken into
account by modifying the corresponding local splines.
The global polynomial

functions

are expressed as

the Kronecker

product of two series,
Q ( x , 2 ) = X(^)-ZCz)
where

X(>：) = {

and

Z(2)

=

, (1+^)/2

{ 1 j2 ， 2 ， . . . ’ 2

(8,2)
, (/-¾)

2 € [0,h]

c)

162

with

x = 2 x / b - 1

义€

c = 1

if r is even

c = ^

if r Is odd

[-1,1]

x e [0,b]

q = (r+1)(s+1)
= t o t a l number of terms in the global function
The series in X-direction is so chosen that the kinematic boundary
conditions at the two ends of X-axis can be easily treated by removing
appropriate unknown displacement parameters corresponding

to the first

two terms in the series. Note that no kinematic boundary conditions are
imposed on the top and bottom surfaces of the plate and a simple power
series

in

Z-direction

suffices.

By

varying

the

parameters

r and

s,

different orders of polynomials can be used in the X- and 2-directions.
As thick plates are treated as 3-D elastic bodies

Comparisons of the results with the exact solutions are shown in Table
8.2. A good agreement between the two sets of results can be observed;
all stresses are within 5 percent difference from the exact solutions.
An 82 percent

reduction

in

the number of unknowns within

the cross-

section is achieved in this analysis.

士丄）"U'土;）Txy(0,0,±;) TxzCo, ? # 0 ) w(?f ? # 0 )

Exact

0.750

0.482

1. 171

-0.0479

-0.705

-0.504

-0.166

0.0473

0.767

0.494

1.010

-0.0490

-0.723

-0.519

-0.005

0.0485

0.787

0.502

1.016

-0.0503

-0.743

-0.527

-0.015

0.0498

0.801

0.534

-0.755

-0.556

1.0

.210

1.896

•2 4 3

1.926

• 270

-0.0511

.961

2.02

.256

0.0505

Table 8.1 Comparison of results for a three-ply laminate
with span-to-thickness ratio of 4.
*

compared with a three-dimensional f i n i t e element s o l u t i o n due to Barker
[119] . The

same g l o b a l polynomial with orders r=4 and s=7 i s used i n

t h i s example.

In order to obtain an accurate transverse shear stresses

d i s t r i b u t i o n , a f i n e r mesh of 9x9 elements, with each p l y layer divided
i n t o three element l a y e r s , i s used to d i s c r e t i s e the cross-section. By
r e f e r r i n g to F i g . 8 . 3 ， i t is apparent that the agreement between the two
sets

of

results

is reasonably good.

In

this example,

the

number

of

unknowns within the cross-section is reduced from 840 to 96 compared to
a full 8-node spline finite prism analysis.
A Laminated composite plate with five layers
A five-ply laminate( 5

=4 )

with ply angles 0°/ 90°/ 0 / 90 /

0 0 i s analyzed. Unlike the previous examples, the total thickness of the
0° and 90° layers are the same, whereas layers at the same orientation
have equal thicknesses. Under these conditions, the effective laminate
stiffness in the X- and Y- directions are the same. The cross-section is
discretised into a 10x10 mesh with each layer of materials divided into
two element layers. Using a global polynomial of orders r=4 and s=7 in

C, is the layerwise coordinate ranging from -1 to 1 for z€[z
z 1
i-r i
It is computationally inefficient to analyze the entire plate with

this layerwise model because of the large number of unknowns involved.
To

reduce

the

computational

cost,

a

special

form

of

constraint

is

introduced. In the area where only the global response are concerned, a
set of through-thickness global functions is selected

to transform the

nodal variables of the layerwise cubic shape functions to a small set of
generalized parameters associated with the global functions. In essence,
the layerwise cubic shape function is enforced to follow a small number
of

pre-defined

patterns.

With

a

limited

sacrifice

number of unknowns is thereby reduced with respect

of

accuracy,

the

to a full layerwise

model. The j-th term of the global-local interpolation can be written as

Table 9.1 ďź&#x161; Interpolation across thickness in different regions
p = number of terms used in equation (9.1)
nl = number of layers of material
r = number of terms taken from the global functions
given in expression (9.3) and r ^ 3 nl + 1.

Having

defined

the

displacement

functions

of

the

plate,

the

strain-displacement relation can be written down as

N
N

Z
N Zf

N

Z

N t

xz^

and

(9.4)

Z
N 2'

or

w:

{c} = Y [B].{6}.

the

element

takes the form,

stiffness matrix

for

each

of

the

three

regions

180
[K]

[K]

[K]

[K]

[K]

(9.5)

[K]
with [K]

[b]t[d][b]

tK]

i

dV

[K^]

[K

] [K

J

[K

]

[K

] [K

]

]

[K

] [K

]

11

[K

21
31

p，p

12

22

32

13

23

33

(9.6)

and
D

11

2 2

D

[Ki J
12 i

D

D

13•

[K.

D

[K J

D

21

22 r

[K

[K

[K

D

23"

D

31

32

]

mn

D

D

33 '

Corresponding
above

Z

13

2

m

n

dz

N

dz

N

Z

m

n

n

m

m n

n

dz

dz

2' 2
m

n

m

Z

n

M

dx dy

D

N

dx dy

—

dx dy +

dx dy +

N

dx dy
+
J

dz

N

N

through-thickness

expressions

+

M N

m n

Z

dz

Z

Z

44 m

N

D

dz

n

2' 2

6 6 m

N 7 N dx dy +
»y

Z'Z^ dz

66

一

dx dy + D

D

N' M

dz

Z Z

44 m n

N

dx dy

N

N
dx dy
,y , x

44

44

Z

Z

m

m

N1 N

dz

n

Z

2

n

n

dx dy

,X

dz

N1 M

dz

NT N

dz

N N

»X

, y

dx dy

dx dy

iVT

n

r

m

D

N dx dy +

,x

Z’ Z* dz

2’Z

33

2

N

- - - -

dz

Z

32

V
i

Z

m

dx dy +

, X

NTN dx dy

N

55

23

N

X

Z 21 dz

Z

22

31

n

m n

12

D

N

2’ 2’ dz

66 m

12

D

[K

dz

11 m n

according

to

D
D

D

2' 2

55

m

n

D

55

Z

m

r

，

Z

m

dz

n

Z

n

dx dy

Nl N dx dy

Z Z’ dz

66 m n

55

,y

dz

NT N
»y
NT N

dx dy

dx dy

interpolation
the

region

is

substituted

concerned.

The

into

the

in-plane

181
eight-node shape f u n c t i o n s , denoted by
integrated

using

2x2

Gaussian

points.

In

the

through-thickness

d i r e c t i o n , however, i n t e g r a t i o n can be done a n a l y t i c a l l y f o r each of the
three

digit denotes the number of elements in each of the in-plane direction.
The third digit denotes the number of 1-D cubic element in the thickness
coordinate.

Values

transverse stress.

in

bracket

are

the

maximum

of

the

corresponding

185

Mesh

crx

669

1.274

1.063

0.206

0. 199

-0.0730

-0.876

-0.834

(0.225)

(0.210)

0.0541

1.291

1.078

0.209

0.196

-0.0740

-0.888

-0.846

(0.228)

(0.211)

0.0548

1.368

1.154

0.224

0.213

-0.0795

—0.951

-0.905

(0.245)

(0.225)

0.0587

1.260

1.051

0.204

0.194

-0.0722

-0.866

一 0 . 8 2 4

(0.224)

(0.211)

0.0534

0.692

0.635

0.226

-0.0341

»0.656

-0.619

(0.228)

0.0332

0.701

0.644

0.229

-0,0346

-0.666

-0.627

(0.231)

0.0336

0.684

0.628

0.223

-0,0337

-0.649

-0.612

(0.225)

0.0328

449

229

exact

669

449

exact

#

？

)

cry (；

xxz (0# , 0 )

( 言 ， 0 , 0 ) TxyCoj

.226

.229

.223

w(；

.0)

12.288
12.288
12.280

12.288

.079

.079

.079

Table 9.3. Comparison of results for a nine-ply laminate with s = 2
and 4. Only Layerwise cubic shape function is used.(*) Mesh : first two
digit denotes the number of elements in each of the in-plane direction.
The third digit denotes the number of 1-D cubic element in the thickness
coordinate.

Values

in

bracket

are

the

maximum

of

the

corresponding

transverse stress.

Assessment of global-local

/wodeJ (GLM)

To examine the accuracy of the global-local model, the plates are
first

discretized

into

element per material
from

the

a

mesh

layer.

through-thickness

of

4x4

B-nodes

elements

Increasing number of
global

functions

and

terms
used

and

a

cubic

is then taken
throughout

the

entire plate. The centre deflections of the resulting global-local model
are compared with the solutions of the corresponding layerwise model and
the

exact

solutions,

see

Table

9.4.

Considering

the

results of

the

four-ply laminates obtained by taking 5 terms from the global functions,
one

For each span in either direction of the X-Y plane, the modified beam
vibration functions are numerically integrated using 10 Gaussian points.
Integration of the 1-D linear shape functions has been discussed in the
Appendix of Chapter 5.

functions, w e will first consider the following continuous beams with
various geometric configuration.
(I) A three-span continuous beam with uniform support spacing and
different end supports. Bending rigidity and mass density are constant
along the overall span.
(II) A four-span continuous beam with non-uniform support spacing
and

simply-supported

ends.

Bending

rigidity

and

mass

density

are

constant along the overall span.
(III) A two-span continuous beam with non-uniform support spacing
and simply-supported ends. The beam is non-prismatic in the sense that
it is made up of two prismatic spans with different bending rigidity and
mass density.
In

and mass p e r u n i t l e n g t h a l o n g the span.
In a d d i t i o n to the i n c r e a s i n g number o f

terms used i n each case,

the s t a r t i n g term o f the t r i a l f u n c t i o n s a l s o varies； for a beam with q
intermediate supports, and if the lowest nf frequecies are sought, then
starting term r = 1 or 1+q
nf ^ number of terms taken, nt ^ 2 x nf
For

comparison

with

analytical

solutions

[125],

the

following

normalized frequency (X) is adopted,

o) = (A/L)

2

(EI/p)

1 /2

where L = beam overall length

The present results are summarized in Tables 10.3 to 10.5.

In all

cases, except for the second non—prismatic two-span beam in Table 10.5,
2 x nf terms are sufficient to give an accuracy of less than 1 percent
difference between the present results and the analytical solutions. For
the

second

non-prismatic

two-span

beam,

a

maximum

of

difference is found at mode 4 when 2xnf terms are used

2.7

percent

to compute the

frequencies. This relatively larger difference can be attributed to the
fact

that

the

mass

density

influence
of

the

of

the sharp changes of bending rigidity and

beam

are

not

taken

into

account

in

the

trial

functions•
In general
rapidly
trial

to

one

can

observe

the analytical

function

is.

When

that

the

present

solution whatever
the

frequencies

the

are

solution

starting

computed

converges

term of the

using

only nf

terms, however, the higher frequencies in most cases are more accurately
predicted b y starting
expected

as

the

appropriately

conditions;

higher

the

Unfortunately,

mode

faster

for

the

the

simply-supported-free

trial

modes

function from

of

shapes

a

of

convergence

cases of

the

term

single-span beam
the

corresponding

is not

guaranteed

1+q.

This

is

approximate more
multi-span
for

all

beam.

boundary

three-span beams with free-free ends or

ends， unreasonable

results,

which

are

excluded

from Table 10.3, are obtained by starting the trial functions from the
term

square laminates are computed using one term of the t r i a l functions i n
each of the two i n - p l a n e d i r e c t i o n s ( r=s=m=n=l ). They are compared with
the three-dimensional f i n i t e d i f f e r e n c e s o l u t i o n of Noor[120]. The f i b r e
o r i e n t a t i o n s o f the d i f f e r e n t laminas a l t e r n a t e between 0Oand 90 O with
respect t o the X - a x i s , and i n the symmetric laminates the 0°layers are
the outer surfaces. The total thickness of the aligned layers are the
same，whereas the individual layers are taken to be of equal thickness.
Geometric and material constants are given as follows ：
span—thickness ratio
E

intermediate l i n e supports are shown i n F i g . 10.2. The p l a t e s are made of
isotropic

material

or

laminated

composite

material.

In

case

of

laminates, the m a t e r i a l constants used are the same as those given i n
the preceding s e c t i o n except that

=

i s replaced by

=G

in

the f o l l o w i n g examples. In a l l cases, the r a t i o between the shorter span
and the t h i c k n e s s i s taken e i t h e r as 5 or 100. Support conditions along
the s i d e s of the p l a t e are i n d i c a t e d by four l e t t e r s that define the
type

Fig. 10.2(b), the locations of the intermediate line supports on X- and
Y- axis are replaced by (0.3b,0.65b) and (0.45b.0.85b) respectively. For
isotropic plates, the natural frequencies are normalized as
A2 = u)2 phb4 / D

where D = Eh3/ 12(1-1^2) and v-0.3

while for composite plates,
X

2

= a)2 phb4 / ( D

11

D

22

) 1/2 where D. = 4E h/45
ii

xi

The starting terms for the present example are taken as r^l+q^ and
s=l+q ,
y
In order

to assess

the accuracy

of

the modified

â&#x20AC;˘
beam functions

relative to other methods available in the literatures, the lowest six
frequencies of plates with the shorter span-to-thickness ratio of 100
(although finite layer method is more appropriate for analyzing plates
with small span-to-thickness ratio) are computed and summarized in Table
10.7. Using five finite layers and five terms of trial functions in each
of the two in-plane directions, the present results for the 2 way-2 span
plate differ by a maximum of 4 percent from those due to Leissa[126].
Results

and 10.10. These sets of r e s u l t s are presented f o r the f i r s t time i n
open

literature,

and

no

comparison with

other

methods

are

made.

Increasing number of terms was used to obtain a converged solution i n
a l l cases. For square p l a t e s of 2 way-2 span and 2 way-3 span, equal
number of terms of t r i a l f u n c t i o n s are used i n each of the two in-plane
d i r e c t i o n s . For 1-way continuous p l a t e s , however, more terms are taken
i n the continuous-span d i r e c t i o n .
Based on the r e s u l t s shown i n Tables 10.8-10.10, a few remarks can
be made.
( i ) By comparing the frequencies of Table 10.8 w i t h the those i n Table
10.7, one can observe that, because of the increased f l e x i b i l i t y due to
transverse

shear

isotropic plates

and

normal

deformation,

the

frequencies of

thick

( Table 10.8) are lower than those o f t h i n plates

(Table 10.7). I t i s apparent that the amount of d i f f e r e n c e depends on
the types of boundary conditions and the arrangement of intermediate
l i n e support.
( i i ) Results i n Table 10.9 i n d i c a t e that the frequencies of the 3-layers
p l a t e s are higher than those of the 5 - l a y e r s p l a t e s despite what types
of boundary conditions are concerned.
Table

10.10

The

same trend can be seen i n

by comparing the r e s u l t s of

the 2 - l a y e r s and 4-layers

p l a t e s . I t i s a l s o noteworthy that the same amount o f material i s used
i n the 2 - l a y e r s and 4ä¸&#x20AC;layers p l a t e s or the 3 - l a y e r s and 5 - l a y e r s plates.
This implies that the frequencies of a p l a t e w i l l increase by d i v i d i n g
the material i n t o more layers.
( i i i ) I t i s well-known that the bending-stretching coupling of skewsymmetric laminates w i l l , i n general,
buckling

adopted as the spanwise trial, functions for a new flat-shell finite
strip. This finite strip has been successfully applied to the flexural
vibration and buckling of plates with complicated support conditions,
bending of rectangular plates with changing longitudinal ridigity and
vibration and buckling problem of prismatic folded plates. In addition,
2-D plane problems with non-rectangular domains can be treated

in a

standard finite element manner. The present finite strip retains the
advantages of the classical finite strip and the spline finite strip
method such as rapid convergence for point-loaded plates, reduced number
of nodal variables when compared with the finite element method and ease
of modelling. Additionally, the present finite strip method possesses a
few

superior features when compared with

the semi-analytical finite

strip and the spline finite strip such as the capability of dealing with
geometric

discontinuities

problems,

its

simplicity

conditions, and

in
in

plate-bending
handling

problems

boundary

and

and

2-D

Internal

plane
support

its ease of implementation. Furthermore, only simple

shape functions are involved in the calculations and computation of the
stiffness matrix involves no numerical Integration.
The original spline finite strip is restricted to the analysis of
thin

plates.

In

this

work,

a

third-order

spline

finite

strip

is

developed as the counterpart of the original spline finite strip which
is capable of treating moderately thick isotropic or laminated composite
plates. The formulation is done within the framework of a third-order
plate

theory.

Its

accuracy

in

performing

linear/nonlinear

plates with different span-to-thickness ratios
standard examples. For nonlinear vibration,

analysis

is assessed by a few

a. special reduction incthod

225

i s proposed f o r reducing the computational cost of the time-stepping
a n a l y s i s . The success of the reduction method r e l i e s on the b e n e f i c i a l
features of the R i t z vectors, a
pre-transformation technique.
A new s p l i n e f i n i t e prism method i s developed and a p p l i e d to the
a n a l y s i s of rectangular, t h i c k , laminated composite p l a t e s . No kinematic
assumptions are imposed and 3-D e l a s t i c i t y theory i s adopted. The method
i s formulated using B3-spline functions i n the l o n g i t u d i n a l d i r e c t i o n
and a set of g l o b a l - l o c a l polynomials i n the c r o s s - s e c t i o n o f
plates.

The

rectangular

method i s

validated

by

means

of

a

few

examples

the
of

laminated plates w i t h a v a i l a b l e a n a l y t i c a l s o l u t i o n s or

The computed shape functions can be s l i g h t l y modified and be used
i n each of the two orthogonal d i r e c t i o n s of a rectangular t h i n - p l a t e
element(Fig.11.1). The r e s u l t i n g element can be assembled i n standard
finite

element

manner

and

C^-continuity i s

maintained

across

the

boundaries between adjacent elements. S p e c i f i e d displacement constraints
can be imposed on the corresponding nodal displacement v a r i a b l e s . The
same idea can be applied to

2~D e l a s t i c problem and a higher-order

f l a t - s h e l l element can be e a s i l y f i g u r e d out.

Redaction

method for nonlinear

vibration

The high computational cost o f

a nonlinear v i b r a t i o n a n a l y s i s

stimulates the development of an e f f i c i e n t reduction method. Only a
preliminary investigation

circular thin-plate finite strip. This finite strip can be applied
the

analysis

of

circular

bridge

slabs

or

other

sector

to

plates with

continuous line supports in the radial directions (see Fig.11.2).
Adaptive finite strip methods
Error estimation for plate bending problems can be developed for
finite strips

in such a way that mesh refinement can be automatically

carried out in the area of rapid stress gradient.
Combination of finite strip and finite prism methods
In a previous study conducted by Chong et.ai-[40], finite strips
are combined with finite prisms for

the analysis of sandwich panels.

This approach can be extended to the calculation of transverse stresses
in a laminated
material

layers

plate;
where

finite prisms
transverse

are used

stresses

to model each

are

required

of

while

the
the

remaining layers are smeared Into finite strips. When compared with a
complete finite