显式方法:需要很小的时间增量步,仅依赖于模型的最高固有频率,与载荷的类型和持续时间无关。通常的模拟需要取10000~1000000个增量步,每个增量步的计算成本相对较低。
隐式方法:对时间增量步的大小没有内在的限制,增量的大小通常取决于精度和收敛情况。典型的隐式模拟所采用的增量步数目比显式模拟小几个数量级,单由于在每个增量步中必须求解一套全域的方程组,故对于每一个增量步的成本,要远高于显式方法。
ABAQUS/Explicit适用的问题类型
高速动力事件
发展显式动力学方法的最初原因:为了分析哪些用隐式方法分析可能极端费时的高速动力学事件。例如,钢板在短时爆炸载荷下的响应,因为迅速施加的巨大载荷,结构的响应变化非常快。对于捕获动力响应,精确地跟踪板内的应力波是非常重要的。由于应力波与系统的最高阶频率相关,因此,为了得到精确解答需要许多小的时间增量。
复杂的接触问题
应用显式动力学方法建立接触条件的公式,要比应用隐式方法容易得多。ABAQUS/Explicit能够比较容易地分析包括许多独立物体相互作用的复杂接触问题,特别适合于分析受冲击载荷并随后在结构内部发生复杂相互接触作用的结构瞬间动态响应问题。
复杂的后屈曲问题
ABAQUS/Explicit能够比较容易地解决不稳定的后屈曲(postbuckling)问题。在此类问题中,随着载荷的施加,结构刚度会发生剧烈变化。后屈曲响应中常常包括接触相互作用的影响。
高度非线性的准静态问题
ABAQUS/Explicit常常能够有效地解决某些在本质上是静态的问题。准静态(quasi-static)过程模拟问题,包括复杂的接触,一般属于这类问题。
材料退化和失效
材料退化(degradation):如混凝土开裂,其拉伸裂缝导致了材料的刚度成为负值。
材料失效(failure):如金属的延性失效,其材料刚度能够退化并且一直降低到零,在这段时间中,单元从模型中被完全除掉。
动力学显式有限元方法
显式时间积分
ABAQUS/Explicit应用中心差分法对运动方程进行显式的时间积分,应用一个增量步的动力学条件计算下一个增量步的动力学条件。在增量步开始时,程序求解动力学平衡方程,表示为用节点质量矩阵$M$乘以节点加速度$\ddot{u}$等于节点的合力(所施加的外力$P$与单元内力$I$之间的差值):
$$ M\ddot{u}=P-I $$ 在当前增量步开始时(t 时刻),计算加速度为:
$$\ddot{u}|_{(t)} = (M)^{-1} \cdot (P-I)|_{(t)} $$ 由于显式算法总是采用一个对角的或者集中的质量矩阵,所以求解加速度并不复杂,不必同时求解联立方程。任何节点的加速度完全取决于节点质量和作用在节点上的合力,使得计算成本非常低廉。
对加速度在时间上进行积分,采用中心差分方法,在计算速度的变化时假定加速度为常数。应用这个速度的变化值加上前一个增量步中点的速度来确定当前增量步中点的速度:
$$ \dot{u}|_{(t+\frac{\Delta t}{2})} = \dot{u}|_{(t-\frac{\Delta t}{2})} + \frac{\Delta t|_{(t+\Delta t)}+\Delta t|_{(t)}}{2} \ddot{u}|_{(t)} $$ 速度对时间的积分加上在增量步开始时的位移,以确定增量步结束时的位移:
$$ u|_{(t+\Delta t)} = u|_{(t)} + \Delta t|_{(t+\Delta t)} \dot{u}|_{(t+\frac{\Delta t}{2})}$$ 这样,在增量步开始时提供了满足动力学平衡条件的加速度。得到了加速度,在时间上“显式地”前推速度和位移。所谓“显式”,是指在增量步结束时的状态仅依赖于该增量步开始时的位移、速度和加速度,这种方法精确地积分常值的加速度。为了使该方法产生精确的结果,时间增量必须相当小,这样在增量步中加速度几乎为常数。
显式动力学方法总结
(1)节点计算
动力学平衡方程
$$\ddot{u}_{(t)} = (M)^{-1} (P_{(t)}-I_{(t)}) $$ 对时间显式积分
$$ \dot{u}_{(t+\frac{\Delta t}{2})} = \dot{u}_{(t-\frac{\Delta t}{2})} + \frac{\Delta t_{(t+\Delta t)}+\Delta t_{(t)}}{2} \ddot{u}_{(t)} $$ $$ u_{(t+\Delta t)} = u_{(t)} + \Delta t_{(t+\Delta t)} \dot{u}_{(t+\frac{\Delta t}{2})} $$ (2)单元计算
根据应变速率$\dot{\varepsilon}$,计算单元应变增量$d\varepsilon$
根据本构关系计算应力$\sigma$
$$\sigma_{(t+\Delta t)}=f(\sigma_{(t)},d\varepsilon)$$
集成节点内力$ I_{(t+\Delta t)} $
(3)设置时间 $t$ 为 $t+\Delta t$,回到步骤(1)。
比较隐式和显式积分程序
共同点:均以所施加的外力、单元内力和节点加速度的形式定义平衡:
$$ M\ddot{u}=P-I $$ 两个程序求解节点加速度,并应用同样的单元以获得单元内力。
不同点:最大的不同在于求解节点加速度的方式。
隐式时间积分方法:Newton迭代求解,使用自动增量步。在时刻$t+\Delta t$增量步结束时,Newton方法寻求满足动力学平衡方程,并计算出同一时刻的位移。由于隐式算法是无条件稳定的,所以时间增量$\Delta t$比应用显式算法相对大一些。对于非线性问题,每一个典型的增量步需要经过几次迭代才能获得满足给定容许误差的解答。每次Newton迭代都会得到对于位移增量$\Delta u_{j}$的修正值$c_{j}$,每次迭代需要求解的一组瞬时方程为:
$$\hat{K}_{j}c_{j}=P_{j}-I_{j}-M_{j}\ddot u_{j}$$ 有效刚度矩阵是$\hat{K}_{j}$是关于本次迭代的切向刚度矩阵和质量矩阵的线性组合。对于一个光滑的非线性响应,Newton方法以二次速率收敛,但如果模型包含高度的非连续过程,则有可能失去二次收敛,并需要大量的迭代。对于大型问题,对这些方程求解器的需求由于对单元和材料的计算需求,这同样适用于ABAQUS\Explicit。隐式分析的最大尺度常常取决于给定计算机中磁盘空间的大小和可用内存的数量,而不是取决于需要的计算时间。
显式时间积分方法:特别适用于求解告诉动力学事件,需要许多小的时间增量来获得高精度的解答,如果事件持续的时间非常短,则可能得到高效率的解答。可以很容易地模拟接触条件和其他一些极度不连续的情况,并且能够一个节点一个节点地求解而不必迭代。最显著的特点是没有隐式方法中所需要的整体切向刚度矩阵,因为显式地前推模型的状态,所以不需要迭代和收敛准则。
自动时间增量和稳定性
稳定性限制了ABAQUS\Explicit求解器的最大时间步长,影响稳定性限制的有关模型设计参数包括模型的质量、材料、网格划分。
显式方法的条件稳定性
条件稳定性的来源:应用显式方法,基于在增量步开始时刻的模型状态,通过时间增量前推到当前时刻的模型状态,这个是的状态能够前推并仍能保持对问题的精确描述的时间是非常短的。如果时间增量大于这个最长的时间步长,则此时间增量已经超出了稳定性限制(stability limit)。
超过稳定性限制的后果:数值不稳定,解答不收敛。
解决办法:一般不可能精确地确定稳定性限制,因而采用保守的估计值。为了提高计算效率,ABAQUS\Explicit选择时间增量,使其尽可能接近但又不超过稳定性限制。
稳定性限制的定义
稳定性限制:以在系统中的最高频率$\omega_{max}$的形式定义。
无阻尼的稳定极限:
$$\Delta t_{stable}=\frac{2}{\omega_{max}}$$ 有阻尼的稳定极限:
$$\Delta t_{stable}=\frac{2}{\omega_{max}}(\sqrt{1+\xi^{2}}-\xi)$$ 式中:$\xi$是最高频率模态的临界阻尼部分,临界阻尼定义了在自由的和有阻尼的振动关系中,有振荡运动和无振荡运动之间的限制。为了控制高频振荡,ABAQUS\Explicit总是以体积粘性的形式引入一个小量的阻尼。
在系统中的实际最高频率基于一组复杂的相互作用,不大可能计算出确切的值。代替的方法:应用一个有效的保守的简单估算,不考虑模型整体,估算在模型中每个个体单元的最高频率,它总是与膨胀模态有关。可以证明,以逐个单元为基础确定的最高单元频率总是高于有限元组合模型的最高频率。
基于逐个单元的估算,稳定极限可以用单元长度 $L^{e}$ 和材料波速 $c_{d}$ 重新定义:
$$\Delta t_{stable}=\frac{L^{e}}{c_{d}}$$ 因为没有明确如何确定单元的长度,所以对于大多数单元类型,上述方程只是关于实际的逐个单元稳定极限的估算。作为近似值,可以采用最短的单元尺寸,但是估算的结果不一定保守。波速是材料的一个特性,对于泊松比为零的线弹性材料:
$$c_{d}=\sqrt{\frac{E}{\rho}}$$ 材料的刚度越大,波速越高,导致越小的稳定极限;密度越高,波速越低,导致越大的稳定极限。
这种简单的稳定极限定义提供了某些直觉上的理解,稳定极限是当膨胀波通过由单元特征定义的距离时所需的时间。
ABAQUS/Explicit中的完全自动时间增量和固定时间增量
在显式分析中所采用是时间增量必须小于中心差分算子的稳定极限,如果未能使用足够小的时间增量,则会导致不稳定的解答。当解答成为不稳定时,求解变量(如位移)的历史响应一般会出现振幅不断增加的振荡,总的能量平衡也将发生显著变化。
在具有大变形和/或非线性材料响应的非线性问题中,模型的最高频率将连续变化,导致稳定极限的变化。对于时间增量的控制,ABAQUS/Explicit有两种方案:
完全自动时间增量:程序中考虑了稳定极限的变化
固定的时间增量:确定固定时间增量的值可以采用在分析步中初始的逐个单元稳定性估算法,或者采用由用户直接指定的时间增量。当要求更精确地表达问题的高阶模态响应时,固定时间增量算法更有用。如果在分析步中应用了固定时间增量,ABAQUS/Explicit将不再检查计算的响应是否稳定。
应用两种估算方法确定稳定极限:逐个单元法和整体法,在分析开始时总是使用逐个单元估算法,并在一定条件下转变成整体估算法。
逐个单元估算法:保守,与基于整体模型最高频率的真正稳定极限相比,它将给出一个更小的稳定时间增量。一般来说,约束(如边界条件)和动力学接触具有压缩特征值响应谱的效果,而逐个单元估算法没有考虑这个效果。
整体估算法:利用当前的膨胀波波速确定整个模型的最高阶频率。为了得到最高频率,将连续更新估算值,一般允许时间增量超出逐个单元估算法得到的值。
质量缩放以控制时间增量
由于质量密度影响稳定极限,所以在某些情况下,缩放质量密度能潜在地提高分析的效率。例如,许多模型需要复杂的离散,因此有些区域常常包含控制稳定极限的非常小或者形状极差的单元。这些控制单元很少,而且可能只存在于局部区域,通过仅增加这些控制单元的质量,就能显著地增加稳定极限,而对模型的整体动力学行为的影响是可以忽略的。
在ABAQUS/Explicit中的自动质量缩放功能,可以阻止这些有缺陷的单元稳定极限的影响。质量缩放可采用两种基本方法:(1)直接定义一个缩放因子;(2)给那些质量需要缩放的单元逐个定义所需要的稳定时间增量。
稳定极限的影响因素
材料:通过对膨胀波波速的限制作用来影响稳定极限。
线性材料:波速是常数,稳定极限的唯一变化来自于最小单元尺寸的变化。
非线性材料:当材料屈服和材料刚度变化时波速发生变化。在整个分析过程中,ABAQUS/Explicit监督在模型中材料的有效波速,并应用在每个单元中的当前材料状态估算稳定性。在屈服之后刚度下降,减小了波速,相应地增加了稳定极限。
网格:稳定极限大致与最短的单元尺寸成比例,应采用尽可能均匀的网格。ABAQUS/Explicit在状态文件(.sta)中提供了网格中具有最低稳定极限的10个单元的清单。
数值不稳定
在大多数情况下,ABAQUS/Explicit对于大多数单元保持了稳定性。但是,如果定义了弹簧和减振器单元,在分析过程中有可能成为不稳定。如果发生了数值不稳定,典型的情况是结果变得无界,没有物理意义,而且解常常是振荡的。
小结
- ABAQUS/Explicit应用中心差分方法对时间进行动力学显式积分。
- 显式方法需要许多小的时间增量,但因为不必同时求解联立方程,所以每个增量计算成本很低。
- 随着模型尺寸的增加,显式方法比隐式方法能够节省大量的计算成本。
- 稳定极限是能够用来前推动力学状态并仍保持精度的最大时间增量。
- 整个分析过程中,ABAQUS/Explicit自动地控制时间增量值以保持稳定性。
- 随着材料刚度的增加,稳定极限降低;随着材料密度的增加,稳定极限提高。
- 对于单一材料的网格,稳定极限大致与最小单元的尺寸成比例。