共旋框架下二維柔性多體動力學分析策略

中圖分類號：0313.7 文獻標志碼：A DOI：10.16385/j.cnki.issn.1004-4523.202310039

Abstract：This paperproposes anaccrateandeficient solutionstrategyforanalyzing dynamic responsesofflexiblemultibodysys tems.Intheproposedstrategyflexiblestructuresaremodeledinthecorotationalframe，thenthediscretemathematicalmodelis solvedbyanoptimizedcompositemethod.Duetotheintroductionofthecorotationalframe，someadvancedlinearelementsanbe directlyemployed，dramaticallydecreasingcomputationalcosts.Foraccuratelycalculatingdnamicresponses，anoptimizedtee sub-stepcompositemethod isdeveloped whereinalgorithmic parametersareoptimizedforminimizing local truncationerrors.The optimizedcompositemethdaievessecondderaccracyuncondiionalstabilityandcontrollblestabilitySomeclacalfle ibledynamicsystemsaresolved inthispaper，and numericalresultsshowthatcomparedtothecurrntlypopular solutionstrategy based on the absolute nodal coordinate formulation and the Generalized- α method，under the same computational accuracy，our strategy has great superiorities in efficiency.

Keywords：multibodydynamic;flexible multibodysystems;corotational frame;parameters optimization;compositemethod

具有強非線性、大變形、大轉動和大量自由度的柔性多體系統在航空航天、機械制造等領域廣泛存在[12]。由于此類動力學系統的復雜性，通常無法找到對應的解析解，目前基于有限單元法和時間積分方法的數值求解策略占主導地位。

一般來說，多體系統中的柔性部件可以在浮動坐標系、慣性坐標系和共旋坐標系下進行建模。浮動坐標系源于剛體動力學，隨后擴展到小變形問題[3，目前在對柔性多體的建模上仍具有較大的局限性。對于柔性多體建模，慣性坐標系具有廣泛的應用，可以通過幾何精確公式（geometricexactformulation，GEF）或絕對節點坐標公式（absolutenod-alcoordinateformulation，ANCF）來實現。兒何精確公式最初是由SIMO4提出的，在過去幾十年間，一些基于幾何精確公式的先進梁5、板/殼單元被陸續提出。但是，由于在截面大轉動上描述的冗余性限制了幾何精確法在柔性多體系統建模領域的發展。為了解決這些限制，絕對節點坐標法被提出。絕對節點坐標公式最早出現在SHABANA的工作中，與幾何精確法不同的是，絕對節點坐標法不再將有限旋轉作為節點坐標，而是使用絕對位移和全局斜率作為單元坐標。利用單元形狀函數和節點坐標在全局坐標系中定義了有限元上各質點的位置和變形。總之，幾何精確法和絕對節點坐標法都能夠準確地分析柔性結構的大范圍運動。但是，當問題局限于小應變和有限轉動時，慣性坐標系不如共旋坐標系簡便。

共旋坐標系[8-1]來源于連續介質力學把剛體運動從整體運動中分離出來，該概念首先因WEMPNER、BELYTSCHKO等9提出。共旋坐標系的關鍵思想是利用一個局部坐標系從總運動中提取純變形，該坐標系隨每個單元平移和旋轉，但不隨單元變形。由于共旋坐標系中單元的純變形較小，因此可以使用幾何線性有限元公式計算共旋坐標系中的單元剛度矩陣和內力矢量，進一步通過將單元剛度矩陣和內力矢量從共旋坐標系到全局坐標系的一致變換來考慮幾何非線性效應。由于共旋坐標系不依賴于幾何線性單元類型，因此可以使用大量現有的具有魯棒性和高精度的幾何線性單元，并將其擴展到共旋坐標系的幾何非線性分析中。

對于空間離散后系統的動力學響應分析，時間積分方法是一種非常有用的數值工具[12-13]，其中一些方法，比如Newmark方法和后向差分公式，已經在許多商業軟件得到廣泛使用。時間積分方法最初是用來求解由常微分方程（ordinarydifferentialequation，ODE）控制的結構動力學系統。為了追求更高的精度和效率，在經典時間積分方法的基礎上發展了一些更先進的時間積分方法，包括參數法[14-15]、保能量法[16-17]、配點法[18-19]等。由于參數法具有二階精度、耗散可控等優點，且易于實現，因此很自然地被推廣到多體系統領域。目前，General-ized- α 方法[14，20]作為參數方法的一種，已經被廣泛地用于求解多體系統。此外，為了精確地計算轉動變量，學者們還構造了Lie-Generalized ?α 方法[21]由于參數方法的動力學平衡方程僅滿足于廣義的時間點，而不是時間離散點，因此他們的加速度精度只有一階。

為了解決這一問題，學者們提出了復合方法。這類方法首次出現在BANK等[22]的工作中，然后由BATHE等[23]進行概念化。復合方法的思想是將不同的時間積分方法有效地組合在一個時間步長中，同時提高了精度、效率和穩定性。更重要地是，復合方法在離散的時間點上嚴格滿足動力學平衡方程。目前已經發展出一些性能較好的復合方法，包括WEN方法[24]、KIM方法[25]、TTBIF[26]等。在這些復合方法中，文獻[27-28]提出的TTBIF在精度上具有顯著的優勢。

在此背景下，本文為柔性多體系統開發了一種基于共旋坐標系和TTBIF的動力學響應求解策略。數值試驗表明，相比當前流行的基于絕對節點坐標系和Generalized 方法建立的求解策略，本文策略在精度和效率上均具有明顯優勢。

1共旋框架下柔性多體動力學方程

柔性多體動力學建模主要研究梁、板、殼的柔性多體模型，而梁作為最簡單、最常見的一種柔性多體模型在工程領域廣泛存在。在此背景下，本文主要研究共旋框架下基于二維梁單元的柔性多體動力學方程。

無論采取哪種坐標系對柔性多體進行建模，對應的動力學行為都可以利用如下微分代數方程進行描述：

式中， q 為廣義坐標向量，包括平動和轉動變量； F_I 為慣性力； F_E 為彈性力； ?（q，t） 為包含運動和幾何約束的向量； 為向量 ?（q，t） 關于 q 求導的矩陣；為拉格朗日乘子； Q（t） 為外部激勵項。采用不同坐標系，式（1）中 F_I 和 F_E 的構造方式是不同的，下面介紹共旋框架下基于一般二維梁單元柔性多體系統動力學方程的建立。

共旋坐標系的核心是利用一個局部坐標系從總運動中提取純變形。如圖1所示，為了描述一個梁模型在平面的運動，首先對每個單元建立一個隨單元進行平動和轉動的局部坐標系 x_LO_Ly_L ，并規定：坐標原點 O_L 位于梁單元左節點1； x_L 沿著梁單元兩節點連線，由節點1指向節點2。 β 和 β₀ 為當前構型和初始構型下 x_L 軸與 x_G 軸的夾角； α=β-β₀ 。

從式（33）和（39）中可以看到，系數 和 的引入可以確保雅可比矩陣的所有元素都是 O（Δt⁰） ，從而避免了由于時間步長取值過小導致的雅可比矩陣奇異性問題。此外，為了方便讀者執行本文提出的策略，圖3提供了Co-TTBIF的計算流程。

圖3Co-TTBIF的計算流程 Fig.3 Step-by-step procedure of Co-TTBIF

3數值算例

本節利用2個算例來驗證本文策略（記為Co-TTBIF）在柔性多體系統動力學仿真中的優勢。作為對比，基于絕對節點坐標單元和共旋坐標單元與Generalized-α方法的組合求解策略也被考慮在本節，分別記為ANCF-Gα和 Co^-Gα 。考慮到TTBIF在一個時間步長內包含3個分步，而Generalized 方法是單分步的。因此，在所有算例中，本文策略的時間步長取為 ANCF-Gα 和 Co^-Gα 的3倍。此外，為了過濾柔性體空間離散引入的虛假振蕩，所有策略采用 ρ_∞=0 。

3.1自由下落單擺模型

首先考慮受重力作用下自由下落的柔性單擺模型[29]，如圖4所示。物理參數假設為：擺長 L=1.2m 、橫截面積 A=a×b=0.015×0.01=0.00015m²， 楊氏模量 E=10MPa 、泊松比 u=0.3， 密度 ρ= 5540kg/m³ 。在本算例中，考慮3種求解策略，分別是：（1）4個考慮剪切變形的二維ANCF梁單元（12個自由度）與Generalized- α 方法；（2）4個考慮剪切變形的二維共旋梁單元（6個自由度）與General-ized- 方法；（3）4個考慮剪切變形的二維共旋梁單元與TTBIF。假設 Δt（ANCF-Gα）=10^-4 S、 Δt （Co-Gα）=10^-4 s和 Δt（Co^-TTBIF）=3×10^-4， S。參考解由取更小 Δt 的ANCF- ΔGα 提供。

圖4共旋坐標系下自由下落的柔性梁模型 Fig.4Free falling flexible beam model in the corotational coordinate system

對于該保守系統，在當前初始條件下，系統的總能量理論上恒為0。圖5～7提供了該模型的能量和位移結果。從圖5中可以看到，ANCF- Gα 、 Co^-Gα 和Co-TTBIF都較好地保守系統的能量。進一步，從圖6和7中可以看到：對于具有大轉動的柔性多體系統，ANCF σ_Gα Co^-Gα 和Co-TTBIF對動力學響應都給出了準確的預測。

為了展示Co-TTBIF在柔性多體系統動力學仿真中的優勢，圖8對ANCF- σ_Gα Co-Ga和Co-TTBIF的Newton迭代次數和CPU時間進行了比較。從中發現：（1）相同單元數下，共旋單元用更少的自由度數實現了與ANCF單元相同的精度；（2）因為TT-BIF的動力學平衡方程嚴格滿足時間離散點，而Generalized 方法的動力學平衡方程只在廣義的時

間點滿足，因此TTBIF更容易收斂，Newton迭代次數更少，計算量更小；（3）相同精度下， Co TTBIF具有顯著的效率優勢，相比當前流行的 ANCF-Gα ，大約可以節約 95% 的計算量。

3.2受迫激勵下的柔性懸臂梁模型

第二個算例考慮一個集中力作用下的柔性懸臂梁模型[30]，如圖9所示，物理參數為：梁長 L₀=10m 橫截面積 A=a×b=0.25×0.5=0.125m²， 楊氏模量 E= 210GPa 密度 ρ=7850kg/m³ 。在懸臂梁的自由端受到正弦激勵 P=P₀sin（ωt） ，其中 P₀=10MN， ω= 50rad/s_o 該懸臂梁模型利用20個共旋梁單元進行空間離散。 Co^-Gα 和Co-TTBIF的時間步長分別取為10^-4 s和 3×10^-4 s，參考解由取更小 ?Δt 的Co-Ga提供。

圖10繪制了 Co^-Gα 和Co-TTBIF在自由端的位移隨時間的演化規律，Ref表示參考解。圖中給出了兩種求解策略的平均絕對誤差（averageabsoluteerrors，AAE）。可以發現：這兩種求解策略都能較好地預測動力學響應，但Co-TTBIF的精度更高一些。圖11比較了 Co^-Gα 和Co-TTBIF在自由端的速度隨時間的演化規律，可以發現由于柔性體的空間離散導致速度中存在數值振蕩，特別是在轉角的計算中。隨著計算時間增加，相比Co-TTBIF，Co^-Gα 的數值振蕩更加劇烈。

進一步，為了討論人工阻尼量對柔性多體系統計算結果的影響，圖12對比了Co-TTBIF和 Co^-Gα 在不同 ρ_∞ 取值下自由端速度的結果。數值結果顯示：隨著耗散程度的逐漸降低，Co-TTBIF和 Co^-Gα 的數值震蕩越來越劇烈。因此，在柔性連續體的動力學仿真中，一定程度的數值耗散量是必備的，其可以有效地過濾掉不需要的虛假信息。

Co-Ga Co-TTBIF 100 50 -50 w 100 （s.ppi）/0 100 Mmmn 5000050° 50 w 100 0 0.2 0.4 0.6 0.8 1.0 時間/s

為了進一步比較 Co^-Gα 和Co-TTBIF的效率，圖13給出了兩種求解策略的Newton迭代次數和CPU時間。結果表明：在Co-TTBIF比 Co^-Gα 略有精度優勢時，其計算量僅為 Co-Gα 的1/8。Co-TTBIF比 Co^-Gα 針對柔性多體采用了完全相同的建模方式，但由于TTBIF的動力學平衡方程嚴格滿足時間離散點，導致其Newton迭代次數更少，計算量顯著低于Generalized 方法。

綜上所述，本節通過對兩個具有大旋轉和大變形的柔性多體系統的動力學分析發現：與當前柔性多體系統動力學分析領域廣泛使用的求解策略ANCF-Gα相比，本文提出的Co-TTBIF在求解同類動力學系統時具有可觀的精度和效率優勢。

4結論

本文針對柔性多體系統提出了一種共旋框架下基于TTBIF的動力學分析策略，記為Co-TTBIF。在本文策略中，柔性體的動力學模型被建立在共旋坐標系下，因此任意先進的幾何線性單元均可以被使用，進而降低動力學方程的非線性程度。隨后，利用優化的復合方法TTBIF對建立的動力學模型進行時域分析，在有效地過濾由于柔性體空間離散引入的數值振蕩的同時可以高精度地計算重要的低頻響應。與已有工作相比，本文的創新之處在于利用共旋理論在保證動力學模型精度的前提下顯著降低模型自由度，同時利用復合型方法求解動力學模型，在保證計算精度和人工可控阻尼的前提下有效降低迭代次數，從而實現柔性多體系統動力學響應的精確和快速仿真。

利用兩個經典的柔性多體動力學模型分析了本文策略Co-TTBIF在處理這類具有強非線性和大范圍運動的動力學系統的表現。計算結果顯示：與當前流行的絕對坐標框架下基于Generalized- 方法的動力學求解策略（ANCF-Gα）相比，在精度相近的前提下，本文策略具有顯著的效率優勢。

參考文獻：

[1] 杜超凡，周曉婷，章定國，等.無網格法在中心剛體-旋 轉柔性梁系統動力學分析中的應用[J].振動工程學 報，2022，35（1）：178-187. DU Chaofan，ZHOU Xiaoting，ZHANG Dingguo，et al.Meshless methods applied to the dynamic analysis forhub-beamsystem[J].JournalofVibrationEngineer ing，2022，35（1）：178-187.

[2] 孫加亮，孫佳昊，金棟平，等.充氣自旋結構的動力學特 性研究[J].振動工程學報，2023，36（3）：662-670. SUN Jialiang，SUN Jiahao， JIN Dongping，et al. Study ondynamic characteristics of inflatable spin structure[J]. Journal ofVibrationEngineering，2023，36（3）：662-670.

[3] DEVEUBEKEBF.Thedynamics of flexible bodies [J].International Journal of Engineering Science， 1976，14（10）：895-913.

[4] SIMOJC.A finite strainbeam formulation.Thethreedimensional dynamic problem.Part I[J].Computer Methods in Applied Mechanics and Engineering，1985， 49（1）：55-70.

[5]MOURID E，MAMOURI S，KOULIR， et al. An efficient implicit time stepping scheme for nonlinear dynamics analysis of 3D geometrically exact beam using exponential mapping[J].Mechanics of Advanced Materials and Structures，2022，29（28）：7189-7203.

[6]ZHANG T，LIU C，TANG H Y.Geometrically exact shell with drilling rotations formulated on the special Euclidean group SE（3）[J]. International Journal for Numerical Methods in Engineering，2021，122（18）： 4886- 4921.

[7]SHABANA A A.Flexible multibody dynamics：review of past and recent developments[J].Multibody System Dynamics，1997，1（2）：189-222.

[8]WEMPNER G.Finite elements， finite rotations and small strains of flexible shels[J]. International Journal of Solids and Structures，1969，5（2）：117-153.

[9]BELYTSCHKO T， SCHWER L， KLEIN M J. Large displacement， transient analysis of space frames[J]. International Journal for Numerical Methods in Engineering，1977，11（1）：65-84.

[10]WANG G，QI Z H，XU JS.A high-precision corotational formulation of 3D beam elements for dynamic analysis of flexible multibody systems[J].Computer Methods in Applied Mechanics and Engineering，2020， 360： 112701.

[11] SHI JB，LIU Z Y，HONG J Z. Multibody dynamic analysis using a rotation-free shell element with corotational frame[J].Acta Mechanica Sinica，2O18，34（4）： 769-780.

[12］邢譽峰，張慧敏，季奕.動力學常微分方程的時間積分 方法[M].北京：科學出版社，2022.

[13]HUGHES TJR. The Finite Element Method：Linear Staticand DynamicFiniteElementAnalysis[M].New Jersey ： Prentice Hall Press，1987.

[14] CHUNG J，HULBERT G M.A time integration algo rithm for structural dynamics with improved numerical dissipation：the generalized- α method[J]. Journal of Applied Mechanics，1993，60（2）：371-375.

[15] ZHANG H M， XING Y F.A threeparameter singlestep time integration method for structural dynamic analysis[J].Acta Mechanica Sinica，2019，35（1）： 112-128.

[16] ZHANG H M， XING Y F， JI Y. An energy-conservingand decaying time integration method for general nonlinear dynamics[J]. International Journal for Numerical Methods in Engineering，2020，121（5）：925-944.

[17]LIU C Y，LI JY，YANG ZQ，et al. Two high-order energy-preserving and symmetric Gauss collocation integrators for solving the hyperbolic Hamiltonian systems [J].Mathematics and Computers in Simulation，2023， 205：19-32.

[18]KIMW，REDDYJN.A new familyofhigherorder time integration algorithms for the analysis of structural dynamics[J]. Journal ofApplied Mechanics，2O17，84 （7）：071008.

[19]FUNG T C.Unconditionally stable collocation algo rithmsforsecond order initial valueproblems[J].Jour nal of Sound andVibration，2001，247（2）：343-365.

[20]邵慧萍，蔡承文.結構動力學方程數值積分的三參數 算法[J].應用力學學報，1988，5（4）：76-81. SHAO Huiping，CAI Chengwen.A three parameters algorithm fornumerical integration of structural dynamic equations[J]. Chinese Journal of Applied Mechanics， 1988，5（4）：76-81.

[21]BRULSO，CARDONA A，ARNOLD M.Lie group generalized- α timeintegration of constrained flexible multibody systems[J].Mechanism and Machine Theory，2012，48：121-137.

[22]BANKRE，COUGHRANWM，FICHTNERW，et al.Transient simulation of silicon devices and circuits [J].IEEE Transactions on Computer-Aided Design of Integrated Circuitsand Systems，1985，4（4）： 436-451.

[23]BATHEKJ，BAIG MMI.On a composite implicit time integration procedure for nonlinear dynamics[J]. Computersamp;.Structures，2005，83（31-32）：2513-2524.

[24]WENWB，LEIM，WANGP，etal.An efficient hybrid implicit time integration method with high accuracy for linear and nonlinear dynamics[J]. European Journal ofMechanics-A，Solids，2023，97：104811.

[25]KIM W.An improved implicit method with dissipation control capability： the simple generalized composite timeintegration algorithm[J].Applied Mathematical Modelling，2020，81：910-930.

[26]JIY，XING YF.An optimized three-sub-step compositetime integration method with controllable numerical dissipation[J].Computers amp;Structures，2O2O，231： 106210.

[27]LIUTH，HUANGFL，WENWB，etal.Further insights ofacomposite implicit timeintegration scheme and its performance on linear seismic response analysis [J].Engineering Structures，2O21，241：112490.

[28]JI Y，XING Y F.A three-sub-step composite method fortheanalysis of rigid body rotations with Euler parameters[J].NonlinearDynamics，2023，111（15）：14309- 14333.

[29]BERZERI M，SHABANA A A.Developmentof simple models for the elastic forces in the absolute nodal coordinate formulation[J].Journal of Sound and Vibration，2000，235（4）：539-565.

[30]LETN，BATTINIJM，HJIAJM.Efficientformulation fordynamicsof corotational 2D beams[J].Compu tationalMechanics，201l，48（2）：153-161.