- Email: [email protected]
A new method in static structural reliability
Probabilistic Engineering Mechanics 17 (2002) 317–325 www.elsevier.com/locate/probengmech
A new method in static structural reliability Kiran B. Chilakamarri1 CK & Associates, 1393 Tabor Avenue, Dayton, OH 45420, USA
Abstract A new method to determine the static structural reliability is presented in this paper. This method takes advantage of the fact that isoprobabilistic transformation, in particular the Rosenblatt’s transform, maps infinite Euclidean space into a unit hypercube and in doing so, reduces the problem of computing probability content of infinite region to a problem of computing a finite volume. The method is applied to several examples including a well-known example with multiple design points. q 2002 Published by Elsevier Science Ltd. Keywords: Reliability; Failure probability; Structures
1. Introduction In structural reliability the safety of a structure is the probability content of a set of safe points determined by a performance function in the parameter space representing the uncertain physical parameters of the structure. The performance function or the limit-state function may or may not have a closed form. Here we assume a limit-state function with a closed form. A fundamental problem in evaluating the probability of failure is the computation of the resulting multi-fold integral. The number of computations increases exponentially with the number of variables. Very frequently an isoprobabilistic transformation or probability transformation is used to standardize the model to convert the joint pdf of the parameters to a joint pdf of independent standard normal variables. A key observation made here is that the probability transformation T (Rosenblatt’s transform in particular) with k random variable maps R k (k-dimensional real space) into a unit hypercube Qk in such a way that for any measurable set S in R k, the volume of the image set TðSÞ is the probability content of S. A simple consequence of this observation is that it is possible to evaluate Pf (the probability of failure) directly in Qk and avoid integration over infinite sets. If in addition the limit-state function admits separation of variables into non-interacting subgroups then the problem of computing multi-fold integral can be split into several sub-problems of computing lower dimensional integrals. The term ‘pulverizer’ is used here to indicate this process. If E-mail address: [email protected] (K.B. Chilakamarri). Now a Visiting Associate Professor in the Department of Mathematics and Statistics, Wright State University, Dayton, OH 45435, USA.
a limit-state function can be pulverized perfectly then the multi-fold integral can be evaluated by calculating k 2 1 one-dimensional integrals where k is the number of random variables. On the basis of these observations an algorithm to evaluate Pf and cdf of performance function is developed. In Section 2, Rosenblatt’s transformation [1] is discussed and a formula for Pf is developed. In the same section the formula is applied to two examples, the first example is a wellknown problem with multiple design points and the second example is an artificial example with highly non-linear limit-state function. In Section 3, the algorithm is described and applied to a simple problem with three variables. In Section 4, the method is applied to three problems each with more than three variables. Example 4 is a well-known problem of one plastic collapse mechanism of one-bay frame. In Example 5, we find the cdf of maximum radial stress of rotating disk and in Example 6 we find the cdf of the natural frequency of a plate. Throughout this paper the results are compared with the existing Monte Carlo simulation results whenever possible.
2. The probability transformation Let ðX1 ; X2 ; …; Xk Þ be a random vector with distribution function Fðx1 ; x2 ; …; xk Þ ¼ Prob½X1 # x1 ; X2 # x2 ; …; Xk # xk : Let Qk denotes the k-dimensional unit hypercube, i.e. Qk ¼ {ðz1 ; z2 ; …; zk Þ;
0 # z1 # 1; …; 0 # zk # 1}:
1
0266-8920/02/$ - see front matter q 2002 Published by Elsevier Science Ltd. PII: S 0 2 6 6 - 8 9 2 0 ( 0 2 ) 0 0 0 1 5 - 2
The probability transformation T: Rk ! Qk also known as
318
K.B. Chilakamarri / Probabilistic Engineering Mechanics 17 (2002) 317–325
Rosenblatt transformation is defined as follows
Carlo simulation as reported by Madsen et al. [3] has a value of 0.00294.
Tðx1 ; x2 ; …; xk Þ ¼ z ¼ ðz1 ; z2 ; …; zk Þ
Example 1. Limit-state function is gðx1 ; x2 Þ ¼ 18 2 3x1 2 2x2 and the joint cumulative function is given by
where z1 ¼ F1 ðx1 Þ ¼ Prob½X1 # x1 ; z2 ¼ F2 ðx2 lx1 Þ ¼ Prob½X2
Fððx1 ; x2 Þ ¼ 1 2 expð2x1 Þ 2 expð2x2 Þ
# x2 lX1 ¼ x1 ; z3 ¼ F3 ðx3 lx2 ; x1 Þ ¼ Prob½X3 # x3 lX1
þ expð2ðx1 þ x2 þ x1 x2 ÞÞ;
¼ x1 ; X2 ¼ x2 ; …; zk ¼ Fk ðxk lxk21 ; …; x1 Þ ¼ Prob½Xk
x1 ; x2 . 0:
ð1Þ
As described by Hohenbichler and Rackwitz [2], there are 2! ¼ 2 distinct transformations. Note that the random variables are not independent for this example.
The random variables Z1 ; Z2 ; …; Zk are uniformly and independently distributed over ½0; 1 and the random vector ðZ1 ; Z2 ; …; Zk Þ is uniformly distributed over the k-dimensional unit hypercube Qk [1]. If in addition, the random variables X1 ; X2 ; …; Xk are independent then
Transformation I: Tðx1 ; x2 Þ ¼ ðz1 ; z2 Þ where, z1 ¼ 1 2 e2x1 and z2 ¼ 1 2 ð1 þ x2 Þ e2x2 ð1þx1 Þ : The image of the curve g ¼ 0 in the Z-square (a unit square) is given by
ðz1 ; z2 ; …; zk Þ ¼ Tðx1 ; x2 ; …; xk Þ
z2 ¼ 1 2 ½10 þ ð3=2Þ lnð1 2 zÞ exp½2{9 þ ð3=2Þ
# xk lX1 ¼ x1 ; X2 ¼ x2 ; …; Xk21 ¼ xk21 :
lnð1 2 zÞ}{1 2 lnð1 2 zÞ}
¼ ðFX1 ðx1 Þ; FX2 ðx2 Þ; …; FXk ðxk ÞÞ: Note that different permutations of X1 ; X2 ; …; Xk can produce different transformations. For the method proposed in paper it is necessary to be able to invert the Rosenblatt’s transformation. Given z1 ; z2 ; …; zk ; we first find x1, after x1 is evaluated both z2 and x1 together yield x2, after x2 is evaluated z3, x1 and x2 can be used to find x3, and so on. The reader is urged to note that the independence of random variables is not required for this process. We assume the set S of failure points is characterized by non-negative values of function gðxÞ the limit-state function. We may define a corresponding failure set Sp ¼ TðSÞ in Qk. Accordingly, the failure surface ›S, the boundary of S, is mapped onto the boundary ›S p. Thus we can calculate the probability of failure either in R k or in the hypercube Qk. Since z is uniformly distributed in Qk ð ð Pf ¼ dFðxÞ ¼ dz ¼ VolumesðSp Þ: ð2Þ Sp
S
Once the boundary ›S p is identified in Qk the problem of evaluating Pf is reduced to a problem of evaluating the volume of a finite domain in Z-space. If the failure region can be characterized by xk , hðx1 ; x2 ; …; xk21 Þ and Tðx1 ; x2 ; …; xk21 Þ denotes the restriction of the map T to the first k 2 1 coordinates then Pf ¼
ð1 ð1 0
0
where z denotes z1 and the failure region in Z-square corresponding to g , 0 is the region above this curve. Thus Pf ¼
ð3Þ Example 1 with multiple design points is a well-known problem in FORM and SORM and originally introduced by Hohenbichler and Rackwitz [2] and discussed by Madsen et al. [3], Ditlevsen and Madsen [4] and Der Kiureghian and Dakessian [5]. The exact solution of the problem by Monte
½10 þ ð3=2Þ lnð1 2 zÞ exp½2{9 þ ð3=2Þ
0
lnð1 2 zÞ}{1 2 lnð1 2 zÞ} dz þ e26 ¼ 0:00294485974: The value of Pf is accurate up to 10 digits. To see how this integral is written, note that we are interested in finding the area of the image of the region g , 0 under the probability transformation in Z-square. For a given z1 ð0 , z1 , 1Þ; we first find the corresponding x1 value from z1 ¼ 1 2 e2x1 : Having found x1, the x2 value is calculated from g ¼ 0: Once x2 is evaluated we find z2 corresponding to z1 from z2 ¼ 1 2 ð1 þ x2 Þ expð2x2 ð1 þ x1 ÞÞ: Note that the support of the image of g ¼ 0 is the interval ½0; 1 2 e26 Þ and the area of interest is the area above the image curve. Transformation II: Tðx1 ; x2 Þ ¼ ðz1 ; z2 Þ; where z1 ¼ 1 2 e2x2 and z2 ¼ 1 2 ð1 þ x1 Þ e2x1 ð1þx2 Þ : Image of the curve g ¼ 0 in Z-square is given by z2 ¼ 1 2 ½7 þ ð2=3Þ lnð1 2 zÞ exp½2{6 þ ð2=3Þ lnð1 2 zÞ}{1 2 lnð1 2 zÞ}
ð1 · · · T h T 21 ðz1 ; z2 ; …; zk21 Þ dz1 dz2 · · ·dzk21 : 0
ð12e26
and again Pf ¼
ð12e29
½7 þ ð2=3Þ lnð1 2 zÞexp½2{6 þ ð2=3Þ
0
lnð1 2 zÞ}{1 2 lnð1 2 zÞ} dz þ e29 ¼ 0:00294485974:
K.B. Chilakamarri / Probabilistic Engineering Mechanics 17 (2002) 317–325
Example 2. Let X and Y be independent standard normal variables. Let the failure region be given by
admit such decomposition. For instance, if
4 ; y$ x{1 þ sin2 ðxÞ}
g¼ x . 0:
Probability transformation is given by ðz1 ; z2 Þ ¼ Tðx; yÞ ¼ ðFðxÞ; FðyÞÞ: The g-function is gðx; yÞ ¼ y 2
4 x{1 þ sin2 ðxÞ}
and Pf ¼ Prob½g $ 0: The failure curve is described by yðxÞ ¼
4 ; x{1 þ sin2 ðxÞ}
then Pf ¼
ð1 0
12F
4 dz: F21 ðzÞð1 þ sin2 ðF21 ðzÞÞÞ
To see why, for a point z in the unit interval, F21 ðzÞ is the corresponding value of x, the y value from the failure curve is yðF21 ðzÞÞ and the area above the curve is Pf. The integral for Pf can be approximated by a finite integral using ^ 5 sigma bounds, further the region of interest can be restricted to {ðx; yÞ : x [ ½0:479177; 5 and y [ ½0; 5} where xp ¼ 0:479177 is the x-value when y ¼ 5 on the curve gðx; yÞ ¼ 0:Thus ðFð5Þ 4 Pf , 12F dz F21 ðzÞð1 þ sin2 ðF21 ðzÞÞÞ Fðxp Þ ¼ 0:0121736: This result is accurate up to six digits. For more than two variables we develop an algorithm described in Section 3.
3. The algorithm We need to introduce new terms to describe the algorithm. Given an expression gðX1 ; X2 ; …; Xk Þ in k variables divide the set of variables {X1 ; X2 ; …; Xk } into two subgroups say {X1 ; X2 ; …; Xi } and {Xiþ1 ; Xiþ2 ; …; Xk } so that g can be expressed as combination of two expressions g 1 and g2 , where g1 involves only X1 ; X2 ; …; Xi and g2 involves only Xiþ1 ; Xiþ2 ; …; Xk : Simply we want to find two expressions g1 ðX1 ; …; Xi Þ and g2 ðXiþ1 ; …; Xk Þ so that gðX1 ; X2 ; …; Xk Þ ¼ Fðg1 ; g2 Þ for some function f. For instance if g ¼ X1 ðX2 þ X3 Þ; then g1 ¼ X1 and g2 ¼ X2 þ X3 : If g ¼ ðX1 þ 3ÞX2 X32 ðX42 2 X52 Þ many choices of g1 and g2 are possible. Some expressions may not
319
X1 ; X2 þ X3 ðX1 þ 10Þ
then no two variables can be separated from the third. Repeat the process with g1 and g2, and so on until the process is terminated. I call this method pulverizer. The end product of the pulverizer method is a record in the form of a Pulverizer Tree. This is a rooted tree in which the root is g, has two children g1 and g2, and the nodes g1 and g2 contain their children and so on (see Fig. 1, 3, 7 and 9). The pulverizer tree provides a quick reference for the reduction of g into smaller expressions. Several pulverizer trees may exist for a given expression. The number of integrals to be evaluated is exactly the number of nodes with degree more than one. Next we introduce the concept of Working Area. Given a real number 1 . 0 and the joint distribution F of k random variables the working area in R k is defined as a region W with Prob½W $ 1 2 1: Similarly, the corresponding working area in Z-space (the hypercube) is TðWÞ ¼ W 0 with VolumeðW 0 Þ $ 1 2 1: As the name indicates, when solving a problem we do not venture outside the working area, so the probability content of the complement of W makes a contribution of at most 1 to the error. For example with two independent normal distributions, X , Nð10; 4Þ and Y , Nð2; 0:1Þ; using ^ 3 sigma bounds, the working area W in R 2 is given by 22 # x # 22 and 1:7 # y # 2:3; and in Zsquare the working area is simply ½Fð23Þ; Fð3Þ £ ½Fð23Þ; Fð3Þ: The error from the ignorance of outside of W is exactly 1 2 ½Fð3Þ 2 Fð23Þ2 : If we wish to use ^ 20 sigma bounds then similar analysis produces a working area with associated error ¼ 1 2 ½Fð20Þ 2 Fð220Þ2 : For an additional example consider two independent variables, X , Nð10; 4Þ and Y , Uð10; 121Þ a uniform distribution. Suppose again for X we use ^ 3 sigma bound, then the working W area in R 2 is, 22 # x # 22 and 10 # y # 121: In Z-square, the working area is ½Fð23Þ; Fð3Þ £ ð0; 1Þ: The error associated with W is 1 2 ½Fð3Þ 2 Fð23Þ: Note that the working area can be extended easily since the concept is additive. Also note that any information about the outside of working area may be used to reduce this type of error. The algorithm I. II.
Write the probability transformation. Decide on a suitable pulverizer tree for the given limitstate function. (Here we want to avoid nodes of high
Fig. 1. Pulverizer tree for Example 3.
320
K.B. Chilakamarri / Probabilistic Engineering Mechanics 17 (2002) 317–325
degree. We also want to choose a short and broad tree). III. To find the cdf of a node labeled hðu; vÞ with a pair of nodes u and v attached, in the working area study the level curves of hðu; vÞ ¼ c (or u ¼ Fðu; cÞ) to write an expression for Fh ðcÞ ¼ Prob½u , Fðu; cÞ using the formula ð1 ð4Þ Fh ðcÞ ¼ Fu f Fv21 ðzÞ; c 0
Any integration technique can be used so long as it also provides error estimates, i.e. functions UðcÞ and LðcÞ so that LðcÞ # Fh ðcÞ # UðcÞ for all c. IV. Start at the bottommost nodes and keep evaluating the cdf of next higher level nodes until we reach the top.
Remark 1. There are two sources of error in the method presented here. As already mentioned, restricting the working area results in a predetermined amount of error. Another source of error is from numerical integration. However, both types of errors can be controlled. Remark 2. Suppose the number of variables is k. If a pulverizer tree has root with degree d p and has d3 number of vertices of degree 3, d4 number of vertices of degree 4,…, and dR number of vertices of degree R. The number of vertices with degree 1 is k, same as the number of variables. Using elementary counting techniques we can show that k ¼ dp þ d3 þ 2d4 þ · · · þ ðR 2 2ÞdR : A non-root vertex of degree iði . 2Þ requires a ði 2 2Þdimensional integration to compute one value of the cdf associated with that vertex and the root requires a ðdp 2 1Þdimensional integration. If N is the grid size, then the total number of computations required is of the order p O d3 N þ d4 N 2 þ · · · þ dR N R22 þ N d 21 : In general if a limit-state function of k variable admits complete pulverization, i.e. d4 ¼ d5 ¼ · · · ¼ 0; then d3 ¼ k 2 2 and dp ¼ 2: Thus Oððk 2 1ÞNÞ computations are required. To reduce the number of computations, trees with small degrees are preferred. Computational errors would compound along a long branch, thus short pulverizer trees are preferred.
EðSÞ ¼ 2500 kg/cm 2, EðAÞ ¼ 2 cm 2, EðFÞ ¼ 4166 kg, sðSÞ ¼ 75 kg/cm2, sðAÞ ¼ 0:04 cm2, and sðFÞ ¼ 125 kg. Failure occurs if S 2 ðF=AÞ # 0: Let A ¼ X1 , Nð2; 0:04Þ; S ¼ X2 , Nð2500; 75Þ and F ¼ X3 , Nð4166; 125Þ and the limit-state gðX1 ; X2 ; X3 Þ ¼ X2 2 ðX3 =X1 Þ: Probability transformation is given by x 22 ðZ1 ; Z2 ; Z3 Þ ¼ FX1 ; FX2 ; FX3 ¼ F 1 ; 0:04
F
x2 2 2500 ; 75
x 2 4166 F 3 125
Writing m ¼ X3 =X1 and g ¼ X2 2 m the pulverizer tree is shown in Fig. 1. First we find Fm, the cdf of m, from FX1 and FX3 and then evaluate the cdf of g from Fm andFX2 : It is not hard to show that if x and y are independent standard normal variables, then Ay þ B ,c Prob ax þ b is the probability content of the region R in R 2 shown in Fig. 2. This probability is exactly equal to the area of the region R in the Z-square. Using this result in our context, we have
125y þ 4166 ,d ; Fm ðdÞ ¼ Prob½ðX3 =X1 Þ , d ¼ Prob 0:04x þ 2 where y and x are standard normal variables. For this problem, ð2b=aÞ ¼ 250; and Fð250Þ , 0 and y ¼ ðad=AÞx þ ðdb 2 BÞ=A ¼ 0:00032 dx þ ð2d 2 4166Þ=125; thus we have ð1 2d 2 4166 21 Fm ðdÞ ¼ F 0:00032dF ðzÞ þ dz 125 0 and using ^ 5 sigma bounds the relevant range of d is
4. Examples We begin with a simple example with three variables [6]. Example 3. Consider a bar of cross-section A and a yield stress S is subjected to a tensile force F where A, S and F are independently distributed normal random variables with
Fig. 2. The area of the marked region in Z-square is the probability content of the region R in x –y plane.
K.B. Chilakamarri / Probabilistic Engineering Mechanics 17 (2002) 317–325
321
below the diagonal of the square working area ð144 # c # 264Þ and (ii) those that are above the diagonal ð269 # c # 394Þ: We can now write the integral for the cdf of m Fm ðcÞ ¼
¼12
ðFX1 ðxp Þ 0
ð1 FX1 ðxp Þ
Fig. 3. Pulverizer tree for Example 4.
½1609; 2661:69:We are now ready to find Fg where g ¼ X2 2 m and FX21 ðzÞ ¼ 75F21 ðzÞ þ 2500 2 Fg ðcÞ ¼ Prob½g # c ¼ Prob½m $ X2 2 c ¼
ð1 0
¼
ð1 0
1 2 Fm FX21 ðzÞ 2 c dz 2
1 2 Fm ½75F21 ðzÞ þ 2500 2 c dz
¼ ð5:7143Þ £ 1025 The FORM analysis [6] yields a value of (5.583) £ 1025. Example 4 is a well-known example [3,7]. Example 4. The failure region of one plastic collapse mechanism of one-bay frame is given by X1 þ 2X2 þ 2X3 þ X4 2 5X5 2 5X6 , 0: The variables X1 ; …; X6 are statistically independent and log-normally distributed with mX1 ¼ mX2 ¼ mX3 ¼ mX4 ¼ 120; mX5 ¼ 50; mX6 ¼ 40 and sX1 ¼ sX2 ¼ sX3 ¼ sX4 ¼ 12; sX5 ¼ 15; and sX6 ¼ 12: The probability transformation is described below ln t 2 4:78252 FX1 ðtÞ ¼ FX2 ðtÞ ¼ FX3 ðtÞ ¼ FX4 ðtÞ ¼ F ; 0:0997513 ðzÞ ¼ FX21 ðzÞ ¼ FX21 ðzÞ ¼ FX21 ðzÞ FX21 1 1 3 4 21
FX 5
¼ expð4:78252 þ 0:0997513F ðzÞÞ; ln x5 2 3:86893 ¼F ; 0:29356
ðzÞ ¼ expð3:86893 þ 0:29356F21 ðzÞÞ; FX21 5 ln x5 2 3:64579 FX 6 ¼ F ; 0:29356
FX4 c 2 FX21 ðzÞ dz 1
for 144 # c # 264
1 2 FX4 c 2 FX21 ðzÞ dz for 264 # c # 394 1
where xp ¼ c 2 72 and xp ¼ c 2 197: Next step involves evaluating F› but Fm ¼ F› since X1 ; X2 X3 ; X4 are all identically distributed. Now we can evaluate Fb the cdf of b ¼ m þ 2›: For the computational efficiency a table of values of Fm is computed, stored and used to approximate the integrals that appear in the formula for Fb. The significant domain for Fm ðcÞ turns out to be a smaller interval ½170; 300 than that given by working area bounds. The interval ½170; 300 is partitioned so that variation of Fm is sufficiently small. Here, a partition 170 ¼ c0 , c1 , · · · , c300 ¼ 340 is used so that lFm ðci Þ 2 Fm ðciþ1 Þl , 0:005192: Once again using ^ 5 sigma bounds the working area for b ¼ m þ 2› in the m – › plane is ½170; 340 £ ½170; 340 and the level curves are shown in Fig. 4. The integral for Fb is easy to write Fb ðcÞ ¼ 8 ! ðFm ðc2340Þ > c 2 Fm21 ðzÞ > > dz Fm > > 2 > 0 > > > ! > < ð1 c 2 Fm21 ðzÞ Fm dz > 2 0 > > > ! > > ð1 > c 2 Fm21 ðzÞ > > dz > : 1 2 F ðc2680Þ 1 2 Fm 2
if 510 # c # 680
if 680 # c # 850
if 850 # c # 1020
m
Once again we construct a table of values for Fb for computational purposes. Here Fb values were stored as {ðci ; zi Þ : 1 # i # 355} with Fb ðci Þ ¼ zi with Maxlziþ1 2 zi l , 0:00538458: Next, we evaluate FD the cdf of D ¼ x5 þ x6 from the distributions of X5 and X6. Using ^ 5 sigma bounds the working area in the x6 2 x5 plane is restricted to
ðzÞ ¼ expð3:64579 þ 0:29356F21 ðzÞÞ: FX21 6 Writing m ¼ X1 þ X4 ; › ¼ X2 þ X3 ; b ¼ m þ 2›; D ¼ X5 þ X6 ; and g ¼ b 2 5D; we choose a pulverizer tree as shown in Fig. 3. We first evaluate Fm ðcÞ the cdf of m. Assuming ^ 5 sigma bounds, the working area in x1 2 x4 plane is a square ½72; 197 £ ½72; 197; and the level curves of m (i.e. x1 þ x4 ¼ c) are straight lines for 144 # c # 394: We distinguish between the two kinds of lines (i) those that are
Fig. 4. Level curves of b in the working area.
322
K.B. Chilakamarri / Probabilistic Engineering Mechanics 17 (2002) 317–325
Examples 5 and 6 are from the article ‘Accuracy of probabilistic methods for HCF’ by Eric Fox in the Proceedings of the 5th National Turbine Engine High Cycle Fatigue Conference HCF 2000, 7 –9 March 2000, Chandler, Arizona. We evaluate the cdf of the expressions involved in each example and compare the results with Monte Carlo results from the above mentioned article.
Fig. 5. Level curves of D in the working area.
½8; 167 £ ½11; 208: The level curves of D are shown in Fig. 5. Since FD ðcÞ ¼ Prob½D # c ¼ Prob½x5 þ x6 # c; we need to find the area in Z-square of the region under the image curve of x5 # c 2 x6 : FD ðcÞ ¼ 8 ðF ðc211Þ X6 > 21 > F c 2 F ðzÞ dz > X5 X6 > > 0 > > > < ð1 FX5 c 2 FX21 ðzÞ dz 6 > 0 > > > > ð1 > > 21 > F c 2 F ðzÞ dz : X5 X6
for 19 # c # 178 for 178 # c # 216 for 216 # c # 375
FX6 ðc2208Þ
where the parameters are assumed to be independent and are distributed as follows: n ¼ Poisson’s ratio ¼ X1 , Uð0:29; 0:31Þ; r ¼ density ¼ X2 , Uð0:279; 0:289Þ; v ¼ rotor speed ¼ X3 ¼, Uð10 000; 11 000Þ; ro ¼ outer disk radius ¼ X4 , Nð8; 0:02Þ and ri ¼ inner disk radius ¼ X5 , Nð2; 0:01Þ where Uða; bÞ denotes uniform distribution over the interval ½a; b: Writing m ¼ ð3 þ X1 ÞX2 ; f ¼ mX32 ; D ¼ X42 2 X52 and g ¼ c0 f D where c0 is a constant comprising of all constants involved in s, we choose the pulverizer tree shown in Fig. 7. The probability transformation is T ¼ ðz1 ; z2 ; …; z5 Þ ¼ ðFX1 ; FX2 ; …; FX5 Þ where
To complete the problem we find Pf ¼ Fg ð0Þ ¼ Fg ðb , 5DÞ ¼ 1 2
Example 5. The maximum radial stress of a rotating disk, s, in psi, is given by 3þn r 2p 2 2 s¼ v ro 2 ri2 ; 60 8 ð9:81Þð39:37Þ
ð1 0
!
FD
Fb21 ðzÞ dz 5
The graph of g ¼ 0; i.e. b 2 5D ¼ 0 and the region g # 0 are shown in Fig. 6. Mid-point Reimann sum is used to approximate the integral and we find Pf ¼ 0:012224 and from the lower and upper Riemann sums we find the lower and upper bounds PLf ¼ 0:012136; PU f ¼ 0:0123135; respectively. Using more refined partition for Fb with Maxlziþ1 2 zi l , 0:00222917; we find Pf ¼ 0:0122188 and the bounds PLf ¼ 0:0121782 and PU f ¼ 0:0122599:
FX21 ðzÞ ¼ 1 ¼
z z þ 0:29; FX21 þ 0:279; FX21 ðzÞ ¼ ðzÞ 2 3 50 100 z þ 10 000; FX21 ðzÞ 4 0:001
¼ 0:02F21 ðzÞ þ 8; FX21 ðzÞ ¼ 0:01F21 ðzÞ þ 2 5 for 0 , z , 1: We can write the integral for Fm ðcÞ ! ð1 c Fm ðcÞ ¼ FX2 dz: 3 þ FX21 ðzÞ 0 1 Since FX2 is not continuous, some caution must be exercised in evaluating the above integral. A study of the level curves of m ¼ ð3 þ X1 ÞX2 help in writing the integrand explicitly
Fig. 6. The failure region g , 0:
Fig. 7. Pulverizer tree for Example 5.
K.B. Chilakamarri / Probabilistic Engineering Mechanics 17 (2002) 317–325
and the integral can be evaluated exactly 8 0 > > > > > 4589:55 2 3:0514:5c þ 5000c ln½179:211c > > < Fm ðcÞ ¼ 227:9 þ 30:3031c > > > > 24781:95 þ 30 544:8c 2 5000c ln½173:01c > > > : 1
323
if c , 0:91791 if 0:91791 , c , 0:923492 if 0:923492 , c , 0:95081 if 0:95081 , c , 0:95659 if c $ 0:95659
Next, Ff the cdf of f ¼ mX32 is obtained from the integral 0 1 ð1 c B C Ff ðcÞ ¼ Fm @ 2 Adz: 21 0 FX3 ðzÞ
image of the area under the curve of type I and FD ðcÞ is obtained by integrating (rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi) 2 FX21 ðzÞ þc FX 4 5
Once again a quick study of the level curves of f in the pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi working area helps in writing the explicit integrals for Ff as from z ¼ FX21 ð ð7:906Þ2 2 cÞ to z ¼ 1: 5 follows Now we can complete the problem by finding Fs the cdf of 8 0 if c , 9:1791 £ 107 > > > p > ðv > c > > Fm dz if 9:1791 £ 107 , c , 9:5659 £ 107 > > 2 > ð1000z þ 10 000Þ 0 > > > < ðvp c Ff ðcÞ ¼ F dz þ up if 9:5659 £ 107 , c , 1:11067 £ 107 m 2 > p ð1000z þ 10 000Þ u > > > 2 > > ð1 > c > > F dz þ up if 1:11067 £ 107 , c , 1:115477 £ 108 > > up m ð1000z þ 10 000Þ2 > > : 1 otherwise where vp ¼ 0:001
qffiffiffiffiffiffiffiffiffiffiffiffiffi c=0:91791 2 10 000 and
qffiffiffiffiffiffiffiffiffiffiffiffiffi u ¼ 0:001 c=0:95659 2 10 000 p
These integrals can be evaluated to provide an exact formula for Ff. Next we evaluate the FD the cdf of D ¼ X42 2 X52 : Using ^ 10 sigma bounds we can find the working area, a rectangle, in the x5 2 x4 plane as shown in Fig. 8. For a given c qffiffiffiffiffiffiffiffiffi D , c $ x24 , x25 þ c $ x4 , x25 þ c:
s ¼ g ¼ c0 f D from the following integral ! ð1 c dz: P½g , c ¼ Fg ðcÞ ¼ Ff c0 FD21 ðzÞ 0 Once again a study of the level curves of g in f 2 D space provides a valuable service in restricting the range for numerical integration. Riemann sums were used to approximate the integral and for that purpose a table for FD was developed so that the Max½Dz , 0:00399814 in evaluation of Fg. For Ff the exact formula was used hence there are no errors from Ff. A sample of results of this analysis and the corresponding Monte Carlo simulation results from the referred article are shown in Table 1. Example 6. The natural frequency of a plate, f, in Hz is
The cdf FD ðcÞ is calculated by integrating (rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi) 2 FX 4 FX21 ðzÞ þc 5
over the interval ð0; 1Þ: We notice that the level curves of D are of three types and accordingly some care must be taken for numerical integration process. For instance if 58:3136 , c , 58:6909; then the area in Z-square corresponds to the
Fig. 8. Level curves of D.
324
K.B. Chilakamarri / Probabilistic Engineering Mechanics 17 (2002) 317–325
Table 1 The cdf of radial stress s of rotating disk
Table 2 The cdf of the natural frequency f of a plate
Percentile
Current new method s
Monte Carlo simulation s
Percentile
Current new method f
Monte Carlo simulation f
0.001 0.01 0.1 0.5 0.9 0.99 0.999
19 19 20 22 23 24 24
19 19 20 22 23 24 24
0.001 0.01 0.1 0.5 0.9 0.99 0.999
1751.4 1874.78 2050.05 2274.12 2508.51 2707.46 2858.2
1752.6 1876.9 2049.8 2274.2 2507.1 2706.3 2858.7
598.5 814 355 002 715 336 594
595.5 810.5 349.8 005.4 714.6 330.9 601
calculated as 0 f ¼
34:03 B @ 2pa2 12
3
11=2
Eh C A g 2 ð1 2 n Þ 386:1
where the parameters are assumed to be independent and are distributed as follows, a ¼ length of plate ¼ X1 , N ð12; 0:1Þ; E ¼ modulus of elasticity ¼ X2 , Nð30 £ 106 , 9 £ 105 Þ; h ¼ plate thickness ¼ X3 , Nð1; 0:05Þ; g ¼ mass per unit area ¼ X4 , Nð0:29; 0:004Þ and n ¼ Poisson’s ratio ¼ X5 , Uð0:29; 0:31Þ; where Uða; bÞ denotes uniform distribution over the interval ½a; b: Writing c ¼ ðX33 =X14 Þ; m ¼ cX2 and w ¼ X4 ð1 2 X52 Þ and combining all the constants appearing in the expression for the frequency, we can write f ¼ 30:7214
rffiffiffiffi m w
Next we need to evaluate ! ð1 d dz; Fm ðdÞ ¼ FX2 Fc21 ðzÞ 0 where FX2 ðtÞ ¼ Fððt 2 30 £ 106 Þ=9 £ 105 Þ: Here again we use ^ 5 sigma bounds. The range of d is the interval ð637:5; 3243Þ: To find the cdf of w we use ^ 5 sigma bounds for X4 the only normal variable appearing in w. The working area is {ðx5 ; x4 Þ : 0:29 # x5 # 0:31 and 0:27 # x4 # 0:31}: A quick study of level curves of w shows that the support of the integrand in the formula for the cdf of w, changes between the intervals ð0:244053; 0:247293Þ; ð0:247293, 0:280209Þ; and ð0:280209; 0:283929Þ. From this information it is not hard to evaluate Fw ðcÞ; where 0 1 c ð1 B 20:29 þ 1 2 ððz=50Þ þ 0:29Þ2 C Cdz: Fw ðcÞ ¼ wB @ A 0:004 0 Finally we can write
and a pulverizer tree chosen for this problem shown in Fig. 9. Using ^ 5 sigma bounds Fc ðcÞ ¼
ð1 0
FX 3
! rffiffiffiffiffiffiffiffiffiffiffiffiffi 4ffi 3 c FX21 ðzÞ dz 1
p3 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! ðFð5Þ c½0:1F21 ðzÞ þ 124 2 1 dz ¼ F 0:05 Fð25Þ ð0:0001728 # c # 0:000111671Þ:
ð30:7214Þ2 m Ff ðcÞ ¼ Prob½F , c ¼ Prob w . c2
¼
ð1 0
(
1 2 Fw
30:7214 c
2
!
!) Fm21 ðzÞ
dz
Again, the nature of the support of the integrand changes at c ¼ 1440:22; 1553.43, 3132.52, and 3378. The values of Ff are computed using Riemann sum approximation for the integral. The integral formula for Fw and a table of discrete values for Fm are used in computing the Riemann sum. Thus the only source of error in computation of Ff is from discretization of Fm and using upper and lower Riemann sums we can find FfU ðcÞ and FfL ðcÞ upper and lower bounds for Ff ðcÞ: In Table 2, the sample values of Ff the mid-point Riemann sum (with Dz ¼ 0:0994814) are compared with Monte Carlo results.
5. Conclusions Fig. 9. Pulverizer tree for Example 6.
The research presented here shows that it is possible to
K.B. Chilakamarri / Probabilistic Engineering Mechanics 17 (2002) 317–325
evaluate reliability with arbitrary accuracy for those limitstate functions with closed form and in most cases this can be done efficiently. The algorithm presented here allows for parallel computations. The accuracy of results is purely a reflection of the accuracy of numerical integration and hence the error can be decreased arbitrarily without increasing the cost exponentially. The need to integrate high dimensional integrals increases with decreasing pulverizability. Even in such cases conditional probability concepts can be used to reduce the computational complexity.
Acknowledgments A preliminary version of this research was presented at NASA Glenn Research Center. I express my thanks to Dr Shantaram S. Pai at NASA GRC for his helpful
325
comments. I express my special thanks to the anonymous reviewers for their valuable and insightful comments.
References [1] Rosenblatt M. Remarks on a multivariate transformation. Ann Math Stat 1952;23:470–2. [2] Hohenbichler M, Rackwitz R. Nonnormal dependent vectors in structural reliability. J Engng Mech Div,ASCE 1981;107:1127–238. [3] Madsen HO, Krenk S, Lind NC. Methods of structural safety. Englewood Cliffs, NJ: Prentice-Hall; 1986. [4] Ditlevsen O, Madsen HO. Structural reliability methods. New York: Wiley; 1996. [5] Der Kiureghian A, Dakessian T. Multiple design points in first and second-order reliability. Struct Safety 1998;20:37– 49. [6] Maymon G. Some engineering applications in random vibration and random structures. Prog Astronaut Aeronaut 1998;1. [7] Tvedt L. Distribution of quadratic forms in normal space—application to structural reliability. J Engng Mech 1990;116(6):1183–97.
A new method in static structural reliability Kiran B. Chilakamarri1 CK & Associates, 1393 Tabor Avenue, Dayton, OH 45420, USA
Abstract A new method to determine the static structural reliability is presented in this paper. This method takes advantage of the fact that isoprobabilistic transformation, in particular the Rosenblatt’s transform, maps infinite Euclidean space into a unit hypercube and in doing so, reduces the problem of computing probability content of infinite region to a problem of computing a finite volume. The method is applied to several examples including a well-known example with multiple design points. q 2002 Published by Elsevier Science Ltd. Keywords: Reliability; Failure probability; Structures
1. Introduction In structural reliability the safety of a structure is the probability content of a set of safe points determined by a performance function in the parameter space representing the uncertain physical parameters of the structure. The performance function or the limit-state function may or may not have a closed form. Here we assume a limit-state function with a closed form. A fundamental problem in evaluating the probability of failure is the computation of the resulting multi-fold integral. The number of computations increases exponentially with the number of variables. Very frequently an isoprobabilistic transformation or probability transformation is used to standardize the model to convert the joint pdf of the parameters to a joint pdf of independent standard normal variables. A key observation made here is that the probability transformation T (Rosenblatt’s transform in particular) with k random variable maps R k (k-dimensional real space) into a unit hypercube Qk in such a way that for any measurable set S in R k, the volume of the image set TðSÞ is the probability content of S. A simple consequence of this observation is that it is possible to evaluate Pf (the probability of failure) directly in Qk and avoid integration over infinite sets. If in addition the limit-state function admits separation of variables into non-interacting subgroups then the problem of computing multi-fold integral can be split into several sub-problems of computing lower dimensional integrals. The term ‘pulverizer’ is used here to indicate this process. If E-mail address: [email protected] (K.B. Chilakamarri). Now a Visiting Associate Professor in the Department of Mathematics and Statistics, Wright State University, Dayton, OH 45435, USA.
a limit-state function can be pulverized perfectly then the multi-fold integral can be evaluated by calculating k 2 1 one-dimensional integrals where k is the number of random variables. On the basis of these observations an algorithm to evaluate Pf and cdf of performance function is developed. In Section 2, Rosenblatt’s transformation [1] is discussed and a formula for Pf is developed. In the same section the formula is applied to two examples, the first example is a wellknown problem with multiple design points and the second example is an artificial example with highly non-linear limit-state function. In Section 3, the algorithm is described and applied to a simple problem with three variables. In Section 4, the method is applied to three problems each with more than three variables. Example 4 is a well-known problem of one plastic collapse mechanism of one-bay frame. In Example 5, we find the cdf of maximum radial stress of rotating disk and in Example 6 we find the cdf of the natural frequency of a plate. Throughout this paper the results are compared with the existing Monte Carlo simulation results whenever possible.
2. The probability transformation Let ðX1 ; X2 ; …; Xk Þ be a random vector with distribution function Fðx1 ; x2 ; …; xk Þ ¼ Prob½X1 # x1 ; X2 # x2 ; …; Xk # xk : Let Qk denotes the k-dimensional unit hypercube, i.e. Qk ¼ {ðz1 ; z2 ; …; zk Þ;
0 # z1 # 1; …; 0 # zk # 1}:
1
0266-8920/02/$ - see front matter q 2002 Published by Elsevier Science Ltd. PII: S 0 2 6 6 - 8 9 2 0 ( 0 2 ) 0 0 0 1 5 - 2
The probability transformation T: Rk ! Qk also known as
318
K.B. Chilakamarri / Probabilistic Engineering Mechanics 17 (2002) 317–325
Rosenblatt transformation is defined as follows
Carlo simulation as reported by Madsen et al. [3] has a value of 0.00294.
Tðx1 ; x2 ; …; xk Þ ¼ z ¼ ðz1 ; z2 ; …; zk Þ
Example 1. Limit-state function is gðx1 ; x2 Þ ¼ 18 2 3x1 2 2x2 and the joint cumulative function is given by
where z1 ¼ F1 ðx1 Þ ¼ Prob½X1 # x1 ; z2 ¼ F2 ðx2 lx1 Þ ¼ Prob½X2
Fððx1 ; x2 Þ ¼ 1 2 expð2x1 Þ 2 expð2x2 Þ
# x2 lX1 ¼ x1 ; z3 ¼ F3 ðx3 lx2 ; x1 Þ ¼ Prob½X3 # x3 lX1
þ expð2ðx1 þ x2 þ x1 x2 ÞÞ;
¼ x1 ; X2 ¼ x2 ; …; zk ¼ Fk ðxk lxk21 ; …; x1 Þ ¼ Prob½Xk
x1 ; x2 . 0:
ð1Þ
As described by Hohenbichler and Rackwitz [2], there are 2! ¼ 2 distinct transformations. Note that the random variables are not independent for this example.
The random variables Z1 ; Z2 ; …; Zk are uniformly and independently distributed over ½0; 1 and the random vector ðZ1 ; Z2 ; …; Zk Þ is uniformly distributed over the k-dimensional unit hypercube Qk [1]. If in addition, the random variables X1 ; X2 ; …; Xk are independent then
Transformation I: Tðx1 ; x2 Þ ¼ ðz1 ; z2 Þ where, z1 ¼ 1 2 e2x1 and z2 ¼ 1 2 ð1 þ x2 Þ e2x2 ð1þx1 Þ : The image of the curve g ¼ 0 in the Z-square (a unit square) is given by
ðz1 ; z2 ; …; zk Þ ¼ Tðx1 ; x2 ; …; xk Þ
z2 ¼ 1 2 ½10 þ ð3=2Þ lnð1 2 zÞ exp½2{9 þ ð3=2Þ
# xk lX1 ¼ x1 ; X2 ¼ x2 ; …; Xk21 ¼ xk21 :
lnð1 2 zÞ}{1 2 lnð1 2 zÞ}
¼ ðFX1 ðx1 Þ; FX2 ðx2 Þ; …; FXk ðxk ÞÞ: Note that different permutations of X1 ; X2 ; …; Xk can produce different transformations. For the method proposed in paper it is necessary to be able to invert the Rosenblatt’s transformation. Given z1 ; z2 ; …; zk ; we first find x1, after x1 is evaluated both z2 and x1 together yield x2, after x2 is evaluated z3, x1 and x2 can be used to find x3, and so on. The reader is urged to note that the independence of random variables is not required for this process. We assume the set S of failure points is characterized by non-negative values of function gðxÞ the limit-state function. We may define a corresponding failure set Sp ¼ TðSÞ in Qk. Accordingly, the failure surface ›S, the boundary of S, is mapped onto the boundary ›S p. Thus we can calculate the probability of failure either in R k or in the hypercube Qk. Since z is uniformly distributed in Qk ð ð Pf ¼ dFðxÞ ¼ dz ¼ VolumesðSp Þ: ð2Þ Sp
S
Once the boundary ›S p is identified in Qk the problem of evaluating Pf is reduced to a problem of evaluating the volume of a finite domain in Z-space. If the failure region can be characterized by xk , hðx1 ; x2 ; …; xk21 Þ and Tðx1 ; x2 ; …; xk21 Þ denotes the restriction of the map T to the first k 2 1 coordinates then Pf ¼
ð1 ð1 0
0
where z denotes z1 and the failure region in Z-square corresponding to g , 0 is the region above this curve. Thus Pf ¼
ð3Þ Example 1 with multiple design points is a well-known problem in FORM and SORM and originally introduced by Hohenbichler and Rackwitz [2] and discussed by Madsen et al. [3], Ditlevsen and Madsen [4] and Der Kiureghian and Dakessian [5]. The exact solution of the problem by Monte
½10 þ ð3=2Þ lnð1 2 zÞ exp½2{9 þ ð3=2Þ
0
lnð1 2 zÞ}{1 2 lnð1 2 zÞ} dz þ e26 ¼ 0:00294485974: The value of Pf is accurate up to 10 digits. To see how this integral is written, note that we are interested in finding the area of the image of the region g , 0 under the probability transformation in Z-square. For a given z1 ð0 , z1 , 1Þ; we first find the corresponding x1 value from z1 ¼ 1 2 e2x1 : Having found x1, the x2 value is calculated from g ¼ 0: Once x2 is evaluated we find z2 corresponding to z1 from z2 ¼ 1 2 ð1 þ x2 Þ expð2x2 ð1 þ x1 ÞÞ: Note that the support of the image of g ¼ 0 is the interval ½0; 1 2 e26 Þ and the area of interest is the area above the image curve. Transformation II: Tðx1 ; x2 Þ ¼ ðz1 ; z2 Þ; where z1 ¼ 1 2 e2x2 and z2 ¼ 1 2 ð1 þ x1 Þ e2x1 ð1þx2 Þ : Image of the curve g ¼ 0 in Z-square is given by z2 ¼ 1 2 ½7 þ ð2=3Þ lnð1 2 zÞ exp½2{6 þ ð2=3Þ lnð1 2 zÞ}{1 2 lnð1 2 zÞ}
ð1 · · · T h T 21 ðz1 ; z2 ; …; zk21 Þ dz1 dz2 · · ·dzk21 : 0
ð12e26
and again Pf ¼
ð12e29
½7 þ ð2=3Þ lnð1 2 zÞexp½2{6 þ ð2=3Þ
0
lnð1 2 zÞ}{1 2 lnð1 2 zÞ} dz þ e29 ¼ 0:00294485974:
K.B. Chilakamarri / Probabilistic Engineering Mechanics 17 (2002) 317–325
Example 2. Let X and Y be independent standard normal variables. Let the failure region be given by
admit such decomposition. For instance, if
4 ; y$ x{1 þ sin2 ðxÞ}
g¼ x . 0:
Probability transformation is given by ðz1 ; z2 Þ ¼ Tðx; yÞ ¼ ðFðxÞ; FðyÞÞ: The g-function is gðx; yÞ ¼ y 2
4 x{1 þ sin2 ðxÞ}
and Pf ¼ Prob½g $ 0: The failure curve is described by yðxÞ ¼
4 ; x{1 þ sin2 ðxÞ}
then Pf ¼
ð1 0
12F
4 dz: F21 ðzÞð1 þ sin2 ðF21 ðzÞÞÞ
To see why, for a point z in the unit interval, F21 ðzÞ is the corresponding value of x, the y value from the failure curve is yðF21 ðzÞÞ and the area above the curve is Pf. The integral for Pf can be approximated by a finite integral using ^ 5 sigma bounds, further the region of interest can be restricted to {ðx; yÞ : x [ ½0:479177; 5 and y [ ½0; 5} where xp ¼ 0:479177 is the x-value when y ¼ 5 on the curve gðx; yÞ ¼ 0:Thus ðFð5Þ 4 Pf , 12F dz F21 ðzÞð1 þ sin2 ðF21 ðzÞÞÞ Fðxp Þ ¼ 0:0121736: This result is accurate up to six digits. For more than two variables we develop an algorithm described in Section 3.
3. The algorithm We need to introduce new terms to describe the algorithm. Given an expression gðX1 ; X2 ; …; Xk Þ in k variables divide the set of variables {X1 ; X2 ; …; Xk } into two subgroups say {X1 ; X2 ; …; Xi } and {Xiþ1 ; Xiþ2 ; …; Xk } so that g can be expressed as combination of two expressions g 1 and g2 , where g1 involves only X1 ; X2 ; …; Xi and g2 involves only Xiþ1 ; Xiþ2 ; …; Xk : Simply we want to find two expressions g1 ðX1 ; …; Xi Þ and g2 ðXiþ1 ; …; Xk Þ so that gðX1 ; X2 ; …; Xk Þ ¼ Fðg1 ; g2 Þ for some function f. For instance if g ¼ X1 ðX2 þ X3 Þ; then g1 ¼ X1 and g2 ¼ X2 þ X3 : If g ¼ ðX1 þ 3ÞX2 X32 ðX42 2 X52 Þ many choices of g1 and g2 are possible. Some expressions may not
319
X1 ; X2 þ X3 ðX1 þ 10Þ
then no two variables can be separated from the third. Repeat the process with g1 and g2, and so on until the process is terminated. I call this method pulverizer. The end product of the pulverizer method is a record in the form of a Pulverizer Tree. This is a rooted tree in which the root is g, has two children g1 and g2, and the nodes g1 and g2 contain their children and so on (see Fig. 1, 3, 7 and 9). The pulverizer tree provides a quick reference for the reduction of g into smaller expressions. Several pulverizer trees may exist for a given expression. The number of integrals to be evaluated is exactly the number of nodes with degree more than one. Next we introduce the concept of Working Area. Given a real number 1 . 0 and the joint distribution F of k random variables the working area in R k is defined as a region W with Prob½W $ 1 2 1: Similarly, the corresponding working area in Z-space (the hypercube) is TðWÞ ¼ W 0 with VolumeðW 0 Þ $ 1 2 1: As the name indicates, when solving a problem we do not venture outside the working area, so the probability content of the complement of W makes a contribution of at most 1 to the error. For example with two independent normal distributions, X , Nð10; 4Þ and Y , Nð2; 0:1Þ; using ^ 3 sigma bounds, the working area W in R 2 is given by 22 # x # 22 and 1:7 # y # 2:3; and in Zsquare the working area is simply ½Fð23Þ; Fð3Þ £ ½Fð23Þ; Fð3Þ: The error from the ignorance of outside of W is exactly 1 2 ½Fð3Þ 2 Fð23Þ2 : If we wish to use ^ 20 sigma bounds then similar analysis produces a working area with associated error ¼ 1 2 ½Fð20Þ 2 Fð220Þ2 : For an additional example consider two independent variables, X , Nð10; 4Þ and Y , Uð10; 121Þ a uniform distribution. Suppose again for X we use ^ 3 sigma bound, then the working W area in R 2 is, 22 # x # 22 and 10 # y # 121: In Z-square, the working area is ½Fð23Þ; Fð3Þ £ ð0; 1Þ: The error associated with W is 1 2 ½Fð3Þ 2 Fð23Þ: Note that the working area can be extended easily since the concept is additive. Also note that any information about the outside of working area may be used to reduce this type of error. The algorithm I. II.
Write the probability transformation. Decide on a suitable pulverizer tree for the given limitstate function. (Here we want to avoid nodes of high
Fig. 1. Pulverizer tree for Example 3.
320
K.B. Chilakamarri / Probabilistic Engineering Mechanics 17 (2002) 317–325
degree. We also want to choose a short and broad tree). III. To find the cdf of a node labeled hðu; vÞ with a pair of nodes u and v attached, in the working area study the level curves of hðu; vÞ ¼ c (or u ¼ Fðu; cÞ) to write an expression for Fh ðcÞ ¼ Prob½u , Fðu; cÞ using the formula ð1 ð4Þ Fh ðcÞ ¼ Fu f Fv21 ðzÞ; c 0
Any integration technique can be used so long as it also provides error estimates, i.e. functions UðcÞ and LðcÞ so that LðcÞ # Fh ðcÞ # UðcÞ for all c. IV. Start at the bottommost nodes and keep evaluating the cdf of next higher level nodes until we reach the top.
Remark 1. There are two sources of error in the method presented here. As already mentioned, restricting the working area results in a predetermined amount of error. Another source of error is from numerical integration. However, both types of errors can be controlled. Remark 2. Suppose the number of variables is k. If a pulverizer tree has root with degree d p and has d3 number of vertices of degree 3, d4 number of vertices of degree 4,…, and dR number of vertices of degree R. The number of vertices with degree 1 is k, same as the number of variables. Using elementary counting techniques we can show that k ¼ dp þ d3 þ 2d4 þ · · · þ ðR 2 2ÞdR : A non-root vertex of degree iði . 2Þ requires a ði 2 2Þdimensional integration to compute one value of the cdf associated with that vertex and the root requires a ðdp 2 1Þdimensional integration. If N is the grid size, then the total number of computations required is of the order p O d3 N þ d4 N 2 þ · · · þ dR N R22 þ N d 21 : In general if a limit-state function of k variable admits complete pulverization, i.e. d4 ¼ d5 ¼ · · · ¼ 0; then d3 ¼ k 2 2 and dp ¼ 2: Thus Oððk 2 1ÞNÞ computations are required. To reduce the number of computations, trees with small degrees are preferred. Computational errors would compound along a long branch, thus short pulverizer trees are preferred.
EðSÞ ¼ 2500 kg/cm 2, EðAÞ ¼ 2 cm 2, EðFÞ ¼ 4166 kg, sðSÞ ¼ 75 kg/cm2, sðAÞ ¼ 0:04 cm2, and sðFÞ ¼ 125 kg. Failure occurs if S 2 ðF=AÞ # 0: Let A ¼ X1 , Nð2; 0:04Þ; S ¼ X2 , Nð2500; 75Þ and F ¼ X3 , Nð4166; 125Þ and the limit-state gðX1 ; X2 ; X3 Þ ¼ X2 2 ðX3 =X1 Þ: Probability transformation is given by x 22 ðZ1 ; Z2 ; Z3 Þ ¼ FX1 ; FX2 ; FX3 ¼ F 1 ; 0:04
F
x2 2 2500 ; 75
x 2 4166 F 3 125
Writing m ¼ X3 =X1 and g ¼ X2 2 m the pulverizer tree is shown in Fig. 1. First we find Fm, the cdf of m, from FX1 and FX3 and then evaluate the cdf of g from Fm andFX2 : It is not hard to show that if x and y are independent standard normal variables, then Ay þ B ,c Prob ax þ b is the probability content of the region R in R 2 shown in Fig. 2. This probability is exactly equal to the area of the region R in the Z-square. Using this result in our context, we have
125y þ 4166 ,d ; Fm ðdÞ ¼ Prob½ðX3 =X1 Þ , d ¼ Prob 0:04x þ 2 where y and x are standard normal variables. For this problem, ð2b=aÞ ¼ 250; and Fð250Þ , 0 and y ¼ ðad=AÞx þ ðdb 2 BÞ=A ¼ 0:00032 dx þ ð2d 2 4166Þ=125; thus we have ð1 2d 2 4166 21 Fm ðdÞ ¼ F 0:00032dF ðzÞ þ dz 125 0 and using ^ 5 sigma bounds the relevant range of d is
4. Examples We begin with a simple example with three variables [6]. Example 3. Consider a bar of cross-section A and a yield stress S is subjected to a tensile force F where A, S and F are independently distributed normal random variables with
Fig. 2. The area of the marked region in Z-square is the probability content of the region R in x –y plane.
K.B. Chilakamarri / Probabilistic Engineering Mechanics 17 (2002) 317–325
321
below the diagonal of the square working area ð144 # c # 264Þ and (ii) those that are above the diagonal ð269 # c # 394Þ: We can now write the integral for the cdf of m Fm ðcÞ ¼
¼12
ðFX1 ðxp Þ 0
ð1 FX1 ðxp Þ
Fig. 3. Pulverizer tree for Example 4.
½1609; 2661:69:We are now ready to find Fg where g ¼ X2 2 m and FX21 ðzÞ ¼ 75F21 ðzÞ þ 2500 2 Fg ðcÞ ¼ Prob½g # c ¼ Prob½m $ X2 2 c ¼
ð1 0
¼
ð1 0
1 2 Fm FX21 ðzÞ 2 c dz 2
1 2 Fm ½75F21 ðzÞ þ 2500 2 c dz
¼ ð5:7143Þ £ 1025 The FORM analysis [6] yields a value of (5.583) £ 1025. Example 4 is a well-known example [3,7]. Example 4. The failure region of one plastic collapse mechanism of one-bay frame is given by X1 þ 2X2 þ 2X3 þ X4 2 5X5 2 5X6 , 0: The variables X1 ; …; X6 are statistically independent and log-normally distributed with mX1 ¼ mX2 ¼ mX3 ¼ mX4 ¼ 120; mX5 ¼ 50; mX6 ¼ 40 and sX1 ¼ sX2 ¼ sX3 ¼ sX4 ¼ 12; sX5 ¼ 15; and sX6 ¼ 12: The probability transformation is described below ln t 2 4:78252 FX1 ðtÞ ¼ FX2 ðtÞ ¼ FX3 ðtÞ ¼ FX4 ðtÞ ¼ F ; 0:0997513 ðzÞ ¼ FX21 ðzÞ ¼ FX21 ðzÞ ¼ FX21 ðzÞ FX21 1 1 3 4 21
FX 5
¼ expð4:78252 þ 0:0997513F ðzÞÞ; ln x5 2 3:86893 ¼F ; 0:29356
ðzÞ ¼ expð3:86893 þ 0:29356F21 ðzÞÞ; FX21 5 ln x5 2 3:64579 FX 6 ¼ F ; 0:29356
FX4 c 2 FX21 ðzÞ dz 1
for 144 # c # 264
1 2 FX4 c 2 FX21 ðzÞ dz for 264 # c # 394 1
where xp ¼ c 2 72 and xp ¼ c 2 197: Next step involves evaluating F› but Fm ¼ F› since X1 ; X2 X3 ; X4 are all identically distributed. Now we can evaluate Fb the cdf of b ¼ m þ 2›: For the computational efficiency a table of values of Fm is computed, stored and used to approximate the integrals that appear in the formula for Fb. The significant domain for Fm ðcÞ turns out to be a smaller interval ½170; 300 than that given by working area bounds. The interval ½170; 300 is partitioned so that variation of Fm is sufficiently small. Here, a partition 170 ¼ c0 , c1 , · · · , c300 ¼ 340 is used so that lFm ðci Þ 2 Fm ðciþ1 Þl , 0:005192: Once again using ^ 5 sigma bounds the working area for b ¼ m þ 2› in the m – › plane is ½170; 340 £ ½170; 340 and the level curves are shown in Fig. 4. The integral for Fb is easy to write Fb ðcÞ ¼ 8 ! ðFm ðc2340Þ > c 2 Fm21 ðzÞ > > dz Fm > > 2 > 0 > > > ! > < ð1 c 2 Fm21 ðzÞ Fm dz > 2 0 > > > ! > > ð1 > c 2 Fm21 ðzÞ > > dz > : 1 2 F ðc2680Þ 1 2 Fm 2
if 510 # c # 680
if 680 # c # 850
if 850 # c # 1020
m
Once again we construct a table of values for Fb for computational purposes. Here Fb values were stored as {ðci ; zi Þ : 1 # i # 355} with Fb ðci Þ ¼ zi with Maxlziþ1 2 zi l , 0:00538458: Next, we evaluate FD the cdf of D ¼ x5 þ x6 from the distributions of X5 and X6. Using ^ 5 sigma bounds the working area in the x6 2 x5 plane is restricted to
ðzÞ ¼ expð3:64579 þ 0:29356F21 ðzÞÞ: FX21 6 Writing m ¼ X1 þ X4 ; › ¼ X2 þ X3 ; b ¼ m þ 2›; D ¼ X5 þ X6 ; and g ¼ b 2 5D; we choose a pulverizer tree as shown in Fig. 3. We first evaluate Fm ðcÞ the cdf of m. Assuming ^ 5 sigma bounds, the working area in x1 2 x4 plane is a square ½72; 197 £ ½72; 197; and the level curves of m (i.e. x1 þ x4 ¼ c) are straight lines for 144 # c # 394: We distinguish between the two kinds of lines (i) those that are
Fig. 4. Level curves of b in the working area.
322
K.B. Chilakamarri / Probabilistic Engineering Mechanics 17 (2002) 317–325
Examples 5 and 6 are from the article ‘Accuracy of probabilistic methods for HCF’ by Eric Fox in the Proceedings of the 5th National Turbine Engine High Cycle Fatigue Conference HCF 2000, 7 –9 March 2000, Chandler, Arizona. We evaluate the cdf of the expressions involved in each example and compare the results with Monte Carlo results from the above mentioned article.
Fig. 5. Level curves of D in the working area.
½8; 167 £ ½11; 208: The level curves of D are shown in Fig. 5. Since FD ðcÞ ¼ Prob½D # c ¼ Prob½x5 þ x6 # c; we need to find the area in Z-square of the region under the image curve of x5 # c 2 x6 : FD ðcÞ ¼ 8 ðF ðc211Þ X6 > 21 > F c 2 F ðzÞ dz > X5 X6 > > 0 > > > < ð1 FX5 c 2 FX21 ðzÞ dz 6 > 0 > > > > ð1 > > 21 > F c 2 F ðzÞ dz : X5 X6
for 19 # c # 178 for 178 # c # 216 for 216 # c # 375
FX6 ðc2208Þ
where the parameters are assumed to be independent and are distributed as follows: n ¼ Poisson’s ratio ¼ X1 , Uð0:29; 0:31Þ; r ¼ density ¼ X2 , Uð0:279; 0:289Þ; v ¼ rotor speed ¼ X3 ¼, Uð10 000; 11 000Þ; ro ¼ outer disk radius ¼ X4 , Nð8; 0:02Þ and ri ¼ inner disk radius ¼ X5 , Nð2; 0:01Þ where Uða; bÞ denotes uniform distribution over the interval ½a; b: Writing m ¼ ð3 þ X1 ÞX2 ; f ¼ mX32 ; D ¼ X42 2 X52 and g ¼ c0 f D where c0 is a constant comprising of all constants involved in s, we choose the pulverizer tree shown in Fig. 7. The probability transformation is T ¼ ðz1 ; z2 ; …; z5 Þ ¼ ðFX1 ; FX2 ; …; FX5 Þ where
To complete the problem we find Pf ¼ Fg ð0Þ ¼ Fg ðb , 5DÞ ¼ 1 2
Example 5. The maximum radial stress of a rotating disk, s, in psi, is given by 3þn r 2p 2 2 s¼ v ro 2 ri2 ; 60 8 ð9:81Þð39:37Þ
ð1 0
!
FD
Fb21 ðzÞ dz 5
The graph of g ¼ 0; i.e. b 2 5D ¼ 0 and the region g # 0 are shown in Fig. 6. Mid-point Reimann sum is used to approximate the integral and we find Pf ¼ 0:012224 and from the lower and upper Riemann sums we find the lower and upper bounds PLf ¼ 0:012136; PU f ¼ 0:0123135; respectively. Using more refined partition for Fb with Maxlziþ1 2 zi l , 0:00222917; we find Pf ¼ 0:0122188 and the bounds PLf ¼ 0:0121782 and PU f ¼ 0:0122599:
FX21 ðzÞ ¼ 1 ¼
z z þ 0:29; FX21 þ 0:279; FX21 ðzÞ ¼ ðzÞ 2 3 50 100 z þ 10 000; FX21 ðzÞ 4 0:001
¼ 0:02F21 ðzÞ þ 8; FX21 ðzÞ ¼ 0:01F21 ðzÞ þ 2 5 for 0 , z , 1: We can write the integral for Fm ðcÞ ! ð1 c Fm ðcÞ ¼ FX2 dz: 3 þ FX21 ðzÞ 0 1 Since FX2 is not continuous, some caution must be exercised in evaluating the above integral. A study of the level curves of m ¼ ð3 þ X1 ÞX2 help in writing the integrand explicitly
Fig. 6. The failure region g , 0:
Fig. 7. Pulverizer tree for Example 5.
K.B. Chilakamarri / Probabilistic Engineering Mechanics 17 (2002) 317–325
and the integral can be evaluated exactly 8 0 > > > > > 4589:55 2 3:0514:5c þ 5000c ln½179:211c > > < Fm ðcÞ ¼ 227:9 þ 30:3031c > > > > 24781:95 þ 30 544:8c 2 5000c ln½173:01c > > > : 1
323
if c , 0:91791 if 0:91791 , c , 0:923492 if 0:923492 , c , 0:95081 if 0:95081 , c , 0:95659 if c $ 0:95659
Next, Ff the cdf of f ¼ mX32 is obtained from the integral 0 1 ð1 c B C Ff ðcÞ ¼ Fm @ 2 Adz: 21 0 FX3 ðzÞ
image of the area under the curve of type I and FD ðcÞ is obtained by integrating (rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi) 2 FX21 ðzÞ þc FX 4 5
Once again a quick study of the level curves of f in the pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi working area helps in writing the explicit integrals for Ff as from z ¼ FX21 ð ð7:906Þ2 2 cÞ to z ¼ 1: 5 follows Now we can complete the problem by finding Fs the cdf of 8 0 if c , 9:1791 £ 107 > > > p > ðv > c > > Fm dz if 9:1791 £ 107 , c , 9:5659 £ 107 > > 2 > ð1000z þ 10 000Þ 0 > > > < ðvp c Ff ðcÞ ¼ F dz þ up if 9:5659 £ 107 , c , 1:11067 £ 107 m 2 > p ð1000z þ 10 000Þ u > > > 2 > > ð1 > c > > F dz þ up if 1:11067 £ 107 , c , 1:115477 £ 108 > > up m ð1000z þ 10 000Þ2 > > : 1 otherwise where vp ¼ 0:001
qffiffiffiffiffiffiffiffiffiffiffiffiffi c=0:91791 2 10 000 and
qffiffiffiffiffiffiffiffiffiffiffiffiffi u ¼ 0:001 c=0:95659 2 10 000 p
These integrals can be evaluated to provide an exact formula for Ff. Next we evaluate the FD the cdf of D ¼ X42 2 X52 : Using ^ 10 sigma bounds we can find the working area, a rectangle, in the x5 2 x4 plane as shown in Fig. 8. For a given c qffiffiffiffiffiffiffiffiffi D , c $ x24 , x25 þ c $ x4 , x25 þ c:
s ¼ g ¼ c0 f D from the following integral ! ð1 c dz: P½g , c ¼ Fg ðcÞ ¼ Ff c0 FD21 ðzÞ 0 Once again a study of the level curves of g in f 2 D space provides a valuable service in restricting the range for numerical integration. Riemann sums were used to approximate the integral and for that purpose a table for FD was developed so that the Max½Dz , 0:00399814 in evaluation of Fg. For Ff the exact formula was used hence there are no errors from Ff. A sample of results of this analysis and the corresponding Monte Carlo simulation results from the referred article are shown in Table 1. Example 6. The natural frequency of a plate, f, in Hz is
The cdf FD ðcÞ is calculated by integrating (rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi) 2 FX 4 FX21 ðzÞ þc 5
over the interval ð0; 1Þ: We notice that the level curves of D are of three types and accordingly some care must be taken for numerical integration process. For instance if 58:3136 , c , 58:6909; then the area in Z-square corresponds to the
Fig. 8. Level curves of D.
324
K.B. Chilakamarri / Probabilistic Engineering Mechanics 17 (2002) 317–325
Table 1 The cdf of radial stress s of rotating disk
Table 2 The cdf of the natural frequency f of a plate
Percentile
Current new method s
Monte Carlo simulation s
Percentile
Current new method f
Monte Carlo simulation f
0.001 0.01 0.1 0.5 0.9 0.99 0.999
19 19 20 22 23 24 24
19 19 20 22 23 24 24
0.001 0.01 0.1 0.5 0.9 0.99 0.999
1751.4 1874.78 2050.05 2274.12 2508.51 2707.46 2858.2
1752.6 1876.9 2049.8 2274.2 2507.1 2706.3 2858.7
598.5 814 355 002 715 336 594
595.5 810.5 349.8 005.4 714.6 330.9 601
calculated as 0 f ¼
34:03 B @ 2pa2 12
3
11=2
Eh C A g 2 ð1 2 n Þ 386:1
where the parameters are assumed to be independent and are distributed as follows, a ¼ length of plate ¼ X1 , N ð12; 0:1Þ; E ¼ modulus of elasticity ¼ X2 , Nð30 £ 106 , 9 £ 105 Þ; h ¼ plate thickness ¼ X3 , Nð1; 0:05Þ; g ¼ mass per unit area ¼ X4 , Nð0:29; 0:004Þ and n ¼ Poisson’s ratio ¼ X5 , Uð0:29; 0:31Þ; where Uða; bÞ denotes uniform distribution over the interval ½a; b: Writing c ¼ ðX33 =X14 Þ; m ¼ cX2 and w ¼ X4 ð1 2 X52 Þ and combining all the constants appearing in the expression for the frequency, we can write f ¼ 30:7214
rffiffiffiffi m w
Next we need to evaluate ! ð1 d dz; Fm ðdÞ ¼ FX2 Fc21 ðzÞ 0 where FX2 ðtÞ ¼ Fððt 2 30 £ 106 Þ=9 £ 105 Þ: Here again we use ^ 5 sigma bounds. The range of d is the interval ð637:5; 3243Þ: To find the cdf of w we use ^ 5 sigma bounds for X4 the only normal variable appearing in w. The working area is {ðx5 ; x4 Þ : 0:29 # x5 # 0:31 and 0:27 # x4 # 0:31}: A quick study of level curves of w shows that the support of the integrand in the formula for the cdf of w, changes between the intervals ð0:244053; 0:247293Þ; ð0:247293, 0:280209Þ; and ð0:280209; 0:283929Þ. From this information it is not hard to evaluate Fw ðcÞ; where 0 1 c ð1 B 20:29 þ 1 2 ððz=50Þ þ 0:29Þ2 C Cdz: Fw ðcÞ ¼ wB @ A 0:004 0 Finally we can write
and a pulverizer tree chosen for this problem shown in Fig. 9. Using ^ 5 sigma bounds Fc ðcÞ ¼
ð1 0
FX 3
! rffiffiffiffiffiffiffiffiffiffiffiffiffi 4ffi 3 c FX21 ðzÞ dz 1
p3 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! ðFð5Þ c½0:1F21 ðzÞ þ 124 2 1 dz ¼ F 0:05 Fð25Þ ð0:0001728 # c # 0:000111671Þ:
ð30:7214Þ2 m Ff ðcÞ ¼ Prob½F , c ¼ Prob w . c2
¼
ð1 0
(
1 2 Fw
30:7214 c
2
!
!) Fm21 ðzÞ
dz
Again, the nature of the support of the integrand changes at c ¼ 1440:22; 1553.43, 3132.52, and 3378. The values of Ff are computed using Riemann sum approximation for the integral. The integral formula for Fw and a table of discrete values for Fm are used in computing the Riemann sum. Thus the only source of error in computation of Ff is from discretization of Fm and using upper and lower Riemann sums we can find FfU ðcÞ and FfL ðcÞ upper and lower bounds for Ff ðcÞ: In Table 2, the sample values of Ff the mid-point Riemann sum (with Dz ¼ 0:0994814) are compared with Monte Carlo results.
5. Conclusions Fig. 9. Pulverizer tree for Example 6.
The research presented here shows that it is possible to
K.B. Chilakamarri / Probabilistic Engineering Mechanics 17 (2002) 317–325
evaluate reliability with arbitrary accuracy for those limitstate functions with closed form and in most cases this can be done efficiently. The algorithm presented here allows for parallel computations. The accuracy of results is purely a reflection of the accuracy of numerical integration and hence the error can be decreased arbitrarily without increasing the cost exponentially. The need to integrate high dimensional integrals increases with decreasing pulverizability. Even in such cases conditional probability concepts can be used to reduce the computational complexity.
Acknowledgments A preliminary version of this research was presented at NASA Glenn Research Center. I express my thanks to Dr Shantaram S. Pai at NASA GRC for his helpful
325
comments. I express my special thanks to the anonymous reviewers for their valuable and insightful comments.
References [1] Rosenblatt M. Remarks on a multivariate transformation. Ann Math Stat 1952;23:470–2. [2] Hohenbichler M, Rackwitz R. Nonnormal dependent vectors in structural reliability. J Engng Mech Div,ASCE 1981;107:1127–238. [3] Madsen HO, Krenk S, Lind NC. Methods of structural safety. Englewood Cliffs, NJ: Prentice-Hall; 1986. [4] Ditlevsen O, Madsen HO. Structural reliability methods. New York: Wiley; 1996. [5] Der Kiureghian A, Dakessian T. Multiple design points in first and second-order reliability. Struct Safety 1998;20:37– 49. [6] Maymon G. Some engineering applications in random vibration and random structures. Prog Astronaut Aeronaut 1998;1. [7] Tvedt L. Distribution of quadratic forms in normal space—application to structural reliability. J Engng Mech 1990;116(6):1183–97.