Stochastic finite element simulation of uncertain structures subjected to earthquake

In present study, the stochastic finite element simulation based on the efficient Neumann expansion technique is extended for the analysis of uncertain structures under seismically induced random ground motion. The basic objective is to investigate the possibility of applying the Neumann expansion technique coupled with the Monte Carlo simulation for dynamic stochastic systems upto that extent of parameter variation after which the method is no longer gives accurate results compared to that of the direct Monte carlo simulation. The stochastic structural parameters are discretized by the local averaging method and then simulated by Cholesky decomposition of the respective covariance matrix. The earthquake induced ground motion is treated as stationary random process defined by respective power spectral density function. Finally, the finite element solution has been obtained in frequency domain utilizing the advantage of Neumann expansion technique.


Introduction
The evaluation of seismic response of structures involves two kinds of stochasticity: one is in the loading conditions and the other is in the structural system itself.The ground motion observed during earthquake is a random process varying both in space and time.Details of such inputs and its effect on structures are studied extensively and are well documented in the field of random vibration.However, there always exist uncertainties in defining the properties of structural parameters.These are evolved from uncertainties in material properties, geometric characteristics and boundary conditions etc. Generally, these parameters are assumed to be described exactly.However the uncertainties associated with these system parameters are also likely to have considerable influence on the response.
The recent development of the stochastic finite element method (SFEM) [1][2][3][4][5] has changed the outlook of alleviating the difficulties in modelling such uncertain system and extensive reviews regarding this are found in Benaroya and Rehak [6], Vanmarcke et al. [7], Manohar and Ibrahim [8].A rational framework for the analysis of uncertain system subjected to stochastic excitation has become feasible [9][10][11].The purpose of the present work is to go beyond the random vibration analysis and to include the uncertainties in structural parameters as well.The application of the stochastic finite element simulation based on the Neumann expansion technique is extended for the analysis of uncertain structure under dynamic random load due to earthquake.

Local averaging method of discretization
In order to discretize the random field, the local averaging method [1,2] has been used.In this method the field variable over an element is approximated by spatial averaging over the element.A homogeneous random scalar field α(x, y) defined in two dimensional plane (x, y) is characterized by its mean, variance and normalized variance function ρ(r x , r y ), where r x = x − x and r y = y − y .The local average of the field over a rectangle A i centred at (x i , y i ) with sides L xi and L yi parallel to the co-ordinate axis x and y respectively is defined as It has been shown [5] that if the correlation of the field is quadrant-symmetric i.e. ρ(r x , r y ) = ρ(−r x , r y ) = ρ(r x , −r y ) = ρ(−r x , −r y ), then the mean vector, variance and covariance matrix of the local average α(x, y) over rectangles A i and A j can be obtained as follows: Where σ = diag[σ 1 , σ 2 , . . .σ N ] and σ i is the standard deviation of the random parameters for i-th discretized element and Here, L xk and L yl (k, l = 0, 1, 2, 3) are the distances characterizing relative positions of the rectangles A i and A j ; and γ(L xk , L yl ) is the normalized variance reduction function of the local average of α(x, y) over the rectangle with sides L xk and L yl as shown in Fig. 1.This function depicts the dependence of the variance of spatial averages on the size of averaging element.Using Equation ( 5), the normalized variance reduction function corresponding to exponential correlation function ρ(r) = exp[−(r/b) 2 ] can be obtained as where b is the correlation parameter and φ(.) is the error function.The value of the error function increases from zero to one as its argument increases from zero to infinity.If α(x, y) is separable, i.e. ρ(r x , r y ) = ρ(r x )ρ(r y ), the two-dimensional reduction factor becomes simply the product of two one dimensional factors i.e. γ(L xk , L yl ) = γ(L xk )γ(L yl ).This one dimensional result can be easily used for computing the variance reduction function of the following exponential covariance function, which has been taken for present numerical study.For a particular analytic expression corresponding normalized variance reduction function can be computed.Various form of analytic expressions i.e. triangular, exponential etc. characterized by correlation parameter have been reported [1].However suitability of any model can be justified by fitting actual experimental data.As no experimental data is available to ascertain the relative merits of alternative models, exponential model has been selected with the purpose of illustrat-ing the analytical procedure.Note that if the field is not homogeneous and or quadrant symmetric or the domain is non rectangular, then Gaussian quadrature can be used [5] to compute the covariance matrix of the random vector field using Equation (5).

Simulation of sample function
If there are 'n' finite elements in the total structure, 'n' property values are then associated with each zero mean homogeneous random field α(x, y).The values α i = α(x i , y i ), i = 1, 2, 3, . . ., n are random with zero mean but correlated.Thus the correlated random vector {α 1 , α 2 , α 3 , . . ., α n } can be obtained as: where, {Z} is a vector consisting of 'n' independent Gaussian random variables with zero mean.
[L] is a lower triangular matrix obtained by Cholesky decomposition of the respective covariance matrix cov(α i , α j ) obtained through local averaging technique.The random vector generated from Equation ( 8), satisfies the original covariance matrix as Once the Cholesky decomposition is over, the sample vector can be easily obtained [12,13] for any number of simulation desired.However, fairly large sample size is required for the simulated covariance matrix approaching the target covariance matrix.In general , the sample size is guided by the degree of accuracy of the resulting ensemble such as mean values and standard deviation.For determining the number of simulation, the standard deviation of desired output can be plotted with number of simulation and sample size may be fixed when the fluctuation is negligible.

Dynamic equilibrium equation
The relative equation of motion of a structure supported on ground and subjected to dynamic forces due to earthquake can be written as [14], In which [M ], [C] and [K] are the global mass, damping and stiffness matrix respectively obtained by assem-bling the corresponding element matrixes and {1} represents the unit column vector.This vector expresses that a unit static translation of the base of a structure produces directly a unit displacement at all the degrees of freedom.{U } is the relative displacement vector, first and second time derivative indicates the velocity and acceleration vector respectively and Üg is the acceleration generated from support motion due to earthquake.

Transfer function
The time domain version of the dynamic equilibrium equation has been transferred to a corresponding frequency domain.The displacement of the linear structural model subjected to a unit amplitude ground acceleration Üg (= e iwt ) can be assumed as: Where H(w) is the frequency response function.Now taking first and second time derivative of Equation (10) and substituting into Equation ( 9), the equation of motion leads to The complex matrix in the bracket term on the left hand side of Equation ( 11) is the dynamic stiffness matrix for the complete structural system.Equation (11) can be represented in a simple form as: where [D] is the dynamic stiffness matrix, {H} is the complex frequency response transfer function and {F } is the forcing vector.The algorithm for solution of the dynamic equation is same as for the static finite element problems except that [D] matrix is complex in nature.So suitable modification is required only for complex matrix operations and the equation has to be solved for each discrete frequency interval to obtain the complex frequency response vector {H}.

Formulation of stochastic element stiffness matrix
The rigidity matrix [D e ] for uncertain modulus of elasticity with deterministic Poisson's ratio can be written separately in the following form: where [C 0 ] is a deterministic matrix and E(x, y) is a random field representing the modulus of elasticity over the domain.It is worth mentioning that the rigidity matrix depends on the nature of spatial variation of the elastic modulus.Various forms of theoretical models are possible [15].In the present study, modulus of elasticity involving homogeneous spatial variation is assumed in the following simple form: where E e is the discretized version of E(x, y) for the eth element which has two parts, the first part is the mean (E e ) and other one is the deviatoric portion (E e α e ).Now the random element stiffness matrix can be spilt up into two parts as where, Here [B] is the standard deterministic strain displacement matrix.Once the mean and deviatoric part of element stiffness matrices are generated, these can be assembled in usual procedure as done in determinis-tic finite element method to obtain the global stiffness matrices.Though, deterministic [B] has been taken in present analysis, it is more elegant to introduce the concept of stochastic shape function [16] in the derivation of strain displacement relationships, particularly for larger stochastic variation.

Formulation of stochastic element mass matrix
The consistent mass matrix of a structural element with random mass density m(x, y) can be written as where N (x, y) is the vector of shape functions.The spatially varying random mass density can be expressed in the following form: Here m e is the discretized version of m(x, y) of e-th element.The random element mass matrix can be also spilt up into two parts as earlier: Where, These element mass matrices are assembled into global co-ordinate system as usual.

Damping matrix
The damping in a structure, in general will be a combination of structural and viscous type.In present study, for simplicity, damping is restricted to structural kind only.The structural damping matrix can be written as [17]: Where i = √ −1 and 'g' is a small constant known as the structural damping factor.Once the stiffness matrix is known, the damping matrix is obtained simply by multiplying the stiffness matrix with structural damping factor.Due to stochasticity in the stiffness matrix, the damping matrix is also uncertain.This can be decomposed into its mean and deviatoric part same as done for stiffness and mass matrices i.e.
[C] = where, It is to mention here that in earthquake analysis, the added mass concept [18] leading to computation of effective mass density of structure is much more realistic.In that case, the natural frequency and mode shape of the structure can be used iteratively for computing new effective mass density.Also the viscous type of damping model which directly use the damping ratio is much more elegant then simple structural type damping modeling.But as the structural parameters are uncertain, the frequency and mode shape of the structure will be uncertain.Thus, to construct the damping matrix and effective mass density computation in such situation the parameter uncertainty has to be included which is another subject matter needs further study.

Direct Monte Carlo solution
In direct Monte Carlo simulation, using the procedure outlined in Sections 3, a set of uncertain sample structures with random elasticity and mass densities are generated.Each of these sample structures is then subjected to the same random harmonic loading as input and respective dynamic responses are evaluated in frequency domain.Finally, the statistics of the responses are evaluated.A direct Monte Carlo simulation is carried out to check the accuracy and efficiency of the proposed Neumann expansion based solution.

Neumann expansion solution
For solving the dynamic equilibrium equation in frequency domain, the dynamic stiffness matrix has to be inverted.The dynamic stiffness matrix defined by Equation ( 12) is random in nature due to spatial variability in modulus of elasticity and mass density.In direct simulation method, the [D] matrix is inverted each time to obtain a sample solution vector {H} which leads to a large amount of CPU time.The Neumann ex-pansion technique has been adopted to circumvent this problem.To apply the Neumann expansion method, the random matrix is decomposed into its respective mean (obtained by replacing the spatial variable parameters with their respective mean values) and fluctuating deviatoric part as: where, The inversion of random dynamic stiffness matrix adopting the Neumann expansion technique [12,19] takes the following form: where, [P ] = Rewriting Equation ( 12) as: the random response vector {H} is obtained by substituting Equation ( 26) in (28), where, The above series solution is equivalent to the solution of the following recursive equation: where r denotes the r-th simulation.Once the inversion of the mean dynamic stiffness matrix [D] is obtained, Equation ( 29) can be used iteratively to obtain the ran-dom response for each simulated structure without further matrix decomposition.The expansion series in Equation ( 29) may be terminated after few terms depending on the convergency rate and accuracy requirements for a particular problem.The solution procedure using the Neumann expansion is same as static problems [19].However, here the dynamic stiffness matrix D is complex in nature.But the static problems involve only stiffness matrix K which is real.Hence the essence of the present work is to investigate the applicability and efficiency of the Neumann expansion based simulation technique for stochastic dynamic systems i.e. when the dynamic stiffness matrix D is complex in nature.

Random ground motion
The seismically induced ground motion is considered to be a stationary random process (with constant frequency content) truncated for a finite duration.The following power spectral density(psd) function suggested by Kanai and Tajimi is used to characterize a strong earthquake ground motion: Where ω is the component frequency, S 0 is the scale factor, ω g and ξ g are the equivalent natural frequency and damping ratio of the ground characterized by a single degree of freedom system.

Frequency response function
Ensemble of frequency response functions for various simulated sample structures are obtained from Equation (29).Now, the expected value and covariance of the frequency response function can be computed from the following equations:

Spectral density function
For linear system, with known transfer function of response, the spectral density function S u for any response variable u can be readily obtained as: Where H represents the transfer function vector for the response u and H * T is its complex conjugate transpose.The expected value and covariance of S u are obtained same as frequency response function.

Mean square displacement
The i-th moments of the spectral density function S u are For i = 0, the above integral gives the value of mean square displacement i.e.
The approximate integration has been performed using Simpson's rule up to the cut-off frequency: Where 'n' is odd and ∆ω is the frequency interval.

Numerical example
The finite element model of a plate fixed at base is shown in Fig. 2. The nodes along the lower edge are fixed and the base is excited due to random ground motion transverse to the plane of the plate.The details of the geometric and material properties of the plate are: The most important aspect of stochastic analysis is the representation of uncertain elastic constants, mass density etc.It is frequently suggested that in absence of any data for fitting probability distribution function for structural parameters, it is acceptable in design to assume Gaussian distribution [20], though non-Gaussian field models for material properties can also be found [21].The Gaussian models have limitations for cases where the parameters experience large variations or when reliability issues are being examined.Again assumption of Gaussian distribution implies the possibility of generating negative values of elastic properties or mass density.However, this difficulty can be alleviated by using truncated Gaussian distribution [12,22].For illustrating the analytical procedure following Gaussian models of elasticity and mass density are assumed in the present work is not too large.In general, c.o.v.larger than 0.2 shows a large deviation and needs to consider too many terms in the Neumann expansion series and advantage of the proposed method is lost.Figure 9 shows the normalized CPU time in estimating the response quantities by Neumann expansion with increasing number of terms in expansion series.However, time saving depends on the order of expansion, numbers of degrees of freedom in the system, variability of the input and required degree of accuracy.It has been experinced that if the order of expansion is restricted and not too large, the method will be advantageous in terms of CPU time, as the number of degree of freedom increases.

Conclusions
By expanding the uncertain dynamic stiffness matrix about its reference values, the Neumann expansion method can be readily introduced in the finite element procedure to solve the system stochasticity problems.The factorization of the stiffness matrix is to be performed only once for entire simulation in the Neumann expansion method leads to a considerable savings in CPU time.In static analysis the convergency of the Neumann expansion series solution [15,19] is guaranteed for each simulated sample structure and can be confirmed numerically irrespective of c.o.v. of the random input parameters, but this is not the case for dynamic systems with large stochastic dispersion.In particular, around the natural frequency of the structure, the psd of response results slow up leading to erroneous results even with too many number of terms in Neumann expansion series.In general, the Neumann expansion solution gives converged solution upto 0.2 c.o.v of uncertain input parameters.However, beyond 15% variatition of random input parameters needs too many terms to consider in the Neumann expansion series and advantage of the proposed method is lost.

Fig. 1 .
Fig. 1.Definition of distances characterizing relative positions of rectangles A i and A j .

Fig. 2 .
Fig. 2. Finite element model of the plate fixed at ground subjected to earthquake loading transverse to the plane of the plate.

Fig. 7 .
Fig. 7. Comparison of expected value of rms deflection at node 46 with varying c.o.v. of uncertain parameters.
L x = 9.0 m, L y = 9.0 m, h = 0.2 m, mean value of the modulus of elasticity, E = 2.0 × 10 10 N/m 2 , mean value of mass density, m = 2400.0Kg/m 3 , Poisson's ratio ν = 0.15, shear correction factor (to consider shear deformation effect) = 5/6.The responses have been computed at a frequency interval of 0.5 and cut off frequency is taken as 50.0 rad/sec.

Fig. 8 .
Fig. 8.Comparison of standard deviation of rms deflection at node 46 of the plate with varying c.o.v. of uncertain parameters.E(x, y) = E[1 + C e E 1 (x, y)] (39) where, C m and C e are the coefficient of variation (c.o.v.) of mass density and modulus of elasticity respectively.The stochastic fields m 1 (x, y) and E 1 (x, y) represents the deviatoric component of the mass density and elasticity respectively assumed as zero mean and unit variance Gaussian process with normalized covariance function as described by Equation (7).The correlation parameter b x and b y are taken as 3.0.The spectral density function of the random ground motion is defined by Equation (32) with S 0 = 3.73×10 −3 m 2 /s 3 , ω g = 15.56 rad/sec, ξ g = 0.64 correspond to a strong earthquake motion.The numerical results obtained by the Neumann expansion method are compared with those from the direct Monte Carlo simulation with respect to accuracy and computational efficiency.It is found that the fluctuation of standard deviation (s.d.) of response is negligible after 100 simulations and with this sample size, nature and trend of various results can be well demonstrated.The expected value and standard deviation (s.d.) of power spectral density (psd) of maximum deflection shown in Figs 3 to 6 as a function of frequency

Fig. 9 .
Fig. 9. Comparison of normalized CPU time of the plate problem.