sectioninput_use.m

薄壁结构有限元计算

% 截面基本参数的输入与基本常数的计算

% 截面端点数
Nv=48;

% 截面边数
Ne=54;

% 截面关联矩阵
B=zeros(Nv,Ne);
Be=[1	3	1	4	5	6	7	5	7	9	9	11	13	11	14	15	16	17	15	16	19	20	19	20	21	22	23	24	22	23	25	27	28	28	30	30	31	32	33	34	34	35	36	38	32	39	40	40	42	43	43	44	45	47
2	4	6	7	6	7	8	15	9	10	14	12	14	17	18	16	17	18	19	20	20	21	22	23	25	23	24	25	30	26	27	28	29	31	31	32	33	33	34	35	36	39	37	39	42	40	41	44	43	44	45	48	46	48
];
for i=1:Ne
    B(Be(1,i),i)=1;
    B(Be(2,i),i)=-1;
end

% 截面始端关联矩阵、终端关联矩阵
Ba=zeros(Nv,Ne);
Bm=zeros(Nv,Ne);
for i=1:Nv
    for j=1:Ne
        if B(i,j)>0
            Ba(i,j)=B(i,j);
        elseif B(i,j)<0
            Bm(i,j)=B(i,j);
        end
    end
end

% 截面基本回路矩阵
Cf=zeros(Ne+1-Nv,Ne);
Cf1=[1	1	1	1	1	2	2	3	3	4	4	5	5	5	5	5	5	6	6	7	7	7	7	7	7
5	6	9	11	15	16	20	21	24	22	25	26	27	28	31	32	34	35	37	38	39	40	42	46	48
];
Cf2=[1	1	1	1	2	2	3	3	3	4	4	5	5	6	6	7	7	7
8	16	17	18	19	21	23	26	27	24	28	29	35	36	38	45	49	50
];
for i=1:size(Cf1,2)
    Cf(Cf1(1,i),Cf1(2,i))=1;
end
for i=1:size(Cf2,2)
    Cf(Cf2(1,i),Cf2(2,i))=-1;
end

% 截面端点坐标矢量与厚度矩阵
X=[7.75	7.75	7.75	7.75	9.6	9.6	9.6	9.6	14.6	14.6	17.32	17.32	17.32	17.32	19.2	19.2	19.2	19.2	23.85	23.85	23.85	27.61	27.61	27.61	27.61	29.95	33.65	35.1	35.1	36.15	36.15	39.3	39.3	39.3	39.3	41.17	41.17	41.17	41.17	47.15	47.15	48.9	48.9	48.9	50.77	50.77	50.77	50.77
];X=X';
Y=[7.15	6.1	1.6	0.5	15.5	7.15	0.65	-1	0.65	-1	7.15	6.1	1.6	0.5	15.5	12	7.15	0.5	15.5	12	7.55	15.5	12.7	12	7.55	12.7	8.25	8.25	7	15.5	8.25	15.5	8.25	7.15	0.5	7.15	6.1	1.6	0.5	0.5	-0.81	15.5	7.15	0.5	7.15	6.1	1.6	0.5
];Y=Y';
T=[1.25	1.25	1	1	1	1	1	1	0.7	1	1	1.25	1.25	1	1	1.5	1.5	1	1	1	0.5	1.85	1	0.4	1.1	0.5	1.1	1.1	1	0.25	1.5	 1	0.5	1	0.5	1	1	1	1	1	1	1	1.25	1.25	1	1	2.5	 1	1	1	1	1	1.25	1.25
];T=diag(T);
% 截面边的长度
L=[1.05	0.95	1.88	1.88	8.35	6.5	1.65	9.6	5	1.5	2.72	1.05	1.1	1.88	1.88	3.5	4.85	6.65	4.65	4.65	3.5	4.45	3.76	3.76	3.76	2.8	0.7	4.45	8.54	2.34	6.04	1.45	1.25	1.05	7.25	3.15	3.15	7.25	1.1	6.65	1.88	1.88	1.05	1.1	9.6	5.97	1.31	1.75	8.35	6.65	1.87	1.87	1.05	1.1
];L=diag(L);

% 截面的边到扭转中心的距离
H=[7.75	7.75	7.15	0.5	9.6	9.6	9.6	15.5	0.5	14.6	0.5	17.32	17.32	7.15	0.5	19.2	19.2	19.2	15.5	12	23.85	23.85	15.5	12	7.55	27.61	27.61	27.61	15.5	12.7	8.25	8.25	35.1	8.25	36.15	15.5	8.25	39.3	39.3	39.3	7.15	0.5	41.17	41.17	15.5	0.5	47.15	0.5	48.9	48.9	7.15	0.5	50.77	50.77
];H=H';

% 横竖边的判别
direction=[-1	-1	1	1	-1	-1	-1	1	1	-1	1	-1	-1	1	1	-1	-1	-1	1	1	-1	-1	1	1	1	-1	-1	-1	1	1	1	1	-1	1	-1	1	1	-1	-1	-1	1	1	-1	-1	1	1	-1	1	-1	-1	1	1	-1	-1
];

% 截面回路围绕面积
Ac=[142.56	16.275	13.16	16.732	61.915	22.8375	144
];Ac=Ac';

% 截面基本关联矩阵
Bf=B(2:Nv,:);

% 计算截面的扇性坐标
I1=ones(Ne,1);
Zero=zeros(Ne,1);
Kw=[(Bf*T*inv(L)*B')' (I1'*L*T*(Ba-Bm)')']';
Qw=-[(Bf*T*H)' 0]';
FC=Kw\Qw;        %"FC" means "flabellate coordinate"

% 计算始端扇性静矩值
Bfm=Bm(2:Nv,:);
omega1=L*T*(Ba-Bm)'*FC/2;
omega2=L*L*(2*Ba-Bm)'*FC/6;
Ks=[Bf;Cf*L*inv(T)];
qs=-[Bfm*L*T*(Ba-Bm)'/2;Cf*L*L*(2*Ba-Bm)'/6];
Sw1=Ks\(qs*FC);

% 计算终端扇性静矩值
Sw2=Sw1+L*T*(Ba-Bm)'*FC/2;
Sw3=[Sw1 Sw2];

% 计算截面静矩值
Sx1=Ks\(qs*X);
Sy1=Ks\(qs*Y);
Sx2=Sx1+L*T*(Ba-Bm)'*X/2;
Sy2=Sy1+L*T*(Ba-Bm)'*Y/2;
Sx3=[Sx1 Sx2];
Sy3=[Sy1 Sy2];
% 计算转动惯量、面积、扭心位置、弯心位置
shareborders=zeros(Ne-Nv+1,Ne-Nv+1);corc=zeros(Ne,2);
for i=1:Ne
      if sum(Cf(:,i))>0
          cfp(:,i)=zeros(Ne-Nv+1,1);
      elseif sum(Cf(:,i))<0
          cfp(:,i)=zeros(Ne-Nv+1,1);
      else
          cfp(:,i)=Cf(:,i);
      end
end
for i=1:Ne
       for j=1:(Ne-Nv+1)
           if cfp(j,i)>0
              corc(i,1)=j;
           elseif cfp(j,i)<0
              corc(i,2)=j;
           else
           end
       end
end
for i=1:Ne
    if corc(i,1)>0
        shareborders(corc(i,1),corc(i,2))=L(i,i)/T(i,i);
       shareborders(corc(i,2),corc(i,1))=L(i,i)/T(i,i);
    else
    end
end
cor=zeros(Ne,2);
for i=1:Ne
    for j=1:Nv
        if B(j,i)>0
            cor(i,1)=j;
        elseif B(j,i)<0
            cor(i,2)=j;
        end
    end    
end     
Ix=0;Iy=0;Ixy=0;Ax=0;Ay=0;Sx=0;Sy=0;Ip=0;Ixf=0;Iyf=0;If=0;Ssx=zeros(Ne,1);Ssy=zeros(Ne,1);Sssx=zeros(Ne-Nv+1,Ne);Sssy=zeros(Ne-Nv+1,Ne);Ssssx=zeros(Ne-Nv+1,1);Ssssy=zeros(Ne-Nv+1,1);St=zeros(Ne-Nv+1,1);
for i=1:Ne
    Ip=Ip+(H(i)^2)*L(i,i)*T(i,i);
    If=If+(FC(cor(i,1))*(2*FC(cor(i,1))+FC(cor(i,2)))+FC(cor(i,2))*(2*FC(cor(i,2))+FC(cor(i,1))))*L(i,i)*T(i,i)/6;
    if direction(i)>0
        Iy=Iy+abs((X(cor(i,1))^3-X(cor(i,2))^3))*T(i,i)/3;
        Ix=Ix+L(i,i)*T(i,i)*(Y(cor(i,1))^2);
        Ixy=Ixy+(X(cor(i,1))+X(cor(i,2)))*L(i,i)*Y(cor(i,1))*T(i,i)/2;
        Ax=Ax+L(i,i)*T(i,i);
        Sy=Sy+(X(cor(i,1))+X(cor(i,2)))*L(i,i)*T(i,i)/2;
        Sx=Sx+Y(cor(i,1))*L(i,i)*T(i,i);
        Iyf=Iyf+(FC(cor(i,1))+FC(cor(i,2)))*L(i,i)*T(i,i)*Y(cor(i,1))/2;
        Ixf=Ixf+(FC(cor(i,1))*(2*X(cor(i,1))+X(cor(i,2)))+FC(cor(i,2))*(2*X(cor(i,2))+X(cor(i,1))))*L(i,i)*T(i,i)/6;
        Ssx(i)=(Sx1(i)+Sx2(i))*L(i,i)/T(i,i)/2;
        if sum(Cf(:,i))>0
            Ssy(i)=(Sy1(i)+Sy2(i))*L(i,i)/T(i,i)/2-L(i,i)^3/12;
        elseif sum(Cf(:,i))<0
            Ssy(i)=(Sy1(i)+Sy2(i))*L(i,i)/T(i,i)/2+L(i,i)^3/12;
        else
            for j=1:(Ne-Nv+1)
                if Cf(j,i)<0
                    Sssy(j,i)=(Sy1(i)+Sy2(i))*L(i,i)/T(i,i)/2+L(i,i)^3/12;
                elseif Cf(j,i)>0
                    Sssy(j,i)=-(Sy1(i)+Sy2(i))*L(i,i)/T(i,i)/2-L(i,i)^3/12;
                end
            end
        end
    else
        Ix=Ix+abs((Y(cor(i,1))^3-Y(cor(i,2))^3))*T(i,i)/3;
        Iy=Iy+L(i,i)*T(i,i)*(X(cor(i,1))^2);
        Ixy=Ixy+(Y(cor(i,1))+Y(cor(i,2)))*L(i,i)*X(cor(i,1))*T(i,i)/2;
        Ay=Ay+L(i,i)*T(i,i);
        Sx=Sx+(Y(cor(i,1))+Y(cor(i,2)))*L(i,i)*T(i,i)/2;
        Sy=Sy+X(cor(i,1))*L(i,i)*T(i,i);
        Ixf=Ixf+(FC(cor(i,1))+FC(cor(i,2)))*L(i,i)*T(i,i)*X(cor(i,1))/2;
        Iyf=Iyf+(FC(cor(i,1))*(2*Y(cor(i,1))+Y(cor(i,2)))+FC(cor(i,2))*(2*Y(cor(i,2))+Y(cor(i,1))))*L(i,i)*T(i,i)/6;
        Ssy(i)=(Sy1(i)+Sy2(i))*L(i,i)/T(i,i)/2;
        if sum(Cf(:,i))>0
            Ssx(i)=(Sx1(i)+Sx2(i))*L(i,i)/T(i,i)/2+L(i,i)^3/12;
        elseif sum(Cf(:,i))<0
            Ssx(i)=(Sx1(i)+Sx2(i))*L(i,i)/T(i,i)/2-L(i,i)^3/12;
        else
            for j=1:(Ne-Nv+1)
                if Cf(j,i)<0
                    Sssx(j,i)=(Sx1(i)+Sx2(i))*L(i,i)/T(i,i)/2-L(i,i)^3/12;
                elseif Cf(j,i)>0
                    Sssx(j,i)=-(Sx1(i)+Sx2(i))*L(i,i)/T(i,i)/2+L(i,i)^3/12;
                end
            end
        end
    end
end
A=Ax+Ay;
Nc=[Sy/A Sx/A];
Nn=[Iyf/Ix -Ixf/Iy];
for i=1:(Ne-Nv+1)
    for j=1:Ne
        if Cf(i,j)>0
            Ssssx(i)=Ssssx(i)+Sssx(i,j)+Ssx(j);
            Ssssy(i)=Ssssy(i)+Sssy(i,j)+Ssy(j);
            St(i)=St(i)+L(j,j)/T(j,j);
        elseif Cf(i,j)<0
            Ssssx(i)=Ssssx(i)+Sssx(i,j)+Ssx(j);
            Ssssy(i)=Ssssy(i)+Sssy(i,j)+Ssy(j);
            St(i)=St(i)+L(j,j)/T(j,j);
        end
    end
end
f0x=-(diag(St)-shareborders)\Ssssx;
f0y=-(diag(St)-shareborders)\Ssssy;
Mx=0;My=0;
for i=1:Ne
    Mx=Mx+Ssx(i)*abs(H(i));
    My=My+Ssy(i)*abs(H(i));
end
xw=-(Mx+2*Ac'*f0x)/Ix;
yw=(My+2*Ac'*f0y)/Iy;
Nw=[xw yw];
Axr=-Ax*yw;
Ayr=Ay*xw;

% 计算扭转函数
Kt=Cf*L*inv(T)*Cf';
TFc=Kt\(2*Ac);       %"TF" means "tortion function"
TFe=Cf'*TFc;    

% 计算与布雷特剪应力对应的扭转常数
JB=2*Ac'*TFc;
Ar=Ip-JB;

% 计算与圣维南剪应力对应的扭转常数
JS=I1'*L*T*T*T*I1/3;
J=JB+JS;