1.5 显示有限元理论简介

结构受冲击载荷作用下的能量吸收性能研究，在数值模拟方面通常采用有限元方法进行求解。结构动力学中的非线性问题常常包括三个方面：几何非线性、材料非线性和边界非线性（或称为状态非线性）。非线性问题与线性问题的求解方法是不一样的，在线性问题中，通过施加全部载荷即可求解得到问题的结果，但是对于非线性问题来说，一般不能一次施加所有的载荷，而要以增量方式施加给定的载荷求解，逐步获得计算结果。结构动力学问题的数值求解方法可以分为振型叠加法和直接积分法两大类。因为振型叠加法只能用来求解线性问题，所以非线性动力学问题采用直接积分法进行求解［212］。直接积分法又分为显式积分法（如中心差分法和精细积分法［213］）和隐式积分法（如Newmark法、Houbolt法和Wilson法等）。因为显式积分方法能够很好地捕捉动态冲击过程中应力波的传播，并且对计算规模不敏感，所以针对结构瞬时非线性响应问题，尤其是较大规模结构动态分析问题的求解，经常采用显式积分方法。LS-DYNA是一款世界著名的通用显式有限元分析程序，是现在很多显式有限元分析程序的基础代码。LS-DYNA能够模拟各种复杂问题，例如高速碰撞、爆炸、金属冲压成型等，其数值分析结果的可靠性已经通过无数的试验得到了验证。本书有关蜂窝结构面内、面外冲击及耐撞性设计都是基于LS-DYNA进行的，下面将对其基本理论进行简要介绍。

1.5.1 弹塑性动力学基本方程

LS-DYNA程序主要运用拉格朗日的增量格式来描述物质的运动过程，其动力学基本控制表达式如下［212］：

平衡方程：

在边界∂b1上满足力边界条件：

在边界∂b2上满足位移边界条件：

在接触截面∂b3上满足边界条件：

式中，σij是柯西应力；ρ是材料的密度；fi是单位体积力；üi为加速度；nij是边界上的单位外法向矢量。

以上平衡方程及边界条件也可以写成如下的积分等效形式：

其中，δui在边界∂b2上满足位移的边界条件。对式（1.5.5）在整个区域内进行积分，根据散度定理，可得：

考虑到：

将式（1.5.6）和式（1.5.7）代入式（1.5.5），可以得到上述微分方程的边界条件等价的变分列式：

对上式进行有限元离散，同时引入材料的本构模型，即可获得动力学方程的矩阵形式：

其中，M是质量矩阵；C是阻尼矩阵；P是外载荷向量；F是节点力向量；是节点速度向量；ü是节点加速度向量。

1.5.2 显示积分算法

采用显式积分法求解碰撞问题的动力学方程时，计算是在单元一级上进行的，并不需要组装总体刚度矩阵和求解整体平衡方程组，所以大部分计算机处理时间被用来计算单元的节点力。采用单点高斯积分既可以大大减少计算量并节省数据存储量，又能避免全积分单元大变形分析过于刚硬的问题。然而单点积分极易引起物理上无法实现但在数学上存在的零能模式，也可以称为沙漏模式（hourglass modes）。沙漏现象主要是由于单元刚度矩阵计算中积分点的不足而导致的单元刚度阵秩不足。对于动力学分析过程中沙漏模式的判断，直观上是观察结构变形是否呈现出锯齿形网格，数值上则一般运用沙漏能与结构总内能的比值。因为沙漏会导致计算结果失真甚至发散，所以有必要对沙漏模式进行抑制。只有沙漏能与结构能量的比值在很小的范围内，计算结果才是有效的。一般这个比值应该在10%以内。沙漏变形能够通过附加很小的刚度或者黏性阻尼来抑制，对于高速冲击问题，则一般选用增加人工黏度的方法，而低速冲击问题中一般选用附加刚度法。引入了抑制沙漏的附加力H后，运动方程可改写为：

在计算力学中，中心差分法是最常用的显示积分算法［213-215］。时间增量可定义为：

中心差分法不需要对刚度矩阵进行分解，在用其求解有限元系统动力学方程时，也不需要组装系统的总体刚度矩阵，所以极大地减少了存储量。

中心差分法最大的问题在于它是条件稳定的，临界积分时间步长一般都很小，所以一般只适合用来求解短时间瞬态响应问题。对于率无关材料且采用常应变单元离散的结构，无阻尼系统的稳定积分时间步长为［214］：

其中，ωmax是系统最大固有频率；α是Courant数，显然，对于稳定的积分时间步长，有0＜α＜1。在LS-DYNA中，低速冲击下α的默认值为0.9，而高速冲击下α的默认值为0.67［211］。

在有限差分法中，式（1.5.12）被称为Courant条件，由Courant、Friedrichs和Lewy于1928年首先发现。由式（1.5.12）可知，临界时间步长取决于最小单元尺寸，所以采用中心差分法进行有限元动力分析求解时，单元划分应尽可能尺寸一致，避免存在很小尺寸的单元，否则将会大大减少临界时间步长，大幅度增加计算量。对于有阻尼系统，中心差分法的临界时间步长为：

其中，ξl是阻尼比。当ξl＝0时，对应于无阻尼的情况，式（1.5.13）右端圆括号中的项等于1；当ξl＞0时，对应于有阻尼情况，式（1.5.13）右端圆括号中的项小于1。因此，有阻尼系统的临界时间步长低于同等无阻尼系统，阻尼作用会降低临界时间步长。

1.5.3 接触－碰撞界面算法

整个碰撞压缩过程中的接触状态很复杂，接触状态的改变也非常剧烈，这样的接触问题属于一种典型的状态非线性问题。接触算法的性能好坏是衡量一个非线性有限元分析软件的重要指标。LS-DYNA［216］提供多达数十种接触类型来处理不同界面之间的接触情况，丰富且高效率的接触类型正是该软件的重要功能和特色。运用这些接触类型可以充分处理各物体之间、可变形体自身各部分之间的接触及多胞结构内部的自接触等。LS-DYNA程序的接触类型主要采用了节点约束法、对称罚函数法及分配参数法这3种不同的算法来处理这些接触-碰撞界面［215］。其中，对称罚函数法是目前最主要的一种算法，另外两种算法仅适用于某些特殊类型的接触。

对称罚函数法自1982年被引入LS-DYNA程序后慢慢成为该程序处理接触-碰撞界面的主要算法。对称罚函数法通过惩罚接触界面之间的穿透深度而得名。一旦程序检测到接触对中一个界面上有节点穿透到另一个界面后，便会在该节点和被穿透界面之间引入一个足够大的力，从而将节点拉回到被穿透的界面上。接触对的两个接触面的地位是对等的，因此叫对称罚函数法。因为对称罚函数法不需要进行烦琐的碰撞和释放条件检测及判断，所以程序实现起来简单，并具有对称性、动量守恒准确性、无噪声、很少引起沙漏等优点。

碰撞界面间的摩擦是一个重要的研究问题，摩擦使得接触问题变得更复杂。在LS-DYNA中采用库仑（Coulomb）摩擦力公式来计算接触界面间的摩擦力。动摩擦系数（FD）应该不大于静摩擦系数（FS），如果两个系数取值不相等，则需要同时指定非零指数衰减系数DC来平滑过渡，此时摩擦系数的计算如下式所示：

其中，μc是摩擦系数；vr是接触截面的相对运动速度。根据库伦摩擦力公式计算得到的摩擦力有可能带来很高的剪应力，以至于超出接触面材料的剪切强度极限值。因此，需要引入黏性摩擦系数来限制最大摩擦力。通过这种方法限制后的最大摩擦力为Flim＝VC × Acont，Acont是接触对真实接触面积；VC一般取一对接触面中强度较低的材料的剪切屈服应力，，其中σ0是材料单轴拉伸屈服应力。