无网格法在中心刚体一旋转柔性梁系统动力学分析中的应用

杜超凡 周晓婷 章定国 高祥

摘要：将无网格点插值法、径向基点插值法、光滑节点插值法用于中心刚体旋转柔性梁的动力学分析。基于浮动坐标系方法，考虑梁的纵向拉伸变形和横向弯曲变形，并计入横向弯曲变形引起的纵向缩短，即非线性耦合项，运用第二类Lagrange方程推导得到作大范围运动的中心刚体旋转柔性梁系统的动力学方程。将无网格法的仿真结果与有限元法和假设模态法进行比较分析，表明其作为一种柔性体离散方法在中心刚体旋转柔性梁的刚柔耦合多体系统动力学的研究中具有可推广性。

关键词：多体系统;动力学;柔性梁;刚柔耦合;无网格法

中图分类号：0313.7;0322

文献标志码：A

文章编号：10044523（2022）01-0178-10

DOI： 10.1638 5/j .cnki.issn.10044523.2022.01.019

引 言

工程中的复杂系统通常由多个独立物体组成，是通过一系列的几何约束连接起来的，并能完成预期动作的一个整体。其中很多都属于柔性附件搭载于刚性主体的刚一柔耦合结构，如空间机械臂、人造卫星天线、直升机旋翼等。上述实际工程问题都可以简化为中心刚体一旋转柔性梁这类典型的刚一柔耦合系统进行动力学分析，而能否正确合理地描述柔性体的变形场将对仿真结果的精度及计算效率产生重要的影响。因此，对变形场离散方法的研究一直是刚一柔耦合系统动力学领域的热点和难点，也是工程和科学领域的迫切需求。目前在该领域被广泛应用的主要是假设模态法（Assumed Mode Method，AMM）和有限元法（Finite Element Method，FEM）[1-2]。假设模态法通常截取较少的模态即可获得较好的数值结果，因此系统白由度少，建立的动力学方程规模也较小，计算效率很高，且便于计算机编程，在数值仿真与实时控制方面具有一定的优势。但当柔性体形状不规则或系统结构复杂时，选取模态函数将变得非常困难[3]。有限元法是将无限白由度的连续体离散为有限白由度的单元集合，以单元节点的弹性位移及一阶导数为广义坐标，在每个单元内建立关于节点坐标的形函数，通过对单元矩阵的组装获得系统的动力学方程[4-5]。发展至今，已有多款商业化的软件可对实际工程中涉及固体及结构的问题进行求解。在柔性多体系统领域，有限元法同样是最常用的离散化手段，与假设模态法相比，其自由度数目往往比较庞大，导致计算效率较低。此外，受单元限制，不易构造高阶连续的形函数，且应力在单元之间并不连续。国内外已有学者开始寻找新的变形场离散方法并应用于旋转柔性梁的动力学研究中，如B样条插值[6-7]及Bezier插值方法[8]。

近年来，无网格法（ Meshless Method）作为一种较新的离散方法得到了迅速发展，成为了研究的热点[9]。无网格法在建立整个问题域的系统代数方程时，只需节点信息，无需划分网格，克服了有限元法前处理复杂的缺点。在构造形函数的过程中，采用更多的节点插值，通常具有高阶连续性，从而提高了计算精度。现有的无网格法有多种，如无网格点插值法（ Point Interpolation Method，PIM）[10-11]、无网格径向基点插值法（Radial Point Interpolation Method，RPIM）”、无单元迦辽金法（Element- free （JalerkinMethod，EFG）”、无网格局部Petrov - Galerkin法（ Meshless Local Petrov-Galerkin Method .MLPG）[14]、再生核粒子法（Reproducing Kernel Particle Meth-od.RKPM）[15]和光滑节点插值法（Node -basedSmoothed Point Interpolation Method， NS_PIM）[16-17]等。目前，将无网格法应用于中心刚体一柔性梁系统的研究鲜有报道[18-21]。

本文采用无网格点插值法（ PIM）、径向基点插值法（RPIM）和光滑节点插值法（NS-PIM）描述柔性梁变形，并在此基础上对旋转柔性梁进行动力学分析。考虑柔性梁的纵向拉伸变形和横向弯曲变形，并计人由横向弯曲变形引起的纵向缩短，即非线性耦合项。采用浮动坐標系描述系统运动，运用第二类Lagrange方程建立系统的动力学方程，编制了相应的仿真软件，通过大范围运动已知和大范围运动未知时的算例，将仿真结果与假设模态法、有限元法等传统离散法所得结果进行对比，表明无网格法应用于该领域的正确性。

1 旋转柔性梁动力学模型

1.1 系统的动能与势能

图1为水平面内运动的中心刚体一旋转柔性梁系统，中心刚体在平面内绕固定转轴旋转，其上以悬臂方式固结柔性梁。以刚体转动中心O为原点建立惯性坐标系Oij，在柔性梁上建立浮动坐标系 。中心刚体的转动惯量和半径分别为Joh和α;柔性梁的弹性模量、长度、密度、横截面积、截面惯性矩分别为E，L，p，S，I。r为刚体上合外力关于转动中心O的主矩。

2 中心刚体一旋转柔性梁系统横向弯曲

固有频率分析

当中心刚体一旋转柔性梁系统作大范围运动时，其纵向振动相比于横向振动，往往忽略不计。采用一次近似模型，忽略梁的纵向振动效应，则其横向振动方程可由式（15）得到：

通过求解式（30）即可求得旋转柔性梁横向弯曲振动的各阶无量纲固有频率。

表1～3为中心刚体无量纲半径比率δ=0，大范围转动速度恒定的情况下，5种不同方法下的旋转柔性梁横向弯曲前三阶无量纲固有频率，假设模态法的模态截断数取3;有限元法取10个单元;无网格点插值法和径向基点插值法将柔性梁离散为11个节点，其中径向基点插值法中形状参数取q=1.03，ac= 4.0;光滑节点插值法将梁离散为81个节点。从表中可知，5种离散方法的旋转柔性梁横向弯曲振动前三阶无量纲固有频率基本一致。随着旋转角速度的增加，前三阶固有频率随之增大，且转速越快，增幅越大。从理论上分析，是由于矩阵θ2（aC+D -M33）+K2产生动力刚化效应，随着旋转角速度的增加，系统的总体刚度矩阵变大，因而其固有频率也随之增大。光滑节点插值法结果总是比其他方法的小，说明其能提供固有频率下界值的特性。随着旋转角速度的增大，有限元法、无网格点插值法、径向基点插值法和光滑节点插值法的结果基本一致，而假设模态法的结果与其他4种方法的误差越来越大。究其原因是因为转速的增大导致柔性梁变形变大，而假设模态法采用的是基于小变形假设下悬臂梁的模态函数，因而在变形较大的情况下会导致误差变大，同时也说明假设模态法不适用于高转速的情况。gzslib202204041803

3 大范围运动已知的动力学仿真

当大范围运动规律已知，即柔性梁的转角规律为已知，式（15）可转化为如下形式：

假设柔性梁由静止开始作大范围旋转运动，展开角速度规律为：式中 T=15 s，15 s后转速到达Ω0，然后以该速度匀速旋转。Ω0分别取为Ω0=4 rad/s和Ω0=20 rad/s。梁的参数取值与文献[23]中相同，分别为：长度L=8m，横截面积S=7.2968×10 -5 m2，截面惯性矩I =8.2189×10 -9 m4，体积密度P=2766.7 kg/m3，弹性模量E= 68.952 GPa。

图2（a）表示Ω0=4 rad/s时，5种离散方法计算的柔性梁末端的轴向拉伸量ω1;图2（b）表示柔性梁末端由横向弯曲变形引起的纵向缩短量即二次耦合变形量ωc;图2（c）为柔性梁末端的纵向变形u1。从图中可知，5种离散方法所得结果相当吻合。梁末端的纵向变形u1为负值，说明梁的轴向拉伸量比横向弯曲变形引起的纵向缩短量小，比较图2（a）和（b）的数值可知，轴向拉伸量为小量，可忽略不计，纵向变形主要由二次耦合变形量ωc导致。

图3（a）表示Ω0=4 rad/s时柔性梁末端横向弯曲变形。图3（b）为15～20 s匀速转动时的局部变形放大图。如图所示，5种离散方法的振幅及频率基本一致，且与文献[8]的仿真结果吻合。

图4（a）表示Ω0=4 rad/s时柔性梁末端横向变形速度，图4（b）为15--20 s匀速转动时的局部变形放大图。如图所示，5种离散方法的仿真结果基本一致。

图5表示Ω0=10 rad/s时柔性梁末端纵向变形，从图中可知，假设模态法的仿真结果与其余4种方法的结果稍有不同，说明随着转速的增大，变形增大，基于小变形假设的假设模态法的误差越来越大。图6表示Ω0=10 rad/s时柔性梁末端横向弯曲变形，图7表示Ω0=10 rad/s時柔性梁末端横向弯曲变形速度，如图所示，假设模态法的仿真结果同样与其余4种方法的结果稍有不同。其余4种离散方法的仿真结果同样基本一致，但相对而言，有限元法、无网格点插值法、径向基点插值法的仿真结果更为接近，从15～20 s的局部放大图中可以看出，光滑节点插值法的振动频率比其他方法略小，这是由于其能提供频率下界值的特性造成的。

图8表示Ω0= 20 rad/s时柔性梁末端纵向变形，图9表示Ω0=20 rad/s时柔性梁末端横向弯曲变形，图10表示Ω0=20 rad/s时柔性梁末端横向弯曲变形速度。从图中可更明显地看出，有限元法、无网格点插值法、径向基点插值法和光滑节点插值法的仿真结果基本一致，但假设模态法的仿真结果已经与其他4种方法出现较大的差别，说明随着转速的增大，假设模态法的精度逐渐降低，进一步说明基于小变形假设的假设模态法不适用于高转速的情况。从15--20 s的局部放大图中可以看出，有限元法、无网格点插值法和径向基点插值法的仿真结果更为接近，光滑节点插值法由于其偏柔性的特性，相对于其他方法，其振幅偏大，频率偏小，利用该特性，结合有限元等方法可从上下界最大程度地获得最优解。

由以上分析可知，随着旋转角速度的增大，假设模态法与其余4种方法的误差越来越大。由式（15）可知，转速的增大导致广义力变大，引起梁变形的增大。而随着梁变形的增大，基于小变形假设的假设模态法误差也越来越大，因此推测假设模态法并不适用于梁变形较大的情形。为了验证这一结论，将梁的弹性模量减小10倍，观察此时5种离散方法梁末端的横向弯曲变形仿真结果。图11表示Ω0=4 rad/s，E=6.8952 GPa时柔性梁末端横向弯曲变形。如图所示，梁的最大变形超过了5m，属于大变形，有限元法、无网格点插值法、径向基点插值法和光滑节点插值法的仿真结果仍然基本一致，而假设模态法很快就发散，说明源于结构力学中固有振型的假设模态法，适用范围仅局限于小变形情况，不能处理大变形问题，而其余4种离散方法均适用于大变形问题。

表4和5分别表示Ω0=10 rad/s和Ω0=20 rad/s时5种离散方法的计算相对时间、大范围旋转角速度恒定时的响应振幅及相对误差。各方法取相同的时间步计算，计算相对时间以3阶模态截断数的假设模态法为标准，相对误差以有限元法为标准。其中，假设模态法纵、横模态截断数分别各取3～7阶。从表中可以看出，模态截断数为3阶的假设模态法计算效率最高，但随着转动速度的增加，大范围旋转角速度恒定时的振幅与有限元法对比误差越来越大，当Ω0=20 rad/s时，误差达到13.86%;同一转动角速度下，假设模态法的模态截断数越多，计算效率越低，但精度并没有明显提高，因此通常认为的通过增加模态截断数来增加假设模态法的精度并不可行;无网格点插值法、径向基点插值法和光滑节点插值法在同一精度误差下计算效率均略高于有限元法，其中光滑节点插值法由于其不包含转角白由度，因而效率最高;精度方面，大范围恒定时的振幅与有限元法基本一致，且光滑节点插值法由于其偏柔性的特点，振幅值较其他方法略大，因而可从上下界获得振幅的最优解。

4 大范围运动未知的动力学仿真

当大范围运动规律未知，即有外驱动力作用于中心刚体时，系统的转动角位移和角速度变为未知，求解式（15）可得柔性梁的动力学响应。梁的参数与大范围运动已知时相同，作用于中心刚体的外驱动力规律为：式中 T= 10 s.10 s后撤掉外驱动力矩。τ0为外驱动力矩的最大值。

图12和13分别表示τ0=1 N.m和τ0=10 N.m时，大范围运动的角位移和角速度以及柔性梁末端的纵向变形u1和横向变形uy。从图中可看出，5种离散方法的仿真结果基本一致。但在去掉外力矩作用时，各方法的响应频率出现差别，这是由于去掉外力矩作用的时刻，各方法仿真的变形状态不同造成的。由图12（a）与13（a）大范围运动角位移可知，当去掉外力矩作用时，柔性梁做小幅摆动，因此梁末端纵向变形和横向变形图与大范围运动已知时的相比很不光滑，出现小锯齿的情形，表明变形的小幅振荡，说明当有外力矩作用时，柔性梁的响应更加复杂。其中光滑节点插值法的梁末端变形稍大于其他方法，体现其偏柔性的特点。比较图12和13中的梁末端纵向变形和横向变形的数值可知，两者不是一个数量级的，因此可忽略其纵向变形。但随着外力矩的增大，梁末端纵向变形与横向变形的数量级差在减小，说明在外力矩增大到一定程度时，则必须考虑梁纵向变形的影响。图14为τ0=50 N.m时，4种离散方法在同样的计算条件下得到的柔性梁末端的横向弯曲变形。由图中可看出FEM、PIM和RPIM的结果很快发散，而NS-PIM仍然收敛，说明了NS-PIM计算上的优势。同时说明PIM和RPIM在中心刚体一柔性梁系统的动力学响应计算上并没有显示出无网格法相对于有限元法的优势，因为对于一维梁单元的离散形式，两者无区别。gzslib202204041803

5 结 论

（1）在转速较低时，5种离散方法所得到的柔性梁横向弯曲前三阶无量纲固有频率基本一致。随着转速的增加，假设模态法的误差越来越大，而其余4种离散方法结果基本吻合，说明基于小变形条件的假设模态法并不适用于高转速的情形。其中光滑节点插值法的结果总是比其他方法的小，能提供固有频率的下界值。

（2）对于假设模态法而言，增加其模态截断数对精度的影响并不大，但计算时间却成倍增加。无网格点插值法、径向基点插值法和光滑节点插值法与有限元法结果基本一致。计算效率方面，三种新方法均比有限元法略高。

（3）假设模态法适用范围仅局限于小变形情况，不能处理大变形问题。外力矩较小时，梁的变形较小，因而可忽略梁纵向变形的影响;随着外力矩的增大，梁的变形增大，且梁的纵向与横向变形数量级在减小，此时应考虑梁纵向变形的影响。仿真结果也说明NS-PIM的适用范围最广。

参考文献：

[1]洪嘉振，尤超蓝.刚柔耦合系统动力学研究进展[Jl.动力学与控制学报，2004，2（2）：1-6.

Hong Jiazhen， You Chaolan. Advances in dynamics ofrigid-flexible coupling system[J]. Journal of Dynamicsand Control， 2004，2（2）：1-6.

[2]Dwivedy S K， Eberhard P.Dynamic analysis of flexiblemanipulators，a literature review[Jl. Mechanism andMachine Theory， 2006， 41（7）：749-777.

[3]吳胜宝，章定国，康新.刚体微梁系统的动力学特性[J].机械工程学报，2010，46（3）：7682.

Wu Shengbao， Zhang Dingguo， Kang Xin. Dynamicproperties of hub-microbeam system[J].Journal of Me-chanical Engineering， 2010， 46（3）：7682.

[4]Chung J，Yoo H H. Dynamic analysis of a rotating can-tilever beam by using the finite element method[Jl.Journal of Sound and Vibration， 2002， 249 （1）：147-164.

[5]Du H， Lim M K， Liew K M.A nonlinear finite elementmodel for dynamics of flexible manipulators[J].Mecha-nism and Machine Theory， 1996， 31（8）：1109-1119.

[6]Lan P， Shabana A A. Integration of B-spline geometryand ANCF finite element analysis[Jl. Nonlinear Dy-namics， 2010， 61： 193206.

[7]Liu Y N， Sun L， Liu Y H， et al. Multiscale Bsplinemethod for 2D elastic problems[J].Applied Mathemat-ical Modelling， 2011， 35： 36853697

[8]范纪华，章定国.旋转柔性悬臂梁动力学的Bezier插值离散方法研究[J].物理学报，2014. 63（15）：154501.

Fan Jihua， Zhang Dingguo. Bezier interpolation methodfor the dynamics of rotating flexible cantilever beam[J].Acta Physica Sinica， 2014， 63（15）： 154501.

[9]Belytschko T， Krongauz Y， Organ D， et al. Meshlessmethods： an overview and recent developments[Jl.Computer Methods in Applied Mechanics and Engineer-ing， 1996， 139：347.

[10]Liu G R， Gu Y T.A point interpolation method fortwodimensional solids[J].International Joumal for Nu-merical Methods in Engineering， 2001. 50 （4） ：937-951.

[11]Liu G R， Dai K Y， Lim K M， et al. A point interpolation mesh free method for static and frequency analysisof twodimensional piezoelectric structures [ J] . Compu-tational Mechanics ， 2002 . 29（ 6） ： 510519.

[12]Wang J G， Liu G R. A point interpolation meshlessmethod based on radial basis functions [J] . InternationalJournal for Numerical Methods in Engineering， 2002，54（ 11） ： 1623-1648.gzslib202204041804

[13]Belytschko T， Lu Y Y， Gu L. Element free Galerkinmethods [Jl. Intemational Journal for Numerical Meth-ods in Engineering ， 1994 ， 37（ 2） ： 229256.

[14]Atluri S N. Zhu T. A new meshless local Petrov-Galerkin （ MLPG） approach in computational mechanics [Jl.Computational Mechanics ， 1998 ， 22（ 2） ： 117-127.

[15]Liu W K， Jun S， Zhang Y F. Reproducing kernel parti-cle methods [Jl. International Journal for NumericalMethods in Fluids ， 1995 ， 20 ： 1081-1106.

[16]Liu G R， Zhang G Y， Dai K Y. A linearly conformingpoint interpolation method （ LC-PIM） for 2D solid mechanics problems [ J] . International Journal of Computational Methods ， 2005 ， 2（ 4） ： 645665.

[17]Liu G R， Zhang G Y. Upper bound solution to elasticityproblems ： A unique property of the linearly conformingpoint interpolation method （ LC-PIM） [J] . InternationalJournal for Numerical Methods in Engineering， 2008，74 ： 1128-1161.

[18]Chen Yuanzhao， Zhang dingguo， Li liang. Dynamicanalysis of rotating curved beams by using absolute nod-al coordinate formulation based on radial point interpolation method [ J] . Journal of Sound and Vibration ， 2019 ，441： 6368.

[19]Xie Dan. Jian Kailin. Wen Weibin. An elementfreeGalerkin approach for rigid-flexible coupling dynamicsin 2D state J].Applied Mathematics and Computation，2017.310：149-168.

[20]謝丹，蹇开林.改进EFG法用于旋转梁的刚柔耦合动力学研究[J].振动工程学报，2017， 30（4）：527-534.

Xie Dan. Jian Kailin. An improved EFG approach forrigid-flexible coupling dynamics for the rotating hubbeam system[Jl. Journal of Vibration Engineering，2017， 30（4）：527-534.

[21]陈渊钊，章定国，黎亮.平面细长梁基于无网格径向基点插值的绝对节点坐标法[J].振动T程学报，2018， 31（2）：245254.

Chen Yuanzhao， Zhang Dingguo， Li Liang. An abso-lute nodal coordinate formulation based on radial pointinterpolation method for planar slender beams[J]. Jour-nal of Vibration Engineering， 2018， 31（2）：245254.

[22]杜超凡.基于无网格法的刚柔耦合系统的动力学建模与仿真[D].南京：南京理工大学，2017.

Du Chaofan.A study on the dynamic modeling and simulation for the rigid-flexible coupled system based onmeshless methods[D].Nanjing： Nanjing University ofScience and Technology， 2017.

[23]和兴锁，李雪华，邓峰岩.平面柔性梁的刚柔耦合动力学特性分析与仿真[J].物理学报，2011， 60 （2）：377-382.