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 ， 这样既能很好地描述材料的不可压缩性又可以进行高效的计算。