基于OpenFOAM 的楔形體入水問題研究

魏姍姍，孫士艷

(江蘇科技大學 船舶與海洋工程學院，江蘇 鎮江 212100)

0 引言

穿浪雙體船最早于20 世紀80 年代初提出，其航行時船首頻繁出水入水，由于船體與波浪之間劇烈的相對運動，船首處容易發生砰擊現象。在入水瞬間，結構的運動引起周圍流體介質運動，流體介質同時對結構物施加反作用力，此時結構物將遭受巨大的流體砰擊載荷，這會對結構體的結構強度造成嚴重威脅。本文將穿浪雙體船船首剖面簡化為2 個完全相同的二維楔形體。

Zhao 和Ma 等[1–2]研究了無重力影響時二維楔形體壓力分布的峰值位置，表明底升角大約大于40°時峰值位置由射流根部轉移到楔形體尖端。Chen[3]建立的二維直升機V 形浮標的數值模擬結果也驗證了Zhao 和Ma 等[2]的結論。陳震[4]表明二維楔形體入水壓力峰值的傳播速度隨著入水角度的增大而減小。張健[5]、王玲等[6]、Chen[7]和陳光茂等[8]分別對二維楔形體的入水砰擊問題進行仿真分析，研究壓力峰值與底升角、入水角度和入水速度的關系。Ma[9]模擬了二維楔體垂直和斜入的入水和水動力沖擊過程。Wu 等[10]表明有限高度非對稱二維楔形體斜入水速度增大時，一段較長時間內的重力效應可以忽略。Sun 等[11]首次分析了二維楔形體在波浪中的入水過程，并研究了重力對楔形體入水過程和物面壓力分布的影響。陳翔等[12]表明了楔形體入水底升角越小流體沿物面爬升速度越快，射流越明顯。Zhang 等[13]表明對稱楔形體入水過程中兩側受到的壓力也是對稱的。

本文采用OpenFOAM 開源軟件，建立具有不同底升角單楔形體和雙楔形體在靜水中入水的CFD 數值模擬模型，采用有限體積法（FVM）離散平均雷諾數方程和連續性方程組成的控制方程組，采用流體體積法（VOF）追蹤自由表面。將基于粘流理論的CFD 數值結果與前人基于勢流理論的數值結果進行對比分析，研究粘流方法與勢流方法在處理入水問題過程中的差別，并重點分析重力和粘性在此過程中所起的作用。

1 數值過程

對于不可壓縮粘性流體，可采用連續性方程和雷諾平均方程來建立流場的控制方程組。連續性方程形式為：

式中，u為速度矢量。

雷諾平均方程可以表達為：

式中：ν為流體運動粘性系數；F為質量力；ρ是流體密度；p為壓力。

采用VOF 方法逼近空氣與水的交界面，在模型中不同的流體狀態共用著一套控制方程組，通過引進相體積分數這一變量，實現對每一個計算單元相界面的追蹤。對于波物相互作用流場，存在空氣和水2 種流體流動，用 ρair表示空氣密度，用ρwater表示水密度，流場內流體密度可表達為：

使用SSTk-ω兩方程湍流模型來閉合RANS 方程，SSTk-ω湍流擴散方程為：

式中：k為湍流動能；為速度的水平分量和垂直分量；ω為湍流特征頻率；Γk和 Γω分別為k和 ω的有效擴散項；Gk和Gω分別為k和 ω的生成項；Yk和Yω分別為k和 ω的耗散項；Dω為橫向擴散項；Sk和Sω為自定義源項。

2 數值結果分析

2.1 物理模型和網格劃分

如圖1 所示，笛卡爾坐標系x軸方向水平向右，y軸方向豎直向上，初始時刻楔形體尖端位于原點且靜置于靜水面上。流體以速度W自下而上流動，等價于楔形體以速度W沿y軸負方向運動。楔形體固定高度h為5 m，底升角 γ分別采用45°和30°，邊界條件見圖1。

圖1 坐標原點和邊界條件示意圖Fig.1 Sketch of the origin of coordinates and boundary conditions

使用blockMesh 工具繪制背景網格，計算域尺度為寬度6 m，厚度0.5 m，高度5 m，其中空氣域為2 m，水深為3 m，如圖2 所示。為了較為準確地模擬楔形體的壓力分布，使用snappyHexMesh 網格生成器分裂加密楔形體周邊網格，最小網格尺度為2.6 mm×2 mm，如圖3 所示。

圖2 計算域示意圖圖Fig.2 Sketch of the computational domain

圖3 楔形體周邊網格示意圖Fig.3 Sketch of the grid surrounding the wedge

2.2 數值結果分析

為了便于分析問題，在下面數值結果中，如無特殊說明，本文將采用垂直速度W、重力加速度g和水密度 ρ對數值結果進行無量綱化，將長度、時間t、速度、壓力P的量綱分別用W2/g，W/g，W和ρW2表示。

2.2.1 45°底升角單楔形體入水

γ=45°單楔形體不同工況下的結果如圖4～圖6 所示。其中圖4（d）和圖5（d）給出了無量綱自由液面形狀，圖6 的相似解引自Sun 等[11]的基于勢流理論的45°底升角無限高度楔形體勻速入水時的壓力分布曲線。

圖4 45°楔形體3 m/s 入水自由液面形狀Fig.4 The free surface profile of 45° wedge with 3 m/s

圖5 45°楔形體6 m/s 入水自由液面形狀Fig.5 The free surface profile of 45° wedge with 6 m/s

圖6 重力影響下45°楔形體表面壓力分布Fig.6 The pressure on the 45° wedge surface with the gravity effect

由圖4 和圖5 可知，當入水時間較短時，自由液面沿物面快速爬升。隨著入水時間的增加，射流尖端到達楔形體頂部，射流沿楔形體頂部飛出，形成了射流分離現象。入水時間相同，重力的作用效果也基本相同，而圖5 比圖4 射流更高的直接原因是其對應的楔形體入水速度更高，被楔形體擾動的射流獲得了更高的分離速度。不考慮重力和射流分離影響時，楔形體的入水問題是自相似的，而圖5 中t=0.1s時刻的自由面無量綱曲線與其余2 個時刻產生偏差，圖5 對應的楔形體入水距離更大，重力影響變明顯，射流無法無限期維持。

圖6（a）中t=0.01 s 的數值結果與其余3 條曲線有明顯差別。這是由于采用VOF 方法確定的氣液交界面是一條密度逐步過渡的細帶狀區域，這將導致采用CFD 方法計算的壓力在入水瞬間計算的不夠精準。速度提升了1 倍后，相同時刻楔形體入水距離已經達到一定距離，楔形體表面壓力不再受自由液面影響，見圖6（b）。但同時在0.1 s，物面壓力又急劇降低，這是由于相似解沒有考慮重力，而基于粘流理論的CFD 結果考慮了重力，隨時間增加重力的作用越來越明顯，導致物面壓力急劇下降，并且遠離相似解。此外，所有CFD 壓力數值結果都是小于相似解的，前面描述的重力只是其中的一個次要原因，另外一個重要原因是砰擊發生時物面壓力梯度較大，CFD 方法由于網格尺度限制在處理大壓力梯度問題時存在一定困難。

2.2.2 30°底升角單楔形體入水

γ=30°單楔形體不同工況下的結構如圖7～圖9 所示。其中圖7（d）和圖8（d）給出了無量綱自由液面形狀，圖9 的相似解引自許國冬[14]的基于勢流理論的30°底升角無限高度楔形體勻速入水時的壓力分布曲線。

圖7 30°楔形體3 m/s 入水自由液面形狀Fig.7 The free surface profile of 30° wedge with 3 m/s

圖8 30°楔形體6 m/s 入水自由液面形狀Fig.8 The free surface profile of 30° wedge with 6 m/s

圖9 重力影響下30°楔形體表面壓力分布Fig.9 The pressure on the 30° wedge surface with the gravity effect

相同入水速度相同時刻，γ=30°楔形體入水瞬間引起的射流均比γ=45°楔形體的更加細長，符合許國冬[14]和Sun[11]等的研究成果。當楔形體底升角越小，流體質點流動速度更快，射流更長。圖8（d）在t=0.1 s 時刻分離射流達到較大的高度和體積，由于自由射流獲得的慣性向上速度來源于楔形體，當它脫離楔形體以后所受的作用力只有重力，使得自由射流的慣性速度變小。越接近自由射流尖端，脫體時間越長，受到的重力影響就越明顯，因此射流尖端附近的流體質點速度首先降低，出現向右下方彎曲的趨勢。

除上一個算例對壓力曲線的分析原因，影響壓力另外一個重要原因是射流分離和重力，見圖9（b）的0.1 s 壓力結果，附體射流層的0 壓力線在曲線中消失，且壓力值降低。

2.2.3 雙楔形體以給定速度入水的數值結果分析

如圖10 所示，雙楔形體A和B兩個頂點的距離設置為0.5 m，其余設置與圖1 一致。

圖10 雙楔形體模型坐標原點和邊界條件示意圖Fig.10 Sketch of the origin of coordinates and boundary conditions of the double wedges

γ=45°雙楔形體不同工況下的結果如圖11～圖14所示。其中圖13 和圖14 分別引入圖6 與之對應的單楔形體數值結果。

由圖11 可知，自由面沿物面快速爬升，形成粘附于物面的附體射流，最終脫離物體頂部，形成一小段自由射流，在雙楔形體中間，自由液面連同自由射流形成一個開口氣腔。圖12 中雙楔形體入水速度提高一倍，附體射流較快地演變為自由射流，t=0.1 s時雙楔形體之間的液面形式為一個凸起的水柱。這是由于入水速度足夠大時，楔形體擾動的自由液面范圍擴大，使得雙楔形體內的兩股液面融合，自由液面以水柱形式抬升。

圖11 45°雙楔形體3 m/s 入水自由液面形狀Fig.11 The free surface profile of 45° double wedges with 3 m/s

圖12 45°雙楔形體6 m/s 入水自由液面形狀Fig.12 The free surface profile of 45° double wedges with 6 m/s

由圖13 和圖14 可知，壓力曲線整體趨勢和峰值出現的相對位置幾乎不變，雙楔形體壓力峰值都大于單楔形體壓力峰值。圖13（d）t=0.1 s 的數值結果與其余不同，雙楔形體壓力顯著高于單楔形體壓力。這是因為對于單楔形體入水，一部分流體受到擾動以后以射流形式沿著物面向上運動，一部分流體沿著自由面向遠方流動，而對于雙楔形體，向遠方流動的流體被前方楔形體阻擋，被阻斷的流體只能以更高的速度向上運動來滿足動量守恒定律，因此物面壓力顯著上升。可以預見，在砰擊發生的短暫時間內，入水速度越大，入水時間越長，雙楔形體之間的相互壓力影響便越顯著。當時間足夠大時，重力的影響變得明顯，這一結論不一定適用，但是長時間的流體與結構相互作用問題并不是本課題的研究內容。14（d）的差別更大是由于速度越高，流體擾動速度越大，流體更快達到相鄰楔形體，相鄰楔形體壓力也升高的更快。

γ=30°雙楔形體不同工況下的結果如圖15～圖18所示。其中圖17 和圖18 分別引入圖9 與之對應的單楔形體數值結果進行對比。

不同于圖11 的氣腔開口，圖15 的氣腔基本封閉了，原因是小底升角導致射流速度更高，射流長度變得更長，使得氣腔更容易封閉，所以圖16（c）中的雙楔形體中間的拱形射流也具有更高的高度。與上一個算例類似，在圖17 和圖18 中，由于雙楔形體之間的相互影響需要時間，時間越長，單楔形體和雙楔形體壓力差別越大。與上一個算例不同的是，本算例楔形體底升角更小，受到擾動的流體運動得更快，因此本算例單楔形體和雙楔形體壓力之間的距離較前一個算例要大。由此可以得出另外一個結論，在楔形體間距不變的條件下，小底升角會使楔形體之間的影響時間提前，影響幅度擴大。

圖15 30°雙楔形體3 m/s 入水自由液面形狀Fig.15 The free surface profile of 30° double wedges with 3 m/s

圖16 30°雙楔形體6 m/s 入水自由液面形狀Fig.16 The free surface profile of 30° double wedges with 6 m/s

圖17 3 m/s 勻速入水30°雙楔形體表面壓力分布Fig.17 The pressure on the 30° double wedges surface of entry with constant velocity 3 m/s

圖18 6 m/s 勻速入水30°雙楔形體表面壓力分布Fig.18 The pressure on the 30° double wedges surface of entry with constant velocity 6 m/s

3 結語

本文采用基于粘流理論的CFD 方法，模擬重力作用下不同底升角單楔形體和雙楔形體以給定速度在靜水中的入水過程，得到如下結論：

1）當入水時間較短時，不同時刻的無量綱自由液面曲線基本重合，隨著時間的增加，射流從物體頂部分離，自由射流在重力作用下發生彎曲向下的現象。

2）當入水時間較短時，排除初始階段的不準確壓力，不同時刻的物面無量綱壓力曲線是基本重合的，并且與基于勢流理論的自相似解非常接近，隨著時間的增加，重力的作用變得顯著，壓力在重力作用下降低。

3）在相同入水速度和相同底升角條件下，雙楔形體物面壓力大于單楔形體物面壓力，并且在砰擊發生的短暫時間范圍內，剛開始影響不大，從某一時刻開始，受一個楔形體擾動的流體到達相鄰楔形體以后，相鄰楔形體壓力差值瞬間變得顯著。此外，若提升速度，被楔形體擾動的流體速度加快，相鄰楔形體受到的壓力影響更大。

4）當降低楔形體底升角時，相當于增加楔形體周圍被擾動的流體速度，因此相鄰楔形體受到的壓力也提升的更快。