Stochastic finite element method for spatial distribution of material properties and external loading

0045-7949@4Wo5o4-4

Compurerr & Slrucrures Vol. 55. No. I, pp. 41-45, 1995 Copyright 8 1995 Ekvier Science Ltd Printed in Great Britain. All rights reserved 0045-7949/9539.50 + 0.00

STOCHASTIC FINITE ELEMENT METHOD FOR SPATIAL DISTRIBUTION OF MATERIAL PROPERTIES AND EXTERNAL LOADING S. Chakraborty and S. S. Dey Department

of Civil Engineering,

Indian

Institute

(Received

of Technology,

7 December

Kharagpur,

721302,

India

1993)

Abstract-In this paper, a stochastic finite element method is presented for the analysis of structures having statistical uncertainities in both material properties and externally applied loads, which are modelled as a homogeneous Gaussian stochastic process. The Neumann expansion technique has been used for the inversion of a stochastic stiffness matrix. The digital simulation technique was adopted to generate the deviatoric part of the stiffness matrix and load vector. A beam problem was taken up for comparison of results obtained by Neumann expansion and the direct Monte Carlo simulation technique. The comparison of the results shows that the values approach those obtained by direct Monte Carlo solutions as the order of expansion of the Neumann series is increased. An excellent agreement was achieved for cases when the coefficient of variation is comparatively small.

NOTATION

are usually spatially distributed and correlated over the region of structures and should be modelled as random variables. Though the modern deterministic finite element method (FEM) is quite elaborate and includes sophisticated mathematical modelling and analysis, a highly refined method of computer evaluation cannot be applied to structures with such uncertainties. To deal with this problem, a stochastic finite element method (SFEM) has been developed using the statistical nature of loads and material properties. The current SFEM is essentially a combination of deterministic FEM with probablistic methods. The most widely used technique is the perturbation method [1,2]. The method was systematically developed by Nakagiri and Hisada [3,4]. Shinozuka et al. investigated the probabilistic model for the spatial distribution of material properties using simulation methods [5-71. The Neumann expansion technique was applied by Shinozuka and Yamazaki. Vanmareke and Grigoriu [8] studied the problem of random shear beams, making use of the scale of fluctuation of a random process. Spanos and Ghanem [9, lo] investigated using orthogonal decomposition to expand the random media. In this paper both the stiffness matrix of the structure and the load vector have been taken as random quantities. The stochastic fields are simulated by Cholesky decomposition of the respective covariante matrix of the random fields. The Neumann expansion technique has been used in the finite element solution of the response variability. The effect of shear deformation has also been considered in the derivation of a stiffness matrix for the beam

width, depth and length of the beam mean value of the modulus of elasticity and load standard deviation of the modulus of elasticity and load coefficient of variation of the modulus of elasticity and load reciprocal of correlation length of random fields, i.e. modulus of elasticity and external loading autocorrelation function of the random fields, i.e. modulus of elasticity and external loading separation vector element stiffness matrix global stiffness matrix global mean stiffness matrix deviatoric part of the stiffness matrix global displacement vector global load vector deviatoric part of the load vector modulus of rigidity Poisson’s ratio shear correction factor identity matrix indicates matrix inversion indicates matrix transposition indicates variance indicates covariance indicates mathematical expectation

INTRODUCTION

Most

structural

engineering

problems

involve

ran-

in loading conditions or in material properties. There is always some probability of variation from the nominal values. So uncertainty abounds in structural engineering problems. These uncertainties domness

41

S Chakraborty and S. S. Dey

42

element. The results are compared with those tained by direct Monte Carlo methods.

ob-

The vector generated inal covariance matrix,

by eqn (5) satisfies the orig-

E(cc UT)= E[LZ(LZ)‘] SYSTEM

The value of Young’s modulus and load per unit length in the structural member is considered to vary stochastically along its undeformed axis. The spatial variations are assumed to be one-dimensional homogeneous stochastic processes. These can be described as E(x) = B[l + e(x)]

(l)

F(x) = m +f(x)l,

(2)

where ,!? and F are the respective mean values of Young’s modulus and load intensity, e(x) and f(x) are modelled as zero mean homogeneous isotropic Gaussian fields, characterized by their respective variantes, CJ~ and a:, and autocorrelation functions, R,(5)

and &(O,

where (3)

DIGITAL

GENERATION

= L . E(ZZT)

LT = c,,

DESCRIPTION

OF THE RANDOM

As E(ZZT) = [I], however, a fairly large sample size is required for the simulated covariance matrix to approach the target covariance matrix. FORMULATION OF THE STOCHASTIC ELEMENT METHOD

Following the general FEM, [ll, 121 we have equation:

procedure of deterministic the following equilibrium

(7) where [K], {U} and {F} are the global stiffness matrix, unknown nodal displacement vector and equivalent nodal load vector, respectively. In the derivation of the element stiffness matrix, the effect of shear deformation has also been taken into account. For the present beam problem the element stiffness matrix consists of

where

FIELDS

In order to perform the SFEM analysis using simulation, the random fields must be generated. For this purpose, the structure is divided into an appropriate number of small finite elements so that the property values within each element can be considered approximately constant. So, if there are n elements in the total structures, then there are n property values associated with each zero mean homogeneous random field. If LX(X)is a zero mean homogeneous Gaussian process, then the values cli = u(x,) (i = 1,2,3 . . n) are random with zero means, but correlated. Thus the correlated random vector {E} = [GI,,CI~,a3,. . , ct,IT can be obtained as

(9) O-10

1 L

L 2 L2

L

-1

z L

L2

24

(10) L

where {Z) = a vector consisting of n independent Gaussian random variables with zero mean and unit standard deviation, and [L] = a lower triangular matrix obtained by the Cholesky decomposition of the covariance matrix C,, defined as

1

24

i

L2

24

E = 2G(l + v).

--

L 24

L2

1 (11)

Here one-point quadrature integration has been performed for deriving the shear stiffness part to avoid the stiff element [ 121.

Cov(cc,, r2). ..... C,, = E(cc, aT) =

FINITE

Var(a:) ......

Cov(cr,, cc,) Cov(a,, a,)

(6)

Spatial distribution of material properties and external loading

43

Now splitting the stochastic stiffness matrix and load vector into their corresponding mean and random deviatoric parts we get

The correlation functions for both the modulus of elasticity and load are assumed in the form of an exponential function, as in eqn (3) and (4). The correlation function includes two parameters, [K] = [iI -I- WI (12) the standard deviation and correlation distance. The standard deviation of the deviatoric part represents {F} = {i) + (@I. the coefficient of variation of corresponding random (13) parameters. The correlation length determines the From eqn (7) displacement can be written as rate of decay of the correlation function with distance. This value, both for modulus of elasticity and {V}=[K+AK]-‘{F+AF). external load, i.e. C, and C,, is taken as 1.0 through(14) out the numerical example. Now, expanding [K]-’ by the Neumann expansion The direct Monte Carlo simulation of the beam response is undertaken to compare with the present we get results obtained by the Neumann expansion method. In direct Monte Carlo simulation the random Gaus[K]-’ = [i + AK]-’ = i (-Q)‘(R]-‘, (15) sian vectors for loads and modulus of elasticity are r=O generated through the Cholesky decomposition of the where Q = [&‘[AK]. covariance matrix [using eqns (5) and (611. Upon generating the rigidity profile of the beam, its reSo using eqn (15) in eqn (14), sponse for the generated random load is determined relying on deterministic FEM. The procedure is (U) = i (-Q)‘[Z$‘(f; + AF) repeated several times to produce an ensemble of the ,=O beam response along the span. Now statistical algorithms are used to compute the mean value and + i (-Q)'[&'{AF}. = f. (-QH~I-‘16 standard deviation of the response at different nodes. ,=O The sample size for simulation purpose should be (161 large enough. The statistical fluctuation of standard deviation obtained by simulation is plotted against The expansion series in eqn (16) is terminated sample sizes. It is observed from Figs 2 and 3 that after a few terms, depending on accuracy and converfluctuation is within a tolerable range after 2000 gency rate. However, for the convergence of the samples. The figures also show the worst case Neumann expansion series, all the eigenvalues of (or = 0.25 and er = 0.29, and the statistical stability [Q] = [K]-‘[AK] should be less than unity. As we will be much better for lower values of coefficient of have used the truncated Gaussian distribution for the variation of the random parameters. So a sample size numerical analysis, all the eigenvalues are always less of 2000 is used throughout the numerical example. than one. However caution should be taken in generating the random Gaussian vectors for large values of coNUMERICAL EXAMPLE efficient of variation of the stochastic fields, because The problem involves a cantilever beam with there is a possibility of generating negative values of stochastic modulus of elasticity along its length, the parameters. During the simulation, the small subjected to a distributed stochastic static load. This number of random moduli of elasticity or loads with example has been taken to facilitate comparison of results with those obtained by direct Monte Carlo 3,Obh /lOOL 7 simulation. STANDARD DEVIATION OF MODULAS =o,2S The physical dimensions, material properties and OF ELASTICITY IO, 1 loading data of the cantilever beam (Fig. 1) are STANDARD DEVIATION OF LOAOKJF1=0.25 L = 5.0, b = h = 1 .O, v = 0.3. The modulus of elasticity is assumed to be Gaussian homogeneous process with mean E = 1.0, and the distributed load is also assumed as Gaussian homogeneous process with mean F = 1.O. All the parameters are represented without any units, so that any units can be specified.

llllll

Ill

ll

ll

l l I IllI 20

I=

L

-x LZS.0 bSh=l.O

Fig. 1. Finite element model with 20 elements (beam with stochastic rigidity under random uniform loading).

1.51

0

I

I

1000

2000 SAMPLE

Fig. 2. Fluctuation

I

3000

1

LOO0

5000

J

SIZE

of the deflection at the tip of the beam with sample size.

44

S. Chakraborty 3,ObhIlOOL

and S. S. Dey

-I

:: Ii 8 ::

2.5

,

8 2 L 0

2.0UF’O.25

a0 2 ,z 5

SAMPLE

I L

I 6

ORDER

should be excluded. For this a trundistribution was used to generate the

values

I 2

SIZE

Fig. 3. Fluctuation of the standard deviation of the deflection at the tip of the beam with sample size.

negative

0.31 0

cated Gaussian random samples. The solution based on Neumann expansion uses eqn (16). The advantage of using eqn (16) is that the matrix inversion is required only once (i.e. inversion of the deterministic part of the stiffness matrix) for all samples. The terms in the expansion series of eqn (16) can be terminated after a few terms, depending on convergence rate. RESULTS AND DISCUSSION

The results of numerical example are presented in Figs 2-7. Figures 4 and 5 show the results for the expected value of deflection and standard deviation of deflection at the tip of the beam vs the order of expansion in the Neumann series. Results of direct simulation are also shown in the same figures. It is observed that the results obtained based on the Neumann expansion method approach those obtained by direct Monte Carlo methods as the order of expansion increases. A notable point is that the rate of convergence is slower in case of standard deviation as standard deviation represents only the deviatoric component, whereas deflection represents the mean

Fig,

5. Convergence

OF

I 6

I 10

EXPANSION

of the standard deflection.

deviation

of

tip

,,Obh/lOOL z I= Y J

0.6-

;

0.6-

-

NEUMANN

*

OIRECT

UF

EXPANSION SIMULATION

-0.1

::

z Oo_1 STANDARD

DEVIATION

OF

E

Fig. 6. Standard deviation of tip deflection of the beam as a function of the standard deviation of the modulus of elasticity.

value plus deviatoric components. It is also experienced that convergence is excellent for relatively lower values of the coefficient of variation. Figures 6 and 7 show the variation of standard deviation of tip deflection of the beam for varying standard deviation of modulus of elasticity and load.

bh/lOOL

2 0.5*

o-0

NEUMANN D

DIRECT

EXPANSION MONTE

P L

METHOD

CARLO

5 h 0 a F

MFTHOO

O.L-

% 5

2.2

-

5 :: o a F:

2.1(TF

= 0.10

2.0 -

0

2

L OROER

Fig. 4. Convergence

6 OF

6

10

EXPANSION

of the tip deflection

of the beam.

STANOARO

DEVIATION

OF

F

Fig. 7. Standard deviation of tip deflection of the beam as a function of the standard deviation of the load.

Spatial distribution of material properties and external loading Corresponding values obtained by direct Monte Carlo simulation are also shown in the same figure. The results agree very well for relatively small coefficients of variation in the stochastic fields. CONCLUDING

2. G. C. Hart and J. D. Collins, The treatment

REMARKS

The response variability of a cantilever beam with stochastic stiffness subjected to random loading was examined in this paper. The advantage of using the Neumann expansion technique is that inversion of the stiffness matrix has to be done only once for all the samples. The convergence of the results from series solution is guaranteed. Sufficiently accurate results can be obtained by considering very few terms of the series for relatively small coefficients of variation of the random fields. However more terms have to be taken for larger coefficient of variation.

8.

9.

10. REFERENCES

11. 1. K. Handa

and G. Kiirrholm, Application of finite element methods in the statistical analysis of structures.

Chalmers struktion.

Tekniska 6 (1975).

Hiigskola

Avd. for

byggnadskon -

12.

45

of randomness in finite element modelling. In National Aeronautic and Space Engineering and Manufacturing Meeting, Los Angeles. Society of Automatic Engineers (1970). S. Nakagiri and T. Hisada, A note on stochastic finite element method (pt 1). Seisan-kenkyu 32, 39 (1988). S. Nakagiri and T. Hisada, A note-on stochastic finite element method (nt 2). Seisan-kenkvu 32, 28 (1988). M. Shinozuka, Monte Carlo solution of structural dynamics. Comput. Struct. 2, 855-874 (1972). F. Yamazaki, M. Shinozuka and G. Dasgupta, Neumann expansion for stochastic finite element analysis. J. Enana Mech. ASCE 114 (1988). F. Yamazaki and M. Shinozuka, Simulation of stochastic fields by statistical preconditioning. J. Engng Mech. ASCE. 116, 268-286 (1990). E. Vanmarcke and M. Grigoriu, Stochastic finite element analvsis of simule beams. J. Engng I _ Mech. ASCE 109 (5) 1203-1214. P. D. Spanos and R. Ghanem, Stochastic finite element analysis of random media. J. Engng Mech. ASCE 115, 1035-1053 (1989). R. Ghanem and P. D. Spanos, Stochastic Finite Element Analysis: a Spectral Approach. Springer, Berlin (1990). 0. C. Zienckiewicz, The Finite Element Method. McGraw-Hill, New York (1977). C. S. Krishnamoorthy, Finite Element Analysis Theory and Programming. Tata McGraw-Hill. New Delhi (1988). -