程志友,謝佳宏,李亞玲
(1. 重慶交通大學(xué) 航運與船舶工程學(xué)院,重慶 400074;2. 重慶交通大學(xué) 經(jīng)濟與管理學(xué)院,重慶 400074)
部分破損沉船打撈工程具有較高的風險性。由于船體破損導(dǎo)致結(jié)構(gòu)強度降低,一旦起吊時內(nèi)力過大,導(dǎo)致沉船起吊斷裂,則需要二次打撈,不僅造成經(jīng)濟損失,還可能引發(fā)海洋污染等諸多嚴重后果。因此,通過優(yōu)化起吊力參數(shù)來降低沉船起吊時的內(nèi)力,避免沉船起吊斷裂是提高沉船打撈安全性的一項重要研究內(nèi)容。
王高陽等[1]以某30萬t VLCC實船為例,選取橫艙壁移動距離和貨油質(zhì)量分布曲線控制點坐標為設(shè)計變量,以計算工況中的最大靜水彎矩值最小為目標,建立了優(yōu)化模型,并采用組合優(yōu)化算法進行優(yōu)化計算,有效降低了靜水彎矩;杜尊峰等[2]采用改進的多目標遺傳算法優(yōu)化各壓載水艙調(diào)配后的壓載水量,有效縮短了荷載傳遞時間并減小了總縱向彎矩,具有一定的工程應(yīng)用價值。
現(xiàn)有關(guān)于降低船舶內(nèi)力的研究多以正常船舶為對象,以其靜水彎矩作為優(yōu)化目標且以沉船為對象的研究相對較少。此外,正常船舶的船體結(jié)構(gòu)完整,一般僅考慮彎矩也能滿足要求。但沉船船體結(jié)構(gòu)往往存在嚴重破損,為降低沉船起吊斷裂風險,需要同時降低沉船的剪力和彎矩。
當前,NSGA-Ⅱ算法[3]是解決多目標優(yōu)化問題的經(jīng)典算法之一,它在NSGA算法的基礎(chǔ)上引入了快速非支配排序和擁擠度計算,降低了算法的時間復(fù)雜度,還引入精英策略,提高了子代種群的多樣性。
綜上,為降低沉船起吊時的內(nèi)力,從梁彎曲理論出發(fā),針對沉船起吊受力特點建立了以起吊力參數(shù)為輸入?yún)?shù)的沉船內(nèi)力數(shù)值計算模型;基于該計算模型,以起吊力參數(shù)為優(yōu)化變量,以內(nèi)力峰值為優(yōu)化目標,建立了起吊力參數(shù)優(yōu)化模型;采用改進的NSGA-Ⅱ算法進行優(yōu)化計算,為沉船打撈提供參考。
某工程案例中[4],沉船基本資料如下:沉船主尺度為92.95 m×14.8 m×7.5 m,船體水中重量約為1 607 t(文中提到的重量均指水中重量),重心為舯后7.42 m。此外,該沉船內(nèi)有兩處剩余貨物,其重量分別為143、286 t,其重心分別位于舯前22.65 m和舯后3.14 m。沉船受損情況如下:距尾垂線約25 m處有一個長度約為7 m的破洞,該處下部雙層底已經(jīng)破損,沉船完全進水。
對于資料不足的沉船,其重量分布可使用比雷斯梯形曲線近似表示[5-6],如圖1。
圖1 比雷斯梯形曲線Fig. 1 Bires trapezoidal curve
圖1中,W為重量,L為船長(取93 m),x′g為重心位置,a、b、c為比雷斯系數(shù)。該船為肥型船舶,故取b=1.174。a、c如式(1):
(1)
參考工程經(jīng)驗,使用兩對1 200 t浮筒打撈沉船,浮筒1圍繞破洞中心布置,避免破洞部分的內(nèi)力過大,浮筒2布置在舯前的適當位置。每對浮筒對稱分布在沉船兩側(cè),提供相同大小的力。貨物及浮筒在沉船上的分布示意如圖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]中僅計算了幾處關(guān)鍵彎矩且部分數(shù)據(jù)未展示。為了清楚表達優(yōu)化效果,筆者對部分缺失數(shù)據(jù)進行了合理取值。
沉船與正常船舶的受載存在明顯不同:正常船舶在重力與浮力的共同作用下保持平衡,而沉船的上升力主要由起吊力提供。因此,需要根據(jù)沉船受力特點建立沉船內(nèi)力數(shù)值計算模型。
將該沉船受載分為Fz、Qz兩類。Fz為吊纜對沉船施加的起吊力合力(圖4)。Qz為除起吊力以外其他所有載荷的等效力,包括船體以及貨物的重力、吸附力等。Fz與Qz大小相等,方向相反,g=9.8 N/kg。
常用吊纜的直徑約為0.1~0.2 m,其相對于船長較小,故可視Fz為集中力,其大小由吊纜直接提供。Qz由多個分布力組成,其分布曲線Q(x)可以按照以下方式計算[5]:
1)由沉船基本參數(shù)和比雷斯系數(shù)可以得到表示船體重量分布曲線的三等分比雷斯梯形曲線。
2)使用梯形曲線將沉船剩余貨物的重量分別等效到沉船三等分后的艏段和舯段,得到剩余貨物重量分布曲線。
3)吸附力大小可由打撈重量(船體重量+剩余貨物重量)與吸附力系數(shù)相乘近似得到。海底底質(zhì)為砂石及沙底,且沉船陷入較淺,參考工程經(jīng)驗[7]取吸附力系數(shù)為0.2。此外,假設(shè)海底土壤與船底接觸均勻,吸附力沿船長平均分布。
三者疊加即可得到Q(x),如圖5。
圖5 等效力分布曲線Q(x)Fig. 5 Equal potency distribution curve Q(x)
按一般的20站計算法[8]計算船舶內(nèi)力時,需要將起吊力等效為相鄰站上的分布載荷,而吊纜的直徑相對于站距而言較小,故將其視為集中力,直接計算沉船內(nèi)力更為合理。此外,避免對起吊力進行等效處理還有利于減少計算步驟,提高后續(xù)使用優(yōu)化算法的迭代效率。
計算船舶內(nèi)力時,可將船體視為箱體梁。根據(jù)梁彎曲理論進行計算,邊界條件均為完全自由端[6,8],即兩端的剪力和彎矩均為0。由2.2節(jié)得到Q(x)的函數(shù)表達式為:
(2)
式中:q1(x)、q2(x)、q3(x)分別為艉段、舯段、艏段上等效力分布曲線的函數(shù)表達式。
將Fz、Qz各自視為一個整體,以左端點(尾垂線所在位置)為原點可得內(nèi)力曲線方程:
(3)
(4)
式中:xi為第i個起吊力的位置參數(shù),即第i個起吊力與左端點的距離;fi為第i個起吊力的大小參數(shù),與xi一一對應(yīng),每一個起吊力的位置參數(shù)和大小參數(shù)并稱為一組起吊力參數(shù),多組起吊力參數(shù)統(tǒng)稱為起吊力參數(shù);n為起吊力的個數(shù),即有n組起吊力參數(shù);Fd(x)為0~x段上Qz的合力;Md(x)為0~x段上Qz對x處的合力矩;N(x)為沉船剪力曲線;M(x)為沉船彎矩曲線。
該計算方法的思路是先計算Fz、Qz各自引起的剪力曲線和彎矩曲線,然后疊加得到沉船的內(nèi)力曲線。該計算方法避開了將起吊力等效為分布載荷,使得計算更加合理,同時也減少了計算步驟,提高了優(yōu)化迭代的效率。
為方便編程計算,將式(3)、式(4)進行離散化處理。首先,將沉船沿船長分為m個單元,獲得m+1個截面。然后使用m+1維向量X表示x,則X的每個分量均表示沉船的一個截面,計算出X對應(yīng)的剪力向量N(X)和彎矩向量M(X),即可得到各截面上的剪力和彎矩。其中,首尾兩個端面上的內(nèi)力均為0。設(shè)Q(x)、Fd(x)、Md(x)分別由m+1維向量Q(X)、Fd(X)、Md(X)表示,可得到沉船內(nèi)力數(shù)值計算模型:
N(X)=Fd(X)-Fu(X)
(5)
M(X)=Md(X)-Mu(X)
(6)
式中:Fd(X)為Fz引起的剪力分布;Md(X)為Fz引起的彎矩分布。
若m取值越大,則單元長度|ΔX|=L/m越小,精度越高,但相應(yīng)的計算速度會降低。常用吊纜的直徑約為0.1~0.2 m,故取|ΔX| =0.000 1L(約0.01 m)較為合適。此時,式(5)、式(6)較式(3)、式(4)的相對誤差約為 ±0.05%。
要降低沉船起吊斷裂風險,需要使內(nèi)力峰值盡可能小,而內(nèi)力的正負僅表示方向,不表示大小,故以內(nèi)力峰值最小為優(yōu)化目標即可得到關(guān)于剪力和彎矩的兩個目標函數(shù):
(7)
式中:(xi,fi)為第i組起吊力參數(shù),表示第i個起吊力的位置和大小,共n組;|N(X)|為剪力絕對值向量;|M(X)|為彎矩絕對值向量。
約束條件如下:
1)平衡起吊要求:沉船起吊過程保持平衡,滿足靜平衡方程,忽略橫向載荷影響,故僅有兩個靜平衡方程。
2)破損部位的彎矩方向要求:破損部位的彎矩方向不得利于發(fā)生裂紋擴展,若破損開口向上,則破損部位的彎矩應(yīng)該小于0;反之,破損部位的彎矩應(yīng)該大于0。
3)變量范圍要求:xi理論取值范圍為0~L,而fi理論取值范圍為0~Fz。考慮到局部強度,xi不得取到破損部位。
4)其他要求:沉船起浮方式較多,不同起浮方式的起吊力參數(shù)限制也不同[9]。采用兩對浮筒起浮的限制較多(圖3),設(shè)同一浮筒上的起吊力位置相對固定,同一浮筒上的起吊力大小相等。
綜上,約束條件的數(shù)學(xué)表達式為:
(8)
式中:xg為等效力Qz的重心與沉船左端點的距離;XP∈[21,29]為破損部位對應(yīng)的截面向量;M(XP)為破損部位的彎矩向量。
內(nèi)力優(yōu)化包含剪力和彎矩兩個目標,屬于多目標優(yōu)化問題。NSGA-Ⅱ 算法作為常見的多目標優(yōu)化算法,其尋優(yōu)部分與基本遺傳算法并無差別?;具z傳算法可以不依賴問題信息進行尋優(yōu),但存在早熟問題,故需要結(jié)合內(nèi)力數(shù)值計算模型特點對其做出改進。
3.3.1 劃分種群,獨立交叉
由于起吊力參數(shù)優(yōu)化模型同時包含xi和fi兩種優(yōu)化變量,若兩者進行相互交叉,則容易相互干擾,丟失各自的優(yōu)秀基因,導(dǎo)致早熟。故當算法執(zhí)行到交叉時,將種群分割成兩個子種群,各自獨立進行交叉后,再組合成新的子代種群。
3.3.2 引入模擬退火算法思想
各組起吊力參數(shù)之間的順序關(guān)系并不唯一,故存在最終效果相同的不同解集,導(dǎo)致對種群進化過程產(chǎn)生干擾。為此,引入模擬退火算法思想[10],將變異概率p和交叉概率c視為模擬退火算法中的溫度,其隨著迭代次數(shù)增多而降低,使得算法在前期有較高的p和c,盡可能對整個空間進行搜索,提高種群的多樣性。而后期p和c逐漸降低,種群盡可能保留前面積累的優(yōu)秀個體,重點進行小范圍內(nèi)的局部搜索。設(shè)溫度下降系數(shù)為t,則改進后的算法流程如圖6。
圖6 改進的算法流程Fig. 6 Improved algorithm flow chart
通過MATLAB實現(xiàn)優(yōu)化計算,算法主要參數(shù)設(shè)置如下:種群大小為200,迭代次數(shù)為500,選擇概率為0.8,變異概率為0.5,交叉概率為0.8,溫度下降系數(shù)為0.987。計算出Pareto最優(yōu)解集對應(yīng)的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前沿拐點對應(yīng)的解作為優(yōu)化計算后的起吊力參數(shù)代表,與按照經(jīng)驗布置的起吊力參數(shù)作對比。
按照工程經(jīng)驗布置的起吊力參數(shù),記為經(jīng)驗方案;優(yōu)化計算出的起吊力參數(shù),記為優(yōu)化方案。兩種方案的起吊力參數(shù)及內(nèi)力峰值見表1。
表1 起吊參數(shù)及內(nèi)力峰值對比
由表1可知,優(yōu)化后的剪力峰值下降了約22.72%,彎矩峰值下降了約16.13%。由此可見,改進后的優(yōu)化算法對于降低沉船內(nèi)力有比較明顯的效果,其優(yōu)化計算出的解對破損沉船打撈的起吊力布置具有較大的參考價值。
從梁彎曲理論出發(fā),針對沉船起吊受力特點建立了以起吊力參數(shù)為輸入,以內(nèi)力分布為輸出的沉船內(nèi)力數(shù)值計算模型。該計算模型較20站計算法更為合理,可以作為獨立模塊計算任意起吊力參數(shù)下的沉船內(nèi)力,其計算精度可以根據(jù)需求進行調(diào)整。
基于沉船內(nèi)力數(shù)值計算模型,以起吊力參數(shù)為優(yōu)化變量,以內(nèi)力峰值為優(yōu)化目標,建立了起吊力參數(shù)優(yōu)化模型;采用改進的NSGA-Ⅱ算法進行優(yōu)化計算,計算出了降低沉船內(nèi)力的起吊力參數(shù)。將該優(yōu)化模型應(yīng)用于實際工程時,應(yīng)注意以下兩點:
1)對于涉及到艙室未完全浸水產(chǎn)生內(nèi)浮力的工況,只需要將內(nèi)浮力曲線疊加到分布曲線Q(x)中即可,后續(xù)優(yōu)化流程不變。
2)對于需要考慮其他客觀因素對起吊力參數(shù)的限制時,只需要在優(yōu)化模型的約束條件中添加相應(yīng)的約束即可,后續(xù)優(yōu)化流程不變。
此外,該優(yōu)化模型可以同時優(yōu)化剪力和彎矩,而對于其計算出的Pareto解集,也應(yīng)該綜合沉船的整體狀態(tài),針對性選擇起吊力參數(shù),盡可能降低沉船起吊斷裂的風險。