第05章-结构可靠度计算方法

结构可靠度计算方法

第五章目 录5.1

一次二阶矩法5.2梯度优化法

5.3非正态分布变量和统计相关变量条件下可靠指标的计算

5.4蒙特卡罗（MonteCarlo）法

5.5响应面法

5.6结构的可靠度设计

5.1一次二阶矩法

一次二阶矩法是针对结构功能函数为变量的一次（即线性）函数，以变量的一阶矩和二阶矩为概率特征进行可靠度计算的一种方法。对于非线性功能函数，一般在某点进行泰勒级数展开并近似地取其一次式，使结构功能函数线性化，然后再用一次二阶矩法计算可靠指标。设影响结构可靠度的n个随机变量为Xi(i=1,2,…,n),结构的功能函数为：Z=g(X1,

X2,

…,Xn)(5-1)极限状态方程为：Z=g(X1,

X2,

…,Xn)=0(5-2)一次二阶矩法将结构功能函数在点展开为泰勒级数一次二阶矩法为了得到线性极限状态方程，近似地只取到一次项：式中表示该导数在点X0i(i=1,2,…,n)处取值。(5-3)一次二阶矩法1.均值一次二阶矩法在可靠性设计的初期，各个随机变量的分布规律难以确定，这些变量的一阶矩（均值）和二阶矩（方差）则较容易得到，对于非线性的结构功能函数，则在均值点进行泰勒级数展开并取其一次式，使结构可靠性设计时计算简单，并称之为均值一次二阶矩法。由于结构功能函数是在均值点展开的，故又简称为中心点法。一次二阶矩法在上述计算过程中并没有考虑到变量的实际分布情况，而只考虑了它们的均值和方差，或者说，是将各个随机变量假定为正态分布或者对数正态分布变量进行计算的。理论与实践均证明，对于非线性极限状态方程，均值一次二阶矩法的计算误差较大，以至于选择形式不同但力学意义等效的非线性极限状态方程时，所得的可靠指标值大不相同。

一次二阶矩法例5-1

圆截面直杆如图5-1所示，承受拉力p=100kN，设杆的直径d和材料的应力屈服极限Fy为随机变量，其均值和标准差分别为：md=3cmσd=0.3cm

MFy=29000.0N/cm2

σFy=2500.0N/cm2试计算直杆的抗拉可靠指标。

图5-1一次二阶矩法

解：若建立用内力表示的极限状态方程，可得：图5-1一次二阶矩法得到线性化的极限状态方程为：

一次二阶矩法由式（4-24）和(4-25)得：一次二阶矩法最后得到可靠指标为：

为了说明采用不同形式的极限状态方程将得到不同的值，下面改用应力表示的极限状态方程进行计算。用应力表示的极限状态方程为：将此方程在均值点线性化得：一次二阶矩法简化后得到：变量的均值和标准差为：最后得到可靠指标：一次二阶矩法由上述结果可见，在同一个问题中获得两个不同的值，这是均值一次二阶矩法存在的严重问题。为了克服这个缺点，人们对均值一次二阶矩法进行改进。一次二阶矩法2.改进一次二阶矩法

均值一次二阶矩法是在均值点附近将非线性功能函数线性化，并据此计算可靠指标，由于均值点一般在可靠区，且距失效边界较远，显然求得的可靠指标误差很大（参见图5-2）。改进一次二阶矩法常简称为一次二阶矩法，又称为验算点法。这一方法是将非线性功能函数的线性化点选为设计验算点P*，并据此计算可靠指标，使得到的可靠指标值具有较高的精度。一次二阶矩法一次二阶矩法对于任意一组随机变量，相应的设计验算点为，根据式(5-3)得在验算点处的线性化极限状态方程为：Z的均值为：

(5-5)

(5-4)

一次二阶矩法设计验算点就在失效边界上，所以有：(5-6)

式(5-5)成为：

(5-7)

在变量相互独立的情况下，由式(4-26)求出的标准差：

(5-8)

一次二阶矩法目前一般采用迭代法得到可靠指标和设计验算点的值。采用式(4-45)将上述根式线性化，得：(5-9)

式中αi可由式(4-44)求得，即(5-10)

一次二阶矩法

αi可以反映变量Xi对综合变量Z的标准差的影响，因此被称为灵敏系数。αi的值在±1之间，且满足下列关系式：(5-11)

一次二阶矩法根据可靠指标的定义有：(5-12)

整理上式得到：(5-13)

一次二阶矩法由式(5-13)可得：

(5-14)

于是得到设计验算点的计算公式：(5-15)

一次二阶矩法在已知、的条件下，要求解可靠指标和设计验算点的值，需要(n+1)个方程。式(5-15)表示n个方程，由于设计验算点位于极限状态面上，所以一般再补充方程（5-6）,并采用迭代法求解。一次二阶矩法适于电算编程的计算方法,其迭代计算步骤如下：

(1)假定一个β值。(2)选取设计验算点的初值，一般取。

(3)计算的值。(4)由式(5-10)计算αi值。一次二阶矩法(5)由式(5-15)计算新的验算点值。

(6)重复步骤(3)至步骤(5),直到前后两次算得的的差值在容许范围内为止。(7)将所得的值代入原极限状态方程(5-2)式，计算结构功能函数g的值。

一次二阶矩法(8)检验极限状态方程的条件是否满足，如果不满足，则计算前后两次β和g的各自差值的比值，并由式估算一个新的β值，然后重复步骤(3)到(7)的计算，直到满足为止。

(9)最后由计算失效概率。一次二阶矩法

设某梁的抗力（应力）为f，荷载作用下梁的最大弯矩为M，弯矩转换为梁应力的参数为A。已知f、M、A为随机变量，其统计特性如下：试求该梁的失效概率及设计验算点的值。

例5-2一次二阶矩法

解：梁的极限状态方程为：按照一次二阶矩法的计算公式，可得：一次二阶矩法第一次迭代：

设初值β=2.5；设计验算点的初值取其均值，即。

(2)计算αi的值得：

一次二阶矩法(3)计算新的设计验算点的值为：

第二次迭代得到：

一次二阶矩法第三次迭代得到的计算结果为：

一次二阶矩法

第三次迭代得到的结果与第二次的结果已很接近，所以第一大轮的迭代可以结束。如果所假定的正确，那么算得的设计验算点的值应满足极限状态方程代入该方程进行检验：一次二阶矩法由检验结果可见，算得的设计验算点不满足极限状态方程，必须进行第二大轮的迭代计算。

再次迭代计算假定β=2.0

，计算结果列于表5-1。

一次二阶矩法

迭代次序1230.8850.88060.8802-0.427-0.4333-0.4341-0.171-0.1916-0.1920（初值3927）

451645344532（初值218140）

211570211270211920

（初值0.0247）

0.024490.024550.02455变量

表5-1

第二大轮的迭代计算结果

一次二阶矩法与β=2.0对应的值为：

为了加快迭代收敛的速度，第三大轮迭代时β的假定值由下式估算：

一次二阶矩法按照β=1.6的假设，算得的设计验算点为：

对应的

第四大轮迭代时，由βn+1式再次估算出β=1.58

，算得的设计验算点为：

一次二阶矩法对应的已近似为零，整个迭代过程结束。

最后得到失效概率为：

一次二阶矩法对应的设计验算点值为：

一次二阶矩法用改进一次二阶矩法计算例5-1圆截面直杆的可靠指标。

本例采用内力表示的极限状态方程计算：即

例5-3一次二阶矩法式中Fy和d为随机变量。用改进一次二阶矩法计算，将最后一大轮（β=

2.872）的迭代计算结果列于表5-2。

最后得到的设计验算点值为：

一次二阶矩法将他们代入极限状态方程得：

这就是最终的设计验算点值，β=

2.872即为该杆的可靠指标。

一次二阶矩法表5-2

最后一大轮(β=

2.872)的迭代计算结果

迭代次数

变量

X1(d)(3.0)X2(Fy)(29000)1X1X21366607.0690.918330.395822.20926157.92X1X290755.63.8620.943320.331892.18726617.0一次二阶矩法3X1X291448.23.7570.946080.323932.18526674.24X1X291544.83.7490.946400.322992.18526680.95X1X291556.43.7480.946440.322882.18526681.76X1X291557.63.7840.946440.322872.81526681.8一次二阶矩法本例的极限状态方程是非线性方程，其失效边界线见图5-3。由图示结果可见，改进一次二阶矩法算得的可靠指标是选失效边界上的设计验算点作为线性化点得到的，算得的可靠指标是原极限状态方程上设计验算点到均值点的距离，也是失效边界线到均值点的最短距离，因此设计验算点是结构的最大可能失效点，由此可看出，改进一次二阶矩法的计算精度是较高的。一次二阶矩法均值一次二阶矩法的线性化点取为均值点，而均值点一般距失效边界较远，所以将原非线性极限状态方程在该点通过泰勒级数展开得到线性化方程，并据此计算可靠指标，其计算误差肯定是很大的。正是由于上述原因，同一结构可靠度计算问题列出不同形式的极限状态方程，得到的线性化极限状态方程不一样，算得的可靠指标也就不同，图5-3中给出的计算结果和图形表示就说明了上述分析的正确性。一次二阶矩法图5-3返回首页一次二阶矩法5.2

梯度优化法

计算可靠指标的几何涵义是求均值点到失效边界的最短距离，而设计验算点就是失效边界上距均值点最近的点。若基本随机变量是不相关的，则可首先将其转换为标准正态分布变量，然后在标准正态变量空间内计算可靠指标，这一计算过程可归结为约束优化问题进行求解，使可靠指标的计算更为方便，对于较复杂的可靠指标计算问题，其优点更加明显。这一计算方法称为梯度优化法（也称为几何法）。

根据几何学的知识，在标准正态变量空间内，极限状态面G(Y)=0上任一点Y=(Y1,Y2,…Yn)到坐标原点的距离为：

(5-16)

梯度优化法MinimizedSubjecttoG(Y)=0(5-17)设极限状态面上到原点距离最小的点为则求解坐标原点到极限状态面的最小距离就归结为以下优化问题：梯度优化法设已知结构可靠度分析中的一组基本随机变量X=(X1,X2,…Xn),

描述结构的功能函数为g(X)=g(X1,X2,…Xn)相应的极限状态方程为g(X)=0。用表示抗力的变量R=R(X)和表示荷载效应的变量S=S(X)做成数学式子，得到g(X)=R(X)-S(X)g(R,S)<0表示超过极限状态的情况，定义为失效事件；g(R,S)

>0则表示结构处于安全状态。梯度优化法为简便起见，设基本变量X=(X1,X2,…Xn)为相互独立的正态变量，通过变换

Y=T(X)(5-18)可以得到一组相互独立的标准正态变量Y=(Y1,Y2,…Yn)，于是功能函数也可转换到标准正态空间(5-19)梯度优化法采用迭代的方法可以确定极限状态面G(Y)=0上距原点最近的点

，然后按照求距离的公式(5-16)计算结构的可靠指标：

(5-20)并由式

(5-21)计算结构的失效概率。

梯度优化法采用的迭代格式如下：

(5-22)是梯度向量(5-23)式中：梯度优化法

α是沿着负梯度方向的单位向量，它垂直于极限状态面，指向背离坐标原点。一般经过3-5个循环的迭代计算，序列Yi便收敛于极限状态面上距原点最近的点Y*，并能达到所需的精度。(5-24)梯度优化法

在例5-1中，圆截面直杆承受的拉力P=100KN，设杆的直径d和材料的应力屈服极限Fy均为正态分布变量，其均值和标准差分别为：

例5-4试用梯度优化法计算直杆的抗拉可靠指标。

图5-1梯度优化法解：用内力表示的极限状态方程为：

用应力表示的极限状态方程为：

根据前面导出的梯度优化法计算公式，用内力表示的极限状态方程算得的结果列于表5-3，用应力表示的极限状态方程算得的结果列于表5-4。梯度优化法表

5-3

用内力表示的极限状态方程的迭代结果

迭代次数

X1(Fy)X2(d)Y1Y2βg(X)迭代初值

29000.03.0000.0000.000----10498.9412665.082.352-0.929-2.1592.3501579.12622654.332.200-0.983-2.6742.85068.89532665.802.190-0.937-2.7152.872-0.08742665.102.185-0.928-2.7182.8720.84652668.102.185-0.928-2.7182.872-0.087梯度优化法表

5-4

用应力表示的极限状态方程的迭代结果迭代次数

X1(Fy)X2(d)Y1Y2βg(X)迭代初值

29000.03.0000.0000.000----1485.29012249.002.116-2.604-2.9473.932-594.66022686.872.175-0.853-2.7512.880-4.60132671.302.183-0.917-2.7252.875-1.46542668.752.184-0.925-2.7202.873-0.29552668.052.184-0.928-2.7202.873-0.263梯度优化法与例5-3一次二阶矩法算得的结果相比，二者的计算精度基本相同，但计算工作量却大大减少，这是梯度优化法的优点之一。梯度优化法本质上与一次二阶矩法相同，都是考虑变量的一阶、二阶矩，采用线性近似的方法计算可靠指标，只是迭代计算的格式不同而已。返回首页梯度优化法5.3

非正态分布变量和统计相关变量条件下可靠指标的计算1.将非正态变量转换为等效正态分布变量设Xi为非正态分布变量，可按以下原则将其转换为等效正态分布变量：

(1)在设计验算点

处，等效正态变量（其均值为，标准差为）的概率分布函数值与原变量Xi（其均值为，标准差为）的概率分布函数值相等。

变量条件下可靠指标的计算非正态分布变量和统计相关(2)在设计验算点处，等效正态变量的概率密度函数值与原变量Xi的概率密度函数值相等，见图5-4。

变量条件下可靠指标的计算非正态分布变量和统计相关变量条件下可靠指标的计算非正态分布变量和统计相关根据条件(1)有：

由于

于是可得等效正态分布变量的均值为：

(5-25)变量条件下可靠指标的计算非正态分布变量和统计相关根据条件(2)有：

(5-26)

考虑到：

变量条件下可靠指标的计算非正态分布变量和统计相关并将式(5-25)代入式(5-26)，可得等效正态分布变量的标准差为：

(5-27)变量条件下可靠指标的计算非正态分布变量和统计相关分别为原变量Xi的概率分布函数和概率密度函数，分别为标准正态分布变量的概率分布函数和概率密度函数，其值可查附表1或通过式（2-31）计算获得。

变量条件下可靠指标的计算非正态分布变量和统计相关(5-28)(5-29)(1).由于随机变量Xi服从对数正态分布，所以变量服从正态分布。于是可将变量Xi转换为变量Yi后再继续计算，此时Yi的均值和标准差分别为：

当变量为对数正态分布时，可得到如下两种简便的计算公式：变量条件下可靠指标的计算非正态分布变量和统计相关(2).等效正态分布变量均值和标准差的计算公式如下：

(5-30)(5-31)变量条件下可靠指标的计算非正态分布变量和统计相关将非正态变量转换为等效正态变量以后，即可用前面所述的方法计算结构的可靠指标。由于这一方法已经被国际安全度联合委员会(JCSS)所采用，故称为JC法。我国各工程结构可靠度设计统一标准中都规定采用这一方法进行结构的可靠度计算。变量条件下可靠指标的计算非正态分布变量和统计相关2.将统计相关变量转换为相互独立变量对于一组正态变量，设其均值为，标准差为，相关系数矩阵为：

变量条件下可靠指标的计算非正态分布变量和统计相关(5-32)

其中为变量Xi与Xj的相关系数。

设与X相应的标准正态随机变量为，于是有

(5-32)变量条件下可靠指标的计算非正态分布变量和统计相关式中

由于变量Yi与Yj的协方差为：

变量条件下可靠指标的计算非正态分布变量和统计相关所以标准正态变量Y的协方差矩阵为：

(5-34)

通过正交变换可将矩阵CY化为对角矩阵CZ’，其对角线元素即为CY的特征值。

变量条件下可靠指标的计算非正态分布变量和统计相关令

(5-35)式中A为正交矩阵，其列向量为CY的特征向量，向量即为一组不相关的正态变量，它的协方差矩阵为：

(5-36)变量条件下可靠指标的计算非正态分布变量和统计相关进而再得到一组不相关的标准正态变量，

并且有：

(5-37)(5-38)变量条件下可靠指标的计算非正态分布变量和统计相关其中：

(5-40)而

(5-39)变量条件下可靠指标的计算非正态分布变量和统计相关3.

变量为相关正态时结构可靠指标的计算设已知n维正态变量X的均值mX、标准差和相关矩阵，则计算步骤如下：

(1)求CY的特征值和特征向量，得到A和。

(2)建立矩阵T和B。

(3)确定迭代初值

变量条件下可靠指标的计算非正态分布变量和统计相关(4)计算

其中为结构功能函数，

(5)求单位向量：

变量条件下可靠指标的计算非正态分布变量和统计相关(6)计算新的验算点值：

(7)求可靠指标：

(8)判别g(X)是否满足精度要求，否则重复步骤(4)~(8)。

(9)迭代收敛后，可由值计算失效概率：

变量条件下可靠指标的计算非正态分布变量和统计相关例5-5

钢梁截面的抗力矩为YZ，Y为钢的屈服极限，Z为与截面尺寸有关的参数，已知在荷载作用下某一截面的实际弯矩为M，则钢梁抗弯功能函数为：

相应的极限状态方程为：

变量条件下可靠指标的计算非正态分布变量和统计相关设变量Y与Z的相关系数为，M与Y、Z在统计上独立，Y、Z、M都是正态分布，其统计特性列于表5-5，求梁截面的可靠指标及破坏概率。

表5-5

随机变量的统计特性

均

值

标

准

差

变

异

系

数

变量条件下可靠指标的计算非正态分布变量和统计相关解：引用记号X=(Y,Z,M),相关系数矩阵为

求得

变量条件下可靠指标的计算非正态分布变量和统计相关迭代计算结果列于表5-6。

变量条件下可靠指标的计算非正态分布变量和统计相关迭代次数

Xixi*1YZM40501000292.8782.17-200-0.804-0.2260.549-2.209-0.6211.5082.7462YZM29.05846.22981301.6260.19688.182-200-0.766-0.2590.549-2.2073-0.74631.69732.883表5-6

迭代计算结果

变量条件下可靠指标的计算非正态分布变量和统计相关3YZM28.72246.40481339.46254.86587.7680-200-0.7585-0.26190.5967-2.1738-0.74891.70822.8634YZM28.86346.48451341.64254.86587.793-200-0.7593-0.26160.5959-2.1738-0.74891.70602.863注：上标*表示在验算点取值。

变量条件下可靠指标的计算非正态分布变量和统计相关迭代结束得：

＝2.863

最大可能破坏点（即验算点）值为：

Y*=28.855,Z*=46.4785,M*=1341.2若不考虑Y与Z的相关性，解得＝3.049，Pf＝0.0014。

变量条件下可靠指标的计算非正态分布变量和统计相关4.变量为相关非正态时结构可靠指标的计算

设已知n维正态变量X的均值mX、标准差和相关矩阵。计算步骤如下：

(1)求CY的特征值和特征向量，得到A和。

(2)

确定迭代初值建立矩阵。

变量条件下可靠指标的计算非正态分布变量和统计相关(3)求非正态变量的正态当量均值和当量标准差，得T和B。

(4)同相关正态变量的步骤(4)~(7)

。

(5)判别G(X)是否满足精度要求，否则重复步骤(3)~(5)。

(6)迭代收敛后，可由值计算破坏概率变量条件下可靠指标的计算非正态分布变量和统计相关例5-6

将例5-5中的变量分布改为：Y－对数正态，Z－对数正态，M－极值I型。计算梁截面的可靠指标和破坏概率。

变量条件下可靠指标的计算非正态分布变量和统计相关解：迭代计算结果列于表5-7。

迭代结束得

＝2.664

最大可能破坏点（即验算点）为Y*=33.753,Z*=47.766,M*=1612.998本例若不考虑Y与Z的相关性，解得＝2.746，＝0.00301。

变量条件下可靠指标的计算非正态分布变量和统计相关表5-7

迭代计算结果

迭代次数Xixi*1YZM4050100052.5191.06139.68549.949965.609292.8782.17－191.0612YZM28.097545.95651257.41073.51222.2978287.072937.800149.7849892.828189.069853.0426－287.0729变量条件下可靠指标的计算非正态分布变量和统计相关3YZM32.577147.42821547.38

4.07212.3714510.79

39.007549.8842525.80

226.231763.4670－510.79

4YZM34.945148.25281687.45

4.36812.4126434.023

39.390949.9197615.413

246.894769.7646－434.023

5YZM33.732747.7701612.27

4.21662.3885413.409

39.215349.9006665.56

235.94866.1928－413.409

变量条件下可靠指标的计算非正态分布变量和统计相关Xi1YZM－0.815－0.2290.532－2.3395－0.65741.52722.87022YZM－0.5436－0.15250.8254－1.5016－0.42132.28012.7625变量条件下可靠指标的计算非正态分布变量和统计相关3YZM－0.4024－0.11290.9085－1.0073－0.28262.27422.50334YZM－0.4898－0.13840.8610－1.3066－0.36922.29682.6685YZM－0.4910－0.13770.8602－1.3081－0.36682.29172.664变量条件下可靠指标的计算非正态分布变量和统计相关

这里给出的变换方法只适用于相关正态变量；对于相关非正态变量，需要采用Rosenblatt变换，而这里提出的变换方法可用作近似计算。例5-6采用Rosenblatt变换得到的结果为

＝2.663返回首页变量条件下可靠指标的计算非正态分布变量和统计相关5.4蒙特卡罗（MonteCarlo）法

蒙特卡罗法主要用于求解具有随机性的不确定性问题，但也能求解确定性问题。蒙特卡罗法又称随机抽样法或统计试验法。计算机的发展和计算技术的提高为蒙特卡罗模拟提供了高效的计算手段，使蒙特卡罗法的应用范围越来越广。随着模拟次数的增加，蒙特卡罗法的计算结果将逐渐趋近于精确解，因此，在结构可靠度计算中，蒙特卡罗法被认为是一种准精确计算方法，而其它近似计算方法的精度也常常用蒙特卡罗法进行验证。由概率定义可知，某事件的发生概率可以用大量试验中该事件发生的频率来估算。因此在结构可靠度计算中，可以通过对随机变量进行大量的随机抽样，然后把这些抽样值一组一组地代入结构功能函数，根据计算得到的功能函数值，确定结构安全与否，最后根据事件发生的次数，计算结构的可靠度或失效概率。蒙特卡罗法1.

随机变量的取样

用MonteCarlo法分析结构可靠度问题，关键是产生已知分布变量的随机数。首先在区间(0，1)内产生均匀分布随机数，然后再变换成给定分布变量的随机数。

蒙特卡罗法(1)伪随机数的产生和检验

在区间(0，1)内产生随机数的方法一般有数学方法和物理方法两种。产生随机数有很多方法，例如取中法、加同余法、乘同余法、混合同余法和组合同余法等等。乘同余法以其统计性质好、周期长等优点而被人们广泛采用。蒙特卡罗法乘同余法的算式为

(5-41)

式中a、c和m为正整数。式(5-41)表示以m为模数的同余式，xi+1为以m除(axi+c)后得到的余数。具体计算时，可引入参数ki，即令

(5-42)蒙特卡罗法式中lnt为取整符号，式(5-41)成为：

(5-43)若将xi+1除以模数m，可得标准化的随机数ui+1

，即

(5-44)蒙特卡罗法例5–7

设a＝3，c＝1，m＝5,

试求8个随机数。解：

xi的初值取为x0=1.0利用式(5–42)得

蒙特卡罗法由式(5–43)得：x1=31+1–50=4再由式(5–44)得：u1=4/5=0.8再令i＝1得：

x2

=34+1–52=3

u1

=3/5=0.6蒙特卡罗法将产生的8个随机数列出如下：

xi12345678ui0.80.60.00.20.80.60.00.2

从所得结果可见，随机数出现周期为4的规律，这种随机数不好。产生这种规律的原因在于常数m取得太小了。这里m＝5，周期为4，周期数小于m。为了得到较好的随机数，可以将m的值取大些，并限定使用适当区间内的随机数，以避免随机数出现周期性循环。蒙特卡罗法建议对二进制电子计算机三个常数的取值如下：

a=27+1，c=1，m=235

为了判别所得随机数的优劣，一般还

应对伪随机数进行统计检验，主要是检验

其均匀性和独立性。蒙特卡罗法a.正态分布

设随机数un和un+1是(0，1)区间内的两个均匀随机数，则可通过下列变换得到服从标准正态分布N(0，1)的两个随机数和

(2)给定分布下随机变量数值的产生

蒙特卡罗法(5-45)

如果随机变量X为一般正态分布则其随机数xn和xn+1由下式计算：

(5-46)

蒙特卡罗法b.对数正态分布

设X为对数正态分布变量，其均值为mx，标准差为σx，变异系数为Vx。因为Y＝lnX为正态分布，所以，根据式(2–38)和(2–39)，得其标准差和均值分别为

蒙特卡罗法变量Y的随机数由式(5–30)和(5–31)产生后,可得到对数正态分布变量X的随机数为：

(5–47)蒙特卡罗法c.极值I型分布

极值I型分布变量的随机数一般是通过其概率分布函数得到,因此先讨论一般情况下变量随机数的产生方法。设已知随机变量X的概率分布函数为Fx(x)，则其随机数可以由下式得到：

(5–48)式中ui为区间(0，1)内的均匀随机数。蒙特卡罗法对于极值I型变量，其分布函数为：

FX(xi)=exp{–exp[–a(xi–k)]}式中α、k是与X的均值mx和标准整σx有关的常量。

设已知随机数ui，由式(5–46)得到：

ui=FX(xi)=exp{–exp[–a(xi–k)]}

蒙特卡罗法于是可得：

(5–49)利用式(2–42)，得：

α＝1.2825／σx，k＝mx–0.450σx，

代入上式得到：

xi=mX–0.45σX–0.7797σX

ln(–lnui)(5–50)蒙特卡罗法2.结构失效概率的计算

设已知统计独立的随机变量X1、X2、···、Xn，其概率密度函数分别为、、···、结构的功能函数为：用MonteCarlo法计算结构的失效概率。

蒙特卡罗法(1)首先通过随机抽样获得各变量的分位值x1、x2、···、xn，如图5-5所示。图5–5

(2)计算功能函数值Zi：

Zi=g(x1，x2，···，xn)

蒙特卡罗法(3)设抽样组数为N，每组随机变量计算得到的功能函数值为Zi，Zi≤0的次数为L，则在大批抽样之后，结构的失效概率可由下式计算：

Pf=L/N (5–51)

即：通过大批抽样得到的结构失效次数与总抽样次数之比即为结构的失效概率。蒙特卡罗法用MonteCarlo法计算结构的失效概率时，需要进行大批抽样，需要确定最低抽样数N。显然，抽样数N与计算结果的精度有关。有文献建议，以95％的置信度保证MonteCarlo法计算精度，其允许误差ε为：(5–52)由上式可见，结构计算次数N越大，误差ε越小。因此，要达到一定的精度，N必须取得足够大。蒙特卡罗法为简便起见，有文献建议N必须满足：

N≥100/Pf (5–53)式中Pf为预先估计的失效概率。在工程实际中，一般要求失效概率Pf很小，这就需要计算的次数很多。例如，工程结构的失效概率一般在0.1％以下，要求计算的次数须达十万次以上。这个要求使实际计算不是遇到困难，就是耗时过多，MonteCarlo法的使用也因此而受到限制。蒙特卡罗法例5–8

设某结构的极限状态方程为g＝R–S＝0，R和S分别为正态和极值I型分布的随机变量，其统计特性为：试用MonteCarlo法求其失效概率。

蒙特卡罗法解：计算过程如下：

1.产生区间(0，1)内的均匀分布随机数ui。

2.产生变量S和R的随机数。变量S为极值I型分布，产生随机数的公式为：

蒙特卡罗法R为正态分布，产生一组随机数的公式为：

3.把所得变量的随机数代入功能函数式g＝R–S，求得g值。

4.重复步骤1至3的计算，累积计算出现g≤0的次数L和总次数N。

蒙特卡罗法5.计算次数足够多，如N满足N≥100／Pf(Pf

是预估计的结构失效概率)，则计算结束。由Pf

＝

L／N算得结构的失效概率。前10次模拟计算所得的u、R、S、g和L值列于表5–8。得到g≤0的次数为1，于是可得结构的失效概率为：L／N＝l／10＝10％

蒙特卡罗法表5-8

前10次模拟计算的结果

NuRuSg=R–SL10.08039104.3920.687187.54216.850020.2795131.3900.490475.54355.847030.8376106.2260.427772.25733.969040.7450120.5650.264163.84456.721050.409577.4810.713589.518-12.0371蒙特卡罗法60.617086.8210.481371.76415.057070.369380.7710.189959.70621.065080.3369125.2910.492975.67549.616090.9023107.4130.738891.23416.1790100.3125128.1770.752792.75135.4260蒙特卡罗法

显然，模拟计算的次数太少，所得结果的精度达不到要求。增加模拟计算的次数以提高计算结果的精度，表5-9给出了模拟1000次计算过程中的部分结果。由表中结果可见，算到500次以后，Pf值已趋于稳定，即

Pf

＝

20％

若按此失效概率估算最少模拟次数，得到的结果是500。

本例用JC法算得的β＝0.753，对应的失效概率为：Pf

＝

22％蒙特卡罗法表5-9

模拟计算的部分结果

N

Pf

N

Pf

1000.2006000.1982000.1907000.2043000.2008000.2004000.19759000.1945000.18810000.196返回首页蒙特卡罗法5.5

响应面法

响应面法（ResponseSurfaceMethod）是用一个简单的显式函数逐步逼近实际的隐式（或显式）极限状态函数，使可靠度计算得到简化。由于该方法可以直接应用确定性结构的计算程序，使得可靠度分析工作更加简便易行，为大型复杂结构的可靠度分析展现出良好的应用前景。

1.一次响应面法

设含有两个基本变量的极限状态函数为，它一般是非线性的，取响应面函数为一次多项式，即

(5-54)

以均值为中心，在区间内选取2+1个样本点，是确定取值界限的选择参数，一般取=1。响应面法由样本点可计算得到的值，再由式（5-54）建立三个方程求解系数。响应面函数确定后，即可计算结构的可靠指标和设计验算点XD的值。然后再以XD为中心选取一组新的样本点，重复上述过程，即可得到与极限状态方程相对应的可靠指标和设计验算点的近似值。这一方法的几何表示见图5-6。

响应面法图5-6两变量的线性响应面

X2X1Rsponsesurface=a+ax+ax12021y=g(x,x)12~Trueresponse

y=g（x1，x2）Y’=g’(X1,X2)=a0+a1X1+a2X2响应面真实面X1X2响应面法为了提高响应面法计算结构可靠度的精度，对于如下的极限状态函数：

(5-55)响应面函数可取为随机变量的完全二次式：

(5-56)确定常量矩阵A、B、C需要足够多的样本点。

2.二次响应面法响应面法

在结构可靠度分析中，涉及到的随机变量往往较多，为了确定式(5-56)中的待定系数，需要取很多的样本点进行分析计算，从而影响了响应面法的计算效率。为了在保证计算精度的前提下提高计算效率，取如下形式的响应面函数：

式中的系数需由2n+1个样本点得到足够的方程来确定。

响应面法响应面函数中不含混合项，确定响应面函数所取的样本点减少。计算方法如下：

第一步是以均值为中心点，在区间(内选取样本点，有文献建议在1-3之间取值。由样本点可计算得到2n+1个函数g(X)值，然后确定响应面函数中的待定数。得到响应面函数之后，即可求出极限状态面上设计验算点的近似值XD。响应面法第二步是选取新的中心点，新中心点XM可选在均值点mX与XD的连线上，并保证满足极限方程g(X)=0，即第三步是以XM为中心点选取新的一组样本点，重复第一步的工作，即可得到极限状态面上设计验算点的值和相应的可靠指标。响应面法整个过程需求解4n+3个函数g(X)的值，这一逼近过程的图形表示见图5-7。

图5-7二次响应面法的逼近过程X1XmmX212XM2XXM1X2X1X)=0)=0Xg(=0g(X)=0g’(Xg’(X)响应面法例5-9

图5-8所示的具有随机参数的非线性振动系统，受到持时和幅值均随机变化的矩形脉冲荷载作用，设极限状态方程为：

式中为系统最大的位移响应，r为最大的线性位移。各变量均假设为统计独立的对数正态变量，统计参数列于表5-10，试用响应面法分析系统的动力可靠度。

响应面法图5-8

非线性恢复力和脉冲荷载作用下的非线性振动

响应面法表5-10

变量的统计参数序号

变量

均值

标准差

1M1.000.052C11.000.103C20.100.014r0.500.055F11.000.206t11.000.20响应面法本例用高效的抽样方法直接计算原方程的失效概率，抽样1500次得到的结果是Pf=3.81×10-2，标准误差为4.5%；若用同样的抽样方法，相同的抽样次数计算二次响应面方程的失效概率，得到的结果是Pf=3.04×10-2

，标准误差为4.9%。两种方法计算结果的精度基本一致，但计算机时后者只是前者的1/12。响应面法

图5-9为一两跨六层的平面刚架，各变量的统计特性列于表5-11，E为弹性模量，I为截面惯性矩。设刚架顶部侧向位移的极限状态函数（单位：mm）为：例5-10试用响应面法计算结构的可靠指标。

响应面法作为比较，表5-12给出了MonteCarlo法抽样200,000次得到的结果以及响应面法算得的可靠指标和失效概率。从得到的结果可见，一次、二次响应面法的计算结果和MonteCarlo法得到的结果很接近。响应面法7.5m22PP15P3PP46P3m6D（x）

图5-9平面刚架响应面法表5-11随机变量的统计参数

随机变量

分布

均值

标准差

柱子的弹模E(mpa)对数正态

200002000柱子的惯矩I(m4)

对数正态

1.0×10-2

1.0×10-3

梁的弹模E（MPa）

对数正态

200002000梁的惯矩I（m4）

对数正态

1.5×10-21.5×10-3P1（N）正态250006250P2（N）正态280007000P3（N）正态290007250P4（N）正态300007500P5（N）正态310007750P6（N）正态320008000响应面法表5-12

几种方法算得的可靠指标

方法

MonteCarlo法

2.2030.01380线性响应面法

2.2552.2360.012070.01268二次响应面法

2.2012.1030.013870.01773响应面法由例9、例10的计算结果可以看出，响应面法计算结构可靠度的效率和精度都比较好，而且适用于静力、动力、线性和非线性等各类问题的求解。返回首页响应面法5.6结构的可靠度设计

在结构可靠度设计中，结构极限状态设计应符合下列要求：

(5-59)式中为基本变量