ABAQUS中冲击动力学问题的求解方法

CAEer吴皓 2023年8月31日浏览：2054 收藏：21

冲击载荷随时间迅速变化。当物体的局部位置受到冲击时，所产生的扰动会逐渐传到未扰动的区域去，这种现象称为应力波的传播。当载荷作用时间短、变化快，且受力物体在加载方向的尺寸又足够大时，这种应力波的传播就显得特别重要。

研究动力学问题最终将简化为求解动力学平衡方程式：节点质量矩阵M乘以节点加速度 ABAQUS中冲击动力学问题的求解方法的图1

等于节点的合力（所施加的外力P与单元内力I之间的差值）：

(2-1)

由于考虑了惯性力的影响，动力学平衡方程中出现了质量矩阵，最后得到的求解方程不是代数方程组，而是常微分方程组。

1. 冲击动力学求解方法

如果加载时间过短或者是动态载荷，需要采用动态分析（dynamic analysis）。复合材料的低速冲击就属于动态分析问题。

动态分析又分为隐式分析和显式分析。在隐式分析中，结构的刚度矩阵需要进行多次生成和求逆，这使得分析求解成本大大增加，而且刚度退化和材料失效常常引起计算收敛问题。在显示分析中，能够避免计算收敛，较好地求解这一问题。

1.1 显式与隐式分析的区别

显式与隐式分析的区别在于：

显式分析需要很小的时间增量步，它仅依赖于模型的最高固有频率，而与载荷的类型和持续的时间无关。通常的模拟需要10000~1000000个增量步，每个增量步的计算成本相对较低。它的求解方法是在时间域中以很小的时间增量步向前推出结果，而无需在每一个增量步求解耦合的方程系统，或者生成总体刚度矩阵。

隐式分析对时间增量步的大小没有内在的限制，增量的大小通常取决于精度和收敛情况。典型的隐式模拟所采用的增量步数目要比显式模拟小几个数量级。然而，由于在每个增量步中必须求解一套全域的方程组，所以对于每一增量步的成本，隐式方法远高于显式方法。

1.2 计算方法选择

复合材料层合板低速冲击损伤涉及到复杂的接触问题、材料刚度随着载荷发生变化的问题、材料的退化（degradation）和失效（failure）导致的严重的收敛问题，这些问题在隐式分析中都无法实现或者求解成本比较昂贵。综上所述，本文选取显式分析方法。

1.3 软件介绍

本文的数值模拟采用商业有限元软件ABAQUS，下面对软件进行简单介绍：

ABAQUS是一套功能强大基于有限元方法的工程模拟软件，可以解决从相对简单的线性分析到复杂的非线性模拟等各种问题。ABAQUS具备十分丰富的单元库，可以模拟任意实际形状，也具有很丰富的材料模型库，可以模拟大多数典型工程材料的性能，包括金属、复合材料、橡胶、钢筋混凝土以及地质材料等。作为一种通用的模拟工具，ABAQUS不仅能解决结构分析（应力/位移）问题，而且能够模拟和研究热传导、电子元器件的热控制（热-电耦合）、声学、土壤力学（渗流-应力耦合分析）和压电分析等领域的问题。

ABAQUS为用户提供了广泛的功能，使用起来十分方便，在许多模拟中（包括高度非线性问题），用户只需提供结构的几何形状、材料性质、边界条件和载荷工况等工程数据。在非线性分析中，ABAQUS能自动选择合适载荷增量和收敛准则。并在分析中不断地调整这些参数值，确保获得精确的解答，用户几乎不必去定义任何参数就能控制问题的数值求解过程。ABAQUS由多个模块构成，包括前后处理模块ABAQUS/CAE、主求解器模块ABAQUS/ Standard和ABAQUS / Explicit等。

ABAQUS/CAE（前后处理）

ABAQUS/CAE是ABAQUS有限元分析的前后处理模块，也是建模，分析和仿真的人机交互平台。该模块根据结构的几何图形生成网格，将材料和截面的属性分配到网格上，并施加载荷和边界条件。该模块可以进一步将生成的模型提交到后处理分析模块中运行，对运行情况进行监测，并对计算结果进行后处理。ABAQUS/CAE的后处理支持ABAQUS分析模块的所有功能，并且对计算结果的描述和解释提供了范围很广的选择，除了云图，等值线和动画显示之外，还可以用列表，曲线等其他常用工具的来完成工程显示。

ABAQUS/Standard（隐式求解器）

ABAQUS/Standard是一个通用分析模块，它能够求解广泛领域的线性和非线性问题，包括静态分析、动态分析，以及复杂的非线性耦合物理场分析等。在每一个求解增量步（increment）中，ABAQUS/Standard隐式地求解方程组。

ABAQUS/Explicit（显式求解器）

使用ABAQUS/Explicit可以进行显式动态分析，它适于求解复杂非线性动力学问题和准静态问题，特别是用于模拟短暂、瞬时的动态事件，如冲击和爆炸问题。此外，它对处理接触条件变化的高度非线性问题也非常有效，例如模拟成型问题，它的求解方法是在时间域中以很小的时间增量步向前推出结果，而无需在每一个增量步求解耦合的方程系统，或者生成总体刚度矩阵。

ABAQUS/Explicit不但支持应力/位移分析，而且还支持完全耦合的瞬态温度/位移分析、声固耦合分析。任意的拉格朗日—欧拉自适应网格功能可以有效地模拟大变形非线性问题。将ABAQUS/Standard和ABAQUS/Explicit结合使用，结合二者的隐式和显式求解技术，可以求解更广泛的实际问题。

综上所述，本文应用显示求解器ABAQUS/Explicit进行数值模拟分析。

2. 动力学显式有限元方法

ABAQUS/Explicit是基于显式算法的有限元程序。显式有限元程序可以分析材料由于力学性能退化而引起材料切线刚度矩阵出现非正定的问题。因此ABAQUS/Explicit程序可以分析纤维增强树脂基复合材料结构受冲击这一高度非线性问题。

使用显式有限元分析时，根据材料的本构模型来确定单元的应力、应变和应变率等。一般情况下，ABAQUS/Explicit求解引擎将每个积分点的应变增量和应变率传递给材料本构，然后再根据相应的材料积分点的信息确定当前状态的应力。对于线弹性材料的计算很简单，当前的总应变与材料刚度矩阵的乘积就是当前状态的应力；而对于非线性材料本构（如塑性材料）的计算就很复杂，本构模型模块将计算得到的应力再传递给求解引擎，由其计算出节点的加速度。然后，根据系统的运动方程求解出节点速度和位移，接着确定下一时间步的应变和应变率。

2.1 显式时间积分

ABAQUS/Explicit应用中心差分方法对运动方程进行显式的时间积分，应用一个增量步的动力学条件计算下一个增量步的动力学条件。在增量步开始时，程序求解动力学平衡方程式：节点质量矩阵M乘以节点加速度 ABAQUS中冲击动力学问题的求解方法的图8 等于节点的合力（所施加的外力P与单元内力I之间的差值），即式(2-1)：

ABAQUS中冲击动力学问题的求解方法的图9 在当前增量步开始时（t时刻），计算加速度：

ABAQUS中冲击动力学问题的求解方法的图10 (2-2)

由于显式算法总是采用一个对角的或者集中的质量矩阵，所以求解加速度并不复杂，不必同时求解联立方程。采用节点的加速度完全取决于节点质量和作用在节点上的合力，使得节点计算的成本非常低。

对加速度在时间上进行积分采用中心差分法，在计算速度的变化时假定加速度为常数。应用这个速度的变化值加上前一个增量步中点的速度来确定当前增量步中点的速度：

ABAQUS中冲击动力学问题的求解方法的图11 (2-3)

速度对时间的积分加上在增量步开始时的位移可以确定增量步结束时的位移：

ABAQUS中冲击动力学问题的求解方法的图12 (2-4)

这样，在增量步开始时提供了满足动力学平衡条件的加速度。得到了加速度，在时间上“显式地”前推速度和位移。所谓“显式”是指在增量步结束时的状态仅依赖于该增量步开始时的位移、速度和加速度。这种方法精确地积分常值的加速度。为了使该方法得到精确的结果，时间增量必须相非常小，这样在增量步中加速度几乎为常数。由于时间增量步必须很小，所以一个典型的分析需要成千上万个增量步。因为不必同时求解联立方程组，所以每一个增量步的计算成本很低。大部分计算成本消耗在单元计算上，以此确定作用在节点上的单元内力。单元的计算包括确定单元应变和应用材料本构关系（单元刚度）确定单元应力，从而进一步计算内力。

下面给出显式动力学方法的总结：

(1) 节点计算

1) 动力学平衡方程

ABAQUS中冲击动力学问题的求解方法的图13 (2-5)

2) 对时间显式积分

(2-6)

ABAQUS中冲击动力学问题的求解方法的图15 (2-7)

(2) 单元计算

1) 根据应变速率 ABAQUS中冲击动力学问题的求解方法的图16 ，计算单元应变增量 ABAQUS中冲击动力学问题的求解方法的图17

2) 根据本构关系计算应力 ABAQUS中冲击动力学问题的求解方法的图18

(2-8)

3) 集成节点内力I ABAQUS中冲击动力学问题的求解方法的图20

(3) 设置时间t为 ABAQUS中冲击动力学问题的求解方法的图21 ，返回到步骤(1)。

2.2 比较隐式和显式时间积分程序

对于隐式和显式积分程序，都是以所施加的外力P、单元内力I和节点加速度形式定义平衡：

(2-9)

其中，M是质量矩阵。两个程序求解节点加速度，并应用同样的单元计算以获得单元内力。两个程序最大的区别在于求解节点加速度上。在隐式程序中，通过直接求解的方法求解一组线性方程组，与显式方法节点计算的成本比较，求解这组方程组的计算成本要高得多。

在完全Newton迭代求解方法的基础上，ABAQUS/Standard使用自动增量步。在时刻 ABAQUS中冲击动力学问题的求解方法的图23

增量步结束时，Newton方法寻求满足动力学平衡方程，并计算出同一时刻的位移。由于隐式算法是无条件稳定的，所以时间增量 ABAQUS中冲击动力学问题的求解方法的图24

比应用显式方法的时间增量相对大些。对于非线性问题，每一个典型的增量步需要经过几次迭代才能获得满足给定容许误差的解答。每次Newton迭代都会得到对于位移增量 ABAQUS中冲击动力学问题的求解方法的图25

的修正值

。每次迭代需要求解的一组瞬时方程为

(2-10)

对于较大的模型，这是一个非常大的计算过程。有效刚度矩阵

是关于本次迭代的切向刚度矩阵和质量矩阵的线性组合。直到一些变量满足了给定的容许误差才结束迭代，如位移修正值等。

在隐式分析中，每一次迭代都需要求解大型的线性方程组，这一过程需要相当数量的计算资源、内存和磁盘空间。

2.3 显式时间积分方法的优越性

显式方法特别适用于求解高速动力学问题，它需要许多很小的时间增量来获得高精度的解答。如果事件持续的时间非常短，则可能得到高效率的解答。

在显式方法中可以很容易的模拟接触条件和其他一些极度不连续的情况，并且能够一个节点、一个节点地求解而不必迭代。为了平衡接触时的内力和外力，可以调整节点加速度。

显式方法最显著的特点是没有在隐式方法中所需要的整体切向刚度矩阵。由于是显式地前推模型的状态，所以不需要迭代和收敛准则。

2.4 显式方法的条件稳定性

稳定性限制了ABAQUS/Explicit求解器所能采用的最大时间步长，这是ABAQUS/Explicit 进行计算的一个重要因素。应用显式方法，基于在增量步开始时刻

的模型状态，通过时间增量 ABAQUS中冲击动力学问题的求解方法的图30

前推到当前时刻的模型状态。使得状态能够前推并仍能保持对问题的精确描述的时间非常短。如果时间增量大于这个最大的时间步长，则时间增量已经超出了稳定性限制（stability limit）。超过稳定性限制的后果是数值不稳定，可能导致计算不收敛。由于一般不可能精确地确定稳定性限制，因而采用保守的估计值。为了提高计算效率，ABAQUS/Explicit选择时间增量，使其尽可能地接近而且又不超过稳定性限制。

在系统中，以最高频率（

）的形式定义稳定性限制。无阻尼的稳定极限由下式定义：

(2-11)有阻尼的稳定极限由下面的表达式定义：

(4-14) (2-12)式中， ABAQUS中冲击动力学问题的求解方法的图34

是最高频率模态的临界阻尼部分。由此可以看出，阻尼通常是减小稳定性限制的。

在系统中的实际最高频率基于一组复杂的相互作用因素，而且不可能计算出确切的值。代替的方法是应用一个有效的、保守的简单估算。不是考虑模型整体，而是估算在模型中每个个体单元的最高频率，它总是与膨胀波的波速有关。可以证明，以逐个单元为基础确定的最高单元频率总是高于有限元组合模型的最高频率。

基于逐个单元的估算，稳定极限可以通过单元长度 ABAQUS中冲击动力学问题的求解方法的图35

和材料的波速

重新定义：

(2-13)

因为没有明确如何确定单元长度，所以对于大多数单元类型，例如一个扭曲的四边形单元，上述方程只是关于实际的逐个单元稳定极限的估算。作为近似值，可以采用最小单元的尺寸，但是估算的结果并不一定是保守的。单元长度越短，稳定极限越小。波速是材料的一个特性，一般定义为：

(2-14)

其中，

是杨氏模型；

是材料密度。由上式可以看出，材料的刚度越大，波速越大，导致越小的稳定极限；密度越大，波速越小，导致较大的稳定极限。

稳定极限是当膨胀波通过由单元特征长度定义的距离时需要的时间。如果知道最小的单元尺寸和材料的波速，就能够估算稳定极限。例如，如果最小单元尺寸是5mm，膨胀波速是5000m/s，那么稳定的时间增量就在 ABAQUS中冲击动力学问题的求解方法的图41

的数量级上。

2.5 自动时间增量

在分析过程中，ABAQUS/Explicit应用上述方程调整时间增量值，使得基于模型的当前状态的稳定极限永不越界。时间增量是自动的，并不需要用户干涉，甚至不需要建议初始的时间增量。因为有限元程序包含了所有的相关细节，所以能够确定出一个有效的、保守的稳定极限。

在显式分析中所采用的时间增量必须小于中心差分算子的稳定极限。如果不能使用足够小的时间增量，则会导致计算不稳定。这时，求解变量（如位移）的时间历史响应一般会出现振幅不断增加的振荡。总体的能量平衡也将发生显著的变化。如果模型只包含一种材料，则初始时间增量与网格中的最小单元的尺寸成正比。如果网格单元尺寸均匀却包含了多种材料，则具有最大波速的单元将决定初始时间增量。

在具有大变形或非线性材料响应的问题中，模型的最高频率将连续变化，并导致稳定极限的变化。对于时间增量的控制，ABAQUS/Explicit有两种方案：完全的自动时间增量（程序中考虑了稳定极限的变化）和固定的时间增量。应用两种估算方法确定稳定极限，分为逐个单元法和整体法。在分析开始时总是使用逐个单元估算法，并在一定条件下转为整体估算法。

逐个单元估算法是保守的，与基于整体模型最高频率的稳定极限相比较，它将给出一个更小的稳定时间增量。而整体估算法应用当前的膨胀波速确定整个模型的最高频率。这种算法得到最高频率将持续地更新估算值[4]。

3. 接触

接触问题在工程中处处可见，复合材料的冲击问题也存在冲击头与层合板之间的接触。当两个物体彼此接触时，垂直于接触面的力作用在两个物体上。如果在接触面之间存在摩擦，可能产生剪力以阻止物体的切向运动。接触模拟的目的是确定表面上发生接触的面积和计算所产生的接触力。

在有限元分析中，接触问题属于非线性问题，但又区别于材料非线性和几何非线性，属于边界条件非线性问题。接触条件是特殊的不连续约束，允许力从模型的一部分传递到另一部分。只有两个表面发生接触时才会产生约束，当分开时，就不存在约束了。分析时，需要判断什么时候表面发生接触并采用响应的约束，什么时候表面分开并解除约束。

复合材料冲击接触分为静态接触和动态接触。前者不考虑接触过程随时间的变化，忽略了惯性效应。实际上，复合材料冲击响应一般由瞬态动力响应或者应力波控制的。本文在计算过程中，采用动态接触。

3.1 ABAQUS接触功能描述

ABAQUS/Explicit提供了两种模拟接触相互作用的算法[29]。通用自动接触（general contact）算法和接触对（contact pair）算法。通用自动接触可以非常简单地定义接触，对于接触表面的类型限制很少；接触对算法对于接触表面的类型有比较严格的限制。通常定义一个接触只需简单地指定所采用的接触算法和将会发生接触作用的表面。在特殊情况下，当默认的接触设置不满足需求时，可以指定接触模拟的其它内容，如考虑摩擦的相互作用力学模型。本文采用通用自动接触算法。

3.2 接触属性

接触分析过程中，两个接触面之间的相对位置关系在整个分析过程中随着接触边界条件以及载荷的变化时刻发生变化。确定了接触面在某一接触段发生接触后，需要考虑两个接触面之间的接触属性。

接触属性包括两部分：接触表面之间的切向作用和法向作用。在切向作用中，两个接触面之间的相对滑动关系、摩擦类型是影响接触的主要因素；在法向作用中两个表面之间的间隙是重点考虑的对象。

两个接触面之间相对滑动时，两表面的相对滑动位移大小不同。主要采用的是有限滑移，指的是两个接触面之间可以有任意的相对滑动。在分析过程中，ABAQUS需要不断地判断从面节点和主面节点的哪一部分发生了接触。

计算过程中还需要考虑接触面之间的约束关系。两个表面之间的距离成为间隙，ABAQUS判断两个表面是否接触到的依据就是判断两个表面之间的间隙是否为零，当两个表面之间的间隙为零时，认为两个表面之间发生了接触。这时，接触面之间就会产生接触力。

文章来源：正脉CAE技术平台