FlexPhysics-Gym v1.20 大版本更新！

提到 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 修正后的变化 面内法向应变之和 替换为单元平均值，即本文条件下的中心值 面内法向应变之差 保持不变 工程剪应变 保持不变 面外应变 保持为零 B-bar 并不是将全部应变都取为中心值，而是只替换面内体积部分。7. 单元内力与切线刚度在本文的小变形框架下，采用修正后的应变进行材料响应计算，并使用同一个 构造单元内力与材料切线刚度。7.1 应力与材料切线采用应力排列材料切线矩阵应与本文的应力、应变排列一致：对于增量形式的材料响应，使用这里， 表示高斯积分点编号，与节点编号 区分。7.2 单元内力采用 高斯积分，单元内力为7.3 单元切线刚度在几何线性条件下，相应的单元切线刚度为其中：• ：等参坐标映射的雅可比矩阵；• ：单元厚度；• ：二维高斯积分权重；• ：第 个积分点的材料切线矩阵。有限变形时不能直接照搬上式是小变形框架下的表达。有限变形问题还涉及变形梯度修正、应力与应变度量，以及相应的几何线性化等内容，不能仅将标准矩阵替换为 就认为得到了完整切线。8. 为什么这不是 CPE4R 的单点积分？一个常见误解是：既然体积部分取了中心值，整个单元是不是就只有一个积分点？不是。在上述 B-bar 表达中：1. 面内体积部分使用单元平均值；2. 面内偏量部分仍随积分点位置变化；3. 应力和材料切线仍在四个高斯点计算；4. 内力与刚度仍按四点积分求和。因此，CPE4 的选择性减缩积分，与 CPE4R 的整体减缩积分不是同一回事。[2,3]对比项 CPE4 CPE4R 常规积分点数 ，共 4 个 1 个 主要应变处理 对体积部分进行 B-bar 修正 采用均匀应变形式的减缩积分 是否属于同一种积分形式 否 否 沙漏控制 不采用 CPE4R 式的单点沙漏稳定化 需要相应沙漏控制 记住这个区别体积部分采用中心值，不等于整个单元采用单点积分。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。 来源：大狗子说数值模拟