劉虓 徐磊 樊天慧 陳超核

(華南理工大學(xué) 土木與交通學(xué)院，廣東 廣州 510640)

隨著世界各國尤其是涉海國家正在積極實(shí)施海洋戰(zhàn)略，以浮式海洋平臺為代表的各種超常規(guī)外形的水面浮體不斷涌現(xiàn)。這些浮體的結(jié)構(gòu)形式與常見的水面浮體——單體船的最大區(qū)別在于：前者由多個(gè)凸多面體組成，而后者只有一個(gè)凸多面體(忽略舵、槳及其它細(xì)小附體結(jié)構(gòu))。浮態(tài)是所有浮體性能的基礎(chǔ)性指標(biāo)。單體船[1- 2]的浮態(tài)研究已經(jīng)成熟，其計(jì)算方法深刻影響了以多體船[3]和浮式海洋平臺[4- 5]為代表的凸多面體組合結(jié)構(gòu)的浮態(tài)計(jì)算。近年來隨著計(jì)算技術(shù)的進(jìn)步，浮態(tài)計(jì)算方法尤其是凸多面體組合結(jié)構(gòu)的浮態(tài)計(jì)算出現(xiàn)了新的趨勢和演變。

首先是計(jì)算模型的新趨勢。浮體的研究對象從常規(guī)船舶過渡到多體船、海洋平臺乃至超大型浮式平臺等，浮體外形的描述方法由基于二維切面的常規(guī)船舶型值表，發(fā)展為三維“面元”[5- 6]模型乃至三維“體素”模型[7- 8]。建立三維面元模型的圖形軟件/格式/庫有：CATIA軟件[9]、STL圖形格式[10]、CGAL幾何算法庫[11]、有限元軟件PATRAN[12]等。三維“體素”模型[8]的提出基于“離散組合”的思想：將復(fù)雜形狀的浮體分解為若干簡單的“體素”(圓柱體、長方體、棱錐等)，首先求解各體素的排水體積、浮心坐標(biāo)及水線面要素，然后進(jìn)行組合得到整個(gè)浮體的浮態(tài)要素。四面體是最基本的體素單元，其它所有的體素均可分解為四面體。劉虓[7]據(jù)此提出了基于四面體的計(jì)算模型，陳彥璋[5]也利用四面體計(jì)算海洋平臺的體積要素。但是他們并沒有將基于四面體模型的方法與傳統(tǒng)的基于二維切面模型的算法進(jìn)行比較研究。

其次是計(jì)算方法的新趨勢。浮態(tài)計(jì)算是一個(gè)反復(fù)迭代的過程，傳統(tǒng)方法分為兩類：矩陣迭代法和非線性優(yōu)化算法。趙曉非[13]提出船舶浮態(tài)計(jì)算的矩陣迭代方法。該方法每次迭代需計(jì)算船舶水面要素(投影面積、漂心坐標(biāo)、慣性矩)和水下體積要素(排水體積和浮心坐標(biāo))，比較繁瑣。近期，學(xué)者們又傾向于將浮態(tài)問題歸結(jié)于非線性優(yōu)化問題，避免了水線面要素計(jì)算[14]。但是這些方法的收斂性存在問題：以Powell法和0.618法為代表的直接搜索法容易陷入局部最優(yōu)解；以遺傳算法為代表的現(xiàn)代啟發(fā)式算法又不能保證每次計(jì)算都能收斂至最優(yōu)解。

文中基于四面體網(wǎng)格模型提出了一種新型的雙重迭代法，通過內(nèi)、外雙層迭代，分別模擬浮體的升沉和旋轉(zhuǎn)運(yùn)動(dòng)，具有明確的物理意義。文中還通過實(shí)驗(yàn)證明了算法的可靠性、精度和良好的收斂性，并進(jìn)一步通過數(shù)值實(shí)驗(yàn)比較了傳統(tǒng)算法和文中算法在精度計(jì)算上的差異和產(chǎn)生原因。實(shí)驗(yàn)對照分析表明：文中方法比傳統(tǒng)方法更適合分析“凸多面體”[13]組合結(jié)構(gòu)的浮態(tài)問題。

1 預(yù)備知識

1.1 凸多面體、坐標(biāo)系和包圍盒

如果多面體在它們每一面所決定的平面的同一側(cè)，則稱此多面體為凸多面體[15]。圖1所示的半潛式海洋平臺即是由立柱、橫箱、下浮體等多個(gè)“凸多面體”構(gòu)成的組合結(jié)構(gòu)。

圖1 半潛式海洋平臺

文中采用浮體坐標(biāo)系xyz和地球/全局坐標(biāo)系xGyGzG這兩個(gè)右手坐標(biāo)系來研究浮體平衡問題。包含浮體的最小長方體稱為“包圍盒”(如圖2所示)。

圖2 浮體坐標(biāo)系xyz、全局坐標(biāo)系xGyGzG和包圍盒

1.2 浮態(tài)參數(shù)/水平面方程

文中假定浮體不動(dòng)，其浮態(tài)體現(xiàn)為浮體坐標(biāo)系下的水平面方程：

nx·(x-xp)+ny·(y-yp)+nz·(z-zp)=0

(1)

式中：水平面的法向量N為〈nx，ny，nz〉，面內(nèi)某點(diǎn)p的坐標(biāo)為(xp，yp，zp)。

平面方程一般表達(dá)式為：

Ax+By+Cz+D=0

(2)

對比式(1)、(2)，得：

(3)

1.3 浮體平衡條件

浮體平衡條件如下：

(Ⅰ)V=P/γ；

(Ⅱ)l(B，G，N)= 0。

式中，V為浮體排水體積，P為浮體重量，γ為水的重度，B為浮體浮心坐標(biāo)，G為浮體重心坐標(biāo)，N為水平面法向量(假設(shè)浮體固定，水面法向量可變)，l為浮力作用線與重力作用線的距離。

(4)

1.4 基于四面體網(wǎng)格的浮態(tài)計(jì)算

四面體(三棱錐)是最基本的三維體單元(如圖3所示)。任何三維浮體均可分解為四面體單元(如圖4所示)。文獻(xiàn)[7]中提供了長方體、圓柱體、三棱柱分解為四面體的方法。

四面體體積為

圖3 四面體單元

圖4 利用FEM軟件將船體分解為多個(gè)四面體單元示意圖

(5)

四面體形心坐標(biāo)為

(6)

式中，四面體(如圖3所示)的頂點(diǎn)i坐標(biāo)為(xi，yi，zi)，其中：i=A，B，C，D。

將三維浮體分解為多個(gè)四面體后，只要逐個(gè)計(jì)算每個(gè)四面體的排水量和浮心坐標(biāo)，就可匯總得到整個(gè)浮體的排水量和浮心坐標(biāo)。文獻(xiàn)[7]中列舉了四面體與水平面相對位置的5種可能情況(如圖5所示)，并逐一進(jìn)行了分析。

圖5 四面體與水面的切割關(guān)系(共5種，此處僅列舉2種)

2 雙重迭代法求浮態(tài)

文中提出的浮態(tài)求解方法是一種迭代法，包含內(nèi)、外兩層迭代。

2.1 內(nèi)層迭代(模擬升沉)

假設(shè)浮體只進(jìn)行升沉運(yùn)動(dòng)(暫時(shí)約束其旋轉(zhuǎn)運(yùn)動(dòng))，可用內(nèi)層迭代法來確定水平面方程。展開具體論述前，先觀察圖6所示的浮體包圍盒，并令8個(gè)頂點(diǎn)距池底的距離為

hi=lxi+myi+nzi，i=1，2，…，8

(7)

式中：l、m、n為水平面法向N(圖6中全局坐標(biāo)系的zG軸指向)與浮體坐標(biāo)系x、y、z軸夾角的余弦。

圖6 浮體包圍盒示意圖

找出最小值h0及其對應(yīng)頂點(diǎn)v0；最大值h1及其對應(yīng)頂點(diǎn)v1。如圖6所示，2號點(diǎn)為v0，距離池底的距離22’即為最小值h0，7號點(diǎn)為v1，距離池底的距離77’即為最大值h1。當(dāng)水平面通過v0(2號點(diǎn))時(shí)，浮體完全露出水面，排水量為0；當(dāng)水平面通過v1(7號點(diǎn))時(shí)，浮體完全淹沒，排水量達(dá)到最大值。真實(shí)的水平面必然通過v0v1的連線(2號點(diǎn)和7號點(diǎn)的連線)。在v0和v1的連線上一定存在某點(diǎn)vt(xt，yt，zt)：當(dāng)水平面通過vt，浮體排水體積V滿足浮體平衡條件(Ⅰ)。

不妨令vt點(diǎn)的坐標(biāo)為

(8)

式中：0≤t≤1。

可使用二分法求t，然后代入式(8)和(1)即可求出水平面方程。下面給出內(nèi)層迭代法求水平面方程的流程圖，如圖7所示。

將t0代入式(8)和(3)計(jì)算水平面方程的4個(gè)參數(shù)，然后分析每個(gè)四面體與水平面的相交關(guān)系，計(jì)算每個(gè)四面體的排水量Vi(0)，最后匯總得到t0對應(yīng)的浮體的排水量即得圖7中V0：

(9)

同理，得到t1和tm對應(yīng)的浮體排水量：

(10)

(11)

圖7 內(nèi)層迭代法流程圖(模擬浮體升沉運(yùn)動(dòng))

2.2 外層迭代(模擬旋轉(zhuǎn))

根據(jù)內(nèi)層迭代法固然能夠令浮體浮力=重量(平衡條件I)，但重力作用線與浮力作用線不一定重合(平衡條件II)。如圖8所示，浮力作用線通過浮心B，方向與水平面法向N相同；重力作用線通過重心G，方向與N相反。令兩條作用線之間的距離為l，如果l≠0，則水平面法向N將會繞某個(gè)軸Axis(假設(shè)浮體不動(dòng))旋轉(zhuǎn)?？筛鶕?jù)右手螺旋法則確定Axis：

(12)

圖8 浮體旋轉(zhuǎn)軸

不妨假設(shè)初始狀態(tài)下，浮體在水面處于正浮狀態(tài)。浮體通過升沉運(yùn)動(dòng)(內(nèi)層迭代法模擬)使得浮力=重力，但如果重力和浮力作用線不重合，那么浮力和重力形成的力偶將驅(qū)動(dòng)浮體旋轉(zhuǎn)(外層迭代法模擬)。假定浮體不動(dòng)，水面旋轉(zhuǎn)，則N繞Axis旋轉(zhuǎn)，產(chǎn)生新的法線方向N′。通過內(nèi)層迭代法再次模擬浮體沿N′的升沉運(yùn)動(dòng)，浮心B將移至新的位置B′，浮力和重力作用線之間的距離l將隨之減小。如果l≠0，N′將繼續(xù)旋轉(zhuǎn)，l進(jìn)一步減小……如此循環(huán)往復(fù)，直至l=0(或者小于指定誤差)。

將以上過程進(jìn)行模擬，得如圖9所示的流程圖。注意：這里假定浮體不動(dòng)，N不斷調(diào)整，直至收斂。

圖9 外層迭代法流程圖(模擬浮體旋轉(zhuǎn))

圖9中，1)指定排水體積V和水平面法向N0，利用圖7所示的內(nèi)層迭代法求水平面方程，然后分析每個(gè)四面體和水平面方程的關(guān)系，最后匯總計(jì)算浮心B0；

2)當(dāng)水面法線為N0時(shí)，計(jì)算通過B0的浮力作用線和通過G的重力作用線之間的距離l0；

3)根據(jù)浮心B0、重心G和水平面法向N0，計(jì)算浮體旋轉(zhuǎn)軸Axis=B0G?N0；

4)N1=水平面法向量N0繞軸Axis旋轉(zhuǎn)δθ1后的法向量(假設(shè)浮體不動(dòng)，水平面法向量N0旋轉(zhuǎn))；

5)流程運(yùn)行至此說明浮體繞Axis旋轉(zhuǎn)，無論旋轉(zhuǎn)角多小，也無法減小重力和浮力作用線之間的距離。其原因在于浮體重心過高，正在傾覆。浮體將繞Axis旋轉(zhuǎn)1°。

3 算例

為驗(yàn)證文中方法，對某半潛式海洋平臺模型(圖10)進(jìn)行浮態(tài)試驗(yàn)，然后使用商業(yè)軟件Maxsurf和文中方法進(jìn)行對比。平臺主尺度如表1所示。

圖10 半潛式海洋平臺模型(單位：mm)

表1 半潛式海洋平臺模型主尺度

平臺構(gòu)件包含：2個(gè)下浮體(如圖11所示，每個(gè)下浮體由1個(gè)長方體和2個(gè)三棱柱組成)、2個(gè)橫撐(如圖12所示)和4根立柱(如圖13所示)。6個(gè)長方體、4個(gè)三棱柱和2個(gè)圓柱體又可進(jìn)一步分解為多個(gè)四面體，具體方法參考文獻(xiàn)[7]。

圖11 下浮體(單位：mm)

圖12 橫撐(單位：mm)

圖13 立柱(單位：mm)

3.1 浮態(tài)實(shí)驗(yàn)

圖14所示為模型在實(shí)驗(yàn)現(xiàn)場的初始正浮狀態(tài)。因?yàn)橹谱髡`差，需在模型頂部安裝平衡塊。模型空載質(zhì)量M=12.56 kg包含了平衡塊，即圖14圓圈內(nèi)的片狀壓鐵。平衡塊使模型能夠正浮在水面上，最終保證模型的重心恰好位于中縱剖面和中橫剖面的交線上。將安裝好平衡塊并在水中保持正浮狀態(tài)的模型從水中取出，晾干水后，在岸上通過一定的方法測出平臺的重心高度為187.5 mm。

圖14 初始正浮狀態(tài)

再次將平臺置于水中，將0.232 kg的壓載(如圖15所示)布置在平臺頂部不同的位置：工況1為橫傾(如圖16所示)，壓載位于2、3號柱之間；工況2為縱傾(如圖17所示)，壓載位于3、4號柱之間；工況3為任意傾斜(如圖18所示)，壓載位于1號柱上方。平臺+壓載的重心數(shù)據(jù)列于表2。3種工況下，測量各支柱標(biāo)尺(如圖17所示)位置的吃水，列于表3。根據(jù)文中提出的算法，編程(下載地址：www.huagongchuanhai.cn/tetrahedron)計(jì)算平臺3種工況下各個(gè)支柱的吃水，計(jì)算值也列于表3。

圖15 測量壓載的質(zhì)量

圖16 2、3號柱之間的工況1(橫傾)壓載

圖17 3、4號柱之間的工況2(縱傾)壓載

圖18 1號柱上方的工況3(任意傾斜)壓載

表2 各工況下的重心

表3 立柱標(biāo)尺線(圖17)吃水實(shí)驗(yàn)值和計(jì)算值

如表3所示，實(shí)驗(yàn)值和計(jì)算值的最大偏差=1.9%，最小偏差=0.0%，平均偏差=0.2%。可認(rèn)為文中算法在浮態(tài)計(jì)算方面是可靠的。

此外，為了驗(yàn)證算法的收斂性，將各工況每次外層迭代的l值(重力和浮力作用線之間的距離)列于表4之中，并繪制曲線(如圖19所示)。觀察工況1情況下l變化曲線，可發(fā)現(xiàn)兩點(diǎn)規(guī)律：(1)僅經(jīng)過2次迭代，l迅速從5.6 mm收斂至0.545 mm，收斂幅度達(dá)90.3%；(2)從第11次迭代開始，l即收斂至0.001 mm，達(dá)到工程精度要求。其它兩種工況下的收斂也呈現(xiàn)出類似規(guī)律。

圖19 重力和浮力作用線之間的距離變化曲線(外層迭代)

表4 重力和浮力作用線距離(外層迭代)

3.2 對比測試

為了進(jìn)一步檢測文中算法，針對表2和表3中的工況1(橫傾)，使用商用軟件Maxsurf建模進(jìn)行浮態(tài)計(jì)算。

圖20 Maxsurf模型(約50個(gè)剖面)

Maxsurf計(jì)算得到的浮心高度為81.90 mm，橫傾角=2.0°。

采用文中算法得到水面方程4個(gè)系數(shù)：A= 0，B=-0.031 8，C=0.999 5，D=-222.694，浮心高度為75.35 mm，橫傾角為arcsin(0.031 8)=1.8°。

由以上計(jì)算結(jié)果可知：在浮心高度和橫傾角這兩項(xiàng)指標(biāo)上，Maxsurf與文中方法差異明顯。下面用AutoCAD驗(yàn)證文中得到的水面方程是否滿足精度要求。

如圖21所示，利用文中算法得到的水面方程在AutoCAD中繪出水面，然后對平臺進(jìn)行切割，去除水上部分，最后使用AutoCAD的Massprop命令獲得平臺水下部分的質(zhì)心(即平臺的浮心)坐標(biāo)為(463.000，9.222，75.359)。將文中算法得到的水面法線方向〈0，-0.031 8，0.999 5〉，AutoCAD 得到的浮心坐標(biāo)以及給定的重心坐標(biāo)代入式(4)，可得重力和浮力作用線的距離l=0.009 5 mm。在型寬為800 mm的情況下，不到0.01 mm的誤差說明文中算法得到的水面方程非常接近“真實(shí)解”。

圖21 去除水面以上部分的平臺(AutoCAD建模)

Maxsurf計(jì)算得到的浮心高度為81.9 mm，“真實(shí)”浮心高度約為75.4 mm，8.6%的誤差顯然并不理想。如圖22所示，如將Maxsurf的剖面數(shù)增加到200，計(jì)算結(jié)果75.1 mm趨近于理想值(如表5第2行所示)。不同的剖面數(shù)對Maxsurf精度產(chǎn)生影響的原因如圖23所示：平臺圓柱體是橫向布的，當(dāng)剖面數(shù)=50時(shí)，圓柱體直徑方向上只布置了兩個(gè)剖面(如圖23(a)所示)，顯然不夠；當(dāng)剖面數(shù)=200時(shí)，圓柱體直徑方向上布置了足夠多的剖面，故其體積計(jì)算精度令人滿意(如圖23(b)所示)。但是對于圓柱體附近的長方體而言，沿縱向布置2-3個(gè)剖面即可，沒有必要布置太多剖面(如圖23所示)。另外對于圓柱體，雖然只要沿船長(縱向)布置足夠多的剖面(如圖24所示)就可以達(dá)到必要的精度，但其實(shí)這是一種“暴力”計(jì)算(相當(dāng)精度條件下，Maxsurf計(jì)算時(shí)間為1.92 s，而文中算法只需要0.26 s，計(jì)算機(jī)處理器主頻 3.7 GHz)。更為合理的方法應(yīng)當(dāng)是沿圓周方向布置剖面(如圖25所示)。但是Maxsurf“脫胎”于船舶性能計(jì)算，其浮態(tài)計(jì)算基于“剖面”的數(shù)值積分，弊端在于只能沿一個(gè)“統(tǒng)一”方向(船長方向)設(shè)置積分剖面(如圖24所示)。文中基于四面體網(wǎng)格提出的浮態(tài)算法則不受此限制，可根據(jù)需要對浮體的每個(gè)構(gòu)件進(jìn)行合理地“獨(dú)立”分網(wǎng)(如圖25所示)。因此文中算法更適合分析“凸多面體”組合結(jié)構(gòu)。

圖22 Maxsurf模型(約200剖面)

表5 Maxsurf在不同剖面數(shù)目情況下的計(jì)算結(jié)果

(a)約50個(gè)剖面 (b)約200個(gè)剖面

(a)圓柱體剖切網(wǎng)格 (b)等效棱柱

(a)圓柱體剖切網(wǎng)格 (b)等效棱柱

4 結(jié)語

文中針對傳統(tǒng)浮態(tài)計(jì)算方法應(yīng)用于“凸多面體”組合結(jié)構(gòu)存在的不足，提出了一種新的浮態(tài)計(jì)算方法，編制了程序，開展了物理和數(shù)值實(shí)驗(yàn)研究，并與傳統(tǒng)方法進(jìn)行了初步比較研究。文中完成的主要工作和得到的主要結(jié)論如下：

(1)提出了一種全新的浮態(tài)計(jì)算方法，通過內(nèi)外兩層迭代分別模擬浮體平衡過程中的升沉和旋轉(zhuǎn)運(yùn)動(dòng)，具有明確的物理意義。相對于傳統(tǒng)的矩陣迭代法，避免了繁瑣的水面要素和水下體積要素的計(jì)算；同時(shí)也因?yàn)槠渚哂忻鞔_的物理意義，避免了傳統(tǒng)方法有時(shí)候無法收斂的缺陷。

(2)編制了相關(guān)的程序，同時(shí)針對某半潛浮式平臺進(jìn)行了建模和數(shù)值分析，還展開了實(shí)際的水池物理實(shí)驗(yàn)。為“凸多面體”組合結(jié)構(gòu)的浮態(tài)計(jì)算和測量提供了有價(jià)值的參考。計(jì)算和實(shí)驗(yàn)數(shù)據(jù)表明：文中算法具有良好的精度和收斂性。

(3)利用著名的浮態(tài)計(jì)算軟件Maxsurf對目標(biāo)平臺進(jìn)行了建模和數(shù)值分析，將獲得的關(guān)鍵浮態(tài)參數(shù)與文中方法進(jìn)行了對比分析。對于二者產(chǎn)生分歧的地方，請來第三方CAD軟件作為“裁判”。分析表明，對于“凸多面體”組合結(jié)構(gòu)，文中算法更容易得到較為精確的結(jié)果。另外在達(dá)到同等計(jì)算精度的前提下，文中算法的收斂速度也比Maxsurf要快。

文中工作是對凸多面體組合結(jié)構(gòu)浮態(tài)計(jì)算和測量的初步探索，其計(jì)算效率還有提升空間。利用有限元軟件進(jìn)行四面體網(wǎng)格剖分(圖4)只是權(quán)宜之計(jì)。將來如果能夠仿照“船舶曲面快速構(gòu)建”這一思路，專門針對“浮體四面體網(wǎng)格快速剖分”提出相應(yīng)的算法，相信一定能夠從整體上提升算法的計(jì)算效率。