结构可靠度分析与设计的编程实践

Project #2 of Loads and Structural Design Methods HIT, Lu Dagang, Spring, 2012Undergraduate Course Loads & Structural Design MethodsProject #2结构可靠度分析与设计的编程实践姓名:杨冬冬 学号： 班级：摘要：本文使用MATLAB软件，运用设计点法、RF法、蒙特卡洛方法等几种求可靠度指标的方法对相关习题进行了编程计算，得到了与预期结果十分接近的结果。1. 引言科学实验和力学的发展是土木工程由完全凭直觉和经验走向了科学与经验相结合的道路。土木工程的设计理论和方法的发展经历了几个阶段：容许应力设计法到极限强度设计法到极限状态设计法。极限状态设计法引入可靠性的概念，应用概率论和数理统计的知识进行设计，经历了四个水准：半经验半概率、近似概率、全概率、最优概率，使结构的设计更加科学和可控。今天，国际上结构设计理论的发展趋势是概率极限状态设计法。在结构建造和使用过程中，结构可靠与不可靠是不可预知的，这是因为建造和使用中存在了诸多不确定性。而结构的可靠度关乎人们的生命财产安全，至关重要，因此学习计算结构的可靠度很有必要。2. 结构可靠度分析的基本原理结构可靠度指的是结构在规定的时间，在规定的条件下，完成预定功能的概率，而通过极限状态的概率计算可以得到安全概率和失效概率，而且。而结构可靠度以结构可靠指标计量，定义为。通过结构可靠指标对结构的可靠度进行判别有着较为科学的依据，可靠度指标与失效概率相关，表明其具有概率意义。在取值中，可靠指标取值越大，则结构的失效概率越小；结构的安全概率越大，结构的保证率越大，更加安全。当用概率描述结构的可靠性时，结构可靠度是一个概率计算问题，即计算结构可靠度就是计算结构在规定时间内、规定条件下结构能够完成预定功能的概率。从简单到复杂或精确程度的不同，先后提出的可靠度计算方法有一次二矩法，二次二矩法，蒙特卡洛方法及其他方法。一次二矩法又分为中心点法和设计点法，其中验算点法是目前可靠度分析最常用的方法。对于一般形式的结构功能函数。式中为n个相互独立随机变量，其平均值和标准差是已知的。在计算=/时，需要计算功能函数Z的平均值和标准差，如果进一步已知各随机变量概率密度函数，精确计算这两个值要进行高维数值积分，这又回到了与直接根据随机变量概率分布计算可靠度相同的问题；如果各随机变量的概率密度函数是未知的，则不能给出和的近似结果，然而展开点的选择却是一个值得考虑的问题，一个简单的方法是将展开点选为平均值处，即中心点。功能函数Z的平均值和方差近似为从而，按中心点法计算的可靠指标为由上式可以看出，中心点法使用了结构功能函数的一次泰勒级数展开式和随机变量的前两阶矩（平均值和方差），故称为一次二阶矩方法。但是该方法存在很多缺点：1.功能函数在平均值处展开不尽合理；2.对于力学意义相同但数学表达式不同的结构功能函数，由此计算的可能不同；3.没有考虑随机变量的概率分布。 于是产生了改进的一次二阶矩方法也就是设计点法：在极限状态曲面上寻找验算点，并在此基础上进行泰勒级数展开，应用随机变量的前二阶矩，采用非正态随机变量的当量正态化，迭代求解结构的失效概率的一种方法，该方法简称验算点法。RF法：RF法用于处理随机变量非正态分布的情况，为JCSS（Joint Committee of Structural Safety）推荐的方法，所计算的可靠指标又称RackwitzFiessler可靠指标。JC法的基本原则：按照等效正态化原则将非正态随机变量转化为当量正态化随机变量解决由于非正态随机变量导致的可靠指标与失效概率不一一对应的不足转化之后可以利用验算点法求解蒙特卡洛方法：对功能函数中所有随机变量进行随机抽样，根据随机变量的抽样值进行结构功能函数值Z=g(x)的计算，根据结构可靠度基本理论可知当功能函数Z0时，结构可靠，当功能函数Z=0时，结构处于极限状态，当功能函数Z0.0001 n=n+1; beta(n)=-12895/(3000*a1(n-1)+4000*a2(n-1)+750*beta(n-1)*a1(n-1)*a2(n-1)-124.2*beta(n-1)*a3(n-1); a1(n)=-(3000+750*beta(n)*a2(n-1)/sqrt(3000+750*beta(n)*a2(n-1)2+(4000+750*beta(n)*a1(n-1)2+(-124.2)2); a2(n)=-(4000+750*beta(n)*a1(n-1)/sqrt(3000+750*beta(n)*a2(n-1)2+(4000+750*beta(n)*a1(n-1)2+(-124.2)2); a3(n)=124.2/sqrt(3000+750*beta(n)*a2(n-1)2+(4000+750*beta(n)*a1(n-1)2+(-124.2)2);endbetaa1a2a3作业5.3clc,clearbeta(1)=2;beta(2)=3;n=2;a1(2)=0;a2(2)=0;cov=1,0,0;0,1,0;0,0,1;while abs(beta(n)-beta(n-1)0.0001 n=n+1; fr(n)=4375+a1(n-1)*beta(n-1)*1312.5; s(n)=121.2; g=0,s(n)*1312.5,-39*1137.5; a1(n)=-g*cov(1,2,3,1)/sqrt(g*cov*g); a2(n)=-g*cov(1,2,3,2)/sqrt(g*cov*g); a3(n)=-g*cov(1,2,3,3)/sqrt(g*cov*g); beta(n)=(39000-4375*121.2)/(a2(n)*1312.5*121.2-39*a3(n)*1137.5);enda1a2beta附件2：一次可靠度方法（RF法）的程序例题5.9%JC%g=R-Qmu_Q=100;mu_R=200;sigma_Q=12;sigma_R=20;Rstar1=mu_R;Qstar1=mu_Q;Qstar=0;Rstar=0; while abs(Qstar,Rstar-Qstar1,Rstar1)=0.0001,0.0001 Qstar=Qstar1;Rstar=Rstar1;%R服从对数正态分布mu_lnR=log(mu_R/sqrt(1+(sigma_R/mu_R)2); mu_R1=Rstar*(1-log(Rstar)+mu_lnR);sigma_lnR=sqrt(log(1+(sigma_R/mu_R)2); sigma_R1=Rstar*sigma_lnR; %Q服从极值型分布alpha=pi/(sqrt(6)*sigma_Q); u=mu_Q-0.5772/alpha;F_Q=exp(-exp(-alpha*(Qstar-u);f_Q=alpha*exp(-alpha*(Qstar-u)-exp(-alpha*(Qstar-u);sigma_Q1=normpdf(norminv(F_Q)/f_Q;mu_Q1=Qstar-norminv(F_Q)*sigma_Q1; %正态化后极限状态方程变为g=R1-Q1beta=(Rstar-Qstar+(mu_R1-Rstar)-(mu_Q1-Qstar)/sqrt(sigma_R12+sigma_Q12); alpha_R1=-sigma_R1/sqrt(sigma_R12+sigma_Q12);alpha_Q1=sigma_Q1/sqrt(sigma_R12+sigma_Q12); Rstar1=mu_R1+alpha_R1*sigma_R1*beta;Qstar1=mu_Q1+alpha_Q1*sigma_Q1*beta;betaRstarQstarend例题5.10（3）例题5.11%g=Z*Fy-M% Z normal; Fy lognormal; M extrem type Iclc,clearmu_Z=100; deta_Z=0.04;%;mu_Fy=40;deta_Fy=0.1;mu_M=2000;deta_M=0.1;sigma_Z=mu_Z*deta_Z;sigma_Fy=mu_Fy*deta_Fy;sigma_M=mu_M*deta_M;sigma_Z1=sigma_Z;sigma_Fy1=sigma_Fy;sigma_M1=sigma_M; beta=0; alpha_Z1=0; alpha_Fy1=0; alpha_M1=0;Zstar=mu_Z+alpha_Z1*beta*sigma_Z1;Fystar=mu_Fy+alpha_Fy1*beta*sigma_Fy1;Mstar=mu_M+alpha_M1*beta*sigma_M1; Zstar1=Zstar;Fystar1=Fystar;Mstar1=Mstar;Zstar=0;Fystar=0;Mstar=0; while abs(Zstar Fystar Mstar-Zstar1 Fystar1 Mstar1)=0.001 0.001 0.001 Zstar=Zstar1;Fystar=Fystar1;Mstar=Mstar1; %Z服从正态分布 mu_Z1=mu_Z;sigma_Z1=sigma_Z; %Fy服从对数正态分布 mu_lnFy=log(mu_Fy/sqrt(1+(sigma_Fy/mu_Fy)2); mu_Fy1=Fystar*(1-log(Fystar)+mu_lnFy); sigma_lnFy=sqrt(log(1+(sigma_Fy/mu_Fy)2); sigma_Fy1=Fystar*sigma_lnFy; %M服从极值型分布 alpha=pi/(sqrt(6)*sigma_M); u=mu_M-0.5772/alpha; F_M=exp(-exp(-alpha*(Mstar-u);f_M=alpha*exp(-alpha*(Mstar-u)-exp(-alpha*(Mstar-u); sigma_M1=normpdf(norminv(F_M)/f_M; mu_M1=Mstar-norminv(F_M)*sigma_M1; sqrt1=sqrt(Fystar1*sigma_Z1)2+(Zstar1*sigma_Fy1)2+sigma_M12); alpha_Z1=-sigma_Z1*Fystar/sqrt1; alpha_Fy1=-sigma_Fy1*Zstar/sqrt1; alpha_M1=sigma_M1/sqrt1; b=solve(mu_Z1+alpha_Z1*x*sigma_Z1)*(mu_Fy1+alpha_Fy1*x*sigma_Fy1)-(mu_M1+alpha_M1*x*sigma_M1)=0); c=eval(b); beta=c(1); Zstar1=mu_Z1+alpha_Z1*beta*sigma_Z1; Fystar1=mu_Fy1+alpha_Fy1*beta*sigma_Fy1; Mstar1=mu_M1+alpha_M1*beta*sigma_M1; beta Zstar Fystar Mstarend作业5.4clc,clearmu_Y=12.5;deta_Y=0.15;mu_X=24;deta_X=0.125;sigma_Y=mu_Y*deta_Y;sigma_X=mu_X*deta_X;sigma_Y1=sigma_Y;sigma_X1=sigma_X; beta=0; alpha_Y1=0; alpha_X1=0;Ystar=mu_Y+alpha_Y1*beta*sigma_Y1;Xstar=mu_X+alpha_X1*beta*sigma_X1; Ystar1=Ystar;Xstar1=Xstar;Ystar=0;Xstar=0; while abs(Ystar Xstar-Ystar1 Xstar1)=0.001 0.001 Ystar=Ystar1;Xstar=Xstar1; mu_lnY=log(mu_Y/sqrt(1+(sigma_Y/mu_Y)2); mu_Y1=Ystar*(1-log(Ystar)+mu_lnY); sigma_lnY=sqrt(log(1+(sigma_Y/mu_Y)2); sigma_Y1=Ystar*sigma_lnY; alpha=pi/(sqrt(6)*sigma_X); u=mu_X-0.5772/alpha; F_M=exp(-exp(-alpha*(Xstar-u);f_M=alpha*exp(-alpha*(Xstar-u)-exp(-alpha*(Xstar-u); sigma_X1=normpdf(norminv(F_M)/f_M; mu_M1=Xstar-norminv(F_M)*sigma_X1; sqrt1=sqrt(3*sigma_Y1)2+sigma_X12); alpha_Y1=-sigma_Y1*Ystar/sqrt1; alpha_X1=sigma_X1/sqrt1; b=solve(3*(mu_Y1+alpha_Y1*x*sigma_Y1)-(mu_M1+alpha_X1*x*sigma_X1)=0); c=eval(b); beta=c(1); Ystar1=mu_Y1+alpha_Y1*beta*sigma_Y1; Xstar1=mu_M1+alpha_X1*beta*sigma_X1; beta Ystar Xstarend 附件3：MCS法的计算程序例题5.11clc,clearN = 1e6;%FY服从对数正态分布miu_FY = 40;V_FY = 0.1;sigma_FY = miu_FY*V_FY;miu_lnFY = log(miu_FY/sqrt(1 + V_FY2);sigma_lnFY = sqrt(log(1 + V_FY2);FY = lognrnd(miu_lnFY,sigma_lnFY,N,1);subplot(3,1,1)hist(FY,20)xlabel(FY)ylabel(Frequency (%)title(Histogram of FY)%Z服从正态分布miu_Z = 100;V_Z = 0.04;sigma_Z = miu_Z*V_Z;Z = normrnd(miu_Z,sigma_Z,N,1);subplot(3,1,2)hist(Z,20)xlabel(Z)ylabel(Frequency (%)title(Histogram of Z)%M服从极值I型分布miu_M = 2000;V_M = 0.1;sigma_M = miu_M*V_M;alpha = 1.282/sigma_M;u = miu_M - 0.45*sigma_M;r = rand(N,1);M = u - log(-log(r)/alpha;subplot(3,1,3)hist(M,20)xlabel(M)ylabel(Frequency (%)title(Histogram of M)g =Z.*FY-M; Nf = zeros(N,1);for i=1:N Nf(i) = any(g(i) 0.0);end n = sum(Nf) %计算失效概率Pf = n/N %计算可靠度指标beta = -norminv(Pf)附件4：分项系数校准分析的程序例题6.3clear,clc%g=r-sg-slb=3.2;dr=0.17;usg=50;ssg=3.5;usl0=70;ssl0=20; sg=usg;sl=usl0;r=sg+sl; sr1=sqrt(log(1+dr2); a=3.14/(ssl0*sqrt(6);u=usl0-0.5772/a;t=exp(-a*(sl-u);f=a*exp(-a*(sl-u)-t);y=exp(-t);ssl=normpdf(norminv(y,0,1),0,1)/f;usl=sl-ssl*norminv(y,0,1); for i=1:100 r0=r; sr=r*sr1; a=3.14/(ssl0*sqrt(6); u=usl0-0.5772/a; t=exp(-a*(sl-u); f=a*exp(-a*(sl-u)-t); y=exp(-t); ssl=normpdf(norminv(y,0,1),0,1)/f; usl=sl-ssl*norminv(y,0,1); ar=-sr/sqrt(sr2+ssg2+ssl2); asg=ssg/sqrt(sr2+ssg2+ssl