駱寒冰，楊 宇，謝 芃，季紅葉

（天津大學a.水利工程仿真與安全國家重點實驗室；b.船舶與海洋工程系，天津300072）

0 引 言

合理預報船首砰擊載荷關系船體結(jié)構安全，是船舶結(jié)構設計者及研究人員重點關注的問題之一。砰擊是個非常復雜的水動力現(xiàn)象，涉及到波浪形狀、船體表面幾何形狀、船波相對運動、空氣層和水彈性等等[1]。由于該問題的復雜性，數(shù)值預報時通常要做許多假設和簡化。船首砰擊載荷的預報方法，通?？梢苑譃閮纱箢悾阂皇腔谀筒ㄐ苑椒M船舶在波浪中的運動，得到船波相對運動，再預報所感興趣位置的砰擊壓力載荷響應；二是把所感興趣的位置簡化成二維剖面或者三維艙段結(jié)構，再開展入水砰擊載荷預報研究。前者側(cè)重于工程應用，目的是得到設計砰擊壓力，后者側(cè)重于機理研究，預報砰擊壓力的時空分布特性，本文研究屬于后者。

隨著計算機技術的發(fā)展，顯式有限元技術開始應用于入水砰擊載荷預報。駱寒冰等[2]采用LSDYNA 商業(yè)軟件對二維剖面開展了入水砰擊數(shù)值模擬研究，預報結(jié)果得到了模型試驗的驗證，盡管采用了并行算法，但還是存在計算效率較低的不足。鑒于OpenFOAM[3]軟件的開源性，二次開發(fā)功能強大，近年來，逐漸開始被應用于船舶與海洋工程水動力領域。羅天等[4]基于OpenFOAM 開發(fā)的數(shù)值求解器naoe-FOAM-SJTU，模擬了船模的橫搖運動。李裕龍等[5]數(shù)值模擬了液艙液體的晃蕩現(xiàn)象。不過，國內(nèi)采用OpenFOAM軟件預報砰擊載荷的少見，國際上已經(jīng)開始這方面的嘗試。Vesselin等[6]針對二維楔形體非對稱入水砰擊問題，采用OpenFOAM預報了砰擊壓力及其液面變化歷程。Andrea等[7]針對某三維簡化結(jié)構開展了入水砰擊數(shù)值模擬，預報結(jié)果與采用PIV技術測試的試驗液面升高和PIV重構壓力分布結(jié)果進行了對比。

本文主要針對楔形體艙段入水試驗模型，采用OpenFOAM 軟件，開發(fā)了基于兩相流求解器Inter?DyMFOAM 數(shù)值模擬結(jié)構入水砰擊的程序包，開展入水砰擊數(shù)值模擬工作。無轉(zhuǎn)角、有轉(zhuǎn)角工況分別指的是結(jié)構對稱、非對稱垂直入水情況，用一個模型可以得到不同斜升角情況下的結(jié)果。數(shù)值模擬的砰擊載荷結(jié)果與模型試驗結(jié)果進行了對比分析，討論了算法的計算效率以及計算精度。本文研究的主要目的是為今后預報三維復雜結(jié)構砰擊載荷問題提供合理可靠的程序算法。

1 OpenFOAM 軟件算法

OpenFOAM 是基于C++按照面向?qū)ο缶帉戦_發(fā)的開源軟件包。本文自主開發(fā)實現(xiàn)了基于兩相流求解器InterDyMFOAM計算楔形體砰擊載荷的程序，其核心代碼主要是基于有限體積法求解粘性不可壓N-S方程系統(tǒng)，使用VOF方法對自由液面進行捕捉，預報砰擊動壓力載荷。基于不可壓縮流體N-S方程的氣液兩相流模型基本控制方程如下：

連續(xù)性方程

動量方程

流體體積函數(shù)的輸運方程

式中，ρ 為密度，U 為流速矢量，t 為時間，p 為壓力，μ 是動力粘性系數(shù)，F(xiàn) 為體積力，包括重力和表面張力的作用；C 為流體體積函數(shù)，C=0 表示單元內(nèi)的物質(zhì)為空氣；當C=1 表示單元內(nèi)的物質(zhì)為水；0

關于控制方程求解的離散格式，求解時間項使用歐拉法，梯度項使用高斯線性法，拉普拉斯項使用高斯線性修正法，插值格式使用線性法。離散方程中求解液面時使用3次修正迭代，動壓力采用高次有限元方程的高效AMG 法求解，并利用對角不完全Cholesky 方法光順，使用PIMPLE 算法控制壓力與速度的耦合計算。湍流模型采用RAS雷諾時均模型中的k-ω SST湍流模型。

由于楔形體入水砰擊整個過程物體在運動，引入動網(wǎng)格文件控制物體周圍網(wǎng)格的運動，動網(wǎng)格文件采用六個自由度剛體運動求解器，給楔形體一定初速度自由下落。計算控制文件中計算時間從0到設定時間，每隔一定時間步長輸出一次結(jié)果，并調(diào)用函數(shù)，輸出每個時間步長時候楔形體表面測點的動壓力值和楔形體整體的受力結(jié)果。數(shù)值模擬時利用mpirun將多塊計算域同時進行計算以實現(xiàn)并行計算，提高計算效率。

該程序包主要包括以時間命名的數(shù)據(jù)文件夾，其中初始狀態(tài)就是名為0 的文件夾，里面有alpha.water.txt、p_rgh.txt 和U.txt 等求解器對應所需的參數(shù)文件；其次是constant 文件夾，包括網(wǎng)格信息文件夾、polyMesh 和眾多性能文件；最后是系統(tǒng)控制文件夾system，包括計算控制文件controlDict.txt、離散格式文件fvSchemes.txt和離散方程文件fvSolution.txt等。

網(wǎng)格劃分采用OpenFOAM 自帶的SnappyHexMesh 網(wǎng)格劃分工具，它比Gambit 等網(wǎng)格劃分工具功能更加強大。SnappyHexMesh 可以自動從STL、OBJ 文件生成六面體網(wǎng)格，通過迭代將一個初始網(wǎng)格細化。

2 楔形體艙段計算模型

本文采用的楔形體艙段模型基于駱寒冰等[8]的入水砰擊試驗，模型的主尺寸長×寬×高為1.440 m×1.354 m×0.746 m，模型底部斜升角為20°，模型總重180 kg，模型橫剖面結(jié)構示意圖如圖1。模型底部為加筋板結(jié)構，左右舷各有三根縱骨，兩根肋骨。在加筋板板格中心位置布置了8 個壓力測試點（P1-P8），位置左右對稱；在縱骨和肋骨上布置了26 個應變測試點（S1-S26），在模型的前后對稱位置分別布置了兩個加速度傳感器（A1、A2）。測點布置如圖2。實驗水池的主尺寸長×寬×高為5.00 m×4.40 m×1.30 m，實驗時候水深為1.20 m。模型試驗現(xiàn)場見圖3。采用日本共和KYOWA 公司的壓力傳感器，數(shù)據(jù)采集系統(tǒng)為DH-5922 動態(tài)信號測試分析系統(tǒng)，采樣頻率為10 kHz。

圖1 模型橫剖面結(jié)構示意圖 Fig.1 Mid-section structure of the model

圖2 模型測點布置示意圖Fig.2 Distribution of measuring points on the model

圖3 入水砰擊模型現(xiàn)場 Fig.3 The model for water entry test without roll angle

圖4 無轉(zhuǎn)角楔形體數(shù)值模擬模型示意圖Fig.4 The model in OpenFOAM without roll angle

實驗模型的轉(zhuǎn)角可以變化，無轉(zhuǎn)角對稱入水時，楔形體左右兩側(cè)斜升角都為20°，如果轉(zhuǎn)動5°，那么楔形體左右兩側(cè)非對稱，斜升角分別為15°和25°。根據(jù)該模型的無轉(zhuǎn)角和有轉(zhuǎn)角試驗工況，本文將開展數(shù)值模擬，對比研究入水過程的加速度和砰擊壓力結(jié)果。

3 無轉(zhuǎn)角入水砰擊載荷研究

3.1 數(shù)值模擬網(wǎng)格模型及算法參數(shù)選擇

表1 不同網(wǎng)格數(shù)量的模型參數(shù)Tab.1 Models with different mesh sizes

圖5 無轉(zhuǎn)角楔形體模型內(nèi)部3D網(wǎng)格 Fig.5 3D fluid meshes of the model

圖6 楔形體局部網(wǎng)格Fig.6 Local fluid meshes around the wedge model

圖7 無轉(zhuǎn)角楔形體網(wǎng)格橫剖面圖 Fig.7 Model meshes in transverse direction

圖8 無轉(zhuǎn)角楔形體網(wǎng)格縱剖面圖Fig.8 Model meshes in longitudinal direction

考慮到入水砰擊試驗采用的壓力傳感器直徑為0.006 m，分別建立了三個計算模型，最小網(wǎng)格尺寸為0.01 m、0.006 m、0.005 m，對應的網(wǎng)格總數(shù)分別為157萬、642萬、989萬。針對989萬網(wǎng)格模型，圖5顯示了用中橫剖面、中縱剖面剖出的水、空氣和楔形體模型內(nèi)部三維網(wǎng)格，氣液相交處的網(wǎng)格進行了加密，楔形體周圍區(qū)域的網(wǎng)格尺寸最小，圖6 為楔形體周圍局部網(wǎng)格。圖7-8 顯示了模型網(wǎng)格橫剖面和縱剖面圖。為保證每個網(wǎng)格的庫朗數(shù)小于1，根據(jù)最大速度與最小網(wǎng)格尺寸計算時間步長取0.000 1 s。采用了k - ω SST湍流模型，根據(jù)經(jīng)驗公式，計算參數(shù)k和ω的值，分別取0.02和2。數(shù)值模擬工作站的處理器為Intel Core i7-6700CPU@3.40 GHz，四核八線程，8 MB 三級緩存，內(nèi)存為24 G，操作系統(tǒng)為Ubuntu16.04LTS 64位。表1顯示了三個模型的網(wǎng)格劃分情況，以及整個數(shù)值模擬過程所用時間。

3.2 入水加速度及測點砰擊壓力時域結(jié)果對比分析

對表1中的三個不同網(wǎng)格密度模型數(shù)值模擬結(jié)果與試驗結(jié)果進行了對比（見圖9），考慮到無轉(zhuǎn)角入水實驗左右壓力實驗結(jié)果的對稱性，選取右側(cè)壓力測點結(jié)果進行分析。圖9 中列出了P5、P6、P7 測點位置砰擊壓力峰值對比示意，其中橫坐標“Model 1”、“Model 2”、“Model 3”、“Test”分別表示三個不同網(wǎng)格模型和入水實驗模型?？梢钥闯?，隨著網(wǎng)格尺寸逐漸減小，測點壓力峰值呈現(xiàn)增大趨勢，Model 2 比Model 1 的結(jié)果明顯大，不過，Model 2和Model 3 的模擬結(jié)果差異不大，Model 3 的結(jié)果與實驗值吻合最好。認為Model 3的網(wǎng)格模型可以合理模擬砰擊壓力，接下來，將從時域角度對比分析Model 3模型的模擬結(jié)果。

圖9 不同網(wǎng)格尺寸模型與模型實驗的測點砰擊壓力峰值對比Fig.9 Comparison of slamming pressures peaks between numerical models and tests

圖10 無轉(zhuǎn)角入水加速度時域結(jié)果對比 Fig.10 Comparison of acceleration results

圖11 無轉(zhuǎn)角測點P5點壓力時域結(jié)果對比Fig.11 Comparison of pressure results on P5

圖12 無轉(zhuǎn)角測點P6點壓力時域結(jié)果對比 Fig.12 Comparison of pressure results on P6

圖13 無轉(zhuǎn)角測點P7點壓力時域結(jié)果對比Fig.13 Comparison of pressure results on P7

（2）液面升高沒有到達測點前，砰擊壓力為0，到達測點后，脈沖增大達到峰值，然后再緩慢減少。隨著入水深度增加，測點砰擊壓力峰值呈現(xiàn)減少趨勢。P5 測點數(shù)值模擬、模型實驗壓力峰值分別為101.2 kPa、115.7 kPa，數(shù)值模擬峰值偏小12.5%，不過，峰值出現(xiàn)的時刻都是0.018 3 s；P6 測點數(shù)值模擬、模型實驗壓力峰值分別為72.6 kPa、75.5 kPa，出現(xiàn)時刻分別為0.029 0 s、0.030 6 s，數(shù)值模擬峰值偏小3.8%，出現(xiàn)的時刻略早于實驗結(jié)果；P7測點數(shù)值模擬、模型實驗壓力峰值分別為33.3 kPa、30.9 kPa，出現(xiàn)時刻分別為0.045 7 s、0.049 2 s，數(shù)值模擬峰值偏大7.7%，出現(xiàn)的時刻略早于實驗結(jié)果。

由上述分析可見，入水加速度、各個測點砰擊壓力時域結(jié)果數(shù)值模擬和模型實驗吻合較好。采用本文提出的OpenFOAM 數(shù)值模擬算法，包括參數(shù)選擇、模型網(wǎng)格劃分等，可以合理地模擬無轉(zhuǎn)角楔形體入水砰擊過程。

3.3 入水過程壓力云圖及三維效應分析

對于Model 3數(shù)值模擬的結(jié)果，選取0.02 s、0.04 s、0.06 s三個時刻，給出了壓力云圖（如圖14）。考慮對稱性，選取楔形體長度方向中間位置所在的中橫剖面右半部分，圖（a-1）、（b-1）、（c-1）分別顯示了三個時刻的二維剖面流體壓力云圖及流體射流情況。圖（a-2）、（b-2）、（c-2）分別顯示了三個時刻的楔形體模型下表面的三維壓力云圖?？梢钥闯觯?/p>

（1）隨著入水深度增加，各個二維壓力云圖清楚地顯示了壓力峰值從楔形體下部往上方移動的過程，液面升高處壓力最大，越向上部移動，壓力峰值越小。也可以觀察到局部射流發(fā)展過程，隨著入水深度增加，射流向楔形體上方移動，0.06 s時射流已經(jīng)超過楔形體折角位置。

（2）數(shù)值模擬楔形體表面的砰擊壓力，沿楔形體長度方向具有一定的三維效應，也就是說，楔形體中部位置壓力較大，兩端位置的壓力略小。不過，從三維壓力云圖中可以看出，壓力小的兩端長度占整體長度比例不大。考慮到楔形體長度L = 1.44 m，半寬B = 0.677 m，L/B>2.0，流體流動的三維效應所造成砰擊壓力的三維效應不明顯。實際模型實驗時候，楔形體兩端有首尾封板，有利于減少流體流動三維效應。

圖14 不同時刻砰擊壓力云圖及射流Fig.14 Pressure contours and water jets at different times

為了更清楚地說明三維效應，沿楔形體長度方向取三個橫剖面，位置分別在1/2 L、1/4 L、1/8 L處，對比不同橫剖面的相同高度測點在入水過程中砰擊壓力。入水模型實驗以及第3.2節(jié)數(shù)值模擬的壓力測點位置位于1/2 L 橫剖面處，同樣在1/4 L、1/8 L 橫剖面的相同高度位置取相應測點，各個測點分別標記為P5-1、P5-2、P5-3、P6-1、P6-2、P6-3，測點P5-1、P6-1 對應第3.2 節(jié)中的P5、P6。圖15-16顯示了各個剖面測點的動壓力時域結(jié)果。P5-1、P5-2、P5-3測點數(shù)值模擬峰值分別為101.2 kPa、99.6 kPa、93.7 kPa，P6-1、P6-2、P6-3 測點數(shù)值模擬峰值為72.6 kPa、67.5 kPa、60.5 kPa。可見，在長度方向（1/8 L，7/8 L）區(qū)間范圍，楔形體表面砰擊壓力的三維效應不明顯。

圖15 無轉(zhuǎn)角測點P5點不同剖面時域結(jié)果對比Fig.15 Comparison of pressure results on P5

圖16 無轉(zhuǎn)角測點P6點不同剖面時域結(jié)果對比Fig.16 Comparison of pressure results on P6

4 有轉(zhuǎn)角入水砰擊載荷研究

4.1 有轉(zhuǎn)角入水砰擊數(shù)值模擬網(wǎng)格模型

模型實驗時候，利用一套機構可以控制楔形體的轉(zhuǎn)角。圍繞楔形體模型中上部的轉(zhuǎn)動軸，角度最大可以轉(zhuǎn)動到10°，這樣，不僅可以做無轉(zhuǎn)角對稱入水砰擊模型實驗，也可以做有轉(zhuǎn)角時候的非對稱入水實驗。本部分選取了轉(zhuǎn)動角度5°、下落高度為0.3 m的實驗工況進行數(shù)值模擬研究。

圖17 有轉(zhuǎn)角楔形體數(shù)值模擬示意圖Fig.17 The model in OpenFOAM with roll angle of 5 degrees

圖18 有轉(zhuǎn)角楔形體網(wǎng)格橫剖面圖Fig.18 Model meshes in transverse direction with roll angle

圖19 有轉(zhuǎn)角楔形體網(wǎng)格縱剖面圖Fig.19 Model meshes in longitudinal direction with roll angle

4.2 有轉(zhuǎn)角入水加速度及測點砰擊壓力時域結(jié)果對比分析

（2）隨著入水深度增加，兩側(cè)測點的砰擊壓力峰值各自呈現(xiàn)減少趨勢，即P4>P3>P2，并且P5>P6>p7。由于楔形體左側(cè)斜升角15°，右側(cè)25°，左側(cè)測點壓力峰值大于右側(cè)測點的，P4>P5，P3>P6，P2>P7，并且左側(cè)對應位置測點出現(xiàn)峰值的時間比右側(cè)的早；

圖20 有轉(zhuǎn)角入水加速度時域結(jié)果對比Fig.20 Comparison of acceleration results

（3）左側(cè)壓力測點，P4數(shù)值模擬、模型實驗壓力峰值分別為50.9 kPa、66.8 kPa，數(shù)值模擬峰值偏小23.5%，出現(xiàn)的時刻分別是0.031 7 s、0.031 8 s；P3 測點數(shù)值模擬、模型實驗壓力峰值分別為42.5 kPa、51.2 kPa，出現(xiàn)時刻分別為0.046 2 s、0.048 1 s，數(shù)值模擬峰值偏小16.8%，出現(xiàn)的時刻略早于實驗結(jié)果；P2測點數(shù)值模擬、模型實驗壓力峰值分別為22.8 kPa、25.1 kPa，出現(xiàn)時刻分別為0.065 8 s、0.049 4 s，數(shù)值模擬峰值偏大8.8%，出現(xiàn)的時刻略早于實驗結(jié)果；

圖21 有轉(zhuǎn)角測點P4點壓力結(jié)果對比 Fig.21 Comparison of pressure results on P4

圖22 有轉(zhuǎn)角測點P3點壓力結(jié)果對比Fig.22 Comparison of pressure results on P3

圖23 有轉(zhuǎn)角測點P2點壓力結(jié)果對比 Fig.23 Comparison of pressure results on P2

圖24 有轉(zhuǎn)角測點P5點壓力結(jié)果對比Fig.24 Comparison of pressure results on P5

圖25 有轉(zhuǎn)角測點P6點壓力結(jié)果對比 Fig.25 Comparison of pressure results on P6

圖26 有轉(zhuǎn)角測點P7點壓力結(jié)果對比Fig.26 Comparison of pressure results on P7

（4）右側(cè)壓力測點，P5數(shù)值模擬、模型實驗壓力峰值分別為25.6 kPa、24.9 kPa，數(shù)值模擬峰值偏小2.5%，出現(xiàn)的時刻分別是0.035 1 s、0.035 4 s；P6 測點數(shù)值模擬、模型實驗壓力峰值分別為12.7 kPa、13.2 kPa，出現(xiàn)時刻分別為0.060 6 s、0.063 6 s，數(shù)值模擬峰值偏小4.7%，出現(xiàn)的時刻略早于實驗結(jié)果；P7測點數(shù)值模擬、模型實驗壓力峰值分別為5.4 kPa、6.5 kPa，出現(xiàn)時刻分別為0.102 3 s、0.102 8 s，數(shù)值模擬峰值偏大15.8%，出現(xiàn)的時刻略早于實驗結(jié)果。

由上述分析可見，入水加速度、各個測點砰擊壓力時域結(jié)果數(shù)值模擬和模型實驗吻合較好。采用本文提出的OpenFOAM數(shù)值模擬算法，可以合理地模擬有轉(zhuǎn)角楔形體非對稱入水砰擊過程。

圖27 有轉(zhuǎn)角不同時刻砰擊壓力云圖及射流Fig.27 Pressure contours and water jets at different times with roll angle

4.3有轉(zhuǎn)角入水過程壓力云圖分析

對于有轉(zhuǎn)角數(shù)值模擬的結(jié)果，選取0.04 s、0.06 s、0.08 s 三個時刻，圖27 列出壓力云圖以及砰擊射流情況，各個圖例的含義同3.3節(jié)?？梢钥闯觯?/p>

（1）隨著入水深度增加，二維壓力云圖清楚地顯示了壓力峰值從楔形體下部往上方移動的過程，越向上部，壓力峰值越小；由于左側(cè)斜升角15°小于右側(cè)25°，各個時刻左側(cè)液面升高處砰擊壓力峰值明顯大于右側(cè)的，左側(cè)的射流比右側(cè)的更加明顯，砰擊更加劇烈；

（2）數(shù)值模擬楔形體表面的砰擊壓力，沿楔形體長度方向具有一定的三維效應，不過，流體流動的三維效應所造成砰擊壓力的三維效應不明顯。有轉(zhuǎn)角非對稱入水情況下，砰擊壓力的數(shù)值模擬和模型實驗的結(jié)果具有可比較性。左側(cè)液面升高處的壓力峰值明顯大于右側(cè)的。

5 結(jié) 論

針對無轉(zhuǎn)角和有轉(zhuǎn)角楔形體艙段入水砰擊問題，本文采用OpenFOAM 開源軟件，開發(fā)了入水砰擊的程序包，開展了數(shù)值預報研究工作，并與入水砰擊模型實驗進行了對比分析。得到了以下結(jié)論：

（1）本文提出了有效預報三維艙段模型入水砰擊問題的數(shù)值模擬方法。開發(fā)了基于兩相流求解器InterDyMFOAM 數(shù)值模擬結(jié)構入水砰擊的程序，包括初始狀態(tài)、constant以及系統(tǒng)控制文件夾。采用有限體積法求解粘性不可壓N-S方程系統(tǒng)，使用VOF方法對自由液面進行捕捉，并且引入了動網(wǎng)格技術，使用并行計算方法。采用RAS雷諾時均模型中的k-ω SST湍流模型。

（2）應用該算法數(shù)值模擬了某三維楔形體艙段的入水砰擊過程。分別選擇有轉(zhuǎn)角與無轉(zhuǎn)角的一個典型工況，模型最小網(wǎng)格尺寸選取為0.005 m，時間步長取0.000 1 s。預報了測點壓力以及加速度時域結(jié)果，包括其峰值大小及其出現(xiàn)時刻，與模型實驗結(jié)果吻合較好。討論了入水砰擊的三維效應。數(shù)值預報的砰擊壓力云圖以及射流情況，合理描述了有轉(zhuǎn)角和無轉(zhuǎn)角的入水砰擊過程。反映該算法能夠合理有效地預報三維艙段模型入水砰擊問題，并且該算法具有推廣應用于更加復雜的三維船首砰擊問題的潛力。