ABAQUS C3D10单元质量矩阵计算（附MATLAB代码）
概述
在动力时程分析、频率分析等分析中，质量矩阵的求解是必不可少的。本次帖子以二阶四面体单元的质量矩阵求解为例，讲解有限元法中质量矩阵求解的一些基础知识。采用Matlab编写程序，将计算的质量矩阵与ABAQUS的C3D10单元相比较，最终结算结果是吻合的。帖子末尾给出了求解二阶四面体单元协调质量矩阵的参考代码。
协调/集中质量矩阵
在有限元分析中，质量矩阵的推导通常依赖于以下假设：物体的质量分布是连续的，并且通过离散化单元来逼近实际的质量分布。我们将使用大朗贝尔原理推导单元的质量矩阵。
从运动的结构中取出一个微小部分，根据达朗贝尔原理，在他的单位体积上作用的惯性力为
式中， 为材料的密度。
在对结构进行离散化以后，取出一个单元，并采用如下形式的位移函数
将上面的式子代入到惯性力公式中为
在利用荷载计算公式
即牛顿第二定律
可见，单元的质量矩阵为
特别需要注意的是，在有限元分析中，质量矩阵的形式并不唯一，大致可分为：协调（一致）质量矩阵和集中质量矩阵。
其中，协调（一致）质量矩阵，因为质量的插值公式与位移的插值公式相同，所以被称之为“一致”质量矩阵，又因为该种质量矩阵是满阵，又被称为协调协调质量矩阵。上面的公式即为协调（一致）质量矩阵。
集中质量矩阵，顾名思义，该种质量矩阵只有主对角元素上有非零元素，矩阵的其余位置均为零元素。该种质量矩阵具有很好的数值计算性能，在涉及到矩阵分解的计算中，集中质量矩阵计算非常高效，通常在显式计算中有很广发的运用。
集中质量矩阵的求解方法也是基于协调质量矩阵，比较通用的方法是将协调质量矩阵每一行的质量集中到主对角元素上，同一行的其余元素均人为置零，这只是获得集中质量矩阵的一个简单的方法。
以上只是简单的介绍了一些质量矩阵，上面的公式并不是通用的，在壳等单元中，质量矩阵的求解往往不是一个简单的工作，下面将上述公式用于计算实体单元的质量矩阵。
C3D10单元质量矩阵
这里以二阶四面体单元为例讲解质量矩阵的求解方法，二阶四面体单元即为ABAQUS中的C3D10单元，该单元的几何形状为
二阶单元的形函数为
形函数矩阵为
其中
将上述形函数代入到协调质量矩阵公式中
上面是一个三重积分，获得解析解极为困难，按照有限元的经典流程，下面是采用数值计算方法计算上面的积分。常用的数值计算方法有高斯积分和Hammer积分等等，这里采用Hammer积分进行计算。
在这里要特别注意的是积分点个数的选择，我一开始是按照四个和五个积分点计算，算来算去和abaqus计算的质量矩阵总是对不上，后来翻一下abaqus的官方文档，发现abaqus的C3D10单元质量矩阵采用的是15个积分点。下面是官方文档的原文 For stress/displacement applications the second-order tetrahedron uses 4 integration points for its stiffness matrix and15 integration points for its consistent mass matrix。
文末会给出15积分点的相关程序。采用Hammer积分的公式为
其中， 为Hammer积分权重， 为Jacobian矩阵，Jacobian矩阵的求解可参考我之前发的贴子。
算例
老规矩，还是经典的悬臂梁受动荷载算例，悬臂梁尺寸为 ，密度为2000，泊松比为0.25，弹性模量为 ，模型一端为固定边界，另一端受简谐荷载，荷载表达式为 ，方向为垂直于梁轴线向下，即y轴负方向，计算总时长为10秒，增量步长为0.01秒，共计1000增量步，计算网格为
荷载和边界为
统计悬臂端，即加载面的位移，将abaqus和matlab计算的y方向位移进行对比如下图，结果是一模一样的，数据有震荡是因为没有考虑阻尼。
下面比较一下求解的质量矩阵具体的数值。首先需要修改inp文件，将abaqus计算的质量矩阵输出到外部文件，具体的输出方法参考我之前写的帖子，然后用matlab程序读取质量矩阵并于自编程序计算结果进行比较。
这里的算例一共有696个节点，因此总体质量矩阵大小为2088×2088的方阵，首先我们直接比较质量矩阵对应元素的绝对误差，即 ，并以自由度个数（2088）为横轴绘制图像为
可以发现，绝对误差的数量级在 ，这个误差基本可以认为是不同计算语言或者软件的数据截断误差造成的。
下面按照同样的方法比较对应元素的相对误差，并以自由度个数（2088）为横轴绘制图像为
可以发现，相对误差的数量级来到了 ，基本可以认为，自编的程序计算质量矩阵与abaqus相同。
为了比较两个矩阵的相似程度，下面采用欧几里得范数来比较两个矩阵，matlab提供了norm（v）可以直接调用，我们比较两个矩阵欧几里得范数的相对误差，伪代码为
最终计算的结果为3.061006180118862e-14，可以理解为两个矩阵的欧几里得距离相对误差为3.061006180118862e-14。
代码
下面给出质量矩阵求解过程中用到的关键函数代码。首先是Hammer积分相关代码。再次提醒：ABAQUS计算C3D10单元采用的是Hammer积分的15个积分点。
function[weights,locations] = hammer3D(npt)
if npt==1 % C3D4
locations = [1/4, 1/4, 1/4];
weights = 1/6.;
elseif npt == 4 % C3D10
a=(5+3*sqrt(5))/20;
b=(5-sqrt(5))/20;
locations = [a, b, b;
b, a, b;
b, b, a;
b, b, b];
weights = [1/24, 1/24, 1/24, 1/24];
elseif npt == 5 % 5积分点
a=0.5;
b=1/6;
locations = [0.25,0.25,0.25;
a, b, b;
b, a, b;
b, b, a;
b, b, b];
weights = (1/6)*[-0.8, 9/20, 9/20, 9/20];
elseif npt == 15 % 15积分点
locations =[0.25d0,0.25d0,0.25d0;...
0.091971078052723d0,0.091971078052723d0,0.091971078052723d0;...
0.724086765841831d0,0.091971078052723d0,0.091971078052723d0;...
0.091971078052723d0,0.724086765841831d0,0.091971078052723d0;...
0.091971078052723d0,0.091971078052723d0,0.724086765841831d0;...
0.319793627829630d0,0.319793627829630d0,0.319793627829630d0;...
0.040619116511110d0,0.319793627829630d0,0.319793627829630d0;...
0.319793627829630d0,0.040619116511110d0,0.319793627829630d0;...
0.319793627829630d0,0.319793627829630d0,0.040619116511110d0;...
0.056350832689629d0,0.056350832689629d0,0.443649167310371d0;...
0.443649167310371d0,0.056350832689629d0,0.056350832689629d0;...
0.443649167310371d0,0.443649167310371d0,0.056350832689629d0;...
0.056350832689629d0,0.443649167310371d0,0.443649167310371d0;...
0.056350832689629d0,0.443649167310371d0,0.056350832689629d0;...
0.443649167310371d0,0.056350832689629d0,0.443649167310371d0];
weights =[0.019753086419753d0,0.011989513963170d0,0.011989513963170d0,...
0.011989513963170d0,0.011989513963170d0,0.011511367871045d0,...
0.011511367871045d0,0.011511367871045d0,0.011511367871045d0,...
0.008818342151675d0,0.008818342151675d0,0.008818342151675d0,...
0.008818342151675d0,0.008818342151675d0,0.008818342151675d0];
end
end
下面是C3D10单元的形函数以及形函数对局部坐标的偏导数。
function [shape,naturalDerivatives] = shapeFunction3D(xi,eta,gamma,elenode)
switch elenode
case 10
shape = [
(2*(1-xi-eta-gamma)-1)*(1-xi-eta-gamma); % N1
(2*xi-1)*xi; % N2
(2*eta-1)*eta; % N3
(2*gamma-1)*gamma; % N4
4*(1-xi-eta-gamma)*xi; % N5
4*xi*eta; % N6
4*(1-xi-eta-gamma)*eta; % N7
4*(1-xi-eta-gamma)*gamma; % N8
4*xi*gamma; % N9
4*eta*gamma % N10
];
naturalDerivatives = [
-4*(1-xi-eta-gamma)+1, -4*(1-xi-eta-gamma)+1, -4*(1-xi-eta-gamma)+1; % dN1/dxi, dN1/deta, dN1/dgamma
4*xi-1, 0, 0; % dN2/dxi, dN2/deta, dN2/dgamma
0, 4*eta-1, 0; % dN3/dxi, dN3/deta, dN3/dgamma
0, 0, 4*gamma-1; % dN4/dxi, dN4/deta, dN4/dgamma
4*(1-2*xi-eta-gamma), -4*xi, -4*xi; % dN5/dxi, dN5/deta, dN5/dgamma
4*eta, 4*xi, 0; % dN6/dxi, dN6/deta, dN6/dgamma
-4*eta, 4*(1-xi-2*eta-gamma), -4*eta; % dN7/dxi, dN7/deta, dN7/dgamma
-4*gamma, -4*gamma, 4*(1-xi-eta-2*gamma); % dN8/dxi, dN8/deta, dN8/dgamma
4*gamma, 0, 4*xi; % dN9/dxi, dN9/deta, dN9/dgamma
0, 4*gamma, 4*eta % dN10/dxi, dN10/deta, dN10/dgamma
];
end
end
