%输入 b:坡度 h:坡高 gama:重度 phi0:内摩擦角 c:粘聚力
%输出 Fsmin x0 y0 R0
Fs=2.0;
b=0.5;
h=25;
gama=20;
phi0=19.6;
c=3;
phi0=phi0*pi/180;
alpha=atan(b);
L=h/sin(alpha);%坡长
m=L*cos(alpha);
beta=[500,500];
x=[];
y=[];
xd=[];
yd=[];
qs=[];
T=[];
P=[];
l=[];
P(1)=0;
y2=[];
h1=[];
W1=[];
for x=-2*h:0.1:h %圆心横坐标取值范围
for y=h:0.1:3*h %圆心纵坐标取值范围
R=sqrt(x^2+y^2); %新圆半径
x1=sqrt(R^2-(h-y)^2)+abs(x); %x最大取值
d1=x1/500; %每块土条宽度
for xd=x1:-d1:0 %圆弧曲线的x坐标取值
yd=y-sqrt(R^2-(x-xd)^2); %圆弧所在方程
beta=atan((xd-x)/(y-yd));%土条坡底角度
if xd<=m
y2=tan(alpha)*xd; %斜坡方程
h1=abs(y2-yd); %土条高度
W1=gama*h1*d1; %土条的重力
else
h1=abs(h-yd);
W1=gama*h1*d1; %土条的重力
end
end
s=1;
n=500;
for i=s:n
if i==1
l(i)=d1/cos(beta(i)); %斜底面长度
qs(i)=0;
T(i)=(c*l(i)+W1(i)*cos(beta(i))*tan(phi0)/Fs);%底面剪力
P(i)=W1(i)*sin(beta(i))-T(i); %条间推力
i=i+1;
else
l(i)=d1/cos(beta(i)); %斜底面长度
qs(i)=cos(beta(i-1)-beta(i))-tan(phi0)*sin(beta(i-1)-beta(i))/Fs; %传递系数
T(i)=(c*l(i)+W1(i)*cos(beta(i))*tan(phi0)/Fs);%底面剪力
P(i)=W1(i)*sin(beta(i))+qs(i)*P(i-1)-T(i); %条间推力
end
i=i+1;
end
P;
if P(n)>0.5
Fs=Fs-0.0001;
elseif P(n)<-0.5
Fs=Fs+0.0001;
else
disp(Fs);
end
end
end
以上是我的代码,报错显示
索引超出数组范围。
出错 nnn6666 (line 53)
l(i)=d1/cos(beta(i)); %斜底面长度