2.3 时间步长

动态问题通常使用显式积分法计算，时间步长在显式分析中尤为重要，它既能影响计算精度又能影响计算效率。下面讲解稳定计算的临界时间步长的定义和如何在有限元模型中控制好时间步长。

2.3.1 临界时间步长

在前面讲到的显示和隐式积分法中，可以看到时间步长Δ t 在数值计算时是比较重要的变量，尤其对于显式积分法，需要时间步长尽量小，这样数值计算的结果就越精确。当然时间步长越小需要计算的时间越长，所以选择合适的时间步长以平衡计算精度和计算效率是至关重要的。

那么什么样的时间步长是合适的呢？首先为了稳定计算，数值计算选用的时间步长必须小于临界时间步长（Δ t ≤Δ t _c ）。那这个临界时间步长又是多少呢？为使数值计算稳定，在显示积分法中（如对于一个不考虑阻尼的系统），这个临界时间步长应该与系统的特征频率有关。

这个公式是从系统平衡方程中得来的。以显式积分法所用的单质点无阻尼弹簧为例，平衡方程为

推导可得

式（2-37）也可以转换成矩阵形式：

式中， 对于这个无阻尼弹簧振子的例子来说， A 矩阵即

然后计算这个 A 矩阵的特征值可以得到

令 ， A ₂ =det（ A ） =1 ， 那么特征值为 。 如果根号里面为负， 那么特征值是虚数解； 如果根号里面为零， 则有两个相同解； 如果根号里面为正， 则有两个不同的实数解。 当矩阵 A 的谱 ρ （ A ） =max（｜ λ _i （ A ₁ ， A ₂ ）｜）≤1 时就有稳定解，这样稳定解区域在双特征值的模为1 时， 也就是图2-7所示双曲线和斜线 （边界线） 的交点。因此可以得到两组稳定性条件： ≤ ， - 1≤ A ₂ < 1； - 1 < A ₁ < 1， A ₂ =1。在这个单质点的例子中， A ₂ =1， 所以需要-1< A ₁ <1， 即

对于简单的单质点弹簧振子， 它的固有频率是 ω = ， 将 ω 代入式 （2-41） 可推导出

在Radioss中使用所有单元的最小时间步长，以便满足所有单元稳定计算的要求，即

式中， 称为系统的临界时间步长， 即Δ t _c 。

稳定计算的临界时间步长既可以使用固有频率描述，也可以使用质量和刚度来描述，在许多参考文献中也用 来描述， 这些不同的描述方法实际上是等价的。 比如以一个一维线弹性的连续介质 （如杆单元TRUSS） 为例， 如图2-8a所示。

通过下面的推导可以将 转换为 的形式。

a）一维杆单元示例 b）应力波在一维结构中传播

这个质量均分在构成单元的节点上， 每一个节点上的质量为 。 那么

注意， 在这个例子中系统的固有频率是 ω _max = 。 如图2-8b所示， c 是应力波扩散的速度， 也可称为声音在固体中传播的速度：

e所以在这个例子里临界时间步长也可以描述为 或者 。

这几种临界时间步长的表达方式存在一定的区别。 Δ t _c = 是描述离散点的时间步长， 所以用于Radioss中的节点时间步长控制， 如/DT/NODA， 而 是描述连续介质的时间步长，所以用于Radioss中的单元时间步长控制， 如/DT/BRICK、 /DT/SHELL、 /DT/AIRBAG、 /DT/AMS等。 当然对于复杂模型， 如果有接触定义， 那么 K 就需要包含接触刚度， 可以使用/DT/INTER进行时间步长的控制。 类似的还有/DT/THERM、 /DT/GLOB等不同的时间步长的控制。

由于需要设置合适的接触刚度，所以材料属性中就需要填写符合实际的密度、刚度。另外，超弹性材料一般被认为是不可压缩的，但是数值计算中不能直接设置泊松比 ν =0.5，这会引起时间步长趋向无穷小而导致计算效率极低。 比如实体单元的刚度为 K = ， 当 ν =0. 5 时 K 将无穷大， 进而导致时间步长无穷小。 所以对于不可压缩材料不建议取 ν =0. 5 ， 而是取 ν ≤0. 495 ， 这样既能很好地描述材料的不可压缩性又可以进行高效的计算。

2.3.2 时间步长控制

Radioss中每一步（cycle）使用的时间步长都是程序内部计算的。为了保持数值计算的稳定性和防止发散，Radioss会将计算出的理论临界时间步长乘以时间步长比例因子Δ T _sca 。计算出的理论临界时间步长可以在Radioss的starter输出文件中打印出来：

如果在engine中使用默认的时间步长控制，比如使用以下设置：

那么在engine计算中的最小时间步长为

那么在engine计算中的第一个循环中就可得到以下结果：

除了选用上面的自由时间步长（即Δ t _min 设置为0），还可以通过一些方法来调整和控制时间步长。

一些复杂的结构网格中存在为数不多的尺寸非常小的单元，那么时间步长就会被这些单元所控制，计算效率大大降低，这时推荐使用时间步长控制的方法来解决，比如使用CST（传统的质量缩放技术）来控制，它的形式如下：

使用这种方法时，如果最小时间步长小于定义的最小时间步长，即Δ T _sca ×min（Δ t _mesh ）≤Δ T _min ，那么时间步长控制就开始起作用。

举例如下。

在starter中的最小时间步长显示如下：

而在engine中定义的时间步长如下：

那么engine中初始的时间步长为

当初始的时间步长小于Δ T _min （7.0E-7）时,Radioss就会通过增加微小的质量来提高临界时间步长 （即Δ t _c ↗= ）。

增加的质量使得Δ T _sca ×min（Δ t _mesh ）≤Δ T _min ,也就是节点的临界时间步长为

此时engine中打印的时间步长依然是7E-07，但是在MAS.ERR中会显示质量的增量。

注意，当时间步长比例因子Δ T _sca 设置为0.67而不是0.9时，就需要增加更多的质量才能满足最小的稳定计算时间步长是7E-07，此时的节点稳定计算时间步长为

另外，时间步长控制中的比例因子Δ T _sca 取值在0~1之间。这个参数在Radioss中默认是0.9。某些分析中甚至推荐用户设置更小的值（比如0.66）。那么什么情况下需要这样的设置呢？这其实也是和材料属性有关的。

图2-9中的Ⅰ图所示应力应变曲线为常见的钢的材料曲线图， 在材料进入塑性阶段后， 有效弹性模量比最初的弹性阶段的有效弹性模量小，而且这个硬化阶段的有效弹性模量随着应变的增加是逐渐变小的， 因此临界时间步长是随着变形的增加而变大的， 此时一般不需要用比例因子减小时间步长来满足计算过程中每个时刻的稳定计算， 使用默认的0. 9就完全可以满足要求； 图Ⅱ中显示的有效弹性模量是逐渐增大的； 图Ⅲ中有效弹性模量在某一区域后急剧增大 （比如泡沫橡胶等材料）， 此时临界时间步长随着构件变形的增大而变小， 因此需要通过比例因子来进一步调小Radioss计算中使用的时间步长， 以保证在整个过程中的计算稳定性。

那么究竟如何设置一个合适的时间步长比例因子呢？首先需要保证不能让增加的虚拟质量太大，一般控制在1%以内对整个模型计算精度的影响可以忽略，这样既保证了计算精度又提高了计算效率；其次不能无限制地让Radioss增加质量，也就是不能在CST中设定太大的最小时间步长而导致增加太多的虚拟质量，从而增加太多的动能，毕竟在非线性动态数值分析中需要保证质量守恒、能量守恒和动量守恒，否则数值计算的结果正确性会受影响；此外，当计算出现“负体积”报错时也有可能是因为设置的强制时间步长太大。更多关于时间步长的信息请参见Radioss理论手册（Radioss Theory Manual）。通常对时间步长比例因子除了默认为0.9外，还有如下建议。

1）使用AMS技术时在 /DT/AMS 中建议设置为Δ T _sca =0.67。

2）对于泡沫材料，则Δ T _sca =0.66。

3）模型如果只有一个单元，则Δ T _sca =0.1。

4）模型如果只有两个单元，则Δ T _sca =0.2。

5）模型如果多于三个单元，则Δ T _sca =0.9。

6）时间步长比例因子一定不能大于1.0。

Radioss中控制时间步长的卡片除了CST，还有STOP、DEL等，即一旦当前的时间步长小于用户给定的最小时间步长计算就终止（STOP），或单元就删除（DEL）。另外，局部增加的质量可以通过在engine中使用/ANIM/MASS卡片在HyperView中以云图形式展现。

2.3.3 单元特征长度

单元类型选择和单元网格大小对时间步长是有影响的。那么对于不同的单元，Radioss中特征长度 l _c 的选取也是不同的，同时不同单元类型临界时间步长的计算方法也是不同的，见表2-2。