验算点法在工程结构可靠度编程中的运用

验算点法在工程结构可靠度编程中的运用北京航空航天大学交通科学与工程学院，北京，100191）摘要：用验算点法对可靠度编程中的关键问题进行了探讨并给出了解决方法。这些问题主要包括随机变量服从正态分布的情形且功能函数为非线性和随机变量不服从正态分布时的当量正态化方法。针对这两种情况运用了两个编程算例来说明验算点法在可靠度分析中的运用。关键字：验算点法、可靠指标、正态分布函数、非正态随机变量、当量正态化1.前言结构可靠度为结构在规定的时间内和规定的条件下完成预定功能的概率。结构可靠性理论的研究，起源于对结构设计、施工和使用过程中存在的不确定性的认识，以及结构设计风险决策理论中计算结构失效概率的需要。结构可靠度计算方法有一次二阶矩方法、二次二阶矩方法、蒙特卡洛方法及其他方法。一次二阶矩方法又分为中心点法和验算点法，其中验算点法是目前可靠度分析中最常用的方法。由于这两种方法都是将非线性功能函数作为一次泰勒级数展开，并使用了随机变量的平均值（一阶矩）和方差（二阶矩），故称为一次二阶矩方法。利用验算点法计算结构的可靠指标时，需要预先知道验算点的坐标值，而对于非线性结构功能函数和非正态随机变量的情形，验算点坐标值是不能预先求得的，因此一般需要迭代求解。2.随机变量服从正态分布的情形功能函数为线性函数图1一个随机变量时的可靠指标（左图为正态随机变量，右图为标准正态随机变量）假定存在n个相互独立的随机变量，X,X，…,X，其均值为卩,卩，…，卩，标准差为12n X1X2 Xn◎Q,...,◎结构功能函数为：XX X

Z=g(X,X,X...X)=a+"aXx123n0iii=1其中aC=l,2,3...n)为常数i将随机变量XC=1,2,3...n)变换为标准正态随机变量Y(i=1,2,3...n)ii(i=1,2,3...n(i=1,2,3...n)2）Y=匚XiGXi则由（1）表示的功能函数表示成Z=g&,Y,Y...Y)=a+》a(pZ=g&,Y,Y...Y)=a+》a(p+gY)=a+工Y1 2 3n0 iXi

i=1从而功能函数的平均值和标准差表示为Xi

iap+0 iXii=1 iaGYiXinp=a+yapZ0 iXii=1 i按照严格的可靠度指标定义i=1a+X

卩=匕=丄

G |'Vz，乙 a2G2ix.i=1可靠度指标和结构失效概率存在精确的对应关系P=①（-B）apiXi■i=^-3）Z=y对极限状态方程 Z=ao+i=1ap+XagY=0iXiii=1a2iG2x得到:ii=1yaGYiXyaGYiXii=1+yapiXi込a2iG2x.i=1 ,a0；工a2g2Vi=1ixi=0与公式（与公式（3）比较，有iXii=4—2a2 2X1i=1—卩=05）令a=cos0YY

iiaQiXii=4—2a2 2X1i=1—卩=05）令a=cos0YY

iiaQ(i=1,2,3,...n)i=1公式（5）可以写成：工aY-B=Xcos0Y-P=0Y Yiiii=1 i=16）公式（6）表示的是一法线式的直线方程，cos0为法线与坐标轴夹角余弦Yi=1cos0 =1Yi验算点在Y空间（标准正态空间）表示为:在X空间表示为：x*=（x*,x*,x*...x*）1 2 3ny*=(y*,y*,y*...y*)1 2 3n两者之间的关系为：X*=卩+Qy*i X Xiiii=1,2,3,...n)根据几何关系有：y*=Pa=Pcos0i Y Yiii=1,2,3,...n)在X空间，验算点坐标值：x*=卩+QiXiPa=卩+Pqcos0Y X Xi i i通常表示为：X*=卩+QPa=p+Pqcos0i X XX X X Xi i i i i iYii=1,2,3,...n)=1,2,3,...n)功能函数为线性函数假定随机变量X1，X2,…,Xn服从正态分布’但结构功能函数不再是线性函数’显然,这时精确求解Z的平均值和标准差是非常困然的。同结构功能函数为非线性的情形一样，如果将可靠指标定义为标准正态坐标系中坐标原点到极限状态曲面的距离，垂足为验算点则不管结构极限状态方程的数学表达形式如何，只要具有相同的力学或物理含义，在标准正态坐标系中，所表示的都将是同一个曲面，曲面上与坐标原点距离最近的点也只有一个。因而，所得到的可靠指标是唯一的，不像中心点法那样，随结构极限状态方程数学表达式的形i=1i其均值和标准差为：|lx=EZZLL)+Y£gJx*Laxi=1i其均值和标准差为：|lx=EZZLL)+Y£gJx*Laxi=1=g**,x,x...x)+^agX1 23n aXi=1 i——(7./V ，X1 2 3n-x*)i-x*)Xiio= -EZ)=[另另aggC*)agC*)=-—i=1j=1axiPagC)oi厶——X O\i ax Xi\i=11-i=1—axihE[(X-EX)&-EX)]ax iii：j所以可靠度指标：22卩=戈\LI计aXi实际上验算点不可知，需要补充条件：x*=卩+Pacos0

iXi对比表达式得到i=1,2,3,...n)cos0X

i,■'!=idgC*)xaSX 「dgC*xaG=1,2,3..n)dXj随机变量X1，x随机变量X1，x2的平均值和标准差分别为卩X1=38,aX1=3.8卩=54，a=5.4，X2 X2均服从正态分布。用验算点法编程算例假定结构功能函数为Z=g（X,X）=XX-1000。X2.3.1算法分析就算可靠指标卩和失效概率P，允许误差取£=10-2.3.1算法分析1）假定验算点，一般取x*（o）=（卩，卩，…,卩）X1 X2 Xn所以x*（o）=卩=38，1X1x*(o)=卩=542）计算卩ZL_aZLgC*X1-x*Xiii=1二dg(x*)\o"CurrentDocument"XadX (i -dg(x*)功能函数对X,X的一阶偏导数为］ =X*,1 2 dX 21i=1dg(x*)X x*dX -12所以卩=x*卩+x*卩-x*x*-100054x*+38x*-x*x*-1000cos0X1cos0X21Y 2X 1 23)计算cos0,cos0X1dg(x*)Xc叫jdg(x*)―X cax x.;1- jag(x*)xcax x\o"CurrentDocument"dg(x*)—x cax x.;1- jX24)x*ij=1j=1重新计算验算点=卩+0QCOS0XX

ii2X1x*c2X)2+(x*c13.8x*(i=1,2,3,...n)5.4x*、1一.4x*)1则x*=卩+0QCOS0,1X1 X1 X1x*=卩+0ccos02 X2 X2 X25)若x*(1)-x*(0)<£11£为规定的允许误差，则停止迭代，所求0即为要求的可靠指标；否则，取x*(0)=x*(1),转2)继续迭代，当验算点误差小于£，结束。112.3.2源程序#include<iostream>#include<cmath>usingnamespacestd;//example3-4doubleBeta(double&x1,double&x2);doubleCosOx1(double&x1,double&x2);doubleCosOx2(double&x1,double&x2);doubleX1=38.0;doubleX2=54.0;doublee=1.0e-3;doubleBeta(double&x1,double&x2){return(54.0*x1+38.0*x2-x1*x2-1000.0)/sqrt(pow(3.8*x2,2)+pow(5.4*x1,2));doubleCosOx1(double&x1,double&x2){return-3.8*x2/sqrt(pow(3.8*x2,2)+pow(5.4*x1,2));}doubleCosOx2(double&x1,double&x2){return-5.4*x1/sqrt(pow(3.8*x2,2)+pow(5.4*x1,2));}intmain(){doublex1=X1;doublex2=X2;doublebeta;doublefai;inti=0;do{i++;X1=x1;X2=x2;beta=Beta(X1,X2);x1=38.0+3.8*beta*CosOx1(X1,X2);x2=54.0+5.4*beta*CosOx2(X1,X2);}while(sqrt(pow(x1-X1,2)+pow(x2-X2,2))>e);cout<<"Example3-4"<<'\n'<<"共"<<i<<"次迭代"<<'\n'<<"B=<<beta<<'\n'<<"x1="<<x1<<'\n'<<"x2="<<x2<<'\n'<<endl;return0;}2.3.3运行结果分析3.随机变量不服从正态分布的情形3.1当量正态化在实际工程中，许多随机变量并不一定服从正态分布，如有的变量服从对数正态分布，有的服从极值I型分布。这样，需要研究随机变量不服从正态分布时的可靠指标计算方法所谓当量正态化，就是将不服从正态分布的随机变量X•等效为正态随机变量X',当量正态ii

化的条件是，在验算点处使得非正态随机变量X的概率分布函数值与当量正态随机变量X'ii的概率分布函数值相等，X的概率密度函数值与X'的概率密度函数值相等，用公式表示为:iifC*)=eXii(x*)=-^pXi fC*)=eXii(x*)=-^pXi icx*—卩x'―i icX'丿ix*—px'

―i icv X'丿i=F (x*)X'ii(x*)X'ii7)8)X'i9)10)cos0X'ij=1dg(x*)XcdX x'idg(x*)―X cX'

j6Xjx*=x*=piXi3.2编程算例+pccos0XXii(i=1,2,3,...n)某钢筋混凝土梁，计算跨度l0某钢筋混凝土梁，计算跨度l0=5.65m，假定均布荷载q的平均值卩S119.71KN.m标准差c=23.94KN.m，服从极值I型分布；R服从对数正态分布，卩=284.577KN.mSR变异系数5=0.082KN.m。用当量正态化方法求梁的可靠指标和失效概率R3.2.1算法分析1)假定验算点，一般取x1)假定验算点，一般取x*(0)=(pX1,pX2,...,pXn取r*(0)=p， s*(0)=pRS2)计算p，QX' X'ii对于服从对数正态分布的情况，由式"lnx*—p'7)得：InX(、x*—pi InX'QlnX'iQ、 lnXi由式(8)可得:- ^ln2)计算p，QX' X'ii对于服从对数正态分布的情况，由式"lnx*—p'7)得：InX(、x*—pi InX'QlnX'iQ、 lnXi由式(8)可得:- ^lnx*—p、QlnXiInXix*—p=i InX'QlnX'ix*

lnXii将(11)代入(12)QlnXi得:InXiJ2兀QlnX'i厂x*—piInX'QlnX'i12)Q=x*QX'ilnXXi i将(13)代入13)11)得：p=x*+x*Clnx*+pX'ii iilnXi)=x*i因为R服从对数正态分布，所以plnR=224.577=5.64V1+82 <1+0.0822 •+0.0822)=0.082lnR1一lnx*+lnip=r*G—lnr*+p)=r*(6.648—lnr*),Q=r*Q =0.082r*R' lnR R' lnR荷载效应S服从极值I型分布，其概率分布函数的参数为:兀 3.14cc— 0.5772「CFcc—a= —= =0.054u=p— =119.71—0.054辰 V6x23.94 SaS3.14108.93令t=exp—aC*—u刀=exp[—0.054C*—108.93”则s的概率密度函数和概率分布函数分别为:f(s*)=SF(s*)=exp(—t)Saexp—a(s*—u)—t=0.054exp[—0.054(s*_108.93)—t所以，py=s*-帖LFs$一irfC*)]}

fC*)S当中涉及到求正态分布函数的反函数，采用以下方法e-1(x)=U，则正态分布函数的反函数是：a式中U的近似值为：a匕 1L：0<fl<Q.534取3’l\=y0 当出二山时 ®/3=1・『・心a>0.5时H-1-«10 丄址LvS岭）2y=-In[4卩（1-卩）]式中 式中 b0=O.157O796288xlOi,b3=-0.2250947176xl0-3,b6=-0.1045274970x10-5,b9=0.3657763036x10-i0,b1=0.3706987906x10-i,b4=0.6841218299x10-5,b7=0.8360937017x10-7,b10=0.6936233982xl0-i2.b2=-0.8364353589x10-3 ,b5= 0.5824238515x10-5,b8=-0.3231081277xl0-8,83)计算卩结构功能函数Z=R-S=R—S',从而得到可靠指标的表达式P—卩R S—+a2S'4)计算cos。Ricos4)计算cos。Ricos0Siicos0Ria—R—v;a2+a2卞R S',cos0a2+a2TRS'5)重新计算验算点+cos0paS'S'ir*=+cos0paS'S'iR RRi6)若x*(1)—x*(0)<8，8为规定的允许误差，则停止迭代，所求P即为要求的可靠指标；否则，取x*（0）=x*（1）,转2）继续迭代，当验算点误差小于8，结束113.2.2源程序#include<iostream>#include<cmath>#include<vector>#include<iomanip>usingnamespacestd;#definePI3.14159265358979//example3-5doubleCumulative(doublez);doubleFxi(doublexi);doubleuR1(doubler0);doubleQR1(doubler0);doubleT(doubles0);doublefS(doubles0,doublet);doubleFS(doubles0);doubleQS1(doubles0);doubleuS1(doubles0);doubleBeta(doubleur1,doubleus1,doubleqr1,doubleqs1);doubleAlphaR1(doubleqr1,doubleqs1);doubleAlphaS1(doubleqr1,doubleqs1);doubleLittleFai(doublex);doubleBigFai(doublex);doubleuR=284.577;doubleuS=119.70;doubleCumulative(doublez){if(z>6.0)return1.0;if(z<-6.0)return0.0;doubleb1=0.31938153;doubleb2=-0.356563782;doubleb3=1.781477937;doubleb4=-1.821255978;doubleb5=1.330274429;doublep=0.2316419;doublec2=0.3989423;doublea=fabs(z);doublet=1.0/(1.0+a*p);doubleb=c2*exp(-z*z/2.0);doublen=((((b5*t+b4)*t+b3)*t+b2)*t+b1)*t;n=b*n;returnn;}doubleuR1(doubler0){returnr0*(6.648-log(r0));doubleQR1(doubler0){return0.082*r0;}doubleT(doubles0){returnexp(-0.054*(s0-108.93));}doublefS(doubles0,doublet){return0.054*exp(-0.054*(s0-108.93)-t);}doubleFS(doublet){returnexp(-t);}doubleQS1(doubles0){returnLittleFai(BigFai(T(s0)))/fS(s0,T(s0));}doubleuS1(doubles0){returns0-BigFai(FS(T(s0)))*QS1(s0);}doubleBeta(doubleur1,doubleus1,doubleqr1,doubleqs1){return(ur1-us1)/sqrt(qr1*qr1+qs1*qs1);}doubleAlphaR1(doubleqr1,doubleqs1){return-qr1/sqrt(qr1*qr1+qs1*qs1);}doubleAlphaS1(doubleqr1,doubleqs1){returnqs1/sqrt(qr1*qr1+qs1*qs1);doubleLittleFai(doublex){returnexp(-x*x/2.0)/sqrt(2.0*PI);}doubleBigFai(doublea){doubleb;doubley;doubleU;vector<double>Barray;Barray.push_back(0.1570796288e1);Barray.push_back(0.3706987906e-1);Barray.push_back(-0.8364353589e-3);Barray.push_back(-0.2250947176e-3);Barray.push_back(0.6841218299e-5);Barray.push_back(0.5824238515e-5);Barray.push_back(-0.1045274970e-5);Barray.push_back(0.8360937017e-7);Barray.push_back(-0.3231081277e-8);Barray.push_back(0.3657763036e-10);Barray.push_back(0.6936233982e-12);if(a>0&&a<0.5)b=a;if(a==0.5)b=1-a;if(a>0.5)