架空输电线路设计计算应力和弧垂程序

clc,clear
%% 用向量储存已知数据并约定符号
g=17.8*10^(-6);E=80000;
T1=[10,-20,40,-10,15,15,10,-10,-5,-5,-5,-5,-5]; %按13种气象的顺序构成温度向量T1
%13种气象的顺序依次为：年均气温，最低气温，最高气温，事故气象，外过有风，外过无风，
%内过电压，安装气象，覆冰无风，覆冰有风(强度)，覆冰有风(风偏)，最大风(强度)，最大风(风偏)
R1=[36.5118,36.5118,36.5118,36.5118,36.9323,36.5118,37.6969,36.9323,...
69.6702,70.6750,70.6750,46.9546,42.2139].*10^(-3); %构成比载向量R1,气象顺序与T1相同
T2=[-5,-20,40,-10,15,15,10,-10,-5,-5,10,-5,-5];%按13种气象的顺序构成温度向量T2
%13种气象的顺序依次为：覆冰有风(风偏)，最低气温，最高气温，事故气象，外过有风，外过无风，
%内过电压，安装气象，覆冰无风，覆冰有风(强度)，年均气温，最大风(强度)，最大风(风偏)
R2=[70.6750,36.5118,36.5118,36.5118,36.9323,36.5118,37.6969,36.9323,...
69.6702,70.6750,36.5118,46.9546,42.2139].*10^(-3);%构成比载向量R2,气象顺序与T2相同
A=[83.5935,133.75];

%%以有效临界档距为分界点确定第1状态气象，以其余12种气象为待求气象
for lr=50:50:800
if lr<312.59
a1=A(1);r1=R1(1);t1=T1(1);
for i=2:13
r2=R1(i);t2=T1(i);a(1)=80;iter=0;
for k=1:500
f(k)=a(k)^3-(a1-E*r1^2*lr^2/(24*a1^2)-g*E*(t2-t1))*a(k)^2-E*r2^2*lr^2/24;
df(k)=3*a(k)^2-2*(a1-E*r1^2*lr^2/(24*a1^2)-g*E*(t2-t1))*a(k);
a(k+1)=a(k)-f(k)/df(k);
es=abs(a(k+1)-a(k));
iter=iter+1;
if es<10^(-6)| iter>500
disp('在迭代第iter次得到第i气象档距为lr 时的导线应力a为：')
i,lr,a(k+1)
end
fv=r2*lr^2/(8*a(k+1));
disp('第i气象档距为lr时的导线最大弧垂fv为：')
fv;
fvy=r1*lr^2/(8*a1);
disp('控制气象年均气温的弧垂为： ')
fvy;
FV10(lr/50)=fvy;%年均气温下各档距的弧垂
FV1(lr/50,i-1)=fv;%构造存储弧垂的矩阵M;各气象条件下各档距的弧垂
M1(lr/50,i-1)=a(k+1);%构造存储应力的矩阵M;各气象条件下各档距的弧垂
end
end
elseif lr>312.59;
a1=A(2);r1=R2(1);t1=T2(1)
for i=2:13
r2=R2(i);t2=T2(i);a(1)=60;iter=0;
for k=1:500
f(k)=a(k)^3-(a1-E*r1^2*lr^2/(24*a1^2)-g*E*(t2-t1))*a(k)^2-E*r2^2*lr^2/24;
df(k)=3*a(k)^2-2*(a1-E*r1^2*lr^2/(24*a1^2)-g*E*(t2-t1))*a(k);
a(k+1)=a(k)-f(k)/df(k);
es=abs(a(k+1)-a(k));
iter=iter+1;
if es<10^(-6)
if es<10^(-6) | iter>500%相邻两次迭代结果的残差约等

于0时，终止迭代并输出结果
disp('在迭代第iter次得到第i气象档距为lr 时的导线应力a为：')
i,lr,a(k+1)
end
fv=r2*lr^2/(8*a(k+1));%求解相应气象条件下的最大弧垂
disp('第i气象档距为lr时的导线最大弧垂fv为：')
fv;
fvy=r1*lr^2/(8*a1);
disp('控制气象年均气温的弧垂为： ')
fvy;
FV20(lr/50)=fvy;%最厚覆冰气象条件下各档距的弧垂
FV2(lr/50,i-1)=fv;%各气象条件下各档距的弧垂
M2(lr/50,i-1)=a(k+1);%各气象条件下各档距的应力
end
end
end
end
end

%%求解有效临界档距处的应力和弧垂
lr=312.59;a1=A(1);r1=R1(1);t1=T1(1);%年均气温为第1状态，确定其参数
for i=2:13%求解其余12种气象条件的应力和弧垂
r2=R1(i);t2=T1(i);a(1)=80;iter=0;
for k=1:500
f(k)=a(k)^3-(a1-E*r1^2*lr^2/(24*a1^2)-g*E*(t2-t1))*a(k)^2-E*r2^2*lr^2/24;
df(k)=3*a(k)^2-2*(a1-E*r1^2*lr^2/(24*a1^2)-g*E*(t2-t1))*a(k);
a(k+1)=a(k)-f(k)/df(k);
es=abs(a(k+1)-a(k));
iter=iter+1;
if es<10^(-6)| iter>500
disp('在迭代第iter次得到第i气象档距为lr 时的导线应力a为：')
i,lr,a(k+1)
end
fv=r2*lr^2/(8*a(k+1));
disp('第i气象档距为lr时的导线最大弧垂fv为：')
fv;
fvy=r1*lr^2/(8*a1);
disp('控制气象最厚覆冰的弧垂为： ')
fvy;
FV30(i-1)=fvy;%档距为312.59米时年均气温气象的弧垂
FV3(i-1)=fv;%档距为312.59米时各气象的弧垂
M3(i-1)=a(k+1);%档距为312.59米时气象的应力
end
end