劉亞沖

（中國船舶及海洋工程設(shè)計研究院 基礎(chǔ)研究部，上海200011）

0 引 言

迄今，梅爾尼科夫方法是研究確定性動力系統(tǒng)全局穩(wěn)定性的唯一解析方法。它通過構(gòu)造系統(tǒng)的梅爾尼科夫函數(shù)，求解該函數(shù)的一階簡單零點得到系統(tǒng)Poincare映射具有Smale馬蹄變換意義下混沌閾值，作為檢測混沌閾值的解析方法已被成功應(yīng)用于許多系統(tǒng)中[1-3]。一般情況下，解析求解梅爾尼科夫函數(shù)時，需要先求出系統(tǒng)的同宿或異宿軌道的參數(shù)表達式，然后將其代入到梅爾尼科夫函數(shù)中，再利用積分表或留數(shù)定理計算。但直接對軌道方程進行解析求解較為困難。針對這一問題，本文選取兩種數(shù)值方法對達芬振子系統(tǒng)求解梅爾尼科夫函數(shù)，并對兩種數(shù)值方法進行對比分析。一種是用函數(shù)逼近的方法得到軌道方程的近似表達式，而該近似表達式易于進行積分計算；另一種是通過對相函數(shù)分析，將梅爾尼科夫函數(shù)的積分轉(zhuǎn)化到易于進行的積分域上進行。

梅爾尼科夫函數(shù)有簡單零點只是系統(tǒng)出現(xiàn)混沌運動的必要條件，因此在計算出混沌閾值后還應(yīng)該進行數(shù)值驗證。龐加萊映射、安全池和李雅普諾夫指數(shù)是三種常用的混沌識別方法，本文選取李雅普諾夫指數(shù)譜對諧和激勵作用下的達芬系統(tǒng)進行數(shù)值驗證。

運動的有界性問題通常采用“安全盆”的概念來描述，它被定義為所有有界解吸引域的集合，其中安全盆侵蝕現(xiàn)象通常用來解釋動力系統(tǒng)的全局失穩(wěn)行為。這一概念最早由Thompson[4]等在研究船舶的傾覆問題時提出，并被應(yīng)用到不同的工程領(lǐng)域中。對于船舶運動系統(tǒng)而言，安全池是滿足船舶橫搖角和橫搖角速度的穩(wěn)態(tài)解落在一有界域的所有初值的集合，穩(wěn)態(tài)解一旦超出這個有界域，船舶就無法提供足夠大的復(fù)原力矩來抵抗傾覆力矩，從而發(fā)生傾覆。在本文最后，采用梅爾尼科夫方法對某型船的橫搖混沌閾值進行了計算，并觀察了安全盆從完整到侵蝕的過程。

1 非線性動力系統(tǒng)的梅爾尼科夫函數(shù)

杜芬振子作為一類非線性動力系統(tǒng)，具有廣泛的應(yīng)用背景。典型的杜芬振子的表達式為：

其中：α為阻尼系數(shù)，β和γ為恢復(fù)力系數(shù)，ζ為激振力幅值。

令 β=-1，γ=1，α=εα1，ζ=εζ1；α,ζ=o（ε）。 將上式寫成狀態(tài)向量正則方程的形式為：

當(dāng) α1=ζ1=0 時，（2）式退化到無擾系統(tǒng)（3），其 Hamilton 量和勢函數(shù)分別為（4）式和（5）式。

無擾系統(tǒng)（3）的平衡點為 s（0,0）、cl（-1,0）和cr（1,0），其中s為鞍點，cl和cr為左右中心。當(dāng)H（0,0）時，存在兩條連接鞍點的同宿軌道。過鞍點的同宿軌道滿足方程

同宿軌道有左右對稱的兩支，t∈（-∞,+∞）時，軌線從鞍點出發(fā)繞中心點順時針旋轉(zhuǎn)一周，最終再回到鞍點，同宿軌道和勢函數(shù)見圖1。由圖不難看出，x（t）為偶函數(shù)，y（t）為奇函數(shù)。

圖1 無擾系統(tǒng)的同宿軌道和勢函數(shù)圖Fig.1 Homoclinic orbit and potential function of undisturbed system

當(dāng) α≠0，ζ≠0，且足夠小時，（2）式仍為可積系統(tǒng)。按照Smale-Birkhoff同宿理論[5]，同宿軌道的穩(wěn)定流形與不穩(wěn)定流形橫截相交可能使非線性動力系統(tǒng)出現(xiàn)混沌。

設(shè)無擾系統(tǒng)同宿軌道的參數(shù)方程為xi（t）和 yi（t），受擾系統(tǒng)（2）同宿軌道的梅爾尼科夫函數(shù)可寫為：

其中：

根據(jù)梅爾尼科夫理論，若M（ t0）具有不依賴ε的簡單零點，即存在M（ t0）=0且，則對于充分小的ε，Poincare映射具有Smale馬蹄變換意義下的混沌。因此，出現(xiàn)同宿軌線橫截相交的必要條件是

2 梅爾尼科夫函數(shù)零點的數(shù)值算法

2.1 類帕德逼近方法

對于較復(fù)雜動力系統(tǒng)，理論計算同宿(或異宿)軌道的參數(shù)方程較為困難，可以采用函數(shù)逼近的思想去近似模擬。帕德逼近方法[6]是一種基于高階泰勒展開構(gòu)造設(shè)解函數(shù)的特殊方法，它以有理函數(shù)作為逼近工具，可以很好地克服泰勒展開在實際應(yīng)用中的不足之處。

由圖1可知，x（t）為偶函數(shù)且有以下初值條件：

設(shè) x （t）在不同時間t處的級數(shù)展開式為

將（11）式代入（3）式可得

將（11）式和（12）式寫為類帕德逼近多項式（QPA）的形式為

圖2 同宿軌道的類帕德逼近解（右支）Fig.2 The Quasi-Pade approximation of homoclinic(right branch)

將上式與（11）式和（12）式的泰勒級數(shù)展開進行對應(yīng)比較后可得 α0=0，α1=21.5，β1=1，因此同宿軌道參數(shù)方程的類帕德逼近函數(shù)為：

圖2是類帕德逼近解（QPA）與Runge-Kutta數(shù)值解（R-K）的對比，可以看出類帕德逼近得到的參數(shù)方程與精確解基本一致。

將（15）式代入（7）式即可求出梅爾尼科夫函數(shù)的簡單零點。

2.2 高斯—勒讓德方法

如果只需計算梅爾尼科夫函數(shù)的簡單零點，不關(guān)注軌道參數(shù)方程解的具體形式，可以繞過參數(shù)方程的求解，通過對Hamilton量做一些變換，將積分轉(zhuǎn)化到易于計算的積分域上[7]。由（4）式可得

在同宿軌道上，對上式分離變量得

其中：x0是同宿軌道與x軸的交點。

將（17）式代入A和B的表達式可得

其中：xi為鞍點。為了計算積分A和B，將積分域 [x0,xi]均勻劃分N個子區(qū)間。設(shè)其中某個子區(qū)間為[a,b]，則在[a,b]中選取三個高斯點可按下式進行：

三個高斯點對應(yīng)的權(quán)重系數(shù)分別為 w1=0.555 6，w2=0.888 9，w3=0.555 6。取 α=0.5，β=-1，γ=1，外激勵頻率變換區(qū)間為0～1.5 rad/s。在計算B時，內(nèi)外積分域應(yīng)分別劃分子區(qū)間，本文中內(nèi)部積分域劃分50個子區(qū)間，外部積分域劃分100個子區(qū)間。最終得到的計算值見圖3。采用Gauss-Legendre積分的優(yōu)點在于，無需重新構(gòu)造差值函數(shù)再進行數(shù)值積分，直接利用Gauss-Legendre公式的離散形式即可。此外，高斯積分點本身就是不等分積分點，因此具有較高的計算效率和精度。

圖3 混沌閾值曲線（類Pade逼近解與Gauss-Legendre解）Fig.3 Chaos threshold curve(QPA versus Gauss-Legendre)

由于類Pade逼近得到的軌道參數(shù)方程與精確解幾乎無差異，圖3中類Pade逼近解可作為精確解。Gauss-Legendre解比類Pade逼近解稍大。由于Gauss-Legendre方法在計算梅爾尼科夫函數(shù)零點時易于實施，在不需要軌道方程的參數(shù)解時可以采用Gauss-Legendre方法。

2.3 數(shù)值驗證

由梅爾尼科夫函數(shù)計算出一階零點只能說明系統(tǒng)存在不變集Λ，這只是系統(tǒng)出現(xiàn)混沌的必要條件[8]。下面將通過數(shù)值計算方法即Lyapunov指數(shù)譜對結(jié)果進行進一步分析。

Lyapunov指數(shù)表示系統(tǒng)在經(jīng)過充分長時間的演變后，單位時間內(nèi)所引起的指數(shù)分離[5]。對于本文中連續(xù)的微分動力系統(tǒng)，選取RHR算法[9]，并對ω=1 rad/s這一激勵頻率進行驗證。由上文計算可得，ζ＞0.811 1α即ζ＞0.406時系統(tǒng)可能出現(xiàn)混沌。分別計算ζ1=0.35和ζ2=0.45時系統(tǒng)的Lyapunov指數(shù)，數(shù)值計算結(jié)果見圖4。

圖4 不同激勵幅值下Lyapunov指數(shù)譜Fig.4 Lyapunov exponents of different excitation range

從圖4中可以看出，當(dāng)激勵幅值小于混沌閾值時，系統(tǒng)的最大Lyapunov指數(shù)小于零，軌道會收縮并最終趨于穩(wěn)定值。當(dāng)激勵幅值大于混沌閾值時，系統(tǒng)的最大Lyapunov指數(shù)大于零，表示軌道分離，運動已不再穩(wěn)定。

3 船舶非線性橫搖動力系統(tǒng)

作為一單自由度的動力系統(tǒng)，船舶橫搖的非線性體現(xiàn)在阻尼力矩的非線性和恢復(fù)力矩的非線性。該橫搖方程的表達式見下式：

其中：Jφφ和ΔJφφ為轉(zhuǎn)動慣量和附加轉(zhuǎn)動慣量，為阻尼項，

因船舶橫搖運動的左右對稱性，非線性恢復(fù)力矩中只有奇次方項?？赏ㄟ^對消滅曲線擬合得到。（21）式兩邊同除以 （Jφφ+ΔJφφ）可得

在（21）式引入小參數(shù) 0＜ε＜＜1，則有

表1 某型船主要參數(shù)列表Tab.1 List of parameters of ship

選取某型船[10]，船型參數(shù)見表1。

該橫搖系統(tǒng)的正則方程形式為：

系統(tǒng)的相圖和勢函數(shù)見圖5。圖中sl是左鞍點，sr是右鞍點，c是中心點。sl和sr即對應(yīng)橫搖的穩(wěn)性消失角，其中sl=-sr=-1.543 rad。

圖5 船舶橫搖系統(tǒng)的相圖和勢函數(shù)Fig.5 Phase plan and potential function of ship-roll system

圖6 船舶橫搖系統(tǒng)的混沌閾值曲線Fig.6 Chaos threshold curve of ship roll system

該橫搖系統(tǒng)的梅爾尼科夫函數(shù)為：

將表1中參數(shù)代入上式，波浪頻率范圍取Ω∈ [0～2.0 rad/s]，將積分域從時間域轉(zhuǎn)到空間域并設(shè)置適當(dāng)?shù)淖訁^(qū)間和高斯積分點可得到該系統(tǒng)的混沌閾值曲線，見圖6。

選取波浪激勵頻率Ω=1.0 rad/s下的混沌閾值進行數(shù)值驗證，此時激勵幅值的閾值為0.031。選取不同的激勵幅值計算得到船舶橫搖的安全池見圖7-9。

從圖中可以看到，當(dāng)激勵幅值小于閾值時，橫搖系統(tǒng)的安全池完整，隨著激勵幅值的增大，安全池逐步發(fā)生破損。在破損域中選取a、b兩點作為初值計算得到的相軌跡見圖10。從a、b兩點出發(fā)的相軌跡繞中心點逐漸振蕩發(fā)散，當(dāng)橫搖角達到左鞍點即穩(wěn)性消失角時，橫搖角和橫搖角速度將持續(xù)增大，此時傾覆發(fā)生。

圖7 橫搖安全池（f=0.025）Fig.7 Safe basin(f=0.025)

圖8 橫搖安全池（f=0.04）Fig.8 Safe basin(f=0.04)

圖9 橫搖安全池（f=0.08）Fig.9 Safe basin(f=0.08)

圖10 破損點的相軌跡Fig.10 The phase trajectory of breakage points

4 結(jié) 語

梅爾尼科夫方法是預(yù)測非線性動力系統(tǒng)出現(xiàn)混沌閾值強有力的解析方法。針對在求解梅爾尼科夫函數(shù)積分時遇到的困難，本文采用了類Pade逼近和Gauss-Legendre積分兩種處理方法，并對比分析了兩種方法的計算結(jié)果及各自的優(yōu)缺點。在不需要軌道參數(shù)方程的具體形式時，可以采用易于工程計算的Gauss-Legendre積分。作為驗證，采用時域仿真計算了系統(tǒng)的Lyapunov指數(shù)譜。最后采用梅爾尼科夫方法對某船舶的橫搖動力系統(tǒng)進行分析，得到了橫搖混沌閾值曲線，觀察了安全池隨激勵幅值增大而逐漸破損的過程，并選取破損域中的點計算了相軌跡。本文工作對于研究船舶傾覆機理，進行穩(wěn)性評估有一定意義。