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

魏姍姍，孫士艷

(江蘇科技大學(xué) 船舶與海洋工程學(xué)院，江蘇 鎮(zhèn)江 212100)

0 引言

穿浪雙體船最早于20 世紀(jì)80 年代初提出，其航行時(shí)船首頻繁出水入水，由于船體與波浪之間劇烈的相對(duì)運(yùn)動(dòng)，船首處容易發(fā)生砰擊現(xiàn)象。在入水瞬間，結(jié)構(gòu)的運(yùn)動(dòng)引起周圍流體介質(zhì)運(yùn)動(dòng)，流體介質(zhì)同時(shí)對(duì)結(jié)構(gòu)物施加反作用力，此時(shí)結(jié)構(gòu)物將遭受巨大的流體砰擊載荷，這會(huì)對(duì)結(jié)構(gòu)體的結(jié)構(gòu)強(qiáng)度造成嚴(yán)重威脅。本文將穿浪雙體船船首剖面簡(jiǎn)化為2 個(gè)完全相同的二維楔形體。

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

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

1 數(shù)值過程

對(duì)于不可壓縮粘性流體，可采用連續(xù)性方程和雷諾平均方程來建立流場(chǎng)的控制方程組。連續(xù)性方程形式為：

式中，u為速度矢量。

雷諾平均方程可以表達(dá)為：

式中：ν為流體運(yùn)動(dòng)粘性系數(shù)；F為質(zhì)量力；ρ是流體密度；p為壓力。

采用VOF 方法逼近空氣與水的交界面，在模型中不同的流體狀態(tài)共用著一套控制方程組，通過引進(jìn)相體積分?jǐn)?shù)這一變量，實(shí)現(xiàn)對(duì)每一個(gè)計(jì)算單元相界面的追蹤。對(duì)于波物相互作用流場(chǎng)，存在空氣和水2 種流體流動(dòng)，用 ρa(bǔ)ir表示空氣密度，用ρwater表示水密度，流場(chǎng)內(nèi)流體密度可表達(dá)為：

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

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

2 數(shù)值結(jié)果分析

2.1 物理模型和網(wǎng)格劃分

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

圖1 坐標(biāo)原點(diǎn)和邊界條件示意圖Fig.1 Sketch of the origin of coordinates and boundary conditions

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

圖2 計(jì)算域示意圖圖Fig.2 Sketch of the computational domain

圖3 楔形體周邊網(wǎng)格示意圖Fig.3 Sketch of the grid surrounding the wedge

2.2 數(shù)值結(jié)果分析

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

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

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

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

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

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

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

圖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

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

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

2.2.3 雙楔形體以給定速度入水的數(shù)值結(jié)果分析

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

圖10 雙楔形體模型坐標(biāo)原點(diǎn)和邊界條件示意圖Fig.10 Sketch of the origin of coordinates and boundary conditions of the double wedges

γ=45°雙楔形體不同工況下的結(jié)果如圖11～圖14所示。其中圖13 和圖14 分別引入圖6 與之對(duì)應(yīng)的單楔形體數(shù)值結(jié)果。

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

圖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 可知，壓力曲線整體趨勢(shì)和峰值出現(xiàn)的相對(duì)位置幾乎不變，雙楔形體壓力峰值都大于單楔形體壓力峰值。圖13（d）t=0.1 s 的數(shù)值結(jié)果與其余不同，雙楔形體壓力顯著高于單楔形體壓力。這是因?yàn)閷?duì)于單楔形體入水，一部分流體受到擾動(dòng)以后以射流形式沿著物面向上運(yùn)動(dòng)，一部分流體沿著自由面向遠(yuǎn)方流動(dòng)，而對(duì)于雙楔形體，向遠(yuǎn)方流動(dòng)的流體被前方楔形體阻擋，被阻斷的流體只能以更高的速度向上運(yùn)動(dòng)來滿足動(dòng)量守恒定律，因此物面壓力顯著上升?？梢灶A(yù)見，在砰擊發(fā)生的短暫時(shí)間內(nèi)，入水速度越大，入水時(shí)間越長，雙楔形體之間的相互壓力影響便越顯著。當(dāng)時(shí)間足夠大時(shí)，重力的影響變得明顯，這一結(jié)論不一定適用，但是長時(shí)間的流體與結(jié)構(gòu)相互作用問題并不是本課題的研究?jī)?nèi)容。14（d）的差別更大是由于速度越高，流體擾動(dòng)速度越大，流體更快達(dá)到相鄰楔形體，相鄰楔形體壓力也升高的更快。

γ=30°雙楔形體不同工況下的結(jié)果如圖15～圖18所示。其中圖17 和圖18 分別引入圖9 與之對(duì)應(yīng)的單楔形體數(shù)值結(jié)果進(jìn)行對(duì)比。

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

圖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 結(jié)語

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

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

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

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

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