中心刚体-柔性梁刚柔耦合动力学模型降阶研究*

李莉 刘铸永 洪嘉振

（上海交通大学工程力学系，上海 200240）

引言

浮动坐标系方法是柔性多体系统动力学最常用的方法，其中柔性体的运动被分成两部分：用浮动坐标系描述的大范围刚体运动和相对于浮动坐标系的小的弹性变形.有限单元法被广泛的采用来描述柔性体的弹性变形，然而有限元节点坐标数目庞大，将会给动力学方程求解带来巨大的计算负担.如何降低柔性体的自由度，是当前柔性多体系统动力学研究的一个重要命题.因此，为了提高柔性多体动力学仿真的计算效率，便于控制设计和实施，就必须要对柔性多体系统动力学模型降阶进行研究.降阶后的模型既要能真实地反映出系统的动力学特性，阶数也要足够低.

模型降阶一般可以从两方面予以考虑［1］：一方面是从建模的角度进行降阶，即根据经典的假设模态法，选择具有较好展开收敛性的模态集，以便用较少的模态对系统进行建模.由于假设模态方法没有考虑作用在物体上的外部载荷的空间特性，其收敛速度很慢.为了提高收敛性，人们提出模态综合法等方法如来改进假设模态方法.但是，如何选取一组恰当的由特征模态和假设模态构成的基矢量集是一件非常困难的事情，因为这取决于工程人员的经验和对此问题的认识程度.模型降阶的另一方面是在模型建立之后，对所建立的模型通过提供合适的降阶方法，进行合理的自由度减缩，从而进一步对模型降阶.从结构动力学或者有限元理论角度出发，目前结构模型的降阶方法主要有：频率截断法、惯性完备性准则和参数匹配法等传统降阶方法.

这些方法仅考虑了模型本身的属性，而忽略了控制系统设计要求以及外界干扰等因素.在现代控制理论中，从系统的可控性、可观性出发，利用一些降阶方法将那些不可观、不可控的模态剔除.Skelto［2］提出的模态价值分析方法和 Moore［3］提出来的平衡准则等，综合考虑了系统本身属性、系统的可观可控性以及外界干扰等因素，是目前比较全面实用的降阶方法.国内学者曲广吉等［1，9］等对柔性航天器动力学模型降阶问题进行了深入的研究，取得了非常良好的结果.章敏、蔡国平［4］研究表明对与频率密集的柔性板，内平衡降阶方法比模态价值分析方法更为有效.在状态空间系统，近年来提出的Krylov子空间模型降阶方法，能够有效的处理大规模动力系统的模型降阶问题［5］.Eberhard［6，7］等对Krylov子空间模型降阶技术进行改进，对柔性体的动力学模型降价进行了很好的研究.总之，当前对柔性结构的模型降阶研究较多，对柔性多体系统动力学的研究相对较少.本文采用Krylov模型降阶方法［8］，对柔性梁刚柔耦合动力学模型降阶问题进行研究.首先采用有限单元法离散弹性变形，建立系统的刚柔耦合动力学方程.然后将柔性梁作为一个输入输出系统，采用Krylov模型降阶方法缩减系统自由度，建立系统降阶的刚柔耦合动力学方程.再次采用传统的模态降阶方法，缩减系统自由度，建立系统降阶的刚柔耦合动力学方程.最后，分别采用有限元全模型、Krylov降阶模型和模态降阶模型进行系统的刚柔耦合动力学仿真.

1 刚柔耦合动力学的模型降阶

做大范围运动的任意一个柔性体的动力学方程可以表示为

其中q表示大范围刚体运动广义坐标，p表示有限元节点坐标，Mr，Me，Mer分别代表刚性质量阵、弹性质量阵、刚柔耦合质量阵，Ke代表弹性刚度阵，Qr，Qe分别为刚体运动和弹性运动的广义力阵.

为了提高柔性多体动力学仿真的计算效率，便于控制设计和实施，就必须要对柔性多体系统动力学模型降阶进行研究.如果有限元节点坐标p（t）∈RN，能够用一个低维的坐标来近似.他们满足如下关系式

其中转换矩阵 V∈RN×s，s≪N.把式（2）带入式（1）降阶的动力学方程可以表示为

降阶的转换矩阵V∈RN×s可以是模态矩阵Φ∈RN×s，也可以通过其它方法如Krylov模型降阶方法得到.

模型降阶过程中暂不考虑大范围刚体运动的影响，仅对柔性体的弹性变形坐标进行降阶.柔性体弹性振动的结构动力学方程为

其中 Be∈RN×m为输入矩阵，Ce∈Rn×N为输出矩阵，u∈Rm×1为载荷向量，y∈Rn×1为观测量，系统的传递函数表示为

不失一般性，下文仅讨论Be＝CTe的工况.引入转换矩阵V∈RN×s，降阶后系统的结构动力学方程为

降阶后系统的传递函数表示为

在相同的输入条件下，降阶前后系统的观测量相同，即降阶前后的传递函数满足以下关系式

假设Ke为非奇异矩阵，传递函数H（s）在σ＝0处泰勒展开为

其中mi∈Rn×m为第i次矩，Krylov模态降阶方法又称为矩匹配法，即降阶后系统传递函数的第 i次矩满足以下关系式

转换矩阵的二阶Krylov形式为［8］

其中数学符号 κr（A1，A2，G）定义为

H（s）在σ≠0处泰勒展开与 H（s+σ）在0处的泰勒展开相等

系统转换矩阵形式为

其中

2 刚柔耦合动力学方程

本节建立了典型的刚柔耦合动力学模型，由一个中心刚体—柔性梁组成.以此为例，说明如何实现柔性多体动力学的模型降阶过程.基于以下基本假设：梁是均匀的，各向同性的，细长的，剪切影响可以忽略；柔性梁的变形为小变形；中心刚体大范围转动的角速度较低，即零次近似耦合模型即可适用.

2.1 运动学描述

如图1所示，中心刚体的半径为r0.柔性梁在固定在刚体的O1点，梁长度为 L.系统在XY平面内转动，重力影响忽略.为了描述系统的运动，建立了两个坐标系：惯性坐标系和浮动坐标系是点O1的位置矢量，是未变形梁上任意一点P0的位置矢量，是变形矢量，梁上任意一点P的位置矢量在坐标系上的位置坐标为

其中 r0＝［r0，0］T，ρ0＝［x，y］T，u＝［u1，u2］T是在坐标系上的位置坐标，A是从到的转换矩阵，给出如下：

其中θ是梁大范围运动的转角.

对（16）式求导可得

其中的变分可以表示为

对（16）式求t的二次导数可得

图1 动力学模型Fig.1 Dynamic Model

2.2 变形描述

采用非笛卡尔变形坐标描述，柔性梁的变形场可用线性梁单元插值函数描述：

此即为单元变形在浮动坐标系上的插值公式，其中pk为单元节点坐标，为平面梁的插值函数.使用布尔定位阵BK，（21）式可以表示成全部节点坐标表示的形式.

其中 p（t）∈RN为梁的整体节点坐标，为整体函数.

2.3 动力学方程

采用Hamilton原理来建立作大范围运动平面柔性梁的动力学方程：

其中，T为系统的动能由中心刚体动能TH和柔性梁动能TB组成，U为系统的势能，W为外力所作的功.

中心刚体和动能有关的变分项可表示为

柔性梁和动能有关的变分项可表示为

柔性梁变形能的变分为

假设作用在柔性梁上的外力为分布参数形式的外力（例如重力），其在连体坐标系下的坐标阵为在惯性坐标系下的坐标阵为f＝Af′.中心刚体受外力矩为Me.那么外力所做的虚功可表示为

将 （24）、（25）、（26）和（27）式代入（23）式，可以得到作大范围运动中心刚体-柔性梁的刚柔耦合动力学方程为

其中 θ，p广义坐标，Mθ，Meθ，Me是分块质量阵，Ke是弹刚度阵，Qθ，Qe是广义外力阵.篇幅关系它们的具体形式不再给出.

降阶后，式（28）可表示为

通过Krylov降阶时W＝V，通过模态降阶时W＝Φ.

3 仿真算例

带末端质量的中心刚体-柔性梁系统的几何和材料参数给出如下：质量密度 ρ＝2.767×103kg／m3，梁长 L＝1.8m，截面惯性矩 I＝1.302×10-10m4，截面积 A＝2.500×10-4m2，弹性模量 E＝6.895×1010N／m2，中心刚体的转动惯量 JH＝0.25kg·m2，中心刚体的半径r0＝0.1m；梁末端有一外力F＝［40.2］T作用；初始时刻梁未变形.

图2 w0＝1rad／s时系统动力学响应Fig.2 Response of the system when w0＝1rad／s

图3 w0＝2rad／s时系统动力学响应Fig.3 Response of the system when w0＝2rad／s

图2，图3给出了中心刚体初始角速度w0＝1rad／s以及w0＝2rad／s时分别采用有限元全模型、Krylov降阶模型和模态降阶模型，进行系统刚柔耦合动力学仿真结果.数值仿真结果表明，与采用模态降阶方法（需要25阶模态）相比，采用Krylov模型降阶方法（需要5阶模态）只需要较低的自由度，就可以得到与采用有限元方法吻合较好的结果.

4 结论

浮动坐标系方法是柔性多体系统动力学最常用的方法，其中柔性体的运动被分成两部分：用浮动坐标系描述的大范围刚体运动和相对于浮动坐标系的小的弹性变形，因此又称之为刚柔耦合动力学.有限单元法被广泛的采用来描述柔性体的弹性变形，然而有限元节点坐标数目庞大，将会给动力学方程求解带来巨大的计算负担.如何降低柔性体的自由度，是当前柔性多体系统动力学研究的一个重要命题.本文以中心刚体-柔性梁系统为例，分别采用Krylov方法和模态方法进行降价.然后分别采用有限元全模型、Krylov降阶模型和模态降阶模型，对中心刚体-柔性梁进行刚-柔耦合动力学仿真.仿真结果表明，与采用模态降阶方法相比，采用Krylov模型降阶方法只需要较低的自由度，就可以得到与采用有限元方法完全一致的结果.说明Krylov模型降阶方法能够有效的用于柔性多体系统的模型降价研究.

1 缪炳祺，曲广吉，夏邃勤等.关于柔性航天器动力学模型降阶问题.中国工程科学技术，2001，3（11）：60～64（Miao BQ，Qu G J，Xia SQ，et al.On the order reduction of dynamicsmodels of flexible spacecrafts.Engineering Science，2001，11（3）：60～64（in Chinese））

2 Skelton R E，Gregory C Z.Measurement feedback and model reduction bymodal costanalysis.In：Proceedings of the Joint Automatic Control Conference.USA：Denver，1979：211～218

3 Moore B C.Principal component analysis in linear system：controllability，observability and model reduction.IEEE Transaction on Automatic Control，1981，26（1）：17～31

4 章敏，蔡国平.密集频率柔性板的内平衡模型降阶及其主动控制研究.工程力学，2009，26（11）：161～167（Zhang M，CaiG P.Balanced reduction and controlof flexible plate with dense frequencies.Engineering Mechanics，2009，26（11）：161～167（in Chinese））

5 Bai Z.Krylov subspace techniques for reduced-order modeling of large-scale dynamical systems.Applied Numerical Mathematics，2002，43（1）：9～44

6 Michael L，Eberhard P.A two-step approach formodel reduction in flexible multibody dynamics.Multibody System Dynamics，2007，17：157～176

7 Fehr J，Eberhard P.Simulation process of flexible multibody systems with non-modal model order reduction techniques.Multibody System Dynamics，2011，25（3）：313～334.

8 Salimbahrami B，Lohmann B.Order reduction of large scale second-order systems using Krylov subspacemethods.Linear Algebra and its Applications，2006，415（2）：385～405

9 次永伟，邱大芦，付乐平等.航天器振动试验控制技术进展.动力学与控制学报，2014，03：193～200（Ci YW，Qiu D L，Fu L P，et al.Progress in spacecraft vibration testing control technology.Journal of Dynamics and Control，2014，03：193～200（in Chinese） ）