船舶橫搖混沌閾值的數(shù)值方法研究

劉亞沖

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

0 引 言

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

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

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

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

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

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

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

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

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

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

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

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

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

其中：

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

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

2.1 類帕德逼近方法

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

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

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

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

將（11）式和（12）式寫為類帕德逼近多項(xiàng)式（QPA）的形式為

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

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

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

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

2.2 高斯—勒讓德方法

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

在同宿軌道上，對(duì)上式分離變量得

其中：x0是同宿軌道與x軸的交點(diǎn)。

將（17）式代入A和B的表達(dá)式可得

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

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

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

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

2.3 數(shù)值驗(yàn)證

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

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

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

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

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

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

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

因船舶橫搖運(yùn)動(dòng)的左右對(duì)稱性，非線性恢復(fù)力矩中只有奇次方項(xiàng)。可通過對(duì)消滅曲線擬合得到。（21）式兩邊同除以 （Jφφ+ΔJφφ）可得

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

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

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

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

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

圖5 船舶橫搖系統(tǒng)的相圖和勢(shì)函數(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]，將積分域從時(shí)間域轉(zhuǎn)到空間域并設(shè)置適當(dāng)?shù)淖訁^(qū)間和高斯積分點(diǎn)可得到該系統(tǒng)的混沌閾值曲線，見圖6。

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

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

4 結(jié) 語

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