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

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种气象为待求气象forlr=50:50:800iflr312.59a1=A(1);r1=R1(1);t1=T1(1);fori=2:13r2=R1(i);t2=T1(i);a(1)=80;iter=0;fork=1:500f(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;ifes10^(-6)|iter500disp('在迭代第iter次得到第i气象档距为lr时的导线应力a为：')i,lr,a(k+1)endfv=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;各气象条件下各档距的弧垂endendelseiflr312.59;a1=A(2);r1=R2(1);t1=T2(1)fori=2:13r2=R2(i);t2=T2(i);a(1)=60;iter=0;fork=1:500f(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;ifes10^(-6)ifes10^(-6)|iter500%相邻两次迭代结果的残差约等于0时，终止迭代并输出结果disp('在迭代第iter次得到第i气象档距为lr时的导线应力a为：')i,lr,a(k+1)endfv=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);%各气象条件下各档距的应力endendendendend%%求解有效临界档距处的应力和弧垂lr=312.59;a1=A(1);r1=R1(1);t1=T1(1);%年均气温为第1状态，确定其参数fori=2:13%求解其余12种气象条件的应力和弧垂r2=R1(i);t2=T1(i);a(1)=80;iter=0;fork=1:500f(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;ifes10^(-6)|iter500disp('在迭代第iter次得到第i气象档距为lr时的导线应力a为：')i,lr,a(k+1)endfv=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米时气象的应力endend