《有限元基础编程百科全书》| C3D8单元（理论、编程、数值验证）

C3D8单元的构造可先回顾平面等参单元的构造方式：《有限元基础编程百科全书》| 平面等参四边形单元，只是单纯的在二维的基础上增加一个维度，本节介绍的是未做任何修正的C3D8单元，对于Abaqus内置的C3D8单元将会另起一小节单独介绍。

单元介绍

如下图所示的六面体单元（等参坐标下），每个节点三个平动自由度，共8个节点，也就是24个自由度，单元的节点位移列阵    和节点力列阵    为：

 

有限元格式

形函数

位移场

单元内某一点的位移，可通过形函数与节点位移描述。

单元内位移场可描述为：

几何方程

单元内某一点的应变，可通过几何矩阵     与节点位移描述。

其中，

本构矩阵

单元刚度方程

C3D8单元的高斯积分点位置与平面单元一致，只是在原有两个方向上再增加一个维度，一维情况下高斯积分点及权重系数汇总见下表：

在对单元刚度数值积分时可对上表的积分点进行三重遍历，下面代码是我的处理方式，仅供参考：

for ig = 1:ngp
    xi_Gauss = samp(ig,1); % 第 ig 个积分点的 xi 坐标
    wi = samp(ig,2);       % 第 ig 个积分点的权重
        for jg = 1:ngp
            eta_Gauss = samp(jg,1); % 第 jg 个积分点的 eta 坐标
            wj = samp(jg,2);        % 第 jg 个积分点的权重
            for kg = 1:ngp
            gama_Gauss = samp(kg,1); % 第 kg 个积分点的 gamma 坐标
            wk = samp(kg,2);         % 第 kg 个积分点的权重
            
            % 计算 B 矩阵和 Jacobian 矩阵
            [B, Jacob] = Bmatrix3D(elenode, elemNodeCoordinate, DOF, xi_Gauss, eta_Gauss, gama_Gauss);
            
            % 使用高斯积分公式更新刚度矩阵
            k = k + B.' * D * B * det(Jacob) * wi * wj * wk;
        end
    end
end

数值案例

考虑如下的数值算例模型，底部固定三个方向的自由度，内壁受到大小为100的表面压强，弹性模量取1000，泊松比取0.3，单元类型为C3D8。分别使用Abaqus和MFEAOOP进行计算，对比位移场和应力场，验证其程序的数值精度。

 

精度验证

将位移场和应力场云图罗列如下：

从上面的结果可以得知，MFEAOOP的C3D8单元与Abaqus的内置C3D8单元大致趋势相同，但是数值方面有些许差距，因为大多数的商软都会对线性单元做一些各方面的修正以适应更多复杂的工况，在后续小节将会分享Abaqus对C3D8单元究竟做了哪些修正，敬请期待！