POLY-WRI 1569-89

Review of Probabillstic Analysis of

Dynamic Response of Systems

with Random Parameters //J - 2 _

by

F. Kozin and J.M. Klosner ,P_qT(_/

Final :Report for

NASA Grant No. NAG- 1- 986

NASA Langley Research Center

Polytechnic University

December 1989
Table of Contents

I,

II.

III.

IV.

V.

VI°

VII.

Section

Introduction

L iouville Equation

Perturbation Methods

Mean Square Approximate Systems

Statistical Linearization of Nonlinear Systems

1. Systems with Known Linearity

2. Systems with Stochastic Parameters

Closure

References

Page

1

4

12

19

26

26

28

31

32
-1-

I. Int-roduction

The objective of this report is to review the various methods that have been

studied in the past to allow probabillstic analysis of dynamic response for

systems with random parameters. In general, the mechanical parameters (i.e.,

spring, damping, Joint parameters, dead zones, etc.) may not be known exactly.

If, for example, the variations about the nominal values are very small, then

the dynamic response would be adequately obtained deterministlcally.

However, for space structures which require precise pointing, it appears that the

variations or uncertainties about the nominal values of the structural details and

of the environmental conditions may be too large to be considered as

negligible.

Thus, these uncertainties must be accounted for on some rational basis

which we shall assume to be defined in terms of probability distributions about

their nominal values. The quantities of concern for describing the response of

the structure includes displacements and velocities, as well as the distributions

of natural frequencies. The exact statistical characterization of the response

would yield Joint probability distributions for the response variables. Since the

random quantities will appear as coefficients, determining the exact

distributions will be difficult at best. Thus, certain approximations will have to

be made. There are a number of techniques that we shall discuss that are

available even in the non-linear case.

In the general case, the n-mass linear structural system possesses the

dynamical description through the differential equation

My + C_r + Ky = f(t) (1.1)

where M, C, K are the mass, damping and stiffness matrices, and f(t) denotes

the external excitation on the structure.
-2-

We shall define the vector equations (1.1) through the vector

X

which represents the system (1.1) as

= Ax + Bf(t), x o =

where

i [o_M__K _M_ICi} , B= {0} M- 1

The solution of (1.3) for Xo(0)=x o ls

We assume

deterministic.

t

x(t) = e-ZtXo + f e-Z(_-_)Bf(r)dr

o

that the matrix B, determined by the mass

Thus, the random quantities appear in matrix

(1.2)

(1.3)

(1.4)

(1.5)

constants, is

A. We shall

denote the external the random variable in A, as X I.... ,X m. Further,

excitations f(t) may or may not be random.

The most important topic for engineering systems is how uncertain

parameter values influence the accuracy of system response prediction. It often

suffices to know how these uncertainities influence the accuracy in estimating

the values of the national frequencies and their corresponding normal modes of

motion in a conservative system (C= 0).

Since linear system response prediction depends upon frequency response or

impulsive admittance, our interest will center on natural frequencies, normal

modes, frequency response, as well as impulse response. The methods of

techniques that we will describe in order to pursue the various subjects are: (1)

Liouville's equation; (2) perturbation methods; (3) mean square approximate
-3-

systems; and (4) non-linear systems, with approximation by linear systems.
-4-

II. Liouville Equation

In this section we derive the method based upon the LIouvllle equation for

the time evolution of the Joint probability distribution function of the state

space (2nx 1) column vector x and the system parameters.

The use of the Liouville equation In mechanics and statistical mechanics is

of long standing and goes back to Maxwell (see for example [1,2,3]). These

references do not consider random system parameters, and average quantities

under equilibrium conditions is of main Interest. While not of direct interest to

us, it is possible to adapt these early methods to our needs. We derive the

needed form of the Liouville equation following a procedure suggested by

Kozin [4] for systems with random parameters and random or deterministic

initial values.

We are interested in the linear equations of motion in the form given by

(1.3) with f= 0:

= Ax , (2.1)

where x is the (2nxl) colwumn vector whose transpose x T has the form x T =

{X1 ..... Xn;X 1 ..... ]in} and A is the (2nx2n) matrix glven bY the first of

(1.4) The vector X Is the state space form for representing the system response;

the components of x will be denoted by Xk(t ), k=l ..... 2n. The random

variables in A are denoted by X 1..... X m. However, since the Liouville

approach applies to general nonlinear as well as linear equations, we consider

the general system

xk = gk(xl, .... X2n; X t..... Xm;t) , k=l ..... 2n (2.2)

Let p(x 1..... X2n;X 1..... Xn;t ) be the Joint probability distribution of the

random quantities ix, ..... X2n; X, ..... Xn). We define the characteristic

function ¢ as
-5-

¢ =E expi _ 8kXk(t ) + j ,

k--_l j

The differentiation of (2.3) with respect to time gives

i=vc? (2.3)

2n (

CO--'_¢= E i E 8k:Xk(t) expl

0t k_--1

The use of (2.2) in (2.4) yields

E 0kxk(t) + ECrXj m i].

k_l i

(2.4)

Since (2.3)

CO--'_'¢ = 1 E 8kE ;k expI E 8kXk(t) + ¢]Xj (2.5)

cot

k=l k=l i--

is essentially the Fourier transform of the Joint density function

p(x 1..... X2n, X 1.... Xm; t), the Inverse Fourier transform of (2.5) produces

The

cop 2n co(gip )

co-?=- E axj (2.6) j=l

solution of (2.6) for p is given by a suitable function of the

independent integrals of the Lagranglan system

dt - dp dXl

1 Icogl + + cog2n I

gl

Let u I ..... U2n be 2n-independent Integrals of (2.7).

dX2n

(2.7)

g2n

Then we know that the

general solution of (2.6) is

p (X 1 ..... Xin; X 1 ..... Xm;t } = h {u 1 .... ,U2n; XI,... ,Xm;t ) (2.8)

where h is an arbitrary function whose form is determined by the initial

conditions on x.

In particular, consider the case where X 1..... X m are

mechanical parameters which are independent of the Integrals.

be written as

explicit random

Then (2.8) can

p(x 1..... X2n; X 1..... Xm;t ) = hl(Ul, .... u2n; t)h2(X 1..... Xm) , (2.9)

where h 2 is the joint density function of the parameters, and h 1 is the

conditional density of the randomness of the initial conditions. To illustrate the
_

form, consider as an example the simple linear system defined as

_'(t) + 02_x(t) = o ]

(2.10) X(to) = Xlo , _(to) = X2oJ

In this case, we would have (Xio,Xeo) random initial values, and w, a random

parameter. The Liouville equation is simply written as

8p

_P ----- --X _ID ._. 022X 1 (2.11)

cot 2 0x 1 COx2

where (2.7) is simply

dt dx I dx 2 dp

...... (2.12)

1 x 2 022Xl 0

In this case, we would easily find

P(Xl,X2;02;t ) = h 1 {Ul(Xlo,X20,t, to;02), u2(Xlo,x20,t,to;02))h2(02)

= h X2 °I XloCOS02(t- to) -+-_sin02(t- to)

02

- 02Xlosin w(t-,to) q- X2oCOS 02(t-- to)}h2(w)

Upon utilizing the fact that the initial values can be defined as

(2.13)

XIO

X 2

= xlcosw(t o- t) + msin02(t o- t) _--- Ul(Xl,X2,to,t,02)

O9

X2o = - 02XxSin02(t o- t) + x2cosw(t o- t) _ u2(xl,x2,to,t,02)

the final probability form for the simple oscillator (2.10) becomes

(2.14)

p(xl,x2;w;t ) = hl{ul(xl,X2,to,t;02), u2(x1,x2,to,t;02)}h2(02 ) (2.15)

For the higher order system, the exact form of p in (2.8) would be obtained as

in (2.15).

We note that the initial conditions may be deterministic so that h I is a

product of impulses at the origin and at unity
-7-

hl(u 1..... u2n;t) = _(Ul)...5(Un)_(Un+ 1- 1)..._(U2n)

where 5(.) is the delta function. Thus (2.15) would become

(2.16)

P(XI ..... X2n; Xl ..... Xm; t)=S(Ul)...5(Un)S(un+l- 1)...5(U2n)h2(X1 ..... X m) •

(2.17)

Let us illustrate thls process by the simplification of the example (2.10)-(2.15).

Consider again the undamped one degree of freedom linear oscillator. Let w 2

be a random parametric value. For to=0,

X2

U 1 = XlCOSCJt- --slnwt ,

_d

u 2 = wx 1 slnwt + x 2 coswt

Assume that w 2 has a dlscrete distribution given by

(2.1s)

m

h2(w2) = _ piS(W2_ %2)

1

For xl=0, x2=l at t= 0, (2.17) becomes

(2.19)

in

P(Xl,X2,W2; t) = 5(ul)5(U 2- 1)_]PiS(w 2- w_)

1

in

+ x2coswt- 1}EPi_(_ 2- _g)

1

(2.20)

Let us determlne the mean of x 1 to illustrate a possible use for (2.20); we have

F "1

Straight integration of (2.21) (see [7]) yields from (2.20)

(2.21)

m sinwit

E{Xl} = EPi1 Wi

Other illustrations, including damping, are given in [4,5].

(2.22)

We are frequently

concerned with the moments of x. Let us show how (2.6) can be employed to

obtain them.

To keep the details simple, consider the linear damped one degree of

freedom system with equation of motion,
-8-

X1 = X2 (= _1)

k

X2 = -- --Xl --

m

Equation (2.6) now takes the form

C

--X2 (= g2)

m

(2.23)

Assume the m, k,

Op 0(glP) c3(g2P)

-- + + = 0 (2.24)

0t 0x 1 0x 2

c are independent of the Initial vector, with h2(m,k,c) the

probability density function of these parameters, and write

p(xl,x2t;m,k,c ) = h I (Ul(Xl,X2,t,m,k,c), u_(Xl,Xvt, m,k,c)} h2(m,k,c)

In this case, we could simply rewrite the function h i as

(2.25)

hl(Ul(Xl,X2,t,m,k,c)) = pl(x,,x2,t[m,k,c) (2.26)

since the parameters (m,k,c) are conditional for h I. Upon inserting (2.26) into

(2.25) and then into (2.24), we obtain

where h_(m,k,c)

e x pe c ratio ns

0Pl 0(glPl) 0(g2Pl)

-- + + = 0 (2.27)

O_t 0X 1 C_ 2

has been factored out. Let us evaluate the conditional

E(x_ ]m,k,c} = ml,o(t) = f f xxpldxxdx2

E{x 2 [m,k,c} = mo.l(t) = f f x2p:dx:dx2

Differentiation (p{_r_ial) of these equations with respect to time produces

(2.2s)

0P 1

I111. 0 = f f x_--_- dXldX 2

But, from (2.23),

0Pl

= I f dXldX 

OP 1 0(glPl) O(g2Pl)

Ot Ox 1 Ox2

The substitution of (2.30) into (2.29), l_he employment of (2.23)

(2.2o )

(2.30)

and simple
-9-

integration by parts of the resulting terms on the right-hand side of (2.29)

finally yields

I_ll, 0 = mo, 1 ,

(2.31)

• k C

mo, 1 = - --ml0- --mo, 1 m m

The same procedure will produce the equations for the conditional moments

E{x_ Im,k,c}, E{x_ ]m,k,c}, etc. We note that for the first conditional moments

we could have taken the conditional expectation of (2.23) to produce (2.31);

however, this procedure only applies to the first moments.

We integrate the moment equations (2.31) to obtain the conditional

moments as a function of time. On multiplying these moments by h2(m,k,c )

and integrating over m,k, and c, we obtain the moments of x 1 and x 2.

It is clear from the above discussion that the Liouville equation will provide

the exact solution for the Joint probability density function p(x 1..... x2a;

X 1..... Xm;t ) in the absence of external forces provided the integrals

u I ..... U2n can be obtained. Further, it provides a straightforward method for

determining the moments of x from which means and variances of x can be

obtained.

The Llouvllle equatlon applies when there are no external forces. We are

interested in the case when external forces are present, of course. Let us see

what can be done along these lines.

The Fokker-Planck equation is the natural extension of the Liouville

equation (see [6,7]). We confine our attention to the case in which the

external force vector f can be obtained by passing gaussian white noise through

a stable linear damped system. We have as equations of motion, conditional on

M= m, K= k, and C= c,
- 10-

dx 1 = x2dt ,

k c

dx 2 =- --xldt- mx2dt + X3dt , (2.32)

m m m

dx 3 = - _x3dt + dB , x 3 = 0 at t=0

where we have employed the differential notation in this case, set f=x3, and

were dB is the Brownian motion increment with

E{dB} = 0 , E{(dB) 2} = a2dt . (2.33)

The last equation of (2.32) represents the fact that the excitation ls obtained by

passing a gaussian white noise through a linear first order stable filter. We

notice that for the Ito system (2.32) xT={xl,x2,x3} is a vector Markoff process

that generates a Fokker-Planck equation.

It can be shown that in this case the Fokker-Planck equation for the

conditional probability density function p, is

I

....0t 0x, (X_pl) 0x 2 -_-X * _ Rm x2 + _-" s P, 0x 3 (- _x3Pl)

a _ 02Pl

+

2 Oxg

(2.34)

We observe that all but the last term on the right are the same as would have

occurred in the Liouville equation in the absence of f. Let the conditional

moments be

mkl,k2,k8 = E{Xl kl x 2k2 x: _} (2.35)

k2 k3

= f f fx k' xe x3 pl(xl,x_,x3) dxldx2dx3

Then, proceeding as in the development of (2.31), we find
11-

ml,o,o _ mo, l,O

°-c m o i o+ 1 m,o i(2.36)

r}10,0,1 _-_ _ _mo,0,1

On multiplying the solutions of (2.36) by h2(m,k,c ) and integrating out the

condition in these three conditional moments, we finally obtain the moments of

x as a function of time. We obtain in analogous fashion the differential

equations for the second conditional moments; we do not do this as the steps

are of a mechanical nature and not of direct interest. The maln point to notice

is that differential equations for the conditional moments of x can be obtained

when an external force is present in the equations of motion provided this force

is produced by passing white noise through a suitable filter.

It ls important to point out that for any gausslan external excitation the

solution vector is gaussian conditioned on the random parameters. Therefore,

all conditional moments can be obtained but not as easily as above [8].

The Llouvllle equation enabled us to obtain, in a straightforward manner,

the exact expression for the conditional probability density function p.

Reference to (2.54) suggests that it will be much more difficult to obtain Pl

from this equation and we shall not pursue this line of thought further.
- 12-

III. Perturbation Methods

The references [9-17] address the eigenvalue (natural frequency) and

eigenvector (normal mode) problem in structural systems by perturbation

methods. Before discussing methods or techniques involved, it is important to

understand at the outset that the geometry of the structure, how its equations

of motion are assembled, the final mathematical form of the equations of

motion, and how randomness In parameters Is introduced have a profound

influence on the nature of the results obtained.

A structure's geometry can be in the form of a linear array (chain) of

elements that may, for example, consist of simple harmonic oscillators strung

together in a line, beam segments continuously connected at a sequence of

supports in a line, etc. The geometry Is the simplest possible in such

arrangements. Plate or shell type structures have a two-dimensional grid-llke

geometry and are next In order of complexity. Finally, we have the general

case in which one or two-dimensional geometries are interconnected in a

complex manner.

The equations of motion depend on the coordinate choice, particularly when

the fact that mass is always distributed is taken into account. Reference [18]

discusses methods of making this choice and illustrates the substantial

difference in response that can occur due to different choices. Reference [16]

also discusses a component mode synthesis method for selecting coordinates

and assembling the equations of motion. A coordinate transformation of the

equations of motion is sometimes employed as in [10,19] and the altered form

of the equations may be advantageous.

Let us briefly present here a typical perturbation procedure. We consider

the free motion of a conservative system governed by the equations
- 13-

I_;+ Ky=0

where y Is the (nx 1) column vector with transpose

yT = {Yl..... Yn}

(3.1)

(3.2)

Now the elements In the symmetric stiffness matrix K are determined by the

bars, beams, columns, joints, etc., making up the structure, and the

uncertainties in the structure reside In these elements. Let there be m

structural elements, and let the stiffness matrix of the im structural element be

K i= (I +Xt)_ , i = I ..... m , (3.3)

which produces the (nx n) random stiffnessmatrix

In

K = E Ki = {Kjk} (3.4)

I

The random variables X I ..... X m are regarded as small perturbation terms

describing the uncertainty present in the structural elements and we assume

E {xi}= o , V_r {Xj = o_2 , (3.5)

K__ is the mean stiffness matrix of the Ith element, K=K T, I.e., K is symmetric

in the Elk, Elk lS the random stiffness element corresponding to Yl and Yk, and

we assume masses of the elements do not change. We note also that we can

write (3.4) as

which gives also

K= K-t- ZXiK i , K = EK I (3.6)

0K

E {K}= K and _ = K__. t , (3.7)

0X_

where K is the stiffness matrix of the structure with each member taking its

mean stiffness.

Assume normal mode motion

y = a cos(wt + ¢)

with c_ the (nx l) column vector defined by

(3.s)
- 14-

C_T __-- {0_1,... , O_n} •

Then, substituting (3.8) into (3.1), we obtain

(K-co2I)_ = 0 ,

where again I is the (nx n) unlt matrix.

(3.0)

(3.1o)

The squared natural frequencies cot2 are determined by the n roots of the

equation

det _K- ¢o2II=o , (3.11)

revealing that the C0r2 and cor are random variables since K contains random

variables. Let the random mode corresponding to w r be the (nxl) column

vector o_r. Then we can write

(K-Wr2I)C_r = 0 , (3.12)

wlth the usual orthogonality relations

c_TIc_s = 0 , _TK_ s = 0 if S _ r

(3.13)

o_TIOer = 1 , c_TKO_r = COr2 ,

where "T" denotes transpose, as before.

We are now interested in expressing the random variables w r and c_r In

terms of a power series in the random variables X i. Consider, for example, the

r _ natural frequency w r of the system expressed in the form

m m m

% = + XiX, + X,pC,Xj+ .... (3.14)

i-_-I i-_-I j_l

where w__r represents the r th natural frequency of the mean system, and the

Xi,XI] .... are to be determined. Once we know the Xi,Xij ..... we can obtain

statistical properties of °Jr or any other quantity of interest. Let us consider a

general formulation of this problem, considering natural frequencies and

normal modes. We follow the method suggested by Zarghame [14] which

appears well suited to computation.
- 15-

Differentiate (3.12) with respect to Xi:

0w r ) Oar _- 2Wr'_l I ar+(K-Wr 2!) -_ii ---0

Next premultlply (3.15) by a T obtaining

(3.1s)

oqw r i

IJ a r = 0

since by the symmetry of K(K=I< T) and (3.12)

(3.16)

aT(K- Wr2I) =0

Thus, with the third of (3.13)

0% 1 arT KLar

= 20j--7

This is to be evaluated at X 1.... =Xn=0 (I.e., X= 0); we obtain

(3.1_)

{a r] = TK, , (3.1S) 1_2_

°_Xi Jo 22 --

where the underb&rred quantities are to be evaluated for the system with mean

stiffness. We note that (3.18) gives the sensitivity coefficients [20,21] of co r

with respect to the X i. The importance of the sensitivity coefficients resides In

the fact that they reveal by their magnitudes those w r that are either sensitive

or insensitive to uncertainity in members values.

The a k, k=l ...... n span the coordinate space; hence, we may write

c_a r

3X:i

We substitute (3.19) into (3.15):

-- = E £(i)% (3.1g)

000) r )

I___{i - 2w r -_iiI a r -{- (K- ¢dr2I) E _r(iJ)aJ = o

J

Now premultiply (3.20) by oekT, obtaining

0w r ) a T K:,- 2w r -_iI a r + a T (K-Wr2I} E]_(i j)a)= 0 ,

l

(3.20)
16-

or, for k_Ar and with the use of (3.13)

Differentiating the next to last of (3.13) with respect to X i gives

G0_ r

_r_,-_-7=0 ,

which on premultiplylng (3.19) by c_TI then demonstrates that _i r) = 0. Thus,

a T K__4ar

Z[k_= , k_r, (3.21)

_- _r_

and hence from (3.19)

where the prime on y_/

0o_ r o_T _,iOt r

= E' _ , (3.22) _r_-_

means that k does not take the value r. When

evaluated at X= O, we have

[a_rI _-_, _?m_

_i lo j ,_2.r 2 - 0)i2 _ (3.23)

where agaln the underbarred quantities are evaluated when members take on

their mean stiffnesses. We note that (3.23) gives the sensitivity coefficients of

the mode shapes with respect to Xj. Without showing the detailed deviation, we

simply state that it can be shown that

_%3r = 1_.. f_(k)(4(r).j_(4(k)_(r)) (021 _ _ _._.¢2)_ 2 (3.24)

G_XiOXj 0 2 t_r| /'_kj m_--rj _'ki [ 0_i J 0 _ 0

and

where

020tr } = _ _(k) _r, ij _ • (3.25)

OXtOXj o k
- 17-

_-_ Lax_Jo [ _XjJo t

rOO rI [0%/

oi' k_r

(3.26)

_!__=- loxjjo _ IoX,/o

Summarizing our results up to this point we have for the random variable

°Jr----°3r+Y] {, (0Xi,oXi-{- 1 i_? 0X_j} XiXJ+o ....

where the partial derivatives are supplied by (3.18) and (3.24).

for the random variable

(3.27)

(3.2s)

We also have

_r -- O_r -+- _ [ O_:i ] Xi "+" "_-1Y]_-'_ 1 2°r} o_2C|o_X:j XiXj + .... (3.29) i o i j o

where (3.23) supplies the first partial derivative, and (3.26) and (3.27) supply

the derivatives in the double sum.

Let us now consider the statistics of w r, etc. Consider Eq. (3.28) first. We

have, on taking expectation,

E{Wr} :W--r + Y 0Xi0Xj o

Even if the X's are independent E{Wr}_A__ _, since the E{X_}_A0 terms are still

present. Now square (3.28) and take expectation

E{_r_}= _+___ [o_Xi0Xj 0% } E {XiXj }

z i ] o

+ i_j_" IOq°Jr/ la_:l E{X,Xj} [ Ox,/o[ OxJ Jo

l O°Jr [ 02¢_r }

+2-- _, ,k I-_-T}o OXjOXk o E{X'XrXk}

02Wr 02Wr

+ 5]Y1,5]_ [

OXkOX_/ i j k _ t o 0

We can now approximate Var Wr; it Is defined as

(3.31)
- 18-

Var Wr2 = E{Wr 2} - [E{Wr}] 2 . (3.32)

Thus, it is a straightforward task to approximate the first two moments of w r.

If we extend (3.28) to cubic, quartrlc .... terms In the X i, then (3.30) and

(3.31) would contain additional terms. How far we should continue this

process will depend on the relative slze of the terms containing E{XiXj},

E{XiXiXk}, etc. and what information we have that would enable us to evaluate

these expectations. It is not usual that we can evaluate any more than E{X_Xj}.

We note that Zarghame's method described above gives sensitivity

coefficients for natural frequencies and corresponding normal modes plus series

expansions for these quantites In terms of the random variables X 1, .... X m

that define the uncertainty present in the stiffness matrix K. Moments of the

quantities are easily obtained, but it Is practically Impossible to obtain

distributions for the natural frequencies and corresponding normal modes. For

confidence interval location and size for a natural frequency, for example, we

must approximate using

E{Wr} 4- 3 _ _te (3.33)

as a rough indication of a 99°/o confidence Interval. This interval gives us some

Idea of the spread in a natural frequency and it could be employed to make

reasonably sure that no steady excitation frequencies were contained therein for

all w r. Alternatively, we might employ the signal to noise ratio

E{Wr} (3.34)

Vv/V r

to obtain an idea of how important stiffness uncertainty is for natural

frequency; if (3.34) is greater than 20 or 30 say, we would regard the location

of °Jr as deterministic; on the other hand, if (3.34) Is less than 5-10, it might be

unwise to ignore this level of variability in the location of Wr, depending, of

course, on the consquences of such uncertainty.
- 19-

IV. Mean Square Approximate Systems

We consider, in this section, a technique for Including disorder or parameter

uncertainty that follows a different line than taken In previous sections.

Specifically, mean square systems are employed. We begin by introducing this

type of system [22,23].

Let us begin with a very simple example in which there is no disorder and

no damping. Let the coordinates ql ..... qn" Then, with

2T = mjkCljClk, 2V = kjkqjq k, _W = fj(t)6qj , (4.1)

where summation is on multiple subscripts. Then, with mass coefficients

included, (4.1) can be rewritten as

mjk_ k + kjkq k = fj(t)

Let, with foj constant,

Then, the forced motion

satisfies

%(t) = fojCOS( t + ¢)

qk = ukcos( wt + ¢)

(kjk -- w2mjk) uk = fo]

(4.2)

(4.3)

(4.4)

Suppose we pick an w and a set of u k which may not be one of the natural

frequencies and normal modes. Then the right of (4.5) will not be zero and we

need force amplitudes Ej to produce this motion:

These equations state that given the foj and _o, the u k are determined by the

solution of this linear system of equations. Further, if w is the natural

frequency w r and the u k define the r th mode shape _rk, then the foj must

vanish. Let us look at the natural frequency problem In an unorthodox

m anne r.

(4.5)
- 20-

(kjk- _2mjk ) Uk = _j (4.8)

The cj are the amplltudes required to maintain the assumed motion; we regard

the Ej as the amplitudes of the constraint forces required to produce the

motion.

Consider next a summation of second order amplitudes

n

I(n,w) = _ e_ > 0 . (4.7)

l

For a fixed w, this is a positive definitive quadratic function of the u's. We can

use this equation to find the natural frequencies w r and corresponding normal

modes O_rk. Assume the u% are normalized in some manner (for example,

Un=l, or better mjkujnk=l). For fixed w, we find the minimum of I(u,w)>O.

Notice that if w=w r the u% that produce a minimum are the O_rk and

I(ark,Wr_)=O, since the El=O, J=l, . . . ,n, in this case. It follows that if for a

specified w we find the minimum of I(u,w 2) and this minimum equals zero,

then this w is a natural frequency and the u that produces this zero minimum

are proportional to the corresponding normal mode. Let us consider another

interesting aspect of this method.

Conslder a frequency window g(w) with the following properties:

g(w) > 0 , oa'< oa < oa" ,

(4.8)

0211 W ;t

f g(w)dw = 1 , f W2g(w)da_ < co

¢d/ Ol

Replace (4.7) with

W II

I(u,g) = f (4.9)

Find the the u that makes (4.9) a minimum. The interesting feature of this

method is that if there is a natural frequency of the system in the frequency

interval (w',w"), the u in mlnI(u,g) determine the normal mode of this natural
- :_i j :'_ : ,: > .... _: .... [.._..-::_,::.,:-'_ T_.;:._:::: ,

-21-

frequency. Let these u be In component form {U_ r), .... un(r)};then the

corresponding natural frequency isdetermined by the Rayleigh quotient:

col: kiku r)u r)

mjkui(r)u r), (4.10)

where we assume we have found the r t_ normal mode and its natural

frequency. It follows that if there is concern that an Interval (co',co") contains a

natural frequency, we have the method for determining if this Is the case

without determining all natural rrequencles of the system. References [22,23]

glve details on this matter we will not discuss In this report.

The computational problem of finding the minlmum of I(u,w _) is carried

out using one of a number of computer codes based upon conjugate gradient

techniques, and, hence, Is not a problem.

So far, there has been no disorder in our system; i.e., the parameters mjk

and kjk have been assumed to take definite values. Let us assume at this point

that mass and stiffness contain random variables. We define, In this case,

02H02 ]

I(u,g) = E f g(a) , (4.11)

l

where, in vector-matrix form

n

_ = uT(K- w_2M)T(K-wZM)u (4.12)

1

n

Since E only operates on _-]_] 2 In (4.11), we have

1

E e = E uT(K-CO._2M)T(K-CO_.2M)u , (4.13)

and co Is a fixed parameter In (4.13). We assume the u are parameters to be

determined. Thus, (4.13) takes the form
- 22 -

E {n} _6i 2 = uTE /(K-co2M)T(K-co2M) }u (4.14)

We note that (K- co2M)T=K- w_2M because of the symmetry assumed for K and

M. In all events, means and second moments of K and M are all the

information needed to determine the expectation in (4.14).

We than proceed as in the deterministic case, since I(u,g) has a

de _e rm lnistic form.

To relate (4.7) to (4.13), all we have to do Is assume

g(¢o) = 6(W-¢0) , (4.15)

where 6(.) is the delta function. The substitution of (4.15) into (4.11) yields

I(u,co) =uTE {(K-w2M)T(K- co2M)} u (4.16)

This expression differs from (4.7) because of the assumed random parameters

in K and M. If in (4.7), w Is a natural frequency of the deterministic system,

I(u,co)=0. The I(u,w)>0 In (4.16) because of the random parameters. Use of

this fact has been made In [32] to obtain an estimate of the variance of natural

frequency COr; the formula is

I(Ur,Wr 2)

Var cor = 4C°r2 , (4.17)

where cor is the rth natural frequency and u r is the corresponding normal mode

for the system with mean parameter values. Monte Carlo simulation [24]

reveals that (4.17) can be conservative and a correction Is suggested. Equation

(4.17) ls easy to use since a minimum for I Is not required. Further, (4.17)

provides a much simpler method for estimating the variance of cor than given in

Section III. However, mean square approximate systems do not provide any

information on E{COr} or on variability in mode shape. Let us next consider

how these systems apply to estimating frequency response with parameter
- 23 -

uncertainty present.

We take the equations of motion in the form

The frequency

external force,

reiwt f----Sj

with r fixed at (_jr--0 for jyAr, _rr=l, the component form of (4.18) is

qj = Z]_ l(w)e iwt

which is exact.

Mq + Ccl+ Kq = f (4.18)

response Z-X(w) and Z(co) are defined as follows. For the

Suppose we try to approximate (4.20) with

qj = _jr e|Wt ,

(4.19)

(4.20)

where the _jr are not known in advance.

satisfied and we must introduce constraint

satisfaction as

From

(4.21)

The equations of motion now are not

forces ej to bring about their

(Kjk- w2Mjk 4- lcOCjk } flkr- 6jr ---- _Sjr (4.22)

* (4.23) I(fl,cO) = E EjrEjr ,

t 1 P

where asterisk denotes complex conjugate. This I is just like (4.11) except the

fl have replaced the u. We find the fl that make (4.23) a minimum, denote

this fl by _. Then, _={_1 ..... _n} is the mean square approximate to the

Zj_l(w). It can be shown that if the system is deterministic (i.e., contains no

random parameters) the _ are exactly the Z_ l(w).

The ejr are complex, hence, the right of (4.23), when written out, is

j_E{[[KjR_CO2Mjk_IOJCjk}I_Ie_jr][(KjI_co2MjI-t-ioJCjt}t_Ir-_jr]} (4.24)
- 24-

It follows that the minimum of (4.23) is for the real and imaginary parts of flkr"

This added complication poses no additional computational problem [25,26].

The method also supplies an error criterion that makes it possible to Judge

the accuracy of the flkr"

References [25,26] describe In some detail how the above technique can be

applied to estimating the frequency response In a number of structures with

specific attention being paid to numerical details of the computations.

Reference [1] also describes how these techniques can be applied to the

construction of a sequence of approximants for a complex system by starting

from a highly constrained Initial system and gradually relaxing the constraints.

In these three references, extensive use is made of the error criterion to

determine when the estimated quantities (usually frequency response) are

sufficiently accurate for the purpose in hand. A comment on what mean square

approximate system provide is now in order.

We observe, for example, that these systems enable us to estimate

frequency response Zkrl(cd) by means of flkr(W). The flkr(Cd) are deterministic

numbers that take into account the means and Variances of the statistical

parameters of the structure. Thus, the flkr(W) provide a deterministic estimate

for Z_-rl(w). In the form given above and in [25,26,27], it is not possible to

obtain statistical information concerning the Z-t(w). However, it is possible to

employ Monte Carlo methods to obtain estimates for the flkr(W) given the

parameters are sample values to obtain sample values for the flkr(W) from

which statistical information can be obtained.

The statistical energy approach (SEA) merits mention at this point since it

also employs average energy concepts [28-32]. Basically, SEA estimates the

average flow of energy from one part of a structure to another. For example, if
- 25-

there is energy input into one part of a complex structure, this method provides

an estimate of how this energy flows into another part of the structure. Thus,

it is possible to estimate average vibrational energy present In any part of the

structure. Information of this type is frequently the only type of information it

is possible to obtain about the response in an extremely complex structure

containing a large number of undamped natural frequencies in 1 Hertz. Insofar

as we know, nothing has yet been done to include the influence of statistical

parameters; however, the work given in [31] suggests it might be possible to do

this.
- 26-

V. Statistical Linearization of Nonlinear Systems

1. Systems with Known Nonlinearity

In a number of real systems the non-llnearity property may be known. The

statistical properties of the response of such systems cannot, in general, be

determined. However, the statistical properties of the response of its linearlzed

version may result In an adequate approximation. To obtain the llnearlzed

form, the method of statistical linearizatlon has often been applied [33].

The basic idea is generated as follows. Consider the true system model

_r(t) + f(y(t)) = n(t) (5.1)

where n(t) is the random excitation vector, which is assumed to be known as

well. The Idea is to approximate (5.1) as the linear form

x(t) + Ax(t) = n(t) (5.2)

average of the mean of the difference terms in (5.3) is minimized, i.e.,

mln E {]]f(y(t))- Ay(t)ll _} (5.4) A

The solution Is simply obtained as

A = E {f(y)Y_} E {yy_}- 1 (5.5)

where ( )_ denotes transpose.

Therefore, A is defined in terms of the y-statistics. But, the y-statistics are

unknown! If they were known, we would not have to linearlze the model. The

traditional idea is to determine A via the linear X-system statistics as

A = E {f(x)x'} E {xx'} -1 (5.6)

The statistics In (5.6) can be exactly obtained in general, but when n ls

The approach is to write (5.1) as,

y(t) +Ay(t) + _(y(t)) - Ay(t)] = n(t) (5.3)

The linearization is obtained by choosing matrix A so that the ensemble
- 27-

Gausslan its determination defined by (5.2) becomes simpler. It should be

pointed out that if f(y) Is a vector polynomial then (5.5) can be determined

exactly as statistical values, through stationary solutions of various forms [34].

The following example for the Dufflng oscillator with white noise excitation

ls very typical.

Example ].- Statistical Iinearlzatlon for the Duffing oscillator with white noise

excitation.

Consider the second order system

y(t) + y(t) + y_(t) = B(t) (5.7)

The B process is the gaussian white noise, given through the classic Brownlan

process, which is Gausslan. The classic Ito form of (5.7) is written as

dy I ----y2dt

dy2 = -(y2+y_)dt + dB

(5.s)

From the conditions of B, we have

E{_B}: o ,

We wish to determine the linear form

which we can write as

E {(dB) 2} --- o'2dt

i_(t) + i(t) + kx(t) = fi(t)

(s.g)

(5.1o)

{

Clearly, k Is determined as

dx 1 = x2dt

dx_ = - (x2+kXl)dt + dB

k = E {x_"x_}E {x_}-': E {x_}E {x_}-1

However, since (x 1, x2) Is gausslan, then we find

k = 3E{x_}

For the system (5.11) we can determine the stationary density as

(5.11)

(5.12)

(5.13)
- 28-

p(xx,x2) = __rra 2v_ exp [ (x2+kx_))___ 1

This allows us to obtain the value of k and the x-moments as

(5.14)

k=l.22a; E{XlX2}= 0 ; E{x_}= 0.4080" ; E{x#}= 0.5000 -2 (5.15)

It is interesting to point out that for this particular example the moments of

(Yl,Y2) can also be obtained for the stationary case as

E{yly2} = 0 ; E{y_} = 0.4780- ; E{y_ 2} = 0.5000-

We see that the error in the linear form Is only in E{x_}, and is

2 (5.16)

E{y_}- E{x _}

E{y_}

This is not an unacceptable error.

0.070

=- or 14.4% (5.17)

0.478

2. Systems with Stochastic Parameters

The question that must be considered here is, if the system has not only

nonlinearities, but also possesses parameters that are random as well, in what

way may we apply statistical linearlzation so that the random constant

coefficients of the true system are reflected in the approximate system. We

illustrate the procedure by considering again the D uffing equation

+ cy + by 2 = B(t) (5.18)

where here b is assumed to be a random constant. The linear form is taken as

+ cx + abx = ]3(t) (5.19)

where a, the llnearlzatlon constant, ls determined from

minE{]Iby a_ o_by[ 2} (5.20)

Ot

which leads to

O_= E{b2y_} E{b2Y_} -'

In the terms of conditional moments (5.18) can be expressed as

(5.21)
- 20-

Following

-1

the traditional path, we determine c_ from the

(5.22)

linear x-system

statistics, so that

-1

-1

The conditional second moments are readily obtained as

(5.23)

2

(7

Eixlx=]b}--0 ; E{x_=Ib}= 2_bc

which, when substituted into (5.23), leads to

2

--. E{xgIb}= _

' 2C

(5.24)

O_2 =

when the mean of b, E{b}=_.

3 _2 1 3 o "2

2 e E{b} 2 eb

The x-moments can then be determined

(5.25)

For illustration, assume

constant b with mean _,

distribution function for the

(5.2o)

opening

p(b) =

{ i/2E_ 0

, _(1- e) <b<b'(l+c)

, otherwise

(5.27)

Th us

E0/b} = --=1 _n (1+C)

2cb (l-c)

where c=bandwidth parameters. Hence,

(5.2s)

E{x_}= 0.204(1/_c) 1/2 (a/c)_n (l+c)

(i-c)

The exact solution for this case Is

(5.20)
- 30-

E{y_}= 0.47S(1Fc)'/_(a Ic) [(I +_),I_ _ ( 1- _)'/_1

and the error is

(5.30)

E{y,_}- E{x_}

E{y_}

which for e=O is 14.5°7o.

0.478 [(l"Jr'£) 1/2- (1-- _)I/2] _ 0.204_n(1-t-c/1- e) (5.31)

o.47s[(1+c) _/:- (1- e),j2]
-31 -

V_ ° Cl(:_ure

The requirement for precision pointing of large scale structures, consisting

of trusses, radio reflectors and optical systems, has motivated the investigation

of how parametric uncertainties affect the system response characteristics.

Randomness of these large structures may arise from many sources, for

example: manufacturing processes; space assembly by humans, variations in

environmental conditions that may bear directly upon the material properties

and resulting mechanical behavior. These uncertalnitles must be accounted for

on some rational basis so that the quantities of concern for describing the

response of the structure can be statistically characterized. Such quantities

include: natural frequencies, normal modes, frequency response, etc. Specific

random techniques that are available and discussed in this report are:

Llouvllle's equation, perturbation methods, mean-square approximate systems,

and statistical llnearlzatlon.

These are the techniques that we shall conslder in order to develop

procedures for determining the response of space structures. It Is conceivable,

however, that advances in these techniques will have to be developed as our

research effort unfolds.
- 32 -

"VII. References

1. J.L. Bogdanoff and F. Kozln, "Simple .4pproxlmants for Complex Linear

Systems," .41.4.4 26th Structures, Structural Dynamics and Materials

Conference, Orlando, FL, April 15-17, 1985.

2. T. Igusa and .4. Der Kureghain, "Stochastic Response of Systems with

Uncertain Properties," Trans., 8th Inter. Conf. Structural Mechanics in

Reactor Technology, Brussels, Aug. 19-23, 1985, M14/5, pp. 117-122.

3. ,4.. Chang and J. Rissanen, "Regulation of Incompletely Identified Linear

Systems," SIAM Journ. Cont., Vol. 6, No. 2, 1968, pp. 327-347.

4. F. Kozln, "On the Probability Densities of the Output of Some Random

Systems," Journ..4ppl. Mech., Vol. 28, 1961, pp. 161-169.

5. P.F. Chenea and J.L. Bogdanoff, '_Impedence of Some Disordered

Systems," Colloquium Mechanical Impedance Methods for Mechanical

Vibrations, .4SME, 1958, pp. 125-128.

6. L. Brillouln, "Inevitable Experimental Errors, Determinism, and

Information Theory," Information and Control, Vol. 2, 1959, pp. 45-63.

7. J.L. Bogdanoff and P.F. Chenea, 'T)ynamlcs of Some Disordered Systems,"

Inter. Journ..Mech. Sci., Vol. 3, 1961, pp. 157-169.

8. T.T. Soong, Random Differential Equations in Science and Engineering, .4cad.

Press, 1973, Sec. 8.3.

9. Lord Rayleigh, Theory of Sound, 2rid Ed., Macmilllan, London 1894 (also

Dover), Section 90.

10. W.H. Wittricks, ''Rates of Change of Elgenvalues with Reference to

Buckling and Vibration Problems," Journ. Roy..4ero. Soc., Vol. 66, No. 9,

1962, pp. 590-591.
- 33-

11. T.T. Soong and J.L. Bogdanoff, "On the Natural Frequencies of a

Disordered Linear Chain or N-Degrees of Freedom," Inter. Journ. Mech.

Sci., Vol. 5, 1963, pp. 237-265.

12. R.L. Fox and M.P. Kapoor, ''Rates of Change of Eigenvalues and

Eigenvectors," AIAA Journ., Vol. 6, No. 12, 1968, pp. 2426-2429.

13. J.D. Collins and W.T. Thomson, "The Eigenvalue Problem for Structural

Systems with Statistical Properties, AIAA Journ., Vol. 7, 1969, pp. 642-648.

14. M.S. Zarghame, Personal Communication, 16 Nov. 1970, Dept. Math.

Arya-Mehr University of Technology, Tehran, Iran.

15. L.A. I(Aefling, "Comment on The Elgenvalue Problem for Structural

Systems with Statistical Properties," AIAA Journ., Vol. 8, No. 7, 1970, pp.

1371-1372.

16. T.K. Hasselman and G.C. Hart, "Modal Analysis of Random Structural

Systems," Proc. ASCE, Journ. Engng. Mech. Division, Vol. 3, No. 6, 1972,

pp. 561-579.

17. M. Shinozuka and C.J. Astill, "Random Eigenvalue Problems in Structural

Analysis," AIAA Journ., Vol. 10, No. 4, 1972, pp. 456-462.

18. M. Shinozuka and J.N. Yang, "Random Vibration of Linear Systems," Inter.

Journ. Solids Structures, Vol. 5, 1969, pp. 1005-1036.

19. A.B. Ku, '_Upper and Lower Bound Eigenvalues of a Conservative Dlscrete

System," Journ. Sound Vib., Vol. 53, No. 2, 1977, pp. 183-187.

20. D.C. Reddy, "Eigenfunctlon Sensitivity and the Parameter Variation

Problem," Electronlchetters, Vol. 3, 1967, pp. 205-212.

21. R.C. Reddy, "Eigenfunction Sensitivity and the Parameter Variation

Problem," Inter. Journ. Cont., Vol. 9, No. 5, 1969, pp. 561-568.
- 34 -

22. J.L. Bogdanoff, "Mean Square Approximate Systems and their Application

in Estimating Response in Complex-Disordered Linear Systems," Journ.

Acoust. Soc. Amer., Vol. 38, No. 2, 1965, pp. 244-252.

23. J.L. Bogdanoff, "Uniqueness of Mean Square Approximate Systems," Journ.

Acoust. Soc. Amer., Vol. 39, No. 2, 1966, pp. 261-268.

24. A.J. Schiff and J.L. Bo_danoff, "An Estimator for the Standard D eviation of

a Natural Frequency," Part 1 and 2, Journ. Appl. Mech., Vol. 39, 1972, pp.

535-544.

25. t(.W. Kayser and J.L. Bogdanoff, "A New Method for Predicting Response

in Complex Linear Systems," Journ. Sound Vib., Vol. 38, No. 3, 1975, pp.

373-389.

26. J.L. Bogdanoff, K.W. Kayser and W. Krleger, "A New

Predicting Response in Complex Linear Systems II," Journ.

Vol. 53, No. 4, 1977, pp. 459-469.

27. T.T. Soong and F.A. CozzarelIt, '_Vibration of Disordered Structural

Systems," The Shock and Vibration Digest, Vol. 8, No. 5, 1976, pp. 21-35.

28. _V. Gersch, "Average Energy Exchanged in Oscillators," Journ. Acoust. Soc.

Amer., Vol. 46, No. 5 (Part 2), 1969, pp. 1180-1185.

29. J.L. Zeeman and J.L. Bogdanoff, "A Comment on Complex Structural

Response to Random Vibration," AIAA Journal., Vol. 7, No. 7, 1969, pp.

1225-1231.

30. _vV. Gersch, '_Mean Square Response in Structural Systems," Journ. Acoust.

Soc. Amer., Vol. 48, No. 1 ( Part 2), 1970, pp. 403-413.

31. J. _Voodhouse, "An Approach to the Theoretical Background of Statistical

Energy Analysis Applied to Structural Vibration," Journ. Acoust. Soe.

Method for

Sound Vib.,
- 35-

Amer., Vol. 69, No. 6, 1981, pp. 1695-1709.

32. R.H. Lyon, Statistical Energy Analysis of Dynamical Systems, MIT Press,

1975.