基于OpenFOAM的破損船舶橫搖阻尼

歐 珊， 毛筱菲， 劉祖源， 黃天奇， 余澤爽

(武漢理工大學 交通學院， 武漢 430063)

船舶橫搖阻尼具有非線性和復雜性，對船舶運動的預報及穩(wěn)性的衡準有重要意義，也是一直以來備受關注的問題，一是對其準確預報在理論上非常困難，而且其非線性使問題變得非常復雜；二是阻尼對橫搖運動影響巨大，研究橫搖通常都要對阻尼進行討論，因此對橫搖阻尼的預報是橫搖運動研究的基礎.

以往的研究中解決橫搖阻尼問題，一種方法是采用經(jīng)驗公式進行求解，如貝爾登公式，尼古拉夫公式，米勒公式，渡邊公式等，還有IKEDA[1]總結的阻尼公式.基于經(jīng)驗公式，2011年第26屆國際拖曳水池會議(ITTC)給出了橫搖阻尼的數(shù)值估算[2].另一種方法是通過船模水池試驗進行確定，分為船模靜水中自由橫搖衰減和強迫振蕩運動試驗，應用能量法來處理.其中，基于船模靜水中自由橫搖衰減得到橫搖阻尼是目前較為常用的方法，但自由衰減試驗只能給出在共振頻率處的船舶阻尼，缺乏頻率的相關性；強迫振蕩試驗雖能給出多個頻率下的附加質量和阻尼，試驗結果的分析仍然會受子樣的限制，另外船模的強迫振蕩運動對試驗設備、測試系統(tǒng)要求較高，目前仍有許多學者致力于試驗測試系統(tǒng)的研究.

隨著計算流體動力學(CFD)方法在船舶領域的應用日益廣泛，涌現(xiàn)出基于雷諾時均(RANS)方程求解船舶橫搖阻尼的研究.最初在二維船體剖面上的研究較多，如文獻[3]和[4].但橫搖問題的三維流動特征明顯，對船舶進行三維模擬更有必要.目前對三維船舶的橫搖研究逐漸增多，基于RANS方程求解方法，Wilson等[5]對 DTMB5512 船進行了自由橫搖和強迫橫搖模擬，并應用ITTC關于CFD不確定度分析規(guī)程對模擬結果進行了驗證和確認.國內應用RANS法進行船舶橫搖阻尼的三維數(shù)值模擬的研究[6-9]，涉及船型、網(wǎng)格、航速、強迫橫搖或自由衰減以及阻尼數(shù)據(jù)處理方法等各個方面.可見，CFD方法在研究橫搖阻尼中已經(jīng)比較普遍，但針對破損船的橫搖阻尼研究較少.目前針對破損船舶的運動模擬仍有許多是基于勢流理論，在計及黏性阻尼影響時幾乎都采用經(jīng)驗或半經(jīng)驗公式，本文基于此，從破損角度來探討橫搖阻尼的特性.

破損的情況很復雜，一般將船舶的破損分為3類：第一類為破損進水后海水全部充滿艙室，不存在自由液面；第二類為艙室局部被淹，而且不與舷外海水相通；第三類為艙室局部被淹并與舷外海水相連通.其中第二類和第三類較常見，本文主要研究這2類破損情況.

第二類破損為內外水耦合問題.船舶與水艙的耦合運動是極為復雜的流體動力學問題，船舶的運動影響水艙內水的流動，水的流動反過來影響船舶的運動.目前還沒有從根本上解決這個問題，但是基于各種假設已有很多的處理方法，一般是將船和水艙分開進行水動力數(shù)值模擬，包括勢流方法和黏性方法，也有將兩者相結合的方法，然后將水艙產(chǎn)生的力和力矩作為外力作用于船舶上從而耦合船舶的運動，并從時域范疇步進求解[10].但是分開處理的方法作了較多假設，與實際有所出入，不能保證時域中每個時間步長都準確求解船舶的運動，而且對于船舶的橫搖運動模態(tài)，用勢流方法求解兩者的耦合運動仍有不妥.

第三類破損為內外水相通問題.由于艙內水與舷外水相通，在船舶運動過程中艙室進水具有動態(tài)特征，進水的流入流出問題導致這類破損不適合用傳統(tǒng)的準靜態(tài)方法進行模擬，即假定破損船舶的進水達到穩(wěn)定之后將艙內水視為液面水平且質量集中于一點，再進行動力響應分析.這種模型簡化了復雜條件，在實際應用上有優(yōu)勢，但是其自由液面水平的假定并不符合物理情形，基于這樣假定得到的艙內水的運動學和動力學結果是不夠精確的，進而導致船舶運動預報的不準確.

因此，發(fā)展基于黏性流理論的破損水動力方法將成為今后的研究重點.本文在應用OpenFOAM對完整船舶的橫搖阻尼研究基礎上，進一步針對第二類和第三類破損船開展純黏性的橫搖數(shù)值模擬，重點研究破損后的橫搖阻尼.一方面由于破損船的橫搖阻尼的數(shù)值研究目前較少，多為試驗研究；另一方面國際海事組織(IMO)第二代完整穩(wěn)性完成后對新一代基于性能的破損穩(wěn)性衡準也將提上日程，由此開展破損船舶的橫搖阻尼的數(shù)值研究，并結合試驗進行探索.

1 數(shù)值計算模型

1.1 研究對象

以一艘圍網(wǎng)漁船為研究對象，通過對船體結構布置及規(guī)范要求的考慮，最終以船中處附近5號漁艙右艙破損開口為例，開口尺寸及位置信息見表1(以船底為參考平面).其中，5號漁艙的2種破損類型在模型處理方式上有所不同：

根據(jù)第二類破損定義可知，該類破損實際上與完整船舶艙室存在自由液面的情況類似，如載液艙、減搖水艙.由于涉及艙室內外水耦合的復雜問題，直接的數(shù)值模擬存在諸多困難，所以在模型處理上，通過一方柱孔(尺寸為 0.04 m×0.04 m×0.092 m)連通破損艙室內部與甲板外部，人為使計算域網(wǎng)格連續(xù)化，船體外表面仍保持完整光滑，如圖1所示.在與物理模型接近的情況下，這種方法可以將非連續(xù)的多域問題簡化為連續(xù)的單域問題，而不改變問題屬性[11].

根據(jù)第三類破損定義可知，該類破損艙室的內外海水是連通的，船體外表面出現(xiàn)破損開口，但由于船體內外海水互通，在數(shù)值建模時可以將破損開口處的艙室內壁面作為船體外表面，這樣既使船體外表面保持連續(xù)性，也保證數(shù)值模型的單域特征.所以在模型處理上直接將破損艙室內壁定義為船體外表面再對整船進行建模和網(wǎng)格劃分，如圖2所示.

通過CATIA對破損船進行物理建模，將模型文件以STL的格式導入OpenFOAM中進行六面體網(wǎng)格劃分，采用OpenFOAM自帶的網(wǎng)格繪制命令BlockMesh和SnappyHexMesh.由于漁船的模型存在附體舭龍骨，所以在繪制網(wǎng)格時必須要保證舭龍骨的表面不失真，這就要在不增加網(wǎng)格數(shù)量的情況下合理設置網(wǎng)格大小以及對船體表面進行“面加密”，如圖3所示.

表1 圍網(wǎng)漁船5號漁艙破損位置及開口尺寸Tab.1 The parameters of damaged fishing compartment m

圖1 圍網(wǎng)漁船第二類破損模型Fig.1 The second situation of damage

圖2 圍網(wǎng)漁船第三類破損模型Fig.2 The third situation of damage

圖3 圍網(wǎng)漁船船體表面網(wǎng)格Fig.3 Grid of hull surface

1.2 控制方程及參數(shù)設置

基于黏性流理論，以雷諾平均的N-S方程為控制方程:

(1)

OpenFOAM程序庫中的interDyMFoam求解器是在自由面求解器interFoam基礎上調用了網(wǎng)格運動模塊建立的，它是求解一類不可壓縮不可滲透恒溫的兩相流問題的求解器.本文采用基于求解器interDyMFoam建立的waveDyMFoam程序庫進行橫搖的數(shù)值模擬，采用流體體積(VOF)方法捕捉自由液面，運用動網(wǎng)格技術包括自適應重劃分網(wǎng)格的網(wǎng)格運動和網(wǎng)格拓撲變化，不僅可以實現(xiàn)船舶六自由度運動模擬也能實現(xiàn)數(shù)值造波.

本文主要研究靜水中的自由衰減問題，湍流模型采用SSTk-ω模型，壓力場和速度場的耦合迭代求解采用PIMPLE算法，時間離散采用歐拉線性的離散方式，梯度項采用高斯線性差值方式，計算控制設置庫朗數(shù)為5，最大時間步長為 0.01 s，計算步長采用可變運行時間，過程文件輸出間隔為 0.2 s.

另外，針對自由橫搖衰減運動的模擬需討論初始橫搖角的給定方式.通常計算自由衰減會采用2種方案：第一種是給定初始橫傾角，在數(shù)值模擬中，加載初始角的方法就是把船體三維模型旋轉一定的角度，這種方式跟試驗狀態(tài)接近，但存在的缺點是在初始橫傾時刻網(wǎng)格無變形，而當船舶從初始橫傾側運動到另一側時網(wǎng)格變形過大，小角度情況下問題不大，但對大角度橫搖就不太適合；第二種是給定初始角速度，在初始時刻船體正浮，給船體一個橫搖的角速度，使之達到橫傾的目的，優(yōu)點是省時省力，缺點是無法精確控制初始橫搖角度的大小，需要進行試算總結出適合該算例的規(guī)律.在實際工作中，根據(jù)對2種方案大量的數(shù)據(jù)對比發(fā)現(xiàn)，計算的阻尼結果差別并不大，所以為了減少工作量，自由橫搖衰減的模擬均采用第二種加載初始角速度的方式.

1.3 網(wǎng)格討論

在網(wǎng)格劃分中，破損船相比完整船不同的是內部網(wǎng)格的繪制，需要設置合理的邊界層參數(shù)，同時，在動網(wǎng)格迭代時對于不同的進水量選取合理的迭代因子，否則計算結果容易發(fā)散.對于破艙內水深的設置，利用文件setFieldsdict單獨設置艙內的液體初始條件，并在SnappyHexMesh文件中調整面精度數(shù)、邊緣精度以及其他的網(wǎng)格因素.

圖4 圍網(wǎng)漁船第二類破損開口處網(wǎng)格劃分細節(jié)圖Fig.4 Detail diagram of grid generation for the second situation of damage

圖4和5分別為圍網(wǎng)漁船第二和第三類破損開口處網(wǎng)格劃分細節(jié)圖.船模計算域尺寸為-L～3L(x軸方向)，-L～L(y軸方向)，-2L～L(z軸方向)，L為船長.考慮船舶橫搖輻射興波經(jīng)壁面反射對船體運動計算的影響進行消波處理，采用的消波方法是在Mayer消波法上的拓展，其原理是增加松弛函數(shù)來消波.另外，局部細化自由面及船體附近網(wǎng)格，以精確捕捉自由表面和船體附近流場信息.需要指出的是，為了便于網(wǎng)格的劃分，在破損模型中艙室內部結構做了簡化處理，如圖6所示.

圖5 圍網(wǎng)漁船第三類破損開口處網(wǎng)格劃分細節(jié)圖Fig.5 Detail diagram of grid generation for the third situation of damage

通過對圖7三種不同網(wǎng)格方案(從左至右:細網(wǎng)格 4.08×106；中等網(wǎng)格 2.14×106；粗網(wǎng)格 1.27×106)進行數(shù)值計算后，分析橫搖和垂蕩運動對網(wǎng)格的依賴性，如圖8所示.綜合考慮精度和時間因素，選擇中等網(wǎng)格作為本文最終的網(wǎng)格方案，其中第二類破損總網(wǎng)格數(shù)在 2.1×106，第三類破損總網(wǎng)格數(shù)在 2.14×106，網(wǎng)格最小尺寸為 0.012 m.

圖7 三種網(wǎng)格方案Fig.7 Three schemes of computational grid

圖8 三種網(wǎng)格下的橫搖和垂蕩運動數(shù)值對比Fig.8 Comparison of roll and heave motion under three different schemes of computational grid

圖9 y+分布Fig.9 y+ contour

近壁面流動的準確捕捉對有物面邊界限制的湍流模擬至關重要，目前工程應用中大多采用壁面函數(shù)法或低雷諾數(shù)模型來處理近壁區(qū)域流動的求解問題.本文采用低雷諾數(shù)湍流模型SSTk-ω，通過邊界層中數(shù)值求解RANS方程實現(xiàn)近壁面流動的模擬.近壁面的網(wǎng)格法向尺度如圖9所示，可見：船體表面的壁面函數(shù)值y+值在絕大部分區(qū)域約為1，在船中及尾部的小部分區(qū)域最大值約為5.

2 數(shù)值模擬結果與分析

2.1 完整船靜水自由橫搖衰減數(shù)值模擬

基于以上的理論基礎和數(shù)值模型，首先對船舶靜水自由橫搖衰減進行CFD模擬.為研究橫搖阻尼的非線性特征設置系列初始橫搖角，得到系列自由橫搖衰減運動隨時間t的幅值，如圖10和圖11所示.通過橫搖消滅曲線擬合方法計算分析出線性阻尼和非線性阻尼2種成分阻尼(具體方法見文獻[12])，并與基于實驗流體動力學(EFD)的模型試驗結果進行對比，匯總如表2.從中得到如下結論：

(1) 從與試驗結果對比中可以看出自由衰減的數(shù)值模擬結果能較好吻合模型試驗結果，橫搖固有周期幾乎一致，阻尼系數(shù)非常接近，驗證了該方法在預報橫搖阻尼研究上的可靠度.

圖10 完整船自由橫搖衰減運動數(shù)值與試驗對比圖Fig.10 Comparision of numerical simulation and experiment for intact ship’s roll free-decay motion

圖11 完整船不同初始橫搖角的自由橫搖衰減運動數(shù)值對比圖Fig.11 Comparision of numerical simulation intact ship’s roll free-decay motion with different initial roll angles

Tab.2Rolldampingcoefficientandnatureperiodofintactship

初始橫搖角/(°)方法無因次線性阻尼系數(shù)無因次平方項阻尼系數(shù)橫搖固有周期/s5CFD0.01370.00602.2759CFD0.01780.00902.27514CFD0.02230.01502.27520CFD0.03150.01732.25019EFD0.03190.00972.230

(2) 如圖12所示，隨著初始橫搖角度的增大，橫搖阻尼的非線性成分也在不斷地增大，在初始橫搖角度達到15° 及以上的時候，阻尼的非線性不可忽略.因此，在橫搖運動理論預報時宜采用非線性阻尼模型，橫搖的線性阻尼模型只在小角度的情況下才適用.

(3) 從圖13和14(圖中p′為船體表面壓力)中可以看出，在舭龍骨附加阻尼的作用下，船舶的自由衰減速度相比光體船舶更快，同時船舶橫搖固有周期有所增加，可以明顯看到船體舭部產(chǎn)生較大的渦，與實際吻合，體現(xiàn)了CFD方法將舭龍骨等附體一體化計算對橫搖阻尼的研究的重要性和必要性.

圖12 不同初始橫搖角下完整船舶線性與非線性橫搖阻尼系數(shù)變化Fig.12 Comparision of numerical simulation and experiment for intact ship’s roll damping coefficients

圖13 帶舭龍骨及光體船舶自由橫搖衰減運動數(shù)值對比圖Fig.13 Comparision of numerical simulation for the roll free-decay motion of intact ship with and without bilge keel

圖14 完整船舶重心位置處橫截面速度矢量圖(初始橫搖角20°)Fig.14 Velocity vector gragh at the cross section of intact ship’s gravity(initial roll angle is 20°)

2.2 第二類破損船靜水自由橫搖衰減數(shù)值模擬

在對完整船的自由橫搖衰減模擬完成的基礎上，應用OpenFOAM對第二類破損船舶的自由橫搖衰減運動進行模擬，探討橫搖阻尼的變化規(guī)律.通過設置不同初始橫搖角，得到不同橫搖時歷，由于此處破損屬于不對稱破損，所以自由橫搖衰減曲線偏離平衡位置，如圖15所示.橫搖阻尼系數(shù)和周期匯總結果如表3所示，從中得到如下結論：

圖15 不同初始橫搖角的第二類破損船舶自由橫搖衰減運動數(shù)值對比圖Fig.15 Comparision of roll angle calculated by CFD method for the second situation of damage

Tab.3Rolldampingcoefficientandnatureperiodofthesecondsituationofdamage

初始橫搖角/(°)方法無因次線性阻尼系數(shù)無因次平方項阻尼系數(shù)橫搖固有周期/s8CFD0.02480.00101.90011CFD0.03340.00131.9009EFD0.03160.00121.93414EFD0.03710.00131.93218EFD0.03960.00211.93228EFD0.04810.00261.933

(1) 將阻尼系數(shù)數(shù)值結果與試驗結果進行比較，得出數(shù)值結果與試驗結果比較接近，體現(xiàn)出阻尼的增長特性，隨著初始橫搖角的增加，線性和非線性阻尼系數(shù)均有所增加，如圖16所示.

圖17 破損前后橫搖阻尼系數(shù)變化Fig.17 Comparision of damping coefficients between intact ship and the second situation of dmage ship

圖16 不同初始橫搖角下第二類破損船舶線性與非線性橫搖阻尼系數(shù)變化Fig.16 Numerical simulation of roll free-decay motion for the second situation of damage with different initial roll angles

(2) 從圖17中發(fā)現(xiàn)將第二類破損后橫搖阻尼與破損前船舶的阻尼進行比較后，一是總阻尼反而有所增加，二是非線性成份所占比例下降，線性成份增加.其原因可能是第二類破損后艙室形成減搖水艙效果，疊加了水艙阻尼.與此同時，注意到固有周期發(fā)生改變導致阻尼起明顯作用的區(qū)間也相應改變，也就是說，雖然破損后艙室形成了減搖效果，但是橫搖阻尼起明顯作用的區(qū)間發(fā)生的變化以及不對稱破損下的橫搖單側幅值增大等因素都對船舶的橫搖運動產(chǎn)生不穩(wěn)定影響.

第二類破損后船舶自由衰減運動不同時刻的流場細節(jié)如圖18所示，隨著時間的增加，破損一側的流體速度矢量不再處于對稱狀態(tài)，體現(xiàn)了非對稱破艙水的晃蕩對船體運動的耦合影響.

圖18 第二類破損船舶重心位置處速度矢量圖(初始橫搖角10°)Fig.18 Velocity vector gragh at the cross section of the second situation of damaged ship’s gravity (initial roll angle is 10°)

2.3 第三類破損船靜水自由橫搖衰減數(shù)值模擬

對第三類破損圍網(wǎng)漁船的自由衰減運動設置2種初始條件：一種是破損艙室初始水位與船舶吃水持平(初始條件I)；另一種是破損艙室初始水相為零(初始條件 II).這樣設置的原因是：① 受試驗條件限制，自由橫搖衰減試驗均在破損進水達到平穩(wěn)后進行，無法考慮船舶破損發(fā)生時刻的動態(tài)進水過程，為了與試驗比較，數(shù)值模擬設置第一種初始條件；② 為了考察破損的進水過程中艙內流場瞬時變化以及對船舶橫搖運動的影響，數(shù)值模擬設置第二種初始條件.通過分析得到如下結論：

(1) 破損進水過程對整個橫搖阻尼的影響并不大，線性和非線性阻尼誤差均為9%左右；兩種初始條件對整個橫搖固有周期則幾乎無影響，誤差為 0.8%，如表4和圖19所示.

(2) 破損進水過程對整個垂蕩模態(tài)的影響較大，如圖20所示，由于船舶吃水的瞬時變化造成垂蕩在第一個周期內的劇烈振蕩，第一個周期之后的浮態(tài)變化則與未考慮破損過程的情況趨于一致.

(3) 破損進水過程在一個周期內達到平穩(wěn)，從圖21可觀測到進水噴涌、液面翻卷及波動等瞬時現(xiàn)象，并且注意到破損開口的上緣部位水流速度矢量變化劇烈， 體現(xiàn)了艙內進水對空氣的排擠影響.

表4第三類破損橫搖阻尼系數(shù)及固有周期匯總表

Tab.4Rolldampingcoefficientandnatureperiodofthethirdsituationofdamage

初始條件無因次線性阻尼系數(shù)無因次平方項阻尼系數(shù)橫搖固有周期/sI0.08560.00151.898II0.07760.00131.915誤差/%9.39.60.8

圖19 有無考慮破損進水過程的橫搖運動時歷對比曲線Fig.19 Time history of roll motion with different initial conditions

圖20 有無考慮破損進水過程的垂蕩運動時歷對比曲線Fig.20 Time history of heave motion with different initial conditions

圖21 第三類破損船進水過程Fig.21 Dynamic process of water penetration

3 結論

運用OpenFOAM開源代碼平臺實現(xiàn)黏性流中完整和破損船舶自由橫搖衰減運動模擬，探討不同初始橫搖角、不同破損形式以及破損過程對橫搖阻尼的影響，并開展試驗驗證研究，得到如下結論：

(1) 破損船體的橫搖阻尼特性比較復雜.船舶在遭遇第二類破損后，橫搖阻尼較完整船舶阻尼有所增加，其中非線性成份所占比例下降，線性成份增加.第三類破損與第二類破損相比，船舶橫搖阻尼也有所增大.破損發(fā)生后，船舶的橫搖固有周期也會發(fā)生改變，但是破損形式對橫搖固有周期影響不大.

(2) 分析第三類破損船舶的破損進水過程發(fā)現(xiàn)，是否考慮船舶達到穩(wěn)態(tài)前的破損過程對整個橫搖阻尼的影響并不大；對整個橫搖周期則幾乎無影響；但是，對整個垂蕩模態(tài)的影響較大，由于船舶吃水的突變造成垂蕩在第一個周期內的劇烈振蕩，第一個周期之后則與未考慮破損過程的情況趨于一致.因此，在研究船舶橫搖阻尼的問題上，破損過程在一定程度上可以簡化處理.

(3) 從方法上嘗試應用OpenFOAM純黏性代碼對破損船舶在波浪中的流固耦合以及艙內外水耦合等難點問題進行研究，實現(xiàn)了CFD預報完整船、破損開口船以及帶自由液面船舶黏性橫搖阻尼的工程應用，為進一步處理波浪中破損船舶的癱船穩(wěn)性問題研究奠定基礎，也為IMO進行破損穩(wěn)性的新一代衡準研究提供參考.

(4) 船舶靜水中的自由橫搖衰減運動得到的阻尼僅為共振頻率處的阻尼，未來將進一步考慮波浪中的自由橫搖衰減運動，研究不同波浪頻率下的船舶阻尼.