transientmasselement.m

薄壁结构有限元计算

% 质量阵的计算与集成

% 输入结构物划分层数、单元质量、各转动惯量、质心位置坐标（其中高程方向为变化的数组）
n=input('Please input the number of layers of the construction: ');
Mem=input('Please input the mass of the element: ');
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: ');
Iz=input('Please input the z-direction of inertia-moment of the section: ');
Ixy=input('Please input the xy-direction of inertia-moment of the section: ');
Iyz=input('Please input the yz-direction of inertia-moment of the section: ');
Ixz=input('Please input the xz-direction of inertia-moment of the section: ');
xc=input('Please input the x-coordinate of centroid of the section: ');
yc=input('Please input the y-coordinate of centroid of the section: ');
zc=input('Please input the z-coordinate array of centroid of the section(1×(n+1)): ');

% 先将总质量阵和单元质量阵分配零空间
M=zeros(7*(n+1),7*(n+1));
Me=zeros(7,7,(n+1));

% 计算第一层和最后一层的单元质量阵，并将值赋予总质量阵中
Me(:,:,1)=[Mem/2 0 0 -Mem/2*yc Mem/2*zc(1) 0 0;
            0 Mem/2 0 Mem/2*xc 0 -Mem/2*zc(1) 0;
            0 0 Mem/2 0 Mem/2*xc Mem/2*yc 0;
            -Mem/2*yc Mem/2*xc 0 Ix+Iy Iyz -Ixz 0;
            Mem/2*zc(1) 0 -Mem/2*xc -Iyz -Ix-Iz -Ixy 0;
            0 -Mem/2*zc(1) Mem/2*yc -Ixz Ixy Iy+Iz 0;
            0 0 0 0 0 0 0];
M(1:7,1:7)=Me(:,:,1);
Me(:,:,(n+1))=[Mem/2 0 0 -Mem/2*yc Mem/2*zc(n+1) 0 0;
            0 Mem/2 0 Mem/2*xc 0 -Mem/2*zc(n+1) 0;
            0 0 Mem/2 0 Mem/2*xc Mem/2*yc 0;
            -Mem/2*yc Mem/2*xc 0 Ix+Iy Iyz -Ixz 0;
            Mem/2*zc(n+1) 0 -Mem/2*xc -Iyz -Ix-Iz -Ixy 0;
            0 -Mem/2*zc(n+1) Mem/2*yc -Ixz Ixy Iy+Iz 0;
            0 0 0 0 0 0 0];
M((7*n+1):(7*n+7),(7*n+1):(7*n+7))=Me(:,:,(n+1));

% 计算中间层的单元质量阵，并将值赋予总质量阵中
for i=2:n
 Me(:,:,i)=[Mem 0 0 -Mem*yc Mem*zc(i) 0 0;
            0 Mem 0 Mem*xc 0 -Mem*zc(i) 0;
            0 0 Mem 0 Mem*xc Mem*yc 0;
            -Mem*yc Mem*xc 0 Ix+Iy Iyz -Ixz 0;
            Mem*zc(i) 0 -Mem*xc -Iyz -Ix-Iz -Ixy 0;
            0 -Mem*zc(i) Mem*yc -Ixz Ixy Iy+Iz 0;
            0 0 0 0 0 0 0];
 M((7*(i-1)+1):(7*(i-1)+7),(7*(i-1)+1):(7*(i-1)+7))=Me(:,:,i);
end