刚柔耦合多体系统动力学建模的子系统法

1、 刚柔耦合系统的物理模型
考虑如图 1 所示的在水平面内作运动的 刚柔耦合系统，其中的 Hub(刚体)作定轴转 动，在其上以悬臂的方式连接着一柔性梁， Hub 上受到力矩 M(t)的作用， 过 Hub 圆心建 立一惯性坐标系 O－XY，在柔性梁上建立 浮动坐标系 o－xy，其中 o 为 Hub 和柔性梁 的连接处，x 轴沿未变形时梁的轴线，其反 向延长线通过转轴 O , y 轴在水平面内，该 坐标系既用来描述梁的大位移运动，也作为 参考系用来描述梁的变形运动。将图１所示
L L 0 i =1 0
Nf
(24,25,26)
1020
Lm = ∑ Vim Ai −
i =1
Nf
f f f 1 f f − T V A Eijm Bi B = E B B ， ∑∑ ∑∑ ∑ ijm i j j m im i 2 i =1 i =1 i =1 j =1 i =1
N
N
N
N
N
当我们取 Ω = 2rad / s 时，即稳态速度较低时，用本文的一次近似模型和零次近似模 型仿真计算得到的梁的末端横向、纵向变形位移和变形速度如图 7~图 10 所示
0.00005 0.00000
L
(19,20,21) (22,23)
′′ ′ ′ X ijm = ∫ EA(ϕ i′′ φ ′j + ϕ i′φ ′j′ ) φ m dx ， Z ijm = ∫ EI φiIV φ ′j + φi′φ j φ m dx
L L 0 0
(
)
Eijm = ∫ ρAK ij φ m dx ， F jm = −∑ Eijm Bi ， Pm = ∫ ρAφ m dx
图1 刚柔耦合系统简图
的刚柔耦合系统分割成两个子系统，子系统 I 是 Hub，子系统 II 是柔性梁。Hub 绕转轴 O 的转动惯量为 J oh ，半径为 a , 均匀柔性梁的长度为 L ，单位长度质量 ρA ，梁的抗弯刚度 为 EI ，拉压刚度为 EA ，梁未变形时对转轴 O 的转动惯量为 J ob 。
刚柔耦合多体系统动力学建模的子系统法*
方建士 章定国 1 （南京理工大学理学院，南京 210094） 摘要： 本文以转动刚体和其上固结的柔性梁组成的刚柔耦合系统为对象， 利用子系统法研究 了该系统的动力学问题并建立了耦合的动力学方程，此动力学方程考虑了“动力刚化”项。 文章所提供的子系统法具有建模效率高，系统方程耦合度低等优点。 关键词：刚柔耦合系统，子系统法，动力学
M z ( 0, t ) = − EI
∂ 2u y ( 0, y ) ∂x 2
= − EI ∑ φi′′( 0 ) Bi ( t )
i =1
Nf
(12)
Nf ∂M z ( 0, t ) = − EI ∑ φi′′′( 0 ) Bi ( t ) Vy ( 0, t ) = ∂x i =1
(13)
其中
φi′′(0) 和 φi′′′(0) 分别为 φi 在 x = 0 处的二阶和三阶偏导数。式(2)又可表示为
(7,8)
K ij = ∫ φi′(ξ ) ⋅ φ ′ j (ξ )dξ
x 0
i, j = 1,2, " , N f
(9)
轴力 N ( x ) 取值如下：
f ∂u N ( x ) = EA x1 = EA∑ ϕ i′ ( x ) ⋅ Ai (t ) ∂x i =1
N
(10)
或当忽略 u x1 时 N ( x ) 用下式近似计算：
Tip-ux位移(m)
0.0002 0.0000 -0.0002 -0.0004 -0.0006
本文 模 型 传 统模 型
本 文模型 传统模型
0
5
10
15
20
-0.0008
0
5
10
15
20
t(s)
t(s)
图 3 柔性梁末端 x 方向的变形位移
图 4 柔性梁末端 x 方向的变形速度
1021
0.00 -0.02
3 ∂2ux ∂ ∂ uy ∂uy ∂ 2 + ρ Aθ ⋅ u + ρ Aθ 2 ( x + u ) + 2ρ Aθ ⋅u y (3) ⋅ = N ( x ) + ρ Aaθ ρ A 2 − EI ⋅ y x 3 ∂t ∂x ∂x ∂x ∂x
ρA
∂ 2u y ∂t
1019
N ( x) =
1 ρ Aθ 2 ( L2 − x 2 ) − ρ Aaox ( L − x ) 2
(11)
截断数， φi′(ξ ) 是 φi (ξ ) 对 ξ 的一阶偏导， aox 是 o 点加速度在 x 轴上的投影。 根据柔性梁变形与内力的关系得
ϕ i (x ), φi ( x ) 是柔性梁第 i 个纵向和横向变形模态， Ai (t ), Bi (t ) 是模态坐标， N f 是模态
2
+ EI
∂ 4u y ∂x
4
=
∂u y ∂ 2 x N ( x) ⋅ − ρ Aaθ − ρ Aθ ( x + ux ) + ρ Aθ u y − 2ρ Aθ ⋅ u ∂x ∂x
(4)
式中 N ( x ) 为梁 x 处的轴力， u x 和 u y 分别为柔性梁的变形位移在 x 轴和 y 轴上的投影。 2.3 整个系统的耦合动力学方程
外力矩，它包括系统外力矩 M ( t ) 和柔性梁在 切口处对 Hub 的作用力和作用力矩所产生的 效应，即
e = M ( t ) − M z ( 0, t ) + Vy ( 0, t ) ⋅ a (2) Mo
图2
子系统间的相互作用
其中 M z ( 0, t ) 和 Vy ( 0, t ) 是柔性梁在切口处受到的作用力矩和作用力，如图 2 所示。 2.2 子系统 II 的动力学方程 子系统 II 是作大位移运动的柔性梁，采用 Euler-Bernoulli 梁模型，考虑梁的纵向和横向 变形， 对梁的微元段作受力分析， 根据牛顿第二运动定律可得微元 x 和 y 方向的动力学方程 分别为[7]：
(27,28,29)
H im = ∫ EIφ iⅣφ m dx ， k m = −∑ Dim Bm ， K m = aPm + Rm + Lm
L 0
N
(30,31,32)
i =1
从上述推导过程看出，两个子系统的动力学方程推导可以相互独立，因此，可以分别采 用适当的动力学原理来建模。在这里子系统 I 是一个定轴转动的刚体，可以方便地推广到多 刚体系统，本文提供的方法适用于末端体带柔性梁的一类刚柔耦合多体系统。
体角速度 ω 的时变规律为
2π t Ωt Ω sin − ω = T 2π T Ω
系统的仿真结果如图 3~图 6 所示
0<t ≤T t >T
0.0000 -0.0005 -0.0010 -0.0015 -0.0020 -0.0025
0.0006 0.0004
Tip-vx速度(m/s)
引言
在航天器、 机器人和高速旋转机构等工程技术领域中， 针对柔性系统动力学的控制问题 人们提出了这样一类由多刚体系统 （或一个中心刚体） 并带有数个柔性附件组成的刚柔耦合 系统。对这类系统的动力学建模面临两个问题：一是建立系统的完整动力学方程，二是动力 学方程中包含柔性梁的“动力刚化”项[1~3]。通常的建模思路是以整个系统为研究对象，求 出相关运动学量，然后再代入合适的动力学方程（诸如 Jourdain 变分原理[4]、Hamilton 变分 原理[5]、拉格朗日方程[6]、凯恩方程[1]等） ，再经过一系列非常繁琐的数学推导后获得系统的 动力学方程。这种用同一种动力学原理进行动力学方程推导的方法，缺乏灵活性，推导过程 冗长， 建模效率不高， 所得动力学方程耦合度高， 且往往是刚性的， 不利于数值积分的求解。 如果说这种建模思路在零次耦合建模时尚可行的话，那么在考虑“动力刚化”项(即一次耦 合建模)时就需要进一步的证明。事实上目前尚少见成熟的多体系统的一次耦合动力学模型 的出现。 本文针对末端体带有柔性附件梁的一类刚柔耦合多体系统动力学建模提出了子系统 另一个是大位移运动的柔性梁， 一个是多刚体系统， 具体思路是将系统分成两个子系统， 法， 先对这两个系统分别建模， 然后再将两个子系统的动力学方程进行合并。 用这种方法推导系 统动力学方程，难度下降，效率高，所得动力学方程耦合度低，有利于数值积分求解。
* m i =1
L
Nf Nf
2
i =1 j =1
∑D
i =1
Nf
im i
T ， R = L ρAxφ dx B − 2θ m m m ∫
0
(17,1Hale Waihona Puke Baidu)
Dim = ∫ ρAφ iφ m dx ， Vim = ∫ ρAϕ iφ m dx ， Yim = ∫ EAϕ i′′ φ m dx
0 0 0
L
= M e = M ( t ) + EI φ ′′( 0 ) B ( t ) − a ⋅ EI φ ′′′( 0 ) B ( t ) J ohθ ∑i i ∑i o i
i =1 i =1 Nf Nf
(14)
把式 (5) ～ (11) 代入式 (3) ～ (4) 再在左右两边同乘以 φ m ( x )dx 并对梁在全长 L 上进行积 分，其中 m = 1,2, " , N f ，再和式(14)耦合，经过运算最后得常微分形式的动力学方程为
0 θ FN f 1 A 1 # # FN f N f A Nf DN f 1 B1 # # DN f N f BN f
e Mo Q1 # = QN f * Q1 # * QN f
J oh k 1 # kN f K1 # KN f
其中
Nf Nf
0 V11 # V1N f 0 # 0
" " % " % "
0 VN f 1 # 0 # 0
Nf Nf
0 F11 # F1N f D11 # D1N f
" " % " " % "
Nf
" VN f N f
(15)
2 D B B Qm = ∑∑EijmB θ ∑ im i i j + ∑∑ZijmBi B j + ∑Yim Ai +θ (aP m + Rm + Lm ) + 2
i =1 j =1 Nf i =1 j =1 i =1 i =1
Nf
(16)
Q = −∑Him Bi + ∑∑ X ijm Ai B j +θ
2、 系统动力学方程
* 1
基金资助：教育部留学回国人员科研启动基金和南京理工大学科研发展基金 章定国：通信作者
1018
2.1 子系统 I 的动力学方程 对于子系统 I 这个简单的系统， 其动力学方程 为：
= M e J ohθ o
e
(1)
为转动角加速度， M o 为 Hub 所受到的 其中 θ
3、 数值仿真实例
对图 1 的系统取如下的结构参数：梁长 L=5m ，截面积 A = 7.3 × 10 m ，惯性矩
2 −5
I = 8.218 × 10−9 m 4 ， 体 积 密 度 ρ = 2.76667 ×103 kg / m3 ， 弹 性 模 量 E = 6.895 ×1010 N / m 2 ，中心刚体半径 a = 0 ，加速时间 T = 15s ， Ω = 8rad / s 。中心刚
为了求解上述偏微分形式的动力学方程，设
u x ( x, t ) = u x1 + u x 2 , u y ( x, t ) = ∑ φi ( x )Bi (t )
i =1
Nf
(5,6)
其中
u x1 = ∑ ϕ i ( x )Ai (t ) ， u x 2
i =1
Nf
N N 1 f f 1 x ∂u y (ξ , t ) d = ξ − =− ∫ ∑∑ K ij Bi (t )B j (t ) 2 i =1 j =1 2 0 ∂ξ 2
0.02 0.01
Tip-uy位移(m)
-0.04 -0.06 -0.08 -0.10
Tip-vy速度(m/s)
0.00 -0.01 -0.02 -0.03
本 文模型 传统模型
本 文模型 传统模型
-0.12
0
5
10
15
20
0
5
10
15
20
t(s)
t(s)
图 5 柔性梁末端 y 方向的变形位移
图 6 柔性梁末端 y 方向的变形速度