張 弢

（巴斯大學(xué) 理學(xué)院 英國 巴斯）

引 言

船舶的浮態(tài)平衡問題是靜力學(xué)中最基本的問題，也是任何浮體在設(shè)計(jì)過程中首先考慮的問題之一，中外學(xué)者多年來均對該問題進(jìn)行了深入研究。傳統(tǒng)的船舶靜力學(xué)原理基于等體積假設(shè)，即體積微幅傾斜時(shí)，浮體旋轉(zhuǎn)軸必定通過當(dāng)前水線面的漂心位置。趙曉非［1］根據(jù)靜力學(xué)原理建立平衡方程并使用牛頓法迭代求解，但每次迭代需要計(jì)算包含水線面慣性矩等參數(shù)的雅可比矩陣。桑松［2］在此基礎(chǔ)上利用最小功原理將該方法擴(kuò)展到適用于任何浮體的變軸橫傾穩(wěn)性上，并推導(dǎo)出浮式結(jié)構(gòu)空間穩(wěn)性臂計(jì)算公式。此外，孫承猛［3］使用有限差分法近似計(jì)算雅可比矩陣，避免了復(fù)雜的公式推導(dǎo)。求方程的根其實(shí)是求函數(shù)極值的特殊情況，可以轉(zhuǎn)化為數(shù)值優(yōu)化中求函數(shù)最小值的問題。因此，馬坤［4］基于非線性規(guī)劃法設(shè)定懲罰函數(shù)來求目標(biāo)函數(shù)的最小值，并將排水體積和浮心表示為3個(gè)自變量——吃水、橫傾角和縱傾角的函數(shù)，進(jìn)而避免計(jì)算水線面慣性矩等參數(shù)。之后，陸叢紅［5］使用約束優(yōu)化的遺傳算法來求解浮態(tài)方程；金寧［6］在此基礎(chǔ)上進(jìn)一步優(yōu)化遺傳算法，使目標(biāo)方程收斂更快；張維英［7］則使用差分進(jìn)化法來求解優(yōu)化問題。

其實(shí)，浮態(tài)平衡方程是二階可導(dǎo)的非線性方程組，使用非約束優(yōu)化（比如線搜索法，置信域法）就可以求解；遺傳算法、差分進(jìn)化法等卻更適合求解約束優(yōu)化的全局最小值問題，并且由于不需要計(jì)算梯度，目標(biāo)方程可以是非連續(xù)不可微的函數(shù)，因此被廣泛應(yīng)用在生物和經(jīng)濟(jì)學(xué)領(lǐng)域（若用于解決船舶浮態(tài)平衡問題，則比牛頓法和置信域法的超線性收斂速度要慢很多）。

本文前4節(jié)將建立浮態(tài)平衡的數(shù)學(xué)模型并論述兩種計(jì)算排水量和浮心的通用方法，第5和第6節(jié)介紹置信域Dogleg算法（又稱“狗腿算法”），并通過MATLAB算例驗(yàn)證該方法用于求解浮態(tài)方程的實(shí)際效果。

1 坐標(biāo)系

圖1 正浮狀態(tài)下的坐標(biāo)系

2 歐拉角和線性變換

剛體和隨剛體運(yùn)動的局部坐標(biāo)系OX′Y′Z′相對于靜止的全局坐標(biāo)系OXYZ的方位可以用3個(gè)歐拉角來表示，但同樣的方位可以用不同組合的3個(gè)歐拉角表示。因此，船舶工程中通常將繞Z′軸的轉(zhuǎn)動稱為首搖（γ），繞Y′軸的轉(zhuǎn)動稱為縱搖（β），繞X′軸的轉(zhuǎn)動稱為橫搖（α），并規(guī)定順序?yàn)槭讚u-縱搖-橫搖。在船舶靜力學(xué)中，首搖并不會改變船體的穩(wěn)性，因此歐拉角可以省略為2個(gè)。在本文中，歐拉角的正負(fù)號按照右手螺旋法則定義，即當(dāng)船首下沉?xí)r縱搖為正，右舷下沉?xí)r橫搖為正。

船體從正浮狀態(tài)（OX′Y′Z′與OXYZ重合）縱搖，可以看作是3空間中的線性變換，用矩陣P表示：

船體從正浮狀態(tài)橫搖的線性變換用矩陣R表示為：

根據(jù)定義的歐拉角順序（即先繞Y′軸縱搖再繞X′軸橫搖），合并后的線性變換表示為矩陣T：

式中：下標(biāo)G表示基底為的線性變換，下標(biāo)L表示基底為的線性變換。因Y′軸與Y軸在正浮時(shí)重合，所以［P］G=［P］L。從式（1）中可以看出先以Y′軸縱搖再以X′軸橫搖等效于先以X軸橫搖再以Y軸縱搖。由于此線性變換只發(fā)生在3空間內(nèi)，因此也可以理解為局部坐標(biāo)不變進(jìn)行基底變換，求得的矩陣相同。至此，船體傾斜后的型線坐標(biāo)（x，y，z）可由式（2）求出：

式中：（x0，y0，z0）為船體正浮的型線坐標(biāo)。

3 浮態(tài)平衡方程

船舶靜力學(xué)的模型通常只涉及到垂蕩、橫搖、縱搖3個(gè)自由度，如圖2根據(jù)牛頓定律可以表示為由3個(gè)未知量以及3個(gè)等式組成的方程組參見式（3）：

其中：

水密度ρ、船體排水量m為已知量，水線面的Z軸坐標(biāo)WLz、縱傾角β和橫傾角α是未知量，重心坐標(biāo)CG是β和α的函數(shù)。由此，式（2）便可表示為該式中：（CGx0CGy0CGz0）為船體正浮時(shí)的重心坐標(biāo)；下標(biāo)x、y和z表示全局參考系OXYZ下的坐標(biāo)。除特別說明以外，之后所有的坐標(biāo)均參考全局坐標(biāo)系OXYZ。

圖2 浮態(tài)平衡

排水體積V和浮心坐標(biāo)CB是WLz、β以及α的函數(shù)。本文所用方法不基于等體積假設(shè)，因此坐標(biāo)系原點(diǎn)無需定于水線面漂心處，而是定于如圖1和圖2所示的船尾中心基線處。這樣做的好處是：

（1）和船體型值表的參考系一致，不需進(jìn)行額外的坐標(biāo)轉(zhuǎn)換。

（2）對船長遠(yuǎn)大于船寬的普通船舶來說，縱傾對浮心變化的影響遠(yuǎn)遠(yuǎn)大于橫傾，因此將縱傾的旋轉(zhuǎn)軸設(shè)在船尾可以減小浮心的變化速度。用數(shù)值優(yōu)化方法求解時(shí)，合理的縮放目標(biāo)函數(shù)可以提高收斂的穩(wěn)定性。

（3）旋轉(zhuǎn)軸的平移并不會影響平衡后的橫傾和縱傾角度，只會影響水線面高度WLz。而此模型中，因水線面不受約束，故可在全局坐標(biāo)系Z軸上任意平移，而對最終計(jì)算結(jié)果沒有影響。

4 浮心和排水體積

船體模型通過式（2）線性變換后與水線面相交，水線面切割模型（通常由離散的多邊形組成，見下頁圖3、圖4）后得到水下的部分，其算法的原理是找到水線面和多邊形的交點(diǎn)并舍棄水線以上部分，再按順序重新連接交點(diǎn)和多邊形在水下的頂點(diǎn)。接下來，計(jì)算式（3）中的排水體積V和浮心坐標(biāo)CB主要有面元法［8］和切片法。

4.1 面元法

如下頁圖3所示，船體表面可劃分成很多平面單元，對水下單元的表面積分計(jì)算其浮力、浮力矩和體積矩等，然后用斯托克斯定理（或格林定理）將表面積分轉(zhuǎn)換為單元輪廓的封閉曲線積分，見式（4）。由于單元的輪廓為多邊形，轉(zhuǎn)換為曲線積分以后，只需對每條邊的直線方程積分即可。通過這個(gè)原理對所有水下單元的積分求和，即可得到總的排水體積和浮心。例如，當(dāng)式 （4） 中curl F·n取常數(shù)1時(shí)，表面積分為面積；當(dāng)curl F·n取z·nz時(shí)，表面積分為體積；當(dāng)curl F·n取0.5·z2·nz時(shí)，表面積分為對Z軸的體積矩。

圖3 面元模型（Heimdal小型工程船）

圖4 切片模型（Heimdal小型工程船）

式中：curl F為向量場F的旋度；n（nx，ny，nz）為平面單元的法向量；dr為曲線的切向量。

4.2 切片法

將船體切割成很多與船長方向垂直的切片（見圖4），計(jì)算每個(gè)切片水下區(qū)域的面積Aw和形心dCB，然后使用式（5）和式（6）進(jìn)行數(shù)值積分（如梯形法）求出排水體積和浮心，其中Lwl為水線長。

文獻(xiàn)［4］按照上述面元法的思路，以格林定理來求Aw和dCB，但本文給出一種更簡潔的表述方法，并可直接應(yīng)用于三維平面。在3中有任意兩個(gè)向量q1（xq1，yq1，zq1），q2（xq2，yq2，zq2），其叉積的長度為：

θ是q1逆時(shí)針旋轉(zhuǎn)到q2的角度。式（7）可以看作是q1和q2圍成的三角形q1Oq2面積的2倍。因此，任意多邊形的每個(gè)頂點(diǎn)可以表示為該平面原點(diǎn)處的向量（參見下頁圖5）并讓最后一個(gè)頂點(diǎn)和第一個(gè)頂點(diǎn)重合（q1，q2，…qn，q1），按順序?qū)⑾噜徬蛄坎娉?，則多邊形的面積為：

當(dāng) π＜θ≤ 2π 時(shí)，ε= 1 ；當(dāng) 0 ≤θ≤ π 時(shí)，ε= 2。

dCB的計(jì)算原理類似，根據(jù)相似三角形和中點(diǎn)定理推導(dǎo)出三角形的形心公式（8）。計(jì)算圖5中每個(gè)三角形qi Oqi+1的形心及面積矩，則整個(gè)多邊形的形心為可由式（9）求出。

圖5 求多邊形的面積和形心

面元法和切片法有各自的優(yōu)勢：面元法需要先將船體模型劃分成網(wǎng)格模型，但網(wǎng)格模型可以和后續(xù)的有限元結(jié)構(gòu)分析以及水動力分析共享，減少后續(xù)工作量；切片法可以直接使用型值表或型線圖的數(shù)據(jù)進(jìn)行穩(wěn)性計(jì)算，減少前期工作量。當(dāng)切片數(shù)量或網(wǎng)格數(shù)量足夠多時(shí)，兩種方法都能準(zhǔn)確地反映船體形狀，因此都是可靠穩(wěn)定的方法。

5 求解非線性方程組

式（3）可以寫成向量形式如下：

f是一個(gè)連續(xù)非線性方程組，求解式（10）相當(dāng)于求解連續(xù)函數(shù)l無約束下的極小值問題：

解決式（11）的基本思想是給出一個(gè)初始值x0，然后迭代求使：

直到發(fā)現(xiàn)極值或者找不出新的xk為止。式中：下標(biāo)k為迭代次數(shù)；pk是單位方向向量；τ是步長。

求無約束函數(shù)極值問題主要有兩類方法：線搜索法和置信域法。線搜索法（比如梯度下降法，牛頓法等）基本思想是先找到搜索方向pk滿足式（12），然后再求出合適的步長τ；而置信域法是先確定一個(gè)搜索半徑Δ，然后在Δ內(nèi)同時(shí)找到滿足式（12）的pk和τ。如果找不到，則縮小Δ直到找到滿足式（12）的pk和τ為止。

找出pk和τ的方法有多種，本文使用的置信域狗腿法具體原理［9］如下。

根據(jù)泰勒定理將目標(biāo)函數(shù)l在xk附近近似為二次函數(shù)，如式（13）所示：

當(dāng)Hk為正定矩陣時(shí)，式（14）的解如圖6所示，可分為以下3種情況。

圖6 狗腿算法求p

式中：pB為無約束條件下，式（14）的全局解。

式中：pU為無約束條件下，式（14）在-gk方向的解，

式中：τ為二次方程的解。

這樣做的好處是在搜索半徑內(nèi)找不到式（14）的全局解時(shí)，盡可能利用二次導(dǎo)數(shù)信息Hk，在P U和P B的連線上找式（14）的近似解p，參見圖6（c）。而當(dāng)搜索半徑相對較小時(shí)，式（13）可去除二次項(xiàng)，簡化為線性近似，在梯度的負(fù)方向求式（14）的近似解p，參見圖6（b）。

當(dāng)Hk為非正定矩陣時(shí)，可以將Hk對角化，并對特征值取絕對值，將Hk轉(zhuǎn)為正定矩陣。無法直接求出Hk時(shí)，可用有限差分法法求近似Hk。

搜索半徑Δ的選取通常依據(jù)目標(biāo)函數(shù)l和二次函數(shù)m的擬合比率r：

圖7 置信域狗腿算法偽代碼

6 收斂結(jié)果

為測試上述置信域狗腿法用于求解浮態(tài)方程的實(shí)際收斂效果，用MATLAB編寫船舶穩(wěn)性程序分別以牛頓法（N）和置信域狗腿法（T）計(jì)算兩艘實(shí)船在不同工況下的浮態(tài)（水線面高度WLz，橫傾角α和縱傾角β），并比較最后目標(biāo)方程的冗余值‖f‖2和迭代次數(shù)。船舶的排水量和重心坐標(biāo)為已知，同時(shí)給出迭代的初始值WLz0、α0和β0。

首先以圖3和圖4所示Heimdal小型工程船為例。該船主尺度為： 總長36 m、寬8.2 m、高3.2 m。

從下頁表1中的計(jì)算結(jié)果可以看出：工況1的兩種方法計(jì)算結(jié)果相同，船體處于接近正浮的狀態(tài)，收斂的精度和迭代次數(shù)也基本相當(dāng)；工況2的排水量增加，重心向船首右舷移動并升至2 m。若使用和工況1相同的初始值時(shí)，兩種方法的計(jì)算結(jié)果都很理想（橫傾23°，縱傾4°），但使用第3組初始值（0.3，0，0）時(shí)，牛頓法會收斂到另一個(gè)平衡位置（橫傾-88°，縱傾7°）。雖然這個(gè)解也滿足平衡方程，但與實(shí)際不符。究其原因是牛頓法并沒有限制步長，在某一步迭代中跨過了最近的解，并收斂到另一個(gè)遠(yuǎn)處的解；置信域法由于步長限制在搜索半徑內(nèi)，因此計(jì)算結(jié)果正常。使用最后一組初始值時(shí)，牛頓法在迭代至第4步時(shí)失效，原因是該步算出的雅可比矩陣不可逆，因此無法算出這一步的解；置信域法由于保證在搜索半徑內(nèi)找到解，參見式（14）-式（17），因此計(jì)算結(jié)果正常。

表1 Heimdal小型工程船浮態(tài)計(jì)算結(jié)果

下面再以某大型集裝箱船為例。該船主尺度為：總長288 m、寬32.2 m、高21.6 m。從表2所示計(jì)算結(jié)果可以看出：工況1時(shí)，初始值為（10，0，0），兩種方法計(jì)算結(jié)果相同，但置信域法的迭代次數(shù)較多。這是由于當(dāng)二次函數(shù)與目標(biāo)函數(shù)擬合情況較好時(shí)，置信域狗腿算法會逐漸增大搜索半徑（如圖6增大1倍）；牛頓法由于沒有限制步長，因此收斂更迅速，但同時(shí)也存在風(fēng)險(xiǎn)，如表1中的第3組例子。工況2時(shí)，重心向船尾右舷移動，初始值不變，但牛頓法再一次求解失敗，原因和表1中的第4組例子一樣，雅可比矩陣不可逆。

表2 某80 000 t集裝箱船浮態(tài)計(jì)算結(jié)果

需要注意的是：本文采用的牛頓法沒有限制步長是為保留其特點(diǎn)并突出和置信域法的區(qū)別，實(shí)際中也可以使用回溯法改進(jìn)牛頓法達(dá)到控制步長的目的，但仍然存在雅可比矩陣可能是奇異矩陣的問題。

7 結(jié) 論

從實(shí)例的收斂結(jié)果來看，牛頓法較為依賴初始值的選擇，因?yàn)檠趴杀染仃嚳赡墚a(chǎn)生不可逆的問題，導(dǎo)致求解中斷；而相比于牛頓法，置信域法不僅更穩(wěn)定、對初始值的依賴程度低，且能確保收斂，因此非常適合用于工業(yè)計(jì)算及穩(wěn)性軟件開發(fā)。

此外，求解浮態(tài)平衡方程的原理也可以推廣到其他穩(wěn)性相關(guān)的問題。比如計(jì)算復(fù)原力臂，只需將式（3）中的第3個(gè)方程去掉，已知橫傾角α求解船體在Z軸和X軸方向的平衡方程，然后計(jì)算浮心和重心的Y軸坐標(biāo)差即可，這種方法同時(shí)考慮了自由縱傾的影響。當(dāng)艙室破損時(shí)，可以把進(jìn)水看成附加質(zhì)量，按照第4節(jié)的方法計(jì)算艙室進(jìn)水的質(zhì)量和重心，疊加到船體本身的質(zhì)量、重心，再求解平衡方程即可。再比如計(jì)算海洋平臺的變軸橫傾穩(wěn)定力臂，可以修改第2節(jié)的線性變換矩陣，使船體沿著指定方位角橫傾即可。