详解有限元法中K、M和C矩阵的集成

概述：该帖子重点讲解有限元法中单元刚度矩阵、质量矩阵和阻尼矩阵如何组装为模型的总体刚度矩阵。首先采用MATLAB程序自编三维八节点单元，即对标ABAQUS中的C3D8单元，采用了B-BAR修正以后，计算的结果与ABAQUS保持一致。然后修改INP文件，将ABAQUS的模型总体矩阵导出为可读文件，再读进MATLAB中，然后用ABAQUS的总体矩阵与自己计算的总体矩阵作对比。设计悬臂梁一端受动荷载作用算例，自编程序采用newmark积分法，ABAQUS采用动力隐式分析步，对比模型关键点的位移、速度和加速度数据，最终计算结果保持一致。附件包含帖子中的ABAQUS算例，有导出刚度矩阵、质量矩阵和阻尼矩阵的INP文件。还有MATLAB程序，ABAQUS导出的矩阵不能直接读取使用，因此，MATLAB程序用于处理ABAQUS导出的矩阵。

（）矩阵的组装原理及程序讲解

勘误：右边大矩阵中的“3p-2 3p-1 p”应为：“3p-2 3p-1 3p”。

如上图所示，左边的小矩阵为单元刚度矩阵，假定为k，假定该单元刚度矩阵的节点编号为1~8，则单刚的大小为24x24，第P个节点的元素如图所示，则该节点元素在单元刚度矩阵中的位置是k(i:k,i:k)；右边的矩阵为总体刚度矩阵，假定为K，假定模型的总节点数为N，则总体刚度矩阵的大小为K(3N，3N)，则P节点对应的元素在总体刚度矩阵中的位置为：K(3P-2:3P,3P-2:3P)，这样就建立了单元刚度矩阵中元素对应在总体刚度矩阵中的位置。下面给出了三种组装的方法，分别对应了不同教材中的理论，注释注明了来源。需要注意的是，对于质量矩阵和阻尼矩阵，有相同的组装方式。

对于下列程序，element为整个模型的单元节点编号，如第n个单元的节点编号为：element（n,:）,kk则储存了模型所有单元的单刚，为三维数组，如第n个单元的单刚为：k(:,:,n)。

第一种组装方法依据的理论是将单刚节点对应的3x3矩阵放到总体刚度矩阵中，具体的程序如下：

function KK=kkAssembly(element,kk)
%  组装总刚
max_node=max(max(element));
KK=zeros(3*max_node,3*max_node);
for i=1:size(element,1)
    KK=KK+formKK(element(i,:),kk(:,:,i),max_node);
end
end
%% 
% 按照自由度对应的分块组装方法，计算速度还可以
function KK=formKK(element,kk,max_node)
% 计算单元节点在全局刚度矩阵中的位置
KK=zeros(3*max_node,3*max_node);
for i=1:8
    for ii=1:8
        node1=element(1,i+1);        
            node2=element(1,ii+1);
        KK(3*node1-2:3*node1,3*node2-2:3*node2)=kk(3*i-2:3*i,3*ii-2:3*ii);
    end
end
end

第二种组装方法的理论如下图：

具体的程序如下：

function KK=kkAssembly(element,kk)
%  组装总刚
max_node=max(max(element));
KK=zeros(3*max_node,3*max_node);
for i=1:size(element,1)
    KK=KK+formKK(element(i,:),kk(:,:,i),max_node);
end
end
%% 采用矩阵转换的方法组装，算到猴年马月也算不完！计算速度令人发指！该方法适用于教学
% 王勖成老师的《有限单元法》，清华大学出版社，第69页，式（2.2.53）
% function KK=formKK(element,kk,max_node)
% G=zeros(8,2*max_node); 
% KK=zeros(2*max_node,2*max_node);
% for i=1:4
%     G(2*i-1:2*i,2*element(1,i+1)-1:2*element(1,i+1))=[1 0;0 1];
% end
% KK=G'*kk*G;
% end

第三种组装方法理论为：

具体的程序为：

function KK=kkAssembly(element,kk)
%  组装总刚
max_node=max(max(element));
KK=zeros(3*max_node,3*max_node);
for i=1:size(element,1)
    KK=KK+formKK(element(i,:),kk(:,:,i),max_node);
end
end
%% 
%% 采用编码法，朱院士的《有限单元法原理与应用》，第四版，第55页，表2-1
% function KK=formKK(element,kk,max_node)
% % 计算单元节点在全局刚度矩阵中的位置
% KK=zeros(2*max_node,2*max_node);
% Globalocation=zeros(1,size(element,2)-1);
% for i=1:size(element,2)-1
%     Globalocation(1,2*i-1)=2*element(1,i+1)-1;
%     Globalocation(1,2*i)=2*element(1,i+1);
% end
% for i=1:8
%     for ii=1:8
%         KK(Globalocation(1,i),Globalocation(1,ii))=kk(i,ii);
%     end
% end
% end
%%

（）修改INP文件导出模型总体刚度矩阵、单元刚度矩阵

**include,input=mo.inp
************************************************************************************
***输出单元刚度矩阵
*Step, name=Emmkk
*static
 1.,1.,1e-05
*File Format, ASCII
*Element Matrix Output, Elset=Part-1-1.SET-1, File Name=EMass, Output File=User Defined, mass=yes
*Element Matrix Output, Elset=Part-1-1.SET-1, File Name=EStiffness, Output File=User Defined, stiffness=yes
*End Step
************************************************************************************
***输出总体刚度矩阵
*Step, name=Gkk
*MATRIX GENERATE, STIFFNESS
*MATRIX OUTPUT, STIFFNESS, FORMAT=COORDINATE
*End Step
************************************************************************************
***输出总体质量矩阵，集中质量阵
*Step, name=Gmm
*MATRIX GENERATE, mass
*MATRIX OUTPUT, mass, FORMAT=COORDINATE
*End Step	
************************************************************************************
***输出总体阻尼矩阵，瑞丽阻尼
*Step, name=Gcc
*MATRIX GENERATE, VISCOUS DAMPING
*MATRIX OUTPUT, VISCOUS DAMPING, FORMAT=COORDINATE
*End Step		

（）算例介绍

悬臂梁尺寸为10x10x40，弹性模量为1e10，泊松比为0.25，密度为2400，一端固定，另一端节点施加简谐动荷载，边界条件和荷载示意图为：

简谐荷载时程曲线为：

模型网格为：

对该计算模型设计三种工况作对比，三种工况的信息分别为：

1. ABAQUS计算，采用动力隐式分析步，计算总时长为10，固定增量步长为0.01，总增量步数为1000。记为工况一。
2. 采用MATLAB完全自编C3D8单元，读取ABAQUS输出的单元节点和节点坐标信息，计算单元刚度矩阵（考虑B-BAR修正）、质量矩阵（集中质量矩阵）和阻尼矩阵（比例阻尼）并集成为总体矩阵，采用newmark时程积分法，计算总时长为10，固定增量步长为0.01，总增量步数为1000。记为工况二。
3. 导出ABAQUS的总体矩阵，替换掉自编程序中的总体刚度矩阵，采用newmark时程积分法，计算总时长为10，固定增量步长为0.01，总增量步数为1000。记为工况三。

（）模型总体刚度矩阵、质量矩阵和阻尼矩阵对比

首先是自编程序计算的总体刚度矩阵数值的部分截图如下：

然后是ABAQUS导出的总体刚度矩阵数值的部分截图如下：

可以发现，计算结果是一致的。将两个矩阵做差记为KK，计算KK矩阵的1、2和无穷范数，结果分别为：0.000644803047180176、0.000345101701596316和0.000643551349639893，计算结果说明自编程序组装的总体刚度矩阵和ABAQUA导出的矩阵是吻合的。质量矩阵他同样是相同的，具体见附件。单刚和质量计算无误，则阻尼一定是准确的。

（）位移、速度和加速度数据对比

提取悬臂端1节点编号位置的位移时程曲线如下图：

提取悬臂端1节点编号位置的速度时程曲线如下图：

结果分析：

• 在计算精度方面，三种工况的计算结果是一致的。说明，自编C3D8单元能够达到ABAQUS的计算精度，位移数据吻合，说明后续的应变和应力数据也是吻合的，但是目前没有编制相应的程序做验证。而且说明了单元矩阵组装为模型总体刚度矩阵的程序是完全正确的。
• 在计算时长方面，目前的程序比ABAQUS计算的要快得多，但是这并不意味着自己编制的程序计算效率就一定超越了ABAQUS，这是不可能的。究其原因，是因为自己编制的程序只是简单的计算了最基础的数据，而且采用了最简单的饿newmark时程积分法，在增量步之间没有进行任何的收敛判断计算，这在线弹性小变形问题中是可以的，但是其他情况不适用；但是ABAUQS中的计算流程则极为复杂，包括计算前的矩阵预处理、矩阵求逆以及为了保证计算精度的处理、增量步之间的收敛判断等等等等，这些通用的计算步骤，占据了简单算例的大部分时间。

（）附件