工程结构可靠度

结构可靠度分析与设计的编程实践摘要：本文基于结构可靠度的分析方法（改进一次二阶矩法，JC法，蒙特卡洛法）的基本原理，采取不同的算例利用MATLAB软件编制了相应计算程序分析结构可靠度。同时，考虑随机变量的分布信息，还对结构构件可靠度校准及分项系数校准等进行了简单分析。通过此次编程实践，加深了对结构可靠度重要性的理解。1. 引言 结构的可靠度是结构在规定时间内，规定条件下结构能够完成预定功能的概率。由于材料生产的不均匀性和结构作用荷载的随机性，导致抗力和荷载都属于符合某种情形的分布，通过概率论知识并依托Matlab数学分析软件我们可以计算出结构的可靠指标和失效概率，确定结构的安全程度。从简单到复杂或精确程度的不同，先后提出的可靠度计算方法有一次二阶矩方法，高次二阶矩方法，蒙特卡洛方法及其他方法。一次二阶矩又分为中心点法和验算点法，其中验算点法是目前可靠度分析最常用的方法。2. 结构可靠度分析的基本原理结构的可靠度通常受各种荷载、材料性能、几何参数、计算公式精确性等因素的影响。这些因素一般具有随机性，特别是当结构可能受到多种因素的影响，而且每一种影响不一定完全服从正态分布时，需要预先对它们进行当量正态化处理。用失效概率来描述结构的可靠度，在理论上是合理的，但在具体计算上又非常困难，因此，必须寻求近似的或简化的计算方法，以达到实用的目的。在实际工程中，结构的功能函数往往是由多个随机变量组成的非线性函数，而且这些随机变量并不都服从正态分布或对数正态分布，因此不能直接采用前述公式计算可靠指标，而需要作出某些近似简化后进行计算。下面本文将介绍分析结构可靠度的几种常用方法。2.1 中心点法若结构功能函数为线性功能函数，即式中：为已知常数（i=1,2,3n）。由此可知功能函数的统计参数为根据概率论中心极限定理，当随机变量的数量n较大时，可以认为Z近似服从于正态分布，则可靠指标直接按下式计算然后再求得得结构的失效概率：若结构的功能函数为非线性情况，即将Z在随机变量的平均值（即中心点）处按泰勒级数展开，并仅取线性项，即接下来的算法就参照线性功能函数的情况。2.2 验算点法JC法是国际“结构安全度联合委员会（JCSS）”推荐的一种计算结构可靠指标更为一般的方法。起主要特点是，对于非线性的功能函数，线性化近似不是选在中心点处，而是选在失效边界上，即以通过极限状态方程上得某一点 的切平面做线性近似，以提高可靠指标的计算精度。这里特定的点称为设计验算点，它与结构最大可能的失效概率相对应，因此也叫验算点法。2.2.1 当功能函数为线性函数时：式中：为已知常数（i=1,2,3n）。为进一步在标准正太坐标系中研究可靠指标的几何意义，按下式将随机变量变换为标准正态随机变量则结构功能函数可由表示为：从而功能函数的平均值和标准差为：这样可靠度指标为：由上式计算的可靠指标与结构失效概率的对应关系：2.2.2 当功能函数为非线性函数时： 为计算可靠指标，将非线性功能函数在验算点处泰勒级数展开，则结构功能函数的一次展开式为：则Z的平均值和标准差可分别近似表示为由此可求得可靠指标：2.3 JC法当随机变量不服从正态分布时，则需要将其进行当量正态化将随机变量等效为正态随机变量，当量正态化的条件是，在验算点处使非正态随机变量的概率分布函数值与当量正态随机变量的概率分布函数值相等，的概率密度函数值与当量正态随机变量的概率密度函数值相等，即则：当量正态化后：上述三式构成非正态随机变量情况下可靠指标的迭代计算公式。3. 结构可靠性设计的基本原理结构设计首先要保证结构的安全性，其次是保证结构的适用性、耐久性和偶然作用下的整体稳定性。为实现这些目的，就要对结构进行合理的设计。结构的设计由最初的完全经验设计法，逐步向以科学试验和理论计算为基础的方向发展。从力学计算理论的发展过程看，经历了容许应力设计法，破损阶段设计法和极限状态设计法；从概率理论的应用过程看，经历了定制设计法、半概率设计法和近似概率设计法。现在。国际上结构设计理论的发展趋势是概率极限状态设计法。3.1 结构可靠度校准在现实工程中，按照结构设计规范规定的设计表达式进行可靠度计算，按规范设计的结构或结构构件的可靠度指标。在这种情况下反映的是一种结构或一类结构构件的可靠指标，代表了设计规范的可靠水平，这一过程称为结构可靠度校准。结构可靠度校准法计算步骤如下：（1） 假设一个荷载效应比；（2） 确定构件抗力的特征值：（3） 确定基本变量的均值和标准差： 均值： 标准差： （4） 确定极限状态方程：（5） 利用PF法求解可靠度指标；（6） 变化荷载效应比，计算不同可靠指标的均值3.2 分项系数的校准在结构设计规范的校准方法中，如果给定了结构或构件的设计表达式且已知变量的概率分布和统计参数，即可求得按设计表达式设计的结构或构件的可靠指标，可靠指标至于作用效应标准值的比值有关，而与作用效应标准值本身的大小无关。这样可以通过先设定作用分项系数、抗力分项系数和作用组合系数，然后计算相应的可靠指标，并与目标可靠指标进行比较，从而确定一组使计算的可靠指标与目标可靠指标最接近的系数值，即确定的、和使下式最小下面是用迭代算法求分项系数的过程：（1）确定极限状态方程和设计方程的表达式。确定基本变量的概率分布类型和近似参数最多有两个未知的均值需要求解。一个是抗力的均值，另一个是相应的荷载效应均值，荷载效应比率用于建立荷载与抗力均值之间的联系第一次迭代，利用极限状态方程 获得两个未知变量均值的关系。（2）利用上述均值获得初始的设计点。（3）对于非正态分布对应的设计点取值，利用当量正态化公式确定等效的正态均值和正态标准差。（4）计算n个变量的方向余弦（5）计算n个变量设计点的坐标值（6）通过求解极限状态方程确定两个未知均值新的关系（7）重复3-6步直到收敛（8）当收敛条件满足，计算分项系数4. 例题的程序运行结果展示4.1 Project2.1的结果展示112页的例题5.4迭代次数12345672.57823.28963.21353.18213.18063.18053.1805-0.5998-0.4552-0.2694-0.1991-0.1847-0.1814-0.1807-0.7998-0.8895-0.9621-0.9793-0.9822-0.9828-0.98290.02480.03890.04150.03630.03460.03430.0342142页的作业5.3迭代次数123.15373.1537-0.9997-0.99970.02320.0232 4.2 Project2.2的结果展示 123页的例题5.9迭代次数123454.39033.94763.76743.76333.7633R123.0790156.9136167.5417168.4619168.4997S123.0790156.9136167.5417168.4619168.4997 127页的例题5.11迭代次数1234564.55184.39194.10824.09834.09824.098210010010010010010023.429.832.432.732.732.72342.32977.53.23883271.33273.23273.3142页的作业5.4迭代次数123452.38232.24682.21392.21322.2132X28.255431.503332.118332.176932.1816Y9.418510.501110.706110.725610.72724.3 Project2.3的结果展示127页的例题5.1102510131906.6667e-0061.2500e-0052.0000e-0052.1667e-0052.3750e-0050.70710.44720.31620.27730.22944.4 Project2.4的结果展示迭代次数12345673.21654.02323.83643.80823.80653.80643.80644.5 Project2.5的结果展示迭代次数1234-0.6402-0.7964-0.8000-0.80001.31581.21600.76820.60480.60000.60000.80790.76110.76000.76001.27661.21771.21601.21605. 结论与讨论结构可靠度设计与分析方法，经历了由简单到复杂，有经验到理论的发展过程，其设计原理更加合理，更加贴近实际。本次编程实践加深了我对结构可靠度分析方法与设计原理的认识。也意识到现阶段还有一些方法与理论不尽合理，仍需进一步完善。6. 附件附件1：演算点法程序function fff(varargin) syms X K=length(varargin); Z1=input(请输入功能函数Z=,s); for n=1:K X(n)=varargin(n); A1(n)=input(请输入平均值=); B1(n)=input(请输入标准差=); end frag=input(=); A2=A1; num=1; while num=1 P1=Z1;T1=0; for n=1:K Q1=diff(Z1,X(n)*(A1(n)-A2(n); P1=Q1+P1; T1=T1+(diff(Z1,X(n)*B1(n)2; end W1=P1/sqrt(T1); WW=subs(W1,X,A2); for n=1:K H1(n)=-diff(Z1,X(n)*B1(n)/sqrt(T1); C1(n)=A1(n)+W1*B1(n)*H1(n); end H2=subs(H1,X,A2); CC=subs(C1,X,A2); R1=0; for n=1:K R1=R1+(CC(n)-A2(n)2; end if sqrt(R1)frag num=1; for n=1:K A2(n)=CC(n); end else num=0; end disp(WW); disp(H2); disp(CC); disp(sqrt(R1); disp( ); endJC法程序function fff2(varargin) syms X K=length(varargin); Z1=input(请输入功能函数Z=,s); for n=1:K X(n)=varargin(n); A1(n)=input(请输入平均值=); B1(n)=input(请输入标准差=); TT(n)=input(如果此变量服从正态分布请输入0，对数分布输入1,极值型分布输入2 ); if TT(n)=1 a1(n)=log(A1(n)/sqrt(1+(B1(n)/A1(n)2); b1(n)=sqrt(log(1+(B1(n)/A1(n)2); else if TT(n)=2 aa(n)=pi/B1(n)/sqrt(6); uu(n)=A1(n)-0.5772/aa(n); end end end frag=input(请输入允许误差=); A2=A1; num=1; while num=1 for n=1:K if TT(n)=0 a2(n)=A2(n); b2(n)=B1(n); else if TT(n)=1 a2(n)=A2(n)+A2(n)*(-log(A2(n)+a1(n); b2(n)=A2(n)*b1(n); else if TT(n)=2 tt(n)=exp(-aa(n)*(A2(n)-uu(n); ff(n)=aa(n)*exp(-aa(n)*(A2(n)-uu(n)-tt(n); FF(n)=exp(-tt(n); b2(n)=normpdf(norminv(FF(n),0,1),0,1)/ff(n); a2(n)=A2(n)-norminv(FF(n),0,1)*b2(n); end end end end P1=Z1;T1=0; for n=1:K Q1=diff(Z1,X(n)*(a2(n)-A2(n); P1=Q1+P1; T1=T1+(diff(Z1,X(n)*b2(n)2; end W1=P1/sqrt(T1); WW=subs(W1,X,A2); for n=1:K H1(n)=-diff(Z1,X(n)*b2(n)/sqrt(T1); C1(n)=a2(n)+W1*b2(n)*H1(n); end CC=subs(C1,X,A2); R1=0; for n=1:K R1=R1+(CC(n)-A2(n)2; end if sqrt(R1)frag num=1; for n=1:K A2(n)=CC(n); end else num=0; end disp(WW); disp(CC); disp(sqrt(R1); disp( );end附件2： 蒙特卡洛法程序function fff3(varargin) syms X K=length(varargin); Z1=input(请输入功能函数Z=,s); for n=1:K X(n)=varargin(n); A1(n)=input(请输入平均值=); B1(n)=input(请输入标准差=); TT(n)=input(如果此变量服从正态分布请输入0，对数分布请输入1,极值型分布输入2 ); end frag=input(请输入模拟精度：变异系数=); N=input(估计需要模拟次数N=); M=1; for n=1:K if TT(n)=0 R1=rand(1,N); X1(n,:)=A1(n)+B1(n)*norminv(R1,0,1); else if TT(n)=1 R2=rand(1,N); a1(n)=log(A1(n)/sqrt(1+(B1(n)/A1(n)2); b1(n)=sqrt(log(1+(B1(n)/A1(n)2); X1(n,:)=exp(a1(n)+b1(n)*norminv(R2,0,1); else TT(n)=2; R3=rand(1,N); aa(n)=pi/B1(n)/sqrt(6); uu(n)=A1(n)-0.5772/aa(n); X1(n,:)=uu(n)-log(-log(R3)/aa(n); end end end Nf=0; for n=1:N ww(n)=subs(Z1,X,(X1(:,n); if ww(n)0.0001) zz=0; for i=1:n if u(i)=1 Fq(i)=Extreamcdf(y(i),mu(i),sigma(i); fq(i)=Extreampdf(y(i),mu(i),sigma(i); b1(i)=normpdf(norminv(Fq(i),0,1),0,1)./fq(i); j(i)=y(i)-norminv(Fq(i),0,1)*b1(i); end if u(i)=2 if y(i)0 b1(i)=y(i)*sigma(i); j(i)=y(i)+y(i)*(-log(y(i)+mu(i); else b1(i)=sigma(i)/2; j(i)=mu(i); end end if u(i)=3 b1(i)=sigma(i); end end bz=0; for i=1:n bz=bz+(diff(f,x(i)*b1(i).2; end bz=sqrt(bz); jz=f; for i=1:n jz=jz+diff(f,x(i)*(j(i)-x(i); co(i)=-diff(f,x(i)*b1(i)/bz; end; final=jz/bz; beita=subs(final,varargin,y); for i=1:n cos1(i)=subs(co(i),varargin,y); end for i=1:n y1(i)=j(i)+beita*b1(i)*cos1(i); end for i=1:n zz=(y1(i)-y(i).2+zz; y(i)=y1(i); end zz=double(zz);enddisp(beita);附件4：分项系数校准分析的程序function fff5() syms r s Z1=r-s; A1(1)=1; A1(2)=1; for n=1:2 B1(n)=input(请输入抗力变异系数=); TT(n)=input(如果此变量服从正态分布请输入0，对数分布输入1,极值型分布输入2 ); if TT(n)=1 a1(n)=log(A1(n)/sqrt(1+B1(n)2); b1(n)=sqrt(log(1+B1(n)2); else if TT(n)=2 aa(n)=pi/B1(n)/A1(n)/sqrt(6); uu(n)=A1(n)-0.5772/aa(n); end end end BB=input(请输入目标可靠度=); A2=A1; for n=1:2 if TT(n)=0 a2(n)=A2(n); b2(n)=B1(n)*A2(n); c2(n)=b2(n)/a2(n); else if TT(n)=1 a2(n)=A2(n)+A2(n)*(-log(A2(n)+a1(n); b2(n)=A2(n)*b1(n); c2(n)=b2(n)/a2(n); else if TT(n)=2 tt(n)=exp(-aa(n)*(A2(n)-uu(n); ff(n)=aa(n)*exp(-aa(n)*(A2(n)-uu(n)-tt(n); FF(n)=exp(-tt(n); CC=subs(FF(n); b2(n)=normpdf(norminv(