C3D8单元B-bar方法 | 公式推导+数值实现

上次推文分享了三维单元的选择性减缩积分技术，该技术将应力分解为体积应力部分和偏应力部分，适用于各向同性的弹性材料，难以扩展到弹塑性材料、轴对称、各向异性、大变形等问题上。

为增加普适性，Hughes    将几何矩阵    分解为体应变部分    和偏应变部分    ，然后对体应变部分    做出某种改进（记为    ），称为B-bar方法，即：

本篇推文将为大家分享三维线性单元的B-bar方法。

公式推导

单元的几何矩阵    可写为：

其中，    表示单元节点个数，子块矩阵    可写为：

其中，    ，    表示与节点    关联的形函数，    表示自由度号，体应变部分    可表示为：

偏应变    ： 

为适应不可压缩材料，体应变部分    需要加一个bar：

   的取值可以为单元形心处的值，即：

本节的程序中采用了形心处的值，读者也可采用单元体积内的平均值：

最终的几何矩阵    可转化为    ：

最终的单元刚度矩阵可写为：

写到这里，程序已经可以顺着公式往下编写了，如果需要对矩阵    的各个元素进行具体细化，可以进一步写为：

其中，

数值案例

由Abaqus帮助文档中：Abaqus->Theory->Elements->Continuum elements->Solid isoparametric quadrilaterals and hexahedra可以看出，为了避免mesh locking，对于轴对称单元、平面应变单元、三维线性单元均做了修正。以致于平时按照课本上的C3D8单元格式编制程序的单元刚度矩阵和Abaqus的C3D8单元刚度矩阵不同。

本小节将B-bar方法修正后的C3D8单元与Abaqus内置的C3D8单元进行做对比，对比内容有位移场、应力场、单元刚度矩阵，模型还是延续上几节的板模型。

   

   

   

   

从上面的云图对比中可以看出，与Abaqus的吻合度极高，那再来看一下刚度矩阵是否与Abaqus一致呢？以第1000号单元刚度矩阵为例：

 

 

可能肉眼看的不是很精确，那就写段程序看一下两个矩阵的元素相差精度：

abaqusMat = Element1000Stiffness; % 第1000号单元刚度矩阵（Abaqus内置的C3D8单元）
matlabMat = solver.ElementMatrix{1000,1}; % 第1000号单元刚度矩阵（自研的C3D8Bar单元）

tolerance = 1e-7; % 设置容差
abs_diff = abs(abaqusMat - matlabMat);
max_diff = max(abs_diff(:));

if max_diff <= tolerance
    disp('矩阵在容差范围内相同');
else
    disp('存在超出容差的差异');
    disp(['最大差异: ', num2str(max_diff)]);
end

最终显示容差为    ，基本可以认为两者矩阵一致。

[1] Hughes, Thomas J. R.The finite element method: linear static and dynamic finite element analysis 2000