- Email: [email protected]
Monte carlo finite element method of structure reliability analysis
Reliability Engineering and System Safety 40 (1993) 77-83
Monte Carlo finite element method of structure reliability analysis Jin Guoliang, Chen Lilt & Dong Jiamei Daqing Petroleum Institute, Anda, Heilongjiang, People's Republic of China (Received 18 August 1991; accepted 10 July 1992)
To resolve the problem in reliability of any complicated engineering structures by use of general reliability methodology one often comes across a great deal of limitations due to realistic difficulties. In order to make use of the advantage that the finite element method could numerically resolve the load response of any structure, and the character that Monte Carlo method could carry out stochastic simulation, the Monte Carlo finite element method for any structure reliability analysis was presented. The advantage of high precision, less work in calculation and rapid convergence of the method developed above has been proved by concrete reliability calculation and analysis of the given engineering structure possessing analytic solutions (e.g. a thick wall cylinder). Using the method for rather complicated structures (e.g. high pressure cast steel casing head bearing of oil well head), it was discovered it is necessary to use only a few simulations to obtain the structure reliability value, which took further steps to prove the effectiveness of the Monte Carlo finite element method.
engineering problems, especially for complicated engineering structures. Therefore, many improved methods were put forward in order to avoid directly obtaining the distribution laws of r and s and the integral of the complex function. Among them, the first order second moment method and the design testing point method are only adapted for resolving a regular structure in which the random variables are suited to the normal distribution and the log-normal distribution, and meanwhile the limit state equation is easily set up. The JC method and the Paloheimo method could be used under the condition that the random variables obey any distribution, but both of them are only adapted for a regular structure where the limit state equation is easily set up with an analytic expression. The Monte Carlo method could be used under the condition that the random variable obeys any distribution, but it is necessary to obtain the expression for the limit state equation and it takes a lot of time to conduct the computing, and sometimes, it is difficult to put into practice because of lack of computer capacity. In a word, all of the methods mentioned above are not perfect for many complicated engineering structures, therefore, it is necessary to do research to think of a new way to analyse the reliability of complicated engineering structures.
1 INTRODUCTION The theory of interference is the basis of structural reliability analysis, it has been used widely for analysing structural reliability. If the structure reliability is assumed as the probability that the stress s in the structure control point is less than the material yield strength, that is, when r
R(t)=f_~fs(s)[~fr(r)dr]ds
(1)
where fr(r) is the probability density function of the material resistance fs(s) is the probability density function of the structural stress, R(t) is the structural reliability. Equation (1) could be used only for resolving the problems in which the distribution laws of both the material resistance r and the stress s of the structure being analysed are known in advance. But it is difficult to obtain the distribution laws of r and s for realistic
ReliabilityEngineeringand SystemSafety0951-8320/93/$06.00 © 1993 Elsevier Science Publishers Ltd, England. 77
78
Jin Guoliang, Chen Lin, Dong Jiamei
2 THE PRESENT SITUATION OF THE MONTE CARLO FINITE ELEMENT METHOD One of the important tasks of structural reliability analysis is to find the structural load response (stress). For complicated engineering structures, very often it is a unique feasible way of adopting a numerical method to calculate their load response. The finite element method has become an effective numerical method along with the rapid development of the technique of the electronic computer. The traditional finite element method is built upon the basis of determinacy of structural system which considers that the structural geometrical size and the load it is subjected to are definite values. But, in fact, the geometrical size of the structure, the load, the material resistance, etc., are usually random variables. Structural reliability analysis is used just to synthesize the indeterminacy of each variable to assess the bearing capacity of the structure. It will have an obvious significance if some appropriate improvement of a large number of effective finite element programs which have been presented at the moment could be completed. The Monte Carlo method is a powerful tool to be used to simulate random variables. It can simulate the random variables that are suitable for any distribution. However, it demands a large amount of calculation if a relative precise solution is needed, which will certainly mean a lot of work for processing the data for a realistic and complicated engineering structure. If some improvement in tackling the problem described above could be successfully made, it would be very useful for the reliability analysis of complicated engineering structures. During the whole process of study, special attention to the basic idea of combination of the Monte Carlo method with the finite element method was paid, and the research work was conducted repeatedly, over and over again, and the Monte Carlo finite element method of reliability analysis for any complicated engineering structures was finally put forward.
3 THE MONTE CARLO FINITE ELEMENT METHOD 3.1 Determination of all basic parameter distribution concerning structure reliability When using the Monte Carlo finite element method to resolve structural reliability, first of all, the basic parameter distribution law for the size, load, material resistance, etc., concerning the structure reliability should be given. According to statistical analysis, the distribution of structural size generally is a normal distribution. If we assume that the + A X and trx
express the tolerance of machine manufacturing and the standard deviation respectively, the relationship between both of them could be established as follows: 1
ax = AS~3
(2)
The distribution law of load could be obtained only from the statistical analysis of a large number of data collected from field structures. The distribution law of material resistance should be obtained from tests in principle. If there are no testing conditions, it can be quoted from technique document directly. Generally, r obeys a normal distribution,z it means,/~r, and standard deviation, or, can be obtained by looking up related statistical materials. 3-5 It is acceptable in engineering applications in this way to estimate the preliminary distribution of the resistance of metal materials.
3.2 Producing uniform distribution random data When using the Monte Carlo finite element method to analyse structural reliability, after the distribution law of each random variable is established, it is necessary to do random sampling for each variable, and meanwhile, to obtain the sampling values of the random variables. In order to obtain the random sampling values of the structural size and external load rapidly and with high precision, two steps may be taken: first to produce uniform distribution random data between zero and one, and then on this basis, transforming them into given variable random sampling. In our study, the multiplier common modulus method was adopted to produce uniform distribution random data between zero and one. The expression for the multiplier common modulus method is as follows:
x(i + 1) = (~. × x(i)) x (mod m)
(3)
r(i + 1) = x(i + 1)/m
(4)
where ~. and m are positive integers. The random data produced in this way are puppet random data, 5'6 if the value of )., x(0) and m were appropriately chosen, the puppet random data sequence can be freely passed through parameter testing, uniform testing and independence testing, then a puppet random data sequence with a long period, good uniformity and good independence could be obtained.
3.3 Sampling random variable of given distribution It will be meaningful that until we have to vary the uniform distribution random data between zero and one into any given distribution random variable and
Monte Carlo finite element method of structure reliability analysis
79
3.5 The establishment of the function of reliability
I
q
0
,
Fig. 1. The audio-visual diagram of random variable sampling.
obtain the sampling of a random variable (such as load, geometrical size, etc.) we need an input to finite element programs. The inverse function method is one widely used. It is suitable for both discrete and continuous distributions. Suppose a random variable r/ has a distribution function F(r/). According to the character of the distribution function it is known that it is monotonic increasing and could have a value between zero and one. Suppose another random variable has uniform distribution between zero and one, then we have: r = F(r/) =
f:
f ( r / ) dr/
(5)
oo
When the value of 7/is obtained between (-o% r/i), the value of r must have a value r(i) between (0, F(r/i)). r(i) is any value of the uniform distribution between (0, 1) (Fig. 1). By this means, we can obtain the sampling value of any given distribution random variable through the transition of the uniform distribution random sequence between (0, 1). Like this, so long as the uniform distribution random data between (0, 1) have been produced, they can be transformed into the sampling of a random variable of any given distribution. When there is no apparent function expression for a relationship between r / a n d r, or even though there is an apparent function expression, but the work load of computing is too large, the sampling transformation should be usedJ
3.4 The calculation of finite element through random sampling For the structure size and external load to which it is subjected, the random distribution law which has been obtained, by means of the method described above, some groups of data sampling value distributed under one of the distribution laws could be obtained, and then every group of data sampling value should be input to the interface of the finite element program, and after that, using the finite element program some stress values at the dangerous points of the structure being analysed will be determined.
Supposing through the finite calculation of n groups data sampling value of structural size and external load to which it is subjected, we found that the n stresses are Sl, s 2 , . . . , s, at the dangerous points in the structure. The stress distribution function Fs(s) can be approximately replaced by the experimental distribution function Fn(s) of the n stress subsamples. Fn(s) is a step function. Let s(1) < s(2) < - • • < s(n) be the statistic of subsample (s~, s2 . . . . , s,) and that the sequence of the data from the large one to the small one was assumed. Then
Fs(n)=Fn(s)=
i
s <~s(i) n s(i)
(6)
s >s(n) and
Fs(s) = lim Fn(s) n---->oo
As has been mentioned above, suppose that the material resistance has a normal distribution, its mean /~r and variance ~ can be taken from Ref. 3; according to this, the value of the distribution function of material resistance r in point s can be determined, it is:
Fr(s) =
O[(s
- Ur)/Or]
(7)
By the interference module dR(t) = fs(s) ds
r) dr
(8)
If we let
G=
fr(r) dr = 1 - Fr(s) = 1 - ~[(s - l~)/a~]
(9)
d H = fs(s) ds
(10)
f H = Jfs(s) ds = Fs(s)
(11)
dR(t) = G d H
(12)
Then
So the structural reliability is:
R(t) =
f2
G dH
(13)
where the limiting values of the integral are considered such that the range of values of G and H is between zero and one. Because the structural reliability is R(t)= S~ G dH, it can be shown that if H stands for horizonatal axis, G, stands for vertical axis, then the area under the curve 1-1 is the structural reliability that can be acquired. It is shown in Fig. 2.
Jin Guoliang, Chen Lin, Dong Jiamei
80
6
should be infinitely approximate. It will be proved by the following calculation, and the times n do not need to be too big, generally less than 20 times, so that the result that is satisfactory for engineering and the relative error compared with the analytic solution is less than 5%.
I
I
4 COMPARISON OF THE MONTE CARLO FINITE ELEMENT METHOD WITH THE ANALYTIC METHOD
t/
Fig. 2. The relationship curve of H-G. The integral method to calculate the area under the curve 1-1 is more precise, the key problem of which is that it is necessary to determine the function of the curve 1-1. Therefore, using the polynomial method to fit the H - G curve, 6 suppose n--1
G = ~, b i x H i
(14)
i=0
where bo, bl, bE . . . . . bn-1 are undetermined coefficients. These undetermined coefficients can be determined through some known arrays of (H, G). According to eqns (6) and (11), we have
i s<~s(i) "In s ( i ) < s < - s ( i + l )
H=
(15)
s >s(n) and G = 1 - O[(s -/~,)/or] It is thus clear that using any stress Si computed by the finite element method, a corresponding group value of (Hi, Gi) can be obtained. Like this, all of the n stresses Si ( i = 1 , 2, 3 . . . . , n) correspond to n groups of values (Hi, Gi). Substituting the n groups of values (Hi, Gi) in eqn (14), n undetermined coefficient equations can be obtained. Resulting from the undetermined coefficient equation group, the undetermined coefficients b0, bl, b2 . . . . . b,_~ can be determined, then eqn (14) can be determined too; substituting eqn (14) in eqn (13), we have
R(t)=
f0'"'
~ bi × H ~d H =
i=0
I?i
b × H/(i + l
In order to compare these methods, a high pressure thick walled cylinder made from ZG35CrMo and with an analytic solution was selected (Fig. 3), and both the analytic method and the Monte Carlo finite element method were used to analyse the structural reliability, respectively. Meanwhile the structural size and material strength of the cylinder were assumed to have a normal distribution. It is clear that from Table 1 and from Figs 4-7, no matter which kind of distribution (normal distribution, log-normal distribution, maximum distribution of extremum I type and Weibull distribution) the external load P has, the result of the structural reliability acquired by means of the Monte Carlo finite element method presented in this paper is consistent with the result obtained by the analytic method such as the design testing point method, the JC method, etc. It is illustrative that the structural reliability of some of the usual distributions in reliability analysis can be obtained more precisely by means of the Monte Carlo finite element method. For the regular thick walled cylinder structure, reliability analysis can be carried out by an existing analysis method. We obtain the theoretical value of the structural reliability under the condition that the internal pressure P to which the structure is subjected agrees with different distributions by using the design
,11
and rl--1
R(t) = ~ [b/(i + 1)]
(16)
i=0
This is the function of the reliability that uses the Monte Carlo finite element method to determine the structural reliability. It can be known from the derivation process described above that R(t) is very closely correlated to the simulation testing times n. The bigger the value of n, the more eqn (14) approximates the H - G curve. When the value of n increases infinitely, both of them
DI DO
!
I
Fig. 3. The thick walled cylinder subject to internal pressure.
Monte Carlo finite element method of structure reliability analysis
81
Table 1. Comparison of the degree of reliabmty acquired by the Monte Carlo finite element method with the analytic method Distribution type of external load P Normal distribution
Log-normal distribution
Maximum distribution of e x t r e m u m I type
Weibull distribution
0.9398
0.8391
0.8554
0.8708
Theoretical value No. of simulations
Computing value
Relative error
Computing value
Relative error
Computing value
Relative error
Computing value
Relative error
5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
0.993 7 0.989 8 0.979 6 0.912 7 0-959 9 0.956 6 0.927 4 0-950 0 0.947 6 0.945 8 0.934 0 0.944 9 0.944 0 0-935 7 0-942 8 0.937 2 0.941 6 0.941 2 0.938 7 0-940 7
5.735 5.320 4-235 -2.884 2.139 1-788 - 1.319 1.085 0.830 0.638 -0.617 0-543 0.447 -0.436 0.319 -0.277 0.192 0-149 -0.117 0-096
0-783 4 0.889 2 0.795 1 0.882 5 0-881 9 0.796 9 0.881 1 0.797 5 0.874 8 0-873 8 0.872 3 0.807 2 0.869 5 0.865 8 0-821 4 0.854 5 0.849 5 0-831 3 0.845 6 0.843 6
-6.638 5-971 -5"244 5.172 5.101 -5-029 5.005 -4.958 4.255 4.135 3.957 -3.802 3.623 3.182 -2-109 1.835 1.239 -0.930 0-775 0.536
0.910 8 0.913 8 0.910 4 0-908 2 0.904 4 0.901 8 0.898 5 0.897 1 0.892 0 0.890 1 0.887 5 0-883 2 0.881 4 0.878 8 0.874 6 0.873 3 0.868 5 0-864 7 0.862 4 0.859 2
6-477 6.827 6-430 6.173 5.728 5.424 5-039 4.875 4.274 4.057 3.753 3.250 3.040 2.736 2-245 2.093 1.531 1.087 0-818 0.444
0.838 3 0.899 9 0.896 7 0.888 1 0.854 8 0-885 1 0-883 3 0-859 3 0-881 1 0-879 9 0.879 5 0.878 4 0.863 7 0.877 6 0.877 3 0-875 6 0.867 1 0.873 7 0.872 4 0-871 0
-3.732 3.342 2-974 1.987 - 1.837 1.642 1.435 - 1.321 1.183 1.045 0.999 0.873 -0.815 0.781 0-746 0-551 -0.425 0.333 0.184 0-023
R~ X~ I. oooo
O.~qS
1tiLl{
~t,,i
~lall
............................
1. oooo
O. foo~
o. qooo o.L~54
O, 900o
~rJ(JX~XXX
~ X X "4.X X ~X 1 ~,~ ~ ~ ..............
O.8000
Jo
,~
Fig. 4. T h e load obeys n o r m a l distribution. , the value of R(t) o b t a i n e d b y using t h e design test p o i n t m e t h o d ; x , t h e value of R(t) o b t a i n e d by using t h e M C F E M ; R(t), reliability; n, no. of simulations.
I
I
I
Fig. 6. T h e load obeys t h e m a x i m u m d i s t r i b u t i o n of e x t r e m u m I type. - - , t h e value of R(t) o b t a i n e d by using the J C m e t h o d ; × , t h e value of R(t) o b t a i n e d by using t h e M C F E M ; n, no. of simulations.
I. 0000
1.0000 0.7000 0.8000
xx 0.9000 0.8708 I
~ v,
o. 8ooo] a
I
I
Fig. 5. T h e load o b e y s l o g - n o r m a l distribution. , the value of R(t) o b t a i n e d b y using t h e J C m e t h o d ; x , t h e value of R(t) o b t a i n e d by using t h e M C F E M ; R(T), reliability; n, no. of simulations.
;o Fig. 7. T h e load obeys t h e W e i b u l l distribution. - - , the value of R(t) o b t a i n e d by using t h e J C m e t h o d ; x , the value of R(t) o b t a i n e d b y using t h e M C F E M ; R(t), reliability; n, no. of simulations.
Jin Guoliang, Chen Lin, Dong Jiamei
82
testing point method and the JC method. By comparing with the result acquired by the Monte Carlo finite element method, it is shown that the convergence speed of the Monte Carlo finite element method is very rapid. It can be shown from Table 1 and from Figs 4-7. When the simulation time n reaches 20 or so, the computing value of the structural reliability acquired by the Monte Carlo finite element method is nearly convergent with the theoretical value, the relative error between them is less than 5%, it is suitable for engineering practice. It is obvious that the Monte Carlo finite element method overcomes the disadvantage that the Monte Carlo method needs to simulate more than ten thousand times, and it takes less time in calculation, it combines the fact that the finite element method can calculate any structural load response with the advantage that the Monte Carlo method can simulate random variables effectively, and it provides a powerful tool for the reliability analysis of complicated structures.
5 TEST U S I N G T H E M O N T E C A R L O FINITE ELEMENT METHOD FOR COMPLEX STRUCTURE RELIABILITY ANALYSIS
steel casing head system was analysed as shown below. The second section structure of the casing head is rather complicated, and it bore the action of both axial load and external load. It is very difficult to obtain the reliability of such complicated structures as shown in Fig. 8, using a common reliability analysis method. The maximum external load which the casing head bears is the pressure when the well is closed. The data illustrated in Table 2, collected from 60 oil wells of some domestic oil fields, are well head pressure when the wells are closed and that which the casing head bears. The external pressure obey the maximum distribution of extremum I type after K - S testing, and its accumulation probability density function is:
Fp(p) = e x p { - e x p [ - 0 . 0 7 8 36(p - 24-0849)]}
Table 3 and Fig. 9 show the result after the Monte Carlo finite element method was used. It is thus clear that the reliability of the second section of the casing head converges at 0.99986 when the number of simulations is less than 20 using the Monte Carlo finite element method.
6 CONCLUSION (1) The Monte Carlo finite element method was developed to calculate structural reliability. It is not only suitable for resolving structural reliability when the random variable has any distribution, but also for resolving the reliability of complicated structures, that is very difficult for the analytic expression of the limit state equation. (2) The result of the analysis of the reliability for the thick walled cylinder possessing an analytic solution by use of the Monte Carlo finite element method indicated that the method
Using the Monte Carlo finite element method, the reliability of the second section of high pressure cast ~[
][.t ~ U l I ~ - ~ ~ 7-'ql~ ~
,
(17)
easin9 hanger vertical hole aut0-c0mpressi0n type
casina connection seal
Fig. 8. The second section of the casing head.
Table 2. The internal pressure that the casing head bears when the wells are closed in some domestic oil fields
No.
Pressure
No.
Pressure
No.
Pressure
No.
Pressure
1 2 3 4 5 6 7 8 9 10 11 12 13 14
23.2 27"2 32.2 33.3 23"8 22"6 35.4 27.0 28.0 23.2 21-6 15.2 17.6 25"6
15
25.1
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
20.2 22"3 22"7 20"6 20.2 19.3 20-8 20.2 20.9 32-5 21-0 20.0 20.5 18-4 26"4
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
19-7 19-2 22"0 15.7 24.5 19.1 32.8 22.3 29-4 24.8 24.0 30.6 37-0 29-0 24-3
46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
28.2 40.7 52.9 50.9 39"3 40-0 50-0 60.0 40.0 45.0 50-0 60.0 67.5 75.0 106.0
Monte Carlo finite element method of structure reliability analysis
83
Table 3
No. of simulations n Reliability R(t) n
R(t) n R(t)
1
2
7 0.999 68 13 0.999 84
8 0.999 63 14 0.999 83
3 0.998 35 9 0-999 92 15 0.999 86
4 0.999 90 10 0.999 78 16 0.999 86
5 0-999 07 11 0.999 80 17 0.999 86
6 0.998 46 12 0-999 80 18 0.999 86
~C~) 1.00000 0.9998B
~
"~
.
.
.
.
.
.
.
.
0.99900
O.ggBoo o.8~7oo
~3
18
a
Fig. 9. The relationship of the reliability R(t) and the no. of simulations n of the second section of the casing head. R(t), reliability; n, no. of simulations. converges rapidly, has rather high precession, and needs less work for computing; it is an effective methodology for determining the reliability of realistic engineering structures. (3) The second section of a high pressure cast steel casing head system in an oil well head and the load that it bears are complex. Its reliability was analysed by use of the Monte Carlo finite element method, which proved again that the Monte Carlo finite element method possessed practical engineering application value.
Wiley, New York, 1980. 4. Lin Huiguo & Zhou Renjun, The World Steel Number Handbook. Mechanical Industry Press, China, 1985. 5. Dai Shuhe & Wang Minge, Reliability Engineering and its Application in Chemical Industry Service. Chemistry Industry Press, China, 1987, p. 194. 6. Ajit Kumar Verma & Murty, A. S. R., A reliability design procedure for arbitrary stress strength distributions. Reliability Engineering and System Safety 26 (1989).
BIBLIOGRAPHY REFERENCES
1. Mu Zhizhong, The Reliability Design of Mechanical Parts. Mechanical Industry Press, China, 1958, p. 116. 2. Jin Guoliang, Chen Lin & Liu Changhai, The study of probability distribution of mechanical properties of the steel used in pressure vessels for the petrochemical industry. Journal of Daqing Petroleum Institute, (1) (1991). 3. Haugen, E. B., Probabilbtic Mechanical Design. John
Huang Kezhong & Mao Shanpei, The Application of Random Methods and Fuzzy Mathematics. Tong Ji University Press, China, 1987, pp. 209-11. Lehmer, D. H., Mathematical methods in large scale computing units. Proc. Second Sym. on Large-scale Digital Calculating Machinery, Harvard University Press, Cambridge, MA, 1951, pp. 141-6. Xu Zhongji, Monte Carlo Method. Shanghai Science and Technology Press, China, 1985, p. 30.
Monte Carlo finite element method of structure reliability analysis Jin Guoliang, Chen Lilt & Dong Jiamei Daqing Petroleum Institute, Anda, Heilongjiang, People's Republic of China (Received 18 August 1991; accepted 10 July 1992)
To resolve the problem in reliability of any complicated engineering structures by use of general reliability methodology one often comes across a great deal of limitations due to realistic difficulties. In order to make use of the advantage that the finite element method could numerically resolve the load response of any structure, and the character that Monte Carlo method could carry out stochastic simulation, the Monte Carlo finite element method for any structure reliability analysis was presented. The advantage of high precision, less work in calculation and rapid convergence of the method developed above has been proved by concrete reliability calculation and analysis of the given engineering structure possessing analytic solutions (e.g. a thick wall cylinder). Using the method for rather complicated structures (e.g. high pressure cast steel casing head bearing of oil well head), it was discovered it is necessary to use only a few simulations to obtain the structure reliability value, which took further steps to prove the effectiveness of the Monte Carlo finite element method.
engineering problems, especially for complicated engineering structures. Therefore, many improved methods were put forward in order to avoid directly obtaining the distribution laws of r and s and the integral of the complex function. Among them, the first order second moment method and the design testing point method are only adapted for resolving a regular structure in which the random variables are suited to the normal distribution and the log-normal distribution, and meanwhile the limit state equation is easily set up. The JC method and the Paloheimo method could be used under the condition that the random variables obey any distribution, but both of them are only adapted for a regular structure where the limit state equation is easily set up with an analytic expression. The Monte Carlo method could be used under the condition that the random variable obeys any distribution, but it is necessary to obtain the expression for the limit state equation and it takes a lot of time to conduct the computing, and sometimes, it is difficult to put into practice because of lack of computer capacity. In a word, all of the methods mentioned above are not perfect for many complicated engineering structures, therefore, it is necessary to do research to think of a new way to analyse the reliability of complicated engineering structures.
1 INTRODUCTION The theory of interference is the basis of structural reliability analysis, it has been used widely for analysing structural reliability. If the structure reliability is assumed as the probability that the stress s in the structure control point is less than the material yield strength, that is, when r
R(t)=f_~fs(s)[~fr(r)dr]ds
(1)
where fr(r) is the probability density function of the material resistance fs(s) is the probability density function of the structural stress, R(t) is the structural reliability. Equation (1) could be used only for resolving the problems in which the distribution laws of both the material resistance r and the stress s of the structure being analysed are known in advance. But it is difficult to obtain the distribution laws of r and s for realistic
ReliabilityEngineeringand SystemSafety0951-8320/93/$06.00 © 1993 Elsevier Science Publishers Ltd, England. 77
78
Jin Guoliang, Chen Lin, Dong Jiamei
2 THE PRESENT SITUATION OF THE MONTE CARLO FINITE ELEMENT METHOD One of the important tasks of structural reliability analysis is to find the structural load response (stress). For complicated engineering structures, very often it is a unique feasible way of adopting a numerical method to calculate their load response. The finite element method has become an effective numerical method along with the rapid development of the technique of the electronic computer. The traditional finite element method is built upon the basis of determinacy of structural system which considers that the structural geometrical size and the load it is subjected to are definite values. But, in fact, the geometrical size of the structure, the load, the material resistance, etc., are usually random variables. Structural reliability analysis is used just to synthesize the indeterminacy of each variable to assess the bearing capacity of the structure. It will have an obvious significance if some appropriate improvement of a large number of effective finite element programs which have been presented at the moment could be completed. The Monte Carlo method is a powerful tool to be used to simulate random variables. It can simulate the random variables that are suitable for any distribution. However, it demands a large amount of calculation if a relative precise solution is needed, which will certainly mean a lot of work for processing the data for a realistic and complicated engineering structure. If some improvement in tackling the problem described above could be successfully made, it would be very useful for the reliability analysis of complicated engineering structures. During the whole process of study, special attention to the basic idea of combination of the Monte Carlo method with the finite element method was paid, and the research work was conducted repeatedly, over and over again, and the Monte Carlo finite element method of reliability analysis for any complicated engineering structures was finally put forward.
3 THE MONTE CARLO FINITE ELEMENT METHOD 3.1 Determination of all basic parameter distribution concerning structure reliability When using the Monte Carlo finite element method to resolve structural reliability, first of all, the basic parameter distribution law for the size, load, material resistance, etc., concerning the structure reliability should be given. According to statistical analysis, the distribution of structural size generally is a normal distribution. If we assume that the + A X and trx
express the tolerance of machine manufacturing and the standard deviation respectively, the relationship between both of them could be established as follows: 1
ax = AS~3
(2)
The distribution law of load could be obtained only from the statistical analysis of a large number of data collected from field structures. The distribution law of material resistance should be obtained from tests in principle. If there are no testing conditions, it can be quoted from technique document directly. Generally, r obeys a normal distribution,z it means,/~r, and standard deviation, or, can be obtained by looking up related statistical materials. 3-5 It is acceptable in engineering applications in this way to estimate the preliminary distribution of the resistance of metal materials.
3.2 Producing uniform distribution random data When using the Monte Carlo finite element method to analyse structural reliability, after the distribution law of each random variable is established, it is necessary to do random sampling for each variable, and meanwhile, to obtain the sampling values of the random variables. In order to obtain the random sampling values of the structural size and external load rapidly and with high precision, two steps may be taken: first to produce uniform distribution random data between zero and one, and then on this basis, transforming them into given variable random sampling. In our study, the multiplier common modulus method was adopted to produce uniform distribution random data between zero and one. The expression for the multiplier common modulus method is as follows:
x(i + 1) = (~. × x(i)) x (mod m)
(3)
r(i + 1) = x(i + 1)/m
(4)
where ~. and m are positive integers. The random data produced in this way are puppet random data, 5'6 if the value of )., x(0) and m were appropriately chosen, the puppet random data sequence can be freely passed through parameter testing, uniform testing and independence testing, then a puppet random data sequence with a long period, good uniformity and good independence could be obtained.
3.3 Sampling random variable of given distribution It will be meaningful that until we have to vary the uniform distribution random data between zero and one into any given distribution random variable and
Monte Carlo finite element method of structure reliability analysis
79
3.5 The establishment of the function of reliability
I
q
0
,
Fig. 1. The audio-visual diagram of random variable sampling.
obtain the sampling of a random variable (such as load, geometrical size, etc.) we need an input to finite element programs. The inverse function method is one widely used. It is suitable for both discrete and continuous distributions. Suppose a random variable r/ has a distribution function F(r/). According to the character of the distribution function it is known that it is monotonic increasing and could have a value between zero and one. Suppose another random variable has uniform distribution between zero and one, then we have: r = F(r/) =
f:
f ( r / ) dr/
(5)
oo
When the value of 7/is obtained between (-o% r/i), the value of r must have a value r(i) between (0, F(r/i)). r(i) is any value of the uniform distribution between (0, 1) (Fig. 1). By this means, we can obtain the sampling value of any given distribution random variable through the transition of the uniform distribution random sequence between (0, 1). Like this, so long as the uniform distribution random data between (0, 1) have been produced, they can be transformed into the sampling of a random variable of any given distribution. When there is no apparent function expression for a relationship between r / a n d r, or even though there is an apparent function expression, but the work load of computing is too large, the sampling transformation should be usedJ
3.4 The calculation of finite element through random sampling For the structure size and external load to which it is subjected, the random distribution law which has been obtained, by means of the method described above, some groups of data sampling value distributed under one of the distribution laws could be obtained, and then every group of data sampling value should be input to the interface of the finite element program, and after that, using the finite element program some stress values at the dangerous points of the structure being analysed will be determined.
Supposing through the finite calculation of n groups data sampling value of structural size and external load to which it is subjected, we found that the n stresses are Sl, s 2 , . . . , s, at the dangerous points in the structure. The stress distribution function Fs(s) can be approximately replaced by the experimental distribution function Fn(s) of the n stress subsamples. Fn(s) is a step function. Let s(1) < s(2) < - • • < s(n) be the statistic of subsample (s~, s2 . . . . , s,) and that the sequence of the data from the large one to the small one was assumed. Then
Fs(n)=Fn(s)=
i
s <~s(i) n s(i)
(6)
s >s(n) and
Fs(s) = lim Fn(s) n---->oo
As has been mentioned above, suppose that the material resistance has a normal distribution, its mean /~r and variance ~ can be taken from Ref. 3; according to this, the value of the distribution function of material resistance r in point s can be determined, it is:
Fr(s) =
O[(s
- Ur)/Or]
(7)
By the interference module dR(t) = fs(s) ds
r) dr
(8)
If we let
G=
fr(r) dr = 1 - Fr(s) = 1 - ~[(s - l~)/a~]
(9)
d H = fs(s) ds
(10)
f H = Jfs(s) ds = Fs(s)
(11)
dR(t) = G d H
(12)
Then
So the structural reliability is:
R(t) =
f2
G dH
(13)
where the limiting values of the integral are considered such that the range of values of G and H is between zero and one. Because the structural reliability is R(t)= S~ G dH, it can be shown that if H stands for horizonatal axis, G, stands for vertical axis, then the area under the curve 1-1 is the structural reliability that can be acquired. It is shown in Fig. 2.
Jin Guoliang, Chen Lin, Dong Jiamei
80
6
should be infinitely approximate. It will be proved by the following calculation, and the times n do not need to be too big, generally less than 20 times, so that the result that is satisfactory for engineering and the relative error compared with the analytic solution is less than 5%.
I
I
4 COMPARISON OF THE MONTE CARLO FINITE ELEMENT METHOD WITH THE ANALYTIC METHOD
t/
Fig. 2. The relationship curve of H-G. The integral method to calculate the area under the curve 1-1 is more precise, the key problem of which is that it is necessary to determine the function of the curve 1-1. Therefore, using the polynomial method to fit the H - G curve, 6 suppose n--1
G = ~, b i x H i
(14)
i=0
where bo, bl, bE . . . . . bn-1 are undetermined coefficients. These undetermined coefficients can be determined through some known arrays of (H, G). According to eqns (6) and (11), we have
i s<~s(i) "In s ( i ) < s < - s ( i + l )
H=
(15)
s >s(n) and G = 1 - O[(s -/~,)/or] It is thus clear that using any stress Si computed by the finite element method, a corresponding group value of (Hi, Gi) can be obtained. Like this, all of the n stresses Si ( i = 1 , 2, 3 . . . . , n) correspond to n groups of values (Hi, Gi). Substituting the n groups of values (Hi, Gi) in eqn (14), n undetermined coefficient equations can be obtained. Resulting from the undetermined coefficient equation group, the undetermined coefficients b0, bl, b2 . . . . . b,_~ can be determined, then eqn (14) can be determined too; substituting eqn (14) in eqn (13), we have
R(t)=
f0'"'
~ bi × H ~d H =
i=0
I?i
b × H/(i + l
In order to compare these methods, a high pressure thick walled cylinder made from ZG35CrMo and with an analytic solution was selected (Fig. 3), and both the analytic method and the Monte Carlo finite element method were used to analyse the structural reliability, respectively. Meanwhile the structural size and material strength of the cylinder were assumed to have a normal distribution. It is clear that from Table 1 and from Figs 4-7, no matter which kind of distribution (normal distribution, log-normal distribution, maximum distribution of extremum I type and Weibull distribution) the external load P has, the result of the structural reliability acquired by means of the Monte Carlo finite element method presented in this paper is consistent with the result obtained by the analytic method such as the design testing point method, the JC method, etc. It is illustrative that the structural reliability of some of the usual distributions in reliability analysis can be obtained more precisely by means of the Monte Carlo finite element method. For the regular thick walled cylinder structure, reliability analysis can be carried out by an existing analysis method. We obtain the theoretical value of the structural reliability under the condition that the internal pressure P to which the structure is subjected agrees with different distributions by using the design
,11
and rl--1
R(t) = ~ [b/(i + 1)]
(16)
i=0
This is the function of the reliability that uses the Monte Carlo finite element method to determine the structural reliability. It can be known from the derivation process described above that R(t) is very closely correlated to the simulation testing times n. The bigger the value of n, the more eqn (14) approximates the H - G curve. When the value of n increases infinitely, both of them
DI DO
!
I
Fig. 3. The thick walled cylinder subject to internal pressure.
Monte Carlo finite element method of structure reliability analysis
81
Table 1. Comparison of the degree of reliabmty acquired by the Monte Carlo finite element method with the analytic method Distribution type of external load P Normal distribution
Log-normal distribution
Maximum distribution of e x t r e m u m I type
Weibull distribution
0.9398
0.8391
0.8554
0.8708
Theoretical value No. of simulations
Computing value
Relative error
Computing value
Relative error
Computing value
Relative error
Computing value
Relative error
5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
0.993 7 0.989 8 0.979 6 0.912 7 0-959 9 0.956 6 0.927 4 0-950 0 0.947 6 0.945 8 0.934 0 0.944 9 0.944 0 0-935 7 0-942 8 0.937 2 0.941 6 0.941 2 0.938 7 0-940 7
5.735 5.320 4-235 -2.884 2.139 1-788 - 1.319 1.085 0.830 0.638 -0.617 0-543 0.447 -0.436 0.319 -0.277 0.192 0-149 -0.117 0-096
0-783 4 0.889 2 0.795 1 0.882 5 0-881 9 0.796 9 0.881 1 0.797 5 0.874 8 0-873 8 0.872 3 0.807 2 0.869 5 0.865 8 0-821 4 0.854 5 0.849 5 0-831 3 0.845 6 0.843 6
-6.638 5-971 -5"244 5.172 5.101 -5-029 5.005 -4.958 4.255 4.135 3.957 -3.802 3.623 3.182 -2-109 1.835 1.239 -0.930 0-775 0.536
0.910 8 0.913 8 0.910 4 0-908 2 0.904 4 0.901 8 0.898 5 0.897 1 0.892 0 0.890 1 0.887 5 0-883 2 0.881 4 0.878 8 0.874 6 0.873 3 0.868 5 0-864 7 0.862 4 0.859 2
6-477 6.827 6-430 6.173 5.728 5.424 5-039 4.875 4.274 4.057 3.753 3.250 3.040 2.736 2-245 2.093 1.531 1.087 0-818 0.444
0.838 3 0.899 9 0.896 7 0.888 1 0.854 8 0-885 1 0-883 3 0-859 3 0-881 1 0-879 9 0.879 5 0.878 4 0.863 7 0.877 6 0.877 3 0-875 6 0.867 1 0.873 7 0.872 4 0-871 0
-3.732 3.342 2-974 1.987 - 1.837 1.642 1.435 - 1.321 1.183 1.045 0.999 0.873 -0.815 0.781 0-746 0-551 -0.425 0.333 0.184 0-023
R~ X~ I. oooo
O.~qS
1tiLl{
~t,,i
~lall
............................
1. oooo
O. foo~
o. qooo o.L~54
O, 900o
~rJ(JX~XXX
~ X X "4.X X ~X 1 ~,~ ~ ~ ..............
O.8000
Jo
,~
Fig. 4. T h e load obeys n o r m a l distribution. , the value of R(t) o b t a i n e d b y using t h e design test p o i n t m e t h o d ; x , t h e value of R(t) o b t a i n e d by using t h e M C F E M ; R(t), reliability; n, no. of simulations.
I
I
I
Fig. 6. T h e load obeys t h e m a x i m u m d i s t r i b u t i o n of e x t r e m u m I type. - - , t h e value of R(t) o b t a i n e d by using the J C m e t h o d ; × , t h e value of R(t) o b t a i n e d by using t h e M C F E M ; n, no. of simulations.
I. 0000
1.0000 0.7000 0.8000
xx 0.9000 0.8708 I
~ v,
o. 8ooo] a
I
I
Fig. 5. T h e load o b e y s l o g - n o r m a l distribution. , the value of R(t) o b t a i n e d b y using t h e J C m e t h o d ; x , t h e value of R(t) o b t a i n e d by using t h e M C F E M ; R(T), reliability; n, no. of simulations.
;o Fig. 7. T h e load obeys t h e W e i b u l l distribution. - - , the value of R(t) o b t a i n e d by using t h e J C m e t h o d ; x , the value of R(t) o b t a i n e d b y using t h e M C F E M ; R(t), reliability; n, no. of simulations.
Jin Guoliang, Chen Lin, Dong Jiamei
82
testing point method and the JC method. By comparing with the result acquired by the Monte Carlo finite element method, it is shown that the convergence speed of the Monte Carlo finite element method is very rapid. It can be shown from Table 1 and from Figs 4-7. When the simulation time n reaches 20 or so, the computing value of the structural reliability acquired by the Monte Carlo finite element method is nearly convergent with the theoretical value, the relative error between them is less than 5%, it is suitable for engineering practice. It is obvious that the Monte Carlo finite element method overcomes the disadvantage that the Monte Carlo method needs to simulate more than ten thousand times, and it takes less time in calculation, it combines the fact that the finite element method can calculate any structural load response with the advantage that the Monte Carlo method can simulate random variables effectively, and it provides a powerful tool for the reliability analysis of complicated structures.
5 TEST U S I N G T H E M O N T E C A R L O FINITE ELEMENT METHOD FOR COMPLEX STRUCTURE RELIABILITY ANALYSIS
steel casing head system was analysed as shown below. The second section structure of the casing head is rather complicated, and it bore the action of both axial load and external load. It is very difficult to obtain the reliability of such complicated structures as shown in Fig. 8, using a common reliability analysis method. The maximum external load which the casing head bears is the pressure when the well is closed. The data illustrated in Table 2, collected from 60 oil wells of some domestic oil fields, are well head pressure when the wells are closed and that which the casing head bears. The external pressure obey the maximum distribution of extremum I type after K - S testing, and its accumulation probability density function is:
Fp(p) = e x p { - e x p [ - 0 . 0 7 8 36(p - 24-0849)]}
Table 3 and Fig. 9 show the result after the Monte Carlo finite element method was used. It is thus clear that the reliability of the second section of the casing head converges at 0.99986 when the number of simulations is less than 20 using the Monte Carlo finite element method.
6 CONCLUSION (1) The Monte Carlo finite element method was developed to calculate structural reliability. It is not only suitable for resolving structural reliability when the random variable has any distribution, but also for resolving the reliability of complicated structures, that is very difficult for the analytic expression of the limit state equation. (2) The result of the analysis of the reliability for the thick walled cylinder possessing an analytic solution by use of the Monte Carlo finite element method indicated that the method
Using the Monte Carlo finite element method, the reliability of the second section of high pressure cast ~[
][.t ~ U l I ~ - ~ ~ 7-'ql~ ~
,
(17)
easin9 hanger vertical hole aut0-c0mpressi0n type
casina connection seal
Fig. 8. The second section of the casing head.
Table 2. The internal pressure that the casing head bears when the wells are closed in some domestic oil fields
No.
Pressure
No.
Pressure
No.
Pressure
No.
Pressure
1 2 3 4 5 6 7 8 9 10 11 12 13 14
23.2 27"2 32.2 33.3 23"8 22"6 35.4 27.0 28.0 23.2 21-6 15.2 17.6 25"6
15
25.1
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
20.2 22"3 22"7 20"6 20.2 19.3 20-8 20.2 20.9 32-5 21-0 20.0 20.5 18-4 26"4
31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
19-7 19-2 22"0 15.7 24.5 19.1 32.8 22.3 29-4 24.8 24.0 30.6 37-0 29-0 24-3
46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
28.2 40.7 52.9 50.9 39"3 40-0 50-0 60.0 40.0 45.0 50-0 60.0 67.5 75.0 106.0
Monte Carlo finite element method of structure reliability analysis
83
Table 3
No. of simulations n Reliability R(t) n
R(t) n R(t)
1
2
7 0.999 68 13 0.999 84
8 0.999 63 14 0.999 83
3 0.998 35 9 0-999 92 15 0.999 86
4 0.999 90 10 0.999 78 16 0.999 86
5 0-999 07 11 0.999 80 17 0.999 86
6 0.998 46 12 0-999 80 18 0.999 86
~C~) 1.00000 0.9998B
~
"~
.
.
.
.
.
.
.
.
0.99900
O.ggBoo o.8~7oo
~3
18
a
Fig. 9. The relationship of the reliability R(t) and the no. of simulations n of the second section of the casing head. R(t), reliability; n, no. of simulations. converges rapidly, has rather high precession, and needs less work for computing; it is an effective methodology for determining the reliability of realistic engineering structures. (3) The second section of a high pressure cast steel casing head system in an oil well head and the load that it bears are complex. Its reliability was analysed by use of the Monte Carlo finite element method, which proved again that the Monte Carlo finite element method possessed practical engineering application value.
Wiley, New York, 1980. 4. Lin Huiguo & Zhou Renjun, The World Steel Number Handbook. Mechanical Industry Press, China, 1985. 5. Dai Shuhe & Wang Minge, Reliability Engineering and its Application in Chemical Industry Service. Chemistry Industry Press, China, 1987, p. 194. 6. Ajit Kumar Verma & Murty, A. S. R., A reliability design procedure for arbitrary stress strength distributions. Reliability Engineering and System Safety 26 (1989).
BIBLIOGRAPHY REFERENCES
1. Mu Zhizhong, The Reliability Design of Mechanical Parts. Mechanical Industry Press, China, 1958, p. 116. 2. Jin Guoliang, Chen Lin & Liu Changhai, The study of probability distribution of mechanical properties of the steel used in pressure vessels for the petrochemical industry. Journal of Daqing Petroleum Institute, (1) (1991). 3. Haugen, E. B., Probabilbtic Mechanical Design. John
Huang Kezhong & Mao Shanpei, The Application of Random Methods and Fuzzy Mathematics. Tong Ji University Press, China, 1987, pp. 209-11. Lehmer, D. H., Mathematical methods in large scale computing units. Proc. Second Sym. on Large-scale Digital Calculating Machinery, Harvard University Press, Cambridge, MA, 1951, pp. 141-6. Xu Zhongji, Monte Carlo Method. Shanghai Science and Technology Press, China, 1985, p. 30.