sectioninput.m

薄壁结构有限元计算

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

% 截面端点数
Nv=input('Please input the number of points of the section: ');

% 截面边数
Ne=input('Please input the number of borders of the section: ');

% 截面关联矩阵
B=zeros(Nv,Ne);
Be=input('Please input the conjuction matrix of the section(2*Ne): ');
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=input('Please input the basic loop matrix of the section(+1): ');
Cf2=input('Please input the basic loop matrix of the section(-1): ');
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=input('Please input the x-coordinate of points of the section: ');X=X';
Y=input('Please input the y-coordinate of points of the section: ');Y=Y';
T=input('Please input the thickness of borders of the section: ');T=diag(T);

% 截面边的长度
L=input('Please input the length matrix of borders of the section: ');L=diag(L);

% 截面的边到扭转中心的距离
H=input('Please input the distance array from borders to tortion center of the section: ');H=H';

% 横竖边的判别
direction=input('Please input the array of direction of borders: ');

% 截面回路围绕面积
Ac=input('Please input the area array of loops of the section: ');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;