驗(yàn)算點(diǎn)法在工程結(jié)構(gòu)可靠度編程中的運(yùn)用

驗(yàn)算點(diǎn)法在工程結(jié)構(gòu)可靠度編程中的運(yùn)用（北京航空航天大學(xué)交通科學(xué)與工程學(xué)院，北京，100191）摘要：用驗(yàn)算點(diǎn)法對(duì)可靠度編程中的關(guān)鍵問題進(jìn)行了探討并給出了解決方法。這些問題主要包括隨機(jī)變量服從正態(tài)分布的情形且功能函數(shù)為非線性和隨機(jī)變量不服從正態(tài)分布時(shí)的當(dāng)量正態(tài)化方法。針對(duì)這兩種情況運(yùn)用了兩個(gè)編程算例來說明驗(yàn)算點(diǎn)法在可靠度分析中的運(yùn)用。關(guān)鍵字：驗(yàn)算點(diǎn)法、可靠指標(biāo)、正態(tài)分布函數(shù)、非正態(tài)隨機(jī)變量、當(dāng)量正態(tài)化前言結(jié)構(gòu)可靠度為結(jié)構(gòu)在規(guī)定的時(shí)間內(nèi)和規(guī)定的條件下完成預(yù)定功能的概率。結(jié)構(gòu)可靠性理論的研究，起源于對(duì)結(jié)構(gòu)設(shè)計(jì)、施工和使用過程中存在的不確定性的認(rèn)識(shí)，以及結(jié)構(gòu)設(shè)計(jì)風(fēng)險(xiǎn)決策理論中計(jì)算結(jié)構(gòu)失效概率的需要。結(jié)構(gòu)可靠度計(jì)算方法有一次二階矩方法、二次二階矩方法、蒙特卡洛方法及其他方法。一次二階矩方法又分為中心點(diǎn)法和驗(yàn)算點(diǎn)法，其中驗(yàn)算點(diǎn)法是目前可靠度分析中最常用的方法。由于這兩種方法都是將非線性功能函數(shù)作為一次泰勒級(jí)數(shù)展開，并使用了隨機(jī)變量的平均值（一階矩）和方差（二階矩），故稱為一次二階矩方法。利用驗(yàn)算點(diǎn)法計(jì)算結(jié)構(gòu)的可靠指標(biāo)時(shí)，需要預(yù)先知道驗(yàn)算點(diǎn)的坐標(biāo)值，而對(duì)于非線性結(jié)構(gòu)功能函數(shù)和非正態(tài)隨機(jī)變量的情形，驗(yàn)算點(diǎn)坐標(biāo)值是不能預(yù)先求得的，因此一般需要迭代求解。隨機(jī)變量服從正態(tài)分布的情形功能函數(shù)為線性函數(shù)功能函數(shù)隨機(jī)變量是一個(gè)正態(tài)隨機(jī)變量，其概率密度函數(shù)和U的密度曲線如圖1示圖1一個(gè)隨機(jī)變量時(shí)的可靠指標(biāo)（左圖為正態(tài)隨機(jī)變量，右圖為標(biāo)準(zhǔn)正態(tài)隨機(jī)變量）假定存在n個(gè)相互獨(dú)立的隨機(jī)變量，，其均值為，標(biāo)準(zhǔn)差為結(jié)構(gòu)功能函數(shù)為：（1）其中為常數(shù)將隨機(jī)變量變換為標(biāo)準(zhǔn)正態(tài)隨機(jī)變量（2）則由（1）表示的功能函數(shù)表示成從而功能函數(shù)的平均值和標(biāo)準(zhǔn)差表示為按照嚴(yán)格的可靠度指標(biāo)定義（3）可靠度指標(biāo)和結(jié)構(gòu)失效概率存在精確的對(duì)應(yīng)關(guān)系對(duì)極限狀態(tài)方程兩端同時(shí)除以得到：（4）與公式（3）比較，有（5）令公式（5）可以寫成：（6）公式（6）表示的是一法線式的直線方程，為法線與坐標(biāo)軸夾角余弦圖2可靠度指標(biāo)的幾何意義及驗(yàn)算點(diǎn)驗(yàn)算點(diǎn)在空間（標(biāo)準(zhǔn)正態(tài)空間）表示為：在空間表示為：兩者之間的關(guān)系為：根據(jù)幾何關(guān)系有：在空間，驗(yàn)算點(diǎn)坐標(biāo)值：通常表示為：功能函數(shù)為線性函數(shù)假定隨機(jī)變量服從正態(tài)分布，但結(jié)構(gòu)功能函數(shù)不再是線性函數(shù)，顯然，這時(shí)精確求解的平均值和標(biāo)準(zhǔn)差是非常困然的。同結(jié)構(gòu)功能函數(shù)為非線性的情形一樣，如果將可靠指標(biāo)定義為標(biāo)準(zhǔn)正態(tài)坐標(biāo)系中坐標(biāo)原點(diǎn)到極限狀態(tài)曲面的距離，垂足為驗(yàn)算點(diǎn)，則不管結(jié)構(gòu)極限狀態(tài)方程的數(shù)學(xué)表達(dá)形式如何，只要具有相同的力學(xué)或物理含義，在標(biāo)準(zhǔn)正態(tài)坐標(biāo)系中，所表示的都將是同一個(gè)曲面，曲面上與坐標(biāo)原點(diǎn)距離最近的點(diǎn)也只有一個(gè)。因而，所得到的可靠指標(biāo)是唯一的，不像中心點(diǎn)法那樣，隨結(jié)構(gòu)極限狀態(tài)方程數(shù)學(xué)表達(dá)式的形式而變。圖3驗(yàn)算點(diǎn)取法如果驗(yàn)算點(diǎn)已知可以在該點(diǎn)一次項(xiàng)展開：其均值和標(biāo)準(zhǔn)差為：所以可靠度指標(biāo)：實(shí)際上驗(yàn)算點(diǎn)不可知，需要補(bǔ)充條件：對(duì)比表達(dá)式得到：編程算例假定結(jié)構(gòu)功能函數(shù)為。隨機(jī)變量的平均值和標(biāo)準(zhǔn)差分別為，，，，均服從正態(tài)分布。用驗(yàn)算點(diǎn)法就算可靠指標(biāo)和失效概率，允許誤差取2.3.1算法分析假定驗(yàn)算點(diǎn)，一般取所以，計(jì)算功能函數(shù)對(duì)的一階偏導(dǎo)數(shù)為，所以計(jì)算重新計(jì)算驗(yàn)算點(diǎn)則，若，為規(guī)定的允許誤差，則停止迭代，所求即為要求的可靠指標(biāo)；否則，取，轉(zhuǎn)2）繼續(xù)迭代，當(dāng)驗(yàn)算點(diǎn)誤差小于，結(jié)束。2.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' <<"β="<<beta<<'\n' <<"x1="<<x1<<'\n' <<"x2="<<x2<<'\n' <<endl; return0;}2.3.3運(yùn)行結(jié)果分析隨機(jī)變量不服從正態(tài)分布的情形3.1當(dāng)量正態(tài)化在實(shí)際工程中，許多隨機(jī)變量并不一定服從正態(tài)分布，如有的變量服從對(duì)數(shù)正態(tài)分布，有的服從極值I型分布。這樣，需要研究隨機(jī)變量不服從正態(tài)分布時(shí)的可靠指標(biāo)計(jì)算方法。所謂當(dāng)量正態(tài)化，就是將不服從正態(tài)分布的隨機(jī)變量等效為正態(tài)隨機(jī)變量，當(dāng)量正態(tài)化的條件是，在驗(yàn)算點(diǎn)處使得非正態(tài)隨機(jī)變量的概率分布函數(shù)值與當(dāng)量正態(tài)隨機(jī)變量的概率分布函數(shù)值相等，的概率密度函數(shù)值與的概率密度函數(shù)值相等，用公式表示為：（7）（8）由公式（7）（8）可解得：（9）（10）當(dāng)量正態(tài)化后：3.2編程算例某鋼筋混凝土梁，計(jì)算跨度，假定均布荷載的平均值，標(biāo)準(zhǔn)差，服從極值I型分布；R服從對(duì)數(shù)正態(tài)分布，變異系數(shù)。用當(dāng)量正態(tài)化方法求梁的可靠指標(biāo)和失效概率3.2.1算法分析假定驗(yàn)算點(diǎn)，一般取取，計(jì)算，對(duì)于服從對(duì)數(shù)正態(tài)分布的情況，由式（7）得：即（11）由式（8）可得：（12）將（11）代入（12）得：（13）將（13）代入（11）得：因?yàn)榉膶?duì)數(shù)正態(tài)分布，所以,荷載效應(yīng)服從極值I型分布，其概率分布函數(shù)的參數(shù)為：,令則的概率密度函數(shù)和概率分布函數(shù)分別為：所以，，當(dāng)中涉及到求正態(tài)分布函數(shù)的反函數(shù)，采用以下方法：，則正態(tài)分布函數(shù)的反函數(shù)是：式中的近似值為：式中b0=0.1570796288×101，b1=0.3706987906×10-1，b2=-0.8364353589×10-3，b3=-0.2250947176×10-3，b4=0.6841218299×10-5，b5=0.5824238515×10-5，b6=-0.1045274970×10-5，b7=0.8360937017×10-7，b8=-0.3231081277×10-8，b9=0.3657763036×10-10，b10=0.6936233982×10-12.計(jì)算結(jié)構(gòu)功能函數(shù),從而得到可靠指標(biāo)的表達(dá)式計(jì)算,重新計(jì)算驗(yàn)算點(diǎn),若，為規(guī)定的允許誤差，則停止迭代，所求即為要求的可靠指標(biāo)；否則，取，轉(zhuǎn)2）繼續(xù)迭代，當(dāng)驗(yàn)算點(diǎn)誤差小于，結(jié)束3.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) b=1-a; y=-log(4*b*(1-b)); doubletemp=0.0; for(inti=0;i<=10;i++) temp+=Barray[i]*pow(y,i); U=sqrt(y*temp); returnU;}intmain(){ doubleR0; doubleS0; doubleR1=uR; doubleS1=uS; doublee=1.0e-3; doublebeta; inti=1; do { i++; R0=R1; S0=S1; beta=Beta(uR1(R0),uS1(S0),QR1(R0