多体系统动力学约束Hamilton方程多步块方法

® t2d+n,)<ln]+ [h {B ® I2d+m )qn+l +
('®D
]}，
■^23 = <^,'
1 2d * n M n +\ + ( ^ ®
)^,] >
2d+m ) ^ +1 +
(卵 K
]，
*^22=^d >
^33= 〇
(6)
4 ) 将上述非线性方程组记为G „+1，求下列方程，即
Z ：：；= Z ：+I- G ：+1/ J
1 引言
约 束 H a m i l t o n 系统[1]作为重要的动力系统之
一 ，在天体力学、多体系统动力学、机械制造以及 航空航天[2]等领域广泛应用。目 前 Hamilton系统正 则微分-代数方程主要解法是几何数值积分方法，包 括辛算法、能量方法、变 分 数值积分方法和L i e 群 方法等。1986年 ，冯康W 首次提出辛算法的概念。 文献[4]提出了辛精细积分算法与最优控制理论。文 献 [5]利 用 矩 阵 形 式 的 分 块 算 法 ，提出了树形多体 H a m i l t o n 系 统 的 隐 式 辛 R u n g e -Kutta算 法 。文献[6] 导 出 了 功 能 梯 度 圆 柱 壳 的 H a m i l t o n 正 则 方 程 ，并基 于辛方法分析了其自由振动特性。文献[7]将精细算 法 与 有 限 元 算 法 结 合 起 来 ，研 宄 了 精 细 辛 有 限 元 方 法 及 其 相 位 误 差 。文 献 [8]提 出 了 约 束 多 体 系 统 任 意 精 度 阶 L i e 群积分方法，并提出用约束最小二乘极
通 信 作 者 ： 丁洁玉，女 ， 1978年 生 ，博 士 ，青岛大学，教 授 ：研究方向— 多体系统动力学与控制。 E-mail: djy@qdu.edu.cn
第3期
杜文静，等：多体系统动力学约束Hamilton方程多步块方法
1017
多步块方法求解常微分方程的算法有很多，但 多 步 块 方 法 求 解 H a m i l t o n 系统正则微分-代数方程 的研宄却较少。本文基于等距节点和非等距节点建 立多步块格式，求 解 多体系统动力学Hamilton体系 正则微分-代数方程。
以指标- 2 的 Hamilton系统正则微分-代数方程
为例，用 r= 3 的等距节点构造的多步块格式进行求
解 ，选 取 步 长 ；!=0.005s, 仿 真 时 间 为 20s, 求解过
线性代数方程，求 解 过 程 利 用 N e w t o n 迭代。以双连杆机械臂系统为例，使 用 r-级多步块格式进
行数值仿真，通过改变步长、仿真时间以及节点数分别对指标-1、-2、- 3 的 H a milton系统正则微
分 -代 数 方 程 进 行 数 值 验 证 。数 值 结 果 表 明 ，本 文 方 法 穗 定 性 好 、精 度 高 、计 算 效 率 高 ，能很好保
(7)
直到达到最大收敛误差为1 0 ' 求得匕+，的值。将求
方程级多步块方法
3 . 1 多步块方法的构造 设单步区间H ,]，时 间 步 长 /!=
々，在
得 的 1^+,的值带入式(4),求 得 y„+，的 值 ，进而求得
?„+1、 凡+1的值 。上述过程为一个时间步长的求解 结果 ，循环重复上述过程，从 而 得 到 《个时刻的广 义位置向量和广义动量尸。
2 ) 由式(4)可 得 r„+1的值。
3 ) 求非线性方程wenku.baidu.com的Jacobi矩 阵 ，其 Jacobi矩阵
可以用一个块矩阵来表示，即 Ju Jn
J = J 2l J 22 J 23
(5)
、•/3I </32 </33y
其中 J u =I^ J n = - h M -\B ® I2d+m), 7,3=0, J 2r h {B ® I2d+m ){Vqq[h{B ® I2dtm)qn+,+
函数；F 为系统势能；匕为势 能 F 关于广义坐标9 的导数； A = [ 4 ，/^，"•，人J T 为 Lagrange 乘子。
对约束方程<P 关于时间（求一阶和二阶导数，可 得指标-2、-1的 Hamilton系统正则微分-代数方程为
q = M{q,p) 'p
■ 々=
- 汽丨( 9 ,〇又
(2 )
Fig. 1 A planar two-link manipulator
选 取 状 态 变 里 9 = 1^
X2 少2
，
广 义 质 量 矩 阵 M = diag[Wi w 丨 /丨 w : w 2 ,2]，
由 Legendre 变 换 得 系 统 的 广
Mq =
j 丨 % ^2 夕2 92'
m2
m 2 八」
第38卷 第 3 期 2021 年 6 月
应用力
学学报
CHINESE JOURNAL OF APPLIED MECHANICS
文 章 编 号 ： 1000-4939(2021)03-1016-06
Vol.38N o .3 Jun.2021
多体系统动力学约束Ham ilton方程多步块方法
杜文静1 丁洁玉u
( 1 青岛大学数学与统计学院266071青岛： 2 青岛大学计算力学与工程仿真研宄中心266071青岛）
摘 要 ：针对多体系统动力学Hamilton体系正则微分-代数方程，基于等距节点、C h e b y s h e v 节点、
L e g e n d r e 节 点 等 非 等 距 节 点 ，在 时 间 区 间 上 建 立 r-级 多 步 块 求 解 格 式 ，得到单步区间各节点的非
(2r -2)
Cr
(2r -2)! (c,- l)(2，-l>
(2r-l)!
-1广丨) (2r - l)!
Cr ~ l
(cr - )l (2r-l) (2r-l)!
(9) 以 为 例 ，此时等距节点为c,.= //3 , = 1,2,3)。
由式(9)可 求 得 B 、T V 的值。将式(9)代入式(4)可得基 于等距节点的多步块求解公式。
0 t(q,t)M '(q,p)p + C>,{qJ) = 0
q = M(q,p)~'p
p = -Vq(q,t)~<Pj(q,t)A
\
,
(3)
(q,p)p + [0i{q,t)M (q,p)p\-
^ '(?>p)p + 2<P?,{q,t)M '(q,p)p + 0„(q,t) = 〇
本 文 建 立 的 多 步 块 格 式 可 以 分 别 求 解 指 标 -1、 -2、- 3 的 Hamilton系统正则微分-代数方程。
(8)
1!
2r —\
将 式 (8)代 入 式 (4)可得
[B 7V ]=
01
C1(2 r-2 ) (2r —2)!
J 2 r-\) Cl .. (2r - l)!
^(2r_l)
C2
••
C2
(2r - l)!
J 2r-\) Cr •• (2r -\)\
c丨-1 〇2 一1
ci(2r-2) (2r - 2)!
G < rt+1, 构 造 以 c,.=丄,(/=丨，2，...，〇为等距节点或 r
C h e b y s h e v 节点、L e g e n d r e 节点为非等距节点的多
步块求解格式，此 方 法 是 L 稳定的。以指标- 3 的约
束 Hamilton方程(式⑴)为例。
令
y = [q \ p TX ] T e R 2d+m, y ,=
3 . 3 基于非等距节点的多步块方法
3.3.1基 于 Chebyshev节点的多步块方法
当 节点c,.(/= 1，2,…，r)为 C h e b y shev多项式的零
点 时 ，有
+ cos --2-i---\- 71 卜(，.=1,2，" ，厂-1) U (r - l)
( 10)
4 数值算例
图 1 为双连杆机械臂系统，杆 1 质量和杆长为 m ^ l k g ， Z^ l m ，其质心坐标为(i p % ) 。杆 2 的 质 量 和 杆 长 为 w 2 = 2 k g ， L 2 = V ^m ，其质心坐标 为 〇c2,h ：)。转 角 用 七 6 表 示 ，初 始 时 刻 时 3 = n/3, 沒2 = 兀/6 〇
持 Hamilton能量误差以及各级约束误差，特别是对于长时间仿真，Hamilton能量误差不会产生漂
移 ，适合长时间多体系统动力学仿真。
关键词：多体系统动力学；Ha m i l t o n 系统；微分-代数方程；多步块方法
中图分类号：T P 301.6; 0175.1 文献标识码：A
D O I : 10.11776/cjam.38.03.D 148
⑴ 又
式 中 ：9 为 Hamilton系 统 的 广 义 坐 标 ，且
9 = [ 糾 ) ，心 ( 0 ，- - . 而 ( 0 ] 丁 ； 夕 为 系 统 的 动 量 ， 且
T …，凡⑴] ; 々 为 广 义 速 度 ， 且 e 4 /d/ ;M 为对称正定的质量矩阵，且 M /?〜 ；
少 为 m 个完整 约 束 ，且 <Z> = [々,鹆 • ，少„,]T ; 气 为 约 束 少 关于广义坐标7 的 导 数 ；" ( ？，/>)为 Hamilton
3 Hamilton系 统 正 则 微 分 -代 数
将式(4)与 之 +,代入式(1),得到非线性代数方
程 组 ，该 方 程 组 含 有 K 2J + W )个 未 知 量 、r(2d + w )
个方程。
本 文 通 过 N e w t o n 迭代循环求解出1^+,的值。具
体步骤如下。
1 ) 给定也+1、九+1、人+1的初值。
= y^ ) ,
^ - i), y„+l =[yiT ^ J ,---,^T]T,
t+l=[y： ，y2r，- ，y：]\ K =UT-1^ 2T-„-,^ T,]T
构造级多步块求解格式得到
Yn+t = h(B ® I2d+m)Y n+i+ (yv ® I2d+m )Yn (4)
式中 ： J5、j V e /T "; ® 为 Kronecker 积 。
1018
应用力学学报
第 38卷
3 . 2 基于等距节点的多步块方法
对 于 等 距 节 点 。=/",(/ = 1，2，"_，幻 ， 对 兄,(/ = 1,2,…，r) 在 （处进行泰勒展开，得到
y ( J ) « a 0 + ^ ( t - t n ) + --- + - ^ r ± - ( t - t n f rA
小 化 算 法 解 决 约 束 稳 定 问 题 。文 献 [9]将 投影算法与 广 义 法 结 合 起 来 ，将其用于求解多体系统动力学 微分-代数方程，从而解决约束稳定问题。文献[10] 对压电层合梁的动力学分析进行了深入研究。
由 于 R u n g e -Kutta方法的级精度低于方法精度， 文 献 [11]首 次 提 出 了 块 方 法 。文 献 [12]基于多处理机 系统提出了求解刚性常微分方程的并行块隐式 R u n g e -Kutta方 法 ，并且证明该方法具有A 稳定性。 文 献 [13]提 出 了 求 解 常 微 分 初 值 问 题 的 并 行 变 步 长 块格式，并提出了一种四点全角全隐式的多步块方 法 来 求 解 常 微 分 方 程 。文 献 [14]将 多 个 有 限 差 分 公 式 组 合 成 块 向 量 ，提 出 了 一 种 求 解 四 阶 初 值 问 题 的 线性多步混合块法。文献[15]对一般的线性多步法 进行改 进 ，提出了 H amilton系统的对称多步法，并 研 究 了 其 长 时 间 的 性 质 。 文献[16]基 于 一 般 线 性 方 法 和 Rattle方 法 ，提出了 Hamilton系统微分-代数方 程的辛对称多步方法。
基 金 项 目 ：国家自然科学基金(11772166; 11472M3)
收 稿 日 期 ： 2019-11-25
修 回 日 期 ： 2021-04-15
第 一 作 者 简 介 ： 杜文 静 ，女 ， 1995年 生 ，青岛 大 学 ，研究:生 ；研宄方向— 微 分 方 程 数 值 求 解 。 E-mail: 649031683@qq.com
2 多 体 系 统 动 力 学 Hamilton系统
正则微分-代数方程
对于受完整约束的动力学系统，其指标- 3 的
H a m i l t o n 系统正则微分-代数方程可表示为
<•1=
8H
—
(；q--,-p--)=
M(q,p)
p
单 步 区 间 中 选 取 r 个 节 点 ，其 中 G < c,< c2 < 〜 <