Uncertain linear structural systems in dynamics: Efficient stochastic reliability assessment

Computers and Structures 88 (2010) 74–86

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Uncertain linear structural systems in dynamics: Efficient stochastic reliability assessment H.J. Pradlwarter, G.I. Schuëller * Institute of Engineering Mechanics, University of Innsbruck, A-6020 Innsbruck, Austria

a r t i c l e

i n f o

Article history: Received 2 April 2009 Accepted 14 June 2009 Available online 25 July 2009 Keywords: Finite element models Linear structures Dynamics Gaussian excitation Uncertain structural parameters Reliability

a b s t r a c t A numerical procedure for the reliability assessment of uncertain linear structures subjected to general Gaussian loading is presented. In this work, restricted to linear FE systems and Gaussian excitation, the loading is described quite generally by the Karhunen–Loève expansion, which allows to model general types of non-stationarities with respect to intensity and frequency content. The structural uncertainties are represented by a stochastic approach where all uncertain quantities are described by probability distributions. First, the critical domain within the parameter space of the uncertain structural quantities is identified, which is defined as the region which contributes most to the excursion probability. Each point in the space of uncertain structural parameters is associated with a certain excursion probability caused by the Gaussian excitation. In order to determine the first excursion probability of uncertain linear structures, an integration over the space of uncertain structural parameters is required. An extended procedure of standard Line sampling [P.S. Koutsourelakis, H.J. Pradlwarter, G.I. Schuëller, Reliability of structures in high dimensions, part I: algorithms and applications, Probabilistic Engineering Mechanics 19(4) (2004) 409–417; G.I. Schuëller, H.J. Pradlwarter, P.S. Koutsourelakis, A critical appraisal of reliability estimation procedures for high dimensions, Probabilistic Engineering Mechanics 19(4) (2004) 463–474] is used to perform the conditional integration over the space of uncertain parameters. The suggested approach is applicable to general uncertain linear systems modeled by finite elements of arbitrary size by using modal analysis as exemplified in the numerical example. Special attention is devoted to the efficiency of the proposed approach when dealing with realistic FE models, characterized by a large number of degrees of freedom and also a large number of uncertain parameters. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction The potentials and abilities of modeling the behavior of fluids, solids and complex materials subjected to various forces and boundary conditions by finite element analyses (FEA) provides a key technology for further developments in the industrialized world. Supported by the available relatively inexpensive and continuously growing computer power, the expectations on FEA are shifting towards reliable and robust computational simulations and predictions of physical events. For achieving this goal, answers to many open issues in computational mechanics need to be provided as discussed e.g. in [10]. The need to design, given uncertain material properties, production processes, operating conditions and fidelity of mathematical computational (FE) models to represent reality, leads to the concern of reliability and robustness of computer-generated predictions. Without some confidence in the validity of simulations, their value is obviously diminished. Since * Corresponding author. E-mail address: [email protected] (G.I. Schuëller). 0045-7949/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2009.06.010

uncertainty is always present in non-trivial realistic applications, uncertainty propagation through the FEA is one of the important issues which must be addressed for further developments. In this paper, a computational efficient reliability estimation procedure for uncertain linear systems subjected to dynamic stochastic loading is presented. The approach is designed to cope with large FE-models in terms of degrees of freedom, a large number of uncertain input quantities and high variabilities in terms of coefficient of variation (e.g. P10%). Uncertainties with respect to dynamic loading and of the parameters describing the mechanical properties of the FE-model are propagated by a stochastic approach: Inherent irreducible (aleatory) uncertainties as well as reducible (epistemic) uncertainties due to insufficient knowledge or modeling capabilities are translated into a probability distribution defining the input of the stochastic analysis. The reliability will be assessed in terms of the first excursion probability of critical responses. For this purpose, Line sampling [7,15] is further developed and extended. Efficient solutions of the first excursion probability for deterministic systems are exploited to compute failure probabilities conditioned on discrete realizations of the uncertain

75

H.J. Pradlwarter, G.I. Schuëller / Computers and Structures 88 (2010) 74–86

structural properties. Computationally efficient procedures to determine the stochastic response in terms of the Karhunen–Loève representation of the conditional critical response are suggested by using modal analysis, impulse response functions in combination with Fast Fourier Transforms (FFT). The associated conditional first excursion probabilities are evaluated by Line sampling. The shown procedure is applicable for any Gaussian distributed excitation represented by a Karhunen–Loève expansion and for an arbitrarily large FE-model. In a further step, the domain of uncertain structural parameters which contributes most to the failure probability is identified by a recently developed estimation procedure which allows to identify the most influential uncertain parameters among a large number. Finally, the unconditional total failure probability is estimated by a novel Line sampling procedure, covering the important failure domain and efficiently integrating over the whole parameter uncertainty space. To demonstrate the applicability of the proposed approach for general FE-models, a realistic FE-model for a twelve story building with more than 24,000 DOF’s and 200 uncertain quantities is analyzed.

2. Methods of analysis 2.1. General outline This section, containing the theoretical developments, is subdivided into six subsections. The first subsection addresses the Karhunen–Loève representation of the stochastic excitation. The second describes the treatment of uncertain structural properties by a stochastic approach. The third considers impulse response functions and its use to evaluate the critical response for linear systems conditioned on specific realizations of the uncertain structural properties. The fourth subsection shows how Line sampling can be applied to estimate the conditional first excursion probability. In the fifth subsection, the domain within the parameter space of uncertainties is identified which contributes most to the first excursion probability. In the last and sixth subsection, the integration over the entire high dimensional parameter space is carried out by a novel Line sampling procedure. 2.2. Representation of uncertain excitation Dynamic excitations acting on the structure are in many cases uncertain. However, although it might be impossible to describe these excitations in a deterministic sense, information on the expected range and its variability is usually available. Such information can be described by the mean value, the standard deviation of the fluctuation, and also the correlation in space and time. The dynamic excitation f ðx; a; AðtÞÞ is a function of the spatial coordinates x, the direction a and the amplitude AðtÞ as a function of time t. This complex dependency on eight scalar quantities fx1 ; x2 ; x3 ; a1 ; a2 ; a3 ; A; tg is considerably simplified by the use of finite element analysis (FEA). In FEA, the spatial coordinates x and the direction of actions a are specified by the degrees of freedom and all forces are reduced to nodal forces. The excitation f can therefore be interpreted as vector with m components, where the structure is assumed to be specified by m degrees of freedom. Hence it suffices in FEA, and as implied here, to specify the dynamic excitation by the m-dimensional vector f ðtÞ as a function of time t. The uncertainties of the excitation f ðtÞ are described mathematically by a stochastic process [8,1,16] characterized as a function of independent random variables N ¼ ðN1 ; N2 ; . . . ; Nn Þ and by deterministic (vector) functions f ðiÞ ðtÞ of time t. Most conveniently, each

of the independent random variables is assumed to follow a standard normal distribution

! 1 n2i qNi ðni Þ ¼ pffiffiffiffiffiffiffi exp  2 2p Z ni P½Ni 6 ni  ¼ qNi ðxÞdx ¼ Uðni Þ;

ð1Þ ð2Þ

1

in which UðÞ denotes the cumulative standard normal distribution. Any Gaussian distributed process, is most conveniently described by the so called Karhunen–Loève presentation (see e.g. [9,6,5,14]).

f ðt; nÞ ¼ f

ð0Þ

ðtÞ þ

n X

ni f

ðiÞ

ðtÞ:

ð3Þ

i¼1

In the above, all vectors on the right hand side do have deterministic properties and the independent random variables assume the following relations,

E½n2i  ¼ 1;

E½ni  ¼ 0;

E½ni nj  ¼ 0 for i–j;

ð4Þ

where E½ denotes the mean or expectation. The representation (3) specifies uniquely the mean lf ðtÞ and the variance r2fk ðtÞ or standard deviation rfk ðtÞ for each degree of freedom k, and the correlation of the uncertain excitation.

lf ðtÞ ¼ f ð0Þ ðtÞ n h i2 X ðiÞ r2fk ðtÞ ¼ fk ðtÞ

ð5Þ ð6Þ

i¼1

E½fj ðt1 Þfk ðt 2 Þ ¼

n X

ðiÞ

ðiÞ

fj ðt 1 Þfk ðt2 Þ:

ð7Þ

i¼1

However, the deterministic vector valued functions f ðiÞ ðtÞ are usually not quantified a priori. They need to be determined from the symmetric covariance matrix C f ðtj ; t k Þ. Assuming that the excitation is discretized at equidistant instants 0 P tk ¼ kDt P T, the covariance matrix is defined as given below,

C f ðt j ; tk Þ ¼ E½ðf ðtj Þ  lf ðtj ÞÞðf ðt k Þ  lf ðtk ÞÞ0 ;

0 6 t 1 ; t 2 6 T;

ð8Þ

0

where a prime ‘‘ ‘‘ denotes the transposed vector and T is the considered duration. After solving the eigenvalue problem,

C f U ¼ UK;

ð9Þ

where U contains the eigenvector and K is a diagonal matrix of eigenvalues, the eigen-solution for the n highest eigenvalues k1 P k2 P    P kn are used. The deterministic vector values functions f ðiÞ ðtÞ are then determined by

f

ðjÞ

ðtk Þ ¼

pffiffiffiffi kj /½kj ;

ð10Þ

where ½k denotes the rows associated with t k and j the column of the eigenvector matrix. In practice it would be extremely difficult, if not unfeasible, to establish the covariance matrix C f for all m degrees of freedom simultaneously and all K time steps tk ; k ¼ 1; . . . ; K. This would lead to a quadratic matrix of the size m  K, which might be beyond a manageable size. For a tractable solution, independent excitation processes are considered separately. A typical example for such a separation could be an earthquake excitation in a specified direction,

f ðtÞ ¼ MI  aðtÞ

ð11Þ

in which M is the mass matrix and I is vector with 1’s for all degrees of freedom in the considered direction and zeros elsewhere. Hence the nodal forces f ðtÞ are for any time t fully correlated, and only the correlations at different times t j – t k need to be described. Hence, after a discretization of the time t by t k ¼ kDt; k ¼ 1; . . . ; K,

76

H.J. Pradlwarter, G.I. Schuëller / Computers and Structures 88 (2010) 74–86

the covariance matrix of the scalar acceleration aðt k Þ has only the size K, for which the eigen-solution poses no difficulty, leading to the representation

aðtk Þ ¼ að0Þ ðt k Þ þ pffiffiffiffi a ðt k Þ ¼ ki Uki "

n X

ni aðiÞ ðtk Þ

ð12Þ

i¼1

ðiÞ

f ðt k Þ ¼ MI  að0Þ ðt k Þ þ

n X

ð13Þ

# ni aðiÞ ðt k Þ :

ð14Þ

i¼1

In practical applications, each of the deterministic Karhunen–Loève excitation terms f ðiÞ ðtÞ can be described by a constant vector F ðiÞ defined over the range of degrees of freedom and a scalar function ðiÞ h ðtÞ as function of time t,

f

ðiÞ

ðiÞ

ðtÞ ¼ F ðiÞ  h ðtÞ:

ð15Þ

As shown above, such a representation applies for earthquake excitation, with F ¼ MI. 2.3. Uncertain structural systems

H ¼ fH1 ; H2 ; . . . ; HS g

ð16Þ

such that for any fixed set H ¼ h, these structural matrices are uniquely specified. Hence, for any fixed set (vector) h, the structural matrices MðhÞ; DðhÞ and KðhÞ can be treated deterministically. It is common practice, to assume the random variables Hi to be standard normally distributed, i.e. with zero mean E½Hi  ¼ 0 and unit variance E½H2i  ¼ 1. This assumption, however, does not imply that the structural components follow a Gaussian distribution or that correlations among structural parts do not exist. Non-Gaussian distributions are realized by non-linear relations of these basic standard normally distributed random variables. Correlations among different structural parts are established by representing different parts as a functions of one or several random variables. 2.4. Impulse response functions Since the structural system is assumed to be linear, the law of superposition is valid. This law implies, that the response uðt; hÞ for fixed structural properties h and dynamic excitation f ðtÞ has, analogous to (3), also a Karhunen–Loève representation, n X

ni uðiÞ ðt; hÞ;

ð17Þ

i¼1

where each deterministic term uðiÞ ðt; hÞ; i ¼ 0; 1; . . . ; n, is the solution of a deterministic dynamic analysis, involving the constant symmetric mass matrix MðhÞ, damping matrix DðhÞ and stiffness matrix KðhÞ in the equation of motion.

€ ðiÞ ðt; hÞ þ DðhÞu_ ðiÞ ðt; hÞ þ KðhÞuðiÞ ðt; hÞ ¼ f MðhÞu

80 6 i 6 n:

v_ ðiÞ ð0Þ ¼ M 1  F ðiÞ :

ð19Þ

For illustration, Fig. 1 shows an example of the impulse response function for the velocity and displacement response of a 2-DOF system. Since the system is linear, also efficient modal analysis applies. After solving the eigenvalue problem,

KðhÞUðhÞ ¼ MðhÞUðhÞKðhÞ the impulse response zðt; hÞ

v ðt; hÞ

ð20Þ is represented in modal coordinates

v ðiÞ ðt; hÞ ¼ UðhÞzðiÞ ðt; hÞ:

ð21Þ

Then, the modal coordinates have initial values ðiÞ z_ j ð0; hÞ is specified by the vector

ðiÞ zj ð0; hÞ

¼ 0, and

z_ ðiÞ ð0; hÞ ¼ UT ðhÞF ðiÞ :

ð22Þ

The free motion in modal coordinates has the explicit solution [8]

Linear structural systems are described by a constant symmetric mass matrix M, damping matrix D and stiffness matrix K. These matrices are realistically considered to be associated with uncertainties. In the proposed stochastic parametric approach [1], the uncertainties of these matrices are also conveniently modeled as functions of independent random variables

uðt; h; nÞ ¼ uðhÞð0Þ ðtÞ þ

function v ðiÞ ðtÞ and the Fast Fourier transform for computing the convolution. The impulse response function is computed as free motion with initial zero displacement and the initial velocity:

ðiÞ

ðiÞ

zj ðt; hÞ ¼

ðiÞ z_ j ð0; hÞ

x0j

efj xj t sinðx0j tÞ

ðiÞ ðiÞ z_ j ðt; hÞ ¼ z_ j ð0; hÞefj xj t cosðx0j tÞ 

fj xj t €zðiÞ _ ðiÞ j ðt; hÞ ¼ zj ð0; hÞe

Hence, to specify the variance of the displacement response, n þ 1 deterministic analyses are required. However, the dynamic response might depend considerably on the structural parameters specified by the vector h. One of the efficient ways to compute the associated dynamic displacement response uðiÞ ðtÞ is the use of the impulse response

ð23Þ f j xj

x0j

  sin x0j t

# ð24Þ

  2fj x0   x2j  j 0 ð1  2f Þ sin x t þ cos x0j t ; j j x0j xj ð25Þ

where all the quantities fj ; xj and x depend actually on the parameter h: 0 j

qffiffiffiffiffiffiffiffiffiffiffiffiffi ðiÞ kj ðhÞ; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x0j ¼ xj 1  f2j ðhÞ:

xj ¼

ð26Þ ð27Þ

The damping ratios fj ðhÞ are, most advantageously, not specified by the an explicit damping matrix DðhÞ, but by modal damping ratios with a moderate increase with frequency xj . The damping ratios might be defined directly as function of uncertain quantities, e.g. fj ðhÞ ¼ cj  expðaj hk Þ with given value for cj and aj . Suppose, the critical responses of interest, e.g. stresses, strains, accelerations, displacements, etc., are comprised in the vector yðt; hÞ. Each component yi ðt; hÞ of the vector must fulfill certain conditions in order to be regarded as reliable. It will be assumed that these critical types of structural response can be represented by a linear combination of the displacement or acceleration re€ ðt; hÞ, sponse uðt; hÞ or u

yðt; hÞ ¼ Q ðhÞuðt; hÞ;

ð28Þ

where Q is a constant matrix, independent of time t, and a function of the structural parameters which might depend on h. Similarly, as the displacement response, the variability due to the random excitation can be cast into a Karhunen–Loève representation.

ðtÞ; ð18Þ

"

yðt; h; nÞ ¼ yð0Þ ðt; hÞ þ

n X

ni yðiÞ ðt; hÞ:

ð29Þ

i¼1

The impulse response function associated with F ðiÞ is then:

wðiÞ ðt; hÞ ¼ Q ðhÞUðhÞzðiÞ ðt; hÞ: ðiÞ yk ðt; hÞ

The critical response due to the excitation term f therefore the result of the convolution

ð30Þ ðiÞ

ðtÞ is

77

H.J. Pradlwarter, G.I. Schuëller / Computers and Structures 88 (2010) 74–86

Velocity impulse response

mass 2

0.5

mass 1

0 −0.5

displacement [m]

mass 1 mass 2

1

velocity [m/s]

Displacements impulse response 0.2

mass 1 mass 2

0.1 0 −0.1

−1 0

2

4

6

8

−0.2

10

0

2

time t

4 6 time [s]

8

10

Fig. 1. Impulse response function.

ðiÞ

yk ðt; hÞ ¼

Z 0

t

ðiÞ

ðiÞ

h ðsÞwk ðt  s; hÞ ds:

ð31Þ

Numerically, this integral is most efficiently computed by using Fast Fourier Transforms, ðiÞ

ðiÞ

h ðtÞ ! FFT ! h ðxÞ;

ð32Þ

ðiÞ ðiÞ wk ðtÞ ! FFT ! wk ð Þ; ðiÞ ðiÞ ðiÞ yk ð Þ ¼ h ð Þwk ð Þ; ðiÞ ðiÞ yk ð Þ ! IFFT ! yk ðtÞ;

ð33Þ

x x

x

x x

 pi;f ðhÞ ¼ P½maxfyi ðt; hÞg P b i 06t6T

[ min fyi ðt; hÞg 6 bi  06t6T Z qðnÞdn ¼

where FFT denotes the discrete Fast Fourier Transform and IFFT its inverse. For the efficiency of the presented approach, it should be stressed that a single modal FE analysis is sufficient to compute the Gaussian distributed response. The associated Karhunen–Loève representation (29) provides the basis for the reliability assessment conditional on the structural parameters h. 2.5. Conditional reliability 2.5.1. First excursion probability In this section, the conditional first excursion problem pf ðN; hÞ is discussed for random Gaussian excitation and for the case where a random realization of H assumes a certain set of deterministic parameters h. The first excursion probability is then defined as the probability that the critical response exceeds at least once within the considered time period ½0; T the threshold bi of the ith component of the critical response yðt; hÞ [8]. The threshold bi  as shown might consist of a lower limit bi and an upper limit b i in Fig. 2 for a single critical response.

ð37Þ

g i ðh;nÞ60

  maxfy ðt; h; nÞg; g i ðh; nÞ ¼ minfb i i 06t6T

ð38Þ

min fyi ðt; h; nÞg  bi g:

ð34Þ ð35Þ

ð36Þ

06t6T

This conditional first excursion probability corresponds to the case of stochastic excitation, specified by the random vector N and its probability density function qðnÞ, subjected to a deterministic strucR tural system. In the next section, the integration pf ¼ pf ðhÞqðhÞ dh over the whole domain of the probability density function qðhÞ will be discussed in detail, leading to the unconditional first excursion problem pf . To simplify the notations and the reliability problem, each critical response will be considered separately. Hence, the subscript ‘‘i” is dropped in the following developments. For this case, where the critical response yðtÞ can be represented by a Karhunen–Loève representation, as given in (29), several efficient numerical procedures have been developed recently to calculate the first excursion probability (e.g. [2,13]). 2.5.2. Line sampling in dynamics In this section, a procedure, based on Line Sampling (LS), is presented. Line sampling [7,15] is particularly efficient if a so called important direction can be computed, which points towards the failure domain nearest to the origin (see Fig. 3). Most importantly, it is not required that the vector a (see Fig. 3) points exactly towards the center of this domain, nor are any assumptions made

 Fig. 2. First Excursions of thresholds b and b.

78

H.J. Pradlwarter, G.I. Schuëller / Computers and Structures 88 (2010) 74–86

2

  1 c ðjÞ If n½i ðcÞ pffiffiffiffiffiffiffi exp  dc 2 2p 1     ðjÞ ðjÞ ; ¼ U b½i þ U b ½i

  Z ðjÞ pf n~½i ; h ¼

þ1

ð43Þ

where UðÞ denotes the cumulative standard normal   distribution, ðjÞ n½i ; h denotes the If ðÞ is an indicator function of failure, and pf ~ failure probability conditioned on the jth randomly selected line ðjÞ ðjÞ n½i þ ca½i and on the strucin the ith important direction n½i ðcÞ ¼ ~ tural parameter set h. The safe domain lies within the bounds h i ðjÞ ðjÞ b½i ; b ½i . Since the superposition law is valid, the response for ðjÞ ?ðjÞ n½i ¼ ca½i þ n½i ðcÞ can be specified by the following function:

     ðjÞ ?ðjÞ y t; c; n½i ¼ cy t; a½i þ y t; n½i  y t; a½i ¼

Fig. 3. Line sampling in high dimensions.

n X

ð44Þ

a½ik yðkÞ ½i ðt; hÞ

ð45Þ

n   X ?ðjÞ ?ðjÞ ðkÞ ¼ yð0Þ ðt; hÞ þ nk½i y½i ðt; hÞ: y t; n½i

ð46Þ

k¼1

k¼1

regarding the shape of the limit state function gðh; nÞ ¼ 0 [7]. LS is robust, since unbiased reliability estimates are obtained, irrespectively of the important direction a, and any deviation from the optimal direction merely increases the variance of the estimate. The basic procedure of LS is for completeness shown in the Appendix A. For structures with the properties specified uniquely by h and which are excited by a Gaussian distributed loading process, the important direction is well known to be proportional to the response yðjÞ ðt; hÞ. Let cðt; hÞ be the standard deviation of the response yðt; hÞ

c2 ðt; hÞ ¼ Var½yðt; hÞ ¼

n X ½yðjÞ ðt; hÞ2 :

ð39Þ

h i ðjÞ  ðjÞ along the jth random line in the ith Then, the safe domain b½i ; b ½i important direction is defined as ðjÞ b½i

 ðjÞ b ½i

j¼1

Then, the unit vector aðt; hÞ pointing towards the important direction is specified by

aj ðt; hÞ ¼

yðjÞ ðt; hÞ cðt; hÞ

ð40Þ

It should be noted that the important direction is a continuous function of time t. Moreover, linearity implies that any excitation defined by n and orthogonal to aðt; hÞ, will have a zero response at time t. Since the important direction aðt; hÞ varies within the time interval 0 6 t 6 T, and the reliability index bðt; hÞ might be of similar size over a considerable portion of time, the efficiency and reliability of Line sampling  can L be increased by sampling in several important directions a½i i¼1 . The algorithm to determine these  ½i L important directions a i¼1 is presented in Appendix B. The procedure involves a transformation of the KL vectors yðjÞ ðt; hÞ such  L that a½i ¼ ei i¼1 , where ei are unit vectors of which the ith component contains the value one and zeros elsewhere. In LS, each point n in the standard normal space is decomposed into the one dimensional space ca½i and the ðn  1Þ dimensional subspace n? ða½i Þ orthogonal to the direction a½i ; i ¼ 1; 2; . . . ; L, ½i

?

½i

n ¼ ca ðhÞ þ n ða Þ:

ð41Þ

  8   y t; n?ðjÞ
   H y t; a½i yðt; a½i Þ   9 ?ðjÞ = b  y t; n½i

 ½i  þ  H y t; a ; 0 6 t 6 T ; yðt; a½i Þ   8   y t; n?ðjÞ
  ½i ¼ min  H y t; a½i : yðt; a½i Þ   9 ?ðjÞ = b  y t; n½i

 ½i  þ  H y t; a ; 0 6 t 6 T ½i ; yðt; a Þ ¼ max

:

ðjÞ

ðjÞ

ðjÞ

 ; b½i < nk ða½i Þ < b ½i

k – i and 0 < k 6 L:

^ðjÞ p f ðhÞ ¼

L X

ðjÞ

pf ½i ðhÞ

ð50Þ

i¼1 ðjÞ

pf ½i ðhÞ ¼ ^f ðhÞ ¼ p

Ni 1 X ðkÞ p ð~n ; hÞ Ni k¼1 f ½i

N 1X ^ðjÞ ðhÞ p N j¼1 f

~nðjÞ ¼ nðjÞ  ða½iT  nðjÞ Þa½i : ½i

2.6. Design point for stochastic structural systems

ðjÞ

ð49Þ

For the very unlikely cases in which the above condition is not satisfied, the line will be ignored and another line is generated instead. The independent estimates pf ðhðjÞ Þ for the failure probability allow to compute an unbiased estimate for the conditional failure probability

and the variance

For each independent random realization ~ n½i , the probability of failðjÞ ure conditioned on the value ~ n½i and h can be determined by

ð48Þ

in which HðyÞ ¼ 1 for y P 0 and HðyÞ ¼ 0 for y < 0. To each of the important directions fa½i gLi¼1 , a disjunct failure domain is associated by imposing the condition

 ðjÞ N In the following it is assumed that ~ n j¼1 are independent sample ? ½i points of the subspace n ða Þ drawn by direct Monte Carlo simulation. The random sample ~ nðjÞ is generated by simulating first n independent standard normally distributed components of a vector  nðjÞ and then subtracting the component which points towards the direction of a½i ,

ð42Þ

ð47Þ

Varðpf ðhÞ ¼

N  2 X 1 ðjÞ ^f ðhÞ ; pf ðhÞ  p NðN  1Þ j¼1

ð51Þ ð52Þ

ð53Þ

which is the basis for deriving further confidence intervals.

In this section, the domain which contributes most to the unconditional total failure probability pf will be determined. The-

79

H.J. Pradlwarter, G.I. Schuëller / Computers and Structures 88 (2010) 74–86

oretically, the failure probability pf jh , conditioned on specific realizations h, need to be integrated over the whole domain of pH ðhÞ

@b2 ðn; hÞ ¼ 0; @hk

pf ¼

which leads to the solution

Z

pf ðhÞpH ðhÞ dh:

ð54Þ

Since the number of uncertain structural parameters, characterized by the random quantities fH1 ; H2 ; . . . ; HS g, is in general large, numerical integration is not feasible. Basically the above integral can be estimated by Monte Carlo sampling

^f ¼ p

N 1 X p ðhðjÞ Þ; N j¼1 f

ð55Þ

hk ¼

where N independent realizations h are drawn from the distribution pH ðhÞ and the associated conditional failure probability pf ðhðjÞ Þ is determined as described in the previous section. Such an approach is only efficient in case the variability of pf ðhðjÞ Þ is small, i.e. might be within one order of magnitude. However, in case the structural variability is high, the conditional first excursion probabilities will vary over many orders of magnitude and direct Monte Carlo approaches will be very inefficient. Considering Eq. (54), it is not difficult to recognize that sampling in the neighborhood of the structural design point h , characterized by

pf ðh ÞpH ðh Þ P pf ðhÞpH ðhÞ;

ð56Þ

contributes most to the total (unconditional) failure probability. In order to identify this domain, composed of the standard normal variables fn; hg, the point fn ; h g with the minimal distance b to the failure surface is determined at the most critical time t 1 (see (B.8)). Since the uncertainty in the excitation and of the structural parameters, specified by n and h , respectively, are independent and therefore orthogonal in the standard normal space, this distance is specified according to Pythagoras by

b2 ðn; hÞ ¼ knk2 þ khk2 :

ð57Þ

Fig. 4 shows a sketch of the failure domain for this case, where bðn; hÞ can only be determined in the complete space ðn; hÞ. The threshold b is reached for 2

b knk2 ¼ 2  c ðt ; hÞ   yð0Þ ðt ; hÞ; yð0Þ ðt ; hÞ  b: b ¼ min½b

ð58Þ ð59Þ

The failure point ðh ; n Þ, with the highest probability density and defined as the nearest failure point to the origin in standard normal space, is derived by imposing the necessary condition

2

3 ðt  ; h Þ

c



ð60Þ

@ cðt  ; h Þ : @hk

ð61Þ

An accurate solution of the above relation can only be obtained in an iterative manner, where convergence is reached for s > S with ðSÞ smaller than the tolerance, ðsþ1Þ

hk ðjÞ

b

k ¼ 1; 2; . . . ; K;

ðsÞ

¼ ð1  wÞhk

þw

2

@ cðt ; hðsÞ Þ @hk c3 ðt ; hðsÞ Þ b

ð62Þ



ðsþ1Þ ¼ khðsþ1Þ  hðsÞ k=khðsþ1Þ k;

ð63Þ

where w ¼ 0:5 leads to a stable and fast convergence. However, there is not much gain in evaluating h very accurately. A first or second step ðs ¼ 0; 1Þ is usually sufficiently accurate for the integration procedures as developed in the next section. The above iterative evaluation of the structural design point h requires a gradient computation which might hamper the efficiency of the approach in case the number of uncertain structural parameters is large. For such cases, a gradient estimation procedure – as shown in the Appendix C – might be used to improve the computational efficiency. 2.7. First excursion probability for stochastic systems In this section, the procedure to estimate the total unconditional first excursion probability by integrating over the complete space of uncertain structural parameters is shown (see Eq. (54)). The proposed approach is to limit the application of LS to the subspace of the structural parameters h. This in fact leads to robust results, even for quite large uncertainties of the structural parameters, since the important directions a½i ðhÞ are then always computed in the optimal directions. The procedure is outlined as follows: 1. Determine the design point h with acceptable accuracy and compute

bS ¼ kh k h aS ¼ bS

ð64Þ ð65Þ

where the index ‘‘S” denotes the subspace of the structural uncertainties. 2. Generate samples fh?ðjÞ gNj¼1 using direct MCS in the subspace of h, which are perpendicular to the vector h . 3. Compute for each parallel line for, say, five discrete points hðjÞ ðck Þ; k ¼ 1; . . . ; 5, as shown in Fig. 5

hðjÞ ðck Þ ¼ h?ðjÞ þ ck aS ; ck ¼ bS þ ðk  3ÞD;

ð66Þ k ¼ 1; 2; . . . ; 5

ð67Þ

D  0:6

ð68Þ ðjÞ

the associated conditional failure probability pf ðck Þ. ðjÞ 4. Estimate the conditional failure probability pf along each line,

1 ðjÞ pf ¼ pffiffiffiffiffiffiffi 2p

Z

1

1

ec

2 =2

ðjÞ

pf ðcÞ dc:

ð69Þ ðjÞ

Fig. 4. Structural design point h .

with quadratic interpolation of the function lnðpf ðcÞÞ by using the available discrete values. This estimation will be discussed subsequently in more detail. 5. Estimate the mean and variance of the failure probability pf ðtÞ according to Eqs. (52) and (53).

80

H.J. Pradlwarter, G.I. Schuëller / Computers and Structures 88 (2010) 74–86 ðiÞ

due to random fluctuation of the estimates pf ðcl Þ. Experience showed that the results using either a linear or quadratic approximation are quite similar. However, a linear approximation is less ðiÞ sensitive to random fluctuations of the estimate pf ðcl Þ. Hence the integration over all points along the line is efficiently approximated, requiring only three to five finite element analyses per line. 3. Numerical example 3.1. General remarks

Fig. 5. Line sampling in the uncertain structural paramameter space.

ðiÞ

It should be noted that pf ðcÞ > 0 for 1 < c < 1, since the position along each line specifies uniquely the properties of a strucðiÞ ture subjected to random excitation. Since pf ðcÞ is always a ðiÞ positive quantity, it is proposed to represent pf ðcÞ as the exponential of a linear or quadratic polynomial, ðiÞ

pf ðcÞ ¼ exp½a0 þ a1 c þ a2 c2 =2;

ð70Þ

where the coefficients are obtained by solving a least square problem

h i ðiÞ a0 þ a1 cl þ a2 c2l =2 ¼ ln pf ðcl Þ ;

l ¼ 1; 2; . . . ; 5

ð71Þ

The advantage of this approximation is that the infinite integral in Eq. (69) can be represented in closed form for a2 < 1



1 a21 ðiÞ : pf ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp a0 þ 2ð1  a2 Þ 1  a2

ð72Þ

h i ðiÞ A linear approximation of ln pf ðcÞ in the least square solution should be applied ða2 ¼ 0Þ, whenever a2 > 0:4, which might be

The method, as developed in Section 2, is now exemplified within this section. Usually, proposed methods have been demonstrated by applying them to exceptionally small academic type of examples, where the feasibility of the approach to real world problems remains an open issue, because computational difficulties and many issues of computational efficiency are avoided. Since reliability evaluations in structural dynamics are computationally very demanding, the computational efficiency of the proposed method is naturally in the focus of interest. However, to be of practical value, the efficiency of procedures should be demonstrated by applying them to models which reflect the complexity of relevant engineering structures. By choosing a realistic FE-model, modeling a 12-story building made of reinforced concrete, the requirement of practical applicability is shown. 3.2. Structural system 3.2.1. Geometry A twelve story building with an additional cellar floor made of reinforced concrete is considered. The FE-model consists of 4046 nodes and 5972 elements using shell and 3-D beam elements for modeling the girders, resulting in 24,276 degrees of freedom. Fig. 6 provides a view of the building showing the displacement field according to the first three mode shapes. The building is axis-symmetric and consists of 13 floors of 24:0  24:0 m. The foundation plate is 0.5 m thick and rests on soil modeled by elastic springs. Fig. 7 shows a plan view for all floors, with the exception that the cellar floor is surrounded by cellar walls of 0.3 m thickness. The height of each story is 4.0 m. The weight of the structure is carried by four groups of concrete walls forming a cross. The four groups of supporting concrete walls are connected by four girders of type g-1 the height of which is 1.2 m and a width identical to the

Fig. 6. First three mode shapes of twelve story building.

81

H.J. Pradlwarter, G.I. Schuëller / Computers and Structures 88 (2010) 74–86

lation coefficients in the order of 0.86. The Young’s moduli of the remaining primary reinforced concrete parts are assumed somewhat higher to be 3.2e+10 N/m2 with a coefficient of variation of 0.144. This relatively large value reflects, to a considerable extent, the lack of knowledge on the dynamic stiffness. It implies also a large correlation of the stiffness, which is reflected by a correlation coefficient in the order of 0.7. The live loads (masses) are represented by a function of five random variables per floor, which are assumed to be independent between the floors and correlated within the floors. The mean value is assumed to be 100 kg/m2 and the coefficient of variation of 0.224 with correlation coefficients of 0.8. 3.3. Dynamic excitation

Fig. 7. Floor plan of twelve story building supported by croncrete walls (w) and girders of type g-1 and type g-2. Critical strains are considered at positions p-1 and p-2.

walls. The thickness of the supporting walls varies with respect to the respective story. In the cellar floor and first and second floor the wall thickness is 0.6 m, the next two floors (third and fourth) 0.5 m, and 0.4 for the fifth until the twelfth floor. The concrete plate of 0.2 m thickness is supported also by girders of type g-2 of height 0.5 m and constant width of 0.3 m. Above the concrete plate, a floor construction with the weight of 150 kg/m2 have been assumed for all floors. 3.2.2. Uncertain structural properties Concrete and reinforced concrete are widely used in building constructions. Its material properties, however, are hard to specify by just a few characteristic quantities. The stiffness of reinforced concrete construction parts depends on many factors which are still difficult to predict and to control. Due to inhomogeneities and possible cracking, the stiffness under dynamic loading reveals uncertainties. Moreover, it might depend also on the loading history (see e.g. [3]). At the present state-of-the-art, quantitative models for reinforced concrete structures, which are both feasible and realistic (including bond-slip, cracking, etc.) are not used yet for ordinary structural analysis as they are still topics of active research. Uncertainties of the stiffness is partially the result of inherent random factors and to a large part by lack of knowledge of the actual mechanisms. A practical way to cover all these uncertainties is to take a single quantity, such as the Young’s modulus, to represent the uncertainties in the stiffness, with the understanding that this quantity should cover the uncertainty of the stiffness and not only of the physical elastic constant. All structural uncertainties are assumed to be represented as a function of 244 independent standard normal variables. The stiffness of the soil bedding is represented by 9 random variables and the modal damping ratio by a single random variable. The mean value of the stiffness of the soil is assumed to be 4.44e+08 N/m3 in the horizontal and 8.88e+08 N/m3 in the vertical direction, with a coefficient of variation of 0.224. 169 random variables do represent the Young’s modulus of various parts of the reinforced concrete structure. The mean value for the Young’s modulus for the foundation plate and cellar walls are assumed to be 2.8e+10 N/m2 and its variability is modeled by a coefficient of variation of 0.215 and corre-

In this numerical example, the dynamic excitation is due to earthquake ground motion in terms of ground accelerations. Future ground motions are highly uncertain, regarding its amplitudes, frequency content and durations, respectively. For physical reasons, the average velocity and acceleration must vanish, and the frequency content depends on the unknown source mechanism and is influenced substantially by local soil conditions. The intensity (acceleration amplitudes) depends on the distance to the earthquake source and the duration on the magnitude of the energy release. It is difficult to cover all these uncertainties by a credible excitation model. In engineering practice, the excitation model is selected such that the model covers uncertainties conditioned on the acceptable hazard. There are several options to derive a suitable realistic excitation. In the ideal case, many records from different events are available – which after a suitable scaling – are used further to establish the covariance matrix (see (8)) from which the Karhunen–Loève terms can be deduced in a straightforward manner. A further option is to establish the covariance matrix such that it is compatible to some response spectra requirements. In this work, an approach based on filtered white noise is used. To cover the unknown frequency content and random amplitudes of the acceleration, the model proposed in [4] is applied, where the acceleration is represented as the output of the response of a linear filter excited by white noise

aðtÞ ¼ X2g v 1 ðtÞ þ 2fg Xg v 2 ðtÞ  X2f v 3 ðtÞ  2ff Xf v 4 ðtÞ;

ð73Þ

in which Xg represents the dominant frequency of the ground and Xf ensure that the spectrum of the acceleration tends to zero for frequencies approaching zero. The linear filter is described by the differential equation [4]

8 9 v1 > > > > > > = <

2

0 6 X2 d v2 g 6 ¼6 > 4 0 dt > v 3> > > ; : > 2

v4

Xg

1

0

2fg Xg

0

0

0

2fg Xg

X2f

9 3 8 9 8 0 > v1 > > > > > > > > > > > < = < 7 0 7 wðtÞ = v2 þ 7 1 5 > v3 > > > 0 > > > > > > ; : ; > : > 2ff Xf 0 v4 0

ð74Þ where wðtÞ denotes white noise. The intensity of wðtÞ can either be characterized by its intensity IðtÞ or by the constant (over frequency x, i.e. Sðx; tÞ ¼ S0 ðtÞ) spectral density S0 ðtÞ ¼ IðtÞ=ð2pÞ, which satisfy the following relation [8]

E½wðtÞwðt þ sÞ ¼ IðtÞdðsÞ ¼ 2pS0 ðtÞdðsÞ

ð75Þ

where dðÞ denotes Dirac’s delta function with the property dðsÞ ¼ 0 R for s – 0 and  dðsÞds ¼ 1;  > 0. Since white noise is uncorrelated for any s > 0, white noise excitation can be discretized by a sequence of independent impulses with zero mean and standard deviation rI ðkDtÞ or with linearly interpolated independent force amplitudes with standard deviation rF ðkDtÞ at subsequent time steps kDt

82

H.J. Pradlwarter, G.I. Schuëller / Computers and Structures 88 (2010) 74–86

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rI ðkDtÞ ¼ IðkDtÞDt rffiffiffiffiffiffiffiffiffiffiffiffiffi IðkDtÞ : rF ðkDtÞ ¼ Dt

Variance of horizontal acceleration

ð76Þ ð77Þ 1

The non-stationary properties of earthquake ground motion is described by a time varying white noise intensity. 0.8

for t 6 0 ð78Þ

2 4

for 0 < t < 2 for 2 6 t 6 8 for t P 8:

[m /s ]

8 0 > > > < t=2 IðtÞ ¼ I0 >1 > > : expð0:16ðt  8ÞÞ

0.4

The following constants are used in this work: I0 ¼ 0:16 m2 =s3 ; Xg ¼ 2p rad=s; fg ¼ ff ¼ 0:6, and xf ¼ 0:1Xg . In a next step, the covariance matrix C ¼ ½C ij  of the acceleration is established with Dt ¼ 0:01 s, and 0 < i; j 6 1600 and C ij ¼ E½aðiDtÞaðjDtÞ]. For this purpose, the impulse response function of the linear filter is determined by the free motion of the filter with zero initial displacement, and v_ 1 ð0Þ ¼ v_ 3 ð0Þ ¼ v_ 4 ð0Þ ¼ 0 and v_ 2 ð0Þ ¼ 1. To obtain the impulse response function of the acceleration aIRF ðtÞ, Eq. (73) is used. The result is shown in Fig. 8. To establish the covariance matrix, a single impulse at time kDt is considered first. The associated covariance is

C ij ðkÞ ¼ Dt  IðkDtÞ



aIRF ðði  kÞDtÞaIRF ððj  kÞDtÞ

for i P k and j P k

0

for i < k or j < k ð79Þ

0.2

0 0

C ij ðkÞ:

4

6

8

10

12

14

16

Fig. 9. Variance of horizontal acceleration over time.

Various KHL−terms for horiz. acceleration 0.4 10−th 30−th 60−th 80−th

0.2

ð80Þ 2

k¼0

[m/s ]

C ij ¼

2

time in [s]

Since the impulses at instants kDt are all independent, the coefficients C ij are given by the sum 1 600 X

0.6

The diagonal terms C ii ¼ E½a2 ðtÞ specify the variance over the time t i ¼ iDt and are shown in Fig. 9. Applying further the procedure as described in Section 2.1, the Karhunen–Loève functions aðjÞ ðtÞ (see Eq. (12)) are computed. For illustration, some of these functions are shown in Fig. 10. Finally, it should be stated that also the vertical earthquake acceleration has been modeled. It turned out that for the structure under consideration this vertical acceleration has only a negligible effect on the critical response and hence it is not described further.

0

−0.2

−0.4 0

2

4

6

8 10 time in [s]

12

14

16

Fig. 10. Karhunen–Loève functions aðjÞ ðtÞ; j ¼ 10; 30; 60; 80 of the horizontal acceleration.

Impulse acceleration response of horizontal filter 8

3.4. Critical response 7 6 5

2

[m/s ]

4 3 2 1 0 −1 −2

0

1

2

3

4 time in [s]

5

6

Fig. 8. Acceleration impulse response function of filter.

7

8

Under dynamic earthquake loading, two reinforced structural components are likely to exceed the critical strain. The strain in vertical direction of the walls in the basement, forming a cross, might be critical. These walls support the weight of the entire structure and are, in addition, highly strained by bending due to the horizontal earthquake accelerations. Since the structure is symmetric, and the horizontal accelerations are in both directions independent and identically distributed, it suffices to consider a single position, indicated by p-1 in Fig. 7. Girders form the remaining critical parts, denoted as g-1 in Fig. 7, which connect separated walls. These girders transmit large shear forces from one group of walls to the counterpart, and contribute essentially to the bending stiffness of the twelve story building. The critical part is the curvature of the girder at the connection to the walls. The position is indicated by p-2 in Fig. 7. Table 1 provides an overview of the strains in these components, i.e. for all floors. Accepting as limit state function a vertical

83

H.J. Pradlwarter, G.I. Schuëller / Computers and Structures 88 (2010) 74–86 Table 1 Critical strain due to static weights and maximum standard deviation due to horizontal earthquake motion. Floor

Width

Wall at pos. p-1

Girder g-1 at p-2

(–)

(m)

v. Strain static 106

v. Strain std. dev. 106

Curvature static 106

Curvature std. dev. 106

0 1 2 3 4 5 6 7 8 9 10 11 12

0.6 0.6 0.6 0.5 0.5 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4

143 110 100 108 97 107 93 80 66 53 39 25 18

266 228 181 170 134 125 80 73 60 46 33 20 12

9 27 27 32 32 41 41 42 42 42 42 42 34

281 466 548 643 654 704 675 620 544 457 367 289 187

strain of 0.0016 for the walls without having to expect serious damage and a maximum acceptable curvature of 0.004 in the girders g-1, the curvature in the girder are more likely to be exceeded. Hence, the reliability analysis will focus on the reliability of the girders for not exceeding the limit state function the curvature of 0.004. Taking into account the static solution of floor five,the limit state function has the lower and upper bounds are b ¼ 0:004þ  ¼ 0:004 þ 000041 ¼ 0:004041, respec000041 ¼ 0:003959 and b tively. 3.5. Reliability of critical component 3.5.1. Estimation by direct monte carlo simulation Before evaluating the first excursion probability, the Karhunen– Loève terms (see Eqs. (29) and (35)) of the critical response, i.e. the curvature at girder g-1, are determined. These terms yðjÞ ðt; hÞ depend only on the set of uncertain structural parameters h and their evaluation requires a single FE analysis for each distinct set hðjÞ . Hence, the investigation of the variability of the response due to uncertain structural properties is the computationally expensive part. The efficiency of the procedure will therefore be governed by the required number of FE analyses; the larger the FE-model, the more important this number will be. Fig. 11 shows some of these critical response functions for the nominal solution h ¼ 0. The associated first excursion problem is determined by Line sampling as outlined in Section 2.5, where

N ¼ 4000 lines have been used. Note that the number N is selected such that the computation time is approximately one fifth of the time for the FE analysis to compute the Karhunen–Loève terms. The failure probability assumes for the nominal system ðh ¼ 0Þ ^f ¼ 4:67  107 and a standard deviation of this estithe value p mate of rp^f ¼ 1:46  107 . However, it is well known that the reliability is quite sensitive with respect to the variations of structural properties. This sensitivity is shown by the results of Direct Monte Carlo simulation (see Eq. (55)) by randomly sampling over the uncertain structural parameters fhðjÞ g1600 j¼1 . Direct Monte Carlo sampling leads ^f ¼ 2:32  104 and the standard deviation of to the estimate p rpf ¼ 3:43  105 . Fig. 12 shows the histogram of all independent estimates. It is observed, that the estimate for the failure probability varies over many orders of magnitude. Hence, the procedure of Direct Monte Carlo Simulation cannot be used in an efficient manner. 3.5.2. Critical domain of uncertain structural parameters The domain of the structural parameters which contributes most to the total failure probability has been derived in Section 2.6. Eq. (61) provides a quantitative guidance to determine this domain. Its solution, shown in Table 2, requires the gradient of the standard deviation due to the dynamic excitation with respect to the uncertain structural parameters. To avoid direct differentiation for all 244 uncertain structural parameters, the gradient estimation

KHL-terms of curvature of girder at p-2 0.0003

Histogram of 1600 DMC estimates p_f=10^{-k} 1

1 10 40 80

0.0002

0.1

%

0.0001 0.01

0

-0.0001

0.001

-0.0002 0

2

4

6

8 time t [s]

10

12

14

16

0.0001 0

2

4

6

8

10

12

14

k Fig. 11. Karhunen–Loève terms yðjÞ ðtÞ; j1; 10; 40; 80 of the curvature of girder g-1 at position p-2 in floor five for h ¼ 0.

Fig. 12. 1600 estimates by direct Monte Carlo of the stochastic structure.

16

84

H.J. Pradlwarter, G.I. Schuëller / Computers and Structures 88 (2010) 74–86

Table 2 The first three most important components (k) of the structural design point h computed by three iterations s ¼ 1; 2; 3. s

cðsÞ

0 0.00070

k 26 244 114

hk 0 0 0

1 0.00190

ð0Þ

2 0.00110

ð1Þ

3 0.00111

ð2Þ

hk 1.92 1.52 0.69

ð3Þ

hk 2.07 1.18 0.76

hk 2.05 1.17 0.75

Twenty inpendent estimates of pf by LS 0.001 0.0009 0.0008 0.0007 0.0006 0.0005

procedure is employed [12,11] to compute these components which significantly influence the standard deviation of the critical response. It was observed that the variability of cðt ; hÞ is governed by only three random variables, namely the numbers 26, 244 and 114. Number 26 controls the stiffness of the girder of all floors, number 244 controls the log-normally distributed damping ratio, and the number 114 the local stiffness deviation at the considered girder at floor 5. These three random variables cause 97% of the total variability of cðt ; hÞ, i.e. Eq. (C.4) lead to  ¼ 0:03. The gradient estimation procedure has been used only in the first step, i.e. to identify the most important parameters and to compute the associated accuracy in terms of . This procedure required 44 FE analyses to arrive at the solution as shown. The iteration for s ¼ 1; 2 in Eq. (62) are then performed by finite differences of the three identified random variables. Table 2 summarizes the result of the failure point nearest to the origin in standard normal parameter space. 3.5.3. Line sampling within the uncertain structural parameter space The procedure outlined in Section 2.7 is applied to compute efficiently the total first excursion probability of the uncertain strucðjÞ tural system. Fig. 13 shows the five discrete values pf ðck Þ of the conditional failure probability at five points along 20 random lines by using bS ¼ 2:48 and D ¼ 0:6. The significant increase of the conditional failure probabilities in the direction of the design point demonstrates the high sensitivity of the failure probability with respect to the particular uncertain structural parameters. Fig. 14 shows twenty independent estimates of the total first excursion probability. To obtain these results 100 FE analyses (20 ^f ¼ 0:00046  5) have been carried out. The estimated mean is p with a standard deviation of rp^f ¼ 0:000040. When compared with results of Direct Monte Carlo sampling in Fig. 12 in the uncertain parameter space, one observes that these results are approximately by a factor of two larger. LS, however, explores the impor-

1

0.1

0.0004 0.0003 0.0002 0.0001 2

4

6

8

10

12

14

16

18

20

Fig. 14. Twenty independent estimates of the total first excursion probability of the ^f ¼ 0:00046 with standard deviation uncertain structure lead to the mean p rp^f ¼ 0:000040.

tant domain systematically and is therefore more reliable. This can be also seen by the substantial smaller coefficient of variation of approximately 8% compared with 16% and a 16 times larger number of FE analyses when using Direct MCS. 4. Conclusions The developments as shown here allow to draw the following conclusions: 1. The reliability evaluation of linear structural systems (general FE-models) with uncertain structural parameters subjected to general Gaussian dynamic excitation is feasible. To the authors knowledge, the feasibility of a reliability evaluation in dynamics for a large FE-model and a large number of structural uncertainties is demonstrated the first time in a numerical example. 2. A novel procedure to identify the critical uncertain parameters has been introduced. 3. Efficiency is gained by performing modal analysis, impulse response functions combined with Fast Fourier Transforms (FFT) and a new Line Sampling approach which focuses on the important failure domain within the uncertain parameter space. 4. The presented approach is accurate and robust, also for cases where the uncertainties of the structural parameters are large. 5. The reliability of linear systems is quite sensitive to structural uncertainties – as demonstrated in the numerical example – and therefore must not be ignored.

Acknowledgement

0.01

This work has been supported by the Austrian Science Foundation (FWF) under contract number P19781-N13 (Simulation Strategies for FE Systems under Uncertainties), which is gratefully acknowledged by the authors.

0.001

0.0001

Appendix A. Line sampling in a single direction 1e-05 2

4

6

8

10

12

14

16

18

20

Fig. 13. Five discrete conditional failure probabilities along twenty lines in the direction of the structural design point using bS ¼ 2:48 and D ¼ 0:6.

Line sampling (LS) [7,15] has been developed for estimating low probabilities of failure in high dimensional reliability problems. LS is particularly efficient if a so called important direction can be computed, which points towards the failure domain nearest to

85

H.J. Pradlwarter, G.I. Schuëller / Computers and Structures 88 (2010) 74–86

the origin (see Fig. 3). LS is robust, since unbiased reliability estimates are obtained, irrespectively of the important direction a, and any deviation from the optimal direction merely increases the variance of the estimate. Suppose the important direction a has been estimated by an appropriate procedure, e.g. gradient estimation, stochastic search algorithm, design point, etc. Let’s assume further that the reliability is computed in standard normal space of dimension n. Then, any point n in the n- dimensional standard normal space is represented by the projection on a

c ¼ aT n

ðA:1Þ

c2½i ðt; hÞ ¼

nþ1i X h

ðjÞ

y½i ðt; hÞ

i2

ðB:6Þ

j¼1

The important direction a½i for LS is selected at time ti ; a½i ¼ aðt i Þ, when the excursion probability assumes its highest value, or equivalently when the reliability index bi ðtÞ assumes it smallest value:

a½i ðhÞ ¼ aðti ; hÞ ðB:7Þ " # " # ð0Þ  ð0Þ  ð0Þ ð0Þ   b  y ðt i Þ y ðt i Þ  b b  y ðt; hÞ y ðt; hÞ  b 6 min ; min ; ; c½i ðti Þ c½i ðti Þ c½i ðt; hÞ c½i ðt; hÞ 06t6T

ðB:8Þ

and the ðn  1Þdimensional space n? ðaÞ:

n ¼ ca þ n? ðaÞ

ðA:2Þ f~ n?ðjÞ ðaÞgNj¼1

LS is performed by selecting random points in the space n? ðaÞ by applying direct Monte Carlo Simulation. Random lines are specified by the variable c and ~ n?ðjÞ ðaÞ

nðjÞ ðcÞ ¼ ca þ ~n?ðjÞ ðaÞ:

ðA:3Þ

j is determined, For each random line, the safe domain bj < c < b leading to an estimate for the failure probability ðjÞ j Þ; pf ¼ Uðbj Þ þ Uðb

ðA:4Þ

where UðÞ denotes the cumulative standard normal distribution. Finally, an unbiased estimate for the failure probability is obtained from the average

pf ¼

N 1 X ðjÞ p : N j¼1 f

ðA:5Þ

Appendix B. Evaluation of important directions in linear dynamics To obtain these specific important directions, aðtÞ at discrete instances ft i gLi¼1 are selected. ½i

a ðhÞ ¼

aðti ; hÞ

ðB:1Þ

The instance ti is determined in the ðn þ 1  iÞ-dimensional subspace n½i  n, defined as the subspace orthogonal to all previously selected important directions fa½j gi1 j¼1 and defining the initial condition n½1 ¼ n ½j

½i

ha ; n i ¼ 0 8 j ¼ 1; . . . ; i  1;

i ¼ 2; 3; . . . ; L

ðB:3Þ

½i

where R is a ðn þ 1  iÞ  ðn  iÞ dimensional orthonormal matrix orthogonal to a½i . This matrix is established by populating first all entries of R½i by independent random numbers and applying subsequently the well known Gram-Schmidt orthogonalization algorithm. The deterministic KL functions yðjÞ ðt; hÞ need also to be transformed because of the change of orientation. Define

Y½i ¼

ð1Þ ð2Þ ðnþ1iÞ ½y½i ðt; hÞ; y½i ðt; hÞ; . . . ; y½i ðt; hÞ

with the initial condition transformed by

ðjÞ y½1 ðt; hÞ

Y½i ðt; hÞ ¼ Y½i1 ðt; hÞ  R½i ;

i ¼ 2; 3; . . . ; L

2 ½i ðt; hÞ

In the following, the actual evaluation of the derivatives ðsÞ @ cðt  ; hðsÞ Þ=@hk is discussed, since it might considerably affect the efficiency of the procedure. For the case where only very few uncertain structural quantities are specified, differentiation for all parameters fhk gKk¼1 by finite differences can be carried out in a straight forward manner, requiring additional K finite element analyses (FEA). In general, however, the number K of uncertain structural quantities might be quite large, where K FEA would seriously reduce the efficiency of the procedure. Although, the number of uncertain structural properties might be large, only a limited number of uncertain parameters have usually a significant effect on the standard deviation cðt  ; hÞ. These quantities are sometimes known a priori and often they cannot be specified with sufficient confidence. Suppose an estimate where only I < K uncertain quantities significantly influence the standard deviation cðt  ; hðsÞ Þ. It is ðsÞ then proposed to compute first these derivatives @ cðt ; hðsÞ Þ=@hi2I by finite difference. These directly evaluated derivatives are denoted as

ji2fIg ¼ @ cðt ; hðsÞ Þ=@hðsÞ i2fIg

ðC:1Þ

Since the set fIg is just an estimate, a check whether indeed no important quantities hjRI are missing is necessary. For this purpose, the standard deviation cðt  ; hÞ has to be computed for the set fcðt  ; hðjÞ þ DðkÞ gIk¼0 , with the following properties:

Dð0Þ ¼ 0; i Dðj>0Þ ¼d i



i ¼ 1; 2; . . . ; K

ðC:2Þ

1 for 0 6 ui;j < 0:5

ðC:3Þ

þ1 for 0:5 6 ui;j 6 1

ðB:2Þ

where h; i denotes the inner product. The ðn þ 1  iÞ-dimensional subspace n½i is computed by the linear transformation

n½i ¼ n½i1  R½i ;

Appendix C. Gradient estimation

ðB:4Þ

ðjÞ

¼ y ðt; hÞ, these terms are also

ðB:5Þ

The variance c of the critical response yðt; hÞ in these subspace n½i is determined analogous as in Eq. (39)

with 0:01 6 d 6 0:1 and uniformly distributed independent random numbers 0 < ui;j < 1 drawn by a random number generator. The accuracy achieved by this procedure can be estimated by comparing the Euclidean norm of the vector b and the vector c,

¼

kck kbk

ðC:4Þ

where the components of the vectors are defined as follow:

bk ¼ cðt ; hðsÞ þ DðkÞ Þ  cðt  ; hðsÞ Þ; X c k ¼ bk  ji DðkÞ i

k ¼ 1; 2; . . . ; J

ðC:5Þ ðC:6Þ

i2fIg

In case  < 0:2, the estimate might be regarded as sufficiently accurate, otherwise the estimate requires improvements. For an efficient procedure for improvements it is referred to [12,11]. An abbreviated version is shown in the following: Components hj , which are likely to have a great effect, can be ðkÞ recognized by the correlation of Dj and the component ck . It is suggested to determine for all components fj R fIgg the estimate

86

H.J. Pradlwarter, G.I. Schuëller / Computers and Structures 88 (2010) 74–86

PJ

j^ j ¼ P

ðkÞ k¼1 Dj c k

 2 ¼ ðkÞ J k¼1 Dj

PJ

ðkÞ k¼1 Dj c k 2

Jd

ðC:7Þ

In a further step, the component m with the absolutely largest esti^ jRfIg j is selected and the derivative jm is computed by ^ m j P jj mate jj finite differences. The improvement is measured by  after updating the vector c

ck #ck  jm DðkÞ m :

ðC:8Þ

This procedure is repeated until a satisfactory accuracy is obtained. If the improvements are insignificant, the sample size J needs to be enlarged as shown in [11]. References [1] Adomian G. Stochastic Systems. New York: Academic Press; 1983. [2] Au SK, Beck JL. First excursion probabilities for linear systems by very efficient importance sampling. Probabilistic Engineering Mechanics 2001;16:193–207. [3] Chryssanthopoulos MK, Dymiotis C, Kappos AJ. Probabilistic evaluation of behaviour factors in ec8-designed r/c frames. Engineering Structures 2000;22: 10281041. [4] Clough RW, Penzien J. Dynamics of Structures. Intern. Student ed. Auckland: McGraw-Hill; 1975.

[5] Ghanem R, Spanos P. Stochastic Finite Elements: A Spectral Approach. Springer-Verlag; 1991. [6] Karhunen K. Über lineare Methoden in der Wahrscheinlichkeitsrechnung. American Academy of Science Fennicade Series A 1947;37:3–79. [7] Koutsourelakis PS, Pradlwarter HJ, Schuëller GI. Reliability of structures in high dimensions, part I: algorithms and applications. Probabilistic Engineering Mechanics 2004;19(4):409–17. [8] Lin YK. Probabilistic Theory of Structural Dynamics. New York: McGraw-Hill, Inc., McGraw-Hill Company; 1967. [9] M. Loève, Fonctions aleatoires du second ordre, supplemeent to P. Levy. In Processus Stochastic et Mouvement Brownien. Gauthier Villars, Paris, 1948. [10] W.L. Oberkampf, S.M. DeLand, B.M. Rutherford, K.V. Diegert, K.F. Alvin, Estimation of total uncertainty in computational simulation. Technical Report, Sandia National Laboratories, Albuquerque, NM SAND2000-0824, USA, 2000. [11] Pellissetti MF, Pradlwarter HJ, Schuëller GI. Relative importance of uncertain structural parameters, part II: applications. Computational Mechanics 2007;40(4):637–49. [12] Pradlwarter HJ. Relative importance of uncertain structural parameters part I: algorithm. Computational Mechanics 2007;40(4):627–35. [13] Pradlwarter HJ, Schuëller GI. Excursion probability of non-linear systems. International Journal of Non-Linear Mechanics 2004;39(9):1447–52. [14] Schenk CA, Schuëller GI. Uncertainty Assessment of Large Finite Elements Systems. Berlin/Heidelberg/New York: Springer-Verlag; 2005. [15] Schuëller GI, Pradlwarter HJ, Koutsourelakis PS. A critical appraisal of reliability estimation procedures for high dimensions. Probabilistic Engineering Mechanics 2004;19(4):463–74. [16] Soong TT, Grigoriu M. Random Vibration of Mechanical and Structural Systems. New Jersey: Prentice-Hall; 1993.