(完整)abaqus系列教程-09显式非线性动态分析

9 显式非线性动态分析
在前面的章节中，已经考察了显式动态程序的基本内容；在本章中，将对这个问题进行更详细的讨论。

显式动态程序对于求解广泛的、各种各样的非线性固体和结构力学问题是一种非常有效的工具。

它常常对隐式求解器是一个补充，如
ABAQUS/Standard；从用户的观点来看，显式与隐式方法的区别在于：
•显式方法需要很小的时间增量步，它仅依赖于模型的最高固有频率，而与载荷的类型和持续的时间无关。

通常的模拟需要取10,000至1,000,000个增量步，每个增量步的计算成本相对较低。

•隐式方法对时间增量步的大小没有内在的限制；增量的大小通常取决于精度和收敛情况。

典型的隐式模拟所采用的增量步数目要比显式模拟小几个数量级。

然而，由于在每个增量步中必须求解一套全域的方程组，所以对于每一增量步的成本，隐式方法远高于显式方法。

了解两个程序的这些特性，能够帮助你确定哪一种方法是更适合于你的问题。

9.1 ABAQUS/Explicit适用的问题类型
在讨论显式动态程序如何工作之前，有必要了解ABAQUS/Explicit适合于求解哪些类问题。

贯穿这本手册，我们已经提供了贴切的例题，它们一般是应用
ABAQUS/Explicit求解的如下类型问题：
高速动力学（high-speed dynamic）事件
最初发展显式动力学方法是为了分析那些用隐式方法（如
ABAQUS/Standard）分析起来可能极端费时的高速动力学事件。

作为此类模拟的例子，在第10章“材料”中分析了一块钢板在短时爆炸载荷下的响应。

因为迅速施加的巨大载荷，结构的响应变化的非常快。

对于捕获动力响应，精确地跟踪板内的应力波是非常重要的。

由于应力波与系统的最高阶频率相关联，因此为了得到精确解答需要许多小的时间增量。

复杂的接触（contact）问题
应用显式动力学方法建立接触条件的公式要比应用隐式方法容易得多。

结论是ABAQUS/Explicit能够比较容易地分析包括许多独立物体相互作用的复杂接触问题。

ABAQUS/Explicit是特别适合于分析受冲击载荷并随后在结构内部发生复杂相互接触作用的结构的瞬间动态响应问题。

在第12章“接触”中展示的电路板跌落试验就是这类问题的一个例子。

在这个例子中，一块插入在泡沫封装中的电路板从1m的高度跌落到地板上。

这个问题包括封装与地板之间的冲击，以及在电路板和封装之间的接触条件的迅速变化。

复杂的后屈曲（postbuckling）问题
ABAQUS/Explicit能够比较容易地解决不稳定的后屈曲问题。

在此类问题中，随着载荷的施加，结构的刚度会发生剧烈的变化。

在后屈曲响应中常常包括接触相互作用的影响。

高度非线性的准静态（quasi-static）的问题
由于各种原因，ABAQUS/Explicit常常能够有效的解决某些在本质上是静态的问题。

准静态过程模拟问题包括复杂的接触，如锻造、滚压和薄板成型等过程一般地属于这类问题。

薄板成型问题通常包含非常大的膜变形、褶皱和复杂的摩擦接触条件。

块体成型问题的特征有大扭曲、瞬间变形以及与模具之间的相互接触。

在第13章“ABAQUS/Explicit准静态分析”中，将展示一个准静态成型模拟的例子。

材料退化（degradation）和失效（failure）
在隐式分析程序中，材料的退化和失效常常导致严重的收敛困难，但是ABAQUS/Explicit能够很好地模拟这类材料。

混凝土开裂的模型是一个材料退化的例子，其拉伸裂缝导致了材料的刚度成为负值。

金属的延性失效模型是一个材料失效的例子，其材料刚度能够退化并且一直降低到零，在这段时间中，单元从模型中被完全除掉。

这些类型分析的每一个问题都有可能包含温度和热传导的影响。

9.2 动力学显式有限元方法
这一节包括ABAQUS/Explicit 求解器的算法描述，在隐式和显式时间积分之间进行比较，并讨论了显式方法的优越性。

9.2.1 显式时间积分
ABAQUS/Explicit 应用中心差分方法对运动方程进行显示的时间积分，应用一个增量步的动力学条件计算下一个增量步的动力学条件。

在增量步开始时，程序求解动
力学平衡方程，表示为用节点质量矩阵M 乘以节点加速度等于节点的合力（在所施u 加的外力P 与单元内力I 之间的差值）：
Mu
=P -I 在当前增量步开始时（t 时刻），计算加速度为：
1()()|()()|t t -=⋅-u
M P I 由于显式算法总是采用一个对角的、或者集中的质量矩阵，所以求解加速度并不复杂；不必同时求解联立方程。

任何节点的加速度是完全取决于节点质量和作用在节点上的合力，使得节点计算的成本非常低。

对加速度在时间上进行积分采用中心差分方法，在计算速度的变化时假定加速度为常数。

应用这个速度的变化值加上前一个增量步中点的速度来确定当前增量步中点
的速度：
()()()()()2
2
(||)
|||2
t t t t
t t t t t t +∆∆∆+
-
∆+∆=+
u
u
u
速度对时间的积分并加上在增量步开始时的位移以确定增量步结束时的位移：
()()()()2
||||t t t t t t t t +∆+∆∆+
=+∆u u u
这样，在增量步开始时提供了满足动力学平衡条件的加速度。

得到了加速度，在时间上“显式地”前推速度和位移。

所谓“显式”是指在增量步结束时的状态仅依赖于该增量步开始时的位移、速度和加速度。

这种方法精确地积分常值的加速度。

为了
使该方法产生精确的结果，时间增量必须相当小，这样在增量步中加速度几乎为常数。

由于时间增量步必须很小，一个典型的分析需要成千上万个增量步。

幸运的是，因为不必同时求解联立方程组，所以每一个增量步的计算成本很低。

大部分的计算成本消耗在单元的计算上，以此确定作用在节点上的单元内力。

单元的计算包括确定单元应变和应用材料本构关系（单元刚度）确定单元应力，从而进一步地计算内力。

这里给出了显式动力学方法的总结：
1. 节点计算
a. 动力学平衡方程
1()()()()()t t t -=⋅-u
M P I b. 对时间显式积分
()()()(2
2
()
2
t t t t t
t t t t t +∆∆∆+
-
∆+∆=+
u
u
u
()()()()2
t t t t t t t t +∆+∆∆+
=+∆u u u
2. 单元计算
a. 根据应变速率，计算单元应变增量ε d ε
b. 根据本构关系计算应力σ
()()(,)
t t t f d σσε+∆=c. 集成节点内力()
t t +∆I 3. 设置时间 t 为，返回到步骤1。

t t +∆9.2.2 比较隐式和显式时间积分程序
对于隐式和显式时间积分程序，都是以所施加的外力P 、单元内力I 和节点加速度的形式定义平衡：
Mu
=P -I 其中M 是质量矩阵。

两个程序求解节点加速度，并应用同样的单元计算以获得单元内力。

两个程序之间最大的不同在于求解节点加速度的方式上。

在隐式程序中，通过直接求解的方法求解一组线性方程组，与应用显式方法节点计算的相对较低成本比较，求解这组方程组的计算成本要高得多。

在完全Newton 迭代求解方法的基础上，ABAQUS/Standard 使用自动增量步。

在时刻增量步结束时，Newton 方法寻求满足动力学平衡方程，并计算出同一时刻t t ∆+的位移。

由于隐式算法是无条件稳定的，所以时间增量比应用于显式方法的时间增t ∆量相对地大一些。

对于非线性问题，每一个典型的增量步需要经过几次迭代才能获得满足给定容许误差的解答。

每次Newton 迭代都会得到对于位移增量的修正值。

j u ∆j c 每次迭代需要求解的一组瞬时方程为
j
j j j j j u M I P c K --=ˆ对于较大的模型，这是一个昂贵的计算过程。

有效刚度矩阵是关于本次迭代的切向j K ˆ刚度矩阵和质量矩阵的线性组合。

直到一些量满足了给定的容许误差才结束迭代，如力残差、位移修正值等。

对于一个光滑的非线性响应，Newton 方法以二次速率收敛，描述如下：
迭代 相对误差 1 1
2 10-2
3 10-
4 . . . .
.
.
然而，如果模型包含高度的非连续过程，如接触和滑动摩擦，则有可能失去二次收敛，而是可能需要大量的迭代过程。

为了满足平衡条件，减小时间增量的值可能是必要的。

在极端情况下，在隐式分析中的求解时间增量值可能与在显式分析中的典型
稳定时间增量值在同一量级上，但是仍然承担着隐式迭代的高昂求解成本。

在某些情况下，应用隐式方法甚至可能不会收敛。

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

对于大型问题，对这些方程求解器的需求优于对单元和材料的计算的需求，对于在ABAQUS/Explicit中的分析这是类似的。

随着问题尺度的增加，对方程求解器的需求迅速地增加，因此在实践中，隐式分析的最大尺度常常取决于给定计算机中的磁盘空间的大小和可用内存的数量，而不是取决于需要的计算时间。

9.2.3 显式时间积分方法的优越性
显式方法特别地适用于求解高速动力学事件，它需要许多小的时间增量来获得高精度的解答。

如果事件持续的时间是非常短，则可能得到高效率的解答。

在显式方法中可以很容易地模拟接触条件和其它一些极度不连续的情况，并且能够一个节点一个节点地求解而不必迭代。

为了平衡在接触时的外力和内力，可以调整节点加速度。

显式方法最显著的特点是没有在隐式方法中所需要的整体切线刚度矩阵。

由于是显式地前推模型的状态，所以不需要迭代和收敛准则。

9.3 自动时间增量和稳定性
稳定性限制了ABAQUS/Explicit求解器所能采用的最大时间步长，这是应用ABAQUS/Explicit进行计算的一个重要因素。

下面一节将描述稳定性限制并讨论在ABAQUS/Explicit中如何确定这个值，还将讨论影响稳定性限制的有关模型设计参数的问题，这些模型参数包括模型的质量、材料和网格剖分。

9.3.1 显式方法的条件稳定性
应用显式方法，基于在增量步开始时刻t的模型状态，通过时间增量t 前推到当前时刻的模型状态。

这个使得状态能够前推并仍能够保持对问题的精确描述的时间是非常短的。

如果时间增量是大于这个最大的时间步长，则此时间增量已经超出了稳定
性限制（stability limite ）。

超过稳定性限制的一个可能后果就是数值不稳定，它可能导
致解答不收敛。

由于一般不可能精确地确定稳定性限制，因而采用保守的估计值。

因为稳定性限制对可靠性和精确性有很大的影响，所以必须一致性和保守地确定这个值。

为了提高计算的效率，ABAQUS/Explicit 选择时间增量，使其尽可能地接近而且又不超过稳定性限制。

9.3.2 稳定性限制的定义
以在系统中的最高频率（ωmax ）的形式定义稳定性限制。

无阻尼的稳定性限制由
下式定义
stable max
2
t ω∆=
而有阻尼的稳定性限制由下面的表达式定义
stable max
2
)
t ξω∆=
-式中，是最高频率模态的临界阻尼部分。

（回顾临界阻尼，它定义了在自由的和有ξ阻尼的振动关系中在有振荡运动与无振荡运动之间的限制。

为了控制高频振荡，ABAQUS/Explicit 总是以体积粘性的形式引入一个小量的阻尼。

）这也许与工程上的
直觉相反，阻尼通常是减小稳定性限制的。

在系统中的实际最高频率是基于一组复杂的相互作用因素，而且是不大可能计算出确切的值。

代替的办法是应用一个有效的和保守的简单估算。

我们不是考虑模型整体，而是估算在模型中每个个体单元的最高频率，它总是与膨胀模态有关。

可以证明，由逐个单元为基础确定的最高单元频率总是高于有限元组合模型的最高频率。

基于逐个单元的估算，稳定极限可以用单元长度和材料波速重新定义：
e
L d c
stable
e d
L t c ∆=因为没有明确如何确定单元的长度，对于大多数单元类型，例如一个扭曲的四边形单元，上述方程只是关于实际的逐个单元稳定极限的估算。

作为近似值，可以采用最短的单元尺寸，但是估算的结果并不一定是保守的。

单元长度越短，稳定极限越小。

波
速是材料的一个特性。

对于泊松比为零的线弹性材料
d c =
其中，E 是杨氏模量，ρ是密度。

材料的刚度越大，波速越高，导致越小的稳定极限；密度越高，波速越低，导致越大的稳定极限。

这种简单的稳定极限定义提供了某些直觉上的理解。

稳定极限是当膨胀波通过由单元特征长度定义的距离时所需要的时间。

如果我们知道最小的单元尺寸和材料的波速，我们就能够估算稳定极限。

例如，如果最小单元尺寸是5 mm ，和膨胀波速是
5000 m/s ，稳定的时间增量就是在1×10-6 s 的量级上。

9.3.3 在ABAQUS/Explicit 中的完全自动时间增量与固定时间增量
在分析的过程中，ABAQUS/Explicit 应用在前一节讨论过的那些方程调整时间增量的值，使得基于模型的当前状态的稳定极限永不越界。

时间增量是自动的，并不需用户干涉，甚至不需要建议初始的时间增量。

稳定极限是从数值模型得来的一个数学概念。

因为有限元程序包含了所有的相关细节，所以能够确定出一个有效的和保守的稳定极限。

然而，ABAQUS/Explicit 容许用户不必顾及自动时间增量。

在第9.7节“摘
要”中简要地讨论了人工时间增量控制。

在显式分析中所采用的时间增量必须小于中心差分算子的稳定极限。

如果未能使用足够小的时间增量则会导致不稳定的解答。

当解答成为不稳定时，求解变量（如位移）的时间历史响应一般会随着振幅的增加而振荡。

总体的能量平衡也将发生显著的变化。

如果模型只包含一种材料，则初始时间增量是直接与网格中的最小单元尺寸成正比。

如果网格中包含了均匀尺寸的单元但是却包含有多种材料，那么具有最大波速
的单元将决定初始的时间增量。

在具有大变形和/或非线性材料响应的非线性问题中，模型的最高频率将连续地变化，并因而导致稳定极限的变化。

对于时间增量的控制，ABAQUS/Explicit有两种方案：完全的自动时间增量（程序中考虑了稳定极限的变化）和固定的时间增量。

应用两种估算方法确定稳定极限：逐个单元法和整体法。

在分析开始时总是使用逐个单元估算法，并在一定的条件下转变为整体估算法。

逐个单元估算法是保守的；与基于整体模型最高频率的真正的稳定极限相比较，它将给出一个更小的稳定时间增量。

一般说来，约束（如边界条件）和动力学接触具有压缩特征值响应谱的效果，而逐个单元估算法没有考虑这种效果。

另一方面，整体估算法应用当前的膨胀波波速确定整个模型的最高阶频率。

这种算法为了得到最高频率将连续地更新估算值。

整体估算法一般地将允许时间增量超出逐个单元估算法得到的值。

在ABAQUS/Explicit中也提供了固定时间增量算法。

确定固定时间增量的值或者采用在分析步中初始的逐个单元稳定性估算法，或者采用由用户直接指定的时间增量。

当要求更精确地表达问题的高阶模态响应时，固定时间增量算法可能是更有用的。

在这种情况下，可能采用比逐个单元估算法更小的时间增量值。

当在分析步中应用了固定时间增量，ABAQUS/Explicit将不再检查计算的响应是否稳定。

通过仔细地检查能量历史和其他的响应变量，用户应当确保得到了有效的响应。

9.3.4 质量缩放以控制时间增量
由于质量密度影响稳定极限，在某些情况下，缩放质量密度能够潜在地提高分析的效率。

例如，许多模型需要复杂的离散，因此有些区域常常包含着控制稳定极限的非常小或者形状极差的单元。

这些控制单元常常数量很少并且可能只存在于局部区域。

通过仅增加这些控制单元的质量，就可以显著地增加稳定极限，而对模型的整体动力学行为的影响是可以忽略的。

在ABAQUS/Explicit中的自动质量缩放功能，可以阻止这些有缺陷的单元不影响稳定极限。

质量缩放可以采用两种基本方法：直接地定义一个缩放因子或者给那些质量需要缩放的单元逐个地定义所需要的稳定时间增量。

这两种方法都容许对稳定极限附加用户控制，详细介绍请参考ABAQUS分析用户手册第7.15.1节“Mass scaling”。

然而，当采用质量缩放时也要小心，因为模型质量的显著变化可能会改变问题的物理模型。

9.3.5 材料对稳定极限的影响
材料模型通过它对膨胀制波波速的限制作用来影响稳定极限。

在线性材料中，波速是常数；所以，在分析过程中稳定极限的唯一变化来自于最小单元尺寸的变化。

在非线性材料中，例如产生塑性的金属材料，当材料屈服和材料的刚度变化时波速发生变化。

在整个分析过程中，ABAQUS/Explicit监督在模型中材料的有效波速，并应用在每个单元中的当前材料状态估算稳定性。

在屈服之后刚度下降，减小了波速并因而相应地增加了稳定极限。

9.3.6 网格对稳定极限的影响
因为稳定极限大致与最短的单元尺寸成比例，所以应该优先使单元的尺寸尽可能的大。

遗憾的是，对于精确的分析采用一个细划的网格常常是必要的。

为了在满足网格精度水平要求的前提下，尽可能地获得最高的稳定极限，最好的方法是采用一个尽可能均匀的网格。

由于稳定极限是基于在模型中最小的单元尺寸，甚至一个单独的微小单元或者形状极差的单元都能够迅速地降低稳定极限。

为了便于发现问题，ABAQUS/Explicit在状态文件（.sta）中提供了网格中具有最低稳定极限的10个单元的清单。

如果在模型中包含了一些稳定极限比网格中其它单元小得多的单元，将模型网格重新划分使其更加均匀可能是值得的。

9.3.7 数值不稳定性
在大多数情况下，ABAQUS/Explicit对于大多数单元保持了稳定。

但是，如果定义了弹簧和减振器单元，它们在分析过程中有可能成为不稳定。

因此，能够在你的分析过程中识别是否发生了数值不稳定性是非常有用的。

如果确实发生了数值不稳定，典型的情况是结果是无界的，没有物理意义的，而且常常解是振荡的。

9.4 例题：在棒中的应力波传播
本例题展示了在前面第2章“ABAQUS基础”中所描述过的显式动态分析的一些
基本思想。

它也描述了稳定极限，以及在求解时网格细划和材料的影响。

棒的尺寸如图9-1所示。

主视图侧视图
图9-1 在棒中波传播的问题描述
为了使问题成为一个一维的应变问题，所有四个侧面均由滚轴支撑；这样，三维模型模拟了一个一维问题。

材料为钢材，其性质如图9-1所示。

棒的自由端承受一个量级为1.0×105 Pa的爆炸载荷，如图9-2所示，爆炸载荷的持续时间为3.88×10-5 s 。

幅
振
时间
图9-2 爆炸载荷的幅值－时间曲线
9-11
9.4.1 前处理——用ABAQUS/CAE创建模型
在这一节中，我们将讨论如何应用ABAQUS/CAE创建这个模拟所用的模型。

在
本手册的在线文档第A.7节“Stress wave propagation in a bar”提供了输入文件。

当通过ABAQUS/CAE运行这个输入文件时，将创建关于该问题的完整的分析模型。

根据
下面给出的指导如果你遇到困难，或者如果你希望检查你的工作，则可以运行这个输入文件。

在附录A“Example Files”中，给出了如何提取和运行输入文件的指导。

如果你没有进入ABAQUS/CAE或者其它的前处理器，可以人工创建关于这个问
题的输入文件，关于这方面的讨论，见Getting Started with ABAQUS/Standard：Keywords Version，第3.4节“Example：stress wave propagation in a bar”。

定义模型几何
在这个例子中，应用可拉伸实体的基本特征，将创建一个三维的可变形物体。

首先画一个棒的二维轮廓图然后将它拉伸成型。

创建部件：
1．在Create Part对话框中，创建一个部件并命名为Bar，接受三维的变形体和可拉伸实体的基本特征的默认设置，对于模型采用近似的尺寸为0.50。

2．应用在图9-3
中给出的尺寸画棒的横截面。

9-12
图9-3 矩形
可以采用如下的步骤：
a.应用位于画图工具箱右上角的Create Lines: Connected工具创建一个
0.20 m高×0.20 m宽的矩形。

b.当完成绘制轮廓图后，在提示区点击Done。

显示Edit Base Extrusion（编辑基础拉伸）对话框。

为了完成部件定义，
你必须指定横截面拉伸的距离。

c.在对话框中，输入拉伸深度1.0 m。

3．将模型保存到名为Bar.cae的模型数据库文件中。

定义材料和截面性质
创建一个单一线弹性材料，命名为Steel，采用密度7800kg/m3，杨氏模量为207E9 Pa和泊松比0.3。

创建一个均匀的实体截面定义，命名为BarSection，接受Steel作为材料，接受Plane stress/strain thickness为1。

将截面定义BarSection赋予整个部件。

创建装配件
进入Assembly模块，并创建一个部件Bar的实体。

模型按照默认方向放置，整体的3轴位于棒的长度方向。

创建几何集合和面
创建几何集合TOP、BOT、FRONT、BACK、FIX和OUT，如图9-4所示（集合OUT包含楞边，在图9-4中如黑粗线所示）。

创建面命名LOAD，如图9-5所示。

这些区域将用于施加载荷和边界条件，以及定义需要的输出变量。

9-13
9-14
图9-4
集合
图9-5 面
定义分析步
创建一个单一的动态、显式分析步，命名为BlastLoad 。

键入Apply pressure load pulse 作为分析步的描述，并设置Time period 为2.0E-4 s 。

在Edit
Step 对话框中，点击Other 页。

为了保持应力波尽可能的尖锐，将
Quadrastic bulk viscosity parameter （二次体积粘性参数，将在第9.5.1节“体积粘性”中讨论）设置为0。

设置输出要求
编辑默认的场变量输出要求，这样在分析步BlastLoad 中，将预先选择的场变
量数据以四个相等的空间间隔写入输出数据库。

删除已存在的默认的历史变量输出请求，而创建一个新的历史变量输出请求的集合。

在Create History Output（创建历史变量输出）对话框中，接受默认的名称H-Output-1和选择的分析步BlastLoad，点击Continue。

点击在
Domain（范围）选项框旁边的箭头，选择Set name（集合名称），然后选择OUT。

在Output Varables（输出变量）列表中，点击在Stresses左边的三角形，点击在S, Stress components and invariants（应力分量与不变量）左边的三角形，并选中S33变量，它是在棒的轴向的应力分量，指定在每1.0E-6s保存一次输出。

定义边界条件
创建一个边界条件，命名为Fix right end，并在所有三个方向上约束棒的右端面（几何集合FIX）（见图9-1）。

创建其他的边界条件，在这些面的法线方向约束顶面、底面、前面和后面（集合FRONT和BACK为1方向，集合TOP和BOT为2方向）。

定义载荷历史
爆炸载荷将以它的最大值瞬时地施加并保持为常数，持续时间为3.88×10-5 s。

然后载荷突然全部去除并保持为零值。

创建一个幅值定义，命名为Blast，采用在图9-6中所示的数据。

在本问题中在任意给定时刻的压力载荷值是指压力载荷的给定量级乘以由幅值曲线插值的值。

创建压力载荷，命名为Blast load，并选择BlastLoad作为载荷施加的分析步。

将载荷施加在LOAD面上。

选择Uniform（均匀）分布，指定值为
1.0E5 Pa作为载荷量级，并选择幅值为Blast。

9-15。