transientstiffnesselement.m

薄壁结构有限元计算

% 输入转动惯量（x向、y向、二向混合、翘曲、x向翘曲、y向翘曲）
% Ix=input('Please input the x-direction of inertia-moment of the section: ');
% Iy=input('Please input the y-direction of inertia-moment of the section: ');
% Ixy=input('Please input the xy-direction of inertia-moment of the section: ');
% If=input('Please input the warp of inertia-moment of the section: ');
% Ixf=input('Please input the x-direction warp of inertia-moment of the section: ');
% Iyf=input('Please input the y-direction warp of inertia-moment of the section: ');

% 输入面积（截面、x向剪切、y向剪切、混合剪切、x向翘曲剪切、y向翘曲剪切、翘曲剪切）
% A=input('Please input the area of the section: ');
% Ax=input('Please input the x-direction of shear-area of the section: ');
% Ay=input('Please input the y-direction of shear-area of the section: ');
% Axy=input('Please input the mixed shear-area of the section: ');
Axy=0;
% Axr=input('Please input the x-direction warp of shear-moment of the section: ');
Axr=Ax*abs(abs(Nn(1))-abs(Nc(1)));
% Ayr=input('Please input the y-direction warp of shear-moment of the section: ');
Ayr=Ay*abs(abs(Nn(2))-abs(Nc(2)));
% Ar=input('Please input the warp of shear-moment of the section: ');

% 输入单元长度
l=input('Please input the length of the element: ');

% 输入单元层数
n=input('Please input the number of layers: ');

% 输入x向静矩、y向静矩、弹性模量、剪切模量、扭转常数、密度、固有频率
% Sx=input('Please input the x-direction of shear-moment of the section: ');
% Sy=input('Please input the y-direction of shear-moment of the section: ');
E=input('Please input the elastic modulus of the material: ');
mu=input('Please input the poisson ratio: ');G=0.5*E/(1+mu);
% J=input('Please input the torsion-constant of the section: ');
rho=input('Please input the density of the material: ');
Aa=input('Please input the area of the structure: ');
omega=1.619*(E*min(Ix,Iy)/rho/Aa/(n*l)^4)^0.5;

% 计算各参数与AA矩阵
b2=rho*omega*omega*A;
b3=rho*omega*omega*Sy;
b4=rho*omega*omega*Sx;
b5=G*Ax;
b6=G*Axy;
b7=G*Axr;
b8=-G*Ax+rho*omega*omega*Iy;
b9=G*Axy-rho*omega*omega*Ixy;
b10=G*Axr-rho*omega*omega*Iyf;
b11=G*Ay;
b12=G*Ayr;
b13=-G*Ay+rho*omega*omega*Ixy;
b14=-G*Ayr+rho*omega*omega*Ixf;
b15=rho*omega*omega*If;
b16=G*Ar;
b17=-G*Ar+rho*omega*omega*If;
AA=[0 1 0 0 0 0 0 0 0 0 0 0 0 0;
   b2/E/A 0 0 0 -b3/E/A 0 0 0 b4/E/A 0 0 0 0 0;
   0 0 0 1 0 0 0 0 0 0 0 0 0 0;
   0 0 b2/G/Ax 0 0 -b5/G/Ax 0 0 0 b6/G/Ax -b4/G/Ax 0 0 b7/G/Ax;
   0 0 0 0 0 1 0 0 0 0 0 0 0 0;
   -b3/E/Iy 0 0 b5/E/Iy b8/E/Iy 0 0 b6/E/Iy b9/E/Iy 0 0 b7/E/Iy b10/E/Iy 0;
   0 0 0 0 0 0 0 1 0 0 0 0 0 0;
   0 0 0 0 0 -b6/G/Ay b2/G/Ay 0 0 b11/G/Ay b3/G/Ay 0 0 b12/G/Ay;
   0 0 0 0 0 0 0 0 0 1 0 0 0 0;
   b4/E/Ix 0 0 -b6/E/Ix b9/E/Ix 0 0 -b11/E/Ix b13/E/Ix 0 0 -b12/E/Ix b14/E/Ix 0;
   0 0 0 0 0 0 0 0 0 0 0 1 0 0;
   0 0 -b4/G/(J+Ar) 0 0 -b7/G/(J+Ar) b3/G/(J+Ar) 0 0 b12/G/(J+Ar) b15/G/(J+Ar) 0 0 b16/G/(J+Ar);
   0 0 0 0 0 0 0 0 0 0 0 0 0 1;
   0 0 0 -b7/E/If b10/E/If 0 0 -b12/E/If b14/E/If 0 0 -b16/E/If b17/E/If 0];
for i=1:14
    for j=1:14
        if abs(AA(i,j))<(3.0e-8)*max(max(abs(AA)))
            AA(i,j)=0;
        end
    end
end

% 计算各参数与BB矩阵
k2=-E*A;
k3=E*Sy;
k4=-E*Sx;
k5=-G*Ax;
k6=-G*Axy;
k7=-G*Axr;
k8=-E*Iy;
k9=E*Ixy;
k10=E*Iyf;
k11=-G*Ay;
k12=-G*Ayr;
k13=-E*Ix;
k14=-E*Ixf;
k15=-G*J-G*Ar;
k16=-E*If;
BB=[1 0 0 0 0 0 0 0 0 0 0 0 0 0;
   0 k2/E/A 0 0 0 k3/E/A 0 0 0 k4/E/A 0 0 0 0;
   0 0 1 0 0 0 0 0 0 0 0 0 0 0;
   0 0 0 k5/G/Ax 0 0 0 k6/G/Ax 0 0 0 k7/G/Ax 0 0;
   0 0 0 0 1 0 0 0 0 0 0 0 0 0;
   0 k3/E/Iy 0 0 0 k8/E/Iy 0 0 0 k9/E/Iy 0 0 0 k10/E/Iy;
   0 0 0 0 0 0 1 0 0 0 0 0 0 0;
   0 0 0 k6/G/Ay 0 0 0 k11/G/Ay 0 0 0 k12/G/Ay 0 0;
   0 0 0 0 0 0 0 0 1 0 0 0 0 0;
   0 k4/E/Ix 0 0 0 k9/E/Ix 0 0 0 k13/E/Ix 0 0 0 k14/E/Ix;
   0 0 0 0 0 0 0 0 0 0 1 0 0 0;
   0 0 0 k7/G/(J+Ar) 0 0 0 k12/G/(J+Ar) 0 0 0 k15/G/(J+Ar) 0 0;
   0 0 0 0 0 0 0 0 0 0 0 0 1 0;
   0 0 0 0 0 k10/E/If 0 0 0 k14/E/If 0 0 0 k16/E/If];
for i=1:14
    for j=1:14
        if abs(BB(i,j))<(1.0e-5)*max(max(abs(BB)))
            BB(i,j)=0;
        end
    end
end

% 计算特征值与特征向量
[Z,lambda]=eig(AA,BB,'qz');

% 形成EE矩阵
EE=[Z(1,:);Z(3,:);Z(5,:);Z(7,:);Z(9,:);Z(11,:);Z(13,:);[[Z(1,:);Z(3,:);Z(5,:);Z(7,:);Z(9,:);Z(11,:);Z(13,:)]*diag([exp(l*lambda(1,1)) exp(l*lambda(2,2)) exp(l*lambda(3,3)) exp(l*lambda(4,4)) exp(l*lambda(5,5)) exp(l*lambda(6,6)) exp(l*lambda(7,7)) exp(l*lambda(8,8)) exp(l*lambda(9,9)) exp(l*lambda(10,10)) exp(l*lambda(11,11)) exp(l*lambda(12,12)) exp(l*lambda(13,13)) exp(l*lambda(14,14))])]];
for i=1:14
    for j=1:14
        if abs(EE(i,j))<(1.0e-6)*max(max(abs(EE)))
            EE(i,j)=0;
        end
    end
end

% 计算各参数与S矩阵
s1=E*A;
s2=-E*Sy;
s3=E*Sx;
s4=G*Ax;
s5=-G*Ax;
s6=G*Axy;
s7=G*Axr;
s8=G*Axr;
s9=E*Iy;
s10=-E*Ixy;
s11=0;
s12=-E*Iyf;
s13=G*Ay;
s14=G*Ay;
s15=G*Ayr;
s16=G*Ayr;
s17=E*Ix;
s18=0;
s19=E*Ixf;
s20=-G*Axr;
s21=G*Ayr;
s22=G*J+G*Ar;
s23=G*Ar;
s24=E*If;
S=[0 s1 0 0 0 s2 0 0 0 s3 0 0 0 0;
   0 0 0 s4 s5 0 0 s6 s6 0 0 s7 s8 0;
   0 s2 0 0 0 s9 0 0 0 s10 s11 0 0 s12;
   0 0 0 s6 -s6 0 0 s13 s14 0 0 s15 s16 0;
   0 s3 0 0 0 s10 0 0 0 s17 s18 0 0 s19;
   0 0 0 s7 s20 0 0 s15 s21 0 0 s22 s23 0;
   0 0 0 0 0 s12 0 0 0 s19 0 0 0 s24];
for i=1:7
    for j=1:14
        if abs(S(i,j))<(1.0e-6)*max(max(abs(S)))
            S(i,j)=0;
        end
    end
end

% 计算X0与Xl矩阵
X0=Z;
for i=1:14
    for j=1:14
        if abs(X0(i,j))<(1.0e-9)*max(max(abs(X0)))
            X0(i,j)=0;
        end
    end
end
Xl=[exp(lambda(1,1)*l)*Z(:,1) exp(lambda(2,2)*l)*Z(:,2) exp(lambda(3,3)*l)*Z(:,3) exp(lambda(4,4)*l)*Z(:,4) exp(lambda(5,5)*l)*Z(:,5) exp(lambda(6,6)*l)*Z(:,6) exp(lambda(7,7)*l)*Z(:,7) exp(lambda(8,8)*l)*Z(:,8) exp(lambda(9,9)*l)*Z(:,9) exp(lambda(10,10)*l)*Z(:,10) exp(lambda(11,11)*l)*Z(:,11) exp(lambda(12,12)*l)*Z(:,12) exp(lambda(13,13)*l)*Z(:,13) exp(lambda(14,14)*l)*Z(:,14)];
for i=1:14
    for j=1:14
        if abs(Xl(i,j))<(1.0e-9)*max(max(abs(Xl)))
            Xl(i,j)=0;
        end
    end
end

% 形成单元刚度矩阵
Ke=[-S*X0*inv(EE);S*Xl*inv(EE)];
Ke1=real(Ke);
Ke2=[Ke1(:,1) Ke1(:,2) Ke1(:,4) Ke1(:,6) Ke1(:,5) Ke1(:,3) Ke1(:,7) Ke1(:,8) Ke1(:,9) Ke1(:,11) Ke1(:,13) Ke1(:,12) Ke1(:,10) Ke1(:,14)];
Ke3=[Ke2(1,:);Ke2(2,:);Ke2(4,:);Ke2(6,:);Ke2(5,:);Ke2(3,:);Ke2(7,:);Ke2(8,:);Ke2(9,:);Ke2(11,:);Ke2(13,:);Ke2(12,:);Ke2(10,:);Ke2(14,:)];
for i=1:14
    for j=1:14
        if abs(Ke3(i,j))<(1.0e-8)*max(max(abs(Ke3)))
            Ke3(i,j)=0;
        end
    end
end