無網格法在中心剛體一旋轉柔性梁系統動力學分析中的應用

杜超凡 周曉婷 章定國 高祥

摘要：將無網格點插值法、徑向基點插值法、光滑節點插值法用于中心剛體旋轉柔性梁的動力學分析。基于浮動坐標系方法，考慮梁的縱向拉伸變形和橫向彎曲變形，并計入橫向彎曲變形引起的縱向縮短，即非線性耦合項，運用第二類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.