顯式和隱式力學(xué)計(jì)算方法對比研究

潘科琪

（上海核工程研究設(shè)計(jì)院，上海　200233）

顯式和隱式力學(xué)計(jì)算方法對比研究

潘科琪

（上海核工程研究設(shè)計(jì)院，上海200233）

0　引言

柔性多體系統(tǒng)動(dòng)力響應(yīng)的計(jì)算可以通過建立多體系統(tǒng)的封閉的拉格朗日第一類方程，對于剛?cè)狁詈舷到y(tǒng)，初始時(shí)刻的條件是已知，因此，動(dòng)力學(xué)方程的計(jì)算歸結(jié)于求解數(shù)學(xué)上的微分代數(shù)方程初值問題。

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

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

1　動(dòng)力學(xué)方程

柔性多體系統(tǒng)的封閉的拉格朗日第一類方程

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

2　計(jì)算方法

2.1紐馬克方法計(jì)算步驟

（1）將公式（1）變?yōu)槿缦滦问?/p>

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

其中，h是時(shí)間步長，γ和β分別是積分常數(shù)。

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

2.2龍哥庫塔方法的計(jì)算步驟

對于邊界條件明確、規(guī)則形狀的梁或者板構(gòu)件的柔性多體系統(tǒng)動(dòng)力學(xué)方程，可以優(yōu)先考慮模態(tài)縮減法離散動(dòng)力學(xué)方程，以減少節(jié)點(diǎn)坐標(biāo)的數(shù)量，減少動(dòng)力學(xué)方程的維數(shù)，此時(shí)，傳統(tǒng)的龍格庫塔法還是適用的。也可采用的MATLAB軟件中的函數(shù)ODE命令直接求解，其求解的主要步驟為：

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

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

3　仿真算例

3.1單擺模型

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

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

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

圖1　曲梁單擺

表1　曲梁單擺的計(jì)算時(shí)間比較

3.2曲柄滑塊多體系統(tǒng)模型

為了進(jìn)一步比較各種計(jì)算方法在處理有鉸約束的多體系統(tǒng)問題的計(jì)算效率，本文對曲梁多體系統(tǒng)的仿真時(shí)間做了比較。表2比較了變曲率連桿的曲柄滑塊機(jī)構(gòu)的動(dòng)力學(xué)問題所耗費(fèi)的計(jì)算機(jī)時(shí)間。

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

其中，ω0=π是曲柄穩(wěn)態(tài)運(yùn)動(dòng)階段的角速度，ts=0.5為加速時(shí)間。

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

圖2　曲柄滑塊結(jié)構(gòu)

表2　曲柄滑塊多體系統(tǒng)的計(jì)算時(shí)間比較

4　結(jié)語

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

［1］凌復(fù)華，殷學(xué)剛，何冶奇.常微分方程數(shù)值方法及其在力學(xué)中的應(yīng)用，重慶大學(xué)出版社，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］王文亮.結(jié)構(gòu)動(dòng)力學(xué).復(fù)旦大學(xué)出版社，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-），女，遼寧開原人，博士研究生，工程師，研究方向?yàn)榉磻?yīng)堆結(jié)構(gòu)力學(xué)2015-05-05

2015-07-03

闡述求解剛-柔耦合動(dòng)力學(xué)方程的兩種計(jì)算方法及其研究進(jìn)展?；陲@式算法即龍哥庫塔法和隱式算法方法即紐馬克方法結(jié)合牛頓迭代法，分別求解單擺和曲柄滑塊系統(tǒng)的線性和非線性的動(dòng)力學(xué)方程。通過結(jié)果對比研究發(fā)現(xiàn)，在保證相同的計(jì)算精度下，隱式算法求解線性和非線性動(dòng)力學(xué)方程的計(jì)算效率都要高于顯式算法。指出在工程計(jì)算建議采用隱式算法進(jìn)行求解。

動(dòng)力學(xué)方程；顯式算法；隱式算法

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.