Probabilistic finite elements for nonlinear structural dynamics

COMPUTER METHODS NORTH-HOLLAND

IN APPLIED

MECHANICS

AND ENGINEERING

56 (1986) 61-81

PROBABILISTIC FINITE ELEMENTS FOR NONLINEAR STRUCTURAL DYNAMICS* Wing Kam LIU, Ted BELYTSCHKO Northwestern

and A. MAN1

University, Department of Mechanical and Nuclear Engineering, U.S. A.

Evanston,

IL 60201,

Received 14 March 1985 Revised manuscript received 24 July 1985

A finite element method applicable to truss structures for the determination of the probabilistic distribution of the dynamic response has been developed. Several solutions have been obtained for the mean and variance of displacements and stresses of a truss structure; nonlinearities due to material and geometrical effects have also been included. In addition, to test this method, Monte Carlo simulations have been used and a new method with implicit and/or explicit time integration and Hermite-Gauss quadrature has also been developed and used. All these methodologies have been implemented into a pilot computer code with two-dimensional bar elements.

1. Introduction The theory of stochastic processes is well established and its development is now available in standard texts. In general, the random uncertainties which are included in a stochastic process can be classified into three major categories: (1) physical uncertainty, (2) statistical uncertainty, and (3) uncertainty in the model. A detailed discussion of these topics can be found in, for example, [l-4]. Traditionally, uncertainty analysis in structural mechanics has concentrated on problems of an almost totally stochastic nature. Within this setting, even a single degree of freedom system with nonlinearities poses a formidable challenge and has not been solved satisfactorily. The most commonly employed solution technique is the Monte Carlo simulation (see e.g. [3]). In general, these simulation procedures are computationally repetitive and therefore expensive, even though they are easily applicable to both linear and nonlinear systems. For linear systems, nonstatistical methods such as ‘second moment analysis’, are available [2]. A related secondorder perturbation technique applied to a special class of linear structural vibrations is discussed in [5]. The emphasis is on the modal decoupling of the equations of motion with uncertain damping. The ‘second moment analysis’ has also been extended in [6] to define the mean and the variance of vector functions. This formulation is mathematically elegant, and Kronecker algebra and matrix calculus are employed. While this formulation has also been extended in [7] to linear stochastic systems with colored multiplicative noise, the direct application of this technique to nonlinear structural dynamics is not feasible. This is because in most nonlinear *The support of NASA Lewis Grant No. NAG3-535 for this research and the encouragement Chamis are gratefully acknowledged. 00457825/86/$3.50

0

1986, Elsevier Science Publishers

B.V. (North-Holland)

of Dr. Christos

W. K. Liu et al., Probabilistic elements for nonlinear structural dynamics

62

Physical Problem Wti

Finite Element -_B

Uncertainties

Solution of

Semi-discretization + and Aasembiy

Semi-discretized

-_+

Equations...

Fig. 1. Schematic of probabilistic

.. . repeated over Simulation Points MONTE CARLO SlMUlATiON (KS)

1 Response 4

Olld

Bounds A

methods.

structural analysis, concern lies more with deviations in loads from a deterministic path and with uncertainties in material properties to which a variance can be assigned, than with completely stochastic loads or systems. Hence, the goal of this paper is to provide techniques which are cost-effective and to enable the engineer to evaluate the effect of uncertainties in complex finite element models. It is also important to be able to treat the effects of uncertainties in a reasonably economical manner; standard Monte Carlo procedures are simply too expensive. Furthermore, the methods should be designed so that they can be incorporated into widely used finite element programs in a natural and concise manner. Thus, the approach should be integrable with the elemental discretization and nodal assembly procedures that characterize finite element theory and software implementation. In the next section, the formulation of the probabilistic finite element method (PFEM) is presented. In Section 3, the computational aspects of PFEM are discussed. In Section 4, the analysis of a two degree of freedom spring-mass probabilistic system is given. Results are also presented for a ten-bar probabilistic system with nonlinearities. The proposed PFEM method is compared to (1) Monte Carlo simulations (MCS) and (2) Hermite-Gauss quadrature (HGQ) schemes. All these methods are schematically depicted in Fig. 1, highlighting the major computational steps. In Section 5, the relative performance of the PFEM as compared to the other two methods is discussed. The reason for the limitation of each solution technique is also presented.

2. Formulation

of the probabilistic

finite element method (PFEM)

We consider the structural system to be governed by the following system of nonlinear algebraic equations which arises from a finite element discretization:

W.K. Liu et al., Probabilistic elements for nonlinear structural dynamics

il!fa-+f(b,

d, d) = F(t)

63

(2.1)

)

where M,f, d, and F are the generalized mass, internal force, displacement, and external force, respectively; and a superscript dot represents the material time (t) derivative. While the internal nodal forces are obtained from one variational statement, they are segregated for convenience_ The probabilistic effects are described through the q-dimensional random vector b; this can include the probabilistic distributions of the material properties; the mass M is assumed to be deterministic. Ail these probabilistic distributions, as reflected in the variance of the material properties, the composite load spectra, etc., are represented by the generalized variance vector, Vat(b). We shall denote the expected value operator by E[*j and use second-order expansions, so E[*] is given by

where q is a vector function of the random variables. The superposed bar denotes ‘at the mean value of b’ and the symbol Cov represents the covariance; summations on i and j from I to q are assumed. If b, is uncorrelated to b, for i # j, then

Cov(b,, b,) = 0 ,

fori#j,

Cov(b,, b,) =Var(bi)

CW

no sum on i .

,

Applying the expected value operator

to (2.1) yields

E[Md] + E[f(b, d, d)] = E[F(t)] . Employing

(2.4)

(2.2) and the chain rules:

E[Md] =imi + $M

& ’

Co%, b,) , I _

b,) +

Co@,,

($ g +$ $} Cov(b,, bj) 1

a2v ab,db,

+K

a”2

-db,db,

-

Cov(b,, b,)

f

(2.5)

,

where ci and b have been replaced by Yand a, respectively. C and K are the damping and stiffness matrices, respectively. They are C = af/dv

and

K = af/c?d .

(2.6)

W. K. Liu et al., Probabilistic elements for nonlinear structural dynamics

64

In the case of a linear structure, f is given by f=Cv+Kd. For simplicity, let us assume that (2.3) holds. This assumption is quite suitable for finite element models which are built up from discrete structural elements, such as bars and beams. Using this simplification and applying perturbation techniques on (2.9, equation (2.4) can be shown to yield F(t)

Ma +f=

(2.7a)

)

MAi+CAi+iAd=AF,

(2.7b)

where (2.8a)

q

a*a

Aa = $I?1 z

Wb,) ,

(2.8b)

Wb, >,

(2.k)

Var(b,) .

(2.8d)

I

Av = i

q d*v

I?1z

I

q a’d

Aa= i I:12

J

Once Ali, AV, and Ad are obtained by solving (2.7b), the second-order

means are

E[li] = E[a] = i + Ati, E[d] = E[v] 2: i + Ai,

(2.9)

E[d] = d + Ad. If one is interested in the deviations in response from a deterministic path due to the uncertainties in material properties to which a variance can be assigned, the number of time integrations (simulations) reduces to only two. These two simulations are Mti+f=F,

(2.10a)

MAa+tAv+lfAd=AF,

(2. lob)

where q da -Ab,,

Aa=x ]=I

ab,

q dv

Av=]T,ab

Ab,, I

W. K. Liu et al., Probabilistic elements for nonlinear structural dynamics

Ad=z

65

4 ati --Ab,, !=I db,

and A bj is equal to the ~~e~~~g~e~ valence of b,. It is also necessa~ to obtain the se~s~ti~ty vectors, i.e. a(.)la see Section 3. Finally, once the means and the sensitivity vectors are determined, the variance vectors can be computed easily by the following first-order formulas:

(2.11)

In the case of uncorrelated terms are denoted by Var(a) - $

b, the covariance matrix becomes a diagonal matrix and the diagonal

($-)’

Var(bj) ,

J

v=(v)

var(d)

=

$

($)2

= iI

($)’

J

I

Var(b,) ,

(2.12)

Var(b,) .

Similar procedures can also be developed for the probabilistic distributions of stresses. However, this direct approach can be very expensive, if the number of random variables is greater than the number of requested probabilistic distributions of stresses. For this situation, an alternative approach, termed an ~~jo~~~~robab~l~s~ic stress analysis, is developed. This is described in [9]. ~~~A~~

2.1.

The uncertainties discussed here are described by discrete random variables. Physical parameters, such as material properties, are often continuous functions in space. When there are uncertainties associated with these parameters, we have random fieZds [3]. The probabilistic distributions at any two points can be represented by a ‘correlation function’. One way to adapt the above procedures to ‘random fields’ is to first do a finite element discretization of the correlation function and thus obtain the covariance matrix. Once this matrix is obtained, the PFEM method as developed here could be used with minor modifications [12,13].

W. K. Liu el: al., P~~b~bi~~sc~c elements

56

for

nonlinear structural dynamics

REMARK 2.2 To the authors’ knowledge, (2.7)-(2.12) represent the first consistently derived second moment probabilistic finite element method (PFEM) which can readily be adapted to existing deterministic finite element computer programs. The second-order terms A&, Ai, and Ad are computed directly from the second moment mean (2.7b). Consequently, (2.9) are second-order accurate and (2.11) are first-order accurate. REMARK 2.3, The complete probability distributions are not available for most random variables, except perhaps for the first two moments. Methods such as MCS or HGQ usually require knowledge of probability density functions. The PFEM method requires only the first two moments and is therefore widely applicable.

3. Computational

aspects of PFEM

The computing procedures essentially involve time integrations of the various equations derived in the previous section. In general, the sensitivity vectors can be obtained directly by integrating the sensitivity equations in time. However, this is not possible for some nonlinear systems. In such cases, the usual procedure is to calculate the derivatives by finite differences [8]. Calculating the finite difference derivatives increases the computation for a probabilistic system. However, results obtained are excellent when compared to the solutions obtained by other methods. The computing procedures for linear and nonlinear systems are described separately below. 3.1. L&rear systems For a linear system, equations

(2.7) become (3.la) (3Sb)

The solutions of (2.7a) and (2.7b) are obtained in sequence so that the additional computation due to the latter is minimized. The solution algorithms, such as implicit and/or explicit time integration, used in (3,la) can be applied directly to (3.lb) with the formulation of only one additional vector ~~~&tion A#. If we examine (2.8a) closely, it can be shown that Al;*;can be computed element-wise once G, d, dF/db,, and Jaldb, are given. In addition, the corresponding variation of the efementaE nodal forces can then be assembled into a description of the probabilistic distribution of the elemental nodal forces for the complete finite element model. It can be easily shown that the governing equations for the sensitivity vectors are obtained by differentiating (2.1) with respect to bj. They are (3.2) where

67

W. K. Liu et al., Probabilistic elements for nonlinear structural dynamics CALCULATED QUANTITIES

ITEM

EQUATION NUMBERS

1

The mean values __ i.e. d, __ v and a _

2

The sensltlvlty derivatives for each random variable b,, J = 1. ,,., q -_ ad I.e., i-d_ ab J

3

BV

dd, _

ab J

.4;

_

(3.la)

1

(3.2a)

q

(3.lb)

1

8a L ab J

and

The second order varlatlons I.e.,

NUMBERS OF INTEGRATIONS

and

Aa

Total = q + 2

Fig. 2. Computational

a13 q=-

a@

af ab,

steps in PFEM.

d, v=constant

Or

3

I

=

From (3.1) and (3.2), it can be seen that the whole procedure uses the same effective stiffness matrix, so only one matrix needs to be triangulated. To evaluate the mean and variance from (2.9) and (2.12), the total number of time integrations required is q + 2. These are: one integration to evaluate the displacement, velocity, and acceleration at the mean value of b (3.la); q integrations to evaluate the sensitivity vectors (3.2); and one more integration to evaluate the second-order variations (3.lb). The computational steps involved in PFEM are shown in Fig. 2. Notice that all time integrations employ the same effective stiffness matrix; parallel computation procedures could be employed, thereby increasing the efficiency tremendously. 3.2. Systems with material and geometrical nonlinearities As in the linear case, the displacement, velocity, and acceleration at the mean value of b are obtained by integrating (2.7a). The relative merits between implicit and explicit time integrations are considered here for a probabilistic nonlinear system. By total differentiation of (2.10a) with respect to b,, i.e. dldb,, we have:

Mda+df=() dbj

db,

’

(3.3a)

W. K. Liu et al., Probabilistic elements for nonlinear structural dynamics

68

(3.3b) Equations

(2.7a) and (3.3a) can be written as (3.4a)

M-

%+1 db,

-

a+1

4@~ d,+A

db,

ab,

_=-

+K

(3.4b)

’

where f, +1 and i are the internal force vector and the ‘tangent stiffness matrix’, respectively, evaluated at b, d,,,, , and t,+ 1. Equations (3.4) can be solved by the implicit Newmark-P algorithm [lo]. The ‘mean value’ equation (3.4a) can be solved by Newton-Raphson iteration

where the residual vector is given by

Ma,+J”

--Y

r n+1=

[F,+1

-i+1-

and the, effective stiffness matrix is ii*

=

{M + p

At2

Ii}“.

The symbol v represents the equilibrium iteration counter at time step IZ+ 1, and iterations are repeated until Aaron approaches zero. Similarly, the first-order sensitivity equation (3.4b) can be written as

aa

2

x*$=-pAt I

4fn+l

db+MI

a&+,

db,

’

where

dri,+l

ab,

-

-

=$+Atg+(i-P)At2$& I I

I

It is observed here that the effective stiffness matrixK* is iden_ticalin both (3.5) and (3.6). Since the triangulated K is given during the iteration procedures, dd, + 1la b, can be obtained simply by forward reductions and back-substitutions; therefore, the number of time integrations is still q +2. The main advantage of employing implicit time integration is its unconditional stability. Therefore, the above methods are best suited for structural dynamics problems dominated by low-frequency response. For impulsive and short-duration transient problems, (2.7a), (3.3a), and (3.3b) can alternatively be solved by explicit integrations. Since fib, d) is nonlinear, the sensitivity vectors can be obtained by central differences. Equations (3.3) are approximated by

W. K. Liu et al., Probabilistic elements for nonlinear structural dynamics

69

(3.7a) -2ti+a-

=o

Ab;

7

(3.7b)

where a+ = a(b + Abj) ,

a- =

f + = f(b + Ab,, d-j,

f-

a(b - Abj) ,

=f(b-

Ab,, d-) .

d+ and d- are similarly defined and Ab, is defined by AbI = (O,O,. . . , Ab,, 0,. . . ,O)' , where t denotes the transpose. With this integrations would still be 4 + 2. However, 2q + 1. These are: one integration for differencing for (3.7a) and (3.7b). Apart mixed time implicit-explicit algorithms [ll] of the algorithms can be achieved.

computational procedure, the total number of time the number of internal force calculations would be the mean (2.7a) and 2q integrations with finite from purely implicit or purely explicit algorithms, could also be employed so that the attributes of each

4. Numerical examples EXAMPLE 4.1. A two degree of freedom spring-mass system. The performance of PFEM, the new method developed here, is evaluated via a two degree of freedom spring mass system. The mean is second-order accurate and the variance is first-order accurate in this example. The computed results are compared with those obtained employing (1) Monte Carlo simulation (MCS) and (2) H ermite-Gauss quadrature (HGQ) schemes. The two latter methods as implemented here are reviewed in Appendix A. The problem statement is depicted in Fig. 3. A sinusoidal vector forcing function is used:

F(t)=250 x lo6 Sin(2000t) m,=0.372

k,=24.0 x lo6

m,=0.248

k,=12.0 x lo6

Fig. 3. A two degree of freedom example.

W. K. Liu et al., Probabilistic elements for nonlinear structural dynamics

70

F(t)

= [

O*O

25.0X lo6 sin 2000t

] .

(4.1)

The random spring constants K, and K, are normally distributed with a coefficient of variation (i.e. cr/p) equal to 0.05. The mean spring constants are 24 x lo6 and 12 x lo”, respectively. DISPLACEMENT(NODE 1)

---___--.MC-

3 G

3

I

-1.000~~

-2.ml .ow

.125

950

375

iGE

,875

(s&,&7~;

Fig. 4. Comparison of the mean displacement at node 1 using: (1) probabilistic Hermite-Gauss quadrature (HGQ); (3) Monte Carlo simulation (MCS).

1.000

1.

1.125

50

finite element method (PFEM);

(2)

DISPLACEMENT(NODE 1) A00

_._.-.PFE,J

---HGQ

_-_-____ MCS

f ;

,200

d

s

j_

.lOO

.ow

.oQo

.125

.250

375

875

1.000

1.125

1

50

Fig. 5. Comparison of the variance of displacement at node 1 using: (1) probabilistic finite element method (PFEM); (2) Hermite-Gauss quadrature (HGQ); (3) Monte Carlo simulation (MCS).

71

W. K. Liu et al., Probabilistic elements for nonlinear structural dynamics The

deterministic masses m, and m2 are 0.372 and 0.248, respectively. A stiffness-proportional damping of 3% is included. The probabilistic equations derived earlier are solved by the implicit Newmark-/ method [7]. The mean amplitude d, is depicted in Fig. 4 for all the three numerical methods-PFE~, HGQ, and MCS. The PFEM solution compares very well with the other two methods. For the variance of d, the PFEM solution, plotted in Fig. 5, seems to overshoot the DISPMCEMENT (NODE 2) s.wo

-.-.-.-

pfqT&J

HGQ

-

______--

MCS

Fig. 6. Comparison of the mean displacement at node 2 using: (1) probabilistic Hermite-Gauss quadrature (HGQ); (3) Monte Carlo simulation (MCS).

mo

finite element method (PFEM);

(2)

DIsPt.mh4EN~(NODE 2) _._.-.f3FEM

---HGQ

________ Mf'S

8

4: 2

.400

a

s

j_

.zSS

.wo .m

.12S

.2So

875

r.mo

Fig. 7. Comparison of the variance of displacement at node 2 using: (1) probabilistic (PFEM); (2) Hermite-Gauss quadrature (HGQ); (3) Monte Carlo simulation (MB).

50

finite element

method

72

W.K. Liu et al., Probabilistic elements for nonlinear structural dynamics

variance at large times. The mean and variance of d, are similarly compared and depicted in Figs. 6 and 7. The maximum coefficient of variation of the displacements d, and d, are found to be 0.13 and 0.10, respectively. The +3a bounds for the displacements d, and d, are plotted in Figs. 8 and 9, respectively. DISPLACEMENTBOUNDS AT NODE l(PFEM) 3.ooo

-s.ow .wo

.........._. "PPER

-MEAN

.125

.250

________. LOWER

,675

575 itiE

Fig. 8. +3a bounds of the displacement

l.MK)

1.125

I. %I

&(x&~

at node 1 using probabilistic

finite element method (PFEM).

DISPLACEMENTBOUNDS AT NODE 2(PFEM) . . . . . . . . . . . . . . UPPER

-6.000 .wa

,125

.260

,376

Fig. 9. k3v bounds of the displacement

_---___

_. LOWER

675

at node 2 using probabilistic

l.wO

1.125

1.

finite element method (PFEM).

W. K. Liu et al., Probabilistic elements for nonlinear structural dynamics

73

EXAMPLE 4.2. A ten-bar probabilistic system with material and geometrical nonlinearities. The problem statement is depicted in Fig. 10. The load time function, which is also shown in Fig. 10, is applied at node 3. This particular loadtime history is chosen such that only four of the ten bars, elements 1, 3,7, and 8, will yield. Therefore, the probabilistic model can be simplified by choosing the yield stresses of these four elements as the normal random variables which have the major impact on the response. The coefficient of variation is 0.05. Since the other six elements do not come close to yield, they are considered deterministic variables. With this approach, instead of 59049 analyses, only 81 analyses are required for the Hermite-Gauss quadrature method. The justification for this drastic simplification is explained in detail in Appendix A. For the PFEM method, the finite difference derivatives are evaluated with an interval Ab, equal to O.O5b, and the equations are solved by explicit time integration. The mean is second-order accurate whereas the variance is first-order accurate. The Monte Carlo simulation results are obtained with 400 simulations. The probabilistic displacement and stress solutions at selected locations are given in Figs. 11-14. The maximum coefficient of variation of the displacement of node 1 is found to be 0.13 and that of the stress in element 1 is 0.11. For this example, the three methods (PFEM, HGQ, and MCS) have been employed and they all compare quite well. The bounds of the displacement and stress can be estimated based on the Chebyschev inequality P((x - ~1 Z- na) G l/n* ,

I-

360”

--+c-

n > 0,

360”

(4.2)

__(

E =30.0 x lo6 E,=30.0 x lo4 A =6.0 p =0.30 q=15000.0 P =175.0 x lo3 P*=O.O Fig. 10. Problem

statement

of a ten-bar

nonlinear

structure.

W. K. Liu et al., Probabilistic elements for nonlinear structural dynamics

74

DlSPl.ACEMENT AT NODE 1 -._.-.f.,CQ __---___. MCS

50.001 1 --PFEM

37 SW-

5 I

lii_l

25.om--

12.500~-

.ow

.I60

..I20

300

.610

.BW

.om

1.120

t ,200

1.440

1 00

TIME (-SEC)

Fig. 11. Comparison of the mean y-displacement of node 1 using: (1) probabilistic (2) Hermite-Gauss quadrature (HGQ); (3) Monte Carlo simulation (MCS).

o.wa \

DISPLACEMENT AT NODE 1 _._._.HGQ _______~ MCS

-PFEM

4 500--

i

5

finite element method (PFEM);

_

x.wo--

s j_

1.500 . .

.clon~

.wo

,100

.JZO

,460

340

TN&C)

.ow

1.120

1.280

1.440

Fig. 12. Comparison of the variance of the y-displacement of node 1 using: (1) probabilistic (PFEM); (2) Hermite-Gauss quadrature (HGQ); (3) Monte Carlo simulation (MCS).

I. 00

finite element method

W. K. Liu et al,, Probabilistic elements for nonlinear structural dynamics

75

STRESS IN ELEMENT1

25.wo-

PFEM

-

.WO, .OOO

,320

_._._.HGQ

,840

---_____. MCS

1.120

Fig. 13. Comparison of the mean stress of member 1 using: (1) probabilistic Hermite-Gauss quadrature (HGQ); (3) Monte Carlo simulation (MCS).

30.000 -

PFEM

I .440

I.600

finite element method (PFEM);

(2)

STRESS IN ELEMENT 1 _._._._ HGQ -_______ MCS

Fig. 14. Comparison of the variance of the mean stresses of member 1 using: (1) probabilistic (PFEM); (2) Hermite-Gauss quadrature (HGQ); (3) Monte Carlo simulation (MCS).

finite element method

76

W.K. Liu et al., Probabilistic elements for nonlinear structural dynamics

where p = E[x] and a* =Var(x). The +-3a bounds (i.e., n = 3) for the displacement and stress are plotted in Figs. 15 and 16, and the solutions can be expected to be within these bounds with 89% confidence level.

DISPLACEMENTBOUNDS AT NODE l(PFEM) 60.001

_._._._ UPPER

-MEAN

_____-__ LOWER _/-.-‘-._.___._.--._.__/._../

.’ 45.oiJo -.

.’

.’ .’ .’

9 G

x 3

5

3o.wo--

_____------__________ -----______________

a

2 15.000-~

.ooai

.oQo

,160

,320

,460

Fig. 15. k3a bounds of the y-displacement

.wO

460

1.120

at node 1 using probabilistic

1.280

1.440

1. 00

finite element method (PFEM).

STRESS BOUNDS IN ELEMENT l(PFEM) 3o.ooo

-MEAN

_._.-.UPPER

________ LOWER

22.500-~

A-

0 25 ”

2

v

l.%ao--

x

z P -f

?SOO--

Fig. 16. f 3u bounds of the stress in member

1 using probabilistic

finite element method (PFEM).

W. K. Liu et al., Probabilistic elements for nonlinear structural dynamics

5. Comparisons

77

among the three methods and conclusions

Based on these numerical studies, we have drawn the following tentative conclusions: (1) Although all three methods agree very well and are evidently comparable in accuracy, PFEM is the most efficient solution procedure for small to medium-size problems. The relative computational efficiency of the three methods is summarized in Table 1. The number of time integrations required for a general structure with q random variables can be summarized as follows: (i) PFEM with partial derivatives evaluated directly: q + 2. (ii) PFEM with partial derivatives evaluated by finite difference: 2q + 1. (iii) HGQ with three-point quadrature: 3q. (iv) MCS with simple Monte Carlo simulation of sample size N: N. (2) Although PFEM is expected to be most accurate when the variances are small, it performs quite well even when the response shows a large coefficient of variation (e.g., 0.13 for the displacement at node 1 in the ten-bar structure). This could be attributed partly to the nature of the probabilistic distribution. For most distributions, values of response faraway from the mean are less likely to be found than those near the mean. Hence second moment analysis about the mean turns out to be quite accurate. (3) The three methods are applicable to linear and nonlinear systems. In linear systems the partial derivatives can be obtained directly. In nonlinear system the brute force method is to obtain these derivatives by finite differences. We are currently investigating ways to compute these derivatives efficiently. However, the methods are problem-dependent. (4) A minor drawback of PFEM is that its accuracy deteriorates for large times even with structural damping. An explanation is given in Appendix B. We are currently investigating several ways of improving this. (5) PFEM can be easily incorporated into widely used finite element programs. (6) A PFEM analysis can be obtained with q + 2 simulations if Cov(b,, bj) # 0 for i # j. For this purpose the bi must be transformed into another set of random variables c, through an eigenproblem such that Cov(c,, c,) = 0 for i # j, and in most cases only a few modes are sufficient [ 131. (7) Currently the PFEM is being extended to the transient analysis of nonlinear continua. The details of the method can be found in [13]. Since this method involves only matrix and vector assembly it can be incorporated in a natural and concise manner in general-purpose finite element programs.

Table 1 Relative computer

2-bar structure IO-bar structure

time by the 3 methods PFEM

HGQ

MCS

1

8

400

1

4

60

78

W. K. Liu et al., Probabilistic elements for nonlinear structural dynamics

Appendix A. Review of Monte Carlo simulation methods and Hermite-Gauss quadrature schemes A. 1. Monte Carlo simulation methods

Various Monte Carlo simulation techniques are now available. In our analysis, the ‘simple’ Monte Carlo method is used. The important feature about the Monte Carlo method is its flexibility. In other words, the computational procedures are the same irrespective of whether the mode is linear or nonlinear, so long as the solution can be obtained from the governing equation. The main idea behind Monte Carlo simulation methods is to randomly generate values of the random variables subject to the probability density function and to calculate the output corresponding to these values. From this set of output, the probabilistic distribution properties, such as the mean and variance, are statistically estimated. In the analysis of the two degrees of freedom probabilistic system, a normal random generator is used. This normal random number generator, RANF, is available on the Northwestern University CDC system. It has been well tested and for large sample size, the distribution is close to normal. This ‘closeness’ can also be estimated by the so-called Central Limit Theorem, which states as follows: “if a population has a finite variance a2 and mean value CL,then the distribution of the sample mean approaches the normal distribution with variance a2/n and mean Al.as the sample size II increases.” The sample size used in this analysis is 400. Both the spring constants K, and K2 are randomly generated and the corresponding displacement solutions are calculated using the exact deterministic solutions. A.2.

Implicit time integration with Hermite-Gauss

quadrature scheme

Let us consider a linear two degrees of freedom probabilistic

Ma,,, + K4+1 = Fn+l

7

system

(A4

where

IZis the time step number and the initial conditions d,, vO, and a, are given. At each time step 12, (A.l) is solved by the Newmark-/ algorithm with /? = 0.25, y = 0.5, and At = O.O2T,,,,,, where T,,, is the smallest fundamental period. The finite difference matrix equation can be shown to be Keffdn+l = Feff ,

(A-3)

Keff = M + p At2 K ,

(A4

F eff= pAt’F,+,

(A.3

where

Jn+l= d,+Atv,l+(i Cl+1= V, + At(l-

+ M&+1 , -P)At2a,, y)a, .

(A.71

W.K. Liu et al., Probabilistic elements for nonlinear structural dynamics

Once d,+l is determined

from (A-3), i.e.,

d n-i-l = (p)-‘F”ff a n+l

79

,

w9

andVn+l can be determined as follows: a n+l

=

(4+1- d,+,)ipAt27

V n+l

=

cl+1+ Y Atan+,

(A-9) (A. 10)

’

The solution procedures are then repeated with tz replaced by n + 1, until y1At is greater or equal to a desired time. Since K, and K2 are random variables, d,+l by (A.8) is ‘implicitly’ a function of KY,and K,. Using the basic definitions of mean and variance, the expected value of d n+l is (A.ll) where fK,(K1) and fKz(K2) are the probability density functions for K, and K2, respectively. In writing (A-11), the assumption that K, is uncorrelated to K2 has been employed. Once the expected value is evaluated, the variance of d,,, can then be computed according to

Var(d,+,)= W~+ll - (Wn+11)2. The Hermite-Gauss in (A.ll) by

quadrature

scheme is to approximate

(A. 12) the double integrals which appears

where II, and n2 are the number of integration points for K, and Kz, respectively, and w, and wI are their corresponding weights. As might be figured from the above equation, if the number of random variables is m, the number of simulations N is N=n,xn,x--xn,

(A. 14)

and N grows exponentially. Therefore, unless the number of random variables is small, this method is not recommended. However, if a physical situation dictates that some of the random variables can be excluded in the calculations, N can then be reduced significantly. Under this circumstance, the Hermite-Gauss quadrature scheme can be an efficient and accurate method. An example of this practical situation has been demonstrated in Section 4. In the analysis of the ten-bar structure, it was predetermined that only four of the ten bars will yield and the yield stresses are chosen as normal random variables. If 3 points are used in evaluating each normally distributed K, in (A.13), the weight and quadrature points are

W.K. Liu et al., Probabilistic elements for nonlinear structural dynamics

80

w, =(i,

2, i)

and

(A. 15)

K,=(/.-3a,jJ,CL+3&

respectively. And for this ‘bias’ integration each time step n) is

procedure,

the number of simulations required (for

Whereas, if one cannot observe a priori that six of the ten bars will not yield, the number of simulations required becomes N = 3l” = 59049 , which makes this a handicapped

method.

Appendix B. Resonant excitation of response sensitivities The equation of motion and the sensitivity spring-mass-damper system are

equation

for a single degree

of freedom

uw

Mi + Ci + Kx = F(t) ,

where

M, C, K are assumed to be dependent on the parameter 6; we are interested in the sensitivity of the response x(t) to this parameter b. Let F(t) be such that x(t) is stable. Under this condition, it is shown below that the response sensitivity axlab is resonantly excited. The damped natural frequencies of the system (B.l) and the sensitivity equation (B.2) are the same. The excitation p(t) involves x, 1, and X in (B.2) and is, therefore, a resonant excitation. Thus approximations for x(t), such as

x(b, + Ab, 0 =x@,,, t>+

b=b

Ab2

(B-4)

0

at b, for any small interval Ab, are valid only for a short duration and the accuracy deteriorates rapidly thereafter. Since the PFEM equations (2.7) use the first- and second-order response sensitivities, they are valid for a short duration only. A similar phenomenon is also observed in the transient response of nonlinear structures. A possible explanation for this phenomenon is that the time ‘t’ has a multiplying effect on the interval ‘Ab’ in the second and third terms in (B.4), and this

W. K. Liu et al., Probabilistic elements for nonlinear structural dynamics

81

results in the deteriorating accuracy. However, the PFEM equations developed in Section 2 are suitable for application when one is interested in a short-time history, e.g., the response due to an impulsive load. Although the PFEM method, as presented here, is not applicable for long duration transient analysis of structures, efforts are being made to remedy this [13].

References [l] Y.K. Lin, Probabilistic Theory of Structural Dynamics (McGraw-Hill, New York, 1967). [2] A.H.S. Ang and W.H. Tang, Probability Concepts in Engineering Planning and Design, Volume I, Basic Principles (Wiley, New York, 1975). [3] E. Vanmarcke, Random Fields, Analysis and Synthesis (MIT Press, Cambridge, MA, 1984). [4] T.K. Caughey, Nonlinear theory of random vibrations, in: C.S. Yih, ed., Advances in Applied Mechanics, 11 (Academic Press, New York, 1971) 209-253. [5] S. Nakagiri, T. Hisada and T. Toshimitsu, Stochastic time-history analysis of structural vibration with uncertain damping, in: Probabilistic Structural Analysis, PVP-Vol. 93 (ASME, New York, 1984) 109-120. [6] F. Ma, Extension of second moment analysis to vector-valued and matrix-valued functions, Intemat. J. Nonlinear Mech., to appear. [7] F. Ma, Approximate analysis of linear stochastic systems with colored multiplicative noise, Internat. J. Engrg. Sci., to appear. [8] R.T. Haftka, Techniques for thermal sensitivity analysis, Internat. J. Numer. Meths. Engrg. 17 (1981) 71-80. [9] W.K. Liu and T. Belytschko, Variational approach to probabilistic finite elements, Progress Report to NASA, December 1984. [lo] N.M. Newmark, A method of computation for structural dynamics, ASCE J. Engrg. 85 (1959) 67-94. [ll] W.K. Liu, T. Belytschko and Y.F. Zhang, Implementation and accuracy of mixed-time implicit-explicit methods for structural dynamics, Comput. & Structures 19(4) (1984) 521-530. [12] W.K. Liu, T. Belytschko and A. Mani, Probabilistic finite elements for transient analysis in nonlinear continua, in: O.H. Burnside and C.H. Parr, eds., Advances in Aerospace Structural Analysis, AD-09 (ASME, New York, 1985) 9-24. [I31 W.K. Liu. T. Belytschko and A. Mani, Random field finite elements, Internat. J. Numer. Meths. Engrg., to appear.