顯式和隱式力學計算方法對比研究

潘科琪

（上海核工程研究設計院，上海　200233）

顯式和隱式力學計算方法對比研究

潘科琪

（上海核工程研究設計院，上海200233）

0　引言

柔性多體系統動力響應的計算可以通過建立多體系統的封閉的拉格朗日第一類方程，對于剛柔耦合系統，初始時刻的條件是已知，因此，動力學方程的計算歸結于求解數學上的微分代數方程初值問題。

動力學方程的求解主要有顯式和隱式兩類算法。常用的顯式算法主要有龍格庫塔 （Runge-kutta Method）和中心差分法［1］，為了保證求解的精度，顯式算法要使積分步長調整到很小，導致CPU的計算時間顯著增加，尤其是存在幾何非線性變形時，顯式算法的計算效率較低，但是其程序編寫較為簡單、方便。與顯式算法不同，隱式算法在計算時并不跟蹤那些對計算結果沒有較大影響的高頻振蕩成分，正因如此，步長稍微大一些時也不會對結果準確性有較大的影響，能夠顯著提高計算效率。常用的隱式算法有紐馬克（Newmark）［2］、威爾遜（Wilson）-θ法［3］、HHT（Hilber-Hughes-Taylor）法［4］、廣義（Generalized-α）法［5］。對于邊界條件明確、規則形狀的梁或者板構件的柔性多體系統動力學方程，可以優先考慮模態縮減法離散動力學方程，以減少節點坐標的數量，減少動力學方程的維數。此時，傳統的龍格庫塔法還是適用的。若采用的MATLAB軟件編程計算，可以調用函數ODE直接求解。

本文采用兩種計算方法，分別求解單擺和曲柄滑塊多體系統方程，通過MATLAB仿真編程計算，分別比較了兩種計算方法求解線性和非線性動力學方程的計算效率。

1　動力學方程

柔性多體系統的封閉的拉格朗日第一類方程

如（1），其中，M、Qe及Qd分別是多體系統的廣義質量陣和廣義力陣，Φe約束方程的雅可比陣，λ是約束方程的拉格朗日乘子列陣。上式是典型的微分代數方程（DAEs）。

2　計算方法

2.1紐馬克方法計算步驟

（1）將公式（1）變為如下形式

根據紐馬克教授于1959年提出來的紐馬克方法［6］，在計算中首先作了如下假設：

其中，h是時間步長，γ和β分別是積分常數。

（2）將方程（2）修改為如下形式

2.2龍哥庫塔方法的計算步驟

對于邊界條件明確、規則形狀的梁或者板構件的柔性多體系統動力學方程，可以優先考慮模態縮減法離散動力學方程，以減少節點坐標的數量，減少動力學方程的維數，此時，傳統的龍格庫塔法還是適用的。也可采用的MATLAB軟件中的函數ODE命令直接求解，其求解的主要步驟為：

在MATLAB軟件中，積分過程中的函數dy=f（t，q）可以表示為：

上面的求解過程中沒有涉及到迭代，流程簡單，對計算方法的高效性要求會有所降低。但是由于在大變形的條件下，模態階數選取的多少對計算精度的影響很大，仿真程序的通用性和精確性較差。

3　仿真算例

3.1單擺模型

單擺梁的材料和幾何參數為：密度=27667kg/m3，彈性模量E=6.8952×1010N/m2，橫截面面積 A=8×10-5m2，橫截面慣性矩 I=1.06667×10-8m2。定曲率曲梁，其坐標系上的初始形狀可以用函數描述為：

其中，αs為曲梁的中心角。

表1所示為用不同計算方法曲梁單擺的剛-柔耦合動力學方程所需計算時間。從表中可以看出，在仿真計算曲梁單擺線性模型的動力學方程時，龍格庫塔法耗時3個小時，紐馬克方法耗時76.93秒，增量法僅需要18.03秒，龍格庫塔法計算效率最低，增量法的計算效率最高；在計算非線性模型的動力學方程時，龍格庫塔法的仿真耗時顯著增多。

圖1　曲梁單擺

表1　曲梁單擺的計算時間比較

3.2曲柄滑塊多體系統模型

為了進一步比較各種計算方法在處理有鉸約束的多體系統問題的計算效率，本文對曲梁多體系統的仿真時間做了比較。表2比較了變曲率連桿的曲柄滑塊機構的動力學問題所耗費的計算機時間。

圖2所示為連桿為曲梁的曲柄滑塊機構。曲柄的長度為1m，橫截面面積A=0.02×0.02m2，變?曲率曲梁的形狀與公式（14）中的相同，其中參數d=1m，l=4m，滑塊的質量0.5kg。施加在曲柄的驅動約束為：

其中，ω0=π是曲柄穩態運動階段的角速度，ts=0.5為加速時間。

從表2中可以看出，在仿真計算曲柄滑塊的線性模型時，龍格庫塔法耗時要10多個小時，紐馬克方法耗時957.50秒，但是在計算非線性模型時，龍格庫塔法耗時要100多個小時，紐馬克方法耗時120多個小時，在處理多柔性體系統的非線性動力學方程時，由于鉸約束方程的增加，紐馬克方法相對于龍格庫塔法仍然具有顯著優勢。

圖2　曲柄滑塊結構

表2　曲柄滑塊多體系統的計算時間比較

4　結語

本文分別采用龍哥庫塔法和紐馬克方法，采用MATLAB編程求解單擺和曲柄滑塊多體系統的線性和非線性的剛-柔耦合動力學方程，在相同的結果精度下，紐馬克方法的計算效率要明顯高于龍哥庫塔法，因此，在商用或者自編程序仿真計算時，采用隱式算法能夠顯著提高計算效率，節省CPU的計算時間。

［1］凌復華，殷學剛，何冶奇.常微分方程數值方法及其在力學中的應用，重慶大學出版社，1990.

［2］Newmark N.M.A method of computation for structural dynamics.Journal of the Engineering Mechanics Division ASCE，1959，85（3）：67-94.

［3］Wilson E L.A computer program for the dynamics stress analysis of underground structures，SESM Report No.68-1.Division of structural Engineering and Structural Mchanics.University of California，Berkeley，CA.

［4］Hilber H M，Hughes T J R，Taylor R L.Improved numerical dissipation for time integration algorithms in structural dynamics.Earthquake Engineering and Structural Dynamics，1977，5：283-292.

［5］Chung J，Hulbert G M.A time integration algorithm for structural dynamics with improved numerical dissipation：The Generalized-Method.Journal of Applied Mechanics，1993，60，371-375.

［6］王文亮.結構動力學.復旦大學出版社，1993.

Kinetic Equations；Explicit Algorithm；Implicit Algorithms

Comparative Research on the Explicit and Implicit Mechanical Calculation Method

PAN Ke-qi
（Shanghai Nuclear Engineering Research&Design Institute，Shanghai 200233）

1007-1423（2015）20-0020-04

10.3969/j.issn.1007-1423.2015.20.005

潘科琪（1984-），女，遼寧開原人，博士研究生，工程師，研究方向為反應堆結構力學2015-05-05

2015-07-03

闡述求解剛-柔耦合動力學方程的兩種計算方法及其研究進展。基于顯式算法即龍哥庫塔法和隱式算法方法即紐馬克方法結合牛頓迭代法，分別求解單擺和曲柄滑塊系統的線性和非線性的動力學方程。通過結果對比研究發現，在保證相同的計算精度下，隱式算法求解線性和非線性動力學方程的計算效率都要高于顯式算法。指出在工程計算建議采用隱式算法進行求解。

動力學方程；顯式算法；隱式算法

Shows two kinds of numerical computation method for rigid-flexible dynamics equation and its developments.Based on the explicit numerical method，respectively solves i.e.Runge Kutta method and implicit numerical method，i.e.Newmark method，linear and nonlinear dynamics equations for curved beam pendulum and slider-crank mechanism.In condition of the same computation accuracy，results comparison show that implicit numerical method computation efficiency is obviously higher than the explicit numerical method no matter for solving linear nor the nonlinear dynamics equation.So，the conclusion that implicit numerical method，i.e.Newmark method is recommended for calculating the mechanics problem in engineering.