基于改進NSGA-Ⅱ的破損沉船起吊力優化

程志友,謝佳宏,李亞玲

(1. 重慶交通大學 航運與船舶工程學院,重慶 400074;2. 重慶交通大學 經濟與管理學院,重慶 400074)

0 引 言

部分破損沉船打撈工程具有較高的風險性。由于船體破損導致結構強度降低,一旦起吊時內力過大,導致沉船起吊斷裂,則需要二次打撈,不僅造成經濟損失,還可能引發海洋污染等諸多嚴重后果。因此,通過優化起吊力參數來降低沉船起吊時的內力,避免沉船起吊斷裂是提高沉船打撈安全性的一項重要研究內容。

王高陽等[1]以某30萬t VLCC實船為例,選取橫艙壁移動距離和貨油質量分布曲線控制點坐標為設計變量,以計算工況中的最大靜水彎矩值最小為目標,建立了優化模型,并采用組合優化算法進行優化計算,有效降低了靜水彎矩;杜尊峰等[2]采用改進的多目標遺傳算法優化各壓載水艙調配后的壓載水量,有效縮短了荷載傳遞時間并減小了總縱向彎矩,具有一定的工程應用價值。

現有關于降低船舶內力的研究多以正常船舶為對象,以其靜水彎矩作為優化目標且以沉船為對象的研究相對較少。此外,正常船舶的船體結構完整,一般僅考慮彎矩也能滿足要求。但沉船船體結構往往存在嚴重破損,為降低沉船起吊斷裂風險,需要同時降低沉船的剪力和彎矩。

當前,NSGA-Ⅱ算法[3]是解決多目標優化問題的經典算法之一,它在NSGA算法的基礎上引入了快速非支配排序和擁擠度計算,降低了算法的時間復雜度,還引入精英策略,提高了子代種群的多樣性。

綜上,為降低沉船起吊時的內力,從梁彎曲理論出發,針對沉船起吊受力特點建立了以起吊力參數為輸入參數的沉船內力數值計算模型;基于該計算模型,以起吊力參數為優化變量,以內力峰值為優化目標,建立了起吊力參數優化模型;采用改進的NSGA-Ⅱ算法進行優化計算,為沉船打撈提供參考。

1 工程概況

某工程案例中[4],沉船基本資料如下:沉船主尺度為92.95 m×14.8 m×7.5 m,船體水中重量約為1 607 t(文中提到的重量均指水中重量),重心為舯后7.42 m。此外,該沉船內有兩處剩余貨物,其重量分別為143、286 t,其重心分別位于舯前22.65 m和舯后3.14 m。沉船受損情況如下:距尾垂線約25 m處有一個長度約為7 m的破洞,該處下部雙層底已經破損,沉船完全進水。

對于資料不足的沉船,其重量分布可使用比雷斯梯形曲線近似表示[5-6],如圖1。

圖1 比雷斯梯形曲線Fig. 1 Bires trapezoidal curve

圖1中,W為重量,L為船長(取93 m),x′g為重心位置,a、b、c為比雷斯系數。該船為肥型船舶,故取b=1.174。a、c如式(1):

(1)

參考工程經驗,使用兩對1 200 t浮筒打撈沉船,浮筒1圍繞破洞中心布置,避免破洞部分的內力過大,浮筒2布置在舯前的適當位置。每對浮筒對稱分布在沉船兩側,提供相同大小的力。貨物及浮筒在沉船上的分布示意如圖2。

圖2 貨物及浮筒的分布示意(單位:m)Fig. 2 Distribution diagram of cargo and buoy

每對浮筒通過對稱分布的4根過船底吊纜對沉船施加吊力,如圖3。

圖3 浮筒受力示意(單位:m)Fig. 3 Diagram of the force on the buoy

在水下起吊階段,沉船離底最后時刻的受力達到最大,也最危險,故取該時刻的沉船為研究對象。實際工程中會控制傾角,較小傾角可以忽略,其受力如圖4。

圖 4 沉船受力示意Fig. 4 Diagram of force on shipwreck

需要指出的是,文獻[4]中僅計算了幾處關鍵彎矩且部分數據未展示。為了清楚表達優化效果,筆者對部分缺失數據進行了合理取值。

2 沉船內力數值計算模型

沉船與正常船舶的受載存在明顯不同:正常船舶在重力與浮力的共同作用下保持平衡,而沉船的上升力主要由起吊力提供。因此,需要根據沉船受力特點建立沉船內力數值計算模型。

2.1 受載分類

將該沉船受載分為Fz、Qz兩類。Fz為吊纜對沉船施加的起吊力合力(圖4)。Qz為除起吊力以外其他所有載荷的等效力,包括船體以及貨物的重力、吸附力等。Fz與Qz大小相等,方向相反,g=9.8 N/kg。

2.2 受載計算

常用吊纜的直徑約為0.1～0.2 m,其相對于船長較小,故可視Fz為集中力,其大小由吊纜直接提供。Qz由多個分布力組成,其分布曲線Q(x)可以按照以下方式計算[5]:

1)由沉船基本參數和比雷斯系數可以得到表示船體重量分布曲線的三等分比雷斯梯形曲線。

2)使用梯形曲線將沉船剩余貨物的重量分別等效到沉船三等分后的艏段和舯段,得到剩余貨物重量分布曲線。

3)吸附力大小可由打撈重量(船體重量+剩余貨物重量)與吸附力系數相乘近似得到。海底底質為砂石及沙底,且沉船陷入較淺,參考工程經驗[7]取吸附力系數為0.2。此外,假設海底土壤與船底接觸均勻,吸附力沿船長平均分布。

三者疊加即可得到Q(x),如圖5。

圖5 等效力分布曲線Q(x)Fig. 5 Equal potency distribution curve Q(x)

2.3 內力計算方法

按一般的20站計算法[8]計算船舶內力時,需要將起吊力等效為相鄰站上的分布載荷,而吊纜的直徑相對于站距而言較小,故將其視為集中力,直接計算沉船內力更為合理。此外,避免對起吊力進行等效處理還有利于減少計算步驟,提高后續使用優化算法的迭代效率。

計算船舶內力時,可將船體視為箱體梁。根據梁彎曲理論進行計算,邊界條件均為完全自由端[6,8],即兩端的剪力和彎矩均為0。由2.2節得到Q(x)的函數表達式為:

(2)

式中:q1(x)、q2(x)、q3(x)分別為艉段、舯段、艏段上等效力分布曲線的函數表達式。

將Fz、Qz各自視為一個整體,以左端點(尾垂線所在位置)為原點可得內力曲線方程:

(3)

(4)

式中:xi為第i個起吊力的位置參數,即第i個起吊力與左端點的距離;fi為第i個起吊力的大小參數,與xi一一對應,每一個起吊力的位置參數和大小參數并稱為一組起吊力參數,多組起吊力參數統稱為起吊力參數;n為起吊力的個數,即有n組起吊力參數;Fd(x)為0～x段上Qz的合力;Md(x)為0～x段上Qz對x處的合力矩;N(x)為沉船剪力曲線;M(x)為沉船彎矩曲線。

該計算方法的思路是先計算Fz、Qz各自引起的剪力曲線和彎矩曲線,然后疊加得到沉船的內力曲線。該計算方法避開了將起吊力等效為分布載荷,使得計算更加合理,同時也減少了計算步驟,提高了優化迭代的效率。

2.4 任意站面的內力計算

為方便編程計算,將式(3)、式(4)進行離散化處理。首先,將沉船沿船長分為m個單元,獲得m+1個截面。然后使用m+1維向量X表示x,則X的每個分量均表示沉船的一個截面,計算出X對應的剪力向量N(X)和彎矩向量M(X),即可得到各截面上的剪力和彎矩。其中,首尾兩個端面上的內力均為0。設Q(x)、Fd(x)、Md(x)分別由m+1維向量Q(X)、Fd(X)、Md(X)表示,可得到沉船內力數值計算模型:

N(X)=Fd(X)-Fu(X)

(5)

M(X)=Md(X)-Mu(X)

(6)

式中:Fd(X)為Fz引起的剪力分布;Md(X)為Fz引起的彎矩分布。

若m取值越大,則單元長度|ΔX|=L/m越小,精度越高,但相應的計算速度會降低。常用吊纜的直徑約為0.1～0.2 m,故取|ΔX| =0.000 1L(約0.01 m)較為合適。此時,式(5)、式(6)較式(3)、式(4)的相對誤差約為 ±0.05%。

3 起吊力參數優化模型

3.1 目標函數

要降低沉船起吊斷裂風險,需要使內力峰值盡可能小,而內力的正負僅表示方向,不表示大小,故以內力峰值最小為優化目標即可得到關于剪力和彎矩的兩個目標函數:

(7)

式中:(xi,fi)為第i組起吊力參數,表示第i個起吊力的位置和大小,共n組;|N(X)|為剪力絕對值向量;|M(X)|為彎矩絕對值向量。

3.2 約束條件

約束條件如下:

1)平衡起吊要求:沉船起吊過程保持平衡,滿足靜平衡方程,忽略橫向載荷影響,故僅有兩個靜平衡方程。

2)破損部位的彎矩方向要求:破損部位的彎矩方向不得利于發生裂紋擴展,若破損開口向上,則破損部位的彎矩應該小于0;反之,破損部位的彎矩應該大于0。

3)變量范圍要求:xi理論取值范圍為0～L,而fi理論取值范圍為0～Fz。考慮到局部強度,xi不得取到破損部位。

4)其他要求:沉船起浮方式較多,不同起浮方式的起吊力參數限制也不同[9]。采用兩對浮筒起浮的限制較多(圖3),設同一浮筒上的起吊力位置相對固定,同一浮筒上的起吊力大小相等。

綜上,約束條件的數學表達式為:

(8)

式中:xg為等效力Qz的重心與沉船左端點的距離;XP∈[21,29]為破損部位對應的截面向量;M(XP)為破損部位的彎矩向量。

3.3 NSGA-Ⅱ算法改進

內力優化包含剪力和彎矩兩個目標,屬于多目標優化問題。NSGA-Ⅱ 算法作為常見的多目標優化算法,其尋優部分與基本遺傳算法并無差別。基本遺傳算法可以不依賴問題信息進行尋優,但存在早熟問題,故需要結合內力數值計算模型特點對其做出改進。

3.3.1 劃分種群,獨立交叉

由于起吊力參數優化模型同時包含xi和fi兩種優化變量,若兩者進行相互交叉,則容易相互干擾,丟失各自的優秀基因,導致早熟。故當算法執行到交叉時,將種群分割成兩個子種群,各自獨立進行交叉后,再組合成新的子代種群。

3.3.2 引入模擬退火算法思想

各組起吊力參數之間的順序關系并不唯一,故存在最終效果相同的不同解集,導致對種群進化過程產生干擾。為此,引入模擬退火算法思想[10],將變異概率p和交叉概率c視為模擬退火算法中的溫度,其隨著迭代次數增多而降低,使得算法在前期有較高的p和c,盡可能對整個空間進行搜索,提高種群的多樣性。而后期p和c逐漸降低,種群盡可能保留前面積累的優秀個體,重點進行小范圍內的局部搜索。設溫度下降系數為t,則改進后的算法流程如圖6。

圖6 改進的算法流程Fig. 6 Improved algorithm flow chart

4 計算結果及對比分析

4.1 計算結果

通過MATLAB實現優化計算,算法主要參數設置如下:種群大小為200,迭代次數為500,選擇概率為0.8,變異概率為0.5,交叉概率為0.8,溫度下降系數為0.987。計算出Pareto最優解集對應的Pareto前沿,如圖7。

圖7 Pareto前沿Fig. 7 Pareto front

由圖7可知:剪力峰值最大為6 050.97 kN,此時的彎矩峰值最小,其值為33 425.82 kN·m;剪力峰值最小為4 003.00 kN,此時的彎矩峰值最大,其值為34 990.82 kN·m。綜合考慮剪力和彎矩,取Pareto前沿拐點對應的解作為優化計算后的起吊力參數代表,與按照經驗布置的起吊力參數作對比。

4.2 結果對比分析

按照工程經驗布置的起吊力參數,記為經驗方案;優化計算出的起吊力參數,記為優化方案。兩種方案的起吊力參數及內力峰值見表1。

表1 起吊參數及內力峰值對比

由表1可知,優化后的剪力峰值下降了約22.72%,彎矩峰值下降了約16.13%。由此可見,改進后的優化算法對于降低沉船內力有比較明顯的效果,其優化計算出的解對破損沉船打撈的起吊力布置具有較大的參考價值。

5 結 語

從梁彎曲理論出發,針對沉船起吊受力特點建立了以起吊力參數為輸入,以內力分布為輸出的沉船內力數值計算模型。該計算模型較20站計算法更為合理,可以作為獨立模塊計算任意起吊力參數下的沉船內力,其計算精度可以根據需求進行調整。

基于沉船內力數值計算模型,以起吊力參數為優化變量,以內力峰值為優化目標,建立了起吊力參數優化模型;采用改進的NSGA-Ⅱ算法進行優化計算,計算出了降低沉船內力的起吊力參數。將該優化模型應用于實際工程時,應注意以下兩點:

1)對于涉及到艙室未完全浸水產生內浮力的工況,只需要將內浮力曲線疊加到分布曲線Q(x)中即可,后續優化流程不變。

2)對于需要考慮其他客觀因素對起吊力參數的限制時,只需要在優化模型的約束條件中添加相應的約束即可,后續優化流程不變。

此外,該優化模型可以同時優化剪力和彎矩,而對于其計算出的Pareto解集,也應該綜合沉船的整體狀態,針對性選擇起吊力參數,盡可能降低沉船起吊斷裂的風險。