多体系统动力学(2)
(2012-02-20 16:13:36)
标签:
杂谈 |
分类: MBS |
2.3 多柔体系统动力学建模
2.3.1 柔性体上点的位置向量、速度和加速度
2.3.1.1.1
柔性体系统中的坐标系
柔性体系统中的坐标系如图2.5所示,包括惯性坐标系( )和动坐标系()。前者不随时间而变化,后者是建立在柔性体上,用于描述柔性体的运动。动坐标系可以相对惯性坐标系进行有限的移动和转动。动坐标系在惯性坐标系中的坐标(移动、转动)称为参考坐标。
与刚体不同,柔性体是变形体,体内各点的相对位置时时刻刻都在变化,只靠动坐标系不能准确描述该柔性体在惯性坐标系中的位置,因此,引入弹性坐标来描述柔性体上各点相对动坐标系统的变形。这样柔性体上任一点的运动就是动坐标系的“刚性”运动与弹性变形的合成运动。由于柔体上各点之间有相对运动,所以动坐标系的选择不是采用连体坐标系,而需要采用随着柔性体形变而变化的坐标系,即“浮动坐标系”。
在研究多柔体系统时,合适的坐标系是非常重要的。在确定浮动坐标系时有两点准则:1、便于方程建立求解;2、柔性体刚体运动与变形运动的耦合尽量小。目前常见的浮动坐标系大致有如下5种,局部附着框架、中心惯性主轴框架、蒂斯拉德框架、巴克凯恩斯框架以及刚体模态框架。采用何种需因实际情况而定。
2.2.3.1.2
柔性体上任意点的位置向量、速度和加速度
在分析刚体平面运动的时候,把复杂的刚体平面运动分解为几种简单的运动。在对柔性体的运动,尤其是在小变形的情况下,也可以采用类似的方法。如某柔性体从位置L1运动到位置L2,其间运动可以分解为:刚性移动->刚性转动->变形运动。对于柔性体上任意一点P,其位置向量为:
.3- 1-158)
为P点在惯性坐标系中的向量; 为浮动坐标系原点在惯性坐标系中的向量;根据式2.2- 14,A为方向余弦矩阵;为柔性体未变形时P点在浮动坐标系中的向量; 为相对变形向量, 可以用不同的方法离散化,与讨论平面问题相同,对于点,该单元的变形采用模态坐标来描述,有:
.3- 2-159)
式(2.3-2-159)中, 为点
满足里兹基向量要求的假设变形模态矩阵,为变形的广义坐标。
柔性体上任一点的速度向量及加速度向量可以对式求对时间一阶导数和二阶导数得到:
.3- 3-160)
.3- 4-161)
图2.6 柔性体变形模型 |
1.2.3.2.1
外加载荷
在ADAMS软件中,外加载荷包括单点力与扭矩、分布式载荷以及残余载荷三部分。
(1)单点力与扭矩
施加于柔性体上某一标记点的单点力和扭矩必须投影到系统的广义坐标上才能起作用,力和扭矩以矩阵形式写出,在标记点的局部坐标系下表示为:
,
.3-5-162)
广义力由广义平动力、广义扭矩(以欧拉角表示的广义力)和广义模态力组成,可表示为:
.3-6-163)
平动坐标下的广义力可以通过转换单点力 到全局坐标基 下来获得,即:
.3-7-164)
其中, 为标记点 上的坐标系相对于全局坐标的欧拉角变换矩阵,可以表示为:
,如图2.6所示, 为局部坐标系 点相对于全局坐标系的转换矩阵,即方向余弦阵; 为因节点的小变形引起的标记坐标的方位变化而引入的转换矩阵; 为定义在柔性体上标记点处的坐标系相对于 点坐标系的常值变换矩阵。
作用在柔性体标记点处的合力矩,可用相对于全局坐标的矢量矩阵表达为:
.3-8-165)
.3-9-166)
将物理坐标系下的力矩向基于欧拉角的坐标系转换,并利用式 ,可得广义合力矩为:
.3-10-167)
通过投影单点力和单点力矩到模态坐标上,可得到 点处的广义模态力。
如在标记点 处施加力 和力矩 ,通过转换到全局坐标上,即有: , ,并将
投影到平动模态坐标、将投影到角模态坐标上,可得广义模态
.3-11-168)
其中, 和 对应于节点 处的平动和转动自由度的模态斜方阵。由于模态矩阵只定义在节点,故单点力和单点力矩只能施加于节点处。
(2)分布式载荷
在ADAMS软件上,分布式载荷可以通过MFORCE的方式来创建。通常在FEM软件中有运动方程:
.3-12-169)
其中, 和 为柔性体上有限单元的质量和刚度矩阵,而 和 为物理节点的自由度矢量和载荷矢量。
利用模态矩阵 将式(2.3-12-169)转换到模态坐标 下,有:
.3-13-170)
上式可简化为:
.3-14-171)
其中, 和 分别为广义质量和刚度矩阵, 为模态载荷矢量。
节点力矢量在模态坐标上的投影为:
.3-15-172)
上式中,如果
是时间的函数则求解时计算开销太大,一种替代的方法是假设空间依赖性和时间依赖性可以分开,把载荷看成是一系列依赖于时间的静态载荷的线性组合,即:
.3-16-173)
载荷在模态坐标上的投影计算可在有限元文件MNF过程中完成,而不必在ADAMS仿真中重复进行。如果定义一系列静载荷为载荷矢量,并使其与系统响应显性相关,即表示成的形式,则模态力又可表示为:
.3-17-174)
(3)残余载荷
先前的假设是施加的载荷向模态坐标的投影已实现,即式(2.3-15-172)成立。然而由于模态截断,很多情况下施加的力并不能进行投影,这些力被称为残余载荷,可以将其看为已投影到了被截断的高阶模态坐标上。残余载荷可表示为:
.3-18-175)
与残余载荷相关的是残余矢量,可看成是把残余载荷施加于柔性体时产生的变形。残余矢量可被当作模态振型加入到Craig-Bampton模态基中,增强的模态基完全能够捕捉残余载荷,否则,残余力的丢失不可避免。
2.2.3.2.2
多柔体系统的能量
(1)动能和质量矩阵
考虑节点 变形前后的位置、方向和模态,柔性体的广义坐标可以表示为:
.3-19-176)
速度表达式(2.3-3-160)在系统广义坐标式(2.3-19-176)的时间导数 中表示为:
.3-20-177)
柔性体的动能为:
.3-21-178)
其中, 和 分别为节点 的节点质量和节点惯性张量;
,为点相对于全局坐标基的角速度在局部坐标基中的斜方阵表示。将式(2.3-20)和关系式 代入式(2.3-21-178),得到动能的广义表达式:
.3-22-179)
上式中的质量矩阵 为 维的方阵,表示为:
.3-23-180)
其中下标 分别表示平动、旋转和模态自由度。质量矩阵的六个独立分量分别表示为:
.3-24-181)
其中9个惯性时不变矩阵列表如下:
表2-1 惯性时不变矩阵列表
惯性时不变矩阵 |
模态坐标数 |
维数 |
|
标量 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
式(2.3-24)中可以明显看出质量矩阵与模态坐标显性相关,而且由于引入转换矩阵 和,质量矩阵也与系统的方向坐标显性相关。质量矩阵中的9个惯性时不变矩阵 可通过计算有限元模型的个节点信息在预处理过程中一次性得到,从而简化运动微分方程的求解。节点信息包括:每个节点的质量 ,未变形时的位置矢量 以及模态矩阵。
(2)势能和刚度矩阵
势能一般分为重力势能和弹性势能两部分,可用下列二次项表示:
.3-25-182)
在弹性势能中, 是对应于模态坐标 的结构部件的广义刚度矩阵,通常为常量。重力势能 表示为:
.3-26-183)
其中 表示重力加速度矢量,重力 可对 求导得;
.3-27-184)
(3)能量损失和阻尼矩阵
阻尼力依赖于广义模态速度并可以从下列二次项中推导得出:
.3-28-185)
上式称为Rayleigh能量损耗函数。矩阵 包含阻尼系数 ,它是常值对称阵。当引入正交模态振型时,阻尼矩阵可用对角线为模态阻尼率的对角阵来表示。对于每一个正交模态的阻尼率都可以取不同值,而且还能以该模态的临界阻尼 的比值形式给出。
3.2.3.2.3
多柔体动力学方程
柔性体的运动方程从下列拉格朗日方程导出:
.3-29-186)
其中, 为约束方程;
为如式(2.3-19)定义的广义坐标;
为投影到 上的广义力;
为拉格朗日项,定义为 , 和 分别表示动能和势能;
表示能量损耗函数;。
将求得的 代入式(2.3-29),得到最终的运动微分方程为:
.3-30-187)
其中, 为柔性体的广义坐标及其时间导数;
2.4 多体系统动力学方程的求解
如2.2.2节所述,在建模和求解过程中,涉及到几种类型的运算和求解:初始条件计算、数学模型自动组装、运动学分析、动力学分析、逆动力学分析和静平衡分析。初始条件计算是非线性位置方程的求解;数学建模是系数矩阵操作;运动学分析是非线性的位置方程和线性的速度、加速度方程的求解;动力学分析是二阶微分方程或二阶微分方程和代数方程混合问题的求解;逆向动力学分析是线性代数方程组求解;静平衡分析从理论上是线性方程组的求解。总的来说,计算多体系统动力学涉及的基本运算包括线性方程组求解、非线性方程组求解、常微分方程组(ODE)求解和微分代数方程组(DAE)求解。
线性方程组和常微分方程组的求解是数值分析中的基本内容。线性方程组根据问题规模可以采用消元法或迭代法,在计算多体系统动力学中,一般采用全主元高斯消元法;常微分方程组的求解可以采用线性多步法或单步的龙格-库塔法。
在这里先简要地给出求解非线性方程组的牛顿-拉夫森(Newton-Raphson)方法,再详细地介绍计算多体系统动力学中至关重要的微分代数方程组求解技术。
2.4.1 非线性代数方程组求解
这里介绍应用广泛的求解非线性代数方程组的牛顿-拉夫森方法。考虑含有n个变量的n个非线性方程组,:
.4-1-188)
其中, ,; ,;0为n维列向量。
令 为式(2.3-1)的解,记 为 的近似,将 在 进行一阶泰勒展开为:
.4-2-189)
其中 为雅可比矩阵,见2.2.3之第2小节部分说明。
令 为改进的近似,如果 足够小,则有:
.4-3-190)
令 ,式(2.4-3)化为:
.4-4-191)
如果雅可比阵 非奇异,则式(2.4-4-191)有唯一解,且得到:
.4-5-192)
于是可根据式(2.4-4)与(2.4-5),由某一初始值开始,逐步迭代,直到收敛,得到有效结果,或者不能收敛,初始值不合适或牛顿-拉夫森方法失效。
2.4.2 微分代数方程组求解
在2.1.3节第2小节中,综述了目前微分代数方程求解的常用技术,根据对广义坐标处理方式的不同,可以分为增广法和缩并法。缩并法一般是在处理过程中将微分代数方程转化为二阶常微分方程再进行求解,增广法更多的是将微分代数方程直接转化为代数方程组的求解序列。
在80年代后期至90年代,微分代数方程求解中缩并法占了主要地位,但90年代增广法展现出新的面貌,体现出其独特的优势。这里介绍F. A. Potra提出的基于超定微分代数方程组(ODAEs - Overdetermined Differential Algebraic Equations)方法的解耦ODAE方法,虽然其数学处理比较复杂,但理解起来是比较简单的,它主要采用微分流形的技术,将微分代数方程组转化为代数方程组的求解序列。该方法充分体现了求解DAE问题中所采用预估-校正策略和显式-隐式方法的区别。
1.2.4.2.1 约束多体系统动力学方程的完整形式
式(2.2-80-87)~(2.2-83-90)给出了约束多体系统动力学方程的完整形式,在这里重写为:
.4-6-193)
.4-7-194)
,
.4-8-195)
,
.4-9-196)
其中, 、 、 分别为系统位置、速度、加速度向量, 是拉格朗日乘子向量, 是时间, 为机械系统质量矩阵, 为约束雅可比矩阵,为外力向量, 为位置约束方程, 和 分别称为速度右项和加速度右项。
系统运动方程(2.4-6-193)~(2.4-9-196)被称为欧拉-拉格朗日(Euler-Langrange)方程,它是微分-代数方程(DAE),DAE不是单纯的常微分方程(ODE)问题,其求解的关键在于避免积分过程中的违约现象。
上述的系统运动方程给出的是最一般的情况。虽然是以二维系统的形式给出,但三维系统如果对各系数矩阵进行适当的组装,可以得到与上述几式一致的运动方程,只是其中各量的含义与二维不同。对于质量矩阵,考虑的是与位置和时间相关的质量矩阵,这是最一般的情况;对于外力向量,考虑的是与位置、速度和时间相关的力,即考虑到了弹簧力、阻尼力等情况;对于约束方程,同时考虑到运动学约束和驱动约束。上述各式给出了明确的函数相关量。
对于上述各式,为了解的存在,一般作如2.2.4节讨论的解的存在性定理中所述的假设:
a. , ;
b. 对任意 且 , 。
2.2.4.2.2
超定微分-代数方程组(ODAEsODAEs)
引入速度矢量 和加速度矢量 ,记:
,
.4-10)-197)
则式(2.4-6-193)~(2.4-9-196)可以另写为:
.4-11-198)
.4-12-199)
.4-13-200)
.4-14-201)
.4-15-202)
.4-16-203)
定义增广广义坐标:
.4-17-204)
在方程组(2.4-11-198)~(2.4-16-203)中,增广广义坐标维数为3n+m,方程组维数为3n+3m,这是一个超定问题,即所谓超定微分-代数方程组(ODAEs),但在前述假设条件下它是相容有解的。
由于系统内在的相容一致性,上述方程组中虽然存在冗余,但可以采取某种方式消除这种冗余性。根据解的存在性定理中的假设,可以认为冗余不会存在于式(2.4-13)~(2.4-16)中,只会存在于式(2.4-11-198)和(2.4-12-199)中,对式(2.4-11-198)和(2.4-12-199)采用微分流形“投影”技术以消除超定性。
.4-18-205)
.4-19-206)
式中, 和 。如果取得有效的 和 ,式式(2-198)和(2-199)(2.4-11)和(2.4-12)从n维降到(n-m)维,消除了系统冗余性。可证明,如果取
和满足:
,
.4-20-207)
非奇异,则方程组(2.4-18-205)、(2.4-19-206)、(2.4-13-200)~(2.4-16-203)与方程组(2.4-21-208)~(2.4-26-213)在某点
的某个邻域内同解。的一种较为有效的取法是:
.4-21-208)
其中 满足:
.4-22-209)
作如上处理后,方程组(2.4-11-198)~(2.4-16-203)可写作:
.4-23-210)
此时, , ,消除了系统超定性。
采用线性多步法:
.4-24-211)
令
.4-25-212)
.4-26-213)
在 时,将 、 、 和 分别写作变量 、 、 和 ,则 可化为:
.4-27-214)
3.2.4.2.3
基于ODAEs的显式方法
考虑显式方法。采用显式线性多步法时,式(2-.4-24211)中 。
在 时,取 值为:
.4-28-215)
将式(2-.4-27214)中第6式和第4式组成的位置方程组:
.4-29-216)
应用式(2.4-29-216)第一式和式(2-.4-22209)第一式,引入新的未知因子 ,可有:
.4-30217)
对式(2-.4-30217)应用式(2-.4-28215)和式(2.4-22-209)第二式,可得:
.4-31-218)
式(2.4-31-218)称为式(2-.4-29216)第一式的扩展等价式。将式(2.4-31)替换到式(2.4-29)中,可得:
.4-32-219)
由式(2.4-32-219)计算出 后,再取 值为:
.4-33-220)
将式(2.4-27-214)中第5式和第3式组成的速度方程组:
.4-34-221)
同式(2.4-32-219)的处理,应用式(2.4-22-209)将式(2.4-34-221)改写为扩展的等价式:
.4-35-222)
式中 为新引入的未知因子。
由式(2.4-35-222)计算出 后,根据由式(2.4-27-214)中第1式和第2式组成的方程组:
.4-36-223)
可计算出加速度 和拉格朗日乘子 。
由式(2.4-32-219)求解 时,由于式(2.4-32-219)是非线性方程,运用Newton-Raphson方法。将式(2.4-32-219)应用Newton-Raphson方法:
,
.4-37-224)
其中,:
.4-38-225)
.4-39-226)
迭代初始值
, .4-40-227)
应用式(2.4-37-224)进行迭代,如果收敛则可计算出 。
由式(2.4-35-222)计算
时,同样应用Newton-Raphson方法:
,
.4-41-228)
其中, 表示同式(2.4-38-225);。
.4-42-229)
迭代初始值:
,
.4-43-230)
由式(2.4-36-223)计算
和时,应用Newton-Raphson方法:
,
.4-44-231)
其中, 表示同式(2.4-38)-225);。
.4-45-232)
迭代初始值:
,
.4-46-233)
由上所述,可以写出基于ODAEs的显式算法。
算法2.4-1
基于ODAEs的显式算法
1. 时间 ,给定相容初始位置 和初始速度 ,由式(2.4-36-223)计算初始加速度 和初始拉格朗日乘子 ;
2. 进入时间步 ;
3. 由式(2.4-25-212)和(2.4-26-213)计算 和 ;
4. 由式(2.4-32-219)计算 。利用式(2.4-37-224)~(2.4-40-227)所述的Newton-Raphson迭代方法;
5. 由式(2.4-35-222)计算 。利用式(2.4-41-228)~(2.4-43-230)所述的Newton-Raphson迭代方法;
6. 由式(2.4-36-223)计算 和 。利用式(2.4-44-231)~(2.4-46-233)所述的Newton-Raphson迭代方法;
7. 如果到终止时间步,结束;否则进行下一时间步 ,转3。
4.2.4.2.4
基于ODAEs的预估-校正方法
考虑隐式方法。采用隐式线性多步法时,式(2.4-174-211)中 。
取 和 为:
.4-47-234)
则系统运动方程式(2.4-27-214)可改写为扩展的等价式为:
.4-48-235)
式中 和 为新引入的未知因子。
用式(2.4-22-209)易证明,当且仅当有唯一确定的
使是式(2.4-48)的解时, 为式(2.4-27-214)的解。
式(2.4-48-235)为隐式方程组,现采用预估-校正方法求解。
引入预估式:
.4-49-236)
.4-50-237)
上述两式为预估式,不同于式(2.4-25-212)和(2.4-26-213)所表示的线性多步法的显式部分和。
利用预估式估计 和 ,式(2.4-47-234)则为:
.4-51-238)
如此,形如式(2.4-32-219)、(2.4-35-222)及(2.4-36-223)的系统位置、速度和加速度方程则相应地为:
.4-52-239)
.4-53-240)
.4-54-241)
据此写出基于ODAEs的隐式算法。
算法2.4-2
基于ODAEs的预估-校正算法
1. 时间 ,给定相容初始位置 和初始速度 ,由式(2.4-54-241)计算初始加速度 和初始拉格朗日乘子 ;
2. 进入时间步 ;
3. 由式(2.4-25-212) 、(2.4-26-213) 、(2.4-49-236)和(2.4-50-237)分别计算 、 、 和 ,并设
, , , , , .4-55-242)
4. 对于 ,计算
.4-56-243)
.4-57-244)
.4-58-245)
式中:
.4-59-246)
5. 如果到终止时间步,结束;否则进行下一时间步 ,转3。
2.5 多体系统动力学中的刚性(Stiff)问题
刚性问题存在于多刚体系统动力学某些情形,更普遍地存在于多柔体系统动力学中,是多体系统动力学的一个重要问题。刚性首先是在常微分方程求解理论中提出,并形成了完整的定义和求解理论。常微分方程刚性理论是多体系统动力学中刚性问题的理论基础,这里先介绍常微分方程刚性问题,再讨论多体系统动力学刚性问题。
2.5.1 微分方程刚性(Stiff)问题
1.2.5.1.1
刚性方程与刚性稳定性
为描述刚性方程性质,先考虑线性系统:
,
.5-1-247)
其中 为解向量函数, 为已知向量函数, 为常系数矩阵,设其特征值为 , ,相应的特征向量为 。
定义(2.5-1-247) 如果在线性系统(2.5-1-248)中, 的特征值 满足:
a. , ;
b.
则称(2.5-1-247)为刚性方程,比值 称为刚性比,通常刚性比
达到就认为是刚性的。
对于非线性系统:
.5-2-248)
令 为(2.5-2-248)满足初始条件
的精确解,在解的邻域内对方程(2.5-2)作线性逼近:
.5-3-249)
或
.5-4-250)
其中 是在点 处 的雅可比(Jacobi)矩阵 的值,如果矩阵 的特征值 满足:
a. , ;
b.
则称方程(2.5-22-248)为刚性方程, 称为在 处的局部刚性比。
在刚性方程中,刚性比 ,矩阵 或 是病态的,故刚性方程也称为病态方程或坏条件方程。刚性方程数值积分过程中存在快变分量 和慢变分量的差别,快变分量要求积分步长很小,而慢变分量则使得在该步长条件下计算步数很多,舍入误差较大,这就是刚性方程数值解的困难所在。
考虑到实际问题中可能会出现单个方程情形,或者矩阵的特征值有实部为零或实部为很小正数的情形,可以给出与实际问题更为接近的刚性方程的定义。
定义2.5-2-248 若线性系统(2.5-1-249)满足条件:
1. 矩阵 的所有特征值实部小于不大的正数;
2. 至少有一个特征值的实部是很大的负数;
3. 对应于具有最大负实部的特征值的解分量变化是缓慢的。
则称(2.5-1)是刚性方程。
对于刚性方程数值稳定性的讨论,一般是针对试验方程:
,
.5-5-251)
进行的。
定义2.5-3-249 一个数值方法以定步长h解试验方程(2.5-5-251),得到线性差分方程的解yn,当
时,若,则称该方法对步长h是绝对稳定的。
定义2.5-4-250
一个数值方法称为A稳定的,如果它的绝对稳定域包含整个左半平面 。
对于A稳定,存在如下结论:
1. 任何显式线性多步法(包括显式RK方法)不可能是A稳定的。
2. A稳定的隐式线性多步法的阶不超过2。
3. 具有最小误差常数的2阶A稳定隐式线性多步法是梯形法。
A稳定的条件过于苛刻,满足条件的数值方法太少,为了突破这个限制,可以放宽稳定性条件,给出稳定的定义。再从刚性问题特点出发,给出刚性稳定性定义。
定义2.5-5-251
一个数值方法称为稳定的,如果它的绝对稳定域包含无限的楔形区域 , 。
定义2.5-6-252
一个数值方法称为刚性稳定的,如果它是收敛的,并存在正常数使在区域 是绝对稳定的,而在区域
上具有高精度且是绝对稳定或相对稳定的。
刚性稳定性的定义,是考虑到微分方程(2.5-1)或(2.5-2)的解中具有较大时间常数的项对解有重大影响。
2.2.5.1.2
刚性微分方程的数值方法
对常微分方程组初值问题:
.5-6-252)
其中 为解向量, 为已知向量函数, 为已知初始向量。
对常微分方程组初值问题,常用的方法有三种。
(1).
线性多步法(LMM)
.5-7-253)
当 时为显式公式,当 时为隐式公式。当 时称为单步法,当 时称为多步法。
(2)2. 预估校正法(PECE)
.5-8-254)
.5-9-255)
式(2.5-8-254)为显式预估公式,式(2.5-9-255)为隐式校正公式。
(3)3. 龙格-库塔法(RKM)
.5-1-2560)
.5-11-257)
如果当 时 ,式(2.5-10-256)和(2.5-11-257)是显式RK方法。如果当 时
,而对角元素不全为0,式(2.5-10-256)和(2.5-11-257)称为半隐式RK方法,若此时均相等,则称为对角隐式RK方法,简称为DIRK方法。
在这些求解常微分方程初值问题数值方法的基础上,产生了求解刚性常微分方程的几类方法,分别是以BDF方法为代表的线性多步法和隐式龙格-库塔方法。BDF方法是一类特殊的隐式线性多步法。一般的隐式龙格-库塔方法计算量较大,这里只给出一类特殊的隐式龙格-库塔方法:对角隐式RK方法。
(1).
向后差分公式(BDF)
.5-12-258)
BDF方法是隐式线性多步法,为k步k阶方法。当k=1,2时,BDF方法是A稳定的,当k=3~6时,BDF方法是稳定和刚性稳定的。实用的BDF方法只能取k=1~6。
(2)2. 对角隐式龙格-库塔方法(DIRK)
.5-13-259)
.5-14-260)
.5-15-261)
对角隐式龙格-库塔方法常用的有2级3阶和3级4阶两个公式,都是A稳定的,此外还有A稳定的4级4阶公式,但给不出A稳定的更高阶DIRK公式。DIRK方法解高频振荡的问题比Gear方法(即BDF方法)好,但对高精度问题比不上Gear方法好,且计算量比Gear方法大。
2.5.2 多体系统动力学中Stiff问题
图2.7
图2.8
如图2.7所示k1=10 N/mm,k2=1000000 N/mm,m1=1kg,m2=1kg,计算可知系统的第一阶频率为11.25Hz,第二阶频率为7117.6Hz,为典型的物理刚性系统。如图2.8所示加入阻尼c1=0.1 N/mm/s,c2=100000N/mm/s,系统变为典型的数值刚性系统。物理刚性系统一定存在大刚度弹簧,系统具有高频;数值刚性系统除在刚度方面存在较大差异外,还有一个很重要的特征是对应高频的阻尼较大,使较高频率被基本阻尼掉,而低频则处于未阻尼状态。数值刚性系统用非刚性数值算法会出现数值困难,这时候需采用专用于求解刚性问题的数值方法。
微分-代数方程的求解,无论是缩并法还是增广法,问题还都是归结为常微分方程初值问题的数值求解,只是求解常微分方程的公式或是用于微分-代数方程转化为常微分方程之后,或是用于转化过程之中。为了使求解的数值方法具有普遍性,既可用于求解良性问题,也可用于求解刚性问题,微分-代数方程所用的常微分方程数值方法一般采用的都是上节所述的求解刚性微分方程的方法。除了最常用的BDF方法,隐式RK方法也被考虑用于求解微分-代数方程问题,此外,预估-校正方法也广泛地用于求解微分-代数方程问题之中,它们都可求解存在刚性问题的微分-代数方程。