Random vibrations of nonlinear elastic systems


Material Information

Random vibrations of nonlinear elastic systems
Physical Description:
vii, 85 leaves : ill. ; 28 cm.
Herbert, Richard Edgar, 1938-
Publication Date:


Subjects / Keywords:
Vibration   ( lcsh )
Elasticity   ( lcsh )
Nonlinear theories   ( lcsh )
bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )


Thesis--University of Florida.
Bibliography: leaves 82-84.
Statement of Responsibility:
By Richard Edgar Herbert.
General Note:
Manuscript copy.
General Note:

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 000568325
notis - ACZ5050
oclc - 13645159
System ID:

Full Text






April, 1964


The author wishes to express his sincere gratitude to Dr. William

A. Nash, Chairman, Department of Engineering Science and Mechanics,

for serving as chairman of his supervisory committee and for his

constant advice and encouragement throughout the author's entire

graduate studies program. He would also like to thank Dr. T. S. George,

Professor of Electrical Engineering, Dr. A. Jahanshahi, Assistant

Professor of Engineering Science and Mechanics, Dr. I. Ebcioglu,

Assistant Professor of Engineering Science and Mechanics, and Dr. R.

G. Blake, Associate Professor of Mathematics, for serving on his

supervisory committee and for the various stimulating discussions he

has held with them over the past few years.

Final thanks go to the Air Force Office of Scientific Research for

their sponsorship of this program.








2.1. Analysis of Deformation . .

2.2. Equations of Motion . .

2.3. Boundary Conditions . .


3.1. Stochastic Processes and Probability Theory .

3.2. Response of Linear Systems . .


4.1. Classification of Random Processes .

4.2. The Fokker-Planck Equation . .


5.1. General Theory . . .

5.2. Some Special Cases ... .

a. Simply Supported Beam .

b. Simply Supported Plate . .






















. . .




6.1. Simply Supported Beam

6.2. Simply Supported Plate




. . .

. .

. . .

. . .

. . .

. . .


Figure Page



NONLINEARITIES . .... ... .. 72



Abstract of Dissertation Presented to the Graduate Council
in Partial Fulfillment of the Requirements for the
Degree of Doctor of Philosophy



Richard Edgar Herbert

April, 1964

Chairman: Dr. William A. Nash
Major Department: Engineering Science and Mechanics

The principal objective of this research was to develop a method

for determination of the response to random excitation of structures

having geometric nonlinearities. This would account for large or finite

deflections of such structures whereas work done in the past has accounted

only for small deformations leading to linearized equations.

The method of attack was to expand the transverse deflection of the

structure in terms of the eigenfunctions of the linear problem and the

longitudinal deflections in terms of orthogonal functions satisfying the

boundary conditions. The equations governing the series coefficients of

the expansions are then easily derived from the Euler-Lagrange variational

equations and are shown to be nonlinearly coupled. The assumption of

uncorrelated loading then permitted the identification of the phase space

of the series coefficients with a Markoff process and thereby permitted

derivation of the Fokker-Planck equation governing the first order

probability density function of the series coefficients. A stationary

solution of this equation was obtained. While the solution found is

valid for several structural elements with various boundary conditions,

a simply supported beam and a simply supported plate were investigated

in detail. It was found that the type of nonlinearities considered had

the effect of reducing the mean-squared response. Furthermore, it was

found that the mean-sqUared value of the first mode still represented a

good estimate of the total mean-squared deflection of the structure but

in contradistinction to the linear theory the effects of the higher

modes must be considered in calculating the response of the first mode.




The rocket and jet propulsion systems of modern air and space

vehicles have given rise to a host of new problems in the theory of

vibrations. The pressure fields generated from these systems vary

randomly in both time and space over a wide range of frequencies.

These random pressure fields can cause severe vibrations to the vehicles

and their components. The importance of these types of vibrations is

evidenced by the fact that in 1958 and again in 1963 a special summer

program on random vibrations was held at the Massachusetts Institute of

Technology. The two volumes of the book "Random Vibrations," edited by

Stephen H. Crandall (1, 2) were outgrowths of these programs. In the

first volume it was pointed out that "Some of the strictly mechanical

problems are still incompletely understood and the tools for

handling them are relatively crude." This statement stands in evidence

as to the need for research in this important field.

The first satisfactory treatment of stochastic motion was presented

by Einstein who studied the random motion of a free particle. Smoluchowski

generalized this theory to other types of Brownian motion and since then

many important contributions have been made to the theory, notably by

Fokker, Planck, Ornstein, Uhlenbeck, Chandraseklar, Kramers, and others.

Numbers in parentheses refer to the List of References.


On the purely mathematical side, some of the outstanding contributions

have been made by Wiener, Kolmogoroff, Feller, and Doob. Generally

speaking, the theory of linear lumped parameter systems, i.e., systems

governed by ordinary differential equations, excited by stochastic

driving functions has reached a relatively high level of sophistication

and the technique of spectral analysis developed by Rice (5) is generally

adequate for the solution of problems of this type.

For nonlinear lumped parameter systems, the theory is not quite

as refined. For small nonlinearities an approximation procedure

known as "equivalent linearization" has been developed (6, 7, 8). The

procedure consists of replacing the nonlinear system with an equivalent

linear system. The crux of the procedure is to choose the stiffness

matrix of the linear system so that the mean square error of the governing

equations is a minimum. When the nonlinearities are not small this

method obviously fails.

Another approximate procedure applicable to systems with small

nonlinearities is the "perturbation method" (see reference 2). The

idea in this method is to assume a series solution in powers of a

parameter which represents the nonlinear component of the system. A

set of linear differential equations governing the coefficients of the

expansion can then be generated. These linear equations can then be

handled by more established techniques. For systems with nonlinearities

which are not small an alternate approach must be made.

For a complete set of references see entries (3) and (4) in the
List of References.

Such an alternate approach exists in the identification of the

trajectory in the phase space with a Markoff process. This assumption,

together with that of Gaussian input permits the derivation of a partial

differential equation, known as the Fokker-Planck equation, governing

the probability density function of the response of the system. This

concept will be more fully exploited in Chapters IV and V. We mention

here only its shortcomings. These are that the assumption of a Markoff

process implies white noise input which can be questioned with regard

to realizability and furthermore, solutions to the Fokker-Planck

equations are known only in a few special cases (9, 10, 11). Despite

these shortcomings it appears to be the most fruitful approach for

analyzing nonlinear systems excited by stochastic functions.

The picture is worse for continuous systems, i.e., systems governed

by partial differential equations. The linear theory has been attacked

by Eringen (12) who used the generalized Fourier analysis developed by

Wiener (13) to the problems of vibrating beams and plates. Other special

linear continuous systems have been studied by many authors (14-20).

Chapter III contains an outline of the standard approach to these


As in lumped parameter systems, the method of "equivalent

linearization" has recently been used to study continuous systems

governed by equations containing small nonlinearities (21). Until now

there has been no analysis presented for the problem of continuous

systems which are governed by equations containing large nonlinearities.

It is the purpose of this paper to show how, with the assumption of a

Markoff process and Gaussian input, the Fokker-Planck equation can be


used to determine the response of nonlinear elastic plates excited by

uncorrelated stochastic loadings.

The nonlinearities of the plate which are considered are those

arising from geometrical considerations. We still make the usual

assumptions of small strains, shears, and extensions compared to unity

so that Hooke's law is valid, but we allow the rotations to be moderately

large, i.e., small compared to unity but large compared to the shears

and extensions. This will be the result when the deflections of a

plate are not small relative to its thickness, a situation very common

in structures subjected to vibrations produced from jet and rocket




The theory of the large forced vibrations of elastic plates is

presented in this chapter. The equations governing the large

vibrations of beams are obtained by properly simplifying the expression

for the kinetic potential.

2.1. Analysis of Deformation

We start our study of plates with the geometry of deformation.

Let the position of a point in an undeformed body at time to be

designated by x, y, z referred to some rectangular coordinate system

X, Y, Z. After deformation, at time t, the point has moved to a new

locationF 7 1 referred to the same X, Y, Z coordinate system.

Then we can write

= Z+


The functions JL I- ,UT represent the projections onto the X, Y, Z

axes of the displacement at time t of a point which was at the position

x, y, z at time to.

If in equations (2.1) we set x = xo = constant, y = yo = constant,

we obtain the equations


These are the equations of a line which at time to was parallel to the Z

axis. At time t it is some curvilinear line as indicated in Figure 1.

Thus, while the x, y, z form a rectangular coordinate system in the

undeformed body, they form a curvilinear coordinate system in the deformed

body. Therefore, if we speak of a stress 6jL we mean the stress on a


Fi4. 1.- Deformation of
a coordinate line.

X 1-! C b

Undeformed Deformed

Fig. 2.- Deformation of an
element of volume.

a. Lagrangian coordinate system.

b. Eulerian coordinate system.

surface which was originally perpendicular to the Z axis and acting in
a direction normal to this surface (see Figure 2). On the other hand a
stress 5, will be the stress on a surface which in the deformed body
is perpendicular to the Z axis and acting in the Z direction (see Figure 2).
Having dispensed with these preliminary considerations of our
coordinate system, we may now proceed to our analysis of strain. If we
introduce the following parameters

e =e -. y3

I -L +4

w = e ( I-)

W, = CC-

then the components of strain can be defined as (22)

6^ = e t[e1 +(e(+_w)1 +(e -W^-

S= ez, = e,+ey e,)' +(I eI,- I

t = e+tle^1 +j e2+ w+ci( e, -L]

ee- ) +
equations (2.4) can be written as (22)(2.4)

e, e xz-z+ W,)

z = + = e .L + ze -- We) + e,-- ez+L4)

4-c-2 x -W, )( e,,+ y)

E.x z ei.4 e,,(x ez wzd, ) ^x ,x
+( -I e, ^ e7y + ,(

We now make the assumption that the strains and rotations of a
volume element are small compared to unity. The assumption of small
strains permits us to use Hooke's Law which is generally not valid for
large strains. When we add the assumption of small rotations,
equations (2.4) can be written as (22)

4 e -L + +

e= +-2L + C )


A further consequence of our assumption is that the quantities (W ,

L, ~ T can be interpreted as the components of the rotation vector of

a volume element. The retention of their squares in equations (2.5) is

due to the fact that while these quantities may be small compared to

unity they may be large compared to the strains so that their squares

may be of the same order of magnitude as the strains. For flexible

bodies this is quite often the case.

We are now ready to start our analysis of thin plates. We consider

that the plate is of constant thickness h and that the XY plane forms

For an excellent treatment of these points see reference (23),
p. 47.


the middle surface of the plate before deformation. Thus, the Z axis

is normal to the plate. Therefore, after deformation the xy surface

will form the middle surface and the z axis will be normal to that


As is always the case in a "Strength of Materials" type of analysis

we must assume a form for the displacement field. For a thin plate with

moderately large rotations we take,

M%^ =r ) JU(YLt-,)- rL^

--a___Ur (2.6)

VCrCX, UZ, 't)=T CJ-CL U)fk

As was pointed out by Biot (24) equations (2.6) are those of the

von Karman plate theory. They are tantamount to assuming that a straight

line originally normal to the middle surface of the plate remains

straight and normal to the middle surface after deformation. This is the

assumption of classical plate theory. Furthermore equations (2.6) permit

a stretching of the middle surface of the plate which in classical theory

is omitted.

The use of equations (2.6) yields the following results for the

parameters given by equations (2.3)


p ^1 3 a^2

ev -E
3T C) 19 2.
e~y = : -z Jy
aa~ ia a

- kL

+-F2 -zj.
a-^ a4w


C z =


u-? aur

z- -

The rotations in the plane of the plate can be considered fairly small

so that in using equations (2.7) to calculate the strain components we

neglect U .

Thus, making use of (2.7) the strain components become


a -2Ur
-r- -4 CIVI 5CY

= ,L (in,
EIL^ Z 'a I*-4



&y La



The quantities

E. 1, E y, *- are obviously the components of the

strain of the middle surface.

+ I_(J.. L
f 214/

Fu =

__L -L

E =0


_ t -4-!X
au a-><

+- AA_ __
a )a Y


That E.X2 = Ez = 0 points up the fact that we are neglecting
the deflections due to the transverse shear stresses, an approximation
that is valid in most cases.

2.2. Equations of Motion
To derive the equations of motion governing the forced vibrations
of plates we employ Hamilton's principle in the form

s Kdt =(2.10)

where K is the kinetic potential of the system given by

K T-V (2.11)

Here T is the kinetic energy of the system and V the potential.
Equation (2.10) states that the displacement field assumed by a system
is such that the kinetic potential is an extremum.
For an elastic structure V consists of two parts, the internal
strain energy Vs and the potential energy of the external loads Ve.

\4 J [.(2 +.12)
V (2.12)

cs~l~~b~B, +5~



e ff[ j + ( +i as (2.13)

where the first expression is integrated over the volume of the body

and the second over the entire surface of the body. The quantities

E ,L + are the components of the boundary traction in the X, Y,
Z directions, respectively. Since in our analysis we have assumed small

strains, no distinction has been made between integration over the

deformed and undeformed body nor between the stresses and the pseudo-

stresses d- which actually should be used in equation (2.12).

For a plate with no shear stresses on top or bottom, equation (2.13)

reduces to


Ve -[ff u s -A
-h (2.14)

+ J- -+ ur) d.td


Ve= (2.15)




c%(x.c k H2 (2.16)

and the line integral is taken around the entire cylindrical boundary

of the plate. We consider only the case when the forces on the

cylindrical boundary are the forces of constraint so that the second

integral of (2.15) will vanish and we have

Ve jy^^ )t) c ) Ms (2.17)

The retention of the second integral in (2.15) would lead to the boundary

conditions but since these have been established in other papers (25, 26)

we will not derive them here. They will be briefly discussed at the end

of the chapter.

Returning to equation (2.12), the strain energy, with the help of

equations (2.8), can be written as

s -ly

-r -z- 2 Z(2.18)

+i; ][( +() +dtcs


Integration over z yields

\V + N- I + N, E+ (2.19)
V2Nr --9ur Ur
-M -)0. I -ZM ^s


NO = dtLz. M,%. =fz gCZ

hlz h12.

N,,= az M 7 z I7 (2.20)


N1= f dz M% = fza6 dz
and we have assumed

S z =0

The kinetic energy of an elastic plate is given by

T f( + + d (2.21)



or after substitution of (2.6) and integration over z

T= rh r+ ] +

where p is the mass density of the plate and the dot indicates

differentiation with respect to time. The second contribution to the

integral is due to rotatary inertia. Since this only effects the higher

modes of vibration we will disregard it. Thus we write

T f/( A. Z Z) (2.23)

Combining (2.17), (2.19) and (2.23) with (2.10) gives

-4 S

++ N -M,, (2.24)

M MI+2. -ctwl1dS o

Application of the techniques of the calculus of variations to (2.24)

would yield the equations of motion in terms of the stresses and


displacements. The stress-displacement relations would then yield the

equations strictly in terms of the displacements (see reference (26)).

For reasons which will become apparent later we will, with the use

of the stress-displacement relations, obtain the kinetic potential in

terms of displacement only.

For an isotropic, elastic, material the stress-strain relations




E +V)E .~
Solving the first two of these equations for O.- and 0 gives

CILJL C^+4 V ) + V


Equations (2.26) and the fourth equation of (2.25) can be
substituted into equations (2.20) and with the aid of (2.8) the required
integration can be performed. The result is the force-displacement

N 2.
M~~~ LI EK1'^j


+' l

M =a Y() 2



J iz0- 1)

E + IX
I+ 1) 2

+ U)

d ^)

M C ) -


We have again employed the assumption


and further


jza-dz = 0

We can now combine (2.11), (2.17), (2.19), (2.23), and (2.27) to

finally obtain the kinetic potential strictly in terms of displacements.


K Arr itL +V w+ ]

E- ( +I +)

Here, the first term represents the kinetic energy, the second the

potential energy of the lateral load, the third the membrane energy and

the fourth the bending energy.


Having the kinetic potential we could employ Hamilton's principle

and obtain the three equations of motion governing u, v, w. In the

ensuing analysis, we will be more interested in equation (2.28) than

in the equations of motion so that the latter will not be derived here.

The simplification of equation (2.28) to the beam equations can

be made by setting V / ,a/c equal to zero so that we have for the

kinetic potential of a beam undergoing large deformations the expression

F =i+t (2.29)
L l

This equation can be simplified further by disregarding the

longitudinal inertia and by writing

s, L w

so that we have

ac-,wr ~Lcjtd g f 2 (2.30)
K f mI

for the kinetic potential of a beam undergoing moderately large vibrations.


To obtain the linearized equations of motion of a plate we disregard

the membrane stresses and the longitudinal inertias so that the kinetic

potential becomes

K0+rPtLf+ ['I L. -&T. JS (2.31)

Application of Hamilton's principle yields (27)


for the equation governing the small vibrations of an elastic plate.

A linear viscous damping term 13 has been added to account for the

damping phenomena.

2.3. Boundary Conditions

For equation (2.32), governing the small vibrations of elastic

plates, two conditions are needed on w along the entire boundary.

For example, the boundary conditions at a clamped edge are w = 0 and

l/lIn = 0 where n is the direction normal to the cylindrical boundary

in the xy plane. A simply supported edge has the conditions w = 0 and

ew/ di = 0 along the boundary. Other boundary conditions are given in



When we consider the large deflections of plates as governed by

equation (2.24), additional boundary conditions must be prescribed. In

addition to the two conditions on w we need one condition on u and one

condition on v along the entire boundary. The case of a clamped plate

is the simplest for then we have u = v = 0. For a simply supported

edge we have un = 0, where un is the longitudinal displacement normal

to the cylindrical boundary in the xy plane. The other boundary condition

is obtained from the condition that the shear stress along the boundary

4Cn vanishes. From (2.25), (2.8), and (2.9) we have

_+ _n au-_ar --Z u ) (2.33)

Since un = w = 0 along s then

S( )l(2.34)

along s and the condition on us at the boundary is dos/M = 0.

Other types of boundary conditions can be considered. For example

see (25).



The basic aim in this chapter is to review a method of solution

of equations of the type

L(u+) mw = q+rw<-

where L is a linear, spatial, differential operator, the dot indicates

differentiation with respect to time t, and q(x,t) is a function which

is not completely deterministic.

We first must specify what we mean by solving the equations. Since

the load q is not deterministic, i.e., only certain statistical properties

of it are known, then quite naturally all we can expect to know of w are

certain statistical properties. The process or experiment for which w

is the result is said to be a stochastic process and w itself is called

a random variable.

Before attempting to solve equation (3.1) it seems natural to first

give a brief introduction to the mathematical descriptions of stochastic

processes. More complete and rigorous descriptions of such processes

can be found in any number of books (28, 29, 30).

3.1. Stochastic Processes and Probability Theory

From the mathematical point of view a stochastic or random process

is a collection of functions y (t), y (t), yN(t) for which there



exists a probability measure. Each yi(t) is called a sample function or

record and the entire collection is called an ensemble. We relate this

mathematical model to the following physical model. Consider a certain

experiment which can be repeated under similar conditions a large

number of times, e.g., the thermal noise arising across a set of identical

resistors. The outcome of each experiment is a different function y (t)

so that we say the function y(t) which we are measuring is a random

variable. On any given trial we cannot predict the outcome so that only

certain statistical information concerning the process y(t) can be


The basic functions that define a random process are the following

set of probability density functions:

W1(Yl,tl) dyl = probability of finding y in the range (yl, y1 + dyl) at

time tl.

W2(Y1,tl; Y2,t2) dy1 dy2 = joint probability of finding y in the range

(yl, Yl + dyl) at time t1 and in the range (Y2, Y2 + dY2) at time t2.

W3(Yltl; Y2,t2; Y3,t3) dy1 dy2 dy3 = joint probability of finding y in

the range (yl, yl + dyl) at time tl, in the range (Y2, Y2 + dy2) at time

t2 and in the range (y3, Y3 + dy3) at time t3.

We continue on in this way indefinitely. These set of functions

must fulfill the following conditions.

1. VV >,

The use of the word variable for these types of functions is


3. 0 o
W C. ( ,t; ...^,w =J*m)) f *Wn(y..i; ..** m, ...^8*v3cd,*

since each Wn implies all previous W,.

It may happen that under a shift of the t axis the functions are

unaffected. Such a process is called stationary and we have W1(y1) dyl =

probability of finding y in the range (yl, y1 + dyl). W2(Y1,y2; t2,t1)

dyI dy2 = joint probability of finding a pair of values separated by a

time t2 t1 in the range (yl, y1 + dyl) and (y2, y2 + dy2). And so on.

For experimental work the condition of stationarity is almost a necessity.

It is tantamount to assuming a steady state condition, i.e., all transients

of the system have disappeared.

Quite often it is necessary to deal with more than one random

process. That is, we may be concerned with several random variables

Yl(t), Y2(t), YN(t). Defining the process we then have the
following probability density functions:

W1(y11, Y21, YN1, tl) dyll dy21 dN1 = probability that

yl falls in the range (y11, yll + dyll),' 2 in the range (Y21' Y21 + d21)'
S. YN in the range (YNI, YN1 + dyN1) at time tl.

W2 (Yl, Y21i YN1, tl; Y12' Y22' YN2' t2) dyll dYN1 dY12
N2 = joint probability that yl falls in the range (yll, 11 + dyll)

S. YN in the range (YNI, YN1 + dYN1) all at time t, and that yl falls

in the range (YNl, YN1 + dN1), yN in the range (yN2' yN2 + dyN2)

all at time t2.

And so on again. For simplicity we can use vectorial notation

and treat the N variables as components of an N dimensional


vector. Then in place of the above we may write Wl(yl,t), W2(j,tl;

Y2,t2), etc., with m (Ylm, Y2m' YNm)
One of the commonest and most useful density functions is the

normal or Gaussian function. For a random variable y(t) the Gaussian
probability density function is defined as (31),

W t )C C ... =

rf\^ = A< n > = n)

and where is the covariance matrix of elements

3C= TC,, )tl wYnmt

The symbol < > indicates ensemble average and R(tn,t m) is the
correlation function defined by

If, = < i n n >

If the process is stationary then

R Ut*,t 1 = R C *,) =




Tn= fL ='M

and the process is completely defined by the correlation function R(T)

and the mean m. Furthermore, when ensemble averages may be replaced by

time averages (known as the ergodicity property) then

Tur)= e^ T/ ^(3.4)


The Gaussian density function can be extended to cover two or more

stochastic processes. A discussion of this situation can be found in


When a constant parameter linear system is driven by a Gaussian

random process then the output of the system is also a Gaussian random

process. For such a situation, knowledge of the mean and correlation

function of the output permits us to write, down any multi-variate density

function for the output process. If the input is non-Gaussian then, in

general, the output is non-Gaussian and the mean and correlation function

no longer completely define the process. For such a situation no general

method exists for finding the probability density function of the output.

One other useful function in stochastic processes is the power

spectral density. It can be defined as the Fourier transform of the

correlation function. Thus,

o (3.5)

F Wo)PRt dO


The quantity F(to) dr3is the amount of power in the frequency range

(W, u + dw) and hence the name power spectral density.

The inverse transform is

R f r F (3.6)

Note that

o (O) fFCo ) = J _,"
'R(0)F ,I VC)L + (3.7)

= total power of the process.

3.2. Response of Linear Systems

Since the correlation function for linear Gaussian systems is so

important we will outline here a method of obtaining this function.

This will later permit comparison between linear and nonlinear theories

for mean squared values of displacement.

We start with the equation governing the system in the form of

(3.1). We formally seek a solution in terms of normal modes. That is,

we consider

L Cu += o (3.8)

and seek solutions satisfying the boundary conditions. Let these be

Sw.,.t (3.9)
ur C < eC~)


whereUon is the natural frequency of the nth mode. Substituting (3.9)

into (3.8) yields

C (o)= 0 Mo C O( (3.10)

We formally seek a solution to (3.1) in the form

U(X ,t)=o C,C-(r ) (3.11)

Substitution of this into (3.1) gives

L LC+P o< m or = q Li) -M(3.12)

Invoking (3.10) reduces this to

Mm+. + el + YLr 2= 0 -)L (3.13)

A basic property of normal modes is that of orthogonality,

i LO-( (TI)O<( C -= ^^ < ~(3.14)

where Cmn is the Kronecker delta.

Thus, multiplying (3.13) by O(m(x) and integrating over all x

gives the result


^+, a +GLto M=^-J^.^c<=f ) +(3.15)

The solution to (3.15) may be written in the form

,(*. -jf(^*.-t)dt (3.16)


GJ (3.17)

the inverse Fourier transform of

[(AX:- .eFsV
The function hn(T) is the impulse response function of equation (3.15).
From (3.11), (3.15), and (3.16), we have, as a formal solution
for each sample response w(x,t), the result

O( )=- c )JO,,(x' }, -a7 (3.18)

where the lower limit on hn(!-) has been changed to -cX since this
function is zero for 2' < 0.


The cross correlation function is the statistical average of

w(x,t) and w(x',t +21) which from (3.18) is

-x (L-,) d-c, 8t,-J) (1)

The function

4(' ,e, t z.)= ) < \ ($ ^,t7 (3.20)

is the cross correlation function of the load q and it completely

determines the cross correlation of the response IU via equation (3.20).

An important special case of 4 is

-(L. o,, I = SY^L-c) s(I -rj (3.21)

where S is the Dirac delta function. This permits no correlation in

space nor time and is obviously not physically realizeable. However,

for lightly damped systems it represents a fair approximation to reality.


Substitution of (3.21) into (3.19) gives

For the important case of lightly damped systems we have

for all n so that (3.17) becomes

T4 ..


= 0 Z

Substitution of (3.23) into (3.22) would yield the cross correlation

function of the response. Upon setting x = x', Z" = 0, we obtain the

mean squared response at the point x

i,<< M J-





R~iz4 em -I


For a simply supported beam,

2. 2.

so that


For systems with more than one space variable, e.g., plates or

shells, the method of attack is essentially the same. It hinges on the

method of normal modes. For complicated structures these functions are

sometimes difficult to obtain and approximate techniques must be employed.



The classification of random processes is discussed in this

chapter. It is shown how the assumption of a Markoff process permits

derivation of the Fokker-Planck equation governing the conditional

probability density function of a set of stochastic variables. It is

then shown how this equation has been used to obtain the stationary

first-order probability density function governing some nonlinear

lumped parameter systems.

4.1. Classification of Random Processes

In order to discuss the classification of random processes it is

necessary to introduce the concept of conditional probability density

functions. We define these functions in the following manner (see

reference 31):

P2 (Y1,tl 72,t2) dy2 = probability that, if y has the value y at
time t1 then y will have values in the interval (y2, y2 + d2) at time

t2 (t2_0 tl).

Pn (Y',tl; Y2,t2; n-1 tn- 1 n tn) dn = probability that if
y has the values Yl, Y2, n-1 at the respective times tl, t2, .

tn-1 then will have values in the interval (y, yn + dy') at time

tn (tn tn-1 tl).



With these definitions we will then have

W( ,(,t, ;,,ih43 = W (= .,N z(j).I .


and so on.

The Pn must fulfill the following conditions


which follow from definition.

Some further important properties of the conditional probability

density function are:

= /P ,%


+ -t1-wo

ff -- --. _
ca0 4P0

-00 -i0 >v

-E3 / --0 1


so that

t:z (?,t ,^,^) = ,( ^.) C^ (4.4)

since it is certain that y2 = yl at t2 = ti. Here (Y2 Yl) is the

Dirac delta function.

We also have

te (4.5)

so that

t ^C.,i z ^,t,,t^~ ,(4.6)

since there is no statistical dependence between values of y observed

at times sufficiently separate.

We are now ready to start classification of random processes.

The simplest type of process is the purely random process for which

we have


so that from (4.1) we have


*** .C^ ^n*)


This last equation tells us that the purely random process is one in

which any ym and n for tm # tn are statistically independent.

The next more complicated process is known as a Markoff process.

It is the situation in which all the information is contained in the

second-order probability density function W2 (y',tl; 72,t2). For the

definition of the Markoff process we have

PRrld ,SC^ji-~.,t-, n-n,(4. 9)

This equation tells us that the probability that y has values in the

range (n, Yn + dyn) at time tn given that it takes on the values yl1

Y2, Yn-1 at times tl, t2, tn-l respectively depends only

on the value of' at the previous time tn-_.

Substitution of (4.9) into (4.1) gives

.C0 (4. 0)

so that the process is completely specified by P2 (yn-1 tn-l I n, tn)

since W1 (Y1,tl) is found from the relation


We can continue on in this way for more complex problems. Thus

P3 (Yn-2, tn-2; Yn-l, tn-1 0n, tn) defines the next more complicated


process, P4 ( -3' n-3; Yn-2' t-2; 7n-l, tn- Yn>, tn) the next, etc.

However, in this analysis it is the Markoff process with which we are

concerned. Therefore the higher-order processes will not be discussed.

It might seem that equations (4.2) are the only restrictions on

P2 (71, tl 1 Y2, t2). However, for a Markoff process this is not the
case. P2 (Yi, tl -2, t2) must satisfy the Smoluchowski equation,

which we will now establish (see reference 31). We start from the


W3(2 jVk ^- (4.12)

integrate over yo and employ (4.9) so that we have

7. '.O (4.13)



so that upon using this in (4.13) we arrive at Smoluchowski's equation

in the form

TL_ OQt t 7= (4,.15^ ^ *)


This is the basic equation of a Markoff process. In the next

section we show how, with the proper assumptions, this equation can be

used to derive the Fokker-Planck diffusion equation.

4.2. The Fokker-Planck Equation

Before actually deriving the Fokker-Planck equation, let us look

at the differential equations with which we will ultimately be dealing.

Consider an N degree of freedom system governed by the differential


Fi^ ty (4.16)

where y = (y1, 2, YN, Y, y2, yN)' Ym are the N variables

with which we are concerned and ym are their derivatives with respect to

t. Further F is a deterministic vector valued function. This equation

would produce a deterministic trajectory in the 2N-dimensional phase

space and would be completely determined by equation (4.16) and the

specification of y at some time to.
Let us add to equation (4.16) a stochastic forcing function f(t)

so that we now have a set of equations governing each sample of y which

have the form

SF(tIV)+ 1) (4.17)

The trajectory of the phase space is now a stochastic process. Clearly,

the position y(t2) at the end of any infinitesimal interval of time

depends only on the value" (tl) at the beginning of the interval (t2,tl)


and on the stochastic forcing function f acting during this interval.

We now make the assumptions that the forcing function f is Gaussian

with zero mean. We further assume that the forcing functions acting

on the system at any two small consecutive time intervals are

statistically independent. These assumptions make it necessary for us

to take

< ^0 (4.18)

< +(t) ({tC -t) > =!?z (4.19)

where fm is the mth component of f and Rmn is some function of m and n.

Furthermore, since the position at the end of an infinitesimal time

interval depends only on the value y(tl) at the beginning of the interval,

and on the forcing function acting during the interval, which according

to our assumption is independent of the forcing function acting outside

the interval, then the trajectory of the phase space is a Markoff

process. It is completely defined by the conditional probability density

function P2 (yl,t \ Y 2,t2), which must satisfy the Smoluchowski

equation (4.15).

With the proper assumptions equation (4.15) can be used to derive

the Fokker-Planck equation governing P2 (y1,t1 | 2,t2). The derivation

can be found in many places (4, 10, 31). Here we follow the derivation

given in (10).

To start the derivation of the Fokker-Planck equation we consider

the first and second moments of the displacement of the phase point in


an infinitesimal time. These are

J ^C~, k^ 'J ^-'taf(4.20)

We assume that these are of order A t and that all higher moments are

of higher order of A t. The first assumption insures the existence of

the following limits.

1 L(4.21)

CL M- Jatro-,- t ,,%y% J

It has been pointed out in (4, 10) that these assumptions are tantamount

to the assumption of a Gaussian process for the disturbances.

Having made these preliminary assumptions we consider an arbitrary

scalar function R(y) which vanishes sufficiently fast to zero at infinity.

Multiplying the Smoluchowski equation (4.15) by this function and

integrating over the entire phase space gives


and we have interchanged the order of integration. We now develop R(y)

in a Taylor series in (y x),


atid- aid


Upon using (4.23), (4.21) and the assumptions concerning the higher
moments, the right-hand side of (4.22) becomes


C4)P ^


+ o A-0 ,

Integrating by parts, writing y for x and putting the result in

(4.22) gives

IR p.( .U-,t\ -?-4 0c^ 1


,ZN ZN ,Z la ~ i
n Y\?t +01~
b I- V.0 tV
^M-l 2..

0- oA-

x ane

+t r ^


Taking limits as A t-- 0 yields

I 2(4. 26)

Since R(y) is arbitrary, the bracketed expression must vanish,

leaving us with the Fokker-Planck equation.

-n + Y(4.27)

This is a parabolic diffusion equation. The required solution is

the positive one with


If all transients of the system have died out and a steady state

condition has been reached, then P(Yo0ol ',t) W1(Y,t) so that the

Fokker-Planck equation becomes

^ ^T^-^^ -d~m~E 2-T ^ f(4.29)

It still remains to be shown how equation (4.29) can be used to

determine W1 governing the variables of equation (4.17). The connection,

of course, is through the moments dm and dmn. Indeed we had,

Ct)P M^ 1S'LC_ %k-Ai)4;. (4.30)
4-L40 4t vv ')T-2


Now at the beginning of the time interval xm had the value ym and at

the end the value Ym + AYm so that dxm = d(Aym) and we have

SL-- Aiw YL-UC' (4.31)

Similarly, it can be demonstrated that

d, L" -t' % > (4.32)
&t .0o 0 *'&*

These moments can now be computed from the differential equation

(4.17) and thus the corresponding Fokker-Planck equation is fully

derived. The method is best illustrated by examples. Therefore, to

close this chapter we consider two lumped parameter systems which

have been analyzed by this method.

The first system is a one-degree of freedom, nonlinear oscillator.

The equations of motion are

Also (4.33)


Writing yl = y, Y2 = equation (4.33) is equivalent to the two equations

(4.35)The coefficients of the Fokker-Plnck equation re therefore

The coefficients of the Fokker-Planck equation are therefore


d At 0

.... < I A = 0


Substitution of these coefficients into equation (4.29) yields
Substitution of these coefficients into equation (4.29) yields

Ho ,,

- :a- aw,

The solution to this equation as given in (9) is

T Lil



at;o A-k d' a-to at*' o~t

ei ILd~

-t r3 Ir rsu.trZC~~I1VJ.~ =O
b~L LI~TP IL ~d- ~ ~


where C is a normalizing constant. If k(y) is linear, then this

reduces to the Gaussian density function as it should.

The second illustration of the application of the Fokker-Planck

equation to lumped parameter systems consists of a loaded nonlinear

string as analyzed by Ariaratnam (11). The N equations governing the

system are of the form


S[ALlw~, +* LL 2.LJ + m (Ai=o

To simplify these equations, the deriving functions and the response

are expanded in terms of the eigenfunctions of the linear problem. Thus,


Substitution of (4.40) into (4.39) yields

(A) (4.41)
A,+ Wt,, +_% L at o <4241%

This equation may be used to derive the moments appearing in the Fokker-

Planck equation. The solution of the resulting equation, as obtained

by Ariaratnam is


W, C=e~B crsA+ +~fS-1\N+1A
I m= -


where C is a normalizing constant and

This expression can be used to obtain the mean squared displacement of

the various masses on the strong.

The Fokker-Planck equation has been used in the past to solve

several nonlinear lumped parameter systems. In the next chapter the

responses of some nonlinear continuous structures such as beams and

plates are investigated by this method.



In this chapter it is shown how the Fokker-Planck equation can be

used to investigate the finite responses of plates which have been

subjected to white noise excitation. A general solution to the Fokker-

Planck equation is given which is applicable to plates with any

boundary conditions. Detailed solutions are presented for a simply

supported beam and a simply supported plate.

5.1. General Theory

In Chapter II it was shown that the forced vibrations of an elastic

plate are governed by the following equation



K f f(& ++r 4' ^if-Wfc (5.2)

and Vs is the strain energy.

Instead of applying variational techniques to equations (5.1) and

(5.2) and thus obtaining the three equations of motion governing u, v, w

we proceed as follows: We expand each sample function of w and q in a



series of the eigenfunctions of the linear problem. This is valid as

long as each sample of w and q has continuous derivitives up to fourth

order (reference 32, p. 370), a condition we now assume. Thus, we write


The infinite series has been terminated at some N, which is later to be

specified. Of course this invalidates the equality sign of equations

(5.3) and (5.4). However, as long as the infinite series representing

q and w converge, the finite sum can be made as accurate as desired by

properly selecting N.

Similar to equation (3.10) the eigenfunctions of the plate must

satisfy the equation

( v >- ,, W = o (5.5)

where the \rm are the eigenvalues determined from the frequency

equation. In addition the Wmn must satisfy the appropriate boundary


Each sample function of u and v is also expanded in an infinite

set of functions. We choose this set to be orthogonal, to satisfy the

boundary conditions and to be such that term by term differentiation


of the infinite series is possible. Thus,


These infinite series have also been terminated with the previous argument

concerning convergence still applying. For convenience the same N has

been chosen for all series.

An example of a proper expansion would be

kL Ot) I MnnIL AAM Tlt- (5.7)

for a clamped rectangular plate of sides a and b and with the origin of

the coordinate system at a corner of the plate. This would be a double
mt x nlvy mV x nITy
Fourier series with terms such as sin a cos cos --- sn b

cos m x os n y omitted. These terms can be omitted if we consider
a b
the extension of u onto the intervals (-a,0) and (-b,0) to be an odd

function. Now since (5.7) is a Fourier sine series, which is zero at

x = 0, x = a, y = 0, y = b, then it is continuous throughout the entire

xy plane and may be differentiated term by term (see Theorem I on page

137 of reference 32). Therefore, for a clamped rectangular plate

equation (5.7) is a series expansion of u which satisfies the boundary

conditions and the conditions for at least one differentiation provided

each sample u is integrable in the Lebesgue sense, a condition we now



If we'now substitute equations (5.3), (5.4) and (5.6) into equation

(5.2) and perform the indicated integration, we obtain


where V is the strain energy Vs after integration and

b-w- c)n= j(C U)-iS (5.9)

Here the orthogonality property of the eigenfunctions has been used.

To satisfy equation (5.1) we consider the umn, Vmn, Wmn as

generalized coordinates and apply the Euler-Lagrange variational

equations, which are of the form

d a( .. (5.10)

where the 9m are the generalized coordinates. Application of (5.10)

to (5.8) yields the following set of differential equations


h"" I
^ w ~ ~ p~_ _



xJ- ph /3~n t~


where we have introduced the same linear viscous damping term of the

linear problem (see equations (2.41) and (3.15)).

We now write (5.11) as a set of first-order equations by setting

'r T- (5.12)

Equations (5.11) then become,

4Lrmtr = ",a" Vf

A m M ~

hQa.r 5U~m



0 &

~Tm ~


k b

-I I

These equations constitute a set of stochastic differential equations.

They are of the same form as (4.17). The vector y has the components

umn, umn vmn, Vmn, Wmn, mn for all values of m and n. This is a

total of 6N2 components. The stochastic driving functions fm are

equivalent to the qmn. Therefore, if we assume that the qmn belong

to a Gaussian random process and that

< % > =0

then the arguments of Chapter IV hold and the trajectory in the phase

space of equations (5.13) constitute a Markoff process. As was

demonstrated in Chapter IV, this implies that the stationary probability

density function W1(y) must satisfy the following Fokker-Planck equation

a ..I WI) m1f



Al*-0 A

mn- A~bl-0 o X

Equations (5.14) imply that

'e'- ,i^^c;,y


The moments dm and dmn of equation (5.16) can be calculated from

the differential equations (5.13). For bookkeeping purposes the

following notation is introduced:

cl = Aff. -- L A U,, >

C9 AA-)

VlA X f.1 < >




-tn )I~-

A-t -o At-

c _* A
Cm rs

Atco at

At-O 1-f

LX lff %

'A LmL AV T >

( A ,,kTA1VS >

/'A U);' a Ur ">

Then, with the aid of equations (5.13) these moments can be calculated

to be

^ v --Jt,
Ai-^o A*

Alt-0 At

-1 k o

- I,, a
10 ,.,V." M


c^ Tftn.


dUr -- wn.


CLvt7L r1,5
J 1-

4tht tat
;I v ( J ^ (^>^ W d:LLA >
~c ;k

all other dmnrs = 0.

Substitution of these expressions into equation (5.15) yields the
following Fokker-Planck equation:

-in~ m L i ~~~l

) (5.21)

--a f S ,), ,

I' V
V V.

'V% bC
9e n W

+ a

A-t -p0

- 12"'eS &

S (



plyvr r

+ ~3,

(,-M V,.) -


This is indeed a complex equation for which there exists no

standard technique for obtaining the solution. By a process of trial

and error the author was fortunate in obtaining a solution for the case

in which

S= ZN (5.22)

This imposes the restriction that the load is completely uncorrelated

in time and space, i.e.,

<^G^Z=tN' ^:-x-)t-') t) (5.23)

The solution to (5.21) with the restriction of (5.22) is

S^Jv4 Z < ^z ^(5.24)

where C'' is a constant to be obtained from the normality condition.

The general uniqueness theorem to the Fokker-Planck equation given in

(33) tells us that this is the only solution.

It is to be noted that the distribution of the velocity variables

is Gaussian. If we integrate over these variables from 0 to +00 we


obtain the probability density function of the umn, mn and wmn. Thus,

W,,( ,....,...3,... = C'^fL- ] (5.25)

The constant C' is obtained from the normality condition.

We are ultimately interested in the probability density function

of Wmn since from it we can determine quantities such as mean squared

displacement, mean squared stress, etc. To obtain this we integrate

equation (5.25) over all umn and vmn so that

v^ (3 .. =r .pf l TeL-Jl^ cku. (5.26)
-co -00

This is as far as can be gone in such generality. To proceed further

we need to compute V. Several important cases are presented in the

next section. It should be noted in passing that when the load has some

spatial correlations it would be advisable to employ some approximate

technique to solve equation (5.21) governing the probability density

function of the modal amplitudes.

5.2. Some Special Cases

a. Simply Supported Beam

The kinetic potential of a beam with moderately large vibrations

is given by equation (2.29) so that the strain energy is

O .
T. V^ dW EP (5.27)


The eigenfunctions of the linear problem are sin mrx so that we

UC7^,t = an XJ (5.21
n=r V-
The first four equations of (5.20) are satisfied identically since

u = v = 0 and since the strain energy has been expressed independently

of u and v.

The probability density function of the wm is then given by

W,CC, -94 fc'^ L-,1 (5.2(

with V given by (5.27) after substitution of (5.28). Thus

where .. [- i (5.3

E- NoL T'


N,- N=L





Here 6G is the mean square deflection of the first mode of the linear
problem obtained by setting =0. We have replaced N'/L by No, which

is interpreted as the power spectral density of the average load acting

on the beam, i.e.,

Go L L

-ao 0

CO (5.32)

N'L 9&t) =N'

where q is now to be interpreted as force per unit length and we have

used the fact that

The static deterministic counterpart of No for a uniformly loaded beam
is, of course, q i.e., the square of the constant load per unit length.

It is seen from equation (5.30) that the nonlinearity of the beam

causes the probability density functions of the modal amplitudes to

become non-Gaussian. Furthermore these variables are no longer

statistically independent.

The mean squared response of the beam is given by

NN A (5.34)
< r 2,= = Y < QA inr M Aw AM I.
15 V%'- V I k- L.



cp oo
= f W....J U)W,(j...)d ...dtj (5.35)
-oD -0o

Because equation (5.30) is an even function of the wn we will have,

< bx- W. >= 0 r-*V (5.36)

so that the modal amplitudes are linearly independent.

The mean squared stress in the beam depends upon < ( which

is given by

2 -

In the linear problem, 4 wm 2 is of the order i/m4 so that the

infinite series does not converge. This was discovered by Eringen (12)

and attributed to two factors: (a) The S functions appearing in

equation (5.23) and (b) The inadequacy of Bernoulli-Euler beam theory.

To obtain expressions for mean-squared stresses, Samuels and Eringen (34)

investigated a Timoshenko beam and Crandall and Yildiz (35) later

investigated many different beam models with different damping mechanisms.

It was found that these more refined theories produced finite mean-

squared stresses.


Since we have only an integral representation for wm2 we

cannot rigorously investigate equation (5.37). However, there is

really no reason to believe that the introduction of the membrane

stresses would cause (5.37) to converge. To investigate the mean-

squared stresses in the nonlinear problem, it would be most desirable

to consider a more refined beam model.
An approximate expression for wm > can be obtained by

substituting the expression

[. R \..... L~j (5.38)

into (5.35) and performing the required integration. The result, after

some juggling, valid only for "moderately large deflections" is

150, r -^ (5.39)

b. Simply Supported Plate

We consider a rectangular plate simply supported on all four sides.

For the linear problem the eigenfunctions are

Av w 2L rt- .(5.40)


where a and b are the lengths of the sides of the plate in the x and y

directions respectively.

Earlier in the chapter it was shown how, for a clamped plate, a

double Fourier series could be constructed for the functions u and v,

which satisfied the boundary conditions and the conditions for

differentiability. For a simply supported plate we define the extension

of u on the interval x = (-a,O) to be odd and on the interval y = (-b,0)

to be even. Then we can write

So=o 0-. b

Now by Theorem I of reference (32) it is evident that this series can

be differentiated with respect to x or y at least once. The summing of

the series with n = 0 is necessary since a cosine expansion must include

a constant term. For convenience we have also started the summation

at m = 0.

For the function v we can proceed along the same lines. We may

therefore write for the three displacements

,'-o Co a. "


For uniformity we have started all sums at m = n = 0 and terminated

them at m = n = N-1 so that the total number of terms of um is still


The kinetic potential of a plate undergoing large vibrations is

given by equation (2.28). The strain energy is therefore given by


where the Ex' Ey' xy are given by (2.9). If we now substitute

(5.42) into (5.43), perform the required integration and substitute

the result into (5.25) we obtain, after some rearranging

\l CAL ^ 1 o,' -n ==

n~h V

-vr A& %A. % 'V -*


vM~ (E t^ fftE 2ff.rr m a. (\- vb) b .
x -a?- 0>
+ ltp HEN' E' As) ~1 ev s A

tErrEEt ^.ur^^ K^^.1) x

+1.L ULn I

+___ +ErZZ AE t-s- vLl

- "Ysv- k, + I.L ,
w p r 2 Wi. A- W r '
OL~t ftf ":,b .%Cp~r 'D r V

'of lp 1 Ar C P
ai61- '- ^P- "


L ,pY S= Imp rE, pF, + ImvF4,E' 2 tfsQe "rus
flL 4L I(O

+ per oE .,s F "p
|L~ lz

I to-

^.C\-u ) +


i-rfvwfYpr .Wv

[f,[k r 2*l

,- V%, 15 Ft E I
K nItpos p F= mpt- t TV


E = fcA --


F =. 1 = CMA'r

A -l cra mn Y
Ampri&t-. M7 --

P" 21- rn rIT d"34

aC.UL a.r

b 4s <


AIA, rnlt.
Awn -nA

c C0a4 I L
OL-- o-


0 Qb b Va \^ ,,

104 I CsAu


C"pr^ C^^CC^^. AAfti rA Ah2Sly

~Y=l 09e~OL CL CL~~~: AR~ ~nty ab
es A-r. = C4IE:*- AAMI
b b

a S

co& w


To obtain the probability density function of the w alone we
need to integrate out the umn and the vmn. If we make use of the


which is given on page 64 of reference (36), then

W ,c(v ,,(A. = -"i ,,u ,...., ..1 ,...

becomes, after some manipulation

becomes, after some manipulation


cf1 AI F-. i- x"
C K [ 4

L b


*PJM^1 -iu

+ L 0 9) b VVL
AV.-Z to 20.

a6 .b

LVW2E.n"2EEE WpWr.( Kj,

Lw pys) )I

4N* #A. s r % p % --

~~(ri-LD) 1v.
0. ~


;L L*

& S

7 -

-00 '- *ae brldr:


If we integrate out the vmn, again making use of (5.47), we finally


W, ic (-

t3 E tTO
I ,C- )L

&& U 0
N-;O w~.:o

+CEL EEEE r Tpi o
v (C & j p n. o f%


O. z2.

____ ( Tj-)L_____________ l
,,' (-V) 6 Q,. -4 V

a- z.b

where as before No is the power spectral density of the average load

acting on the plate, i.e.,

/ s


N -00

N0 = N'/ab

jl0 i^O


Again, the effect of the nonlinearity is to cause the modal

amplitudes to become statistically dependent and non-Gaussian. The

mean-squared displacement of the plate is given by

2.QA/y S A^Tl^AA MT_ Xv^ (5.51)

%V b b


--O -"

Having reduced the mean-squared displacement to quadratures, we must

again stop and seek some approximation to finish the job. This is left

to the next chapter.

As was the case in the beam, the mean-squared stresses for the

linear plate do not converge and so will not be investigated here.

The investigation of the large random vibrations of a simply

supported beam and simply supported plate have been reduced to quadratures

in this chapter. We leave to the next chapter the numerical investigation

of specific cases.



In the previous chapter the probability density function of the

modal amplitudes of a plate undergoing large deflections was derived.

Detailed calculations were presented for a simply supported beam and

a simply supported plate. The mean-squared deflection of these

structures was reduced to quadratures. For the beam an approximate

formula was developed.

In this chapter we investigate more closely some of the results

of the previous chapter. For the beam, the linear, approximate

nonlinear and integral representation of the mean-squared displacement

are compared numerically for a range of the parameters. For the plate

the mean-squared displacement of the first mode, which is the

predominant term, is numerically investigated for different aspect

ratios b/a.

6.1. Simply Supported Beam

In Figure 3, the mean-squared deflection at the center of the

beam .SZ as determined by the linear theory, the approximate formula

(5.39), and numerical integration of equation (5.35) is plotted against

rB with R = 1/2. For the range of Id considered it was found, as

in the linear theory, that sufficient accuracy is obtained with N = 1.

It was also found that the approximate formula is valid over a small











o 0




P14 ,


range of 6o As long as the difference between linear and nonlinear

theories is not greater than 10 per cent, this formula gives an excellent

estimate of the true mean-squared deflection. This is to be expected

in view of the approximation made in deriving formula (5.39).

Also plotted in Figure 3 is the mean-squared deflection determined

by numerical integration with N = 1 and R = 1/4 and 1/8. These curves

indicate that for lower values of the radius of gyration the curve of

the nonlinear theory begins to deviate sooner from the straight line of

the linear theory. This is to be expected since with diminishing R the

role of bending dimishes. Furthermore, if the beam were rectangular

then h2 = 12R2 so that

Z -{ (6.1)

Now for R = 1/2, the nonlinear theory begins to appreciably deviate

from the linear theory at B2 = 0.3 or at = 0.1. For R = 1/4 we

have appreciable deviation at Q-g2 = 0.009 or S= 0.11, and for

R = 1/8 at d-B2 = 0.004 or P = 0.14. Since, as is well known, the

linear theory of beams is valid only for w/h

In Figure 4 the mean-squared deflection as determined from the

linear theory and from numerical integration of equation (5.35) is

plotted against larger values of Bo2 for R = 1/2. Of the three curves,

the uppermost represents numerical integration with N = 1. The lowest

curve represents 4 wl2 > as determined by numerical integration with

N = 3. The middle curve represents the total deflection, i.e.,


- c;o CO







0 0




O a



4 l2 + with N = 3. Two observations are evident from

these curves. Firstly, consideration of the beam as a one-degree of

freedom system (i.e., N = 1) is no longer valid for such large values
2 2
of -Bo2. Secondly, as in the linear theory, < wl2 > gives a

fairly close estimate (a few per cent) of the total mean-squared
deflection. However, in computing ( wl > it is now necessary to

consider the effects of the second and third modes. That is to say,

the nonlinear coupling is so strong, for the range of parameters

considered in Figure 4, that the second and third modal coefficients

have a significant effect on the mean-squared value of the first modal


It should be pointed out that for a rectangular beam with R = 1/2

we have h2 = 3, so that with B = 1.0, the highest value in the

graph of Figure 3, we have

C~~ II z ir ,6


This value is not in excess of the applicability of the nonlinear theory

considered nor of the values of practical interest.

6.2. Simply Supported Plate

The first mode is the predominant one in the calculation of the

mean-squared displacement of the plate. Therefore, if in equation (5.50)

we consider only the first mode we have


+ U .
L -


__ ( xyWn%%i

0- +-i. t- Pw),b-rn
-- l

- -b
b14 CL"

1 L + I n
- C (*>**'' ^ T1

Z. o.b
+1- 2.
6't a b~


L =.n -



-j- -


C. I
3tZIn -h g- + +
8b..> a ^ 0 'g.',Q




0. Z.b



+~v-)i + "L ,- %w ...
-. z. b

---- ---

+ &

W1 C,) =


EhrrY hL
_gEhrr' I Y8
2NoC\- ri



and all other Lnllll = Kmnllll = 0 so that (6.3) reduces to

\N, C =C PA

r -



Wit E
+ ~lY8X

Ia b
o,. I=>

I..+ 2- ,

L +




I of = b/a

This is of the form

^^}- ^^C


PELe +,
la^Ho ^-i)" O^

where ) and o( are easily calculable once V, a and b are known.

The mean-squared deflection of the plate is



vt% tII

(I, Cy\~L
o- Lt,


Unfortunately this integral is not tabulated and we must resort to

numerical methods.

In Figure 5 the mean-squared deflection at the center of the plate
T as determined from numerical integration of (6.7) is plotted
against "2 for 1) = 1/3, h = 1, and two difference aspect ratios
of = 1 and 2 = 2. For larger aspect ratios the linear and nonlinear

theories give higher mean-squared deflections. This is what is to be

expected since increasing the aspect ratio is equivalent to moving

apart one pair of supports of the plate.

Also evident from Figure 5 is that the nonlinearity of the plate

causes a reduction in the mean-squared deflection. Furthermore, for

the smaller aspect ratio the percentage deviation of the nonlinear

theory is slightly greater than for the higher aspect ratio.

In this chapter the mean-squared deflection of the beam and plate

has been investigated by numerical integration. The results indicate

reduction of these quantities and a significant coupling of the modes

for sufficiently large deflections.


J 0 Co o o
c 0 0

o 6 6











a '







This paper has been devoted to the study of the random vibrations

of some nonlinear elastic systems. The equations of motion of thin

elastic bodies with large deflections have been presented in variational

form. General truncated series expansions of the middle surface

displacements have been performed and the Euler-Lagrange variational

equations have been used to obtain the nonlinearly coupled differential

equations governing the series' coefficients..

The assumption of uncorrelated, Gaussian loading permitted the

phase-space of the series' coefficients to be identified with a Markoff

process. This, in turn, permitted the derivation of the Fokker-Planck

equation governing the probability density function of the series'

coefficients. The moments for this equation were obtained from the

governing set of differential equations.

A solution to the Fokker-Planck equation for the case of spatially

uncorrelated white noise has been obtained. Several special cases have

been worked out in detail. In particular, the probability density

functions for a simply supported beam and a simply supported plate have

been computed. The mean-squared displacements of these systems were

reduced to quadratures. For the case of the beam, an approximate

formula was developed.



Numerical integration of the exact expression for the mean-squared

displacement of the beam has shown a reduction of this quantity as

compared to the linear theory. Calculations have also shown a limited

range of applicability of the approximate formula. Furthermore, it was

shown that the first mode still gives a good estimate of the mean-squared

response but the coupling effect of the modes is so important for

sufficiently large deflections, that the effect of the higher modes must

be taken into account when computing the mean-squared value of the first

mode. Numerical integration also showed a reduction of the mean-squared

deflection of the plate.

While this paper has presented a method of attacking the random

vibrations of elastic systems with large deflections, it must be

pointed out that the solutions presented are valid only for completely

uncorrelated loadings. This is, of course, physically unrealizeable.

Nevertheless, the results of this analysis should give some insight into

the problems of the random vibrations of nonlinear elastic systems.


1. Crandall, S. H., ed. Random Vibration. Vol. 1, the M.I.T. Press,
Cambridge, Massachusetts, 1948.

2. Crandall, S. H., ed. Random Vibration. Vol. 2, the M.I.T. Press,
Cambridge, Massachusetts, 1963.

3. Unlenbeck, G. E., and Ornstein, L. S. "On the Theory of Brownian
Motion," Phys. Rev., Vol. 36, 1930, pp. 823-841.

4. Wang, M. C., and Unlenbeck, G. E. "On the Theory of Brownian
Motion II," Rev. Mod. Phys., Vol. 17, 1945, pp. 323-342.

5. Rice, S. 0. "Mathematical Analysis of Random Noise," Bell System
Technical Journal, Vol. 23, 1944, pp. 282-332; Vol. 24, 1945,
pp. 46-156.

6. Booton, R. C. "The Analysis of Nonlinear Control Systems with
Random Inputs," Proc. Symposium on Nonlinear Circuit Analysis,
Vol. 2, 1953, pp. 369-391.

7. Caughey, T. K. "Response of Van Der Pol's Oscillator to Random
Excitation," J. Appl. Mech., Vol. 26, 1959, pp. 345-348.

8. Caughey, T. K. "Random Excitation of a Loaded Nonlinear String,"
J. Appl. Mech., Vol. 27, 1960, pp. 575-578.

9. Chuang, K., and Kazda, L. F. "A Study of Nonlinear Systems with
Random Inputs," Trans. Am. Inst. Elec. Engrs., Part II,
Applications and Industry, Vol. 78, 1959, pp. 100-105.

10. Ariaratnam, S. T. "Random Vibrations of Nonlinear Suspensions,"
J. Mech. Engrg. Sci., Vol. 2, 1960, pp. 195-201.

11. Ariaratnam, S. T. "Response of a Loaded Nonlinear String to
Random Excitation," J. Appl. Mech., Vol. 29, 1962, pp. 483-485.

12. Eringen, A. C. "Response of Beams and Plates to Random Loads,"
J. Appl. Mech., Vol. 24, 1957, pp. 46-52.

13. Wiener, N. "Generalized Harmonic Analysis," Acta Mathematica,
Bd. 55, 1930, pp. 117-258.

14. Ornstein, L. S. "Zur Theorie der Brownschen Bewegung fur Systeme,
worin mehre Temperaturen vorkonmen," Zeltschrift fur Physik,
Bd. 41, 1927, pp. 848-856.



15. Van Lear, G. A., and Unlenbeck, G. E. "Brownian Motion of Strings
and Elastic Rods," Phys. Rev., Vol. 38, 1931, pp. 1583-1598.

16. Press, H., and Houboldt, J. C. "Some Applications of Generalized
Harmonic Analysis to Gust Loads on Airplanes," J. Aeronautical
Sci., Vol. 22, 1955, pp. 17-26.

17. Thomson, W. T., and Barton, M. V. "Response of Mechanical Systems
to Random Excitation," J. Appl. Mech., Vol. 24, 1957, pp. 46-52.

18. Lyon, R. H. "Response of Strings to Random Noise Fields," J. Acoust.
Soc. Am., Vol. 28, 1956, pp. 391-398.

19. Nash, W. A. "Response of an Elastic Plate to a Distributed Random
Pressure Characterized by a Separable Cross Correlation,"
Tech. Note No. 1, Contract No. AF 49(638)-328, Engrg. and
Industrial Exp. Sta., Univ. of Fla., Gainesville, 1961.

20. Bogdanoff, J. L., and Goldberg, J. E. "On the Bernoulli-Euler Beam
Theory with Random Excitation," J. Aero-Space Sci., Vol. 27,
1960, pp. 371-376.

21. Caughey, T. K. "Response of a Nonlinear String to Random Loading,"
J. Appl. Mech., Vol. 26, 1959, pp. 341-344.

22. Novozhilov, V. V. Theory of Elasticity. English translation, Israel
Program for Scientific Translations, Jerusalem, 1961.

23. Novozhilov, V. V. Foundations of the Nonlinear Theory of Elasticity.
English translation, Graylock Press, Rochester, N. Y., 1953.

24. Biot, M. A. "Elastizitatstheorie zweiter Ordung mit Anwendungen,"
Z. a. M. M., Bd. 20, 1940, pp. 89-99.

25. Wang, C. "Nonlinear Large Deflection Boundary Valve Problems of
Rectangular Plates," N.A.C.A. TN-1425, 1948.

26. Herrman, G. "Influence of Large Amplitudes on Flexural Motion of
Elastic Plates," N.A.C.A. TN-3578, 1956.

27. Timoshenko, S., and Woinowsky-Kreiger, S. Theory of Plates and
Shells. 2nd edition, McGraw Hill, New York, N. Y., 1959.

28. Doob, J. L. Stochastic Processes. John Wiley and Sons, New York,
N. Y., 1953.

29. Kolmogorov, A. N. Foundations of the Theory of Probability. 2nd
edition of English translation, Chelsea, N. Y., N. Y., 1956.


30. Feller, W. Probability Theory and Its Applications. John Wiley
and Sons, New York, N. Y., 1950.

31. Middleteon, D. An Introduction to Statistical Communication Theory.
McGraw Hill, New York, N. Y., 1960.

32. Tolstov, G. P. Fourier Series. English translation, Prentice Hall,
Englewood Cliffs, N. J., 1962.

33. Caughey, T. K. "Derivation and Application of the Fokker-Planck
Equation to Discrete Nonlinear Dynamic Systems Subjected to
White Random Excitation," J. Acoust. Soc. Am., Vol. 35,
pp. 1683-1692.

34. Samuels, J. C., and Eringen, A. C. "Response of a Simply Supported
Timoshenko Beam to a Purely Random Gaussian Process,"
J. Appl. Mech., Vol. 25, 1958, pp. 496-500.

35. Crandall, S. H., and Yildiz, A. "Random Vibrations of Beams,"
J. Appl. Mech., Vol. 29, 1962, pp. 267-275.

36. Grobner, W., and Hofreiter, N. Integraltafel. Zweiter Teil,
Bestimmte Integrale, Springer-Verlog, Wien, 1958.


Richard Edgar Herbert was born November 11, 1938 at Mount Vernon,

New York. In June 1956 he was graduated from A. B. Davis High School.

He attended the Cooper Union School of Engineering, a privately

endowed, tuition-free college. Upon receiving the degree Bachelor of

Civil Engineering in June 1960, he enrolled as a half-time graduate

student in the Department of Engineering Mechanics at the University of

Florida. He was graduated with the degree Master of Science in

Engineering in February 1962. In September 1962 he embarked upon

doctoral graduate work.

The author was a half-time Instructor from September 1960 to

February 1962 in the Department of Engineering Mechanics except for

the summer of 1961 when he was an Assistant in Research. From February

1962 to September 1962 he was a full-time Research Associate in the same

department. When he started his doctoral studies he was granted an

Engineering College Fellowship and subsequently a Ford Foundation


Richard Edgar Herbert is a member of Sigma Xi and Omega Delta Phi.


This dissertation was prepared under the direction of the chairman

of the candidate's supervisory committee and has been approved by all

members of that committee. It was submitted to the Dean of the College

of Engineering and to the Graduate Council and was approved as partial

fulfillment of the requirements for the degree of Doctor of Philosophy.

April 18, 1964

Dean, College of Engineering
Dean, College of Engineering

Dean, Graduate School



4t / 7, .-- -

-. a,.sJ_ .[,

-- --,, --'- -.


3 1262 08553 79091111111111111 I
3 1262 08553 7909