RANDOM VIBRATIONS OF NONLINEAR
ELASTIC SYSTEMS
By
RICHARD EDGAR HERBERT
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE UNIVERSITY OF FLORIDA
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE
DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
April, 1964
ACKNOWLEDGMENTS
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.
TABLE OF CONTENTS
ACKTIOWLEDGMENTS .
LIST OF FIGURES .
ABSTRACT .
CHAPTER
I. INTRODUCTION . .
II. THEORY OF PLATES . .
2.1. Analysis of Deformation .
2.2. Equations of Motion .
2.3. Boundary Conditions .
III. THE RESPONSE OF LINEAR SYSTEMS TO RANDOM EXCITATION .
3.1. Stochastic Processes and Probability Theory .
3.2. Response of Linear Systems .
IV. THE FOKKER-PLANCK EQUATION AND ITS APPLICATION TO SOME
NONLINEAR LUMPED PARAMETER SYSTEMS .
4.1. Classification of Random Processes .
4.2. The Fokker-Planck Equation .
V. APPLICATION OF THE FOKKER-PLANCK EQUATION TO NONLINEAR
ELASTIC SYSTEMS . .
5.1. General Theory . .
5.2. Some Special Cases ... .
a. Simply Supported Beam .
b. Simply Supported Plate .
iii
Page
ii
v
vi
1
5
5
13
22
24
24
29
35
35
40
49
49
59
59
63
. .
Page
CHAPTER
VI. NUMERICAL INVESTIGATION .
6.1. Simply Supported Beam
6.2. Simply Supported Plate
VII. CONCLUSIONS .
LIST OF REFERENCES .
BIOGRAPHICAL SKETCH .........
. .
.
. .
. .
. .
. .
LIST OF FIGURES
Figure Page
1. DEFORMATION OF A COORDINATE LINE 6
2. DEFORMATION OF AN ELEVEN[ OF VOLUME ... .. 6
3. MEAN-SQUARED DEFLECTION AT CENTER OF BEAM FOR SMALL
NONLINEARITIES .... ... .. 72
4. MEAN-SQUARED DEFLECTION AT CENTER OF BEAM FOR LARGE
NONLINEARITIES . ... 74
5. MEAN-SQUARED DEFLECTION AT CENTER OF PLATE ... 79
Abstract of Dissertation Presented to the Graduate Council
in Partial Fulfillment of the Requirements for the
Degree of Doctor of Philosophy
RANDOM VIBRATIONS OF NONLINEAR ELASTIC SYSTEMS
By
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.
vii
CHAPTER I
INTRODUCTION
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.
-2-
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
problems.
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
-4-
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
engines.
CHAPTER II
THEORY OF PLATES
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+
(2.1)
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
(2.2)
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
Z
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
+
(2.3)
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 )
(2.5)
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.
-10-
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
surface.
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)
-11-
p ^1 3 a^2
ev -E
3T C) 19 2.
e~y = : -z Jy
aa~ ia a
- kL
+-F2 -zj.
a-^ a4w
(2.7)
C z =
a
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
-12-
a -2Ur
-r- -4 CIVI 5CY
= ,L (in,
EIL^ Z 'a I*-4
(2.8)
where
&y La
~%
+2
(2.9)
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
I
+,la"~
b(itI
_ t -4-!X
au a-><
+- AA_ __
a )a Y
-13-
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.
Thus,
\4 J [.(2 +.12)
V (2.12)
cs~l~~b~B, +5~
-14-
and
e ff[ j + ( +i as (2.13)
SL
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
h
h/l.
Ve -[ff u s -A
-h (2.14)
+ J- -+ ur) d.td
or
Ve= (2.15)
kit
a.4
-15-
where
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
-16-
Integration over z yields
\V + N- I + N, E+ (2.19)
S
V2Nr --9ur Ur
-M -)0. I -ZM ^s
where
NO = dtLz. M,%. =fz gCZ
hlz h12.
N,,= az M 7 z I7 (2.20)
hli
N1= f dz M% = fza6 dz
-h/2.
and we have assumed
k/l
S z =0
-h/i
The kinetic energy of an elastic plate is given by
T f( + + d (2.21)
V
-17-
or after substitution of (2.6) and integration over z
T= rh r+ ] +
S(2.22)
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)
S
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
-18-
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
are
=J-
S+U)E
(2.26)
E +V)E .~
Solving the first two of these equations for O.- and 0 gives
CILJL C^+4 V ) + V
(2.26)
-19-
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
relations
N 2.
M~~~ LI EK1'^j
N-
+' l
M =a Y() 2
(2.27)
Here,
J iz0- 1)
E + IX
I+ 1) 2
+ U)
d ^)
M C ) -
-20-
We have again employed the assumption
o
_lNz
and further
h/z
jza-dz = 0
-h/z.
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.
Thus,
K Arr itL +V w+ ]
(2.28)
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.
-21-
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.
-22-
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
S
K0+rPtLf+ ['I L. -&T. JS (2.31)
Application of Hamilton's principle yields (27)
(2.32)
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
(27).
-23-
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).
CHAPTER III
THE RESPONSE OF LINEAR SYSTEMS TO RANDOM EXCITATION
The basic aim in this chapter is to review a method of solution
of equations of the type
(3.1)
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
-24-
-25-
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
determined.
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 >,
i.
The use of the word variable for these types of functions is
traditional.
-26-
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
-27-
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 ... =
where
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 *,) =
(3.3)
=-R(.C)
-28-
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)
-T
The Gaussian density function can be extended to cover two or more
stochastic processes. A discussion of this situation can be found in
(31).
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
FjWJ~~c~-ex?01
-29-
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)
-d0
Note that
o (O) fFCo ) = J _,"
'R(0)F ,I VC)L + (3.7)
-0C
= 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~)
-30-
whereUon is the natural frequency of the nth mode. Substituting (3.9)
into (3.8) yields
2
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)
n=t
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
-31-
^+, a +GLto M=^-J^.^c<=f ) +(3.15)
The solution to (3.15) may be written in the form
,(*. -jf(^*.-t)dt (3.16)
where
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.
-32-
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)
LL
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.
-33-
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 ..
aon
= 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
L
i,<< M J-
(3.24)
ct3
=-
nzi
(3.22)
(3.23)
R~iz4 em -I
-34-
For a simply supported beam,
2. 2.
so that
(3.25)
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.
CHAPTER IV
THE FOKKER-PLANCK EQUATION AND ITS APPLICATION TO SOME
NONLINEAR LUMPED PARAMETER SYSTEMS
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).
-35-
-36-
With these definitions we will then have
W( ,(,t, ;,,ih43 = W (= .,N z(j).I .
(4.1)
and so on.
The Pn must fulfill the following conditions
(4.2)
which follow from definition.
Some further important properties of the conditional probability
density function are:
= /P ,%
(4.3)
+ -t1-wo
ff -- --. _
ca0 4P0
-00 -i0 >v
-E3 / --0 1
-37-
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)
+-t.a-
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
(4.7)
so that from (4.1) we have
(4.8)
*** .C^ ^n*)
-38-
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
(4.11)
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
-39-
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
equation,
W3(2 jVk ^- (4.12)
integrate over yo and employ (4.9) so that we have
7. '.O (4.13)
But,
(4.14)
so that upon using this in (4.13) we arrive at Smoluchowski's equation
in the form
TL_ OQt t 7= (4,.15^ ^ *)
-40-
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
equation
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.
--b
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)
dj;
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)
-41-
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
-42-
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
(4.22)
and we have interchanged the order of integration. We now develop R(y)
in a Taylor series in (y x),
(4.23)
atid- aid
-43-
Upon using (4.23), (4.21) and the assumptions concerning the higher
moments, the right-hand side of (4.22) becomes
JFC^)O^
C4)P ^
(4.24)
+ 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
(4.25)
,ZN ZN ,Z la ~ i
n Y\?t +01~
b I- V.0 tV
^M-l 2..
0- oA-
~
x ane
+t r ^
-44-
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
(4.28)
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
-45-
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)
Also
-46-
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
-'A
d At 0
.... < I A = 0
(4.36)
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
(4.37)
(4.38)
L
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- ~ ~
-47-
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
(4.39)
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,
(4.40)
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
-48-
W, C=e~B crsA+ +~fS-1\N+1A
I m= -
(4.42)
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.
CHAPTER V
APPLICATION OF THE FOKKER-PLANCK EQUATION TO NONLINEAR
ELASTIC SYSTEMS
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
i,
where
K f f(& ++r 4' ^if-Wfc (5.2)
S
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
-49-
-50-
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
NW
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
conditions.
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
-51-
of the infinite series is possible. Thus,
(5.6)
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
assume.
-52-
If we'now substitute equations (5.3), (5.4) and (5.6) into equation
(5.2) and perform the indicated integration, we obtain
(5.8)
where V is the strain energy Vs after integration and
b-w- c)n= j(C U)-iS (5.9)
S
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
-53-
h"" I
^ w ~ ~ p~_ _
-V
(5.11)
xJ- ph /3~n t~
I
(yk
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
C*)
(5.13)
av
0 &
VIVT%_
~Tm ~
-54-
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
(5.14)
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
-55-
where,
Al*-0 A
(5.16)
mn- A~bl-0 o X
Equations (5.14) imply that
'e'- ,i^^c;,y
with
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 < >
(5.19)
-56-
ALr
-tn )I~-
4tk
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
(5.20)
J M
c^ Tftn.
-57-
dUr -- wn.
--4-
Ip,'_
CLvt7L r1,5
J 1-
m^n-vn'~s
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)
^m5r.M
--a f S ,), ,
a^J^1--Q~~'b11-^
I' V
V V.
'V% bC
9e n W
+ a
6wr^v%
A-t -p0
- 12"'eS &
S (
3V
(~'h~L
plyvr r
FL V%
+ ~3,
av,,
(,-M V,.) -
-58-
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)
-Z--
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
-59-
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
1NL FOL3
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)
-60-
The eigenfunctions of the linear problem are sin mrx so that we
L
write
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
where
E- NoL T'
(5.31)
N,- N=L
L.
8)
9)
0)
-61-
Here 6G is the mean square deflection of the first mode of the linear
so
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
2
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.
-62-
where
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.
-63-
2
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.
2
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)
-64-
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. "
-65-
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
N2
The kinetic potential of a plate undergoing large vibrations is
given by equation (2.28). The strain energy is therefore given by
(5.43)
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 -*
-66-
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- "
(5.45)
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 ) +
where
i-rfvwfYpr .Wv
[f,[k r 2*l
,- V%, 15 Ft E I
K nItpos p F= mpt- t TV
-67-
0.
E = fcA --
b
F =. 1 = CMA'r
LQ.
A -l cra mn Y
Ampri&t-. M7 --
0Oo
P" 21- rn rIT d"34
aC.UL a.r
b 4s <
AAML Q'w Y
-r--,,
AIA, rnlt.
Awn -nA
OL
c C0a4 I L
OL-- o-
(5.46)
b
0 Qb b Va \^ ,,
104 I CsAu
ac
a.
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
.KSI
co& w
-68-
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
relation
(5.47)
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
(5.48)
cf1 AI F-. i- x"
C K [ 4
L b
(5.49)
*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. ~
nnpClo~e\l
;L L*
- ZEE
& S
7 -
-00 '- *ae brldr:
-69-
If we integrate out the vmn, again making use of (5.47), we finally
obtain
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%
(5.50)
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
oO
N -00
N0 = N'/ab
jl0 i^O
-70-
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
where
--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.
CHAPTER VI
NUMERICAL INVESTIGATION
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
-71-
-72-
4-1
0
1-I
14.
0C
a4
,0
o4
-4
ow
0
o 0
44
0
(3
cun
44>
P14 ,
-73-
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
reasonable.
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.,
-74-
- c;o CO
*r4
o,
0
4-1
0
-1
0
rFl
0
0 0
0
4J
:J3
a0
O a
i"
-75-
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
2
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
coefficient.
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
q-2.
C~~ II z ir ,6
(6.2)
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
-76-
+ U .
L -
-IL
__ ( xyWn%%i
0- +-i. t- Pw),b-rn
-- l
+2..
- -b
b14 CL"
1 L + I n
- C (*>**'' ^ T1
Z. o.b
+1- 2.
6't a b~
9a-.b
6c4b*
b"Cb''
L =.n -
CeL
+b
8b
801
-j- -
aaO
(6.4)
C. I
3tZIn -h g- + +
8b..> a ^ 0 'g.',Q
(
1X
2.
0. Z.b
(6.3)
Now,
+~v-)i + "L ,- %w ...
-. z. b
---- ---
+ &
W1 C,) =
CRILf
EhrrY hL
_gEhrr' I Y8
2NoC\- ri
11
-77-
and all other Lnllll = Kmnllll = 0 so that (6.3) reduces to
\N, C =C PA
r -
PO
I[c+l1
1,g
Wit E
+ ~lY8X
Ia b
o,. I=>
I..+ 2- ,
Sb
L +
L_ CL
(6.5)
where
L
0O
I of = b/a
This is of the form
^^}- ^^C
(6.6)
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
-CO
(6.7)
4-
vt% tII
~
(I, Cy\~L
o- Lt,
-78-
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
2
T as determined from numerical integration of (6.7) is plotted
2
against "2 for 1) = 1/3, h = 1, and two difference aspect ratios
Po
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.
-79-
J 0 Co o o
c 0 0
o 6 6
l-
4,
4-I
0
4-
cd
r.-l
0.
4-i
o
0)
go
4-1
0
U
-X4
0
0U
1.a
,
Fr
a '
3
04
CO C
gl
0
0
o
o
CHAPTER VII
CONCLUSIONS
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.
-80-
-81-
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.
LIST OF REFERENCES
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.
-82-
-83-
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.
-84-
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.
BIOGRAPHICAL SKETCH
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
Fellowship.
Richard Edgar Herbert is a member of Sigma Xi and Omega Delta Phi.
-85-
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
SUPERVISORY COMMIITTEE
Chairman
4t / 7, .-- -
-. a,.sJ_ .[,
-- --,, --'- -.
!^.^.C
UNIVERSITY OF FLORIDA
3 1262 08553 79091111111111111 I
3 1262 08553 7909