基于FLUENT的內(nèi)孤立波質(zhì)量源造波方法*

姜 海 郭海燕 張 林 王 偉

(中國海洋大學工程學院 青島 266100)

海洋內(nèi)波是指在海水穩(wěn)定層化的海洋中產(chǎn)生的、最大振幅出現(xiàn)在海洋內(nèi)部的波動。海洋內(nèi)波的產(chǎn)生應具備兩個條件: 一是海水密度穩(wěn)定分層, 二是要有擾動能源, 缺一不可。由于海水的密度分布經(jīng)常處于不均勻狀態(tài), 因此海洋內(nèi)波是一種比較普遍的現(xiàn)象(杜濤等, 2001; 柯自明等, 2009; Xu et al, 2012)。其中,內(nèi)孤立波是一種強非線性的潮成內(nèi)波, 非線性與色散效應的平衡使得其在傳播過程中可以保持穩(wěn)定的波形和不變的速度(方欣華等, 2005)。內(nèi)孤立波通常具有很大的振幅和波長, 在傳播過程中攜帶有巨大的能量, 對海洋結(jié)構(gòu)物構(gòu)成很大的威脅(Roberts,1975; Osborne et al, 1980; 張波等, 1998)。

現(xiàn)在廣泛使用的內(nèi)孤立波理論有 KdV理論、mKdV理論、eKdV理論及MCC理論等(Michallet et al,1998; Choi et al, 1996, 1999; Helfrich et al, 2006)。近些年來, 隨著計算機技術(shù)及計算流體力學的發(fā)展, 采用計算流體力學(CFD)的方法研究內(nèi)孤立波的生成傳播及其與結(jié)構(gòu)物的相互作用, 已成為重要的研究方法之一。內(nèi)孤立波數(shù)值模擬常用的方法有兩類——仿物理造波法和純數(shù)值造波法。仿物理造波法包括平板拍擊法、雙推板法、搖板法等, 純數(shù)值造波法包括速度入口法和源項造波法。付東明等(2009)結(jié)合KdV理論發(fā)展了雙推板數(shù)值造波方法; 陳鈺等(2009)基于KdV理論, 采用速度入口法, 建立了可有效模擬弱非線性海洋內(nèi)孤立波的二維分層數(shù)值水槽; 高原雪等(2012)基于MCC理論, 采用速度入口法, 建立二維內(nèi)孤立波數(shù)值水槽, 實現(xiàn)了振幅可控的內(nèi)孤立波數(shù)值模擬; 李水娟等(2012)基于KdV理論, 采用速度入口法建立了三維內(nèi)孤立波分層數(shù)值水槽; 王旭等(2014)結(jié)合 eKdV理論,從理想流體完全非線性歐拉方程出發(fā), 發(fā)展了內(nèi)孤立波質(zhì)量源數(shù)值造波方法, 所得結(jié)果跟理論與實驗吻合較好; 王偉等(2016)基于 KdV、mKdV理論, 采用平板拍擊法進行內(nèi)孤立波數(shù)值模擬,并詳細分析了內(nèi)孤立波波致流場變化。

自從 Brorsen等(1987)基于“水下爆炸”原理, 首次提出基于勢流理論的適合于邊界積分方程法的質(zhì)量源函數(shù)造波方法后, 質(zhì)量源方法已被廣泛應用于表面波數(shù)值模擬(Lin et al, 1999; 劉秀麗等, 2015; 周玲玲等, 2015), 但目前采用質(zhì)量源法模擬內(nèi)孤立波的論文很少。王旭等(2014)采用線源形式的質(zhì)量源, 結(jié)合eKdV理論成功模擬了內(nèi)孤立波。本文采用兩個點源形式的質(zhì)量源, 基于FLUENT軟件, 從黏性不可壓縮流體的Navier-Stokes方程出發(fā), 結(jié)合KdV、eKdV理論建立有限水深兩層流體中的內(nèi)孤立波質(zhì)量源數(shù)值模擬方法, 并將所得結(jié)果與理論和實驗相對比, 驗證方法的可行性。

1 數(shù)學模型

1.1 造波方法

流體力學中的源方法是在流體力學控制方程中加入特定的源項, 以此激發(fā)產(chǎn)生所需要的流場。源項是一個廣義量, 它代表了那些不能包括到控制方程的非穩(wěn)態(tài)項、對流項與擴散項中的所有其它各項之和,包括質(zhì)量源項、動量源項、能量源項和湍流源項。其中, 質(zhì)量源方法是在連續(xù)性方程中加入附加質(zhì)量源項, 通過對附加質(zhì)量源項進行特定的設置, 以此來產(chǎn)生目標流場。Brorsen等(1987)對質(zhì)量源方法進行了詳細的數(shù)學物理推導, Kim等(2010)從雷諾輸運定理出發(fā), 詳細推導了質(zhì)量源表達式。

本文用質(zhì)量源法建立二維內(nèi)孤立波數(shù)值水槽,如圖1所示。上層流體的深度為h1, 密度為ρ1, 下層流體的深度為h2, 密度為ρ2, 建立如圖1中所示的坐標系oxz, ox軸位于兩層流體尚未擾動時的分界面上,oz軸垂直向上。在左邊界處, 上下層流體中各放置一個點源形式的質(zhì)量源, 上層中的源 Ω1的強度為S1( x, z, t), 下層中的源 Ω2的強度為 S2(x, z, t)。源 Ω1不斷釋放質(zhì)量, 源 Ω2不斷吸收質(zhì)量, 從而產(chǎn)生下凹型的內(nèi)孤立波。

圖1 內(nèi)孤立波數(shù)值水槽示意圖Fig.1 Scheme of the numerical wave tank for generating internal solitary wave

在質(zhì)量源區(qū)域, 加入源項后, 流體的控制方程變?yōu)槿缦滦问?

式中, S （ x, z , t）為源強, μ為動力黏度系數(shù)。

假定質(zhì)量源引起的流體質(zhì)量的增減對應入射的孤立波波面的變化, 則對于源Ω1、Ω2有:

式中, C為內(nèi)孤立波相速度, η（ t ）為內(nèi)孤立波波面。

質(zhì)量源為點源, 在質(zhì)量源區(qū)域內(nèi)均勻分布, 可認為質(zhì)量源強度僅為時間 t的函數(shù), 即 Si( x, z, t) = Si( t),則質(zhì)量源Ω1、Ω2的源強的表達式如下:

式中,1A、2A分別為源Ω1、Ω2的面積。

1.2 消波方法

采用海綿層阻尼消波法(韓朋等, 2009)進行消波。該方法的原理是, 通過在特定的計算區(qū)域的動量方程中添加阻尼源項來削弱傳入該區(qū)域的波動。阻尼消波的效果與消波段的長度有關(guān), 該長度應至少為入射波波長的一倍。消波區(qū)加入阻尼源項后的動量方程為:

其中, ν為運動黏度系數(shù), ()xφ為阻尼消波系數(shù), 是在阻尼段起點為零的單調(diào)遞增函數(shù)。阻尼消波系數(shù)取線性分布的形式, 如下式所示:

式中: α為阻尼系數(shù)內(nèi)的經(jīng)驗參數(shù),1x為海綿層的起始坐標,2x為海綿層末端坐標;sL為海綿層的長度, 在數(shù)值計算中通常取1—2倍波長。

2 內(nèi)孤立波數(shù)值模擬

2.1 模型建立與網(wǎng)格劃分

王飛(2015)在中國海洋大學物理海洋實驗室, 利用密度分層水槽, 采用重力塌陷法制造不同振幅的內(nèi)孤立波, 通過實驗研究內(nèi)孤立波的生成及流場特性。水槽長15m, 寬0.35m, 高0.7m, 額定水深0.6m。采用 Oster雙缸法制取密度分層流體, 上層流體高度9.5cm, 密度1035kg/m3, 下層流體高度48.5cm, 密度1054kg/m3; 采用 PIV技術(shù)監(jiān)測流場信息。為與實驗結(jié)果進行對比, 本文選取上述流體參數(shù), 以實驗中得到的三個波高作為設定波高進行數(shù)值模擬。

黃文昊等(2013)通過實驗, 研究了四種內(nèi)孤立波理論的適用范圍, 利用內(nèi)孤立波的非線性參數(shù)ε和色散參數(shù)ξ來定量表征各內(nèi)孤立波理論的適用區(qū)間, 結(jié)果表明:ξ＜0.1,ε≤ξ時, KdV理論適用;ξ＜0.1,ξ＜ε＜時, eKdV理論適用。以此為依據(jù), 對三組工況分別選用適用的內(nèi)孤立波理論作為數(shù)值模擬輸入條件。數(shù)值模擬工況如表1所示:

表1 數(shù)值模擬工況Tab.1 Conditions of the numerical simulations

網(wǎng)格劃分: 數(shù)值水槽采用結(jié)構(gòu)化網(wǎng)格進行劃分,x方向網(wǎng)格尺寸取Lx= 2 cm,z方向網(wǎng)格尺寸取Lz= 0 .5cm, 共 87000個網(wǎng)格。消波區(qū)的長度取為1—2倍波長。

源的設置: 兩質(zhì)量源的大小應滿足點源假設, 即其高度Hs應遠小于數(shù)值水槽的水深H, 寬度LΩ應遠小于內(nèi)孤立波的波長λ, 通常質(zhì)量源寬度應小于波長的5% (Linet al, 1999; Hafsiaet al, 2009),故將源Ω1、Ω2分別放置在左邊界處上、下層流體中的一個網(wǎng)格內(nèi)。Lin等(1999)指出質(zhì)量源位置距離分界面1/3層厚處模擬結(jié)果與理論最相符, 并將此位置稱為“standard elevation”。故參考 Lin 等(1999)的結(jié)論, 將兩源距離兩層流體分界面的距離各設為相應層厚度的 1/3, 即Ω1距離兩層流體分界面 3.5cm,Ω2距離分界面 16cm。從后面的模擬結(jié)果來看, 這樣的設置取得了良好的效果。

2.2 FLUENT設置

邊界條件: 水槽左邊界設為對稱邊界; 水表面采用剛蓋假設(方欣華, 2005), 設為固壁邊界; 水底設為固壁邊界; 右邊界設置為壓力出口。

液面追蹤方法: 采用VOF方法(韓朋等, 2009)來追蹤兩層流體分界面處的內(nèi)孤立波波面。流體體積函數(shù)F定義為單元內(nèi)流體所占有的體積與該單元可容納流體體積之比。若單元體被指定相流體占滿, 則該單元F值為 1; 若單元體被另一相流體占滿, 則單元的F值為0; 若單元體的F值在0與1之間, 則表示該單元為包含兩相流體的交界面單元。

UDF設置: 利用 FLUENT中的 DEFINE_SOURCE宏定義上下兩源的源強, 添加到 FLUENT中造波源區(qū)域的源項加載區(qū), 嵌入到造波區(qū)控制方程中; 利用 DEFINE_PROFILE宏在壓力出口邊界給定計算域靜水壓強, 使得出口邊界處水體能夠自由流出, 保證計算域內(nèi)壓力出口附近上下層水深為定值(張齊焰, 2008; 陳鈺等,2009)。

求解設置: 控制方程的對流項采用一階迎風格式,擴散項采用中心差分格式, 壓力速度耦合方式采用PISO算法, 壓力插值選用體積力加權(quán)法, 兩相界面的構(gòu)造方法選用幾何重構(gòu); 初始時水槽中流場速度取為0,時間步長設為0.01s, 每個時間步最大迭代次數(shù)20次。

3 數(shù)值造波結(jié)果及分析

3.1 質(zhì)量源法中內(nèi)孤立波的生成

兩質(zhì)量源的強度隨模擬時間不斷發(fā)生變化, 使得造波源附近的流場速度發(fā)生變化。圖 2為工況 A中兩造波源附近兩個時刻的速度矢量圖, 圖中白線分別為數(shù)值水槽左邊界和上邊界, 黑色矩形分別為質(zhì)量源Ω1和Ω2。從速度矢量圖中可以看出, 隨著源Ω1不斷釋放質(zhì)量、源Ω2不斷吸收質(zhì)量, 源附近上下層流體的速度場均發(fā)生變化, 上層流體產(chǎn)生沿x正向的水平速度, 下層流體產(chǎn)生沿x負向的水平速度, 這與內(nèi)孤立波波致流場中的速度趨勢是相同的。

圖3為 A工況造波初期四個不同時刻質(zhì)量源區(qū)域附近兩相界面圖, 展示了造波開始后左側(cè)邊界附近兩相界面的變化, 其中灰色部分為上層流體, 黑色部分為下層流體。從圖中可以看出, 兩層流體的界面在源Ω1不斷釋放質(zhì)量、源Ω2不斷吸收質(zhì)量的驅(qū)動下,緩慢而平穩(wěn)地發(fā)生下凹, 并最終形成一個完整的下凹形內(nèi)孤立波。

圖2 質(zhì)量源附近速度矢量圖Fig.2 Velocity vector near the mass source

圖3 質(zhì)量源附近兩相界面變化圖Fig.3 Change at the boundary of two phases near mass source

3.2 數(shù)值造波結(jié)果及對比

圖 4為三個工況下內(nèi)孤立波數(shù)值模擬所得波形(CFD)與理論波形(KDV)和實驗波形的對比。從圖中可以看出, 三個工況下數(shù)值模擬波形均光滑對稱, 且與理論及實驗波形吻合良好; 各工況數(shù)值模擬所得波高分別為-3.988cm、-7.320cm和-8.092cm, 與設定波高吻合良好。從對比中可以看出, 本文所述質(zhì)量源方法模擬內(nèi)孤立波所得波形和波高同理論解和實驗結(jié)果是一致的。

圖4 數(shù)值模擬波形與理論及實驗的對比Fig.4 Comparison in the results of numerical, theoretical and experimental simulations

圖 5為波谷斷面處波致水平流速沿垂向分布數(shù)值模擬、理論解及實驗結(jié)果的對比圖。從圖中可以看出, 工況A中數(shù)值模擬與理論解及實驗均吻合良好。工況 B、C中, 上層流體中波致水平流速較大, 致使實驗監(jiān)測數(shù)據(jù)失效, 但在該層數(shù)值模擬與理論解吻合較好; 在下層流體中, 數(shù)值模擬與理論解及實驗均吻合較好。

圖6為內(nèi)孤立波波致水平流速數(shù)值模擬、理論解及實驗的對比圖。從對比中可以看出, 工況A中模擬峰值比理論和實驗結(jié)果略小, 圖像介于理論解和實驗之間, 整體上與二者吻合良好; 工況 B、C中同樣由于上層流體的流速較大, 超出系統(tǒng)的測量范圍, 造成實驗監(jiān)測數(shù)據(jù)失真, 仍可以看出, 上層流體中數(shù)值模擬與理論解吻合良好, 下層流體中三者的吻合度均良好。

4 結(jié)論

采用兩個點源形式的質(zhì)量源, 分別放置于兩層流體中的上下層中, 上層流體中的質(zhì)量源釋放流體,下層流體中的質(zhì)量源吸收流體, 激發(fā)產(chǎn)生下凹型內(nèi)孤立波。推導質(zhì)量源強度表達式, 并利用商業(yè)軟件FLUENT, 基于其 UDF編譯功能, 以內(nèi)孤立波 KdV和 eKdV理論為依據(jù), 從不可壓縮流體的 Navier-Stokes方程出發(fā), 發(fā)展了一種內(nèi)孤立波數(shù)值造波方法,并將數(shù)值模擬結(jié)果與理論及實驗進行對比。結(jié)果表明:

圖5 波谷斷面水平流速沿垂向分布Fig.5 Profiles of horizontal velocity at trough section

圖6 波致水平流速對比圖Fig.6 The horizontal vector of wave-induced velocity

(1) 造波過程中, 在上層流體中的源釋放質(zhì)量、下層流體中的源吸收質(zhì)量的驅(qū)使下, 上下層流體速度呈現(xiàn)出相反的情形, 使得兩層流體分界面逐漸下凹, 最終形成一個不斷向前傳播的下凹形內(nèi)孤立波。

(2) 三個工況下波形、波高及波致水平速度的數(shù)值模擬結(jié)果均與理論及實驗的吻合較好, 驗證了此方法的有效性, 表明本文所述質(zhì)量源法能夠有效模擬內(nèi)孤立波的生成與傳播。

(3) 該方法作為一種純數(shù)值造波方法, 避免了動網(wǎng)格的使用, 計算效率高, 耗時短, 另外不需要已知速度等邊界條件, 可用于研究速度邊界未知的情況,如波流相互作用等。

(4) 該方法同樣可以模擬上凸型內(nèi)孤立波。對于上凸型內(nèi)孤立波, 其波面 η （ t ）始終為正, 則(5)、(6)兩式中的源強表達式 S1（t）始終為負, S2（t）始終為正,表示上層流體中的源吸收質(zhì)量, 下層流體中的源釋放質(zhì)量, 從而驅(qū)動兩層流體界面發(fā)生上凸, 形成上凸型內(nèi)孤立波。

實際海洋中環(huán)境復雜, 由于表面波、流以及海底地形等因素的影響, 內(nèi)孤立波波形通常不如理論波形和數(shù)值模擬波形那樣光滑對稱(Xu et al, 2012),在以后的內(nèi)孤立波數(shù)值模擬中可以考慮波流相互耦合以及復雜海底地形下的內(nèi)孤立波的數(shù)值模擬, 這樣得到的波形將更趨近于實際海洋中的內(nèi)孤立波的波形。

方欣華, 杜 濤, 2005. 海洋內(nèi)波基礎和中國海內(nèi)波. 青島:中國海洋大學出版社, 101—115, 258—281

王 飛, 2015. 內(nèi)孤立波作用下小尺度豎直圓柱體的水動力特性研究. 青島: 中國海洋大學博士學位論文, 21—32

王 偉, 郭海燕, 王 飛等, 2016. 內(nèi)孤立波波致流場數(shù)值模擬研究. 海洋與湖沼, 47(3): 502—508

王 旭, 林忠義, 尤云祥, 2014. 兩層流體中內(nèi)孤立波質(zhì)量源數(shù)值造波方法. 上海交通大學學報, 48(6): 850—855

付東明, 尤云祥, 李 巍, 2009. 兩層流體中內(nèi)孤立波與潛體相互作用數(shù)值模擬. 海洋工程, 27(3): 38—44

劉秀麗, 段夢蘭, 高 攀等, 2015. 基于 OpenFoam 的數(shù)值波浪水槽研究. 復旦學報(自然科學版), 54(3): 373—378,385

張 波, 黃長穆, 1998. 中國南海流花11-1油田的深水開發(fā)技術(shù). 中國海上油氣(工程), 10(3): 36—44

張齊焰, 2008. 非線性波浪沖擊透空結(jié)構(gòu)物數(shù)值模擬. 杭州:浙江大學碩士學位論文, 31—40

李水娟, 余真真, 羅龍洪等, 2012. 三維數(shù)值波浪水槽內(nèi)孤立波特性分析. 水電能源科學, 30(7): 96—99

杜 濤, 吳 巍, 方欣華, 2001. 海洋內(nèi)波的產(chǎn)生與分布. 海洋科學, 25(4): 25—28

陳 鈺, 朱良生, 2009. 基于Fluent的海洋內(nèi)孤立波數(shù)值水槽模擬. 海洋技術(shù), 28(4): 72—75, 100

周玲玲, 丁全林, 蘭慶琳等, 2015. 基于LB方法的質(zhì)量源造波的數(shù)值波浪水槽. 水電能源科學, 33(3): 99—103

柯自明, 尹寶樹, 徐振華等, 2009. 南海文昌海域內(nèi)孤立波特征觀測研究. 海洋與湖沼, 40(3): 269—274

高原雪, 尤云祥, 王 旭等, 2012. 基于MCC理論的內(nèi)孤立波數(shù)值模擬. 海洋工程, 30(4): 29—36

黃文昊, 尤云祥, 王 旭等, 2013. 有限深兩層流體中內(nèi)孤立波造波實驗及其理論模型. 物理學報, 62(8): 084705

韓 朋, 任 冰, 李雪臨等, 2009. 基于VOF方法的不規(guī)則波數(shù)值波浪水槽的阻尼消波研究. 水道港口, 30(1): 9—13

Brorsen M, Larsen J, 1987. Source generation of nonlinear gravity waves with the boundary integral equation method.Coastal Engineering, 11(2): 93—113

Choi W, Camassa R, 1996. Weakly nonlinear internal waves in a two-fluid system. Journal of Fluid Mechanics, 313: 83—103

Choi W, Camassa R, 1999. Fully nonlinear internal waves in a two-fluid system. Journal of Fluid Mechanics, 396: 1—36

Hafsia Z, Hadj M B, Lamloumi H et al, 2009. Internal inlet for wave generation and absorption treatment. Coastal Engineering, 56(9): 951—959

Helfrich K R, Melville W K, 2006. Long nonlinear internal waves. Annual Review of Fluid Mechanics, 38: 395—425

Kim G, Lee M E, 2010. On the internal generation of waves: control volume approach. Ocean Engineering, 38(8—9): 1027—1029

Michallet H, Barthélemy E, 1998. Experimental study of interfacial solitary waves. Journal of Fluid Mechanics,366(1): 159—177

Osborne A R, Burch T L, 1980. Internal solitons in the Andaman Sea. Science, 208(4443): 451—460

Roberts J, 1975. Internal Gravity Waves in the Ocean. New York:Marcel Dekker Inc, 2—13

Xu Z H, Yin B S, Yang H W et al, 2012. Depression and elevation internal solitary waves in a two-layer fluid and their forces on cylindrical piles. Chinese Journal of Oceanology and Limnology, 30(4): 703—712