FERUM
FERUM
FERUM
UCB/SEMM-2000/08
STRUCTURAL ENGINEERING,
MECHANICS AND MATERIALS
by
Bruno SUDRET
and
Armen DER KIUREGHIAN
November 2000
by
Bruno Sudret and Armen Der Kiureghian
A report on research supported by
Electricit de France
under Award Number D56395-T6L29-RNE861
Report No. UCB/SEMM-2000/08
Structural Engineering, Mechanics and Materials
Department of Civil & Environmental Engineering
University of California, Berkeley
November 2000
Acknowledgements
This research was carried out during the post-doctoral stay of the rst author at the
Department of Civil & Environmental Engineering, University of California, Berkeley.
Contents
Part I : Review of the literature
1 Introduction
. . . . . . . . . . .
Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Generalities
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1
1.2
10
2.1
10
2.2
. . . . . . . . . . . . . . . . . .
10
2.3
. . . . . . . . . . . . . . . . . . . .
11
2.4
. . . . . . . . . . .
11
. . . . . . . . . . . . . . . . . . . . . . .
13
3.1
. . . . . . . . . . . . . . . . . . . . . . . . .
13
3.2
14
15
. . . . . . . . . . . . . . . . . . . . . . . . . .
17
5.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
17
5.2
19
vii
Contents
viii
5.3
5.4
5.2.1
Denition . . . . . . . . . . . . . . . . . . . . . . . . . .
19
5.2.2
Properties . . . . . . . . . . . . . . . . . . . . . . . . . .
20
5.2.3
. . . . . .
21
5.2.4
Conclusion
. . . . . . . . . . . . . . . . . . . . . . . . .
22
. . . . . . . . . . . . . . . . . . . . .
23
. . . . . . . . . . . . . . . . . . . . . . . .
23
5.3.1
Introduction
5.3.2
. . . .
24
25
5.4.1
25
5.4.2
Variance error . . . . . . . . . . . . . . . . . . . . . . . .
26
6.2
. . . . . . . . . . . . . . . . . . .
26
Early results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
26
6.1.1
EOLE
6.1.2
OSE
vs. KL
. . . . . . . . . . . . . . . . . . . . . . . .
26
. . . . . . . . . . . . . . . . . . . . . . . . .
27
28
6.2.1
28
6.2.2
28
6.2.3
28
6.2.4
. . . . . . . .
29
6.2.5
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . .
31
vs. KL
Conclusions
. . . . . . . . . . . . . . . . . . . . . . . . .
32
. . . . . . . . . . . . . . . . . . . . .
32
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
34
37
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
37
38
Contents
3
ix
. . . . . . . . . . . . . . . . . .
39
3.1
40
3.2
. . . . . . . . . . . . . . . . . . . .
40
40
4.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
40
4.2
41
4.3
. . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
41
. . . . . . . . . . . . . . . . . . . . . . . . . . .
42
5.1
43
5.2
. . . . . . . .
43
44
47
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
47
47
2.1
47
2.2
. . . . . . . . . . . . . . . . . . . . . . . . . .
48
2.3
49
2.4
Probabilistic transformation . . . . . . . . . . . . . . . . . . . . .
50
2.5
FORM, SORM
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
52
2.6
54
2.6.1
Early approaches . . . . . . . . . . . . . . . . . . . . . .
54
2.6.2
. . . . . . . . .
56
2.6.3
Conclusion
. . . . . . . . . . . . . . . . . . . . . . . . .
57
. . . . . . . . . . . . . . . . . . . .
57
3.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
57
3.2
58
. . . . . . . . . .
Contents
3.2.1
. . . . . . . . . . . . .
59
3.2.2
60
3.2.3
60
3.2.4
3.2.5
Examples
. . . . . . . . . . . . .
. . . . .
61
. . . . . . . . . . . . . . . . . . . . . . . . . .
62
3.3
62
3.4
63
3.5
. . . . . . . . . . . . . . . . .
64
Sensitivity analysis
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
65
67
5.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
67
5.2
67
5.3
. . . . . . . . . . . . . . . . . . . .
68
5.4
69
5.5
71
5.6
71
5.7
Conclusions
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
72
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
73
Conclusions
. . . . . . . . . . . . . . . . . . . . . . .
75
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
75
77
2.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
77
2.2
. . . . . . . . . . .
77
2.3
. . . . . . . . . . . . . . . . . . .
77
2.4
2.5
. . . . . . .
L2 ( ; F ; P )
. . . . . .
78
79
Contents
2.6
3
xi
. . . . . . . . . . . . . . . . . . . .
81
Computational aspects . . . . . . . . . . . . . . . . . . . . . . . . . . . .
83
3.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
83
3.2
83
3.3
Solution algorithms . . . . . . . . . . . . . . . . . . . . . . . . . .
83
3.4
Hierarchical approach . . . . . . . . . . . . . . . . . . . . . . . . .
84
Extensions of SSFEM . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
85
4.1
85
4.1.1
. . . . . . . . . . . . . . . .
86
4.1.2
. . . . . . . . . . . . . . . . . .
86
4.2
87
4.3
Hybrid SSFEM
88
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.1
. . . . . . . . . . . . . . . . . .
4.3.2
4.3.3
Concluding remarks
88
. . . . . . . . . . . . . . . .
88
. . . . . . . . . . . . . . . . . . . .
89
89
. . . . . . . . . . . . . . . . . . .
90
93
A.1.1 Denition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
93
94
96
6 Conclusions
99
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
99
Contents
xii
103
Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
. . . . . . . . . . . . . . . . 104
107
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
. . . . . . . . . . . . . . . . . . . . . . . . 107
2.1
2.2
Discretization procedure
. . . . . . . . . . . . . . . . . . . . . . . . . . . 109
3.1
3.2
3.2.2
Two-dimensional case
3.2.3
. . . . . . . . . . . . . . . . . . . 111
. . . . . . 111
3.3
3.4
3.5
3.4.1
General formulation
3.4.2
. . . . . . . . . . . . . . . . . . . . 112
. . . . . . . . . . . . . . . . . . . . . . . 115
Visualization tools
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
Contents
xiii
3 Implementation of SSFEM
1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
1.1
Preliminaries
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
1.2
SSFEM pre-processing
. . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
2.1
2.2
Polynomial chaos
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
3.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
3.2
3.3
3.4
. . . . . . . . . . . . . . 125
3.4.2
3.4.3
3.5
. . . . . . . . . . . . 123
3.4.1
variable
119
. . . . . . . . . . . . . . . . . . . . . . . . . . . 125
Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
SSFEM Analysis
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
4.1
. . . . . . . . . . . . . . . . . 127
4.2
4.3
4.4
. . . . . . . . . . . . . . . . . . . . . . . . . . 130
5.2
5.3
Reliability analysis
5.4
Conclusions
. . . . . . . . . . . . . . . . . . . . . . . 131
. . . . . . . . . . . . . . . . . . . . . . . . . . 131
. . . . . . . . 132
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
Contents
xiv
133
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
2.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
2.2
2.3
2.4
. . . . . . . . . . . . . . . . 135
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
3.2
3.3
3.4
. . . . . . . . . . . . . . 137
4.2
4.3
4.4
5
. . . . . . . . . . . . . . . . . . . 134
. . . . . . . . . . . . . 139
. . . . . . . . . . . . . . . . . . 139
. . . . . . . . . . . . . . . . . . . 140
4.2.1
4.2.2
. . 140
. . . . . . . . . . . . . . . . . . . 143
4.3.1
4.3.2
Conclusions
. . . . . . . . . . . . . . . . . . . . . 145
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
147
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
. . . . . . . 148
. . . . . . . . . . . 148
Contents
2.2
3
xv
. . . . . . 148
. . . . . . . . 150
3.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
3.2
. . . . . . . . . . . . . . . . . 150
3.2.1
3.2.2
3.2.3
. . . . . . . . . . . . . . . . . . . 151
3.3
3.4
3.5
. . . . . . . . 156
3.6
One-dimensional
. . . . . . . . 156
3.7
3.8
Introduction
. . . . . . . . . . . . . . . . . . . . . . . . 160
3.8.2
3.9
3.10
Conclusions
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162
. . . . . . . 164
4.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164
4.2
4.3
4.4
4.5
Conclusions
Conclusions
6 Conclusion
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169
171
Part I
Review of the literature
Chapter 1
Introduction
Modeling a mechanical system can be dened as the mathematical idealization of the
basic variables
(describing system geometry, loading, material properties), response variables (displace-
ment, strain, stresses) and the relationships between these various quantities.
For a long time, researchers have focused their attention on improving structural models
(beams, shells, continua, ...) and constitutive laws (elasticity, plasticity, damage theories, ...). With the development of computer science, a great amount of work has been
devoted to numerically evaluate approximated solutions of the boundary value problems
describing the mechanical system. The nite element method is probably nowadays the
most advanced approach for solution of these problems.
However the increasing accuracy of the constitutive models and the constant enhancement of the computational tools does not solve the problem of identication of the
model parameters and the uncertainties associated with their estimation. Moreover, in
most civil engineering applications, the intrinsic randomness of materials (soil, rock,
concrete, ...) or loads (wind, earthquake motion, ...) is such that deterministic models
using average characteristics at best lead to rough representations of the reality.
Accounting for randomness and spatial variability of the mechanical properties of materials is one of the tasks of stochastic or probabilistic mechanics, which has developed fast
in the last ten years. The aim of this report is to present a state-of-the-art review of the
existing methods in this eld. Having industrial applications of these methods in mind,
attention will be mainly focused on nite element approaches. The literature on this
topic has been classied by Matthies
et al.
presents the latest developments of various approaches. These three contributions will
Chapter 1. Introduction
be used as a basis for the present report, where only those parts related to the present
the theories aiming at calculating the rst two statistical moments of the response
the
They are based on the denition of a limit state function. As failure is usually
associated with rare events, the tails of the probability density functions (PDFs)
of response quantities are of interest in this matter.
the
It has been noted in the reviewed literature that these three categories of approaches are
investigated by dierent communities of researchers having few interactions with each
other. We will try to show in this report that the ingredients utilized in these methods
have many common features.
Note that the above classication is somewhat subjective. Indeed results obtained as
byproducts of the main analysis tend to break the walls between these classes, as shown
in the following examples, which will be investigated in detail later on :
1 The book by Haldar and Mahadevan (2000) should be mentioned for the sake of completeness. Due
to its recent publication, it could not be reviewed for the present report.
2. Outline
However, it is expected that methods pertaining to one of the above categories will
not be ecient in the computation of byproducts. As mentioned before, due to the
compartmentalization of the research groups, no signicant comparisons have been made
so far.
2 Outline
A common ingredient of all the methods mentioned above is the need to represent
the spatial variability of the input parameters. This is done by using a random eld
representation. For computational needs, these random elds have to be discretized in
an optimal way. Chapter 2 covers the methods for discretization of random elds.
Chapter 3 deals with second moment methods in the context of nite element analysis.
These methods give results in terms of
The ingredients for reliability approaches are rst introduced in a general context. Then
the specic modications to be introduced in the nite element context are presented.
This approach was rst proposed by Der Kiureghian and Taylor (1983).
Chapter 5 is devoted to the spectral stochastic nite element method (SSFEM). This
method was introduced by Ghanem and Spanos (1991 ). The main concepts will be
presented as well as a summary of the applications found in the literature.
As a conclusion, we will present a scheme for comparing the SSFEM and the reliability
approaches on the problem of evaluating small probabilities of occurrence as well as
the PDF of a given response quantity for a specic example. As noticed before, no
comparison of this type has been reported in the literature so far.
Chapter 2
Methods for discretization of random
elds
The engineering applications in the scope of this report require representation of uncertainties in the mechanical properties of continuous media. The mathematical theory for
this is random elds. For denitions and general properties, the reader is referred to Lin
(1967) and Vanmarcke (1983).
1 Generalities
The introduction of probabilistic approaches in mechanical problems requires advanced
mathematical tools. This section is devoted to the presentation of some of them. Should
the reader be already familiar with this material, the following will give him/her at least
the notation used throughout the report.
trial.
outcomes of a trial form the sample space of the phenomenon, denoted hereinafter by .
An event E is dened as a subset of containing outcomes 2 . Probability theory
aims at associating numbers to events, i.e. their probability of occurrence. Let P denote
this so-called probability measure. The collection of possible events having well-dened
probabilities is called the -algebra associated with , denoted here by F . Finally the
probability space constructed by means of this notions is denoted by ( ; F ; P ).
A real
random variable X
is a mapping
X : ( ; F ; P )
variables, the probability density function (PDF) and cumulative distribution function
fX (x)
and
FX (x),
being possibly
X , the
dependency on the outcomes may be added in some cases as in X ( ). A random vector
(2.1-b)
(2.1-c)
are :
E [X ] =
Z 1
x fX (x) dx
2 = E (X
E [X n ] =
Z 1
)2 =
E [].
Z 1
(x
Cov [X ; Y ] = E [(X
Cov [X ; Y ] =
n-th
)2 fX (x) dx
xn fX (x) dx
X )(Y
and
is :
Y )]
Z 1Z 1
(x X )(y
Y ) fX;Y (x ; y ) dx dy
(E [X 2 ] < 1) is
space.
(2.5)
< f ; g >L2 (
) =
f (x) g (x) d
1. Generalities
Galerkin procedure.
H (x) attached to point x is a random variable or a random vector. It is one- or multidimensional according to the dimension d of x, that is d = 1 or d > 1. For the sake
of simplicity, we consider in the following univariate multidimensional elds. In practical terms, this corresponds to the modeling of mechanical properties including Young's
modulus, Poisson's ratio, yield stress, etc., as statistically independent elds.
Gaussian
if any vector
power spectrum
1
SHH (! ) =
2
(2.6)
Z 1
~(x) e
i!x dx
nite set of random variables fi ; i = 1 ; ::: ng, grouped in a random vector denoted by
:
H (x) Discretization
! H^ (x) = F [x ; ]
(2.7)
The main topic here is to dene the best approximation with respect to some error estimator, that is the one using the minimal number of random variables. The discretization
methods can be divided into three groups :
point discretization, where the random variables fi g are selected values of H ()
at some given points
average discretization, where fi g are weighted integrals of H () over a domain
(2.8)
xi .
i =
H (x) w(x) d
10
Reviews of several discretization methods can be found in Li and Der Kiureghian (1993);
Ditlevsen (1996); Matthies
xc of this element :
H^ (x) = H (xc) ;
(2.9)
H^ () is
= fH (x1c ) ; ::: H (xNc e )g (Ne being
The
approximated
eld
then
x 2
e
entirely
dened
by
the
random
vector
and covariance matrix are obtained from the mean, variance and autocorrelation
^ ()
coecient functions of H () evaluated at the element centroids. Each realization of H
is piecewise constant, the discontinuities being localized at the element boundaries. It has
been shown (Der Kiureghian and Ke, 1988) that the MP method tends to
over-represent
(2.10)
q
X
i=1
Ni (x) H (xi )
x 2
e
q is the number of nodes of element e, xi the coordinates of the i-th node and Ni
^ () is
polynomial shape functions associated with the element. The approximated eld H
obtained in this case from = fH (x1 ) ; ::: H (xN )g, where fxi ; i = 1 ; ::: N g is the set
where
11
(2.11)
E H^ (x)
(2.12)
Each realization of
=
=
q
X
H^ () read :
x)(xi)
Ni (
i=1
q X
q
X
i=1 j =1
midpoint method.
(1995). Assuming that every integration appearing in the nite element resolution scheme
is obtained from integrand evaluation at
discretize the random eld by associating a single random variable to each of these Gauss
points. This gives accurate results for short correlation length. However the total number
of random variables involved increases dramatically with the size of the problem.
see Ditlevsen (1996). In the context of point discretization methods, the approximated
eld
H^ ()
follows :
H^ (x) = a(x) +
(2.13)
where
q
X
i=1
a(x)h and bi (x) areidetermined by minimizing in each point x the variance of the error
Var H (x) H^ (x) subject to H^ (x) being an unbiased estimator of H (x) in the mean.
These conditions write :
(2.14)
(2.15)
8 x 2
;
Minimize
with
E H (x)
H^ (x) = 0
Eq.(2.15) requires :
(2.16)
12
(2.17)
2
H (x) H^ (x)
q
X
i=1
q X
q
X
i=1 j =1
8 i = 1 ; ::: q
(2.19)
q
X
Cov [H (x) ; i ] +
j =1
bj (x)Cov [i ; j ] = 0
H(x) + b(x) = 0
(2.20)
where
(2.21)
(2.22)
H^ (x) = (x)
T
H (x)
1
q
X
i=1
i
1 H(x) i
from which it is seen that OLE is nothing but a shape function discretization method
where, setting the mean function aside, the shape functions read :
(2.23)
N OLE (x)
i
H(x) i =
1
q
X
j =1
(2.24)
TH(x) 1 H(x)
13
(2.25)
Thus requiring the error variance to be minimized is equivalent to requiring the error
and the approximated eld to be uncorrelated. Both statements can be interpreted as
constant
the element :
H (x) d
e
H^ (x) =
e
j
ej
(2.26)
Vector
is
then
dened
= fH e ; e = 1 ; ::: Ne g.
as
the
H e ; x 2
e
collection
of
these
random
variables,
are
that
is
computed from
(1983) gives results for homogeneous elds and two-dimensional rectangular domains.
the local variance of the random eld (Der Kiureghian and Ke, 1988).
Diculties involved in this method are reported by Matthies
et al. (1997) :
the approximation for non rectangular elements (which can be dealt with by a
collection of non overlapping rectangular ones) may lead to a non-positive denite
covariance matrix.
i
is almost impossible
to obtain except for Gaussian random elds. For the sake of exhaustivity, recent
work from Knabe
14
ab
and also investigated by Takada (1990 , ) in the context of stochastic nite elements.
It is claimed not to require any discretization of the random eld and thus seems to be
particularly attractive. In the context of
the
linear elasticity,
element stiness matrices as basic random quantities. More precisely, using standard
nite element notations, the stiness matrix associated with a given element occupying
a volume
reads :
k =
e
(2.27)
B T D B d
e
D denotes the elasticity matrix, and B is a matrix that relates the components
where
D(x ; ) = Do [1 + H (x ; )]
1
where D o is the mean value and H (x ; ) is a zero-mean process. Thus Eq.(2.27) can
(2.28)
be rewritten as :
(2.29)
e ( )
=k
; k
e + e ( )
o
e ( )
H (x ; ) B T D o B d
e
functions with respect to the coordinates. Hence they are polynomials in the latter, say
(2.30)
e ( )
ij
Pij
Pij
Pij (x; y; z ) H (x ; ) d
e
are obtained from those of
B and Do . Let us
as :
(2.31)
Pij (x; y; z ) =
NWI
X
l=1
alij xl y l z l
exponents
(2.32)
el ()
xl y l z
l H (x ; )d
e
15
it follows that :
NWI
X
(2.33)
k
e ( ) =
ij
alij
in a matrix
l=1
k =k
e
(2.34)
e+
o
NWI
X
l=1
kel el
keo and fkel; l = 1 ; ::: NWIg are deterministic matrices and el
mesh-dependent as
et al.
of discretization similar to the shape function approach (see section 2.2). Moreover, if the
correlation length of the random eld is small compared to the size of integration domain
e ,
the accuracy of the method is questionable. Indeed, the shape functions usually
e.g.
Young's modulus and cross-section) may not give good results when these properties
are rapidly varying in the element. The problem of accuracy of the weighted integrals
approach seems not have been addressed in detail in the literature. A comprehensive
study including the denition and computation of error estimators would help clarify
this issue.
Applications of the weighted integral method for evaluating response variability of the
system will be discussed later (Chapter 3, section 4).
sian random elds were considered, with three dierent correlation structures, namely
16
A (x ; x0 ) = exp(
(2.35)
k x x0 k )
a
k x x0 k2 )
B (x ; x0 ) = exp(
a2
2:2 k x x0 k
sin(
)
0
a
C (x ; x ) =
2:2 k x x0 k
a
(2.36)
(2.37)
where
and
l=a :
(2.38)
Applying OLE, four dierent sets of discretization points are used, namely the nodes of
the element under consideration, or the nodes of 3, 5 or 7 adjacent elements respectively.
For type A correlation, the error remains large even for a small element size
this case (because the autocorrelation function is not dierentiable at the origin,
see Vanmarcke (1983))
It is seen that OLE gives better results than SF in all cases. As mentioned before,
OLE is basically a SF approach, where the shape functions are not prescribed
polynomials, but the optimal functions to minimize the variance of the error.
to the initial
one is also given by Li and Der Kiureghian (1993). In all cases, OLE leads to better
accuracy in the discretization than MP, SA and SF.
17
Figure 2.1: Discretization errors for OLE method with varying grid and element size
(after Li and Der Kiureghian (1993))
(2.39)
w(x)
H (x ; ) w(x) d
(2.40)
weighted integrals
1
e (x) =
1
0
if
x 2
e
otherwise
18
Figure 2.2: Comparison of errors for MP, SA, SF and OLE for varying element size (after
Li and Der Kiureghian (1993))
as
a nite summation :
H^ (x ; ) =
(2.41)
N
X
i=1
H^ (x ; o )
L2(
)
point of view, the basis used so far are not optimal (for instance, in case of MP, SA and
SF, because the basis functions
The discretization methods presented in the present section aim at expanding any real-
19
Table 2.1: Weight functions and deterministic basis unifying MP, SF, SA, OLE methods
Method
weight function
(x xc)
1
e (x)
j
ej
MP
SA
(x
SF
w(x)
'i (x)
1
e (x)
1
e (x)
polynomial
xi)
functions
Ni (x)
shape
(x
OLE
xi)
Denition
position
The set of
H (x ; o ) is expanded
is
(2.42)
8 i = 1 ; :::
(2.43)
H (x ; ) = (x) +
1 p
X
i=1
i i () 'i(x)
fi(); i = 1 ; ::: g denotes the coordinates of the realization of the random eld
with respect to the set of deterministic functions f'i g. Taking now into account all
possible realizations of the eld, fi ; i = 1 ; ::: g becomes a numerable set of random
where
variables.
When calculating
20
to
(2.44)
This means that
fi; i
= 1 ; ::: g
(Kronecker symbol)
forms a set of
orthonormal
separation of the
Properties
(2.45)
M p
X
i=1
minimized (with respect to the value it would take when any other complete basis
fhi(x)g is chosen).
Due to the orthonormality of the eigenfunctions, it is easy to get a closed form for
each random variable appearing in the series as the following linear transform :
1
[H (x ; ) (x)] 'i (x) d
(2.46)
i () = p
i
Hence when H () is a Gaussian random eld, each random variable i is Gaussian.
It follows that fi g form in this case a set of independent standard normal variables.
Furthermore, it can be shown (Love, 1977) that the Karhunen-Love expansion
of Gaussian elds is almost surely convergent.
From Eq.(2.45), the error variance obtained when truncating the expansion after
(2.47)
M
X
i=1
21
The righthand side of the above equation is always positive because it is the variance of some quantity. This means that the Karhunen-Love expansion always
5.2.3
Eq.(2.42) can be solved analytically only for few autocovariance functions and geometries
of
. Detailed closed form solutions for triangular and exponential covariance functions
for one-dimensional homogeneous elds can be found in Spanos and Ghanem (1989),
= [ a ; a].
L2(
).
'k (x) =
(2.48)
where
dik
1
X
i=1
dik hi (x)
are the unknown coecients. The Galerkin procedure aims at obtaining the
'k (:) when truncating the above series after the N -th term. This
projecting 'k onto the space HN spanned by fhi (:) ; i = 1 ; ::: N g.
best approximation of
is accomplished by
(2.49)
N (x) =
N
X
i=1
dik
Z
CHH (x ; x0 ) hi (x0 ) d
x0
(2.50)
HN in
< N ; hj >
L2(
).
Z
k hi (x)
This writes :
N (x) hj (x) d
= 0 j = 1 ; ::: N
(2.51)
CD = BD
22
(2.52-a)
(2.52-b)
(2.52-c)
(2.52-d)
Bij
Cij
Dij
ij
=
=
hi (x) hj (x) d
Z
Z
= dij
= ij j
Kronecker symbol
and eigen-
i . This solution scheme can be implemented using the nite element mesh shape
functions as the basis f(hi ()g (see Ghanem and Spanos (1991b , chap. 5.3) for the exam-
values
ple of a curved plate). Other complete sets of deterministic functions can also be chosen,
as described in the next section.
5.2.4
Conclusion
Due to its useful properties, the Karhunen-Love expansion has been widely used in
stochastic nite element approaches. Details and further literature will be given in Chapters 5.
The main issue when using the Karhunen-Love expansion is to solve the eigenvalue
problem (2.42). In most applications found in the literature, the exponential autocovariance function is used in conjunction with square geometries to take advantage of the
closed form solution in this case. This poses a problem in industrial applications (where
complex geometries will be encountered), because :
the scheme presented in Section 5.2.3 for numerically solving(2.42) requires additional computations,
in a
square-shape volume and use the latter to solve in a closed form (when possible) the
eigenvalue problem. Surprisingly, this assertion, earlier made by Li and Der Kiureghian
(1993), did not receive attention in the literature.
23
Introduction
The Karhunen-Love expansion presented in the above section is an ecient representation of random elds. However, it requires solving an integral eigenvalue problem to
determine the complete set of orthogonal functions
ab initio a
complete set of orthogonal functions. A similar idea had been used earlier by Lawrence
(1987).
Let
Z
(2.53)
hi (x) hj (x) d
= ij
orthonormal, i.e.
( Kronecker symbol)
Let
outcomes of the eld, the coecients in the expansion become random variables. Thus
the following expansion holds :
(2.54)
H (x ; ) = (x) +
where
1
X
i=1
i () hi (x)
that :
(2.55-a)
i () =
(2.55-b)
If
H ()
[H (x ; )
Z
Z
variables, possibly
correlated. In
(x)] hi (x) d
fig1i=1
24
fhi(x)g1i=1
(Legendre polynomi-
als were used by Zhang and Ellingwood (1994)) and select the number of terms
retained for the approximation,
Compute
the
covariance
e.g. M .
matrix
of
the
zero-mean
Gaussian
vector
H^ (x ; ) = (x) +
(2.56)
5.3.2
M
X
i=1
i () hi (x)
The discretization of Gaussian random elds using OSE involves correlated Gaussian
random variables
standard
matrix
=
(2.57)
where
is related to by :
= 1=2
(2.58)
(2.59)
i () =
M
X
k=1
ik k k ()
Hence :
H^ (x ; ) = (x) +
(2.60)
M
X
M
X
ik
i=1 k=1
M p
X
)+
k k ()
k=1
= (x
k k () hi (x)
M
X
(2.61)
'k (x) =
M
X
i=1
ik hi (x)
i=1
ik hi (x)
25
M p
X
H^ (x ; ) = (x) +
(2.62)
k=1
By comparing the above developments with the numerical solution of the eigenvalue
problem associated with the autocovariance function (2.42) (see Section 5.2.3), the following important conclusion originally pointed out by Zhang and Ellingwood (1994) can
fhi(x)g1i=1 is strictly
equivalent to the Karhunen-Love expansion in the case when the eigenfunctions 'k (x)
fh (x)g1
i
CHH
i=1 .
Kiureghian (1993). It is an extension of OLE (see section 2.4) using a spectral representation of the vector of nodal variables
Assuming that
() = +
(2.63)
where
.
N p
X
i=1
i i () i
fi ; i = 1 ; ::: N g are independent standard normal variables and (i ; i) are the
verifying :
i = i i
Substituting for (2.63) in (2.13) and solving the OLE problem yields :
(2.65)
H^ (x ; ) = (x) +
N
X
i ()
i=1
p iT H (x)
i
terms, the
26
5.4.2
Variance error
(2.66)
r
X
1
i=1 i
Ti H (x)
2
As in OLE and KL, the second term in the above equation is identical to the variance
H^ (x). Thus EOLE also always under-represents the true variance. Due to the form of
(2.66), the error decreases monotonically with r , the minimal error being obtained when
no truncation is made (r = N ). This allows to dene automatically the cut-o value r
of
Remark
of
technique of
reduction
i in (2.63). This
as :
reducing the number of random variables in the shape functions method (Liu
et al.,
1986 ),
reducing the number of random variables before simulating random eld realizations (Yamazaki and Shinozuka, 1990)
EOLE
vs.
KL
The accuracy of the KL and EOLE methods has been compared by Li and Der Kiureghian (1993) in the case of one-dimensional homogeneous Gaussian random elds.
The error estimator (2.38) was computed for dierent orders of expansion
r. The results
i.e.
for a
given cut-o number r. A deeper analysis shows, as pointed out by Li and Der Kiureghian
3 Estimator (2.38) is dened as a Sup.
27
Figure 2.3: Comparison of errors for KL and EOLE methods with type A correlation
Eq.(2.35) (after Li and Der Kiureghian (1993))
6.1.2
OSE
vs.
KL
Zhang and Ellingwood (1994) applied the OSE method to a one-dimensional Gaussian random eld dened over a nite domain
fh (x)g1
i
(2.67)
where
i=0
hn (x) =
2 n + 1 x
P
2a n a
Pn is the n-th Legendre polynomial. The authors introduced two error estimators
based on the covariance function to evaluate the respective accuracy of KL and OSE
methods. To reach a prescribed tolerance, it appears that the number of terms
included in OSE is one or two more than the number of terms required by KL.
to be
28
6.2.1
The following
(2.68)
H (x) is homo-
geneous (See Part II, Chapter 2, Section 4). In the following numerical application, a
one-dimensional homogeneous Gaussian random eld having the following characteristics
is chosen :
Domain
= [0 ; 10],
Correlation length
6.2.2
a = 5.
Figure 2.4 represents the estimator (2.68) for the three discretization schemes at dierent orders of expansion. On each gure, the mean value of
"rr(x) over
is also given.
As expected from the properties of the Karhunen-Love expansion described in Section 5.2.2, the KL approach provides the lowest mean error. The EOLE error is close
to the KL error while the OSE error is slightly greater (20 points were chosen for the
EOLE discretization, which means that the size of each element in the EOLE-mesh is
LRF
a=10). As already stated by Li and Der Kiureghian (1993), the point-wise vari-
sian random eld under consideration is non dierentiable because of the exponential
autocorrelation function.
6.2.3
The results for exponential square autocorrelation function (see Eq.(2.36)) are presented
in gure 2.5. As there is no analytical solution to the eigenvalue problem associated
29
0.25
Order of expansion : 2
Mean Error over the domain :
KL : 0.231
EOLE : 0.233
OSE : 0.249
0.4
0.35
0.3
Order of expansion : 4
Mean Error over the domain :
KL : 0.112
EOLE : 0.117
OSE : 0.132
KL
EOLE
OSE
0.2
KL
EOLE
OSE
rr(x)
rr(x)
0.15
0.25
0.2
0.1
0.15
0.1
0.05
0.05
10
10
10
x
0.14
Order of expansion : 6
Mean Error over the domain :
KL : 0.0729
EOLE : 0.0794
OSE : 0.0931
0.18
0.16
0.14
KL
EOLE
OSE
Order of expansion : 10
Mean Error over the domain :
KL : 0.0427
EOLE : 0.0521
OSE : 0.0666
0.12
0.1
KL
EOLE
OSE
0.08
rr(x)
rr(x)
0.12
0.1
0.06
0.08
0.06
0.04
0.04
0.02
0.02
10
Figure 2.4: Point-wise estimator for variance error, represented for dierent discretization
schemes and dierent orders of expansion (
with the Karhunen-Love expansion in this case, only EOLE and OSE are considered.
It appears that EOLE gives better accuracy in this case also.
6.2.4
vs.
order of expansion
The mean of
functions.
As expected, at any order of expansion, the smallest mean error is obtained by KL
(if applicable). EOLE is almost always better than OSE. The EOLE-mesh renement
necessary to get a fair representation depends strongly on the autocorrelation function,
30
0.35
0.03
EOLE
OSE
Order of expansion : 2
Mean Error over the domain :
EOLE : 0.081
OSE : 0.105
0.3
0.25
EOLE
OSE
Order of expansion : 4
Mean Error over the domain :
EOLE : 0.0017
OSE : 0.0048
0.025
rr(x)
rr(x)
0.02
0.2
0.015
0.15
0.01
0.1
0.005
0.05
10
10
Figure 2.5: Point-wise estimator for variance error, represented for dierent discretization
schemes and dierent orders of expansion (
0.45
0.4
0.4
0.35
0.25
0.25
0.2
0.2
0.15
0.15
0.1
0.1
0.05
0.05
10
12
14
16
Order of expansion r
a -
EOLE
OSE
0.3
KL
EOLE
OSE
0.3
rr(x) dx / ||
rr(x) dx /||
0.35
b-
5
6
7
Order of expansion r
10
tures
type is considered (gure 2.7-b), then EOLE is more accurate than OSE whatever the
mesh renement.
case of the exponential square autocorrelation function is much smaller than that obtained for the exponential autocorrelation function, whatever the discretization scheme.
For
31
0.1
0.3
EOLE, LRF / a = 1/2
EOLE, LRF / a = 1/4
EOLE, LRF / a = 1/6
EOLE, L / a = 1/8
0.25
EOLE, L / a = 1/2
RF
EOLE, L / a = 1/4
RF
EOLE, LRF / a = 1/6
EOLE, LRF / a = 1/8
0.09
0.08
RF
OSE
OSE
0.07
rr(x) dx / ||
rr(x) dx / ||
0.2
0.06
0.05
0.04
0.15
0.1
0.03
0.02
0.05
0.01
0
a -
8
10
Order of expansion r
12
14
16
b-
5
6
7
Order of expansion r
10
and OSE
6.2.5
Conclusions
The series expansion discretization schemes (KL, EOLE and OSE) all ensure a rather
small variance error as soon as a few terms are included.
When the
gives the best accuracy. EOLE is more accurate than OSE if the underlying mesh is
i.e. LRF =a
suciently rened (
(1993), EOLE is more ecient with a ne mesh and a low order of expansion than with
a rough mesh and a higher order of expansion.
When the
exponential square
in this case. Generally speaking, the variance error computed with an exponential square
autocorrelation function is far smaller than that computed for the exponential case.
Thus in practical applications, if there is no particular evidence of the form of the
autocorrelation function, the exponential square form should be preferred, since it allows
practically an exact discretization (mean variance error <
This result holds for both EOLE and OSE discretization schemes. Furthermore, this
form of autocorrelation function implies a dierentiable process, which would be more
realistic for most physical processes.
32
The case of
the case
of a Gaussian eld :
(2.69)
(2.70)
which are of great importance for modeling material properties due to its non-negative
domain.
Although it has not be used in the literature, the translation procedure could be applied
with any of the series expansion schemes described in the last section including KL and
OSE.
e.g. the MP, SA, SF, OLE, EOLE methods. A critical parameter
for ecient discretization is the typical size of an element or the grid size.
Several authors including Der Kiureghian and Ke (1988) and Mahadevan and Haldar
(1991) have pointed out that the nite element- and the random eld meshes have to be
designed based upon dierent criteria. Namely :
stress gradients
of the
response. Should some singular points exist (crack, edge of a rigid punch, ...), the
mesh would have to be locally rened.
LRF
correlation
33
LRF
(2.71)
a4 to a2
This range was conrmed by Li and Der Kiureghian (1993) (see details in section 4)
by computing the error estimator (2.38) and by Zeldin and Spanos (1998) by
comparing the power spectra of
In the context of reliability analysis (see Chapter 4), Der Kiureghian and Ke (1988)
and Mahadevan and Haldar (1991) reported numerical diculties of the procedure
when the length
LRF
pearing in the discretization are highly correlated and the diagonalization of the
associated covariance matrix leads to numerical instabilities.
As far as the EOLE method is concerned, the short study presented in section 6
a=10 and a=5 for the exponential autocorrelation function and between a=4 and a=2 for other cases. However, in contrast to point- and average discretization methods, the fact that LRF
is rather small does not imply that the number of random variables r used in the
discretization is large, since r is prescribed as an independent parameter.
LRF
structed on a regular pattern (segment, square, cube). However, in the context of reliability analysis, Liu and Liu (1993) showed that the renement of the random eld mesh
should be connected to the
This seems to be a common feature with the nite element mesh : when the response
quantities of interest are localized in a specic subdomain of the system, it is possible
to choose a coarse mesh in the regions far away from this subdomain.
In the applications, some authors simply construct the random eld mesh by
grouping
several elements of the nite element mesh in a single one (see Liu and Der Kiureghian
(1991 ); Zhang and Der Kiureghian (1993, 1997)). This allows to reduce dramatically
the size of the random vector
applications, where the nite element mesh is generally automatically generated, having
variable element size with unprescribed orientation. Indeed, in this case, it would not
be practical to dene the random eld mesh by grouping elements of the nite element
mesh.
34
9 Conclusions
This chapter has presented a review of methods for discretization of random elds that
have been used in conjunction with nite element analysis. Comparisons of the eciency
of these methods found in the literature have been reported, and new results regarding
the series expansion methods have been presented. The question of the design of the random eld mesh has been nally addressed. As a conclusion, advantages and weaknesses
of each method are briey summarized below :
Methods yielding
they provide more accurate representations for the same mesh renement.
The SA method is practically limited to Gaussian elds since the statistics of the
random variables involved in the discretization cannot be determined in any other
case. However, it may be extended to non Gaussian elds obtained by translation
of a Gaussian eld, see Section 7.
former is the most ecient in terms of the number of random variables required
for a given accuracy. However, it requires the solution of an integral eigenvalue
problem. The latter uses
correlated
Although applicable to any kind of eld, both of these methods are mainly ecient
for Gaussian random elds, since the variables involved in the discretization are
Gaussian in this case. As an extension, non Gaussian random elds obtained by
translation can be dealt with, see Section 7.
It is possible to reduce the number of random variables involved in a discretization procedure by a spectral decomposition of their covariance matrix. Only those
eigenvalues with greatest value are retained in the subsequent analysis. This reduction technique has been applied in conjunction with SF. Coupled with OLE, it
yields EOLE.
9. Conclusions
35
eld, involving its autocovariance function or the orthogonal basis functions respectively. These realizations are continuous. In terms of accuracy, EOLE is better
than OSE if the random eld mesh is suciently rened (
LRF
The covariance matrix of the random variables involved in EOLE is readily available, since these variables correspond to selected points in the domain
. Solving
Chapter 3
Second moment approaches
1 Introduction
Historically, probability theory was introduced in mechanics in order to estimate the
response variability
value when the input parameters themselves vary around their means. The aim is to
understand how uncertainties in the input
propagate
For this purpose, second order statistics of the response are to be evaluated.
Suppose the input randomness in geometry, material properties and loads is described
by a set of
random variables, each of them being represented as the sum of its mean
i . The
element analysis, the second moment methods aim at evaluating the statistics of the
nodal displacements, strains and stresses from the mean values of the input variables
and the covariance matrix of
The
perturbation method
with the discretization scheme presented in Chapter 2, Section 3.2 is presented in Section 4. The
quadrature method
38
for non-linear dynamic problems. It uses a Taylor series expansion of the quantities involved in the equilibrium equation of the system around their mean values. Then the
coecients in the expansions of the left- and right-hand sides are identied and evaluated
by perturbation analysis.
In the context of nite element analysis for quasi-static linear problems, the equilibrium
equation obtained after discretizing the geometry generally reads :
K U = F
(3.1)
F are varying around their mean. As a consequence, the three quantities appearing
o
in the above equation will also vary around the values K o ; U ; F o they take for these
vector
read
(3.2)
K = Ko +
(3.3)
U = U
(3.4)
F = Fo +
o+
N
X
i=1
N
X
i=1
N
X
i=1
K
U
I +
i i
I
i i
N X
N
1X
K II + o(k k2)
2 i=1 j =1 ij i j
N X
N
1X
II + o(k k2 )
U
+
2 i=1 j =1 ij i j
I +
i i
N X
N
1X
F II + o(k k2)
2 i=1 j =1 ij i j
(3.5)
(3.6)
K Ii = @@K =0
i
2K
@
K IIij = @ @ =0
i
By substituting (3.2)-(3.4) in (3.1) and identifying the similar order coecients on both
1 For the sake of simplicity, the dependency of the random variables
is not written in the sequel.
39
U o = Ko 1 F o
U Ii = K o 1 F Ii K Ii U o
U IIij = K o 1 F IIij K Ii U Ij K Ij U Ii K IIij U o
(3.7)
(3.8)
(3.9)
E [U ] U
(3.10)
N X
N
1X
+
U II Cov [i ; j ]
2 i=1 j =1 ij
Cov [U ; U ]
N X
N
X
i=1 j =1
U Ii U Ij T Cov [i ; j ] =
@U T
Cov [i ; j ]
@
@
=0
=0
i
j
i=1 j =1
N X
N
X
@
ij =
(3.12)
(i ; j ) :
Cov [i ; j ]
i j
Cov [U ; U ]
(3.13)
@U T
ij i j
@
@
=0
=0
i
j
i=1 j =1
N X
N
X
@
It is seen that each term of the summation involves the sensitivity of the response to the
parameters
(partial derivatives of
i ). Eq.(3.13) thus allows to interpret what quantities are most important with respect
to the response variance. The second-order approximation of the covariance matrix can
also be derived. It involves up to fourth moments of
et al. (1986b ).
40
et al.
settlement of a pile embedded in an elastic soil layer with random Young's modulus.
et al.
ab
(1986 , )
applied the perturbation method to static and dynamic non-linear problems, including :
et al., 1986b );
et al., 1986a );
turbine blade (shell element) with random side load and length (same reference).
The results compare fairly well with Monte Carlo simulations. However the coecient
of variation of the random quantities is limited to 10% in all these examples.
ab
and Takada (1990 , ), couples the perturbation method with the representation of the
stochastic stiness matrix presented in Chapter 2, Section 3.2. This representation can
be put in the following form :
(3.14)
k =k
e
e+
o
NW
XI
l=1
kel el
41
el are zero-mean weighted integrals of the random eld and kel are deterministic
matrices. By assembling these contributions over the N elements of the system, a global
stochastic stiness matrix involving NW I N random variables is obtained.
where
U is used :
N NW
X
XI
@U
U = Uo +
e
(3.15)
e=1 l=1
@el
el =0
U =U
(3.16)
N NW
X
XI
e=1 l=1
el K o 1
@K
@el
Uo
Cov [U ; U ] =
N X
N NWI
X
X NW
XI
e1 =1 e2 =1 l1 =1 l2 =1
Ko
1
@K
@el11
U U
o
oT
@K T
@e2
l2
2
2
The last term in the above expression is obtained from the denition (2.32) of the
weighted integrals. For example, in the one-dimensional case :
(3.18)
Cov el11 ;
el22
E el11 el22 =
e1
e2
where :
(3.19)
SHH (! ) satisfying :
CHH (x1
(3.20)
one nally obtains :
Cov el11 ;
(3.21)
el22
=
=
Z
Z
R
R
x2 )
SHH (! )
Z
R
SHH (! ) ei!(x1
e1
x2 ) d!
xl1 ei!x1 dx
1
SHH (! ) vle11le22 (! ) d!
e2
x2 l2 ei!x2 dx2
42
where
vle11le22 (! ) is dened as :
vle11le22 (! ) =
(3.22)
e1
xl1 ei!x1 dx
e2
x2 l2 ei!x2 dx2
Substituting (3.21) in (3.17) and gathering the diagonal terms into the vector
Var [U ]
leads to :
Var [U ] =
(3.23)
Z
R
V (!) is a vector
SHH (! ) V (! ) d!
having
K o , k
e,
l
Eq.(3.22). These functions can be given closed-form expressions after some algebra in
Deodatis (1990) examined the following upper bound for the variance of
Ui .
Taking
indeed each component of (3.23) separately, the following trivial inequality holds :
Var [Ui ]
(3.24)
where
H2
Z
R
The above equation yields an upper bound on the response variance, which is independent of the correlation structure of the input eld. This is valuable since this kind of data
is dicult to obtain in practice. However, Eq.(3.24) has limited practical use because
the quantities
max jVi(! )j are not easily available. In the application of this method to
a frame structure, Deodatis (1990) computed these maximum values graphically after
plotting the functions
Vi (! ).
tained by a nite element code. It is based on the quadrature of the integrals dening
these moments.
43
(3.25)
The
quadrature of the
E g i (X ) =
[g (x)]i fX (x) dx
E g i(X )
(3.26)
NP
X
k=1
wk [g (xk )]i
(wk ; xk ) are integration weights and points associated with fX respectively, and
NP is the order of the quadrature scheme. For instance, if X has a uniform distribution
over [ 1 ; 1], the integration points in the above equation are the well-known Gauss
points. More generally, it is possible to compute integration weights and points associated with other PDFs fX at any order. Baldeweck (1999) gives tables for normal and
where
Xi 's
(3.27)
E S i (X
1;
::: XN ) =
RN
stress component) is to
'(:) is the PDF of a standard normal variate, and s(x1 ; ::: ; xN ) is usually known in
an algorithmic sense, i.e. through a nite element code. Generalizing Eq.(3.26), Eq.(3.27)
where
can be estimated as :
(3.28)
E S i (X
It is seen that
NPN
1;
::: XN )
NP
X
k1 =1
evaluations of
NP
X
kN =1
S (i.e.
a priori.
44
negligible compared to others, so that some terms in Eq.(3.28) are not computed.
Remark
From Eq.(3.27), it is seen that the method is not a pure second moment
approach, since information about the distributions of the basic random variables is
included in the calculation. However, it is presented in this chapter because it gives
primarily the moments of the response quantities.
Baldeweck applied the quadrature method to compute the rst four moments of the
response quantities. From these rst moments, a probability distribution function can
be tted (so-called Johnson or Pearson distributions). This PDF can be nally used to
get reliability results.
Several examples in structural mechanics, geotechnics and fracture mechanics are presented. In each case, the results obtained with the quadrature method compare well with
e.g.
that non-linear problems can be dealt with as easily as linear problems. However, the
limitation on the number of random variables is a severe restriction of this approach (at
most 4 were used in the applications described by Baldeweck (1999)).
Due to the Taylor series expansion, accurate results are expected only in case of
small variability of the parameters. The upper limit on the coecient of variation
(COV) for which the results are acceptable strongly depends on the degree of
nonlinearity of the system. Fair comparisons with Monte Carlo simulations have
been obtained for COV up to 20%. However the upper limit may dier a lot with
respect to the kind of mechanical problems under consideration. It is important to
i.e.
note that, while the results from perturbation analysis are distribution-free (
they only require the second moments of the input variables), the Monte Carlo
results by necessity must be obtained for specic distributions. In this sense the
comparison is dependent on the assumed distribution in the Monte Carlo analysis.
45
the derivatives of
K and F
i ,
possibly at the second order level. This derivations are usually performed at the
element level before assembling the system. In most situations, they can be done
analytically, leading however to intricate formul. These calculations can be time
consuming, particularly when the second order terms are included. Higher-order
approximations are totally intractable.
The weighted integral method described in Section 4 allows to obtain a) second order
statistics of the response given a prescribed input random eld and b) an upper bound
on the response variance which is independent of the correlation structure of the eld.
It is claimed that the method does not use any particular discretization scheme for the
random eld. However, as stated in Chapter 2, Section 3.2, the accuracy of the method
for problems in which the correlation length is small compared to the size of the structure
may not be good. This has been illustrated by Zhang and Der Kiureghian (1997, Chap. 2)
on the example of an elastic rod in tension, having constant cross-section and random
Young's modulus. Moreover, the following severe limitations of the method have to be
recognized :
it is limited to elastic structures, where the Young's modulus is modeled as a random eld. Although extension of the method to nonlinear problems was claimed
(Deodatis and Shinozuka, 1991), such a derivation could not be found in the literature.
the bounds (3.24) on the response variance are dicult to compute in practice,
due to the complex expression of the response variability functions
V (!).
e.g. NW I
NW I N
Due to its simplicity, the quadrature method is appealing for problems involving a reduced set of random variables. It is neither limited to linear problems nor small coe-
46
Chapter 4
Finite element reliability analysis
1 Introduction
Reliability methods aim at evaluating the
eling takes into account randomness. Classically, the system is decomposed into components and the system failure is dened by various scenarii about the
joint
failure of
e.g. a given structure) describing the randomness in geometry, material parameters and
loading. If needed, suppose that one of the discretization scheme described in Chapter 2
has been applied to represent random elds as functions of a nite set of random vari-
denote a vector of such eects, whose values enter in the denition of the failure of the
system. These two vectors are related through the
(4.1)
S = S ()
mechanical transformation :
48
e.g.
through a
The values of
satisfying
g (S ) = 0
dene the
of the structure.
g(S) = max , if the failure occurs when the displacement at a given point
exceeds a given threshold
max .
g(S) = o J2 (), if the failure occurs when a given point yields (o is the yield
stress and J2 ( ) the equivalent Von Mises stress).
Let us denote now by
(4.2)
Pf =
Z
g(S )0
fS (s)ds
The multi-fold integral (4.2) over the failure domain is not easy to compute.
Thus approximate methods have been developed in the last 25 years. Exhaustive presentation of this domain can be found in the monograph by Ditlevsen and Madsen (1996).
In the sequel, only the main concepts are summarized.
49
margin is dened by :
Z=R S
(4.3)
Z
Z
C =
(4.4)
(4.5)
and
were to be
Z
Z
Z
Z
jointly normal,
( C )
where
so
can
Z = g (S )
(4.7)
S and covariance matrix SS being known. The mean and covariance of Z
are not available in the general case where g (S ) is non-linear. Using the Taylor expansion
around the mean value of S :
(4.8)
Z = g (S ) + (rs g )TS=S (S S ) + o k (S S )2 k
the mean
Z g(S )
Z2 (rs g )TS=S SS (rs g )S=S
(4.11)
MVFOSM = q
g (S )
50
g 3 ()). Variations can be important in some problems (Ditlevsen and Madsen, 1996,
chap. 5).
The problem of invariance was solved by Hasofer and Lind (1974) in a second moment
context by recasting the problem in the standard space using a linear transformation of
random variables. Essentially, the authors showed that the point of linearization should
be selected as the point on the limit state surface nearest to the origin in the standard
space, the distance to this point being the rst-order second moment reliability index
FOSM . Later, a non-linear probability transformation was employed to account for probability distribution of the variables including non Gaussian distributions. This method
is described in the sequel.
Y = Y ()
(4.12)
such that
is a
Gaussian random vector with zero mean and unit covariance matrix.
:
. Several cases sorted by ascending order of diculty are listed in the sequel
as examples :
is a Gaussian random vector with mean and covariance matrix . The
diagonalization of the symmetric positive denite matrix
(4.13)
= A Y +
allows to write :
Y ()
=
Jy ; =
A 1 ( )
A1
(4.16)
51
f (xi )
Jy ; = diag
'(yi )
(4.17)
ij =
(4.18)
Cov [i ; j ]
i j
correlations was solved by Der Kiureghian and Liu (1986). The authors proposed
two dierent models :
the
the closed form expression for the joint PDF and CDF become tedious to
manage when dealing with a large number of random variables.
the
random variables and complies with almost any valid correlation structure.
Due to these characteristics, only the Nataf model is presented now. From the
marginal PDF of
(4.19)
(4.20)
is dened :
fZ (z ) = 'n (z ; Ro )
(2 )n=2 det Ro
exp
1 T
z Ro 1 z
2
(4.21)
then reduces to :
'n(z ; Ro )
'(z1 ) : : : '(zn )
(4.22)
ij =
Z 1 Z 1
x
i
i
xj
j
j
52
Der Kiureghian and Liu (1986), Liu and Der Kiureghian (1986), Ditlevsen and
Madsen (1996).
From a reliability point of view, assuming the basic variables
PDF, vector
writes :
(4.23)
Y () = Lo 1 Z
= Lo 1 1 [F1 (x1 )] ; ::: 1 [Fn (xn )] T
J y ; = Lo
(4.24)
1 diag
f (xi )
'(yi )
as follows :
8
>
y1
>
>
>
<y
=
2 =
>
:::
>
>
>
:
yn =
(4.25)
1 [F (x1 )]
1 [F (x
2 jx1 )]
1 [F (xn jx1 ;
::: xn )]
i (Ditlevsen and
Madsen, 1996).
g (S ) g (S ()) = g (S Y 1 (Y )) = G(Y )
(4.27)
where
(4.28)
Pf =
Z
G(y )0
'(y) dy
1
2
ky
k2
53
This PDF has two interesting properties, namely it is rotationally symmetric and decays
exponentially with the square of the norm
origin of the
gure 4.1
reliability index
(4.29-a)
(4.29-b)
= T y
y = argmin fk y k
j G(y) 0g
This quantity is obviously invariant under changes in parametrization of the limit state
function, since it has an intrinsic denition,
state surface.
The solution
i.e.
the
point or the most likely failure point in the standard normal space. When the limit state
function G(y ) is linear in y , it is easy to show that :
Pf = ( )
(4.30)
where
yn
y
G(y) = 0
y1
When
54
by solving (4.29),
Pf
(4.31)
Pf
= ( )
Geometrically, this is equivalent to replacing the failure domain by the halfspace outside
the hyperplane tangent to the limit state surface at
becomes a better approximation when
is large.
curvature tting
tting where semi-paraboloids interpolate the limit state surface at given points around
the design point (Der Kiureghian et al., 1987).
Recently, higher order approximation methods (HORM) have been proposed by Grandhi
and Wang (1999), where the limit state surface is approximated by higher order polynomials. The amount of computation needed appears to be huge compared to the improvement it yields.
Early approaches
As mentioned earlier, the classical reliability methods (FORM/SORM) require the determination of the
design point, which is dened as the point on the limit state surface
closest to the origin, in the standard normal space. The constrained optimization problem (4.29) is equivalent to :
(4.32)
y = argmin Q(y) = 1
2
ky
k2
j G (y ) 0
(4.33)
L(y ; ) = 12 k y k2 + G(y)
(y ; ) = argmin L(y ; )
55
y + rG(y) = 0
G(y ) = 0
(4.35)
(4.36)
positive Lagrange multiplier can be obtained from (4.35), then substituted in the
same equation. This yields the rst-order optimality conditions :
The
k rG(y) k y + k y k rG(y) = 0
(4.37)
This means that the normal to the limit state surface at the design point should point
towards the origin.
Hasofer and Lind (1974) suggested an iterative algorithm to solve (4.34), which was later
used by Rackwitz and Fiessler (1978) in conjunction with probability transformation
techniques. This algorithm (referred to as HLRF in the sequel) generates a sequence of
points
T
rG(yi)
yi+1 = rG(yki)rG (yyi ) kG(yi) k r
G(y ) k
(4.38)
Eq.(4.38) can be given the following interpretation : at the current iterative point
limit state surface is linearized,
tangent to
G(y)
at
y = yi.
yi, the
orthogonal projection
of
yi
tangent hyperplane.
As the limit state function and its gradient is usually dened in the original space, it
is necessary to make use of a probabilistic transformation such as those described in
Section 2.4. The Jacobian of the transformation is used in the following relationship :
ry G(y) = rg() J ; y
(4.39)
The HLRF algorithm has been widely used due to its simplicity. However it may not
converge in some cases, even for rather simple limit state functions. Der Kiureghian and
de Stefano (1991) have shown that it certainly diverges when a principal curvature of
the limit state surface at the design point satises the condition
modied versions of this algorithm have been developed (Abdo and Rackwitz, 1990;
Liu and Der Kiureghian, 1991 ). The latter reference presents a comprehensive review
of general purpose optimization algorithms, including the gradient projection method
(GP), the augmented lagrangian method (AL), the sequential quadratic programming
method (SQP), the HLRF and the modied HLRF (mHLRF). All these algorithms have
been implemented in the computer program CALREL (Liu
56
purposes, and tested with several limit state functions. It appears that the most robust
as well as ecient methods are SQP and mHLRF.
Although the modied mHLRF was an improvement over the original HLRF, no proof
of its convergence could be derived. Thus further work has been devoted in nding an
yi+1 = yi + i di
(4.40)
(4.41)
with
i = 1
T
rG(yi) y
di = rG(yki)rG (yyi ) kG(yi ) k r
i
G(y i ) k
i
In the above equations, di and i are the search direction and the step size respectively.
The original HLRF can be improved by computing an optimal step size i =
6 1.
For this purpose, a merit function m(y ) is introduced. At each iteration, after comput(4.42)
that is :
i = argmin fm(yi + di )g
(4.43)
This non-linear problem is not easy to solve. It is replaced by the problem of nding a
i such that the merit function is suciently reduced (if not minimal). The so-called
Armijo rule (Luenberger, 1986) is an ecient technique. It reads :
(4.44)
i = max bk j m(yi + bk di ) m(y i ) abk k rm(yi ) k2 ; a; b > 0
value
k 2N
Zhang and Der Kiureghian (1995, 1997) proposed the following merit function :
m(y) =
(4.45)
1
2
k y k2 +c jG(y)j
It attains its minimum at the design point provided the same condition on
is
fullled.
Both
properties
are
sucient
to
ensure
that
the
global
algorithm
dened
by
57
Conclusion
(see Section 3 for a detailed presentation). Thus an ecient optimization algorithm for
determining the design point should call the smallest number of each evaluation. From
this point of view, the iHLRF algorithm is the most ecient algorithm. The reader is
referred to Liu and Der Kiureghian (1991 ) and Zhang and Der Kiureghian (1997) for
detailed cost comparisons.
ry G(y)
in terms of load eects, which are related to the basic random variables
. The chain
In this expression,
is given by
of
sensitivity are useful in various applications such as optimal structural design and determination of importance of parameters.
For our purpose of determining the design point, the evaluation of the gradient has to be
ecient (because of the numerous calls in the iHLRF algorithm) and accurate (because
its value enters an iterative convergent procedure, which is driven by tolerance checking).
The straightforward application of a
may be inappropriate in
the size of the nite variation of the parameters and is thus dicult to x in advance.
A computationally more ecient approach employs the
1986 ). Recalling the formalism developed in Chapter 3, the rst order variation of the
58
U Ii is given by :
U Ii = K o 1 F Ii K Ii U 0
(4.47)
Thus the same
mean value stiness matrix K o is used for evaluating all the components
of the displacement gradient vector. This method requires basically one complete nite
i = 1 ; ::: N .
been proposed
to circumvent the drawbacks of the above methods. It is presented in the sequel rst
for an elastic linear problem. Then the extension to geometrically non-linear structures
and Der Kiureghian, 1993), plane stress elastoplastic damaged structures (Zhang and
Der Kiureghian, 1997; Der Kiureghian and Zhang, 1999) will be briey summarized.
R=
(4.48)
[Z
e
B d
e =
t
[ Z
e
writes :
N obo d
e +
T
Z
@
e
N to dSe = F
T
where
m ; l ; g be those vari-
ables representing uncertain material properties, external loads and the geometry of the
structure respectively. Let
displacements. In case of small strain linear elastic structures, the following relationships
hold :
(4.49)
(4.50)
(4.51)
(4.52)
X
U
R
F
=
=
=
=
X [g ]
U [X (g ) ; m ; l ]
R [X (g ) ; U (X (g ) ; m ; l ) ; m ]
F [X (g ) ; l ]
59
To simplify the presentation, it is assumed in this section that the limit state function
only depends on the displacement vector
U:
g (S ()) g (U ())
(4.54)
Thus its gradient becomes :
(4.55)
3.2.1
m yields :
@R
@R
(4.56)
r
m U + = 0
@U
@ m
@R
we obtain from the above equation :
Introducing the stiness matrix K =
@U
(4.57)
K rm U = @@R
m
Dierentiating (4.53) with respect to
m gives :
@R [
T @ d
=
B
@ m e
e
@ m e
(4.58)
= D(m) B ue
Thus :
(4.60)
@R [
T @ D B u d
=
B
e
e
@ m e
e
@ m
The right hand side of (4.60) is then obtained by evaluating derivative quantities in each
element, then assembling them in a global vector exactly as a regular vector of nodal
forces.
by a random
60
em
of
m
represents the
constant
Young's modulus in
(4.61)
@R [
T D B u d
=
B
o
e
e
@ m e
e
(4.62)
x 2
e can be written as
(4.63)
D(m ; x) = (A1 + A2 m ) Do
thus depending linearly on the random variables (see Eq.(2.65) for the exact exprespoint
@R [
=
B T A2 Do B ue d
e
@ m e
e
(4.64)
3.2.2
l yields :
K rl U = @@F
l
Usually
@F [
=
@ l e
Z
@@o bo d
e + N T @@to dSe
l
l
@
e
l contains load intensity factors for point-wise, surface or body forces. Hence
3.2.3
g yields :
@R @X
@U @F @X
+K =
@X @ g
@ g @ X @ g
(4.67)
which simplies in :
(4.68)
K rg U
@F
=
@X
@R
@X
@@X
@X
@ g
61
g . The dierence inside the brackets should however be paid more attention, since the
domain of integration
e in (4.48) is dependent on X . To carry out these derivatives, it
is necessary to map the integral domains onto a xed conguration. Such a mapping is
a standard scheme for the so-called isoparametric elements. The derivation of all these
quantities for truss and four-node plane elements can be found in Liu and Der Kiureghian
(1989).
3.2.4
can be
written as :
(4.69)
linear systems (
ri U = K
(4.70)
1
rU
@F
@i
@R
@i
rU g(U ) rU
K = rU g
(4.71)
Then the following equality holds :
rg(U ()) = rU g K
T
(4.72)
=
@@F
1
@F
@
@R
@
@R
@
inverse stiness matrix used in solving (4.71) is readily available from the initial
@ F =@
and
@ R=@
62
The assumption (4.54) is now relaxed. If strains or stresses appear in the limit state
function
(S = U ; ; "),
@B
@ @D
=
B
u
+
D
e
@ @
@
(4.73)
ue + D B @@ue
(m ; l ; g ).
rsg rS
3.2.5 Examples
Der Kiureghian and Ke (1988) considered the reliability of elastic structures, namely a
beam with stochastic rigidity and applied load. Both random elds were assumed to be
homogeneous and Gaussian. They were discretized by the MP method on the one hand,
the SA method on the other hand. These methods respectively over- and under-represent
the variance of the original eld. Thus they allow to bound the exact results (within the
framework of FORM approximations).
Two limit state functions were dened, one in terms of midspan deection and another
in terms of bending moment. The reliability index is computed with dierent correlation
lengths and element sizes. It appears that good accuracy is obtained when the ratio
between element size and correlation length is
of MP and
EI . The reliability index was computed for various truncated expansions of the eld involving an increasing number of terms M . The convergence
is attained for M =6-10 depending on the choice of the autocovariance function. Using
i.e. analytical expressions for its eigenfunctions), the convergence was obtained
63
replaced by :
(4.74)
(4.75)
Included in the above reference is the derivation of all partial derivatives of interest for
truss and plane four-node elements. A summary of the procedure can also be found in
Liu and Der Kiureghian (1991 ). It is emphasized that the non-linearities do not make
the gradient computation more expensive, provided the stiness matrix is replaced by
the
any additional iteration. In conjunction with the gradient method, the gradient response
is also obtained by a single forward resolution. However the analytical expressions for
the partial derivatives and their coding is much more cumbersome in the non-linear case.
In the above reference, the gradient operators were coded in FEAP (Zienkiewicz and
Taylor, 1989). The reliability of a square plate with a hole having random geometry was
investigated. The applied load had a random intensity. The elastic material properties
M () U (t ; ) + C () U_ (t ; ) + R [U (t ; ) ; ] = P (t ; )
where the dots denote the time derivatives. The response gradient with respect to
is
denoted by ;
(4.77)
V = @@U
2 For the sake of simplicity of the notation, only one sensitivity parameter
64
M V + C V_ + K (U ) V = @@P
(4.78)
@M
U
@
@C _
U
@
@R
@ U
xed
U n+1 = L1 U n+1 ; U n ; U_ n ; U n
U_ n+1 = L2 U n+1 ; U n ; U_ n ; U n
(4.79)
(4.80)
When substituting for (4.79), (4.80) in (4.76), a non-linear system of equations is obtained, which is usually solved by a Newton-Raphson scheme. When convergence is
achieved, the gradient vector
linear case (see Section 3.2). The matrix used in this substitution is identical to that
used for determining
The quantities appearing in the right hand side contain the derivatives of the internal
forces
(such
as plasticity), this derivative is complicated. Zhang and Der Kiureghian (1993) present
the complete derivation in the case of
J2
hardening.
The rst example presented in the above paper is a plane strain analysis of a strip with
a circular hole. Cyclic loading under quasi-static conditions was applied. The response
gradient with respect to yield stress and hardening parameters was computed using
DDM and nite dierence with various nite variation size. Convergence from the latter
to the former when variation size tends to zero was assessed.
The second example is a truss structure under dynamic loading. Geometrical nonlinearities due to large displacements are taken into account. The material of the truss
follows
In both cases, the CPU time for the gradient computation is a small fraction of the
time required for the response run. This result is applicable to any number of variables
provided the adjoint method is used.
J2
plasticity. In this case, the discretized constitutive laws take a special form
because of the zero stress constraint. A coupling with the damage model by Lemaitre
and Chaboche (1990) is introduced and the sensitivity formulas developed.
4. Sensitivity analysis
65
system
locations of the rst damage threshold crossing. A similar study is presented by Der
Kiureghian and Zhang (1999). It is emphasized that taking into account the spatial
variability in reliability dramatically changes the result,
4 Sensitivity analysis
The determination of the design point requires the computation of the gradient of the
mechanical response. This can lead to tedious analytical developments and coding when
the direct dierentiation method is used, but allows for addressing strongly non-linear
reliability problems. However, the gradient computation contains information which can
be used for
sensitivity analysis.
analysis.
The relative importance of the basic standard normal random variables entering the
reliability analysis can be measured by means of the vector
= k yy k
(4.81)
where
dened as :
y
denotes the coordinates of the design point in the standard normal space.
with respect to parameters entering the denition of the limit state function (g ( ; g )) or the probability
distribution function f ( ; f ) of the basic random variables . In the former case, the
sensitivity of is (Ditlevsen and Madsen, 1996, chap. 8) :
(4.82)
d
=
dg
1
@g ( ; g )
k ry G(y(g ) ; g ) k @g
In the latter case, it involves the partial derivative of the probability transformation
(4.12) and turns out to be :
(4.83)
@ Y ( (f ) ; f )
d
= (f )T
df
@f
66
where
obtained as :
dPf
d
= '( )
d
d
(4.84)
The papers referred to in Section 3 all include sensitivity analysis of the reliability index,
which gives better insight of the problems under consideration.
A special use of the sensitivity analysis is the following. Suppose the limit state function
is dened in terms of one load eect, say, one component of displacement
uio
and a
threshold :
g (U ) = uio
(4.85)
Obviously, the probability of failure associated with this limit state function is
to the cumulative distribution function of
ui
uio
Pf
evaluated at
with respect
identical
reasoning was applied in Liu and Der Kiureghian (1991 ) and Zhang and Der Kiureghian
(1997) to determine the complete PDF of a response quantity.
Sensitivity analysis can also be used to identify random variables whose uncertainty has
insignicant inuence on the reliability index and which can be replaced by deterministic
values (
0: 3
deserve to be modeled as random elds for a better accuracy of the results. The
examples considered to determine this empirical value of 0.3 included a clamped beam,
a portal frame and a two-dimensional plate with a hole.
et al.
(1997) entirely
67
nite element computations have been presented. It has been observed that the two most
important issues making this marriage possible are :
a practical method for computing gradients of the limit state function. The
dierentiation method
direct
(N + 1)
random variables).
As a consequence, when a large number of random variables is used together with the
nite dierence method (for instance, when a commercial nite element code is used, the
source code of which not being accessible), the direct approach may eventually be really
time consuming (Lemaire, 1998). The
oers an alternative in
this case.
Let
(4.86)
g (x) g^(x) = ao +
N
X
i=1
ai xi +
N
X
i=1
aii
x2 +
i
N X
N
X
i=1 j =1;j 6=i
aij xi xj
68
i.e. number
of nite element runs) is required to build the surface. Then the reliability analysis
can be performed by means of the analytical expression (4.86) instead of the true limit
state function. This approach is particularly attractive when simulation methods such
as importance sampling (Bucher and Bourgund, 1990) are used to get the reliability
results.
"rr(a) =
(4.87)
NF
X
yk
k=1
g^(xk ) 2
(4.89)
Find
a = Argmin
( NF
X
k=1
yk
2
V x a
T ( k)
After basic algebra (see for instance Faravelli (1989)), the solution of the above problem
turns out to be :
a = VT V
1
VT y
k
where V is the matrix whose rows are the vectors V (x ) (see Eq.(4.88)) and y is the
k
k
vector whose components are y = g (x ).
(4.90)
The various response surface methods proposed in the literature dier only in the terms
be able to solve (4.90). Furthermore, the tting points have to be chosen in a consistent
way in order to get independent equations,
3 The subscripts
i.e. an invertible V T V .
(i; j ) vary as described in Eq.(4.86). The variation is not explicitly written here for
69
Wong selected values symmetrically around the mean at a distance of one standard
deviation, that is :
xi = i i
(4.91)
The number of tting points increases exponentially with the number of random variables
Bourgund (1990) proposed a simplied quadratic expression without cross terms, which
is dened by only
8
1
>
>
<
2i
>
>
: 2i+1
x
x
x
(4.92)
i-th
i = 1 ; ::: N
i = 1 ; ::: N
ei
i-th basis
vector of the space of parameters, whose coordinates are f0 ; ::: ; 1 ; 0 ; ::: g, and f is an
where
i
= X
= X f i ei ;
= X + f i ei ;
random variable,
is the
x is computed. Then
a new center point xM is obtained as a linear interpolation between X and x , so that
(4.93)
A second response
requires only
out for structural systems involving a great number of random variables. Importance
sampling is nally used to get the reliability results.
Rajashekhar and Ellingwood (1993) later considered the approach by Bucher and Bourgund (1990) as the rst two steps of an iterative procedure they pushed forward until
convergence. They also added cross terms to the response surface denition, obtaining
better results in the numerical examples.
70
Analyzing the three papers presented above, Kim and Na (1997) observed that, in each
e.g.
the basic random vector) and arranged along the axes or the diagonals of the space
of parameters, without any consideration on the orientation of the original limit state
surface. The authors argued that these procedures may not converge to the true design
point in some cases.
Alternatively, they proposed to determine a series of
In each iteration, the tting points used in the previous step are projected onto the
previous response surface, and the obtained projection points (which lie closer to the
actual limit state surface) are used for generating the next response surface. In each
iteration, an approximate reliability index is readily available, since the response surface
is linear. In some sense, this method nds the design point without solving the minimization problem usually associated with FORM. The authors assessed the validity of
this so-called
1,000,000 samples), rst by using an analytical limit state function, then by studying a
frame structure and a truss.
Starting from the paper by Kim and Na (1997), Das and Zheng (2000) recently proposed
to enhance the linear response surface by adding square terms. The tting points dening
the nal linear response surface are reused to produce the quadratic surface. SORM
analysis is then performed.
Lemaire (1997) presents a synthetic summary of the response surface methods (called
adaptive because of successive renement until convergence around the design point)
and draws the following conclusions :
it is better to cast the response surface in the standard normal space rather than
in the original space for reliability problems. All quantities being adimensional,
there is a better control of the regression.
Provided enough tting points are used, the choice of the type of experimental
design is not fundamental.
The quality of the response surface has to be checked. Dierent indicators are
proposed to estimate the accuracy, including :
the back-transformation of the tting points from the standard normal space
to the original space, in order to exclude non physical points,
V T V appearing in Eq.(4.90),
the belonging of the obtained design point to the original limit state surface.
71
In order to reuse at best the nite element results, a data base keeping track of the nite
element runs should be constructed.
et al. (1998) and Pendola et al. (2000c ) proposed a benchmark problem in non-
linear fracture mechanics. Crack initiation in a steel pipe submitted to internal pressure
and axial tension is under consideration. Dierent nite element codes including Ansys
and Code_Aster
(developed
by Lemaire and his colleagues), Comrel (developed by Rackwitz and his colleagues) and
Proban (developed by Det Norske Veritas). As far as accuracy is concerned, the direct
coupling and the response surface method give identical results for probabilities of failure
divide by 10 the number of nite element runs for the specic example. However, a nite
dierence scheme for gradient computation was applied in the direct coupling, which is
not optimal. In this example, an axisymmetric non-linear nite element model was used.
A similar comparison has been carried out by Defaux and Heining (2000) on the problem
of an hyperbolic cooling tower submitted to thermal and wind loading. A linear elastic
three-dimensional nite element model using thin shell elements was used. In this case,
the direct coupling and the response surface method gave the same results for similar
computational cost.
Ryfes
72
interpolation tools, and can thus be used instead of quadratic response functions to
approximate the limit state function. After being trained with a set of input/output
data (here realizations of the vector of basic random variables, and corresponding value
of the limit state function obtained after a nite element calculation), the neural network
can produce reliable output values for any input at low cost.
Hurtado and Alvarez (2000) compare two types of networks called
and
multi-layer perceptrons
radial basis functions networks. After being trained, the networks are used together
with a crude Monte Carlo simulation to get the probability of failure. A system reliability
problem associated with the collapse of a frame is considered. It appears that the radial
basis functions network provide the best results with a rather small number of training
samples.
Pendola
The neural network replaces the quadratic response surface obtained after the iterative
procedure described in section 5.4. Applying this approach to the benchmark problem
described in Pendola
et al.
those obtained by the response surface method, the number of nite element simulations
being however reduced by a factor 2.
5.7 Conclusions
Although the response surface method is an old idea, it seems to have gained new
consideration in recent years. The up-to-date approach consists in generating quadratic
response surfaces iteratively, where the center point converges to the design point. After
convergence, any reliability method can be applied with the response surface,
e.g. FORM,
6. Conclusions
73
As a summary, the response surface method appears to give accurate results for most
problems applied, and may be faster than the direct coupling when a small number
of random variables is considered, and when it is not possible to implement the direct
dierentiation method (for instance, when a commercial nite element code is used).
Otherwise, the direct coupling will probably require less or equal amount of computation. These conclusions can change in the near future due to the introduction of neural
networks in the eld of reliability analysis.
6 Conclusions
In this chapter, methods coupling reliability and nite element analysis have been presented. The classical approach of reliability (FORM/SORM) has been summarized. The
need of computing response gradients was emphasized. For this purpose, the direct dierentiation method has been presented. It allows sensitivity analysis for general problems
including material and geometrical non-linearities and dynamics. Using this approach,
the computational cost of the gradient is a small increment over the cost of the non-linear
response itself.
It has been shown that reliability analysis allows for obtaining PDFs of any response
quantity. It should be noticed that this approach will give accurate results only for the
tails of the PDF. Indeed it is based on FORM, which may be inaccurate for low reliability
indices (large probabilities).
The response surface method has been presented as an alternative to direct coupling. It
is also applicable to the most general problems and does not require the implementation
of gradients in the nite element code. Whether one method is more ecient than the
other depends fundamentally on the number of random variables included in the analysis
and the way gradients are computed.
Chapter 5
Spectral stochastic nite element
method
1 Introduction
spectral stochastic nite element method (SSFEM) was proposed by Ghanem and
Spanos (1990, 1991a ) and presented in a comprehensive monograph by Ghanem and
Spanos (1991b ). It is an extension of the deterministic nite element method (FEM) for
The
with deterministic
geometry, material properties and loading. The evolution of such a system is governed
by a set of partial dierential equations (PDE) and associated boundary conditions and
initial state. When no closed-form solution to these equations exists, a discretization
procedure has to be applied in order to handle the problem numerically. In the usual
nite element method, the geometry
that are the nodes of the nite element mesh. Correspondingly the response of the
system,
ui ; i
If a material property such as the Young's modulus is now modeled as a random eld,
76
random vector of nodal displacements U (), each component being a random variable ui() yet to be characterized.
above paragraph results in approximating the response as a
A random variable is completely determined by the value it takes for all possible out-
2 . Adopting the same kind of discretization as for the spatial part would
result in selecting a nite set of points f1 ; ::: Q g in . The Monte Carlo simulation of
comes
i have to be selected
that an accurate description of the response would require a large value for Q.
SSFEM aims at discretizing the random dimension in a more ecient way using series
expansions. Two dierent procedures are used.
the input random eld is discretized using the truncated Karhunen-Love expansion presented in Chapter 2, Section 5.2.
ui ()
nomial chaos.
Computational issues regarding the peculiar system of equations eventually obtained will be addressed in Section 3.
hybrid
SSFEM.
Some technical developments including the denition of the polynomial chaos and the
additional tools related to the discretization of lognormal random eld are gathered in
an appendix at the end of this chapter.
77
N N (N
K U =F
(5.1)
where the
matrices k
element stiness
k =
e
(5.2)
B T D B d
e
D(x ; ) H (x ; ) Do
(5.3)
where
Do
is a
constant
matrix. The
Karhunen-Love
expansion
(Eq.(2.43)) :
(5.4)
H (x ; ) = (x) +
1 p
X
i=1
i i () 'i(x)
of
H (:)
writes
78
(5.5)
where
=k
e ( )
e+
o
1
X
i=1
kei i()
keo is the mean element stiness matrix and kei are deterministic matrices obtained
by :
k = i
e
i
(5.6)
'i (x) B T D o B d
e
Assembling the above element contributions eventually gives the stochastic counterpart
of the equilibrium equation (5.1) (assuming a deterministic load vector
"
Ko +
(5.7)
1
X
i=1
F) :
K i i() U () = F
no closed-form solution for such an inverse exists. An early strategy adopted by Ghanem
and Spanos (1991 ) consists in using a Neumann series expansion of the inverse stochastic
stiness matrix to get an approximate response. Eq.(5.7) can be rewritten as :
"
Ko I +
(5.8)
which leads to :
"
U ( ) = I +
(5.9)
1
X
i=1
1
X
i=1
Ko
1K
i i ( )
# 1
K o 1 K i i()
U () = F
U 0 ; U 0 = Ko 1 F
The Neumann series expansion of the above equation has the form :
U ( ) =
(5.10)
1
X
k=0
( 1)k
" 1
X
i=1
#k
K o 1 K i i() U 0
"
U () = I
1
X
i=1
K o 1 K i i() +
1 X
1
X
i=1 j =1
K o 1 K i K o 1 K j i()j () + : : : U 0
U ().
i and k in
79
(5.12)
ui ()
where
1
X
j =0
fPj g1j=0
; F ; P)
L2(
The set of
Pj fk ()g1k=1
Pj fk ()g1k=1
(5.13)
uij
in Eq.(5.13) forms a
uij
e.g. :
= i11 i22 : : : ipp
basis
ui () in this
basis.
Referring to the inner product dened in
not orthogonal. For instance, 1 () and 13 () are two basis random variables
4
inner product is E [1 ] = 3. For further exploitation of the response, such as
is however
whose
polynomial chaos3 proposed by Ghanem and Spanos (1991b ) possesses this property.
The details of its construction are quite technical and not essential to the understanding
of SSFEM. Thus they are given in Appendix A.1 at the end of this chapter.
To proceed, let us assume that any random variable
(5.14)
u() =
1
X
j =0
uj j ()
where
(5.15-a)
(5.15-b)
(5.15-c)
o 1
E [ j ] = 0
E [ j () k ()] = 0
j>0
j 6= k
3 Also referred to as Wiener chaos from the name of the mathematician who derived it rst.
4 Eq.(5.14) has been preferred to a more detailed notation such as
the sake of simplicity.
u() =
j=0
uj j fk ()gM
k=1
for
80
the coordinates
rst term
1
X
U () =
(5.16)
Uo
U j j ()
j =0
dierent
0
expansion (5.11). The latter, denoted by U , is that obtained by a perturbation approach
(see Chapter 3, Section 2).
By denoting
(5.17)
i=0
1
X
K i i ( )
U j j ()
j =0
F =0
For computational purposes, the series involved in (5.17) are truncated after a nite
M;P =
(5.18)
M PX1
X
i=0 j =0
K i U j i() j () F
; F ; P)
L2(
minimizing
HP ,
which yields :
(5.19)
E [M;P k ] = 0 k = 0 ; ::: P
cijk = E [i j k ]
F k = E [ k F ]
(5.20)
(5.21)
Note that
(5.22)
M PX1
X
i=0 j =0
cijk K i U j = F k k = 0 ; ::: P
(5.23)
K jk =
M
X
i=0
cijk K i
81
PX1
(5.24)
j =0
K jk U j = F k k = 0 ; ::: P 1
K jk a matrix of size
N N . The P dierent equations can be cast in a linear system of size NP NP as
follows :
2
6
6
6
4
(5.25)
U j is a N
K oo : : :
K 1o : : :
.
..
KP
1;o
:::
K o;P
K 1;P
KP
3 2
Uo
U1
1
6
1 7
7 6
7 6
5 4
1;P 1
.
..
UP
7
7
7
5
6
=6
6
4
Fo
F1
7
7
7
5
.
..
FP
KU =F
.
..
U () =
(5.27)
PX1
j =0
U j j ( )
of
HP
is usually 10-35
The
mean
The
E [U ]
namely U o , since E [ j ( )] = 0 for j > 0.
nodal displacement vector
Cov [U ; U ] =
(5.28)
the coecients
E [ 2 ]
Appendix A.1).
PX1
i=1
E 2i
is :
U i U Ti
i 's (See
82
Ui
Eq.(5.27). In the case when this equation is limited to quadratic terms (second
order polynomial chaos), a closed-form expression for the
of
characteristic function
U has been given by Ghanem (1999a ), which can be then numerically Fourier-
Ghanem and Spanos (1991 ). However no such application could be found in the
literature.
It seems possible to couple the general reliability tools developed in Chapter 4,
Section 2 with SSFEM. Let us consider for instance a limit state function of the
following form :
g (U ()) = u uio
(5.29)
uio is a nodal displacement under consideration and u is a prescribed threshold. Substituting the io -th component of the vectorial equation (5.27) in (5.29)
yields the following analytical polynomial expression of the limit state function :
where
(5.30)
g (U ()) = u
PX1
j =0
This limit state function is already cast in the standard normal space due to the
denition of the polynomials
j (fk ()g1
k=1 ). Moreover, its gradient with respect
to the basic random variables can easily be obtained in closed-form as well. Determining the design point and associated probability of failure should thus be
straightforward.
Of course, this approach requires having solved (5.25) beforehand and it is probably
not ecient when a single reliability problem is to be solved. In contrast, it might
be interesting when the probability density function of a response quantity is to be
determined by sensitivity analysis after repeated FORM analyses (see Chapter 4,
Section 4), or when
sponse. This approach must eventually be compared to the classical nite element
reliability approach developed in chapter 4 in terms of accuracy and eciency.
3. Computational aspects
83
3 Computational aspects
3.1 Introduction
As it can be seen in Eq.(5.25), the size of the linear system resulting from the SSFEM
approach increases rapidly with the series cut-o number
P . Whenever
classical direct
methods are used to solve the system, the computational time may blow up rapidly.
This is the reason why early applications of SSFEM were limited to a small number of
degrees of freedom
of SSFEM addressed (Ghanem and Kruger, 1996; Pellissetti and Ghanem, 2000).The
main results are reported in this section.
cijk
K i and the coecients cijk . Storing K as these building blocks K i along with
Kruger (1996) took the example of a 4-term KL expansion. Using a second (resp. third)
order polynomial chaos, the proposed method requires 11 times (resp. 33 times) less
memory compared to the classical global storage. It turns out that a large number of
coecients
cijk are zero (see the tables in Ghanem and Spanos (1991b , chap. 3)).
It is recalled that
'i (x). Since the mean of these uctuations is zero, and if they are bounded
the storage.
conjugate
84
time they are operated on. Since only a small number of coecient
state-of-the-art review) have been proposed. They essentially replace the original system
K U =F
by :
M 1K U =M 1F
(5.31)
K^ oo 0 : : :
6 0
K^ 11 0
6
0
0
M = diag fK^ jj g 6
4
.
..
:::
K^ P
3
7
7
7
5
.
..
where
inc
K^ jj = cojj Linc
Ko U Ko
1;P 1
N = 264 degrees
of freedom and
P = 5; 15,
the authors showed that the proposed preconditioner allows to divide the number of
iterations by 12-15 compared to the Jacobi preconditioner. Moreover, the former leads
to a number of iterations independent of the coecient of variation of the input eld
hierarchical
i.e. P ) does not change the lower-order basis functions. This leads
to the following solution strategy. Suppose the linear system (5.26) is partitioned as
5 The larger the coecient of variation, the larger the condition number of
4. Extensions of SSFEM
follows :
(5.33)
where
()l
and
85
K ll K lh U l = F l
K hl K hh
Uh
Fh
It is then possible to enhance successively the solution (5.27) starting from a lower-order
solution and using (5.35) by adding one basis polynomial at each time. As the lowerorder coecients are not modied along the procedure, this could lead to a successive
built-up of error.
On an example application, Ghanem and Kruger (1996) found no signicant discrepancy
between the results obtained by this procedure and those obtained by a direct higherorder resolution. However, there is no proof or evidence that this is a general result. It
is likely that the accuracy of the results obtained by the hierarchical approach decays
when the coecient of variation of the input eld increases. A more exhaustive study
should be carried out to assess the validity of this approach.
4 Extensions of SSFEM
4.1 Lognormal input random eld
The use of Gaussian random elds is quite common in the context of probabilistic mechanics. However these elds are not well suited to modeling material properties (Young's
modulus, yield stress, etc.) which are by their nature
positive
coecients of variation, realizations of the eld could include negative outcomes that are
physically meaningless. In contrast, the lognormal eld appears attractive in this sense.
A lognormal eld can be dened by a transformation of a Gaussian eld
l(x) = eg(x)
(5.36)
g (x) as :
SSFEM approach, Ghanem (1999 ) proposed to expand them into the polynomial chaos.
Due to the particular form of (5.36), this leads to closed-form expressions.
86
4.1.1
l = eg +g
(5.37)
where
(5.38)
1
X
i=0
li i ( )
i ( ) is the i-th Hermite polynomial in this case. Due to the orthogonality properties of the i 's, the coecients li can be obtained as :
where
li =
(5.39)
li =
(5.40)
1
E [ i ( + g )]
exp [g + g2 ]
2
E [ i ]
2
gi
i!
of any lognormal random variable into the (one-dimensional) polynomial chaos reduces
to :
l = l
(5.41)
4.1.2
1 i
X
g
i!
i=0
i ( )
(5.42)
M
X
i=1
g (x) :
(5.43)
l(x) =
1
X
i=0
li (x) i ()
4. Extensions of SSFEM
87
To use SSFEM in conjunction with a lognormal input random eld is now straightforward : the procedure described in Section 2 applies, where Eq.(5.4) is replaced by
Eq.(5.38). The stochastic equilibrium equation (see Eq.(5.7)) now writes :
1
X
(5.44)
i=0
K i i() U () = F
by :
dijk = E [ i j k ]
(5.45)
The polynomial chaos expansion of the input random eld introduces a new approximation in SSFEM, which probably decreases the accuracy of the method. This accuracy
has not been stated by Ghanem and his co-workers. Whether a fair accuracy could be
obtained with a manageable number of terms in the series expansion is of crucial impor-
bc
are provided in Ghanem (1999 , ). Regarding reliability problems, the accuracy in the
tails of PDFs is probably also aected by the use of the polynomial chaos expansion of
the input random eld.
f1 ; ::: M g for the rst one, fM +1 ; ::: M 0 g for the
second, etc. All these variables are then merged in a single list, the size of which determines the dimension of the polynomial chaos expansion of the response. This technique
increases
dramatically the size of the polynomial chaos basis (see for instance table 5.2, page 96),
which basically controls the computation time.
6 We suppose here the statistical independence of these elds
88
The SSFEM formalism consists in expanding the response process over a basis of
L2( ; F ; P ), namely the polynomial chaos. If the basis functions i () were Dirac
delta functions (
i ), where i
, a collocation-like
procedure along the random dimension would be obtained. Thus the response process is
now considered as the
a nite set of
the denition of
Let
depending on the matrix storage and solving scheme, there should exist a threshold level
for which one procedure (SSFEM or Monte Carlo simulation) becomes more ecient than
the other.
4.3.2
The
hybrid
(5.47)
U ()
P
X1
j =0
U j j () +
Q
X1
j =0
j () = (
j ) results in :
U j j ()
The above expansion is substituted for in the equilibrium equation, and the obtained
residual is made orthogonal both to the
89
Concluding remarks
Details of the hybrid method can be found in Ghanem (1998 ). However no convincing
application of these ideas has been published so far. Moreover, the delta functions do
not form a numerable set, and their use as a basis of
L2( ; F ; P ) or a subspace of it
is questionable. The decoupling assumption which leads to the iterative procedure mentioned above is not really argued. There is globally a lack of mathematical justication
of the method. Further theoretical research as well as applications are needed to assess
the validity of this approach.
ab
Early applications (Spanos and Ghanem (1989); Ghanem and Spanos (1991 , ))
dealt with one- and two dimensional linear elastic structures : a cantilever beam
with Gaussian exural rigidity
EI
its free end, a square plate clamped along one edge and subjected to a uniform inplane tension along the opposite edge with Gaussian Young's modulus, a clamped
curve plate for which the KL expansion had to be computed numerically. Coecients of variation of the response as well as PDF's were determined and compared
to those obtained by Monte Carlo simulation. The number of nite elements in
these examples was limited to 16. The maximal accuracy adopted in these examples was a 3rd order - 4-dimensional polynomial chaos. For medium COV of the
input (say 15-30%), only these most accurate results compare fairly well with the
Monte Carlo simulation results.
Ghanem and Brzkala (1996) addressed the problem of a two-layer soil mass with
deterministic properties in each layer and Gaussian random interface, subjected
to a constant pressure on part of its free surface. In this case, the random eld
representing the Young's modulus of the material is not Gaussian due to its nonlinear relationship with the Gaussian eld dening the interface. Thus the stiness
matrix had to be expanded over the polynomial chaos as in Eq.(5.44).
dia was addressed by Ghanem (1998 ). The permeability coecients as well as the
diusion coecient are modeled as Gaussian random elds and discretized using
90
Karhunen-Love expansion. The eective head and the contaminant concentration are expanded into the polynomial chaos. The numerical results include the
coecients in the expansion of the concentration as well as the variance of the
latter over the domain. No comparison with other approaches (
e.g.
Monte Carlo
simulation) is given.
The problem of heat conduction was addressed by Ghanem (1999 ). In this case,
both the conductivity and the heat capacity are modeled as Gaussian or lognormal
random elds. As an example, a one-dimensional domain of unit length is subjected
to a constant ux at one end and perfectly insulated at the other one. The initial
temperature of the domain is uniform. It is divided in 10 elements. The COV
of the input is up to 40%. The results are presented in terms of the coecients
of the polynomial chaos expansion. Neither post-processing of these results nor
comparison with Monte Carlo simulation is provided in this paper. Due to the
relatively small number of terms included in the expansions
(M = 2
3),
it is
easy task (see Ghanem (1999 ) for some theoretical developments and references
on this topic). Moreover, it is not clear how to enforce the negativeness of the yield
criterion when the stress is now a random quantity.
Thus the authors simplied the problem introducing two
bounding solids,
whose
mechanical properties allow to bound the stresses. The plastic ow rule is thus
applied with
deterministic
sumptions are questionable, this is the only example of real non-linear application
of SSFEM, which shows that a lot of work remains in this matter.
intrinsic.
91
Formally, Eq.(5.27) can be interpreted as a polynomial response surface for the displacement eld, dened by means of the basic random variables
usual response surface methods such as those described in Chapter 4, Section 5, SSFEM
allows to dene it at any order in a consistent framework.
Note that in all applications found in the literature, only the Karhunen-Love expansion
has been used to discretize the input Gaussian random eld. The use of other schemes
such as OSE or EOLE would however be possible and in some case more practical than
KL (for instance when other correlation structures than that with exponential decay are
dealt with).
The approximate solution Eq.(5.27) is obtained in the context of Galerkin minimization
of residuals. General convergence properties to the exact solution are associated with
this procedure : when the number of terms in the series tends to innity, SSFEM tends
to be exact.
However the following limitations of the method have to be recognized :
it is practically limited to
The amount of computation required for a given problem is much greater than that
of the equivalent deterministic problem. Typically 15-35 coecients are needed to
characterize each nodal displacement. As a consequence a huge amount of output
data is available. The question of whether this data is
really
error estimator
the method has been carried out, except some comparisons with Monte Carlo
simulations presented in early papers by Ghanem and Spanos.
Although it is claimed in dierent papers quoted above that the reliability analysis
is a straightforward post-processing of SSFEM, no application could be found in
the literature. The application of SSFEM to reliability analysis remains broadly
an open problem. Important issues such as the accuracy of SSFEM in representing
the tails of the PDFs of response quantities have to be addressed for this purpose.
When lognormal random elds are used, another accuracy issue comes up. Even
for a single variable, only an innite number of terms in the expansion reproduces
the lognormal characteristic. This means that the input eld dened by using only
a few terms in the polynomial chaos expansion (Eq.(5.38)) can be far from the
actual lognormal eld.
92
As a conclusion, it is noted that SSFEM is a quite new approach. Although limited for
the time being, it deserves further investigation and comparisons with other approaches
to assess its eciency.
Appendix
93
Appendix
A.1 Polynomial chaos expansion
A.1.1 Denition
The polynomial chaos is a particular basis of the space of random variables
L2( ; F ; P )
dn e 2 x2 1 2
hn (x) = ( 1)n
e2x
dxn
(5.48)
Multidimensional Hermite polynomials can be dened as products of Hermite polynomials of independent standard normal variables. To further specify their construction, let
us consider the following integer sequences :
= f1 ; ::: pg j 0
i = fi1 ; ::: ipg ij > 0
(5.50)
(5.51)
i; () =
(5.52)
p
Y
k=1
(i ; ) is :
hk (ik ())
Pp
k=1 k
f i;() j
= pg and by p the space they span. p is a subspace of
; F ; P ), usually called homogeneous chaos of order p. The subspaces p are or2
thogonal to each other in L ( ; F ; P ). This is easily proven by the fact that they are
mials
L2(
1
M
(5.53)
k=0
= L2 ( ; F ; P )
94
where
u() = uo
(5.54)
In this expression
o+
1
X
i1 =1
uo ; ui1 ; ui1 i2
u()
1 X
1
X
i1 =1 i2 =1
and second order homogeneous chaoses respectively. The lower order homogeneous chaos
have the following closed-form expression :
=
1 (i ) =
2 (i1 ; i2 ) =
3 (i1 ; i2 ; i3 ) =
o
(5.55-a)
(5.55-b)
(5.55-c)
(5.55-d)
Remark
1
i
i1 i2 i1 i2
i1 i2 i3 i1 i2 i3
i2 i3 i1
i3 i1 i2
The polynomial chaos can be related to the (non orthogonal) basis associated
with the Neumann series expansion, see Eq.(5.13). For this purpose, let us introduce the
orthogonal projection p
relationship holds
of
L2( ; F ; P ) onto
p.
(5.56)
nite number M
are constructed by
are for instance selected from the Karhunen-Love expansion of the input random eld.
by
p (1 ;
associated :
(5.57)
M
Y
i=1
hi (i )
p (1 ;
M +p 1
p
7 This relationship and other mathematical properties of the polynomial chaos can be found in
Ghanem (1999 ).
Appendix
95
M = 4)
p = 4). As an example,
j
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
p
j -th basis polynomial j
p=0
1
p=1
1
2
2
1 1
p=2
1 2
22 1
13 31
p=3
2 (12 1)
1 (22 1)
23 32
14 612 + 3
2 (13 31 )
p=4
(12 1)(22 1)
1 (23 32 )
24 622 + 3
p,
is
is
given by :
P=
(5.58)
p
X
M
k=0
+k
k
and
p.
It is seen that
u ( which was a single number in the deterministic nite element method) is now
represented by P coecients, it is easily seen that SSFEM will require a large amount of
quantity
computation. This may be worthwhile, considering that the whole probabilistic structure
of
coecients.
input random elds. In the original SSFEM, the Karhunen-Love expansion (see chapter 2, Section 5.2) is used under the assumption that the input eld is Gaussian. The
choice of
is thus directly related to the accuracy requested in this random eld dis-
96
P (M
p=1
p=2
p=3
p=4
10
15
15
35
70
28
83
210
captured in describing the solution process. Typical values used in the applications are
M = 4 and p = 2; 3.
g (x) :
g^(x ; ) = g (x) +
(5.59)
M
X
i=1
gi (x) i()
where
M
1X
1
lo (x) = (x) = exp [g (x) +
gi (x)2 ] = exp [g (x) + g^2 (x)]
2 i=1
2
l
g^ (x) is the standard deviation of g^(x). The other ones simplify after some algebra
to :
(5.63)
li (x) = (x)
l
E [ i ( + g (x))]
E [ 2i ]
Appendix
97
i (),
M
Y
gj (x)j
E [ i ( + g (x))] j =1
= M
E [ 2i ]
Y
(5.64)
j =1
Finally, letting
tend to
j !
be written as :
M
Y
(5.65)
l(x) = (x) +
l
1
X
i=1
gj (x)j
X j =1
M
Y
j =1
j !
()
Chapter 6
Conclusions
1 Summary of the study
This report has presented several techniques using the nite element method coupled
with probabilistic approaches. Methods for discretizing random elds, obtaining second moment statistics of the response, probabilities of failure, or approximations of the
stochastic response process itself have been reviewed. In each case, advantages and limitations have been analyzed and examples of application taken from the literature have
been reported.
So far, these examples deal with simple geometries (beams, square plates, sometimes
plates with a hole) and few elements (up to one hundred). Thus the random eld discretization obtained directly or indirectly from the nite element mesh involves a manageable number of variables. However, some work remains on the topic of treating in a
really independent fashion the random eld- and nite element meshes (both of them
being for instance generated automatically with respect to their respective criteria), and
connect them properly.
Perturbation-based approaches were presented in Chapter 3. From a practical point of
deviation). They require gradient operators at the element level in the nite element
code. For strongly non-linear limit state functions, they are expected to be accurate
only with small coecients of variation of the input variables. This could be a limitation
when geomaterials are involved. The CPU time becomes very large when the number of
random variables is medium (say 20-50). The second order approach may be intractable
in this case.
The nite element reliability approach (see Chapter 4) is based on the
coupling of nite
element calculations and a reliability algorithm determining the design point. It allows to
Chapter 6. Conclusions
100
compute the probability of failure of a system with respect to a given limit state function.
Due to this coupled formulation, it is possible to use state-of-the-art nite element codes
by linking them to the reliability program. It has been applied to general non-linear
problems including plasticity, plane stress plasticity and damage. It is applicable to
industrial problems in its current state of development. Current research on this topic
is related to time- and space-variant reliability.
The spectral stochastic nite element method has been applied to linear problems, and it
is not applicable to general non-linear problems yet. However, it is a rather new approach
and deserves further exploration. The main idea of obtaining an approximation of the
stochastic response process itself is denitely attractive due to the wide spectrum of
byproducts it can yield. The SSFEM method is computationally demanding but, on
the other hand, gives a full characterization of the output quantities. Whether this
information is really needed for practical applications is an open question. So is also the
question of the eciency and accuracy of SSFEM in the context of reliability analysis.
perturbation method and Monte Carlo simulation for second moment analysis,
direct coupling between FORM analysis and a deterministic code for reliability
problems.
The implementation issues and comparison results are presented in Part II of the present
report.
Part II
Comparisons of
Stochastic Finite Element Methods
with
Matlab
Chapter 1
Introduction
1 Aim of the present study
Part I of the present report reviewed methods coupling nite element analysis with a
probabilistic description of the input parameters. Emphasis has been put on taking into
account the spatial variability of material properties. This has been done by introducing
the concept of
Second moment approaches (including the perturbation and the weighted integral methods) have been reviewed as well as the so-called nite element reliability methods. Finally, the spectral stochastic nite element method (SSFEM) has been presented, which
is claimed to provide after post-processing second moment as well as reliability results.
As already stated in Part I, there has been little comparison of SSFEM with the other
approaches, at least no comparison with the perturbation method in the context of
second moment analysis, and no reliability study at all. The current part of this report
aims at making these comparisons and thus evaluating the eciency and accuracy of
SSFEM with respect to more classical approaches.
As already discussed in Part I, Chapter 5, SSFEM is only well established for linear
problems so far. Thus elastic two-dimensional mechanical problems have been chosen for
the present study. The conclusions of the study should be understood only in this context.
It is reminded that both the perturbation method and the nite element reliability
methods can be and have already been applied to general non-linear problems (including
large strains, plasticity) as well as dynamics. These approaches have a much larger scope
than SSFEM, at least in its present stage of development. However, in the case when
Chapter 1. Introduction
104
2 Object-oriented implementation in
Matlab
In order to carry out the comparisons mentioned above, numerical tools had to be
implemented. The Matlab environment was chosen for this purpose. The rst reason
is the ease of implementation due to the numerous toolboxes for numerical analysis
provided by Matlab . The second reason is the ability of developing software within
the
pack information into some kinds of objects. Having adopted this way of programming,
it should not be a hard task to transfer the Matlab code into a true object-oriented
3 Outline
The second part of this report is divided into four chapters. The rst two chapters are
devoted to implementation issues, the last two chapters to the comparisons mentioned
above.
Chapter 2 presents a new
toolbox is later used by the dierent programs required by the present study. It practically implements the spectral discretization schemes discussed in Part I, Chapter 2,
Sections 5-6.
Chapter 3 presents the implementation of the SSFEM method. A detailed description of
the implementation of the polynomial chaos expansion (see Part I, Chapter 5) is given.
Post-processing techniques to get second-moment and reliability results are also detailed.
Chapter 4 is devoted to second-moment approaches. The formulation of the perturbation method is particularized to the situation when the randomness in the system is
limited to a random eld describing the Young's modulus of the material. The problem
of simulating random elds representing material properties is then addressed. Finally
the various methods are compared on the example of computing the settlement of a
foundation over an elastic soil mass with spatially varying Young's modulus.
Chapter 5 is devoted to reliability analysis. The post-processing of SSFEM by FORM
and importance sampling is compared to a direct coupling between FORM and a deter1 The
desired.
Matlab
3. Outline
105
ministic nite element code. The serviceability of a foundation over an elastic soil mass
with spatially varying Young's modulus is investigated.
Chapter 2
Implementation of random eld
discretization schemes
1 Introduction
Taking into account material spatial variability in nite element analysis requires the
introduction of random elds and the implementation of
discretization schemes
such
Linear Estimation (EOLE), and Orthogonal Series Expansion (OSE). In order to get a
versatile tool that can be used by itself, an object-oriented implementation in Matlab
is aimed at. All the input data dening the eld, as well as all the quantities required
to evaluate realizations are gathered in a
The three proposed discretization schemes are basically implemented for Gaussian random elds. As an extension, lognormal elds are dealt with by exponentiation.
(2.1)
H^ (x; ) = +
M
X
i=1
Hi (x) i()
108
fi() ; i
standard
normal
variables
and
related to the autocorrelation function of the eld in case of EOLE (see Part I,
Eq.(2.65)),
fhi(x)g1i=0
e.g.
Legendre
e.g.
RFinput.Type
: its value is
RFinput.Mean
RFinput.Stdv
RFinput.CorrLength
scalar
'Gaussian'
in this case.
.
.
elds.
(2.2)
and
'exp2'
(x ; x0 ) =
8
>
<exp(
>
:exp(
`x
`y
8
>
<exp(
(2.3)
x x0 2
) ) for 1D elds
0
`
(x ; x ) =
x x0 2 y y 0 2
>
:exp( (
) (
) ) for 2D elds
`x
`y
1 In Matlab as in C, components of a
structure type
elds
operator .. In the sequel, the word entry is used instead of eld in order to avoid any confusion
with the random eld under consideration.
3. Discretization procedure
109
RFinput.OrderExp
RFinput.NbPts
involves the denition of a grid. This entry contains the number of points along
each direction, dening a uniform grid over the domain. This is a scalar in case of
1D elds and an array of length 2 in case of 2D elds.
"
^l(x ; ) = exp +
(2.4)
M
X
i=1
Hi (x) i ()
In the context of SSFEM, the latter equation is expanded over the polynomial chaos
basis (Part I, Chapter 5, Section 4.1) as :
^l(x ; ) = +
(2.5)
P
X
i=1
li (x) i ()
The input parameters for such elds are the same as those for a Gaussian eld except
the following entries :
RFinput.Type
: its value is
'Lognormal'
in this case.
.
RFinput.LNMean
RFinput.LNStdv
.
l
3 Discretization procedure
From the random eld input and the geometry of the mechanical system (dened as an
array containing the mesh nodal coordinates,
RFinput
110
is determined from the array of nodal coordomain of discretization of the random eld, and is
RF
RF.Domain.
As an alternative, this entry can be input in the form of the following list :
RF.Domain =
fxmin ; xmax g for one-dimensional elds (resp. RF.Domain = f[xmin ; ymin] ; [xmax ; ymax]g
for two-dimensional elds. This option is useful when the random eld toolbox is used
by itself to numerically compare dierent discretization schemes.
H^ (x; ) = +
(2.6)
where
M
X
i=1
x 2
RF
(2.7)
RF
(x ; x0 ) 'i(x0 ) d
x0 = i 'i (x)
3.2.1
One-dimensional case
Suppose
(2.8)
where
jx x0 j
`
1991 ) :
for
i odd, i 1 :
(2.9-a)
(2.9-b)
i =
2`
1 + !i2 `2
i =
1
sin 2!i a
a+
2!i
3. Discretization procedure
111
!i is the solution of :
1
(2.10)
!i tan !i a = 0
`
where
for
i even, i 2 :
i =
(2.11-a)
in the range
2`
1 + !i2 `2
(2.11-b)
[(i
i =
1) ; (i
a
1
) ]
2 a
1
sin 2!i a
2!i
!i is the solution of :
1
1
(2.12)
tan !i a + !i = 0
in the range
[(i
) ;i ]
`
2 a a
All coecients fi ; !i g are computed for i = 1 ; ::: M and stored as additional entries
where
of
RF.
3.2.2
Two-dimensional case
Following Ghanem and Spanos (1991 ), the solution of the two-dimensional eigenvalue
problem is simply obtained by products of one-dimensional solutions,
i = 1i1D 1i2D
'(x) '(x ; y ) = 'i1 (x) 'i2 (y )
(2.13)
(2.14)
where superscript
e.g. :
In implementation, the products of the 1D eigenvalues are computed and sorted in descending order, and the
subscripts
3.2.3
If
RF
(2.15)
e.g.
RF = [xmin ; xmax ], a shift parameter is computed :
T = xmin +2 xmax
RF =
RF
x
min xmax xmax xmin
T=
;
2
2
(2.17)
H^ (x; ) = +
M p
X
i=1
i 'i (x
T ) i() ; x 2
RF
112
uniform
RF .
The
number of points dening the grid (in each direction for 2D problems) is specied in the
RFinput.NbPts. First the coordinates of all the nodes of this grid are computed,
say fx1 ; ::: xN g. The array containing these coordinates is stored in the entry RF.COORD.
entry
H^ (x; ) = +
(2.18)
where
M
X
i()
p iT C x;xi
i
i=1
C x;xi = f(x xi ) ; i = 1 ; ::: M g, and (i ; i) is the solution of the eigenvalue
problem :
C i = i i
(2.19)
RF.COORD.
greatest eigenvalues
i and corresponding eigenvectors i are then computed and stored as additional entries
of
RF.
It is noted that the full eigenvalue problem does not have to be solved, if an
General formulation
The discretized random eld in this case reads (see Eq.(2.54) in Part I) :
H^ (x; ) = +
(2.21)
where
i=1
hi (x) i ()
is dened by :
(2.22)
M
X
(k ; l) = E [kl] =
RF
RF
3. Discretization procedure
113
=
(2.23)
where
as follows :
= 1=2
can be
1=2
is a diagonal matrix whose terms are the square roots of the diagonal terms
H^ (x; ) = +
(2.25)
3.4.2
M p
X
i=1
i
(M
X
k=1
k i hk (x) i ()
fhn(x)g1n=1 is
equations :
(2.26-a)
(2.26-b)
P0 (x) = 1 ; P1 (x) = x
1
[(2 n + 1) x Pn (x) n Pn 1 (x)]
Pn+1 (x) =
n+1
n1
Pn ( 1) = (8 1)n
>
<0
Pn (0) =
n!
>
: n=2 n
2
2 !
Pn (1) = 1
(2.27-a)
(2.27-b)
(2.27-c)
n odd
if n even
if
Z 1
(2.28)
Pm (x) Pn (x) dx =
8
<0
2
:
2m+1
n 6= m
if n = m
if
114
hn (x) =
(2.29)
2n 1
x T
Pn 1 (
)
2a
a
n = 1 ; 2 ; :::
where :
xmin + xmax
2
xmax xmin
a =
2
T =
(2.30-a)
(2.30-b)
are the
(2.31)
2p
(k ; l) = 2a (2 k
1)(2 l
1)
Z xmax Z xmax
xmin
xmin
x0 T
x T
) Pl 1(
) dx dx0
( x ; x0 ) Pk 1(
a
a
x T
a
x0 T
v =
a
u =
(2.32-a)
(2.32-b)
Eq.(2.31) reduces to
(2.33)
Z 1Z 1
a 2 p
(k ; l) = 2 (2 k 1)(2 l 1)
(a u ; a v )Pk 1(u) Pl 1(v ) du dv
1 1
2p
(k ; l) = a 2 (2 k
where
1)(2 l
1)
NPG
X NPG
X
i=1 j =1
f(wi ; Xi) ; i = 1 ; ::: NPGg are the integration weights and points, respectively.
NPG should be
greater than the number of random variables M . In implementation, NPG = 16 was used
in the computation. To eciently evaluate Eq.(2.34), two matrices are rst computed :
Matrix
(2.35)
of size
aXj )
4. Visualization tools
Matrix
LP
of size
115
M NPG, containing :
LP(k ; l) = Pk 1(Xl )
(2.36)
After
performed using a Matlab built-in procedure. The correlation matrix, the eigenvalues
and eigenvectors are stored as additional entries of the object
RFinput.LNStdv,
(stored in
RFinput.LNMean
of the
ln(1 + 2 =2 )
1 2
= ln
2
=
(2.37)
(2.38)
RF.Mean
and
RF.Stdv
respectively.
Then the underlying Gaussian eld is discretized using one of the three methods described in the preceding sections.
4 Visualization tools
When dealing with Gaussian random elds, the accuracy of the discretization can be
evaluated by means of the following
point estimator :
h
(2.39)
116
KL method
#
" 1
X p
1
"rr(x) = 2 Var
i 'i (x) i()
M +1
(2.40)
1
X
M +1
= 1
EOLE method
M
X
i=1
i '2i (x)
M
X
1
i=1 i
"rr(x) = 1
(2.41)
OSE method
i '2i (x)
iT C x;xi
2
fig1i=1 in the KL
and EOLE expansions, it can be seen from Eqs.(2.40)-(2.41) that the variance of
the error is simply the dierence between the variances of
to :
Var H^ (x)
2
"rr(x) = 1
(2.42)
Such a result does not hold when dealing with correlated variables. In the case
of OSE method, the expression for error estimator (2.39) requires a little more
H^ (x ; ) =
(2.43)
M
X
i=1
hi (x)i ()
i () =
(2.44)
H (y ; ) hi (y ) dy
(2.45)
Var H (x)
H^ (x) = E
2
H (x) H^ (x)
h
= E H 2 (x) + E H^ 2 (x)
(2.46)
E H^ 2 (x) = 2
M X
M
X
i=1 j =1
2 E H (x) H^ (x)
hi (x) hj (x) C (i ; j )
5. Conclusion
where
117
"
(2.47)
M
X
hi (x)
i=1
H (y ) hi(y ) dy
By regrouping the deterministic terms outside the expectation operation, one gets :
E H (x) H^ (x) =
(2.48)
M
X
i=1
= 2
hi (x)
M
X
i=1
hi (x)
hi (y )E [h(x)h(y )] dy
Z
hi (y ) (x ; y ) dy
E [H 2 (x)] = 2
(2.49)
"rr(x) = 1+
M X
M
X
hi (x) hj (x) C (i ; j ) 2
i=1 j =1
M
X
i=1
hi (x)
hi (y ) (x ; y ) dy
where the integral in the above equation is numerically computed using Gaussian
integration.
RF
(2.50)
" = j
1 j
"rr(x) d
RF
RF
yielding a mean
error smaller than a prescribed tolerance. Figure 2.1 shows an example of the output for
dierent discretization schemes, in case of one-dimensional eld. So does Figure 2.2 in
case of two-dimensional elds and EOLE discretization scheme.
5 Conclusion
This chapter has presented the implementation of a random eld discretization toolbox
within Matlab . Due to the object-oriented programing, this toolbox is easily extensible
e.g.
e.g. with
triangular decay-
ing) or three-dimensional elds. This toolbox has been used to compare the accuracy of
the series expansion discretization methods in Part I, Chapter 2, Section 6.
118
0.25
Order of expansion : 4
Mean Error over the domain :
KL : 0.112
EOLE : 0.117
OSE : 0.132
0.2
KL
EOLE
OSE
rr(x)
0.15
0.1
0.05
10
x
Figure 2.1: Error
- one-
dimensional eld
0.4
0.3
0.2
0.1
50
0
30
40
30
Discretization
20
scheme :EOLE
10 Order of expansion :3
20
10
0
Figure 2.2: Error estimator computed for EOLE expansion - two-dimensional eld
Chapter 3
Implementation of SSFEM
1 Introduction
1.1 Preliminaries
As presented in Part I, Chapter 5, the Spectral Stochastic Finite Element Method (SSFEM) aims at representing the mechanical response of a system (
e.g.
the vector of
nodal displacements) through its coecients over a basis of the space of random variables
u() =
(3.1)
PX1
j =0
uj j fk ()gM
k=1
where
variable
linear problems. Thus the current implementation is limited to linear elastic two-dimensional me-
chanical problems. Furthermore, the material Young's modulus will be the only parameter considered as spatially variable, and consequently modeled as a random eld.
120
random eld discretization are provided in the Matlab workspace. In the context
of an object-oriented implementation, this data is gathered in two objects called
Model
and
RF
PC.
section 3.
uk ,
fukj g describing
Remark
for the sake of clarity. The user can of course choose any other name, provided there is
consistency in the input le specifying the data to be put in the Matlab workspace.
2 SSFEM pre-processing
2.1 Mechanical model
To get started with a nite element analysis, the following data describing the mechanical
model has to be provided in the Matlab workspace.
A ag variable (
e.g. CONEC)
A connectivity array (
2. SSFEM pre-processing
e.g. ELData)
(ELData.type) and its constitutive material ((ELData.mat). Each entry is an integer array of size
to
121
NbElts.
ELData.type(i) = 1
for
i = 1 ; ::: NbElts.
e.g. MATS)
A material structure (
stitutive law. For linear elastic isotropic material, the entries corresponding to
material #i are :
MATS{i}.E
: Young's modulus
MATS{i}.nu
: Poisson's ratio
MATS{i}.initialstress
MATS{i}.bodyforces
forces in
x and y directions.
BC
is thus 0).
A load vector (
Model.
the parameters of the mechanical model as a single variable to subroutines. The input
data described so far is sucient to run a deterministic analysis. For each application,
such an analysis is carried out systematically in order to check the data as well as the
quality of the mesh with respect to the output quantity under consideration.
e.g. RFinput),
in a structure (
coordinates
Section 3.
COORD,
122
3 Polynomial chaos
3.1 Introduction
M -th dimensional p-th order polynomial chaos consists in a set of multidimensional
Hermite polynomials in f1 ; ::: M g, whose degree does not exceed p (see Part I, Chap-
The
non-negative integers
f1 ; ::: M g as follows :
=
(3.2)
where
hq (:)
is the
M
Y
i=1
hi (i ) ;
i 0
M
X
i=1
0 ; ::: P
1g1 .
E 2j
in
L2( ; F ; P ).
cijk = E [i j k ]. These coecients are required when the input random eld is
Gaussian.
dijk = E [ i j k ]. These coecients are required when the input random eld is
lognormal.
1 In the actual implementation, subscript
on array indexing.
3. Polynomial chaos
123
h0 (x) = 1
d hq (x)
(3.3-a)
(3.3-b)
dx
hq (0) =
(3.3-c)
Polynomial
= q hq 1 (x)
hq
8
>
<0
>
:(
is stored as an array of
1)q=2
q!
q odd
if q even
if
2q=2 2q !
cell array feature of Matlab . A cell array Z is basically an array whose elements
Z fig ; i = 1 ; ::: can be of any kind (i.e. not necessarily the same for all of them). In the
the
present example, this feature is mandatory, since the arrays representing the Hermite
polynomials have dierent lengths.
i in the sequence, skip i boxes and put a ball in the next box;
conversely, for each sample of ball positions, each integer equals the number of
empty boxes (possibly 0) between two consecutive balls.
of degree@ = q
is
the number
of the
M +q 1
M +q
=
M 1
q
M = 4 ; q = 2)) :
For a given
124
ball sample
Integer sequence
Polynomial basis
1 0 1 0
h1 (1 ) h1 (3 ) = 1 3
h2 (4 ) = 42
0 0 0 2
2)
Ball sample
Integer sequence
(M = 4 ; q =
Polynomial basis
42
0 0 0 2
3 4
0 0 1 1
32
0 0 2 0
0 1 0 1
2 4
0 1 1 0
2 3
22
0 2 0 0
1 0 0 1
1 4
1 0 1 0
1 3
1 1 0 0
1 2
12
2 0 0 0
M = 4 ; q = 2)
From the current sample, the next one is recursively obtained by shifting the
i.e.
the ball is
already in the rightmost box), then the rightmost ball that can be shifted by one
box to the right is found. This ball is shifted, and all the balls lying to its right
are brought back to its immediate right.
For each ball sample, the corresponding integer sequence
as an array of length
M.
3. Polynomial chaos
125
E [hp( ) hq ( )] = pq p!
(3.4)
where
E [ ] =
(3.5)
where
M
Y
i=1
i !
and 0 otherwise.
3.4.2
When dealing with Gaussian random elds within SSFEM, the coecients
(resp.
cijk =
k ) corre-
l6=i
j0 =
6 i, j0
hj0
and
and
j0
hj0 .
Otherwise,
Y
l6=i
l !
(3.8)
E [ hp ( )hq ( )] =
Z
R
p1
2
, one gets :
hp(x) hq (x) x e
1 x2
2
dx
126
dhp ( )
dhq ( )
d
E [ hp ( )hq ( )] = E
f
hp ( )hq ( )g = E
hq ( ) +
h ( )
d
d
d p
(3.9)
(3.10)
cijk .
fdijk
E [ i j k ]g are required. From Eq.(3.2) and the independence of the standard normal
variables f1 ; ::: M g, it follows that :
dijk =
(3.11)
M
Y
l=1
l . Unfor-
tunately, no simple formula similar to Eq.(3.10) can be derived in this case. Thus the
rather inecient following algorithm is used :
(3.12)
r=0
arl l r
Using the linearity of the expectation operator and the closed form solution for
the moments of the standard normal variable
(3.13)
l +X
l +
l
l , one obtains :
l = 1 ; ::: M
l +X
l +
l
r=0
r even
arl
r!
2r=2
r
2 !
Eq.(3.11).
Due to the symmetry in the subscripts of the coecients subscripts, only those
are associated with
0ijkP
dijk that
4. SSFEM Analysis
127
3.5 Conclusion
This technical section has described the practical implementation of the polynomial
chaos. All the basis polynomials and related coecients are eventually gathered in a
single object called
PC,
Ghanem and Spanos (1991 ), no details about the implementation of SSFEM could be
found in the literature. The method proposed by Ghanem and Spanos used symbolic
calculus, which is much more complicated to implement than the approach proposed in
the present report. It is believed that the present chapter provides new solutions to this
problem.
4 SSFEM Analysis
As any nite element software, the core of the SSFEM program consists in computing
element stiness matrices and nodal forces, assembling element contributions and solving
the obtained linear system. These dierent steps are described in detail in the sequel.
(3.14)
where
e ( )
H (x ; ) B T D o B d
e
is the elasticity matrix computed with unit Young's modulus. Substituting for the truncated series expansion (2.1) into (3.14) leads to computing the following
matrices :
the
k =
e
(3.15)
B T Do B d
e
k =
e
i
Hi (x) B T Do B d
e i = 1 ; ::: M
deterministic
128
f =
e
(3.17)
where
B o d
e +
T
N T b d
e
i.e. :
K =
(3.18-a)
Ki =
(3.18-b)
[
e
[
e
k =
e
kei =
[Z
e
e
[Z
e
global stiness
B T Do B d
e
Hi (x) B T D o B d
e
The Galerkin technique associated with SSFEM then leads to the following system of
equations using the above matrices (see Part I, Eq.(5.22)) :
M PX1
X
(3.19)
i=0 j =0
cijk K i U j = F k k = 0 ; ::: P
Uo
U .
i=0
right hand side of Eq.(3.19) are all zero except the rst one
(3.20)
K jk =
M
X
i=0
+
cijk K i = cojk K
M
X
i=1
cijk K i
i.e.,
K o K
and
F k appearing in the
j ; k = 0 ; ::: P
2
(3.21)
2A
6
6
6
4
K oo : : :
K 1o : : :
.
..
KP
1;o
:::
K o;P
K 1;P
KP
.
..
3 2
1
6
1 7
7 6
1;P 1
76
5 4
Uo
U1
.
..
UP
7
7
7
5
=6
6
6
4
Fo
0
.
..
3
7
7
7
5
3 3 scheme was tried but gave practically the same result for larger computation time, and was
thus abandoned.
4. SSFEM Analysis
129
KU =F
(3.22)
from the
large matrix
is
P , where N
i.e.
twice
xed
degrees of freedom
i.e.
deterministic
uj k = 0
(3.23)
where
uj k
is the
8 k 2 I ; 8 j = 0 ; ::: P
in (3.21).
8 k 2 I ; 8 j = 0 ; ::: P
1,
(j N + k)
Fj k
is set equal to 0.
However the same result can be obtained with greater computational eciency by apply-
130
; K i ; K jk ) are sparse, that is, contain a great number of zeros. Thus the sparse
i.e. K
storage option in Matlab is used in their declaration. In this case, only the non zero
terms of are stored together with the indices of their position.
After having declared these matrices as sparse, the user can perform any operation
without worrying about the storage issues, which greatly simplies the implementation.
The solution of the problem,
built-in solver.
By running examples, it appeared that the most time consuming step in the analysis is
the second level of assembly. The solving step itself usually represented a small fraction
of the total computation time. This can be explained by the fact that an interpreted
code is used for the assembly steps, whereas the solver is a compiled Matlab routine.
A dierent behavior would be expected if the program were to be fully compiled.
5 SSFEM post-processing
The crude output of the SSFEM analysis is a set of nodal displacement coecients
PX1
U () =
(3.24)
U j j fk ()gMk=1
j =0
e ( ) =
(3.25)
PX1
j =0
(3.26)
(x ; ) = B (x)
"
(3.27)
(x ) = +
PX1
j =0
M
X
i=1
i ()
Do B (x)
PX1
j =0
5. SSFEM post-processing
131
E [U ] = U
(3.28)
Cov [U ; U ] =
(3.29)
Uo
PX1
i=1
E 2i
U i U Ti
Similar results can be obtained for strain and stress components by means of Eqs.(3.26)(3.27).
(3.30)
where
g (U ()) = u
(3.31)
P 1
X
j =0
failure. Thereafter, importance sampling around the design point can be performed in
order to get a more accurate value of the probability of failure. This is computationally
cheap since the expression of the limit state function is analytical.
It is emphasized that the limit state function is already dened in terms of standard normal variables, which avoids any probabilistic transformation within the
iHLRF algorithm
(3.32)
0 if k = 0
@
=
Q
@k
k hk 1 (k ) l6=k hl (l ) otherwise
132
uio
(3.33)
where
Fio
fio (u) =
(3.34)
uio
uio . Using
is :
dPf
d
dFio
=
= '( (u))
du
du
du
d
=
du
(3.35)
k r
g (U ( (u)) ; u) k
with respect to
@g ( ; u)
@u
u
fio (u) =
(3.36)
'( (u))
k r g(U ( (u)) ; u) k
To compute the entire PDF of a nodal displacement, a FORM analysis is carried out for
dierent thresholds
u,
6 Conclusions
This chapter has presented the implementation of the Spectral Stochastic Finite Element Method (SSFEM) in the Matlab environment. Object-oriented programming
was aimed at, rst to allow a versatile utilization of the code, second to build a base for
later implementation in a true object-oriented language like C++.
An original implementation of the polynomial chaos basis was proposed, which is believed
to be simpler than the only other implementation found in the literature,
i.e.
that by
Monte Carlo simulation (resp. direct coupling between a deterministic nite element code
and the FORM procedure). This is the goal of Chapter 4 (resp. Chapter 5).
Chapter 4
Second moment analysis
1 Introduction
The aim of this chapter is to compare second moment methods for elastic twodimensional mechanical problems involving spatial variability of material properties.
Precisely, modeling the Young's modulus of the material as a random eld, the mean
value and standard deviation of response quantities such as nodal displacements are
computed. Three dierent methods have been implemented.
The crude
eld, then carrying out a deterministic nite element analysis of the mechanical
problem, and nally, statistically treating the response quantities of interest.
The
perturbation method
rst and second order, the specic formulation associated with random elds being
rst discussed.
These three approaches are applied to compute the statistics of the settlement of a
foundation over an elastic layer, whose heterogeneity is accounted for by modeling its
Young's modulus as a random eld.
134
H (x ; o ) of the
random eld modeling the material Young's modulus, the element stiness matrix is
given by (See Eq.3.14) :
(4.1)
e ( )
o
H (x ; o ) B T (x) Do B (x) d
e
Using the isoparametric formulation, this integral is recast over a reference domain
(4.2)
where
e (
o)
(4.3)
where
NPG
e (
o)
NPG
X
i=1
i and related
135
H (x( i ) ; 0 )
follows :
(4.4)
(4.5)
H^ (x(i ) ; 0 ) = +
M
X
i=1
"
Hi (x(i )) i (o )
H^ (x(i ) ; 0 ) = exp +
M
X
i=1
if H () Gaussian
#
(4.6)
(4.7)
EMC [U ] =
VarMC [U ] =
sim
1 NX
U (i )
Nsim i=1
Nsim
"N
sim
X
i=1
2
Nsim is the number of samples considered, U (i ) is the nodal displacement vector
2
associated with sample i , and U (i ) is the vector containing the square values of the
where
nodal displacements.
In implementation, the sums
each nite element run. The accuracy of the Monte-Carlo simulation is estimated by the
coecient of variation of the empirical mean (4.6), which is given by :
(4.8)
COV MC =
VarMC [U ]
pN EMC [U ]
sim
136
e.g. Liu et al. (1986b ,a ); Ghanem and Spanos (1991a ,b ); Deodatis and Shinozuka
(1991); Anders and Hori (1999)) assess the validity of the proposed approach by comparing the results with those obtained by a Monte Carlo simulation. In practice, these
authors use Gaussian random elds.
When material properties are modeled, the use of Gaussian random elds is questionable. Indeed realizations of Gaussian random variables can be negative valued, whereas
the material properties are positive in nature. The results obtained by Monte Carlo simulation in this case are denitely doubtful due to the possible negative outcomes. Either
these negative outcomes have to be discarded, which introduces a bias in the simulation,
Moreover, the general result, which says that Monte Carlo simulation is asymptotically
exact (the more samples, the best result) does no longer hold : in this case indeed, the
more samples, the more non physical outcomes.
This problem has not received much attention in the literature, the authors of the papers
mentioned above do not even bring up the question. From the remarks above, the following strategy will be adopted in the present study : Monte Carlo simulation of Young's
modulus will only be performed with lognormal distributions.
is limited to spatial variability of the material Young's modulus. Thus the basic variables used in the Taylor series expansion are
= f1 ; ::: M g used in the random eld discretization, and the formulation of the perturbation method presented in a general context in Part I, Chapter 3, Section 2 can be
simplied.
(4.9)
K () K o +
(4.10)
U () U o +
N
X
i=1
N
X
i=1
137
K Ii i + 21
U Ii i + 12
N X
N
X
i=1 j =1
N X
N
X
i=1 j =1
K IIij i j
U IIij i j
where the notation used in the above equations has been introduced in Part I, Chapter 3,
Section 2. Moreover, the load vector
K () =
(4.11)
"
[Z
+
M
X
i=1
Hi (x) i
B T D o B d
K with respect to i is :
[Z
@K
I
K i @ =0 =
Hi (x) B T D o B d
second derivatives
From Eqs.(4.9)- (4.10), the Taylor series expansion of the equilibrium equation is :
(4.13)
Ko +
M
X
i=1
K Ii
M X
M
1X
o
I
U IIij ij = F
U + U i i + 2
i=1 j =1
i=1
M
X
i and ij
U o = Ko 1 F
U Ii = K o 1 K i U = Li U where Li = K o 1 K i
U IIij = U IIji = K o 1 K Ii U Ij + K Ij U Ii = Li U Ii + Lj U Ij
= (Li Lj + Lj Li ) U o
138
that :
E [i] = 0
Cov [i ; j ] = ij
(4.17-a)
(4.17-b)
(Kronecker symbol)
Thus the rst- and second-order estimates of the nodal displacement mean values are :
(4.18-a)
EI [U ] =
Uo
(4.18-b)
EII [U ] =
U o + 12
M
X
i=1
L2i U o
CovI [U ;
(4.19)
U] =
M
X
i=1
U is :
U Ii U Ii T
E [i j k ] = 0
E [i j k l ] = ij kl + ik jl + il jk
(4.20-a)
(4.20-b)
(4.21)
CovII [U ;
U ] = CovI [U ;
M X
M
M X
M
1X
1X
T
II
II
U] + 4
U ii U jj + 2
U IIij U IIij T
i=1 j =1
i=1 j =1
some sense, corresponds to selecting a dierent point around which the expansions are
carried out), the global stiness matrix would write :
(4.22)
K ( ) =
"
[Z
e
exp +
M
X
i=1
Hi (x) i
B T Do B d
(4.23)
K
(4.24)
K IIij
I
i
139
[
@K
=
e Hi (x) B T D o B d
@i =0 e
e
[Z
@2K
=
=
e Hi (x) Hj (x) B T D o B d
= e+ 2 ,
2=
where
is the standard
deviation of the approximate underlying Gaussian eld. This means that the expansion
in this case would not be carried out around the mean value
the results is expected to be worse than that obtained in the previous section. Moreover,
the computation of the second order terms
which makes the whole computation much more time consuming. This type of expansion
will not be used in the numerical applications.
2B
A
E;
Due to the symmetry, half of the structure is modeled by nite elements. Rigorously
speaking, there is no more symmetry in the system when random material properties are
introduced. However, it is believed that this simplication does not signicantly inuence
the results. The parameters selected for the computation are gathered in Table 4.1.
140
A rened mesh was rst used to get the exact maximum displacement under the
foundation (point A in Figure 4.1). Then dierent meshes were tried in order to design
an optimal mesh,
smallest number of elements. The mesh displayed in Figure 4.2-a was eventually chosen.
It contains 99 nodes and 80 elements.
Table 4.1: Settlement of a foundation - Parameters of the deterministic model
Soil layer thickness
Foundation width
Applied pressure
Soil Young's modulus
Soil Poisson's ratio
Mesh width
a -
t
2B
P
E
L
Mesh
30 m
10 m
0.2 MPa
50 MPa
0.3
60 m
b -
Deformed shape
Figure 4.2: Settlement of a foundation - Mesh and deformed shape obtained by a deterministic analysis
For the input parameters given in Table 4.1, the maximum displacement obtained
with the most rened mesh is
in Figure 4.2-a is
uexact
= 5:49
A
( ; ) can
variation E = E =E of E by :
141
and coecient of
ln(1 + E2 )
1 2
= ln E
2
=
(4.26)
(4.27)
E
UA
corresponding to any
UA (E ) = uo
(4.28)
E
E
UA (E ) = eln(uo E )
(4.29)
Thus
the
maximum
LN (ln(uo E ) ; ),
settlement
whose
mean
UA
value
follows
and
lognormal
standard
deviation,
distribution
after
some
UA = uo (1 + E2 )
UA = UA E = uo E (1 + E2 )
(4.30)
(4.31)
These analytical expressions are used in the sequel to assess the validity of the rstand second order- perturbation and SSFEM methods. Results are presented for dierent
coecients of variation
4.2.2
E .
Numerical results
The dierent softwares developed for dealing with random elds are used in this example
by setting the correlation length of the random eld equal to
` = 10000
m, and by
random eld.
Perturbation method
ure 4.3. Considering a good approximation as one giving less than 5% discrepancy
from the exact result, it appears that the (constant) rst-order estimate of the mean is
acceptable only for
value, whatever
142
0.9
1.45
Analytical
First Order
Second Order
1.4
1.35
0.7
1.3
0.6
U / uo
1.25
0.5
1.2
U / uo
Analytical
First Order
Second Order
0.8
0.4
1.15
0.3
1.1
0.2
1.05
0.1
1
0.95
a -
0.1
0.2
0.3
0.4
Coefficient of variation E
0.5
0.6
b-
0.1
0.2
0.3
0.4
Coefficient of variation E
0.5
0.6
Perturbation results
0.9
1.45
Analytical
Order PC =1
Order PC =2
Order PC =3
1.4
1.35
0.7
1.3
0.6
U / uo
1.25
0.5
A
1.2
U / uo
Analytical
Order PC =1
Order PC =2
Order PC =3
0.8
0.4
1.15
0.3
1.1
0.2
1.05
0.1
1
0.95
0.1
0.2
0.3
0.4
Coefficient of variation
0.5
0.6
a -
b-
0.1
0.2
0.3
0.4
Coefficient of variation E
0.5
0.6
SSFEM results
exact Taylor series expansion of the stiness matrix has only constant
E .
E , which is
As a conclusion, it is seen on this example that the rst-order perturbation method can-
143
not be expected to give satisfactory results for medium to large coecients of variation
of the input random eld. In contrast, the second-order approach is much more accurate,
whatever the value of
latter is inexpensive to apply due to the fact that the second-order derivatives
K IIij are
all zero.
SSFEM method
Results regarding the SSFEM approach are given in Figure 4.4 for
dierent orders of expansion. As only one random variable is used to discretize the
underlying random eld, the size of the polynomial chaos of order
p is P = p + 1. The
polynomial chaos is used to expand both the global stiness matrix and the vector of
nodal displacements, as described in Part I, Chapter 5, Section 4.1.
It appears that a good accuracy is obtained with
modulus, whereas
p=2
p=3
Figure 4.4-b that the rst order result is closer to the exact solution than the second
order. No satisfactory explanation to this behavior could be found.
Problem statement
To account for the heterogeneity in the soil, the Young's Modulus is now modeled as a
homogeneous
mean value :
E = 50 MPa,
standard deviation :
[0 ; 0:5],
within
E = E E
varies
LRF =` = 1=2 1=10. The mean of the error variance (2.50) has
been computed in each case for dierent orders of expansion M , the results are reported
element size
LRF
satises
in Table 4.2.
If selecting a tolerance of 10% for the accuracy in the discretization, it appears from
Table 4.2 that
144
Eq.(2.50))
M LRF =` = 1=2 LRF =` = 1=3 LRF =` = 1=4 LRF =` = 1=5 LRF =` = 1=10
1
0.467
0.462
0.461
0.460
0.459
0.236
0.228
0.226
0.225
0.224
0.151
0.143
0.140
0.139
0.138
0.085
0.078
0.076
0.075
0.074
0.048
0.041
0.038
0.037
0.036
of the EOLE grid does not signicantly improve the accuracy, for a given
M.
Thus
As there is no closed form solution to the statistics of the response, a Monte Carlo
simulation is performed. For each coecient of variation
obtained simulation error measured by Eq.(4.8) is increasing from 0.1% to 1.5% when
the COV of the input random eld varies from 0.05 to 0.5)
4.3.2
Numerical results
UA and the standard deviation UA of the maximum settlement UA are
0.7
1.3
Monte Carlo
Order PC =1
Order PC =2
Perturbation 1st
Perturbation 2nd
1.25
1.2
Monte Carlo
Order PC =1
Order PC =2
Perturbation 1st
Perturbation 2nd
0.6
0.5
1.15
U / uo
U / uo
0.4
1.1
0.3
1.05
0.2
1
0.1
0.95
0.9
a -
0.1
0.2
0.3
Coefficient of variation E
0.4
0.5
0.1
0.2
0.3
Coefficient of variation
0.4
0.5
b-
Figure 4.5: Settlement of a foundation - Young's modulus modeled by a lognormal random eld
145
UA and UA are smaller than the values they take in case of homogeneous
Young's modulus (See Figures 4.3-4.4). As in the case of homogeneous Young's modulus,
considered. Both the second-order perturbation and the second-order SSFEM methods
give good accuracy for any COV in the range under consideration.
CPT
4508
11.8
21
19.5
1446
Regarding the Monte Carlo simulation, the computation time comes almost only from
the successive deterministic nite element runs. The second order perturbation method
requires about twice as much time as the rst order.
As far as SSFEM is concerned, there is a huge dierence between the computational
cost of the rst- and the second-order methods. Higher order computations could not
be carried out. This can be explained by the fact that the discretization of the random
eld required
M = 4
e.g., p = 3).
large size of the polynomial chaos basis even for small order (
Completely
dierent computation times would have been observed in the rst example (
of lognormal Young's modulus), where only one random variable was used.
i.e in case
146
It is emphasized that the computed times reported in Table 4.3 are related to the discretization scheme employed (here EOLE) and would be dierent if another discretization scheme was chosen. They are of course related to the Matlab implementation,
5 Conclusions
In this chapter, three second moment methods have been presented in the context of
spatial variability of material properties and applied to a geomechanical problem.
The particular formulation of each method in the present context has been derived. It
appears that the perturbation method is inexpensive to apply up to second order, due
to the fact that the second-order derivatives of the stiness matrix with respect to the
basic random variables are zero.
The accuracy of each method has been investigated for dierent values of the coecient of
variation of the input random eld. The cost of each approach has been nally evaluated.
After compiling all the results, it appears that the second-order perturbation method is
the most attractive for problems involving random elds, because it is inexpensive and
accurate even for large coecients of variation of the input. The SSFEM approach also
gives accurate results when applied at second and higher orders. The computation time
may however blow up when more than 2 or 3 random variables are used to discretize
the random eld.
Chapter 5
Reliability analysis of mechanical
problems involving random spatial
variability
1 Introduction
The aim of this chapter is to compare two dierent approaches for solving reliability
problems based on elastic two-dimensional analyses involving random spatial variability
of material properties.
iHLRF
into account the spatial variability, the deterministic code Femrf described in
Chapter 4, Section 2.2 is used. Details are given in Section 2.
Both approaches are applied to compute the reliability index associated with the maximum settlement of a foundation lying on an elastic soil layer.
A parametric study is carried out using a Gaussian (resp. lognormal) random eld in
Section 3 (resp. Section 4). Comparisons of the two approaches described above are made
by varying successively :
the order of expansion in the random eld discretization as well as that of the
polynomial chaos expansion,
148
Remark
As already stated in Chapter 4, Gaussian random elds are not well suited
to model material properties, due to possible negative outcomes. In the context of reliability analysis, even the computed design point could happen to correspond to non
physical values, for instance when large coecients of variation of the input are considered. However, as Gaussian random elds have been used extensively in the literature
together with the Karhunen-Love expansion, they will be used for comparison purposes
in Section 3.
by the
iHLRF
(o) is used in
in that section is now applied to problems for which the randomness is limited to random
elds describing the material's Young's modulus.
The limit state function under consideration in the present study is :
(5.1)
where
Accordingly, the gradient of the limit state function with respect to the basic random
149
(5.2)
(5.3)
1 ; 0 ; ::: ]
B the matrix yielding the strain components from the nodal displacements ue , and D is the elasticity matrix. In case of
where
@i
with respect to
D = H (x) Do ( +
(5.5)
M
X
i=1
Hi (x) i) Do
@ D @H (x)
=
Do = Hi(x) Do
@i
@i
(5.6)
@U
=
@i
(5.7)
1
[ Z
e
Hi (x) B
Do B d
ue
@U
=
@i
(5.8)
where
Ki
is the
i-th
K 1 Ki U
(5.8) in (5.2) nally gives the following closed-form expression for the gradient of the
limit state function :
(5.9)
150
g (U ()) = u UA ()
(5.10)
where
is a given threshold.
In this section, the random eld modeling the Young's modulus of the soil is supposed
to be Gaussian, and has the following properties :
mean value
E
= 50 MPa,
E =E 2 [0 ; 0:5],
exponential
E ,
E =
the literature make use of this kind of correlation structure together with the
Karhunen-Love discretization scheme, this form is assumed in this section. To get
i.e `x = 10000 m, `y = 30 m.
u.
uo
uo = 5:42 cm.
151
Direct coupling
direct
of the in-
put random eld. Results are reported in Table 5.1 together with the accuracy of the
discretization (estimator
Table 5.1: Inuence of the accuracy in the random eld discretization - Direct coupling
approach
3.2.2
"
direct
0.269
2.694
0.129
2.631
0.082
2.627
0.060
2.627
0.048
2.627
SSFEM +FORM
SSF EM
3.2.3
and
p,
the results
For our choice of parameters, it appears that the direct coupling allows to get 2-digit accuracy in the reliability index
to
" 10%.
direct as soon as M
M , SSF EM
converges to
direct
of the polynomial
(M = 2 ; p = 5).
resulting SSFEM system of equations is 198 21 = 4148.
computation times is
P =21.
This corresponds to
152
direct
M
1
2.694
2.631
2.627
2.627
and
p - SSFEM approach
SSF EM
4.665
3.008
2.741
2.685
2.681
4.510
2.904
10
2.656
15
2.611
21
2.614
4.487
10
2.889
20
2.645
4.480
15
2.885
M = 2 (" = 0:129)
variation of the input is 0.2. The reliability index is computed for dierent thresholds
of maximum settlement
polynomial chaos expansion are used in this case). Results are reported in Table 5.3.
When direct coupling is used, it is observed that the number of iterations required by the
iHLRF
out using 3 terms in the Karhunen-Love expansion; the results are equal to those given
in Table 5.3 with less than 1% discrepancy).
i.e less than 5% discrepancy between SSF EM and direct ) are obtained only for u 20 cm ( 4). For
5, a good accuracy using SSFEM would require higher orders of expansion (p > 5),
which becomes intractable in our example. Moreover, the example was cooked up so
that
M =2
use up to 5-th order polynomial chaos expansions. Usually, the number of random vari-
e.g.
153
order SSFEM method would be practically applicable. From Table 5.3, it is seen that,
on the present example, the
154
Table 5.3: Inuence of the threshold in the limit state function - Direct coupling and
SSFEM results
(cm)
10
12
15
20
30
direct
0.553
1.856
2.631
3.143
3.648
4.139
4.601
# Iterations
10
12
13
SSF EM
0.392
0.504
0.564
0.564
0.552
2.451
1.859
1.821
1.842
1.858
4.509
2.904
2.655
2.610
2.614
6.568
3.787
3.298
3.161
3.126
9.656
4.926
4.065
3.782
3.674
14.803
6.523
5.054
4.535
4.304
25.096
9.093
6.498
5.564
5.118
155
Table 5.4: Inuence of the correlation length of the input random eld (
` = 10
m) -
M
1
"
0.550
direct
3.441
2
0.335
0.232
0.175
3.215
3.181
3.180
SSF EM
6.198
3.978
3.583
3.474
3.443
5.778
3.690
3.327
3.232
3.208
5.671
3.625
3.277
5.646
3.608
0.140
3.179
10
0.071
3.179
When comparing column#2 of Table 5.4 with column#2 of Table 5.1 (which corresponds
to
two-digit accuracy on the reliability index is obtained when direct coupling is used.
SSF EM
and
direct )
156
M = 10 or
more. In contrast, it would not be possible to apply SSFEM with p > 2 when M is more
than 10, which means that the obtained reliability index would probably be inaccurate.
variation of the input random eld. Results are reported in Table 5.5.
When the direct coupling is used, convergence of the
(the
, the more iterations). The values obtained are within 1% of those obtained with
M = 3. It is observed that the reliability index strongly decreases when the variability
higher
value induces a relatively large reliability index (see Section 3.3). For larger
however,
the results are not very good either. Some FORM analyses carried out after SSFEM do
not converge, some others converge to a wrong design point, especially when the order
of the polynomial chaos is large. This may be explained by the fact that the polynomial
response surface associated with SSFEM is undulatory in this case (due to higher order
polynomials) and may have several local design points. As an example, for the case of
E =0.4, it is observed that the convergence to the true reliability index is not monotonic
with increasing p. Thus the result obtained for a given order cannot be a priori positioned
with respect to the true value.
From these examples, it appears that SSFEM coupled with FORM cannot be applied
one-dimensional in
discretization error
two-dimensional and isotropic, with a correla` = 30 m. The coecient of variation of the eld is set equal to 0.2 and the
threshold in the limit state function is u =10 cm. The direct coupling and the SSFEM
157
Table 5.5: Inuence of the coecient of variation of the input random eld - Direct
M = 2)
E
0.1
0.2
0.3
0.4
direct
8.277
4.132
2.759
2.069
SSF EM
30.706
13.769
10.702
9.578
9.043
14.803
6.523
5.054
4.535
4.303
9.257
3.925
2.994
2.666
2.467
6.301
2.455
1.708
y
y
2.045
4.380
1.370
3.062
1.592
1.227
0.5
1.655
0.807
y For these values, the iHLRF algorithm applied after SSFEM has not converged after 30
iterations.
method are applied with dierent orders of expansion
and
Table 5.6
It can be seen from column #2 that the discretization error
"
that obtained for a one-dimensional random eld. For instance, even 50 terms in the
Karhunen-Love expansion do not allow to have
The direct coupling can still be applied up to this order of expansion though. In-
158
M
1
"
0.586
direct
4.212
2
0.442
0.362
0.303
3.271
3.239
3.019
SSF EM
7.647
4.924
4.427
4.281
4.232
5.823
3.748
3.387
3.293
3.269
5.736
3.692
3.343
5.303
3.414
3.098
0.273
2.946
10
0.177
2.876
50
0.059
2.826
deed, as it will be explained in Section 3.7, the computation time for direct coupling
is approximately linear with M. The best result obtained with direct coupling here is
direct = 2:826, which is probably a slight over-estimate of the true reliability index.
M = 4 and p = 3 yielding
SSF EM = 3:098, which means at best 10% accuracy in the reliability index.
(M ; p)
159
index.
Table 5.7: Computer processing time required by direct coupling and SSFEM methods Gaussian random elds
CPT y Direct
Coupling (")
20.6
33.6
CPTy
SSFEM (")
2.3
2.5
3.2
4.0
4.8
3.9
8.0
22.6
58.2
129.0
4.7
3
4
43.7
30.3
296.4
1888.7
8.7
127.4
11.4
2
3
53.8
65.8
The CPT for a deterministic nite element run with constant Young's modulus was 0.57".
From column #2 of Table 5.7, one can see that the CPT required by the direct coupling is
of the element stiness matrices. Each of these matrices requires the evaluation of the
random eld realization at four points (the Gauss points), and each evaluation takes a
time exactly proportional to the order of expansion
gradients computed is also proportional to
M.
In contrast, when using SSFEM, the CPT increases extremely fast with the order of
the polynomial chaos expansion. Thus the method can be eciently applied only when
a small number of terms
the reliability index under consideration is suciently small so that the second order
SSFEM already gives a fair estimate.
160
Introduction
The SSFEM approach allows to get an approximation of the random response of the
structure in terms of polynomials in standard normal variables, see Eq.(3.1). In the
context of reliability analysis, this allows to dene
analytical
the analytical polynomial expression of the limit state function contains information that is lost when the linearization resulting from FORM is used.
the limit state surface obtained from SSFEM could be globally accurate, however
not necessarily around the true design point, which means that applying FORM
could give poor results.
Moreover, since the limit state function is inexpensive to evaluate due to its analytical
expression, simulation methods such as importance sampling become attractive.
3.8.2
Numerical results
order of expansion
(resp.
limit state surface is an hyperplane). For all these cases, it can be seen in Table 5.8 that
importance sampling gives exactly the same results as FORM (the last-digit discrepancy
being explained by the fact that only 10,000 samples are used in the simulation).
161
Table 5.8: Post-processing of the SSFEM results - Comparison between FORM and
importance sampling
M
1
F ORM
SSF
EM
IS
SSF
EM
4.665
4.669
3.008
3.012
2.741
2.738
2.685
2.689
4.510
4.515
2.904
2.891
2.656
2.633
2.611
2.578
2.614
2.580
4.487
4.490
2.889
2.872
2.645
2.608
Signicant discrepancies between the two approaches appear only for higher orders of
polynomial chaos expansion,
of the reliability index, which means that the FORM result is satisfactory in all cases.
From this short study, the following conclusions can be drawn :
for the example under consideration, the limit state surface dened analytically
after the SSFEM analysis is suciently smooth so that the rst-order reliability
method gives good results.
after having determined the design point by FORM, importance sampling allows to
evaluate more accurately the probability of failure at low cost, due to the analytical
denition of the limit state function.
F ORM
SSF
EM
ans
IS
SSF
EM
reliability estimates by SSFEM in the previous sections are due to the truncation
of the polynomial chaos expansions and not due to the FORM approximation.
162
dened analytically and thus easy to evaluate. This allows to compute at low cost the
probability density function of a response quantity, as described in Chapter 3, Section 5.4.
are solved . To improve the eciency, the starting point of each analysis is chosen as the
design point of the previous analysis. This allows convergence of the iHLRF algorithm
within 3 iterations.
50
45
40
35
fU (u)
30
25
20
15
10
5
0
0.12
0.1
0.08
0.06
u
0.04
0.02
Figure 5.1: Probability density function of the maximum displacement obtained by multiple FORM analyses after SSFEM
It can be seen that the obtained PDF has its mode close to
uo = 5:42
cm (which is
the value obtained from a deterministic nite element analysis) and that it looks like
a lognormal distribution, in agreement with the results of Chapter 4, Section 4.2. It
should be emphasized that the far tails of the PDF computed by this method may be
inaccurate, as observed in Section 3.3.
3.10 Conclusions
From the above comprehensive parametric study, the following conclusions can be
drawn :
1 This is done in a matter of seconds on a personal computer.
163
a deterministic nite element code converges to a limit when the discretization error
" tends to zero. This convergence is always obtained by upper values. As soon as
" 10%, the method gives a 2-digit accuracy for , whatever its value (at least in
the range
The CPT for each iteration is increasing linearly with the order of expansion
iHLRF
chaos increases. This means that SSFEM may be applicable in some cases to
solve reliability problems. However, this convergence presents an unstable behavior,
which makes the method unreliable.
2
When
3,
the value
when the correlation length of the input random eld is small to medium, because of the large number of terms required in the Karhunen-Love expansion
for a fair discretization (Section 3.6).
when the reliability index is large, because of the high order of the polynomial
chaos expansion required.
E 0:3),
the SSFEM approach followed by FORM may not converge or may converge to a
wrong result (Section 3.5).
Finally, it is noted that importance sampling after FORM analysis is inexpensive
to carry out due to the analytical expression of the limit state function. Thus it
allows to rene the evaluation of the probability of failure of the system at low
cost.
2 Slightly greater values can certainly be obtained in a fully compiled implementation. However, this
does not change fundamentally the conclusions.
164
When both methods are employed and give the same results, the SSFEM analysis
can be post-processed to compute the PDF of the response quantity appearing
in the limit state function (see derivations in Chapter 3, Section 5.4). It can also
be used to perform several FORM analyses with dierent limit-state functions.
This seems to be the only case where SSFEM could give something more than the
direct coupling approach. The direct coupling results are needed however in order
to check the validity of the SSFEM solution.
4 Settlement of a foundation over an elastic soil mass Lognormal input random eld
4.1 Introduction
In this section, the direct coupling and SSFEM methods are applied together with a
one-dimensional lognormal random eld modeling the Young's modulus of the material.
As far as direct coupling is concerned, the only dierence with the preceding section
is the way the random eld realizations are evaluated in Femrf: Eq.(4.5) is now used
instead of Eq.(4.4). As far as SSFEM is concerned, the introduction of lognormal elds
requires the stiness matrices to be expanded into the polynomial chaos, as explained
in Part I, Chapter 5, Section 4.1.
It is emphasized that the discretization of the random eld is not exactly identical for the
two approaches. When using the direct coupling, it corresponds to the exponentiation
of a truncated series expansion of a Gaussian eld. When using SSFEM, it corresponds
to a truncated polynomial chaos expansion such as that described in Part I, Chapter 5,
Section 4.1.
The deterministic problem under consideration is the same as in Section 3. The mean
value and coecient of variation of the Young's modulus are
`x = 10000 m
state function is u = 10 cm.
direction being
and
`y = 30
The parametric study presented in this section is limited to the inuence of the orders
of expansion on the reliability index, as well as the threshold in the limit state function.
Indeed, it is believed that the poor results obtained in Section 3 for small correlation
length and/or large coecient of variation of the eld would not be better when a
lognormal eld is considered.
165
and
direct
M
1
3.528
3.452
3.447
3.447
and
SSF EM
4.717
3.714
3.569
3.561
3.560
4.562
3.617
10
3.474
15
3.467
4.539
10
3.606
4.532
15
3.603
Focusing on column #2, it appears that the direct approach gives a 2-digit accuracy for
the reliability index as soon as
Broadly speaking,
pansion
SSF EM
tends to
direct
instance, for
M =1, SSF EM converges to 3.560 instead of 3.528. This comes from the fact
that the representations of the lognormal eld are not identical in the two approaches,
as mentioned above.
M = 2. The reliability
u by
means of direct
coupling and SSFEM (dierent orders of polynomial chaos expansion are used in this
case). Results are reported in Table 5.10.
166
Table 5.10: Inuence of the threshold in the limit state function - Lognormal input
random eld
direct
(cm)
0.473
2.152
10
3.452
12
4.514
15
5.810
20
7.480
30
9.829
# Iterations
SSF EM
0.401
0.477
0.488
0.488
2.481
2.195
2.165
2.166
4.562
3.617
3.474
3.467
6.642
4.858
4.559
4.534
9.763
6.494
5.918
5.846
14.964
8.830
7.737
7.561
25.367
12.655
10.475
10.044
When direct coupling is used, it is observed that the number of iterations required by
the
iHLRF algorithm to get the design point increases slightly with direct , however not as
much as in the case of Gaussian elds (see Table 5.3). The accuracy of the results does
not depend on the value of
in the Karhunen-Love expansion; the results are equal to those given in Table 5.10 with
less than 1% discrepancy).
167
is better than in the Gaussian case. When using a 4-th order polynomial chaos expansion,
the reliability index is obtained within 1-2% accuracy whatever its value.
In other words, SSFEM applied with lognormal random elds appears to be more reliable
than in the case when it is applied with Gaussian elds. This is an interesting property,
since lognormal elds are more suited to modeling material properties. The fact that the
polynomial chaos expansion has to be used also for representing the input eld seems
not deteriorate the accuracy of the results.
and
p.
Table 5.11: Computer processing time required by direct coupling and SSFEM methods Lognormal random elds
CPT Direct
Coupling (")
22.4
31.2
40.5
49.8
58.8
CPT
SSFEM (")
3.5
5.5
8.8
16.6
36.3
6.6
23.1
324.2
2952.6
8.2
188.2
11.6
829.0
13.4
By comparing the results in Table 5.11 with those in Table 5.7, the following conclusions
can be drawn :
168
As far as direct coupling is concerned, almost the same CPT is observed, whether
the random eld is Gaussian or lognormal. This is explained by the fact that the
only dierence between the two calculations is an exponentiation operation each
time the random eld is evaluated.
As far as SSFEM is concerned, the CPT reported in Table 5.11 are much greater
than those reported in Table 5.7. There are two main reasons explaining this
dierence:
if
Ne
polynomial chaos basis, see Eq.(5.58) in Part I). As it has been mentioned
already,
large.
The second level of assembling requires more computational eort since, for
a given pair
(M ; p), there
and
reliability index, the computer processing time is so large that the SSFEM method is
not ecient at all (see numbers printed in bold characters in Table 5.11).
4.5 Conclusions
The parametric study carried out using lognormal input random elds has shown that :
the direct coupling gives accurate results, whatever the value of the reliability
index, at a cost similar to that obtained when using Gaussian elds.
the SSFEM method gives better results with lognormal elds than with Gaussian
random elds. Using a 4-th order polynomial chaos expansion allows to get 1-2%
accuracy on the reliability index. However, the computation time in this case is
huge compared to that of the direct coupling (about 100 times when
M = 2).
Practically the method could be applied only when the correlation length is large,
so that the Karhunen-Love expansion with
M =1
5. Conclusions
169
5 Conclusions
In this chapter, reliability problems have been solved using two dierent methods, namely
the direct coupling between the
and the SSFEM method post-processed by the same algorithm. Both methods have been
applied to assess the serviceability of a foundation lying on an elastic heterogeneous soil
layer. The Young's modulus of the soil was successively modeled as a Gaussian and
lognormal random eld.
The case of Gaussian random elds with exponential autocorrelation function has been
investigated rst, because this type of elds has been extensively used in the literature,
however without convincing comparisons or appreciation of the results. It appears that
a fair discretization of the eld may require more than a few terms, even when the
correlation length is not small. The accuracy in the discretization turns out to be a key
issue. It is noted that this point is never discussed in the papers making use of this kind
of expansions together with SSFEM.
Whatever the parameters, the direct coupling appears robust and fast, the cost of the
analysis increasing linearly with the order of expansion of the input random eld.
As far as SSFEM is concerned, fair results can sometimes be obtained, usually using a
high order polynomial chaos expansion (
the random eld discretization, the cost becomes rapidly prohibitive. Consequently, only
results obtained with a low order polynomial chaos expansion are available in this case.
They appear poor compared to those obtained by direct coupling. In some cases, the
computed reliability index may not even be correct (for instance when large coecients
of variation of the input are considered).
The case of lognormal random elds has been investigated as well. The direct coupling
provides reliable results, approximately at the same cost as in the Gaussian case. The SSFEM approach appears more stable than in the Gaussian case. However the computation
As a conclusion, the direct coupling appears far better than the SSFEM approach for
solving reliability problems, because it is robust and fast. The SSFEM approach could
however be applied
together with the direct coupling in some cases (i.e. when it has proven
Chapter 6
Conclusion
The goal of the second part of the present study was to compare dierent methods
taking into account spatial variability of the material properties in the mechanical analysis. In order to be able to compare a broad spectrum of methods, attention has been
focused on elastic two-dimensional problems. For this purpose, dierent routines have
been developed in the Matlab environment, namely :
the iHLRF algorithm for nding the design point in FORM analysis,
additional routines for perturbation analysis, Monte Carlo simulation and importance sampling.
The detailed implementation of the programs has been presented in Chapters 2-3. An
object-oriented programing was aimed at, in order to get easily extensible code.
The dierent programs were applied to compute the moments of the response of a
foundation lying on an elastic soil layer (up to second order), as well as to assess its
serviceability with respect to a maximum settlement criterion.
As far as second moment analysis is concerned, both the second-order perturbation and
SSFEM methods provide good results, whatever the coecient of variation of the input
Chapter 6. Conclusion
172
random eld. However, the former turns out to be computationally more ecient because
of the special form it takes for the considered application.
As far as reliability analysis is concerned, the direct coupling turns out to be far better
than SSFEM, because it provides excellent accuracy whatever the type of random eld
and the selected parameters. SSFEM can give fair results in some cases, but usually at
a cost much greater than that of the direct coupling. Unfortunately, in some cases, it
converges to a wrong solution, which makes it unreliable.
It is noted that the perturbation method and the direct coupling approach have a far
larger scope than SSFEM, since they have been applied to all kinds of non-linear problems, whereas SSFEM is still more or less limited to linear problems.
As a conclusion, it is noted that the present study is the rst attempt to compare a
broad spectrum of stochastic nite element methods on a given application. Throughout
the description of the implementation, it has been seen that these methods have more
in common than what the dierent research communities involved in their development
sometimes think, at least from a computational point of view. Of course, a single example
(
i.e.,
conclusions of the superiority of some methods over others, but it gives at least a new
light on their respective advantages and shortcomings.
Bibliography
-point algorithm for large time invariant and
time-variant reliability problems, in Reliability and Optimization of Structures, Proc.
3rd WG 7.5 IFIP conference, Berkeley.
Anders, M. and Hori, M., 1999, Stochastic nite element method for elasto-plastic body,
18971916.
Baecher, G.-B. and Ingra, T.-S., 1981, Stochastic nite element method in settlement
predictions,
449463.
ASCE, 110, 3,
J. Eng. Mech.,
357366.
Brenner, C. and Bucher, C., 1995, A contribution to the SFE-based reliability assessment
of non linear structures under dynamic loading,
265273.
Bucher, C.-G. and Bourgund, U., 1990, A fast and ecient response surface approach
for structural reliability problems,
Das, P.-K. and Zheng, Y., 2000, Cumulative formation of response surface and its use
in reliability analysis,
309315.
Proc.
8th ASCE specialty Conference on Probabilistic Mechanics and Structural Reliability.
Bibliography
174
Deodatis, G., 1990, Bounds on response variability of stochastic nite element systems,
565585.
Deodatis, G., 1991, The weighted integral method, I : stochastic stiness matrix,
Mech., 117, 8,
J. Eng.
18511864.
Deodatis, G. and Shinozuka, M., 1991, The weighted integral method, II : response
variability and reliability,
18651877.
Der Kiureghian, A. and de Stefano, M., 1991, Ecient algorithms for second order
reliability analysis,
29062923.
Der Kiureghian, A. and Ke, J.-B., 1985, Finite element based reliability analysis of frame
structures, in
Der Kiureghian, A. and Ke, J.-B., 1988, The stochastic nite element method in structural reliability,
8391.
Der Kiureghian, A., Lin, H.-Z., and Hwang, S.-J., 1987, Second order reliability approximations,
12081225.
Der Kiureghian, A. and Liu, P.-L., 1986, Structural reliability under incomplete probability information,
85104.
Der Kiureghian, A. and Taylor, R.-L., 1983, Numerical methods in structural reliability,
Proc. 4th Int. Conf. on Appl. of Statistics and Probability in Soil and Structural
Engineering, 769784.
in
Der Kiureghian, A. and Zhang, Y., 1999, Space-variant nite element reliability analysis,
173183.
Chichester.
Faravelli, L., 1989, Response surface approach for reliability analysis,
115, 12,
J. Eng. Mech.,
27632781.
a
simulation, J. Appl. Mech., ASME, 65,
Ghanem, R.-G., 1998 , Hybrid stochastic nite elements and generalized Monte Carlo
10041009.
b
dia, Comp. Meth. Appl. Mech. Eng., 158,
Bibliography
175
a
plementation, Comp. Meth. Appl. Mech. Eng., 168,
Ghanem, R.-G., 1999 , Ingredients for a general purpose stochastic nite elements im-
1934.
Ghanem, R.-G., 1999 , The nonlinear gaussian spectrum of log-normal stochastic processes and variables,
c
properties, J. Eng. Mech., ASCE, 125, 1,
964973.
Ghanem, R.-G., 1999 , Stochastic nite elements with multiple random non-Gaussian
2640.
Ghanem, R.-G. and Brzkala, V., 1996, Stochastic nite element analysis of randomly
layered media,
361369.
Ghanem, R.-G. and Kruger, R., 1996, Numerical solution of spectral stochastic nite
element systems,
289303.
Ghanem, R.-G. and Spanos, P.-D., 1990, Polynomial chaos in stochastic nite elements,
Ghanem, R.-G. and Spanos, P.-D., 1991 , Spectral stochastic nite-element formulation
for reliability analysis,
23512372.
Grandhi, R.-V. and Wang, L., 1999, Higher-order failure probability calculation using
non linear approximations,
185206.
Handa, K. and Andersson, K., 1981, Application of nite element methods in the statistical analysis of structures, in
ICOSSAR'81, 409420.
Hasofer, A.-M. and Lind, N.-C., 1974, Exact and invariant second moment code format,
111121.
Hisada, T. and Nakagiri, S., 1981, Stochastic nite element method developed for structural safety and reliability, in
SAR'81, 395408.
Hisada, T. and Nakagiri, S., 1985, Role of stochastic nite element method in structural
safety and reliability, in
SAR'85, 385395.
Hohenbichler, M. and Rackwitz, R., 1981, Non normal dependent vectors in structural
safety,
12271238.
Bibliography
176
Hornet, P., Pendola, M., and Lemaire, M., 1998, Failure probability calculation of an
axisymmetrically cracked pipe under pressure and tension using a nite element code,
Hurtado, J.-E. and Alvarez, D.-A., 2000, Reliability assessment of structural systems
paper #290.
Kim, S.-H. and Na, S.-W., 1997, Response surface method using vector projected sampling points,
319.
Kleiber, M., Antnez, H., Hien, T.-D., and Kowalczyk, P., 1997,
Parameter sensitivity
z_
two-parameter random elds, Prob. Eng. Mech., 13, 3,
Knabe, W., Przewlcki, J., and R yski, 1998, Spatial averages for linear elements for
147167.
Int. J. Num.
18491863.
Lemaire, M., 1997, Finite element and reliability : combined methods by surface response,
lishers.
Lemaire, M., 1998, Elments nis et abilit : un mariage la mode, in A. Mebarki,
D. Boissier, and D. Breysse, Editors,
FIAB'98.
Press.
Li, C.-C. and Der Kiureghian, A., 1993, Optimal discretization of random elds,
Mech., 119, 6,
J. Eng.
11361154.
Liu, P. and Liu, K., 1993, Selection of random eld mesh in nite element reliability
analysis,
667680.
Liu, P.-L. and Der Kiureghian, A., 1986, Multivariate distribution models with prescribed marginals and covariances,
105112.
Bibliography
177
Liu, P.-L. and Der Kiureghian, A., 1989, Finite element reliability methods for geomet-
UCB/SEMM/89-05, University
of California, Berkeley.
a
linear uncertain structures, J. Eng. Mech., ASCE, 117, 8,
Liu, P.-L. and Der Kiureghian, A., 1991 , Finite element reliability of geometrically non
18061825.
Liu, P.-L. and Der Kiureghian, A., 1991 , Optimization algorithms for structural reliability,
Struct. Safety, 9,
161177.
Liu, P.-L., Lin, H.-Z., and Der Kiureghian, A., 1989, Calrel user's manual, Tech. Rep.
a
linear structural dynamics, Comp. Meth. App. Mech. Eng., 56,
Liu, W.-K., Belytschko, T., and Mani, A., 1986 , Probabilistic nite elements for non
6186.
Liu, W.-K., Belytschko, T., and Mani, A., 1986 , Random eld nite elements,
Int. J.
18311845.
Addison &
Struct. Safety, 5,
3545.
Mahadevan, S. and Haldar, A., 1991, Practical random eld discretization in stochastic
nite element analysis,
Struct. Safety, 9,
283304.
Matthies, G., Brenner, C., Bucher, C., and Guedes Soares, C., 1997, Uncertainties in
probabilistic numerical analysis of structures and solids - stochastic nite elements,
283336.
Pellissetti, M. and Ghanem, R.-G., 2000, Iterative solution of systems of linear equations
arising in the context of stochastic nite elements,
31, 8-9,
607616.
Pendola, M., Chtelet, E., and Lemaire, M., 2000 , Finite element and reliability: combined method by response surface using neural networks,
Struct. Safety,
submitted
for publication.
Pendola, M., Hornet, P., and Mohamed, A., 2000 , Parametric study of the reliability
of a cracked pipe using a non-linear nite element analysis, in R.E Melchers and M.G
Bibliography
178
Pendola, M., Mohamed, A., Lemaire, M., and Hornet, P., 2000 , Combination of nite element and reliability methods in nonlinear fracture mechanics,
Reliability Engineering
1527.
Phoon, K., Quek, S., Chow, Y., and Lee, S., 1990, Reliability analysis of pile settlements,
17171735.
Rackwitz, R. and Fiessler, B., 1978, Structural reliability under combined load sequences,
Comput. Struct., 9,
489494.
Rajashekhar, M.-R. and Ellingwood, B.-R., 1993, A new look at the response surface
approach for reliability analysis,
Spanos, P.-D. and Ghanem, R.-G., 1989, Stochastic nite element expansion for random
media,
a
ment analysis, Prob. Eng. Mech., 5, 4,
10351053.
Takada, T., 1990 , Weighted integral method in multidimensional stochastic nite ele158166.
Takada, T., 1990 , Weighted integral method in stochastic nite element analysis,
Eng. Mech., 5, 3,
Prob.
146156.
Massachussets.
Vanmarcke, E.-H. and Grigoriu, M., 1983, Stochastic nite element analysis of simple
beams,
12031214.
Waubke, H., 1996, Dynamische Berechnungen fr den Halbraum mit streuenden Parametern mittels orthogonaler Polynome, Ph.D. thesis, Technische Universitt, Mnchen.
Wong, F.-A., 1985, Slope reliability and response surface method,
111, 1,
J. Geotech., ASCE,
3253.
268287.
Zeldin, B.-A. and Spanos, P.-D., 1998, On random eld discretization in stochastic nite
elements,
320327.
Zhang, J. and Ellingwood, B., 1994, Orthogonal series expansion of random elds in
reliability analysis,
26602677.
Bibliography
179
Zhang, Y. and Der Kiureghian, A., 1993, Dynamic response sensitivity of inelastic structures,
2336.
Zhang, Y. and Der Kiureghian, A., 1995, Two improved algorithms for reliability analysis,
in R. Rackwitz, G. Augusti, and A. Borr, Editors,
Zhang, Y. and Der Kiureghian, A., 1997, Finite element reliability methods for inelastic
McGraw-Hill,