楊壯滔, 張 濤, 段 浩, 朱 敏
?
基于數值波浪水槽的圓柱浮標運動仿真與分析
楊壯滔, 張 濤, 段 浩, 朱 敏
(中國船舶重工集團公司 第705研究所昆明分部, 云南 昆明, 650106)
研究了不同海況條件下外形及重浮心距對圓柱形浮囊浮標運動響應的影響規律。使用雷諾平均的N-S方程(RANS)和流體體積模型(VOF), 基于勢流粘流相結合的方法, 建立了數值波浪水槽。對比分析了不同物理時間下所造波的波形與理論波形, 結果表明, 該方法的造波精度滿足工程應用需求。同時結合6自由度(6DOF)模型和重疊網格方法, 在不同海況作用下對2種類型浮標的運動特性進行仿真。結果表明, 文中所建立的浮標運動仿真分析方法是正確有效的, 文中所做工作可為圓柱形浮囊浮標設計提供參考。
圓柱形浮標; 數值波浪水槽; 6自由度模型; 運動特性
小型浮標在軍用和民用領域的應用前景日益廣闊。但由于尺寸小, 質量輕, 容易受海浪海流的影響, 產生較為劇烈的運動。浮標在海浪的作用下會產生6自由度運動(six degrees of freed- om, 6DOF), 分別為縱蕩(surge)、橫蕩(sway)、垂蕩(heave)、橫搖(roll)、艏搖(yaw)、縱搖(pitch)運動, 從而對其工作造成不利影響, 因此需要對浮標運動響應進行研究, 以指導其外形設計[1]。
與其他形狀浮標相比, 圓柱型浮標受到波浪激勵的影響最小。因此小型浮標的形體多選用圓柱形[2]。工程上常用Morison方程對浮標在波浪激勵下的運動特性進行仿真分析[3]。1977年, Berteaux等[4]使用該方法研究分析了圓柱形浮標的運動特性。2010年, 曲少春等[2]使用該方法對圓柱形浮標進行運動特性分析得出對該型浮標改進的建議。Morison方程是一種基于勢流理論的方法, 無法求解波浪遇到浮標后發生的繞射、越浪、破碎、爬高和渦旋等現象帶來的影響[5]。
隨著計算機技術的發展, 基于雷諾平均的N-S方程(reynolds averaged navier-stokes, RANS)方法的數值波浪水槽能有效求解波浪遇到物體后發生的繞射等現象帶來的影響, 在船舶耐波性分析領域應用較廣, 也可用于研究圓柱形浮標在海浪激勵下的運動特性研究[6]。在相同波浪激勵作用下, 圓柱形浮標的運動特性受其外形和重浮心距影響較大, 但暫無相關研究報道。因此, 文中將使用數值波浪水槽對2種浮囊類型和不同重浮心距的圓柱形浮標在波浪作用下的運動特性進行對比分析, 總結垂蕩和縱搖規律, 為圓柱形浮標提供參考設計。
對于不可壓非定常粘性流體, RANS方法是目前工程上常用的計算流體力學(computatio- nal fluid dynamics, CFD)方法之一。其中, 空間離散采用有限體積法, 時間推進采用雙時間法。該方法的連續性方程和動量方程可寫為


在數值造波時, 需要跟蹤氣體和液體的界面。因而選用流體體積(volume of fluid, VOF)方法處理氣體和液體的界面跟蹤問題, 該方法跟蹤界面是通過求解相連續方程實現的[8]。且

文中使用CFD方法求解波浪作用下浮標運動響應, 在計算中, 浮標姿態隨時間變化而變化。由于不同姿態下所計算出的壓力分布不同, 為了準確求解浮標在不同姿態下所受到的力和力矩, 需要使用6DOF模型來求解浮標的姿態變化。文中選用重疊網格法來實現計算中的浮標運動[9]。
如圖1所示, 該方法的網格分為重疊區域和背景區域, 背景區域為整個計算域, 重疊區域為包含求解對象的運動區域。在計算時, 2個區域將進行數據傳遞和數值差分, 利用重疊區域相對背景區域發生的運動來反應求解對象發生的運動。在劃分網格時, 需對2個區域交界處網格進行控制, 使得在交界處2個區域的網格尺寸保持一致。
建立運動區域后, 需要使用6DOF模型來實現運動區域在流體作用力下的運動計算。該模型將求解對象視為剛體。將其動量方程和動量矩方程對時間求導, 得


將式(4)和式(5)展開成3個方向可得

求解流程如圖2所示。
文中的數值仿真在一個3D數值水池中進行, 該數值水池以RANS方法和多相流模型為基礎。對小型浮標在波浪作用下的運動進行仿真時, 首先需要生成滿足計算要求的波浪。在數值仿真中, 從機理上常見的造波方法可分為2類: 源造波法和邊界造波法。源造波法主要包括質量源造波和動量源造波, 邊界造波法主要包括推板造波、搖板造波、楔形沖箱造波等模擬物理模型的方法以及邊界流體流動速度函數的直接輸入法[10]。
文中使用的造波方法為邊界流體流動速度函數的直接輸入法。根據線性理論, 波面方程和速度場表達式分別為


若采用純粘性流的計算方法, 數值水槽出口處需要增加阻尼區域以防止波浪在出口邊界處發生反射。由于仿真所需物理時間較長, 出口處的阻力會引起波浪幅值的衰減。因此使用不需要阻力區域的勢流與粘流相結合的方法來進行造波, 該方法計算區域分為勢流方法區域和粘流方法區域2部分。如圖3所示, 勢流區域將粘流區域包裹住, 粘流區域在波浪傳播方向的長度等于2個波長, 寬為1個波長, 浮標模型位于此區域。勢流區域的外邊界到粘流區域與勢流區域的交界的距離也分別為2個波長。因此, 文中所建立的數值造波水池沿波浪傳播方向長度等于6個波長, 寬度等于5個波長。由于出口邊界為流體流動的速度函數, 因此不需要消波阻尼區域來防止出口反射, 且生成的波浪幅值不會因阻尼的存在而隨物理時間衰減。
仿真中使用線性規則波作為波浪輸入。為真實反映海況, 采用有意波高作為規則波的波高輸入。2~5級海況波高周期如表1所示。

表1 不同海況下波浪參數
分別選取以2級海況為目標所造的波形圖與理論值進行比較, 波面方程見式(7), 波長為6.1 m, 波高為0.366 m, 周期約為1.98 s。為驗證在反應浮標運動響應的時間范圍內所造波浪的精度, 分別取第0 s、10 s、30 s和50 s的波形進行對比。
圖4中, 剛開始造波時2種方法所造的波形與理論波形一致, 隨著計算物理時間的推移(見圖5~圖7), 由于純粘流造波方法需要在數值水槽出口處加阻尼, 導致出口處波幅開始衰減。
計算的物理時間越長, 波幅的衰減越大, 甚至出現波長發生改變的狀況。對于勢流與粘流相結合的方法, 由于在出口處不需要消波阻尼, 因此不會使波幅發生衰減(見圖4~圖7), 只會因數值計算帶來較小誤差, 對浮標在波浪作用下運動響應的求解影響較小。因此使用勢流和粘流結合的方法進行數值造波, 精度滿足工程實際需求。
對質量和浮囊排水體積相同, 浮囊外形不同的2種浮標進行仿真。為了控制仿真輸入值, 選用單方向規則波作為輸入波, 由于浮標和水槽是對稱的, 因此選用垂蕩幅值、縱搖幅值和縱搖極值作為浮標響應輸出值。其中: 垂蕩幅值能反映浮標豎直方向運動的位移量和速度; 縱搖幅值能反應浮標定軸轉動的角速度; 縱搖極值能反映浮標定軸轉動的旋轉角度。浮標主視圖如圖8所示, 其中左圖浮標浮囊為扁平形, 簡稱扁平型浮標, 右邊浮標浮囊為細長形, 簡稱細長型浮標。浮標參數如表2所示。建立3套坐標系, 如圖9所示, 分別為浮標坐標系和大地坐標系。規定靜水面與浮標軸線相交的交點為大地坐標系的原點,軸垂直于水平面方向向上,軸方向指向波浪傳播方向,軸垂直于、軸向外。

表2 浮標的物理參數

圖8 浮標主視圖
Fig. 8 Main view of buoy
由于浮標在水面運動使得浮標排水實時變化, 導致其浮心位置隨著浮標運動而發生變化。在仿真中, 將浮標視為剛體且質量和質量分布不發生變化, 因此選用浮標質心作為浮標坐標系的坐標原點O,z軸垂沿浮標軸線方向向上, 當浮標漂浮在靜止水面時,x軸和y軸與大地坐標系的軸和軸指向一致。在運動過程中, 浮標的平移響應由浮標坐標系和大地坐標系的相對位移來表示, 浮標的轉動響應由2套坐標系坐標軸的相對夾角來表示。

輸入波形條件, 對不同浮標進行仿真, 并監測波面圖, 如圖10和圖11所示。直到輸出運動響應值呈周期性變化, 則取30 s的穩定周期內的輸出值作為結果。由于該仿真得到的結果是時勵曲線圖, 為了較為明顯地反映出浮標在波浪影響下的響應情況, 取浮標運動穩定后位移幅值的平均值作為垂蕩響應幅值。
如圖12所示, 浮標的垂蕩運動幅值受海況和浮囊外形影響較大, 總體上隨著海況等級的增加而增加。在相同的海況和重浮心距的條件下, 扁平型浮囊浮標的垂蕩幅度大于細長型浮囊浮標的垂蕩幅度。
圖13中, 扁平型浮標的幅值響應算子(respo- nse amplitude operator, RAO)值更接近1, 而細長型浮標的RAO值小于1, 說明扁平型浮標在垂直方向受波面影響比細長型浮標大。使用理論方法對浮標在規則波作用下的運動進行分析, 進一步驗證仿真所得結論。由于浮標的尺寸長度遠小于波浪波長, 重量較輕, 因此采用Froude-Krylov理論將浮標所受流體的拖曳力和慣性力線性化[4]。
在垂直方向上, 浮標運動方程


如圖14所示, 浮標的縱搖運動幅值受海況、浮囊外形和重浮心距的影響較大, 總體上也隨海況等級的增加而增加。在相同的海況和重浮心距的條件下, 若海況等級較低, 細長型浮標的縱搖幅值低于扁平型的; 若海況等級較高, 則細長型浮標的縱搖幅值高于扁平型的。在相同浮標外形和海況的情況下, 浮標的縱搖幅值不隨重浮心距增大而減小。圖14中, 50 mm重浮心距細長型浮標在2級海況作用下較扁平型浮標小很多。
如圖15所示, 當浮標在1個波峰的作用下發生偏轉后姿態開始回復, 若還未回復到平衡位置時下一個波峰再次使浮標發生偏轉。使得浮標運動幅度較小而極值較高。運動幅值反映了縱搖運動的速度大小, 為反映縱搖時位移量的大小需要通過縱搖運動的極值。如圖16所示, 可以看出,在相同海況和重浮心距的條件下, 扁平型浮標的縱搖響應極值要小于細長型浮標的極值。
同理, 浮標繞質心縱搖的運動方程


扁平型浮囊豎直橫截面積較大, 當浮標偏轉角度一致時, 扁平型浮囊的排水體積變化量大于細長型, 使得扁平型浮標的回復力大于細長型浮標。而扁平型浮囊的直徑大, 相同排水體積下, 沾濕面積小于細長型浮囊, 使得扁平型浮標受到水平方向波浪的拖曳力較小。而當海況等級較小時, 波長較短, 相對于浮囊液面不平, 扁平型浮囊直徑大, 使得浮囊左右排水體積相差較大, 造成的擾動力矩也大于細長型浮標。因此, 海況等級較低時, 細長型浮標縱搖運動的幅值較小。當海況等級較高時, 波長較長, 相對于浮囊液面較平, 浮囊左右排水體積相差較小, 造成的擾動力矩較小, 由于扁平型浮標的回復力大于細長型浮標, 受水流拖曳力小于細長型浮標, 因此海況等級較大時, 扁平型浮標縱搖運動的幅值較小。而且扁平型浮標的回復力對于細長型浮標, 受水流拖曳力小于細長型浮標, 使得其更容易回復到平衡位置, 因此扁平型浮標的縱搖極限值小于細長型浮標。
對于相同浮囊外形的浮標, 在運動中排水體積變化較小, 浮心位置變化較小, 只有通過降低重心來提高重浮心距。若重浮心距越大, 受到的回復力矩越大, 同時拖曳力對重心的力矩也越大, 故增大重浮心距不能有效減小浮標縱搖的幅度。
文中使用勢流粘流相結合的方法進行規則波的數值水池造波, 并選取不同物理時間的波形圖與純粘流方法所造的波形圖進行對比, 得出勢流粘流相結合的方法所造波與理論值對比誤差較小, 可以用于浮標運動的仿真計算。使用該方法造波, 對不同外形和不同重浮心距的圓柱形浮標在不同海況作用下的運動進行仿真, 得到浮標在波浪影響下運動仿真的結果, 總結垂蕩和縱搖規律, 使用理論方法進行分析, 驗證仿真結論可信。結論如下:
1) 浮標的垂蕩運動幅值、縱搖運動幅值和極值整體上隨海況等級的增加而增加;
2) 在相同海況和重浮心距條件下, 細長型浮標的垂蕩運動幅值小于扁平型浮標。
3) 在相同海況和重浮心距條件下, 細長型浮標的縱搖運動姿態角極限值大于扁平型浮標;
4) 相同重浮心距條件下, 低等級海況下細長型浮標的縱搖運動幅值小于扁平型浮標, 而高等級海況下細長型浮標的縱搖運動幅值大于扁平型浮標;
5) 在相同海況和外形的條件下, 重浮心距越大, 浮標的橫搖幅值不一定越大。
文中所得結論可以在浮標耐波性設計中起指導性作用, 探究如何削減浮標運動響應帶來的不利影響將是下一步研究方向。
[1] 戴洪磊, 牟乃夏, 王春玉, 等. 我國海洋浮標發展現狀及趨勢[J]. 氣象水文海洋儀器, 2014, 31(2): 118-121, 125.Dai Hong-lei, Mou Nai-xia, Wang Chun-yu, et al. Deve- lopment Status and Trend of Ocean Buoy in China[J]. Me- teorological, Hydrological and Marine Instruments, 2014, 31(2): 118-121, 125.
[2] 曲少春, 鄭琨, 王英民. 圓柱形浮標運動分析與仿真[J].計算機仿真, 2010, 27(6): 363-367.Qu Shao-chun, Zheng Kun, Wang Ying-min. Analysis and Simulation of Spar Buoy Motion[J]. Computer Simulation, 2010, 27(6): 363-367.
[3] 王樹青, 梁丙臣. 海洋工程波浪力學[M]. 青島: 中國海洋大學出版社, 2013.
[4] Berteaux H O, Goldsmith R A, Schott III W E. Heave and Roll Response Of Free Floating Bodies Of Cylindrical Shape: WHOI-77-12[R]. Leipzig, Germany: International Transport Fouum, 1977.
[5] 孫斌. 波浪作用下玻璃鋼浮標水動力特性的數值分析[D]. 長沙: 長沙理工大學, 2011.
[6] 唐歆. 海洋資料浮標水動力分析及結構研究[D]. 上海: 上海海洋大學, 2012.
[7] 馬崢, 黃少鋒, 朱德祥. 湍流模型在船舶計算流體力學中的適用性研究[J]. 水動力學研究與進展, 2009, 24(2): 207-216.
[8] Ma Zheng, Huang Shao-feng, Zhu De-xiang. Study on Applicability of Turbulence Model in Ship Computational Fluid Dynamics[J]. Chinese Journal of Hydrodynamics, 2009, 24(2): 207-216.
[9] Hirt C W, Nichols B D. Volume of Fluid(VOF) Method for the Dynamics of Free Boundary[J]. Journal of Comp- utational Physics, 1981, 39(1): 201-225.
[10] 趙發明, 高成君, 夏瓊. 重疊網格在船舶CFD中的應用研究[J]. 船舶力學, 2011, 15(4): 332-341.Zhao Fa-ming, Gao Cheng-jun, Xia Qiong. Overlap Grid Research on the Application of Ship DFD[J]. Journal of Ship Mechanics, 2011, 15(4): 332-341.
[11] 方昭昭, 朱仁傳, 繆國平, 等. 基于數值波浪水池的波浪中船舶水動力計算[J].水動力學研究與進展, 2012, 27(5): 515-524.Fang Zhao-zhao, Zhu Ren-chuan, Miao Guo-ping, et al. Numerical Calculation of Hydrodynamic Forces for a Ship in Regular Waves Bsed on Numerical Wave Tank[J]. Chinese Journal of Hydrodynamics, 2012, 27(5): 515-524.
(責任編輯: 楊力軍)
Simulation and Analysis of Cylindrical Buoy Motion Based on Numerical Wave Flume
YANG Zhuang-tao, ZHANG Tao, DUAN Hao, ZHU Min
(Kunming Branch of the 705 Research Institute, China Shipbuilding Industry Corporation, Kunming 650106, China)
The influences of a buoy’s shape and gravity center distance on its motion response in different sea conditions are studied. Using the Reynolds-averaged Navier-Stokes equations(RANS) and the volume of fluid(VOF) model, a numerical wave flume is established via combination of potential flow and viscous flow. The waveform generated by this numerical wave flume and the theoretical waveforms are analyzed for different physical time, and the results show that the accuracy of the generated waveform satisfies engineering application. In addition, the six degrees of freedom (6DOF) model and the overlap grid method are employed to simulate the motion characteristics of two kinds of buoys in different sea conditions. It is shown that the present buoy motion simulation and analysis method is correct and effective. This study may facilitate the design of cylindrical buoys.
cylindrical buoy; numerical wave flume; six degrees of freedom(6DOF) model; motion characteristic
楊壯滔, 張濤, 段浩, 等. 基于數值波浪水槽的圓柱浮標運動仿真與分析[J]. 水下無人系統學報, 2018, 26(3): 207-213.
TJ67; TB71.2
A
2096-3920(2018)03-0207-07
10.11993/j.issn.2096-3920.2018.03.004
2017-10-18;
2017-12-19.
楊壯滔(1993-), 男, 在讀碩士, 主要研究方向為水中兵器總體設計.