基于重疊網(wǎng)格方法的楔形體入水抨擊數(shù)值模擬

蔡 威，張金鳳

(天津大學(xué) 水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室，天津 300072)

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

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

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

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

1 數(shù)學(xué)模型

1.1 控制方程

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

(1)

(2)

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

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

(3)

(4)

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

(5)

使用Menter[16]提出的SSTk-ω兩方程紊流模型來閉合RANS方程，該模型綜合了標(biāo)準(zhǔn)k-ω模型在近壁面區(qū)的計(jì)算優(yōu)點(diǎn)和標(biāo)準(zhǔn)k-ε在遠(yuǎn)場(chǎng)的計(jì)算優(yōu)點(diǎn)，其紊動(dòng)能方程和紊流耗散率方程的表達(dá)如下

(6)

(7)

式中：k為紊動(dòng)能；ω為紊動(dòng)能耗散率。其他各項(xiàng)參數(shù)的表達(dá)式如下

(8)

(9)

(10)

(11)

(12)

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

φ=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 重疊網(wǎng)格

圖1 重疊網(wǎng)格示意圖Fig.1 Schematic diagram of overset grids

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

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

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

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

圖2 重疊網(wǎng)格方法求解器的計(jì)算流程Fig.2 Flow chart of the overset grids solver

(3)根據(jù)邊緣(邊界)節(jié)點(diǎn)和插值貢獻(xiàn)節(jié)點(diǎn)的位置關(guān)系，求解插值系數(shù)；

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

(14)

式中：φL為邊緣(邊界)節(jié)點(diǎn)上的某一流場(chǎng)計(jì)算信息，如壓強(qiáng)、相分?jǐn)?shù)等；αk是其對(duì)應(yīng)的第k個(gè)插值貢獻(xiàn)節(jié)點(diǎn)的插值系數(shù)；φk是其對(duì)應(yīng)的第k個(gè)插值貢獻(xiàn)節(jié)點(diǎn)的流場(chǎng)計(jì)算信息值。

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

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

2 數(shù)值模擬參數(shù)設(shè)置

表1 計(jì)算分組Tab.1 The groups of calculation

2.1 計(jì)算分組設(shè)置

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

2.2 網(wǎng)格劃分

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

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

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

2.3 邊界條件

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

3 結(jié)果分析

二維水槽的水深設(shè)置為8 m，物塊初始時(shí)刻下端點(diǎn)位于水面處且為水槽正中央的位置。

3.1 垂向入水

為了便于分析結(jié)果，避免液體密度和速度的影響，引入無量綱化的附加壓強(qiáng)參數(shù)Cp，其表達(dá)式為

(15)

式中：ρ為水的密度，取1 000 kg/m3；v為物塊入水的速度，取2 m/s；p為物塊入水產(chǎn)生的附加壓強(qiáng)，其表達(dá)式為

p=P總-p0

(16)

式中：P總為總壓強(qiáng)；p0為靜水時(shí)此處的壓強(qiáng)。

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

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

5-a 算例15-b 算例2圖5 垂向入水楔形體底部壁面的附加壓強(qiáng)參數(shù)分布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時(shí)的附加壓強(qiáng)場(chǎng)分布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 斜向入水楔形體底部壁面附加壓強(qiáng)參數(shù)分布Fig.7 Distribution of additional pressure parameter on the bottom wall of the oblique water entry wedge

3.2 斜向入水

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

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

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

4 結(jié)論

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