有限元仿真中的两大“锁死”难题：剪切自锁与体积自锁的机理与应对策略!

7月前浏览878

在有限元仿真分析中，我们常常会遇到计算结果与理论预期严重不符的情况，其中剪切自锁（Shear Locking）和体积自锁（Volumetric Locking）是两类典型的数值问题。这两种在仿真分析中的“锁死”现象会导致模型表现出异常的刚硬特性，使变形计算结果远小于实际值，常常造成非线性分析不收敛。本文将解析这两种锁死现象的产生机理、典型特征以及工程实践中常用的解决方案。

—01—

剪切自锁：弯曲变形“冻结”之手

剪切自锁是有限元分析中最常见的数值问题之一，我们将理论上或实际上不应该产生剪切变形的场景中，由于单元原因在仿真中产生了剪切变形，或寄生了剪切能量，最终导致位移结果偏小的这种现象称为剪切自锁。

这种现象的本质是单元位移场插值函数无法准确描述弯曲变形，导致本应不存在的“伪剪切应变”出现，从而吸收了大部分外力功，使结构表现出异常的过高刚度。

从力学原理看，当材料受纯弯曲时，根据材料力学中的平截面假定，变形前后截面应保持平面且始终垂直于中性轴，理论上剪应力应为零。但在有限元离散模型中，完全积分的线性单元每条边仅由两个节点定义，其边的位移模式只支持线性变化，不能呈现单元边界弯曲形式，如图1所示。

(a) 单元正常弯曲变形

(b) 单元剪切自锁变形

图1 单元正常弯曲与剪切自锁

当应用完全积分线性单元来模拟受纯弯曲的小块材料时，单元上表面拉伸，下表面压缩，这些都与纯受弯的小块材料一致，但是在每个积分点处，中间直线被迫“倾斜”，表现为图1所示的竖直虚线和水平虚线之间的夹角不再为90°，这说明了这点的剪应力不为零，即在积分点处产生了伴生剪切，这与纯受弯的情况不符。

这种伪剪切变形消耗了变形能，而不产生所希望的弯曲变形，使得实际弯曲变形减小，即表现为单元“过硬”，计算得到的挠度可能比理论值低30%以上。

理论和实践表明，剪切自锁的产生通常与三个因素密切相关：

① 单元阶次：线性单元比二次单元更容易锁死；

② 受力状态：纯弯或弯主导工况；

③ 网格形态：单元长宽比越大，锁死越明显。

从以上剪切自锁的三个密切相关因素可以看出，剪力自锁主要影响受弯曲载荷的完全积分线性单元的行为。而二次单元(如下图)的边界可以弯曲，故它没有严重的剪力自锁问题。

图2 完全积分的线性单元与二次单元

但是如果二次单元发生扭曲或它的弯曲应力有梯度(复杂受力状态)，那么也将有可能出现某种程度的剪切自锁，这两种情况在实际问题中是可能发生的。

因此，对于完全积分的线性单元，我们只有当确信载荷只会在模型中产生很小的弯曲时，才可以采用完全积分的线性单元。如果对载荷产生的变形类型有所怀疑，可应采用不同类型的单元分析对比。

对于完全积分的二次单元，在复杂应力下也有可能发生自锁，因此如果在模型中应用这类单元，也应仔细检查计算结果。

那么在工程实践中，我们如何缓解剪切自锁呢？通常可通过以下四种方法缓解剪切自锁：

① 采用缩减积分单元

缩减积分的概念如图3所示，上图为完全积分，下图为缩减积分（注意在ANSYS软件中，Full Integration 实际是选择缩减积分，即将求解应变的矩阵分为体积项和偏差项，对体积项采用缩减积分，对偏差项采用完全积分）。此处所讲的缩减积分（Reduced Integration）是对体积项和偏差项均采用缩减积分。但该方法会出现沙漏(变形的零能模式)。

图3 完全积分单元与缩减积分单元

值得注意的是，剪切自锁现象还受到单元长宽比的较大影响。如下表所示，有分析表明：当简支梁网格长宽比达到10:1（200mm×20mm单元）时，即使用缩减积分，跨中位移误差仍高达30%；而当采用近正方形单元（20mm×33.33mm）时，结果与理论解（1.042mm）基本吻合。

② 采用增强应变单元（非协调单元）

该方式适用于四边形或六面体低阶单元，该方法是在单元模型中增加附加自由度，以克服剪切自锁。单元质量越好，计算精度越高。计算原理如图4所示，在单元中附加了自由度相当于提高了形函数的阶数；同时允许网格中产生缝隙和重叠，故也称为非协调模式。

图4 增强应变单元

③ 采用高阶单元

如前所述，剪切自锁主要是因为线性单元的边不能弯曲，因此采用二次或更高阶单元（如PLANE183），可大为缓解剪切自锁现象，但这会造成计算成本增加。

④ 网格优化

控制长宽比（建议<3）或在弯曲方向布置足够的单元层数（一般不小于3层）。

—02—

体积自锁：不可压缩材料的“刚性枷锁”

与剪切自锁不同，体积自锁主要发生在处理近似不可压缩材料（如橡胶、生物软组织）的仿真中，当泊松比υ接近0.5时，完全积分单元会出现体积应变被过度约束的现象。从本构关系看，体积模量K的计算公式为：

当υ→0.5时，分母趋近零，K→∞，材料表现出无限大的体积刚度，导致单元无法发生合理的体积变形。在超弹性材料的大变形分析中，这种锁死会使静水压应力呈现不真实的“棋盘式”振荡分布，甚至造成求解中止。

体积自锁的产生需要两个条件：

① 使用完全积分单元；

② 材料泊松比接近0.5。

值得注意的是，平面应力问题不会出现体积自锁，因为平面外应变可补偿体积变化。对于非平面应力问题，如何缓解体积自锁呢？通常可采取如下措施:

①采用混合U-P单元(杂交单元)

混合U-P（位移-压力）公式是解决体积自锁最有效的方法。其核心思想是将静水压应力(图5)作为独立自由度引入单元，与位移场分开求解。

静水压力是一个点处应力张量的均匀部分，也可以理解为所有主应力的平均值，其数学表达为 $ρ = - \frac{1}{3} (σ_{1} + σ_{2} + σ_{3})$ 。其中 $σ_{1}$ 、 $σ_{2}$ 、 $σ_{3}$ 是主应力。它是一种等向性应力状态，作用于材料的各个方向，主要引起体积变化。

图5 静水压力下的单元

在ANSYS中，采用U-P杂交单元可通过插入命令流实现。设置KEYOPT(6)=1启用该技术，同时建议配合体积相容公差（VTOL，默认1e-5）调整，如下:

例如对SOLID185单元，可插入command添加如下命令流：

① 在橡胶材料下启用U-P杂交单元:

KEYOPT,MATID,6,1

② 在求解栏中插入一条 Command：

/SOLU

SOLC,,,,1E-3

② 此外缓解体积自锁现象的其他措施还包括：

B-Bar方法：对体积应变项单独降阶积分；

增强应变单元技术：适合六面体/四边形单元；

材料调整：υ设为0.495～0.499平衡精度与收敛；

缩减积分：但需警惕沙漏效应。

—03—

结语:精准仿真的“解锁”之道

在实际工程中，应对剪切自锁和体积自锁的“抗锁死”技术需权衡计算精度、稳定性和效率。

需要特别注意的是，应对剪切自锁或体积自锁的这些技术往往会相互影响甚至冲突。例如同时采用缩减积分和增强应变可能导致过度软化。在应对这些问题时建议：

① 先尝试网格细化，厚度方向应至少3层单元；

② 对于橡胶密封件，首选混合U-P单元技术；

③ 关注沙漏能比率，尽量控制在<5%，不应超过10%。

剪切自锁和体积自锁作为有限元分析中的典型数值问题，其本质都是单元形函数与实际变形模式不匹配导致的过度约束现象。理解其背后的力学原理（如平截面假定、体积模量发散）和数值机理（积分方案、形函数阶次），才能合理选择单元技术。现代商业软件（如ANSYS、ABAQUS）虽然提供了多种高级单元选项，但作为工程师仍需掌握“问题诊断—方法匹配—结果验证”的完整闭环，特别是在处理非线性材料或复杂接触时，建议通过小规模试算验证单元技术的适用性。只有正确“解锁”这些数值陷阱，才能使仿真真正成为设计优化的可靠指南。

来源：薛定谔的Cube

你真的了解阻尼吗？一文了解不同阻尼的物理机制及软件实现方法！

前段时间，笔者在《模态分析理论干货，从动力学方程解构模态分析的核心！》一文中我们对模态分析中的阻尼做了一定介绍。为进一步了解阻尼作用，本文对此再做一些的探讨。01 结构中的阻尼我们知道，如果一个无阻尼结构做自由振动，则它的振幅会是一个常数。然而在现实世界中，结构中蕴含的能量会被结构运动耗散，结构振动的幅度将逐渐减小，直至停止。这种能量的耗散我们称之为阻尼。通常我们假定阻尼是粘性阻尼，即指振动系统中运动受力大小与运动速度成正比且方向相反的阻力所引起的能量损耗。这是一种基于流体粘性特性的阻尼现象。当粘性阻尼包含在结构运动方程中时，动力学方程可写为:其中：[C]——为结构的阻尼矩阵。结构能量耗散来自于许多因素。其中包括结构连接的摩擦(摩擦生热和做功)、局部材料的迟滞效应，即材料在受力变形过程中，应力与应变之间存在相位差，导致在加载和卸载过程中形成滞回环，从而产生能量损耗的现象。阻尼是一种表征结构能量耗散的方法，它包含了重要的能量吸收耗散功能而又无需模拟具体的效果。我们在进行结构运动分析时，模态分析中特征模态的计算是关于无阻尼系统的，但大多数工程中都包含了某种阻尼，尽管这种阻尼可能很小。对于每个模态，有阻尼和无阻尼的固有频率之间的关系如下: 可以看到，对于临界阻尼比较小(ξ＜0.1)时，有阻尼系统的特征频率非常接近无阻尼系统。当ξ增大时，用无阻尼系统表征有阻尼系统的特征频率就变得不太准确了。而当ξ接近于1时，采用无阻尼系统表征有阻尼系统的特征频率就失效了。如果结构处于临界阻尼状态，在任何扰动后，结构都不会有摆动而是尽可能迅速恢复到它的初始构形，如下图。在《模态分析理论干货，从动力学方程解构模态分析的核心！》一文中，我们知道阻尼包括:粘性阻尼、摩擦(库伦)阻尼、材料阻尼(滞后阻尼)、人工阻尼等。这几种阻尼均有累积叠加的特性。那在仿真软件中是如何实现对阻尼的定义和描述的呢？02 不同阻尼的物理机制和描述不同类型的阻尼有不同的物理机制描述。我们知道，对于材料阻尼(材料内部阻尼)，其物理产生机制是材料内部颗粒摩擦、缺陷变化或分子链运动导致机械能转化为热能。可采用粘性阻尼、滞后阻尼对其进行描述。粘性阻尼的阻尼力与速度成正比(F=C×V)，耗能与频率显著相关。而滞后阻尼的耗能与材料应变幅值相关，与频率无关，可迟滞回线模型描述。关于迟滞回线模型大家可以看看《什么是迟滞回线（Hysteresis loop）?》这片文章，本文不做讲述。对于摩擦阻尼(一种结构阻尼)，其产生的物理机制是结构连接处（螺栓、轴承）的干摩擦（库伦阻尼）或声辐射导致的能量损耗。其描述常采用库伦摩擦模型，即F=μN。仿真中，流体阻尼往往不可忽视，其产生的物理机制是物体在流体中运动时，流体黏性阻力消耗能量。物体在流体中的阻尼数学描述通常需要根据运动情况加以区分。例如，在空气中，低速运动时，产生的阻尼可以采用粘性阻尼描述。而在高速运动时，阻尼力与运动速度不能通过简单的线性关系描述，其阻尼力通常通过如下公式描述: 其中:ρ为流体的密度；υ为物体相对流体的速度；Cd为阻力系数，由实验或仿真确定，受形状、表面粗糙度及雷诺数影响；A为特征面积（m²），通常为迎流投影面积。各因素对流体阻尼的影响如下:因素对阻尼的影响流体密度ρ 密度越大，阻力呈线性增加速度υ 阻力与v^2成正比，高速时非线性显著物体形状流线型设计可降低阻力系数30%-50%表面粗糙度粗糙表面加剧湍流，增大摩擦阻力温度温度升高降低流体粘度可能减小摩擦阻力(但对气体密度影响复杂)03 仿真软件如何处理阻尼？在仿真软件中，软件通过数学模型将各类阻尼效应集成到动力学方程中，主要有以下几种方法。① 线性阻尼模型假定阻尼为粘性阻尼，直接将阻尼效应嵌入运动方程，即: 通过求解带有阻尼项的运动方程，可以得到结构在阻尼作用下的响应。这种方式最为直接，也是我们最熟悉的方式。 ② 复模量模型是描述粘弹性材料在动态载荷下应力-应变关系的核心理论工具，其核心是通过复数形式的模量统一表征材料的弹性储能和粘性能量耗散特性。复模量定义为：其中:M'为储能模量：表征材料储存可恢复弹性形变能量的能力，对应应力-应变同相分量。M''为损耗模量：反映材料因内部摩擦（如分子链运动摩擦）导致的能量耗散，对应应力-应变正交分量。在复模量分析中，量化阻尼效率的关键参数是损耗因子tanδ，tanδ = M''/M'。其值越大，材料耗能能力越强（如橡胶tanδ≈0.5）。复模量模型根据受力形式分为三类：类型表达式适用变形模式应用复数杨氏模量拉伸/压缩橡胶减震垫、高分子薄膜复数剪切模量剪切粘弹性阻尼层、隔振支座复数体积模量体积压缩密封材料、液压系统 ③ 瑞利阻尼模型瑞利阻尼模型假设阻尼矩阵是质量和刚度矩阵的线性组合，即: 上式中α和β是需要定义的常数。需要明确：阻尼正比于质量和刚度的假设没有严格的物理基础。实际上我们对阻尼的分布知之甚少，也就不能保证其他更复杂的模型的正确性。一般来说瑞利阻尼模型对超过临界阻尼10%大阻尼系统不再可靠。在许多仿真软件中我们可以精确定义系统每阶模态的瑞利阻尼。在动力学仿真分析中，恰当地定义阻尼对于获得精确的仿真结果十分重要。但正如本文前面所讲述的那样，阻尼只是近似地模拟了结构吸收耗散能量的特性，各类软件并非试图去模拟引起这种效果的物理机制。因此，我们在进行仿真模拟中确定所需要的阻尼数据是困难的。有时，我们可以通过动态试验获取所需的数据，但许多情况我们不得不通过查阅资料或通过经验获取一些数据，在这种情况下，我们必须谨慎地对待模拟的结果，并通过参数分析研究和评估模拟对阻尼值得敏感性。来源：薛定谔的Cube