蔡 威，張金鳳

(天津大學 水利工程仿真與安全國家重點實驗室，天津 300072)

在近些年來海洋資源的開發(fā)利用中，浮式結構由于其建造成本低、可預制、可拆卸、施工安裝方便等諸多優(yōu)點，越來越受到工程應用領域的重視。浮式結構在海洋環(huán)境中的水動力特性的相關研究方法主要有實驗研究和數值模擬研究，其中實驗研究具有成本高、時間長、操作復雜的特點，而數值模擬具有效率高的特點。浮式結構和流體相互作用的數值模擬研究問題的本質是一個具有自由表面的流固耦合問題，而且需要模擬物體的大幅度運動。這一系列特點與經典的物塊入水抨擊數值模擬相同，因此本文通過對物塊入水抨擊的數值模擬來驗證數值模型的準確性。

對于兩相流的自由表面追蹤問題，研究者目前使用的主要有Level set方法[1]和VOF方法[2]。Gu等[3]使用Level set方法追蹤自由表面，模擬固體的垂向和斜向入水過程；王平等[4]使用VOF方法對圓柱體的入水和下潛過程進行了模擬。對于流固耦合問題，數值模擬方面有很多的研究成果。胥飛等[5]使用SPH(Smoothed Particle Hydrodynamics)方法對二維船體入水砰擊過程進行了模擬；Khayyer等[6]利用MPS(Moving Particle Semi-implicit)方法模擬了二維波浪對海岸結構物的沖擊力；Johnson等[7]基于IBM(Immersed Boundary Method)方法模擬流動問題中的移動固壁邊界；任安祿等[8]利用ALE(Arbitrary Lagrangian-Eulerian)方法來模擬圓柱繞流渦致振動；Rahman等[9]使用多孔體模型研究了半潛式浮式防波堤的運動；Shen Z[10]等使用重疊網格(Overset grids)方法研究了船只的自推進和機動；Quon等[11]使用重疊網格方法分析了波浪能發(fā)電裝置的性能；Cheng等[12]使用重疊網格方法研究了浮式風機的氣動性能。無網格的SPH方法和MPS方法都存在計算量較大的問題，而目前IBM方法不能較好地解決物體大幅度運動的問題，重疊網格方法可以模擬物體的大幅度運動，計算效率較高且對原有求解器的繼承性好。

入水抨擊是海洋工程和海軍工程中的一個重要課題，通常做法是將物體假設為二維楔形體。Von Karman[13]最早采用附加質量近似代替流體作用來分析入水抨擊問題，提出附加質量法計算水上飛機入水抨擊荷載； Dobrovol[14]假設流體具有自相似性，忽略重力影響，使用自相似法，推導出了定速度入水時抨擊壓強分布的理論解。Mei[15]在Dobrovol的基礎上使用邊界有限元(Boundary Element Method，BEM)離散和求解方程，得到了更為一般情況下的解析解。

本文求解Navier-Stokes方程，在OpenFOAM開源程序的基礎上進行二次開發(fā)， 利用VOF方法追蹤自由表面，使用重疊網格技術，建立了流體(波浪)和浮體結構相互作用的二維數值模型，通過模擬楔形物塊垂向和斜向入水過程，驗證數值模型的準確性，并簡單分析楔形物塊入水規(guī)律。

1 數學模型

1.1 控制方程

本文通過求解有限體積法離散的RANS方程，采用VOF方法追蹤自由表面，來模擬氣液兩相流的運動。流體基本控制方程RANS為

(1)

(2)

式(1)和式(2)分別為連續(xù)方程和動量方程，其中：u為速度；ρ為流體密度；v為運動粘滯系數；S為源項，其表達式為

S=Fσ+ρg=Cκα+ρg

(3)

(4)

式中：Fσ為作用在曲面微元上的表面張力之和；C為表面張力系數；κ為界面處的曲率；α為VOF方法中定義的體積分數，并且滿足相方程

(5)

使用Menter[16]提出的SSTk-ω兩方程紊流模型來閉合RANS方程，該模型綜合了標準k-ω模型在近壁面區(qū)的計算優(yōu)點和標準k-ε在遠場的計算優(yōu)點，其紊動能方程和紊流耗散率方程的表達如下

(6)

(7)

式中：k為紊動能；ω為紊動能耗散率。其他各項參數的表達式如下

(8)

(9)

(10)

(11)

(12)

模型中的系數α、β、σk和σω的計算方法為(用φ代表)

φ=F1φ1+F2φ2

(13)

式中：α1=5/9，α2=0.44；β1=3/40，β2=0.082 8；σk1=0.85，σk2=1；σω1=0.5，σω2=0.856；β*=0.09；a1=0.31。

1.2 重疊網格

圖1 重疊網格示意圖Fig.1 Schematic diagram of overset grids

楔形體的運動模擬采用重疊網格法。在這種方法中，計算區(qū)域網格根據所覆蓋區(qū)域的特點，劃分為多塊具有重疊或者嵌套部分的子網格，如圖1所示流場劃分為背景網格(虛線)、固體劃分為貼體網格(實線)。當物體運動時，貼體網格隨之做無變形的剛體運動。數值計算在各個子網格中分別進行，同時在重疊的區(qū)域通過網格節(jié)點插值來實現計算信息的傳遞。

重疊網格方法的關鍵在于建立域連接信息(DCI， domain connectivity information)，用于子網格間的計算信息傳遞，其處理流程為：

(1)搜索洞單元，找出不參與計算的節(jié)點單元；

(2)為邊緣(邊界)節(jié)點從重疊區(qū)域的另外一套網格中尋找合適的插值貢獻節(jié)點；

圖2 重疊網格方法求解器的計算流程Fig.2 Flow chart of the overset grids solver

(3)根據邊緣(邊界)節(jié)點和插值貢獻節(jié)點的位置關系，求解插值系數；

(4)優(yōu)化重疊區(qū)域，減少計算量。

(14)

式中：φL為邊緣(邊界)節(jié)點上的某一流場計算信息，如壓強、相分數等；αk是其對應的第k個插值貢獻節(jié)點的插值系數；φk是其對應的第k個插值貢獻節(jié)點的流場計算信息值。

如圖1所示為重疊網格示意圖，其中方形空心點為洞點，方形實心點為背景網格的邊緣(邊界)節(jié)點，圓形空心點為物體貼體網格的邊緣(邊界)節(jié)點。

對于定常邊界問題，在求解物理問題之前，只需對DCI的建立進行一次處理，并將其存儲在計算機內存中。對于移動邊界問題，需要在流場計算期間重復這些步驟。圖2給出了所開發(fā)的重疊網格方法求解器的整個計算流程。它與常規(guī)的CFD程序非常相似，區(qū)別在于包含了挖洞、DCI建立和插值。

2 數值模擬參數設置

表1 計算分組Tab.1 The groups of calculation

2.1 計算分組設置

圖3-a所示為二維對稱楔形體定速度入水的示意圖，其中α和β分別為物塊左側和右側與水平面所成夾角，θ為物塊對稱軸與豎直方向所成夾角，l為物塊頂部長度，v和u分別為物塊在Z方向和X方向的分速度。算例設置如表1所示。

2.2 網格劃分

二維模型的背景網格為均分網格，尺寸為X=16 m、Z=16 m，X方向和Z方向的網格數均設為160，單個網格尺寸為0.1 m×0.1 m。楔形體部分的貼體網格設置為漸變網格，尺寸為X=6.2 m、Z=5.4 m，楔形體的頂部長度取2.4 m。為了較為準確地模擬楔形體底部的壓力分布，楔形體底部使用snappyHexMesh進行五級漸變加密，最小網格尺寸為3.125 mm×3.125 mm。同時使用extrudeMesh只對Y方向的網格進行處理，將網格變?yōu)閱螌泳W格，可以顯著減少網格數量，3個算例的計算總網格數分別為194 608、131 800和186 168。圖4給出了算例一網格劃分示意圖。

3-a 二維對稱楔形體3-b 模型示意圖圖3 二維對稱楔形物塊入水示意圖Fig.3 Schematic diagram of two-dimensional symmetrical wedge block entering water

圖4 算例一楔形體周圍嵌套網格示意圖Fig.4 Sketch of nested grid around wedge of case 1

2.3 邊界條件

如圖3-b所示，二維模型邊界設置為：左右壁面采用固壁邊界條件，法向梯度為零；前后面采用側壁邊界條件，法向始終不參與求解；楔形體采用物塊邊界條件，邊界上的通量恒定為零；大氣邊界為自由表面邊界條件；底邊采用底床邊界條件，邊界上的速度恒定為零。

3 結果分析

二維水槽的水深設置為8 m，物塊初始時刻下端點位于水面處且為水槽正中央的位置。

3.1 垂向入水

為了便于分析結果，避免液體密度和速度的影響，引入無量綱化的附加壓強參數Cp，其表達式為

(15)

式中：ρ為水的密度，取1 000 kg/m3；v為物塊入水的速度，取2 m/s；p為物塊入水產生的附加壓強，其表達式為

p=P總-p0

(16)

式中：P總為總壓強；p0為靜水時此處的壓強。

圖5給出了算例1和算例2的垂向入水時楔形體底部壁面上的附加壓強分布。分別給出了0.05 s、0.06 s、0.07 s時的無量綱化壓強分布曲線，結果表明不同時刻的數值模擬結果幾乎一致，滿足自相似性。作為對比，給出了Mei[15]的解析解和Dobrovol[14]的理論解，結果與Mei[15]的解析解吻合良好，與Dobrovol[14]的理論解存在明顯的偏差。Mei[15]的解析解和Dobrovol[14]的理論解均基于勢流理論和自相似性求解，但是Mei[15]的解析解進一步考慮了射流對壁面壓強分布影響，更符合實際。

圖6給出了算例1和算例2的垂向入水時楔形體底部壁面上的附加壓強場分布。結合圖5和圖6，可以看出：(1)最大附加壓強并不一定出現在楔形體下端點處；(2)高于靜水面位置處的附加壓強依然較大。算例1中的楔形體，最大附加壓強位置在下端點處，而且附近的壓強變化較小，壓強大的區(qū)域分布較為集中。算例2中的楔形體，最大附加壓強位置分布在距離下端點一段距離的底部壁面兩側，位置高于靜水面，附近的壓強變化較大。

5-a 算例15-b 算例2圖5 垂向入水楔形體底部壁面的附加壓強參數分布Fig.5 Distribution of additional pressure parameter on the bottom wall of the vertical water entry wedge

6-a 算例16-b 算例2圖6 垂向入水楔形體在0.05 s時的附加壓強場分布Fig.6 Distribution of additional pressure field around the vertical water entry wedge at 0.05 s

7-a u=0.6 m/s,v=-2 m/s7-b u=-0.6 m/s,v=-2 m/s圖7 斜向入水楔形體底部壁面附加壓強參數分布Fig.7 Distribution of additional pressure parameter on the bottom wall of the oblique water entry wedge

3.2 斜向入水

實際情況中，物塊入水會有一定的斜向速度，斜向入水速度會對底部壁面上的壓強分布產生影響。Xu[17]在Mei[15]的基礎上，進一步考慮了入水時的射流對壓強分布的影響，使用邊界元分析了楔形體斜向入水問題。本文模擬了水平向速度分別為0.6 m/s 和 -0.6 m/s，垂直速度為-2 m/s，即水平速度和垂直速度絕對值的比值u/|v|為0.3和-0.3時，算例3的底部壁面在0.05 s時刻的附加壓強分布，結果與Mei[15]的BEM結果吻合良好(如圖7所示)。同時給出了水平與垂向速度絕對值之比u/|v|為0.1、0.5、-0.1和-0.5的情況，其結果如圖8所示。

對比分析不同水平與垂向速度之比的情況，可以看出：(1)附加壓強的極值出現的相對位置幾乎不變；(2)附加壓強的最大值隨著水平與垂向速度之比的變化而出現明顯的變化。對于算例3的楔形體，附加壓強的最大值對應的相對位置X/vt≈-4.14。由于算例3中鍥形體入水角度α=20°和β=40°，楔形體入水時左側更貼近水面，因此入水時，左側水面受到擠壓更大，壓強極大值會出現在左側。當水平速度為0時，無量綱化壓強參數極大值為16.87，位于楔形體底部壁面左側。當楔形體斜向左下方入水時，由于左側的擠壓效應增強，隨著水平向的速度增大，附加壓強的最大值增大。當水平與垂向速度絕對值之比為-0.5時，無量綱化附加壓強參數的極大值為21.46。當楔形體斜向右下方入水時，由于左側的擠壓效應減弱，隨著水平向的速度增大，壓強的最大值減小。當水平與垂向速度絕對值之比為0.5時，無量綱化壓強參數極大值為13.10。

8-a 楔形體斜向左下方入水8-b 楔形體斜向右下方入水圖8 不同u/|v|時斜向入水楔形體底部壁面壓強分布對比Fig.8 Distribution of additional pressure parameter on the bottom wall of the oblique water entry wedge in different u/|v|

4 結論

本文采用開源計算流體力學軟件OpenFOAM的兩相流求解器和重疊網格方法，建立了可以模擬物體大幅度位移運動的流固耦合模型。通過對經典的物塊入水問題的模擬，分析和驗證了楔形體物塊垂向和斜向入水狀況下，物塊底部壁面上的壓強及其極值位置的分布。數值模擬結果與理論分析或已有數值模擬結果吻合良好，同時闡明了楔形體垂向和斜向入水的壓強分布規(guī)律，認為水平速度的增大會改變附加壓強的極值大小，但是不會改變其出現的相對位置。