结构可靠度作业及答案.doc

(完整word版)结构可靠度作业及答案
《结构可靠度理论与应用》
1、如图所示圆截面直杆，承受拉力P=120KN，已知材料的强度设计值fy的均值μfy=310MPa，标准差为σfy=25MPa，杆直径d的均值μd=30mm，标准差为σd=3mm，在功能函数为：（1） ；（2） ，在这两种情况下，试用中心点法求其可靠度指标和可靠度。
程序：
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.0977
Pr1 =0.9820
betaC2 =3.3259
Pr2 = 0.9996

2、粒状土承受剪切应力τ=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.9429
Pf =0.1729

3、某钢梁承受确定性弯矩，抗弯截面模量，服从正态分布；钢材强度f服从对数正态分布（），极限状态方程为=0。试用中心点法和验算点法求可靠指标及梁的失效概率，并比较其计算结果。
中心点法：
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.6509
Pf1 =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;
while abs(norm(x)-normX)/normX>1e-6
normX=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;
end
beta2=beta
Pf2=normcdf(-beta2)
结果：
Beta2 =3.9421
Pf2 = 4.0391e-005
验算点法得到的可靠度指标大一些，验算点法考虑了随机变量的分布，使用验算点法得到的结果更加准确。

4、已知某钢筋混凝土受压短柱的极限状态方程为，
抗力R服从对数正态分布；恒载，服从正态分布；活载服从极值I型分布，。试用JC法求当目标可靠指标[]=3.7时，构件截面的抗力平均值
程序：
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;
while abs(R1-R0)>1e-6
R1=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);
end
RR=R0*sqrt(1+deltaR^2)*exp(-beta*sigmaLnR*cosR);
muR=RR
结果：
muR=320.0119

5、设某构件正截面强度计算的极限状态方程为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;
while abs(norm(x)-normX)/normX>1e-6
normX=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;
end
Pf1=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(g<0));
Pf2=nF/nS
结果：
Pf2=0.2387

6、设构件的极限状态方程为：。式中，，服从对数正态分布；，为正态分布；，为正态分布；，为正态分布。试用蒙特卡洛法计算该结构构件的可靠度。
程序：
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.5012

7、设构件的极限状态方程为，，为对数正态分布；，为对数正态分布；，为正态分布；，为极值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);
syms x alpha k;
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