三维非协调元C3D8I理论介绍+数值实现

形函数修正

在4节点平面非协调元的基础上，还可以继续推广至8节点和20节点的三维单元以及16节点的厚壳单元中。对于8节点三维非协调单元(C3D8I)的形函数可扩充为：

位移模式变为：

式中的    为附加的自由度，新引入的形函数在单元边界上呈二次抛物线变化，附加的位移项在单元与单元的交界面处是不连续的，也叫做不协调的，称非协调模式（incompatible model），相应的单元称非协调单元。从数学上看，通过引入    、    、    项，使得形函数中的二次项趋于完备，从而提高计算精度。非协调单元不满足有限元解的收敛性，可以通过分片试验来检验解的收敛性。

刚度方程

在求解单元刚度矩阵时，应变关系矩阵只需在原有格式上扩充三列即可：

式中， 

非协调模式不额外增加单元节点力，单元刚度矩阵可写为：

由静力凝聚技术可得非协调元的单元刚度矩阵为：

单元改进

以上单元刚度矩阵的形式，在进行分片试验时，发现雅可比行列式为常数时，即单元体积为一个定值，单元是收敛的，反之不收敛。为了解决这一问题，Taylor等人建议在计算    时，雅可比矩阵取单元形心处（    ）的值。具体操作呢可点击跳转平面非协调单元详解，也可以结合具体代码C3D8I.m进行学习。

数值案例

考虑如下图所示的悬臂梁模型，模型左端节点全部自由度固定，在模型上表面施加大小为10的表面压强，弹性模量    ，泊松比    ，厚度方向布置一层单元。

 

将自研C3D8I单元程序与Abaqus自带的C3D8I单元做位移精度对比：

   

   

从上图可看出自研的程序与Abaqus自带的C3D8I单元精度吻合较好，用户可在此基础上开展进一步研究。