非線性(時滯)動力學問題的一種新暫態(tài)時程積分解法

顏廷浩, 任傳波, 周繼磊， 吳凱偉， 曹軍帥， 孫志釧

(山東理工大學交通與車輛工程學院, 淄博 255000)

在時域中,數(shù)值解法可直觀地確定非線性系統(tǒng)的運動狀態(tài)(位移、速度、加速度)與系統(tǒng)參數(shù)和初始條件的關系,對于線性和非線性系統(tǒng),已提出大量數(shù)值計算方法[1-4]。對于一般的動力學系統(tǒng),其動力學微分方程可以寫為時間離散化的方程形式,這種方程的數(shù)值解越來越受到關注。

時程積分方法能夠有效地求解微分方程的數(shù)值解,Subbaraj等[5]和Dokainish等[6]研究了直接時程積分在結構動力學系統(tǒng)中的顯式和隱式計算方法,證明這種逐步積分方法是一種強有力的計算手段。Katsikadelis[7]提出了一種直接時間積分方法,用于求解結構線性和非線性多自由度系統(tǒng)動態(tài)響應的運動方程。Chen等[8]對非線性結構直接積分算法的穩(wěn)定性進行了相關的研究。Keierleber等[9]研究了高階隱式動態(tài)時間積分方法。Zuijlen等[10]研究了結構動力學和流固耦合計算的顯式和隱式高階時間積分。為提高時間積分求解方法的精度,鐘萬勰[11-12]提出暫態(tài)歷程的精細計算方法,以其高精度、無條件穩(wěn)定等優(yōu)點得到了廣泛的應用。對于線性定常結構動力系統(tǒng),精細積分法可以得到逼近于精確解的數(shù)值結果。在此基礎上,王立峰等[13]研究了求解非齊次線性常微分方程的精細積分方法。張慶云等[14]研究了一種提高增維精細積分法計算精度的方法。顧元憲等[15]、張靖姝等[16]、任傳波等[17]研究了非齊次動力方程,剛柔耦合系統(tǒng)和結構動力學系統(tǒng)等問題。

時滯在各類系統(tǒng)中廣泛存在,甚至不可避免,包含時滯的系統(tǒng)可以用時滯動力學微分方程來表示。Bellen等[18]通過一種均勻修正的隱式Runge-Kutta方法計算了時滯微分方程數(shù)值解。Shakeri等[19]通過同倫攝動法求解了時滯微分方程,算例顯示該方法可以取得令人滿意的效果。包含時滯的非線性系統(tǒng)越來越多地被研究,這類問題可以利用非線性微分方程的近似解析法來計算。實際研究中,主要利用MATLAB中函數(shù)對時滯動力學微分方程進行數(shù)值求解。

現(xiàn)提出一種求解非線性微分方程的新暫態(tài)時程積分方法。在此基礎上,對包含時滯的非線性動力學微分方程進行了數(shù)值求解研究。在暫態(tài)時程積分過程中,將非線性項看作非齊次項,在瞬態(tài)區(qū)間起始時刻處進行Taylor展開并結合Romberg數(shù)值積分進行計算。Taylor展開時,將系統(tǒng)狀態(tài)方程連續(xù)引入到非線性項導數(shù)的求解過程中,來計算高階導數(shù)。對于含有時滯的非線性微分方程,將時滯項同樣看作非齊次項,通過分析時滯項暫態(tài)時程積分的上下限來確定時滯項的值。暫態(tài)時間內未知的時滯項數(shù)值利用線性插值進行計算,然后通過 Romberg 積分求得時滯積分項的數(shù)值。以強非線性振子方程為例,研究新暫態(tài)時程積分方法求解強非線性微分方程的有效性和精確度,又以包含時滯的Duffing方程和van der pol方程為例,研究新暫態(tài)時程積分方法對含有時滯的非線性微分方程的求解效果。

1 非線性動力學問題的暫態(tài)時程積分解法

對于一般的非線性問題可由以下微分方程表示：

(1)

(2)

式(2)中：H、A、B是與M、C、K有關的常系數(shù)矩陣。將式(2)右側外力項和非線性項均視為非齊次項,根據(jù)常微分方程理論,式(2)有一般解：

s)]{AF(s)+Bf[s,y(s)]}ds

(3)

取步長為Δt=tk+1-tk,對式(3)進行離散得：

s)]{AF(s)+Bf[s,y(s)]}ds=

s)]Bf[s,y(s)]ds

(4)

式(4)可寫為

(5)

式(5)中：

T=exp(HΔt)

(6)

(7)

(8)

對于式(4)右端的第1項中T可以根據(jù)文獻[12]方法精細得出,式(4)右側的第2項可以根據(jù)文獻[17]方法精細得出。

對于式(4)右側第3項通過Romberg積分方法求解得,由文獻[20]可知Romberg積分格式為

m=1,2,3…

(9)

其中：

(10)

(11)

式(9)中m的取值要視問題的求解精度而定,m的取值越大精度越高。本文中取m=2進行計算。

根據(jù)式(9)～式(11)所給格式,對式(4)中第3項進行計算,步驟如下：

(12)

tk+Δt時刻所對應的y值。

離散方程(3)的時間步長為Δt,在時間段[tk,tk+1]中,tk時刻所對應的f[tk,y(tk)]已知,因此可對非線性項f(t,y)在tk處進行泰勒展開,以求得時間段內非線性項的未知值f[tk+1/4,y(tk+1/4)]、f[tk+1/2,y(tk+1/2)]、f[tk+3/4,y(tk+3/4)]、f[tk+1，y(tk+1)]。

泰勒展開的形式如式(13)所示,展開階次為4階。

(13)

式(13)中：

(14)

將非線性系統(tǒng)狀態(tài)方程(2)代入到式(14)中,則式(14)中的各高階導數(shù)項可寫為如下形式的方程：

(15)

通過重復將系統(tǒng)狀態(tài)方程代入到式(15)中,高階導數(shù)被有效地降階。當求得Taylor展開式(13)中的系數(shù)D、E、F后,利用該式即可求得時間段[tk,tk+1]中未知時刻的非線性項值,如式(16)所示。

(16)

將式(16)所得系統(tǒng)狀態(tài)值代入Romberg積分式(12)中,則非線性項得解。利用逐步遞推的暫態(tài)時積公式(4),則可求得非線性微分式(1)時域內的數(shù)值解。

下面求解三次強非線性振子方程。三次強非線性振子方程形式為

(17)

式(17)中的參數(shù)分別取為：a=1.5,b=3。

(18)

圖1 4階Runge-Kutta法和新暫態(tài)積分法求解強非線性振子方程的結果對比

對于非線性項f=-ay13-bsiny1,按照式(13)～式(14)在t=tk處進行Taylor展開,結果如下：

(19)

將式(19)代入式(16)中,則時間區(qū)間[tk,tk+1]內未知值f[tk+1/4,y(tk+1/4)]、f[tk+1/2,y(tk+1/2)]、f[tk+3/4,y(tk+3/4)]、f[tk+1,y(tk+1)]得求,利用Romberg積分式(12),可求得非線性項的值,通過逐步遞推積分式(4),即可求得強非線性振子式(17)時域內的數(shù)值解。

利用上述方法對式(17)進行求解,仿真時長為50 s。作為對比,選擇4階Runge-Kutta法求解式(17),結果如圖1所示。由圖1可知,新暫態(tài)時程積分法可以有效求得非線性微分方程的數(shù)值解。時間步長相同時,暫態(tài)時程積分法和4階Runge-Kutte法可以取得相近的結果,如表1所示,誤差量級在10-3左右,這證明了暫態(tài)時程積分解法的有效性。

表1 新暫態(tài)時程積分法(ITTI)和Runge-Kutta法的結果和誤差

2 含時滯的非線性動力學問題的暫態(tài)時程積分解法

對于一般包含時滯的結構動力學系統(tǒng),可以用以下動力學微分方程表示為

(20)

暫態(tài)時程積分法宜于處理一階方程,因此將式(20)寫為狀態(tài)方程形式為

(21)

式(21)中：H、A、B是與M、C、K有關的常系數(shù)矩陣。將式(17)右側后3項均看作非齊次項,根據(jù)常微分理論,式(21)有一般解：

{AF(s)+Bf[s,y(s)]+Cy(s-τ)}ds

(22)

取時間步長Δt=tk+1-tk,對式(18)進行離散，即

s)]{AF(s)+Bf[s,y(s)]+

Cy(s-τ)}ds=exp(HΔt)yk+

(23)

式(23)右側前3項根據(jù)上文方法可解,時滯項采用第一部分式(9)～式(12)所給Romberg積分進行求解。

Romberg積分需要求解的方程為

(24)

進行式(12)的計算時,[tk,tk+1]時間區(qū)間內時刻

所對應的y(s-τ)值可由線性插值法計算。

(25)

則通過Romberg積分式(9)～式(12),方程中的時滯項得解,然后利用式(19)逐步遞推的積分公式即可求得含時滯非線性微分方程的數(shù)值解。

2.1 求解包含阻尼時滯擾動項的Duffing方程

包含阻尼時滯擾動項的Duffing方程可以描述為

(26)

(27)

對非線性項f=-cy13進行Taylor展開可得

(28)

通過式(9)～式(16)即可求得暫態(tài)時程積分過程中每一步內非線性項的暫態(tài)時程積分值。

式(26)～式(28)中參數(shù)值分別為：a=0.04,b=-0.2,c=0.53,A=0.4,ω=0.32,g為時滯項增益系數(shù),τ為時滯量。

當阻尼時滯擾動項參數(shù)分別為g=-1.51 N/m,τ=1.325 s時。時滯積分項中：當tk<1.325 s,y2(t-τ)|tk=0,當tk≥1.325 s時y2(t-τ)|tk=y2(tk-13 250),其中

250。

通過線性插值式(24)、式(25)和Romberg積分式(9)～式(12)即可求得暫態(tài)時程積分中每一步中時滯項的積分值。

利用逐步遞推的暫態(tài)時程積分式(23),即可解得含時滯的非線性微分方程的時域數(shù)值解。

利用MATLAB中時滯微分方程求解函數(shù)和本文方法分別對方程進行計算,可得結果如圖2所示。

時滯擾動項參數(shù)為：g=-1.51 N/m,τ=1.325 s

當阻尼時滯擾動項參數(shù)為g=-1.51 N/m,τ=0.5 s時,可得結果如圖3所示。

觀察仿真結果圖3可得,本文方法解得結果和MATLAB所解結果除開始幾秒鐘存在差異外,其余仿真結果幾乎一樣,MATLAB所繪制仿真圖像有明顯的轉折點,本文解法所繪制圖像更為連續(xù)光滑。當仿真時間為150 s,MATLAB所解數(shù)值結果共有 3 674 個,并且步長是不均勻的。本文方法時間步長為0.000 1 s,因此150 s仿真時間內共有150萬個仿真結果點。由于本文方法是一種逐步遞推的暫態(tài)時程積分方法,當仿真時長較長時的結果與MATLAB結果相同時,足以說明本文方法的正確性以及仿真時長較短時結果的精確性。因此兩種方法仿真初始時的差異可以歸結為二者的步長不同,本文方法更為精細,因此可以表現(xiàn)出MATLAB所沒有表現(xiàn)的振動特征。

2.2 求解包含阻尼時滯作用的Van Der Pol受迫振動系統(tǒng)

Van Der Pol方程揭示了由于系統(tǒng)受到非線性阻尼作用而導致系統(tǒng)出現(xiàn)孤立周期振蕩行為。包含阻尼時滯作用的Van Der Pol受迫振動系統(tǒng),可由以下方程表示：

(29)

(30)

式(29)、式(30)中的參數(shù)值分別為：a=5,b=1,c=1,A=5,ω=2.466,g=3.5 N/m,τ=1.5 s。

通過和算例2相同的方法對非線性項和時滯項進行處理,利用暫態(tài)時積分求得包含阻尼時滯項的Van Der Pol方程的時域數(shù)值解。

分別采用MATLAB中時滯微分方程求解函數(shù)和本文方法對系統(tǒng)進行仿真,可得結果如圖4所示。

時滯擾動項參數(shù)為：g=-1.51 N/m,τ=0.5 s

時滯擾動項參數(shù)為：g=-1.51 N/m,τ=0.5 s

仿真時間為600 s時,MATLAB dde23函數(shù)仿真數(shù)值結果共11 794個,步長間隔并不均勻。本文方法步長為0.000 1 s,相對于MATLAB方法計算步長更加精細。由圖4可知,在仿真開始幾秒鐘,兩種解法的結果存在小差異,原因與算例2相同,在幾秒種后兩種方法的結果幾近一致。通過兩種方法所解得的曲線,認為可取得相近的數(shù)值結果,并且本文方法更為精細,因此所得曲線更加光滑。

3 結論

(1)采用逐步遞推的暫態(tài)時程積分方法對非線性微分方程進行求解,將非線性方程中的非線性項看作非齊次項。在每一段暫態(tài)時間區(qū)間內,方程中的齊次和非齊次項分別進行暫態(tài)時程積分計算,二者疊加可作為這一段時間歷程內方程的實際狀態(tài)。暫態(tài)時間逐步向前遞推即可連續(xù)求得微分方程時域數(shù)值解。通過這種方法解得的非線性微分方程數(shù)值解和4階Runge-Kutta法可以取得相近的數(shù)值結果,二者的誤差在10-3左右。

(2) 對于包含時滯的非線性微分方程,將時滯項同樣看作是非齊次項。在暫態(tài)時間區(qū)間內,通過分析時滯項的積分上下限確定時滯項的具體值,然后結合Romberg積分方法對時滯項進行積分計算。在Romberg積分過程中,對時間段內未知的時滯項值,通過線性插值的方法進行計算。通過這種逐步遞推的暫態(tài)時程積分方法,求得的含時滯非線性微分方程的數(shù)值解和MATLAB函數(shù)的解大致相同,可以證明本文求解方法的有效性。并且本文方法相較于MATLAB解法更為精細,可控性更好,對方程中參數(shù)可以更方便地進行控制,為含時滯非線性微分方程的研究提供了新的路徑。