船舶橫搖混沌閾值的數值方法研究

劉亞沖

（中國船舶及海洋工程設計研究院 基礎研究部，上海200011）

0 引 言

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

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

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

1 非線性動力系統的梅爾尼科夫函數

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

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

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

當 α1=ζ1=0 時，（2）式退化到無擾系統（3），其 Hamilton 量和勢函數分別為（4）式和（5）式。

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

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

圖1 無擾系統的同宿軌道和勢函數圖Fig.1 Homoclinic orbit and potential function of undisturbed system

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

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

其中：

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

2 梅爾尼科夫函數零點的數值算法

2.1 類帕德逼近方法

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

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

設 x （t）在不同時間t處的級數展開式為

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

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

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

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

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

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

2.2 高斯—勒讓德方法

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

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

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

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

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

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

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

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

2.3 數值驗證

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

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

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

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

3 船舶非線性橫搖動力系統

作為一單自由度的動力系統，船舶橫搖的非線性體現在阻尼力矩的非線性和恢復力矩的非線性。該橫搖方程的表達式見下式：

其中：Jφφ和ΔJφφ為轉動慣量和附加轉動慣量，為阻尼項，

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

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

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

選取某型船[10]，船型參數見表1。

該橫搖系統的正則方程形式為：

系統的相圖和勢函數見圖5。圖中sl是左鞍點，sr是右鞍點，c是中心點。sl和sr即對應橫搖的穩性消失角，其中sl=-sr=-1.543 rad。

圖5 船舶橫搖系統的相圖和勢函數Fig.5 Phase plan and potential function of ship-roll system

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

該橫搖系統的梅爾尼科夫函數為：

將表1中參數代入上式，波浪頻率范圍取Ω∈ [0～2.0 rad/s]，將積分域從時間域轉到空間域并設置適當的子區間和高斯積分點可得到該系統的混沌閾值曲線，見圖6。

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

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

圖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 結 語

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