![Radioss 基础理论与工程高级应用](https://wfqqreader-1252317822.image.myqcloud.com/cover/906/41309906/b_41309906.jpg)
2.3.1 临界时间步长
在前面讲到的显示和隐式积分法中,可以看到时间步长Δt在数值计算时是比较重要的变量,尤其对于显式积分法,需要时间步长尽量小,这样数值计算的结果就越精确。当然时间步长越小需要计算的时间越长,所以选择合适的时间步长以平衡计算精度和计算效率是至关重要的。
那么什么样的时间步长是合适的呢?首先为了稳定计算,数值计算选用的时间步长必须小于临界时间步长(Δt≤Δtc)。那这个临界时间步长又是多少呢?为使数值计算稳定,在显示积分法中(如对于一个不考虑阻尼的系统),这个临界时间步长应该与系统的特征频率有关。
![](https://epubservercos.yuewen.com/143706/21511157501519806/epubprivate/OEBPS/Images/23_07.jpg?sign=1739284359-Xvx0BkYFj6qU4sY0JdLov43P3Ar6RAYF-0-a593b5a1d13ec6c9588f7237929eaa3e)
这个公式是从系统平衡方程中得来的。以显式积分法所用的单质点无阻尼弹簧为例,平衡方程为
![](https://epubservercos.yuewen.com/143706/21511157501519806/epubprivate/OEBPS/Images/23_08.jpg?sign=1739284359-S8qAjzjwZ8IamVInk9E45eZneaNeT6Nx-0-7f383b9531d42b9a028f6e463e0a0387)
推导可得
![](https://epubservercos.yuewen.com/143706/21511157501519806/epubprivate/OEBPS/Images/24_01.jpg?sign=1739284359-lZVuy1iPkipRB1PkJbYuEfsOlQctX4n5-0-f15b553191e78a7bdb42ff42940c3613)
式(2-37)也可以转换成矩阵形式:
![](https://epubservercos.yuewen.com/143706/21511157501519806/epubprivate/OEBPS/Images/24_02.jpg?sign=1739284359-VkRijgVe9m2OlztVI97tlV9HOcKclHG0-0-8c1f07326c0265d2cd7b0a4947c019c0)
式中, 对于这个无阻尼弹簧振子的例子来说,A矩阵即
![](https://epubservercos.yuewen.com/143706/21511157501519806/epubprivate/OEBPS/Images/24_03.jpg?sign=1739284359-0hQwvXDslhH9kNlzzzkJs8Z70llZ4Ixe-0-2f99dc1ccc99c6844e9d2d7eb829e3db)
然后计算这个A矩阵的特征值可以得到
![](https://epubservercos.yuewen.com/143706/21511157501519806/epubprivate/OEBPS/Images/24_04.jpg?sign=1739284359-200Odyzj2zh6VLB0depBsKUMUN3WjU8U-0-af45f13b6e976cf1c4af72a28c7e4126)
令, 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, 即
![](https://epubservercos.yuewen.com/143706/21511157501519806/epubprivate/OEBPS/Images/24_09.jpg?sign=1739284359-QrGrJ550QjMFYUSaZB69ZjSrVv2zNwaX-0-de56549347a3d57d2ca432f96d9396d0)
图2-7 稳定解区域
![](https://epubservercos.yuewen.com/143706/21511157501519806/epubprivate/OEBPS/Images/24_10.jpg?sign=1739284359-WRaH4ZVMBgB8TIcLg1p9hQQwsZIJ99ai-0-e35e72db86bded3c40f282dda2ac38e6)
对于简单的单质点弹簧振子, 它的固有频率是ω=, 将ω代入式 (2-41) 可推导出
![](https://epubservercos.yuewen.com/143706/21511157501519806/epubprivate/OEBPS/Images/24_12.jpg?sign=1739284359-vyikhn6Dx5NdjqRV2RHULGnayxu2S08e-0-96947606bf4c03312fc21ebf62926d08)
在Radioss中使用所有单元的最小时间步长,以便满足所有单元稳定计算的要求,即
![](https://epubservercos.yuewen.com/143706/21511157501519806/epubprivate/OEBPS/Images/24_13.jpg?sign=1739284359-rrHfV3QOT60auqWTkj8NpSUo499lnHIA-0-2ca220cfa8a099c0d03ae101165e2886)
式中,称为系统的临界时间步长, 即Δtc。
稳定计算的临界时间步长既可以使用固有频率描述,也可以使用质量和刚度来描述,在许多参考文献中也用来描述, 这些不同的描述方法实际上是等价的。 比如以一个一维线弹性的连续介质 (如杆单元TRUSS) 为例, 如图2-8a所示。
通过下面的推导可以将转换为
的形式。
![](https://epubservercos.yuewen.com/143706/21511157501519806/epubprivate/OEBPS/Images/25_01.jpg?sign=1739284359-XVuOPRpRMDgLFrCcxOWFcSG9HPYytAM9-0-06deeeba65de54d9f19967d2bc4478dc)
图2-8 一维结构
a)一维杆单元示例 b)应力波在一维结构中传播
![](https://epubservercos.yuewen.com/143706/21511157501519806/epubprivate/OEBPS/Images/25_02.jpg?sign=1739284359-D4clXZ0vi5BEVltArCVids6NLxBxiUft-0-9d816edb2d85c44b7a41327179145a28)
这个质量均分在构成单元的节点上, 每一个节点上的质量为。 那么
![](https://epubservercos.yuewen.com/143706/21511157501519806/epubprivate/OEBPS/Images/25_04.jpg?sign=1739284359-55wBt7jLnbg4cM1nwDduTIFk4lldonIx-0-661ca917a27b684787af81d4df834bd5)
注意, 在这个例子中系统的固有频率是ωmax =。 如图2-8b所示, c是应力波扩散的速度, 也可称为声音在固体中传播的速度:
![](https://epubservercos.yuewen.com/143706/21511157501519806/epubprivate/OEBPS/Images/25_06.jpg?sign=1739284359-US0ERb50yTv73pEngmqSmIEAf57Lg9D5-0-191a23fa2c682123a3421999f506ad43)
e所以在这个例子里临界时间步长也可以描述为或者
。
这几种临界时间步长的表达方式存在一定的区别。 Δtc =是描述离散点的时间步长, 所以用于Radioss中的节点时间步长控制, 如/DT/NODA, 而
是描述连续介质的时间步长,所以用于Radioss中的单元时间步长控制, 如/DT/BRICK、 /DT/SHELL、 /DT/AIRBAG、 /DT/AMS等。 当然对于复杂模型, 如果有接触定义, 那么K就需要包含接触刚度, 可以使用/DT/INTER进行时间步长的控制。 类似的还有/DT/THERM、 /DT/GLOB等不同的时间步长的控制。
![](https://epubservercos.yuewen.com/143706/21511157501519806/epubprivate/OEBPS/Images/25_11.jpg?sign=1739284359-sR0JBXOLybzmTAIHnoaX0V3UKlgzGisF-0-538ee9e5dc7c9c7a8c0be52c012fad8c)
由于需要设置合适的接触刚度,所以材料属性中就需要填写符合实际的密度、刚度。另外,超弹性材料一般被认为是不可压缩的,但是数值计算中不能直接设置泊松比ν=0.5,这会引起时间步长趋向无穷小而导致计算效率极低。 比如实体单元的刚度为K=, 当ν=0. 5 时K将无穷大, 进而导致时间步长无穷小。 所以对于不可压缩材料不建议取 ν =0. 5 , 而是取 ν≤0. 495 , 这样既能很好地描述材料的不可压缩性又可以进行高效的计算。