stiffness.m

薄壁结构有限元计算

% 集成刚度矩阵
K=zeros(7*n,7*n);
for i=1:(n-1)
    K=BendTorsionStiffnessAssemble(K,Ke3,i,i+1);
end
for j=1:7
    for i=1:7
        K(i,j)=K(i,j)*2;
    end
end
for i=1:7*n
    for j=1:7*n
        if abs(K(i,j))<(1.0e-8)*max(max(abs(K)))
            K(i,j)=0;
        end
    end
end

% 集成质量矩阵
M=zeros(7*n,7*n);
mem=rho*A*l;
Ixx=0;Iyy=0;
for i=1:Ne
    if direction(i)>0
        Iyy=Iyy+L(i,i)*T(i,i)*(L(i,i)^2+4*l^2);
        Ixx=Ixx+L(i,i)*T(i,i)*(12*(H(i)^2)+l^2);
    else
        Ixx=Ixx+L(i,i)*T(i,i)*(L(i,i)^2+4*l^2);
        Iyy=Iyy+L(i,i)*T(i,i)*(12*(H(i)^2)+l^2);
    end
end
Ixx=Ixx*rho*l/12;
Iyy=Iyy*rho*l/12;
for k=1:n
    M(7*k-6,7*k-6)=mem;
    M(7*k-5,7*k-5)=mem;
    M(7*k-4,7*k-4)=mem;
    M(7*k-3,7*k-3)=(Ix+Iy)*rho*l;
    M(7*k-2,7*k-2)=Ixx;
    M(7*k-1,7*k-1)=Iyy;
    M(7*k,7*k)=If*rho*l;
end