潘科琪

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

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

潘科琪

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

0　引言

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

動(dòng)力學(xué)方程的求解主要有顯式和隱式兩類(lèi)算法。常用的顯式算法主要有龍格庫(kù)塔 （Runge-kutta Method）和中心差分法［1］，為了保證求解的精度，顯式算法要使積分步長(zhǎng)調(diào)整到很小，導(dǎo)致CPU的計(jì)算時(shí)間顯著增加，尤其是存在幾何非線性變形時(shí)，顯式算法的計(jì)算效率較低，但是其程序編寫(xiě)較為簡(jiǎn)單、方便。與顯式算法不同，隱式算法在計(jì)算時(shí)并不跟蹤那些對(duì)計(jì)算結(jié)果沒(méi)有較大影響的高頻振蕩成分，正因如此，步長(zhǎng)稍微大一些時(shí)也不會(huì)對(duì)結(jié)果準(zhǔn)確性有較大的影響，能夠顯著提高計(jì)算效率。常用的隱式算法有紐馬克（Newmark）［2］、威爾遜（Wilson）-θ法［3］、HHT（Hilber-Hughes-Taylor）法［4］、廣義（Generalized-α）法［5］。對(duì)于邊界條件明確、規(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)的龍格庫(kù)塔法還是適用的。若采用的MATLAB軟件編程計(jì)算，可以調(diào)用函數(shù)ODE直接求解。

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

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

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

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

2　計(jì)算方法

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

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

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

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

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

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

對(duì)于邊界條件明確、規(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)的龍格庫(kù)塔法還是適用的。也可采用的MATLAB軟件中的函數(shù)ODE命令直接求解，其求解的主要步驟為：

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

上面的求解過(guò)程中沒(méi)有涉及到迭代，流程簡(jiǎn)單，對(duì)計(jì)算方法的高效性要求會(huì)有所降低。但是由于在大變形的條件下，模態(tài)階數(shù)選取的多少對(duì)計(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ì)算方法曲梁?jiǎn)螖[的剛-柔耦合動(dòng)力學(xué)方程所需計(jì)算時(shí)間。從表中可以看出，在仿真計(jì)算曲梁?jiǎn)螖[線性模型的動(dòng)力學(xué)方程時(shí)，龍格庫(kù)塔法耗時(shí)3個(gè)小時(shí)，紐馬克方法耗時(shí)76.93秒，增量法僅需要18.03秒，龍格庫(kù)塔法計(jì)算效率最低，增量法的計(jì)算效率最高；在計(jì)算非線性模型的動(dòng)力學(xué)方程時(shí)，龍格庫(kù)塔法的仿真耗時(shí)顯著增多。

圖1　曲梁?jiǎn)螖[

表1　曲梁?jiǎn)螖[的計(jì)算時(shí)間比較

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

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

圖2所示為連桿為曲梁的曲柄滑塊機(jī)構(gòu)。曲柄的長(zhǎng)度為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í)，龍格庫(kù)塔法耗時(shí)要10多個(gè)小時(shí)，紐馬克方法耗時(shí)957.50秒，但是在計(jì)算非線性模型時(shí)，龍格庫(kù)塔法耗時(shí)要100多個(gè)小時(shí)，紐馬克方法耗時(shí)120多個(gè)小時(shí)，在處理多柔性體系統(tǒng)的非線性動(dòng)力學(xué)方程時(shí)，由于鉸約束方程的增加，紐馬克方法相對(duì)于龍格庫(kù)塔法仍然具有顯著優(yōu)勢(shì)。

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

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

4　結(jié)語(yǔ)

本文分別采用龍哥庫(kù)塔法和紐馬克方法，采用MATLAB編程求解單擺和曲柄滑塊多體系統(tǒng)的線性和非線性的剛-柔耦合動(dòng)力學(xué)方程，在相同的結(jié)果精度下，紐馬克方法的計(jì)算效率要明顯高于龍哥庫(kù)塔法，因此，在商用或者自編程序仿真計(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-），女，遼寧開(kāi)原人，博士研究生，工程師，研究方向?yàn)榉磻?yīng)堆結(jié)構(gòu)力學(xué)2015-05-05

2015-07-03

闡述求解剛-柔耦合動(dòng)力學(xué)方程的兩種計(jì)算方法及其研究進(jìn)展。基于顯式算法即龍哥庫(kù)塔法和隱式算法方法即紐馬克方法結(jié)合牛頓迭代法，分別求解單擺和曲柄滑塊系統(tǒng)的線性和非線性的動(dòng)力學(xué)方程。通過(guò)結(jié)果對(duì)比研究發(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.