針對凸多面體組合結構的二重迭代浮態(tài)計算方法

劉虓 徐磊 樊天慧 陳超核

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

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

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

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

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

1 預備知識

1.1 凸多面體、坐標系和包圍盒

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

圖1 半潛式海洋平臺

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

圖2 浮體坐標系xyz、全局坐標系xGyGzG和包圍盒

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

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

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

(1)

式中：水平面的法向量N為〈nx，ny，nz〉，面內某點p的坐標為(xp，yp，zp)。

平面方程一般表達式為：

Ax+By+Cz+D=0

(2)

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

(3)

1.3 浮體平衡條件

浮體平衡條件如下：

(Ⅰ)V=P/γ；

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

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

(4)

1.4 基于四面體網格的浮態(tài)計算

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

四面體體積為

圖3 四面體單元

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

(5)

四面體形心坐標為

(6)

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

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

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

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

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

2.1 內層迭代(模擬升沉)

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

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

(7)

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

圖6 浮體包圍盒示意圖

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

不妨令vt點的坐標為

(8)

式中：0≤t≤1。

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

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

(9)

同理，得到t1和tm對應的浮體排水量：

(10)

(11)

圖7 內層迭代法流程圖(模擬浮體升沉運動)

2.2 外層迭代(模擬旋轉)

根據內層迭代法固然能夠令浮體浮力=重量(平衡條件I)，但重力作用線與浮力作用線不一定重合(平衡條件II)。如圖8所示，浮力作用線通過浮心B，方向與水平面法向N相同；重力作用線通過重心G，方向與N相反。令兩條作用線之間的距離為l，如果l≠0，則水平面法向N將會繞某個軸Axis(假設浮體不動)旋轉?？筛鶕沂致菪▌t確定Axis：

(12)

圖8 浮體旋轉軸

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

將以上過程進行模擬，得如圖9所示的流程圖。注意：這里假定浮體不動，N不斷調整，直至收斂。

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

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

2)當水面法線為N0時，計算通過B0的浮力作用線和通過G的重力作用線之間的距離l0；

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

4)N1=水平面法向量N0繞軸Axis旋轉δθ1后的法向量(假設浮體不動，水平面法向量N0旋轉)；

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

3 算例

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

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

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

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

圖11 下浮體(單位：mm)

圖12 橫撐(單位：mm)

圖13 立柱(單位：mm)

3.1 浮態(tài)實驗

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

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

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

圖15 測量壓載的質量

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

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

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

表2 各工況下的重心

表3 立柱標尺線(圖17)吃水實驗值和計算值

如表3所示，實驗值和計算值的最大偏差=1.9%，最小偏差=0.0%，平均偏差=0.2%?？烧J為文中算法在浮態(tài)計算方面是可靠的。

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

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

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

3.2 對比測試

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

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

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

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

由以上計算結果可知：在浮心高度和橫傾角這兩項指標上，Maxsurf與文中方法差異明顯。下面用AutoCAD驗證文中得到的水面方程是否滿足精度要求。

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

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

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

圖22 Maxsurf模型(約200剖面)

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

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

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

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

4 結語

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

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

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

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

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