结构可靠度作业及答案

《结构可靠度理论与应用》1、如图所示圆截面直杆，承受拉力P=120KN，已知材料的强度设计值fy的均值μfy=310MPa，标准差为σfy=25MPa，杆直径d的均值μd=30mm，标准差为σd=3mm，在功能函数为：（1）2(/4)ZdrF；（2）24/ZrFd，在这两种情况下，试用中心点法求其可靠度指标和可靠度。程序：clear;clc;muX=[310;30];sigmaX=[25;3];Z=pi/4*muX(2)^2*muX(1)-120e3;Zx=[pi/4*muX(2)^2;pi/2*muX(1)*muX(2)];betaC1=Z/norm(Zx.*sigmaX)Pr1=normcdf(betaC1)Z=muX(1)-4/pi*120e3/muX(2)^2;Zx=[1;8*120e3/pi/muX(2)^3];betaC2=Z/norm(Zx.*sigmaX)Pr2=normcdf(betaC2)运行结果：betaC1=2.0977Pr1=0.9820betaC2=3.3259Pr2=0.99962、粒状土承受剪切应力τ=52KPa，其剪切面法向应力w服从正态分布，均值为100KPa，标准差为20KPa，土的磨擦角φ服从正态分布，均值为35º，标准差为5º(=0.0873弧度)。w和φ相互独立，极限状态方程为：Z=wtanφ-τ=0，用中心点法计算β值和失效概率pf。程序：clear;clc;muX=[100;35*pi/180];sigmaX=[20;5*pi/180];Z=muX(1)*tan(muX(2))-52;Zx=[tan(muX(2));muX(1)/cos(muX(2))^2];betaC=Z/norm(Zx.*sigmaX)Pf=normcdf(-betaC)运行结果：betaC=0.9429Pf=0.17293、某钢梁承受确定性弯矩138MkNm，抗弯截面模量63(89010,0.05)，服从正态分布；钢材强度f服从对数正态分布（262,fMPa0.1f），极限状态方程为ZfWM=0。试用中心点法和验算点法求可靠指标及梁的失效概率fP，并比较其计算结果。中心点法：clear;clc;muX=[262e6;8.9e-4];cvX=[0.1;0.05];sigmaX=cvX.*muX;Z=muX(1)*muX(2)-1.38e5;Zx=[muX(2);muX(1)];beta1=Z/norm(Zx.*sigmaX)Pf1=normcdf(-beta1)结果：beta1=3.6509Pf1=1.3066e-004验算点法：clear;clc;muX=[262e6;890e-6];cvX=[0.1;0.05];sigmaX=cvX.*muX;sLn=sqrt(log(1+sigmaX(1)/muX(1)^2));mLn=log(muX(1))-sLn*2/2;muX1=muX;sigmaX1=sigmaX;x=muX;normX=eps;whileabs(norm(x)-normX)/normX1e-6normX=norm(x);Z=x(1)*x(2)-1.38e5;Zx=[x(2);x(1)];cdfX=logncdf(x(1),mLn,sLn);pdfX=lognpdf(x(1),mLn,sLn);nc=norminv(cdfX);sigmaX1(1)=normpdf(nc)/pdfX;muX1(1)=x(1)-nc*sigmaX1(1);Zs=Zx.*sigmaX;aX=-Zs/norm(Zs);beta=(Z+Zx'*(muX1-x))/norm(Zs)x=muX1+beta*sigmaX1.*aX;endbeta2=betaPf2=normcdf(-beta2)结果：Beta2=3.9421Pf2=4.0391e-005验算点法得到的可靠度指标大一些，验算点法考虑了随机变量的分布，使用验算点法得到的结果更加准确。4、已知某钢筋混凝土受压短柱的极限状态方程为(,,)0ZgRGQRGQ，抗力R服从对数正态分布0.17R；恒载(GGN53,3.71)GkNkN，服从正态分布；活载Q服从极值I型分布，70,QkN20.31QkN。试用JC法求当目标可靠指标[]=3.7时，构件截面的抗力平均值R？程序：clear;clc;deltaR=0.17;muG=53;sigmaG=3.71;muQ=70;sigmaQ=20.31;beta=3.7;sigmaLnR=sqrt(log(1+deltaR^2));aEv=sqrt(6)*sigmaQ/pi;uEv=-psi(1)*aEv-muQ;S0=[muG,muQ];R0=muG+muQ;R1=0;cosR=0;whileabs(R1-R0)1e-6R1=R0;cdfQ=1-evcdf(-S0(2),uEv,aEv);pdfQ=evpdf(-S0(2),uEv,aEv);sigmaQ1=normpdf(norminv(cdfQ))/pdfQ;muQ1=S0(2)-norminv(cdfQ)*sigmaQ1;muS=[muG,muQ1];sigmaS=[sigmaG,sigmaQ1];sigmaR1=sigmaLnR*R0;cosS=-([-1,-1].*sigmaS)./norm([1,-1,-1].*[sigmaR1,sigmaS]);cosR=-1*sigmaR1/norm([1,-1,-1].*[sigmaR1,sigmaS]);S0=muS+cosS.*sigmaS*beta;R0=S0(1)+S0(2);endRR=R0*sqrt(1+deltaR^2)*exp(-beta*sigmaLnR*cosR);muR=RR结果：muR=320.01195、设某构件正截面强度计算的极限状态方程为Z=R-S=0。其中R和S分别为正态和极值I型分布的随机变量，其统计量为R(100,20)和S(80,24)，20和24为标准差。试用JC法和蒙特卡罗模拟分别求解构件失效概率。JC法：clear;clc;muX=[100;80];sigmaX=[20;24];g=muX(1)-muX(2);gX=[muX(2);muX(1)];aEv=sqrt(6)*sigmaX(2)/pi;uEv=-psi(1)*aEv-muX(2);muX1=muX;sigmaX1=sigmaX;x=muX;normX=eps;whileabs(norm(x)-normX)/normX1e-6normX=norm(x);g=x(1)-x(2);gX=[1;-1];cdfX=1-evcdf(-x(2),uEv,aEv);pdfX=evpdf(-x(2),uEv,aEv);nc=norminv(cdfX);sigmaX1(2)=normpdf(nc)/pdfX;muX1(2)=x(2)-nc*sigmaX1(2);gs=gX.*sigmaX1;alphaX=-gs/norm(gs);beta1=(g+gX'*(muX1-x))/norm(gs);x=muX1+beta1*sigmaX1.*alphaX;endPf1=normcdf(-beta1)结果：Pf1=0.2221蒙特卡罗法：clear;clc;muX=[100;80];sigmaX=[20;24];aEv=sqrt(6)*sigmaX(2)/pi;uEv=-psi(1)*aEv-muX(2);g=muX(1)-muX(2);gX=[muX(2);muX(1)];nS=1e6;ig=ones(nS,1);x=[normrnd(muX(1),sigmaX(1),nS,1),-evrnd(uEv,aEv,nS,1)];g=x(:,1)-x(:,2);nF=sum(ig(g0));Pf2=nF/nS结果：Pf2=0.23876、设构件的极限状态方程为：21213410Zxxxxx。式中，111(,)(25,0.23)xxx，服从对数正态分布；222(,)(0.0113,0.3)xxx，为正态分布；333(,)(0.0006,0.3)xxx，为正态分布；444(,)(0,0.1)xxx，为正态分布。试用蒙特卡洛法计算该结构构件的可靠度。程序：clear;clc;muX=[25;0.0113;0.0006;0];sigmaX=[5.75;0.3;0.3;0.1];sLn=sqrt(log(1+sigmaX(1)/muX(1)^2));mLn=log(muX(1))-sLn*2/2;nS=1e7;ig=ones(nS,1);x1=lognrnd(mLn,sLn,1,nS);x2=normrnd(muX(2),sigmaX(2),1,nS);x3=normrnd(muX(3),sigmaX(3),1,nS);x4=normrnd(muX(4),sigmaX(4),1,nS);g=1-x1.*x2-x1.*x1.*x3-x4;nR=sum(ig(g=0));Pr=nR/nS结果：Pr=0.50127、设构件的极限状态方程为1234Zxxxx，111(,)(2234.32,0.1)xxx，为对数正态分布；222(,)(949.59,0.1)xxx，为对数正态分布；333(,)xxx(1521.9,0.109)，为正态分布；444(,)(496.1,0.292)xxx，为极值I型分布。试用蒙特卡洛法计算该结构构件的可靠度。程序：clear;clc;muX=[2234.32;949.59;1521.9;496.1];sigmaX=[0.1;0.1;0.109;0.292];sLn1=sqrt(log(1+sigmaX(1)/muX(1)^2));mLn1=log(muX(1))-sLn1*2/2;sLn2=sqrt(log(1+sigmaX(2)/muX(2)^2));mLn2=log(muX(2))-sLn2*2/2;nS=1e6;ig=ones(nS,1);x1=lognrnd(mLn1,sLn1,1,nS);x2=lognrnd(mLn2,sLn2,1,nS);x3=normrnd(muX(3),sigmaX(3),1,nS);symsxalphak;EVIpdf='alpha*exp(-alpha*(x-k)-exp(-alpha*(x-k)))';EVIcdf='exp(-exp(-alpha*(x-k)))';alpha=1.2825/sigmaX(4);k=muX(4)-0.5772/alpha;EVIpdfStar=eval(vpa(subs(EVIpdf,x,muX(4))));EVIcdfStar=eval(vpa(subs(EVIcdf,x,muX(4))));aEv=normpdf(norminv(EVIcdfStar))/EVIpdfStar;uEv=muX(4)-norminv(EVIcdfStar)*sigmaX(4);x4=evrnd(uEv,aEv,1,nS);g=x1+x2-x3-x4;nR=sum(ig(g=0));Pr=nR/nS结果：Pr=1