2.3.1 临界时间步长
在前面讲到的显示和隐式积分法中,可以看到时间步长Δt在数值计算时是比较重要的变量,尤其对于显式积分法,需要时间步长尽量小,这样数值计算的结果就越精确。当然时间步长越小需要计算的时间越长,所以选择合适的时间步长以平衡计算精度和计算效率是至关重要的。
那么什么样的时间步长是合适的呢?首先为了稳定计算,数值计算选用的时间步长必须小于临界时间步长(Δt≤Δtc)。那这个临界时间步长又是多少呢?为使数值计算稳定,在显示积分法中(如对于一个不考虑阻尼的系统),这个临界时间步长应该与系统的特征频率有关。
这个公式是从系统平衡方程中得来的。以显式积分法所用的单质点无阻尼弹簧为例,平衡方程为
推导可得
式(2-37)也可以转换成矩阵形式:
式中, 对于这个无阻尼弹簧振子的例子来说,A矩阵即
然后计算这个A矩阵的特征值可以得到
令, A2 =det( A) =1 , 那么特征值为。 如果根号里面为负, 那么特征值是虚数解; 如果根号里面为零, 则有两个相同解; 如果根号里面为正, 则有两个不同的实数解。 当矩阵 A 的谱ρ( A) =max(|λi( A1 , A2 )|)≤1 时就有稳定解,这样稳定解区域在双特征值的模为1 时, 也就是图2-7所示双曲线和斜线 (边界线) 的交点。因此可以得到两组稳定性条件:≤, - 1≤A2 < 1; - 1 < A1 < 1,A2 =1。在这个单质点的例子中, A2 =1, 所以需要-1<A1 <1, 即
图2-7 稳定解区域
对于简单的单质点弹簧振子, 它的固有频率是ω=, 将ω代入式 (2-41) 可推导出
在Radioss中使用所有单元的最小时间步长,以便满足所有单元稳定计算的要求,即
式中,称为系统的临界时间步长, 即Δtc。
稳定计算的临界时间步长既可以使用固有频率描述,也可以使用质量和刚度来描述,在许多参考文献中也用来描述, 这些不同的描述方法实际上是等价的。 比如以一个一维线弹性的连续介质 (如杆单元TRUSS) 为例, 如图2-8a所示。
通过下面的推导可以将转换为的形式。
图2-8 一维结构
a)一维杆单元示例 b)应力波在一维结构中传播
这个质量均分在构成单元的节点上, 每一个节点上的质量为。 那么
注意, 在这个例子中系统的固有频率是ωmax =。 如图2-8b所示, c是应力波扩散的速度, 也可称为声音在固体中传播的速度:
e所以在这个例子里临界时间步长也可以描述为或者。
这几种临界时间步长的表达方式存在一定的区别。 Δtc =是描述离散点的时间步长, 所以用于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 , 这样既能很好地描述材料的不可压缩性又可以进行高效的计算。