ABAQUS CPE4 单元的默认 B-bar method
提到 CPE4 单元，我们通常会想到：四节点、双线性、平面应变、完全积分。
但这里的“完全积分”，并不意味着直接将标准应变—位移矩阵代入所有积分项，不作任何修正。
实际上，Abaqus 对这类一阶完全积分等参单元采用了选择性减缩积分技术，也称为 B-bar 方法。[1,2]
它的核心思想可以概括为：
保留面内偏量部分，将各积分点的面内体积应变替换为单元平均值。
本文从应变分解出发，推导 CPE4 单元的 B-bar 矩阵，并解释三个容易混淆的问题：
• 为什么修正系数是 ，而不是 ？ • 为什么修正后仍然满足 ？ • 为什么使用中心值，并不意味着整个单元采用单点积分？
本文的讨论范围
以下采用小变形、小应变框架，讨论常规 CPE4 单元的 B-bar 应变修正。
这些表达式用于说明其基本原理，不应直接视为有限变形情况下完整的单元切线公式。
1. 从标准应变—位移关系开始
1.1 节点位移与应变排列
用大写 表示节点编号：
单元节点位移向量写为
采用如下应变排列：
其中，
为工程剪应变。
平面应变条件要求
平面应变不等于平面应力
平面应变条件约束的是面外应变，而不是面外应力。因此， 并不意味着 。
1.2 标准节点子矩阵
对于节点 ，记形函数的空间导数为
与该节点对应的标准应变—位移子矩阵为
将四个节点的子矩阵拼接，可得
因此，标准应变—位移关系为
这里， 的尺寸为 ，完整矩阵 的尺寸为 。
2. B-bar 方法修正的是什么？
定义面内法向应变之和
由于平面应变条件下 ，它也等于小应变意义下的体积应变：
对于位移插值得到的应变场，
不同积分点处的形函数导数通常不同，因此，各积分点处的 也可能不同。
B-bar 方法不直接使用这些各自变化的体积应变，而是将它们替换为同一个单元平均值。[1]
这相当于放松单元内部对体积变形施加的过多约束，从而缓解体积锁死。[1]
修正的是体积变形的表达方式，而不是把整个应变场都变成常量。
3. 面内分解：为什么系数是 1/2？
将标准节点子矩阵分解为
这里， 和 分别表示本文定义的面内偏量部分与面内体积部分。
3.1 面内体积部分
将 平均分配到 和 两个方向，每个方向的贡献为 。
相应的节点子矩阵为
它只影响两个面内法向应变，不影响面外应变和工程剪应变。
3.2 面内偏量部分
由
可得
这部分包含面内法向应变差与剪切变形。
3.3 不要与三维球量—偏量分解混淆
三维应变张量的球量部分为
其中 为三维单位张量。
但本文采用的是二维面内分解：修正只分配给 和 两个方向，因此系数为 。
Abaqus 理论文档明确区分了这两种情况：二维单元的修正系数为 ，三维单元为 。[1]
1/2 的真正含义
这里采用 ，不是说三维球量分解中的 不再成立，而是因为 CPE4 的这一应变修正仅作用于面内两个方向。
因此，不能直接把这里的面内偏量部分等同于三维应变张量的完整偏量部分。
4. 从单元平均值到中心值
4.1 单元平均的定义
记单元面积为 。对于恒定厚度的平面单元，体积平均可以写成面积平均：
相应的平均体积应变为
将位移插值代入，可得
4.2 Q4 单元的中心值表达
对于常规直边、双线性等参 Q4 平面单元，在恒定厚度、有效映射的条件下，形函数空间导数的单元面积平均值等于其自然坐标中心值。[1]
记
于是有
以及
因此，对于本文讨论的 Q4 平面单元，用单元平均值替换体积部分，可以写成用自然坐标中心值替换体积部分。
“中心”指哪里？
这里的中心特指自然坐标中的 ，不应笼统地理解为任意几何意义下的形心。
“平均值等于中心值”也不宜不加条件地推广到三维六面体或其他类型的单元。
5. B-bar 矩阵的构造
B-bar 修正保留当前位置的面内偏量部分，并用单元平均值替换面内体积部分：
其中，
对于本文讨论的 Q4 单元，也可以写成
等价地，
将各项展开，得到
其中，不带上标 的导数在当前积分点计算，带上标 的导数在自然坐标中心计算。
完整的单元 B-bar 矩阵为
修正后的应变由
给出。
6. 从应变角度理解这一步修正
上面的矩阵表达，可以转化为更直观的应变表达。
令
则
也就是说，B-bar 方法在两个面内法向应变上，分别加上同一个修正量：
由此可以直接得到四个结论。
6.1 面内体积应变变为中心值
将前两式相加：
因此，各积分点采用相同的面内体积应变。
6.2 面内法向应变差保持不变
将前两式相减：
两个方向增加了相同的修正量，所以应变差不会改变。
6.3 工程剪应变保持不变
由于矩阵第四行不作修正，
6.4 平面应变条件保持不变
由于矩阵第三行始终为零，
将上述结果汇总如下：
B-bar 并不是将全部应变都取为中心值，而是只替换面内体积部分。
7. 单元内力与切线刚度
在本文的小变形框架下，采用修正后的应变进行材料响应计算，并使用同一个 构造单元内力与材料切线刚度。
7.1 应力与材料切线
采用应力排列
材料切线矩阵应与本文的应力、应变排列一致：
对于增量形式的材料响应，使用
这里， 表示高斯积分点编号，与节点编号 区分。
7.2 单元内力
采用 高斯积分，单元内力为
7.3 单元切线刚度
在几何线性条件下，相应的单元切线刚度为
其中：
• ：等参坐标映射的雅可比矩阵； • ：单元厚度； • ：二维高斯积分权重； • ：第 个积分点的材料切线矩阵。
有限变形时不能直接照搬
上式是小变形框架下的表达。有限变形问题还涉及变形梯度修正、应力与应变度量，以及相应的几何线性化等内容，不能仅将标准矩阵替换为 就认为得到了完整切线。
8. 为什么这不是 CPE4R 的单点积分？
一个常见误解是：
既然体积部分取了中心值，整个单元是不是就只有一个积分点？
不是。
在上述 B-bar 表达中：
1. 面内体积部分使用单元平均值； 2. 面内偏量部分仍随积分点位置变化； 3. 应力和材料切线仍在四个高斯点计算； 4. 内力与刚度仍按四点积分求和。
因此，CPE4 的选择性减缩积分，与 CPE4R 的整体减缩积分不是同一回事。[2,3]
| | ||
记住这个区别
体积部分采用中心值，不等于整个单元采用单点积分。
9. B-bar 方法能解决所有锁死问题吗？
不能。
B-bar 方法主要针对的是体积锁死，不应将其理解为所有锁死问题的通用解决方案。
对于一阶四边形单元，弯曲问题中的剪切锁死仍需另外关注。Abaqus 文档也明确提醒，一阶完全积分单元在容易发生剪切锁死的弯曲问题中存在局限。[1,3]
因此：
• 缓解体积锁死，不等于消除剪切锁死； • 采用 B-bar，不等于任意粗网格都能得到可靠结果； • 采用 B-bar，也不等于采用了混合压力—位移单元。
10. 小结
在本文的小应变框架下，CPE4 单元的 B-bar 修正可以概括为
对于本文讨论的恒定厚度 Q4 平面单元，又可写为
其关键特征是：
• 修正系数为 ： 修正只分配到两个面内法向方向； • 体积应变采用平均值： 在本文条件下等于自然坐标中心值； • 面内偏量保持不变： 法向应变差与工程剪应变不受修正； • 平面应变条件保持不变： ； • 仍采用四点积分： 并非将整个单元改为单点积分。
最后，用一句话记住它：
CPE4 的 B-bar 方法，是在保留面内偏量变化的同时，用单元平均体积变形替换积分点处的局部体积变形。
参考资料
[1] Dassault Systèmes. Abaqus 2025 Theory Guide: Solid isoparametric quadrilaterals and hexahedra，尤其是 “Fully integrated first-order isoparametric elements” 一节。
[2] Dassault Systèmes. Abaqus 2025 User Subroutines Guide: UMAT，尤其是 “Deformation Gradient” 一节。
[3] Dassault Systèmes. Abaqus 2025 Benchmarks Guide: Performance of continuum and shell elements for linear a nalysis of bending problems。
