數值計算求解非線性動力方程

郭 沖

（山西省公路局長治分局勘測設計所，山西 長治 046000）

1 引言

如果外荷載或者結構的剛度、質量、阻尼等任何一項隨時間變化，那么這樣的結構體系就是非線性的，這種結構體系的動力方程也是非線性的。對于此類非線性動力方程，其解析解往往較難由理論推導得到，這類問題可以通過數值計算方法對動力響應方程進行求解。常用的數值計算方法包括激勵線性插值法、中心差分法、Newmark的平均加速度法和線性加速度法等。

2 各數值算法的基本原理

2.1 激勵線性插值法

對于線性體系，通過在每個時間間隔里對激勵進行插值，并利用線性體系響應解析解，能夠推導出一種有效的數值計算方法。只要保證插值的時間間隔較短，就可以得到令人滿意的求解結果。

對于欠阻尼體系（ξ＜1），有非線性動力方程

在時間間隔0≤τ≤△ti內，反應u（τ）由3部分組成：

（1）τ=0時刻的初位移ui和初速度u觶i引起的自由振動響應；

（2）零初始條件下對階躍力pi的反應；

（3）零初始條件下對斜坡力（△pi/△ti）τ的反應。

經推導可得出第i+1個時間步的位移與速度公式如下：

ω：外荷載頻率；

ωn：結構固有頻率。

給定初始條件，代入式（1）即可得到結構在任意時刻的響應。

2.2 中心差分法

該方法是基于對位移時間導數的有限差分近似進行的。取固定的時間步長△t，則時刻i的速度和加速度的中心差分格式為：

將上式的速度加速度代入動力方程中得：

為了計算ui+1，需要知道位移ui和ui-1。若已知初始位移u0，則u-1可表示為：

至此，所有參數均可以通過初始條件與公式進行遞推求得，則任意時刻的結構響應可以通過數值計算得到。

2.3 Newmark法

N.M.Newmark發展了一類時間步進法，它們基于下面的公式：

當 γ=0.5，β=0.25時，上式為平均加速度法；當 γ=0.5， β=1/6時，上式為線性加速度法。對式（2）以增量格式重寫如下：

3 程序實現

上文對激勵線性插值法、中心差分法與Newmark法的基本原理進行了介紹，借助數學計算軟件Matlab，對這幾種方法編制了相應的計算程序。現將Matlab中編制的函數展示如下。

3.1 激勵線性插值法

function y=excitation（m，k，Tn，e，dt，t，F，te）

a=[0：dt：t]；n=length（a）；

for j=1：n

if a（j）

p（j）=F.*sin（pi.*a（j）./te）；

else

p（j）=0；

end

end

Wn=2.*pi./Tn；Wd=Wn.*（1-e.^2）.^0.5；

Si=sin（Wd.*dt）；Co=cos（Wd.*dt）；

E=exp（-e.*Wn.*dt）；

A=E.*（（e./（1-e.^2）.^0.5）.*Si+Co）；

B=E.*（Si./Wd）；

C=（2.*e./（Wn*dt）+E.*（（（1-2.*e.^2）./（Wd.*dt）-（e./（1-e.^2）.^0.5））.*Si-（1+2*e./（Wn*dt））.*Co））/k；

D=（1-2.*e./（Wn*dt）+E.*（（2.*e^2-1）./（Wd.*dt）.*Si+2.*e.*Co./（Wn.*dt）））./k；

At=-E.*（Wn./（（1-e.^2）.^0.5）.*Si）；Bt=E.*（Co-（e./（1-e.^2）.^0.5）.*Si）；

Ct=（-1./dt+E.*（（Wn./（（1-e.^2）.^0.5）+（e./（dt.*（1-e.^2）.^0.5））.*Si+Co./dt））/k；

Dt=（1-E.*（（e./（1-e.^2）.^0.5）.*Si+Co））./（k.*dt）；u（1）=0；ut（1）=0；

for i=2：n

u（i）=A.*u（i-1）+B.*ut（i-1）+C.*p（i-1）+D.*p（i）；

ut（i）=At.*u（i-1）+Bt.*ut（i-1）+Ct.*p（i-1）+Dt.*p（i）；

end

x1=[0：dt：t]；m=length（x1）；

y1=zeros（m，1）；y2=zeros（m，1）；

for i=1：n

y1（i）=u（i）；

y2（i）=ut（i）；

end

3.2 中心差分法

function difference（m，k，c，u0，ut0，dt，t，te，Tn）

if dt<（Tn./pi）

u（2）=u0；ut（2）=ut0；x=[0：dt：t]；n=length（x）；

for j=2：n+1

if x（j-1）

p（j）=10.*sin（pi.*x（j-1）./te）；

else

p（j）=0；

end

end

utt（2）=（p（2）-c.*ut（2）-k.*u（2））./m；

u（1）=u（2）-（dt）.*ut（2）+（dt）.^2./2.*utt（2）；

kt=m./（dt）^2+c./（2.*dt）；

a=m./（dt）^2-c./（2*dt）；b=k-2.*m./（dt）^2；

for i=2：n+1

pt（i）=p（i）-a.*u（i-1）-b.*u（i）；

u（i+1）=pt（i）./kt；

end

x1=[0：dt：t]；y1=zeros（n，1）；

for j=1：n

y1（j）=u（j+1）；end

3.3 Newmark法（線性加速度法）

function liner（m，k，c，u0，ut0，dt，t，te）

bet=1/6；gama=1/2；u（1）=u0；ut（1）=ut0；

x=[0：dt：t]；n=length（x）；

for j=1：n

if x（j）

p（j）=10.*sin（pi.*x（j）./te）；

else

p（j）=0；

end

end

utt（1）=（p（1）-c.*u（1）-k.*u（1））./m

kt=k+（gama./bet）./dt.*c+（1./bet）./（dt）^2.*m

a=1./（bet.*dt）.*m+（gama./bet）.*c

b=1./（2.*bet）.*m+dt.*（gama./（2.*bet）-1）.*c

for i=1：n-1

dp（i）=p（i+1）-p（i）；

dpt（i）=dp（i）+a.*ut（i）+b.*utt（i）；

du（i）=dpt（i）./kt；

dut（i）=（gama./bet）.*du（i）./dt-（gama./bet）.*ut（i）+dt.*（1-（gama./bet./2））.*utt（i）；

dutt（i）=（1./bet）./（dt）^2.*du（i）-（1./bet）./dt.*ut（i）-（1./bet./2）.*utt（i）；

u（i+1）=u（i）+du（i）；

ut（i+1）=ut（i）+dut（i）；

utt（i+1）=utt（i）+dutt（i）；

end

x1=[0：dt：t]；m=length（x1）；y2=zeros（m，1）；

圖1 理論解與數值計算結果對比

for i=1：n

y2（i）=u（i）；end

4 求解實例

計算實例如下：單自由度體系具有如下特性：m=0.253 3千磅力·秒2/英寸、k=10千磅力/英寸、Tn=1 s（ωn=6.283弧度/s），ξ=0.05。試確定體系在正弦脈沖力p（t）=10sin（πt/0.6）作用下1s內的反應（△t=0.05 s），圖1給出了理論解與數值計算結果對比。

為驗證求解結果的正確性，首先給出該體系在正弦脈沖作用時間1 s內的理論解：

[1]唐友剛.高等結構動力學[M].天津：天津大學出版社，2002.

[2]呂同富，康兆敏，方秀男.數值計算方法[M]，北京：清華大學出版社，2009.

[3]邢棉.MATLAB數值計算在高等數學教學中的應用探討[J].中國科教創新導刊，2010.

[4]李初曄，王增新.結構動力學方程的顯式與隱式數值計算[J].航空計算科學，2010，01.