張周昊，張洪生，王宇鑫，馬婷婷

(上海海事大學(xué) 海洋科學(xué)與工程學(xué)院，上海 201306)

海洋內(nèi)波是最大振幅位于海洋內(nèi)部的波動(dòng)，海水穩(wěn)定層化和外力擾動(dòng)是內(nèi)波發(fā)生的必要條件。內(nèi)波作為能量和動(dòng)量垂向運(yùn)輸?shù)闹匾d體，有利于海洋生態(tài)的保護(hù)，但會(huì)對(duì)海洋工程和水下艦只的安全造成很大的威脅，故深入研究海洋內(nèi)波的傳播特性具有重要意義。研究海洋內(nèi)波主要有理論研究、物理模型試驗(yàn)、衛(wèi)星遙感和數(shù)值模擬方法。其中數(shù)值模擬方法不僅操作相對(duì)簡便而且經(jīng)濟(jì)高效，還可以與理論研究等方法互為對(duì)照。隨著計(jì)算機(jī)技術(shù)的發(fā)展，數(shù)值模擬方法的應(yīng)用越來越廣泛[1-2]。

Zhang 等[3]基于Euler 方程，采用改進(jìn)的SIMPLE 算法，模擬了密度為連續(xù)分層流體中內(nèi)波的傳播和變形，并將等水深水域的數(shù)值解和理論解進(jìn)行了對(duì)比，定性分析了存在潛堤的變水深水域中不同時(shí)刻的密度分布，研究了變水深水域中內(nèi)波傳播的特性。李景遠(yuǎn)等[4]基于OpenFOAM 平臺(tái)，引入密度輸運(yùn)方程，基于雷諾時(shí)均的N-S 方程，模擬了密度連續(xù)分層流體中內(nèi)孤立波的傳播。

雖然密度連續(xù)分層更貼近海洋的真實(shí)情況，但把實(shí)際海洋概化為兩層密度不同的流體不失為一種有效的簡化手段。高原雪等[5]利用FLUENT 軟件，以Euler 方程為控制方程，模擬了密度為兩層的水域中內(nèi)孤立波的生成與傳播。陳鈺等[6]基于FLUENT 軟件，通過給定水位和入口速度的造波方法，數(shù)值模擬了弱非線性內(nèi)孤立波的傳播。Terez 等[7]基于二維不可壓的Boussinesq 近似的N-S 方程，用重力塌陷的造波方法數(shù)值模擬了大波幅內(nèi)孤立波，通過追蹤分布在混合區(qū)和密躍層的拉格朗日粒子觀察運(yùn)動(dòng)。韓鵬等[8]對(duì)長非線性內(nèi)孤立波在斜坡地形上的傳播和演化進(jìn)行了理論和數(shù)值模擬研究。Baie 等[9]應(yīng)用k-ε模型，實(shí)現(xiàn)了潛艇生成內(nèi)波的數(shù)值模擬。在利用FLUENT 軟件進(jìn)行計(jì)算時(shí)，不僅在內(nèi)波的數(shù)值模擬中經(jīng)常使用紊流模型，而且在表面波的數(shù)值模擬中也經(jīng)常使用紊流模型[10]。

海洋表面邊界條件可分為剛蓋假定和自由表面假定兩種邊界條件[11-12]，雖然將剛蓋假定視為海面近似邊界條件具有一定的合理性，但實(shí)際海洋中海面并非如同剛性表面一般靜止不動(dòng)，因此模擬在自由表面邊界條件下內(nèi)波的傳播特性具有重要的研究意義。本文利用FLUENT 軟件，對(duì)兩層密度或深度不同的流體，上邊界條件分別作剛蓋假定和自由表面假定，對(duì)微幅內(nèi)波的傳播進(jìn)行數(shù)值模擬，研究在上表面兩種假定下，密度場和速度場的異同及其產(chǎn)生原因。

1 兩層流體中小振幅內(nèi)波傳播的理論解

小振幅內(nèi)波問題可歸結(jié)為二維的線性問題。假定海水中沒有其他流動(dòng)，x軸取在靜止的兩層流體交界面上，z軸取向上為正，如圖1 所示，則不可壓縮流體的運(yùn)動(dòng)方程和連續(xù)性方程可表示為[11]：

圖1 穩(wěn)定層化兩層流體中內(nèi)波傳播示意Fig.1 Sketch of propagation of internal waves in density stratified two-layer fluid

式中：u和w分別為流體質(zhì)點(diǎn)水平和垂向速度；t為時(shí)間；p為壓強(qiáng)；ρ0=ρ0(z)為靜止時(shí)分層流體穩(wěn)定的密度分布，假定其只隨深度變化；ρ′為流體的擾動(dòng)密度；g為重力加速度。

假定上、下層流體交界面處內(nèi)波頻率為 ω、水平波數(shù)為k，文獻(xiàn)[11]推導(dǎo)了剛蓋假定下內(nèi)波的頻散關(guān)系式和速度場理論解；文獻(xiàn)[12-14]推導(dǎo)了自由表面邊界條件下內(nèi)波的頻散關(guān)系式和速度場理論解。剛蓋假定下的頻散關(guān)系式為：

式中：ρ1、ρ2分別為上、下層流體密度；b1、b2分別為上、下層流體深度。

自由表面假定下的頻散關(guān)系式為：

式中：k0=ω2/g；γ=ρ1/ρ2；t1=th(kb1)，t2=th(kb2)。

2 兩層流體中內(nèi)波傳播的數(shù)值模型

2.1 模型的建立及邊界條件的設(shè)置

數(shù)值水槽如圖2 所示，原點(diǎn)取在模型左側(cè)的兩層流體交界處。數(shù)值水槽中內(nèi)波傳播的有效區(qū)域?yàn)镠JEF，其長度取各工況下7 倍波長L，消波區(qū)域?yàn)镕EKI，消波區(qū)長度取為1.5L，設(shè)定總水深為100 m。給定振幅A為1.0 m，周期T為400 s。

圖2 數(shù)值水槽示意Fig.2 Sketch of numerical wave flume

本文采用FLUENT 軟件中的k-ε模型進(jìn)行數(shù)值模擬，并利用VOF(Volume Of Fluid)方法追蹤自由液面位置，模型的控制方程和VOF 法的介紹可參見文獻(xiàn)[6,15]，方程的離散采用了有限體積法。通過改變HI處的邊界條件來實(shí)現(xiàn)兩種邊界條件下內(nèi)波傳播的數(shù)值模擬。在模擬剛蓋假定下內(nèi)波的傳播時(shí)，HI處為對(duì)稱（symmetry）邊界條件；模擬自由表面假定下內(nèi)波的傳播時(shí)，增加高度為50 m 的空氣層，空氣層上表面HI處為壓力進(jìn)口（Pressure inlet）邊界條件。設(shè)HJ、JK和KI均為固壁邊界條件。

模型采用結(jié)構(gòu)化網(wǎng)格，傳播區(qū)x方向網(wǎng)格大小為各工況波長L的2%；在兩層流體交界面以及上層流體的上表面及其鄰近處，z方向采用加密網(wǎng)格，其余位置z方向網(wǎng)格均取1 m。在消波區(qū)域的x方向采用1.05 倍網(wǎng)格大小遞增的漸變網(wǎng)格。控制流場全局庫朗數(shù)（Courant number）為0.02，讓FLUENT 自動(dòng)選擇合適的時(shí)間步長，總模擬時(shí)間為6 000 s。

2.2 數(shù)值造波與消波

采用平板拍擊的方法來實(shí)現(xiàn)造波[16]。造波板上下拍動(dòng)，帶動(dòng)流體運(yùn)動(dòng)以達(dá)到造波的目的，其平衡位置在兩層密度不同的流體交界面處。根據(jù)體積守恒定律，造波板上下拍擊時(shí)所排開的流體體積，與波面從造波板邊緣沿內(nèi)波傳播方向至無窮遠(yuǎn)處的積分相等，即：

式中：D為造波板長；z(t)為某時(shí)刻造波板偏離平衡位置的距離；η(x,t)為兩層流體交界面處波面位移函數(shù)。

由于內(nèi)波沿x軸正向傳播，因此：

式中：c為內(nèi)波相速度。將式（5）兩側(cè)對(duì)時(shí)間求導(dǎo)，再結(jié)合式（6），設(shè)η(x,t)=Acos(kx-ωt)，推導(dǎo)可得：

式中：um(t)為造波板垂向運(yùn)動(dòng)速度。由式（7）可見運(yùn)動(dòng)函數(shù)由需要生成的內(nèi)波波速、振幅、波數(shù)、頻率和造波板長度確定。

采用FLUENT 用戶自定義函數(shù)（UDF）中的剛體運(yùn)動(dòng)DEFINE_CG_MOTION 宏，結(jié)合推導(dǎo)的運(yùn)動(dòng)方程，控制造波板上下平移。

在水槽的末端，即在FEKI區(qū)域設(shè)置消波層進(jìn)行消波。消波的方法是在消波區(qū)內(nèi)，在每個(gè)迭代步內(nèi)速度u和w都要乘以消波系數(shù)μ(x)。μ(x)有線性和指數(shù)等多種形式?，F(xiàn)選取線性形式：

式中：a為 常數(shù)，本文取為0.02；x0為消波區(qū)起始位置水平坐標(biāo)；x為消波區(qū)內(nèi)任一點(diǎn)的水平坐標(biāo)；B為消波區(qū)長度，常取為1～2 倍波長L。

2.3 敏感性分析

為了驗(yàn)證計(jì)算結(jié)果的正確性，需對(duì)采用不同Δz大小網(wǎng)格的計(jì)算結(jié)果進(jìn)行敏感性分析。采用ρ1=1 020 kg/m3、ρ2=1 028 kg/m3、b1=10 m、b2=90 m 的剛蓋假定數(shù)值水槽作為對(duì)比分析算例，該算例比其他算例（下文論及）的數(shù)值解與相應(yīng)理論解之間的差別要大。Δz分別取0.5 和1.0 m。圖3 是x=1 200m、z=0 m 處流體質(zhì)點(diǎn)垂向速度隨時(shí)間的變化曲線。由圖3 可見，兩種網(wǎng)格的計(jì)算結(jié)果差別很小，與理論解的變化趨勢(shì)一致。鑒于采用Δz=1.0 m時(shí)能減少計(jì)算時(shí)間，因此本文采用Δz=1.0 m的網(wǎng)格進(jìn)行計(jì)算。

圖3 (1 200,0)處垂向速度隨時(shí)間的變化過程線Fig.3 Time series of vertical velocity at point (1 200,0)

3 不同情況下速度場數(shù)值解與相應(yīng)理論解的比較

3.1 上下層流體不同密度差的影響

工況設(shè)置如表1所示。假定ρ1、ρ2、b1、b2四個(gè)參數(shù)及周期T，用兩種假定下的色散關(guān)系求得各工況下波數(shù)k（同種工況兩種假定下波數(shù)k均取3 位有效數(shù)字后近似相等）。由表1 可知，在上下層水深及波周期不變的情況下，上下層流體的密度差越大，波數(shù)k就越小，內(nèi)波傳播速度越快。

表1 上下層不同密度工況設(shè)置Tab.1 Cases with different densities for upper and lower layers

在x方向傳播區(qū)長度的一半處，設(shè)置上下兩個(gè)測(cè)點(diǎn)。下測(cè)點(diǎn)取在兩層流體交界面處（z=0 m），監(jiān)測(cè)該處流體質(zhì)點(diǎn)的垂向速度w；上測(cè)點(diǎn)取在上層流體的上表面（z=50m），監(jiān)測(cè)該處流體質(zhì)點(diǎn)的垂向速度w和水平速度u。計(jì)算的水氣交界面處水平速度結(jié)果見圖4～5。由于篇幅所限，只給出Case 1 及Case 4 的部分結(jié)果。

圖4 Case 1(750,50)位置水平速度歷時(shí)曲線Fig.4 Time series of horizontal velocity at point (750,50) for Case 1

圖5 Case 4(1 750,50)位置水平速度歷時(shí)曲線Fig.5 Time series of horizontal velocity at point (1 750,50) for Case 4

由圖4～5 可知，各工況中的速度在兩種邊界條件下，數(shù)值解與各自的理論解吻合良好；剛蓋與自由表面兩種邊界條件下的數(shù)值解基本無差別。對(duì)比圖4 和5 可知，如其他條件不變，上下層流體之間密度差越大，上層流體表面的水平速度振幅就越大，水平運(yùn)動(dòng)越劇烈。鑒于上層流體表面的垂向速度，在剛蓋條件下為0，在自由表面條件下盡管不等于0，但幅值很小，因此沒有給出其計(jì)算結(jié)果圖。

計(jì)算了波峰和波谷時(shí)垂向速度和水平速度沿水深的分布，波峰和波谷均取在波型充分穩(wěn)定的時(shí)間段。限于篇幅，文中僅給出了水平速度沿水深分布的計(jì)算結(jié)果，見圖6～7。在兩種邊界條件下，數(shù)值解與理論解基本一致。在波峰時(shí)刻，下層流體的數(shù)值解和理論解吻合得更好；在波谷時(shí)刻，上層流體的數(shù)值解和理論解吻合得更好。

圖6 自由表面假定下波峰時(shí)水平速度沿水深的分布Fig.6 Distribution of horizontal velocity along water depth at the moment for wave crest under the assumption of free surface

圖7 剛蓋假定下波谷時(shí)水平速度沿水深的分布Fig.7 Distribution of horizontal velocity along water depth at the moment for wave trough under the assumption of rigid lid

3.2 上下層流體不同水深比的影響

工況設(shè)置如表2 所示。控制上下層流體的密度不變，改變深度比。由于篇幅所限，只給出工況較為極端的Case 5 和Case 8 的計(jì)算結(jié)果，不同位置處流體質(zhì)點(diǎn)的速度歷時(shí)曲線見圖8～12。

圖8 Case 5(1 200,0)位置垂向速度歷時(shí)曲線Fig.8 Time series of vertical velocity at point (1 200,0) for Case 5

表2 上下層不同水深工況設(shè)置Tab.2 Cases with different water depths of upper and lower layers

由表2 可以明顯看出，當(dāng)上下層流體深度b1、b2互換時(shí)，波數(shù)k近似相等，這體現(xiàn)了色散方程的“對(duì)偶性”[17]。

圖9 Case 5(1 200,10)位置垂向速度數(shù)值解歷時(shí)曲線Fig.9 Time series of numerical vertical velocity at point (1 200,10) for Case 5

圖10 Case 5(1 200,10)位置水平速度數(shù)值解歷時(shí)曲線Fig.10 Time series of numerical horizontal velocity at point (1 200,10) for Case 5

圖11 Case 8(1 200,0)位置垂向速度歷時(shí)曲線Fig.11 Time series of vertical velocity at point (1 200,0) for Case 8

由圖8 和11 可見，兩層流體交界面處的垂向速度，無論在哪種邊界條件下，數(shù)值解與理論解相比均已有較大區(qū)別，體現(xiàn)出一定的非線性影響；當(dāng)b1小于b2時(shí)，數(shù)值解的波峰處幅值小于波谷處幅值；當(dāng)b1大于b2時(shí)，數(shù)值解的波峰處幅值大于波谷處幅值。當(dāng)b1小于b2時(shí)，剛蓋假定下所計(jì)算的垂向速度幅值大于自由表面假定下所計(jì)算的垂向速度幅值；當(dāng)b1大于b2時(shí)則相反。

圖12 Case 8(1 200,90)位置水平速度數(shù)值解歷時(shí)曲線Fig.12 Time series of numerical horizontal velocity at point (1 200,90) for Case 8

由圖9 可見，當(dāng)b1很小時(shí)，自由表面假定下的數(shù)值解有明顯的波動(dòng)。而當(dāng)b1大于b2時(shí)，波動(dòng)幾乎為0，因此文中沒有給出其計(jì)算結(jié)果圖。圖10 和12 是上層流體表面的水平速度隨時(shí)間變化曲線?？梢?，當(dāng)b1小于b2時(shí)，數(shù)值解的波峰尖陡、波谷平坦，呈現(xiàn)出非線性影響；而當(dāng)b1大于b2時(shí)，則非線性影響微弱。在入射條件相同的情況下，b1越小上層流體的表面速度越大。

鑒于波谷時(shí)刻與波峰時(shí)刻的速度沿水深分布圖相類似及兩種假定下數(shù)值計(jì)算結(jié)果類似，文中僅給出了波峰時(shí)的垂向和水平速度沿水深分布，見圖13～14，R 代表剛蓋，F(xiàn) 代表自由表面。由圖13 可見數(shù)值解垂向速度的極大值不再位于兩層流體交界面處。b1小于b2時(shí)，垂向速度極大值的位置向下偏移，導(dǎo)致監(jiān)測(cè)點(diǎn)記錄的數(shù)據(jù)更小；b1大于b2時(shí)，極大值的位置向上偏移，監(jiān)測(cè)點(diǎn)記錄到的數(shù)據(jù)也更小。在圖8 和11 中也體現(xiàn)出波峰時(shí)數(shù)值解的幅值相對(duì)于理論解有所偏小。圖14 是波峰時(shí)刻水平速度沿水深的分布。由圖14 可見，密躍層越靠近上層表面，上層流體內(nèi)數(shù)值解和理論解差別越大，下層流體內(nèi)差別越小，反之亦然。同種工況內(nèi)，兩種邊界條件下的流體質(zhì)點(diǎn)水平速度沿水深的分布，沒有明顯的區(qū)別。

圖13 波峰時(shí)垂向速度沿水深的分布圖Fig.13 Distribution of vertical velocity along water depth at the moment for wave crest

圖14 波峰時(shí)水平速度沿水深的分布圖Fig.14 Distribution of horizontal velocity along water depth at the moment for wave crest

各工況在波峰波谷時(shí)刻，水平速度與垂向速度的數(shù)值解和理論解之間的一致程度可用指標(biāo)d來衡量[18-19]：

式中：x(j) 為理論解；為 理論解均值；y(j)為數(shù)值解。d=0時(shí)表示數(shù)值解與理論解完全不一致，而d=1表示數(shù)值解與理論解完美匹配。

由表3 及圖4～7 可見，Case 1～Case 4 盡管下層流體密度不同，但上下兩層深度相等，密躍層剛好位于水深的1/2 處，此時(shí)總體說來，數(shù)值解與理論解相差較小。Case 5～Case 8 上下兩層流體密度固定，但上下兩層深度差別較大，由表3 及圖8～14 可見，此時(shí)數(shù)值解與理論解的差別較大。特別是當(dāng)上下兩層水深差值較大時(shí)，數(shù)值解與理論解的差別更大，其中又以密躍層接近靜水面時(shí)表現(xiàn)得更為明顯，這體現(xiàn)了非線性的影響。這是因?yàn)閿?shù)值解是非線性解，而理論解是線性解。在上層水深很小時(shí)，在自由表面假定下水氣交界面處出現(xiàn)了較為明顯的垂向速度。鑒于在實(shí)際海洋中上層水深遠(yuǎn)小于下層水深，由此可以推斷，對(duì)于實(shí)際的海洋分層流體，采用更為真實(shí)的自由表面假定應(yīng)更為合理。

表3 各工況水平和垂向速度的數(shù)值解和理論解一致性分析Tab.3 Consistency analysis between the numerical and analytical results for the horizontal and vertical velocities of different cases

4 結(jié)語

本文基于k-ε模型和VOF 方法，以FLUENT 軟件為計(jì)算平臺(tái)，采用平板造波的方法，實(shí)現(xiàn)了剛蓋假定和自由表面假定下小振幅內(nèi)波傳播的數(shù)值模擬，通過比較z方向不同網(wǎng)格大小的計(jì)算結(jié)果，檢驗(yàn)了計(jì)算結(jié)果的合理性。

設(shè)置上下兩層流體不同的密度差和深度，數(shù)值模擬了剛蓋假定和自由表面假定下多種情況內(nèi)波的傳播；通過與各自理論解的比較，闡述了數(shù)值解與理論解之間的異同；通過對(duì)兩種假定下計(jì)算結(jié)果的比較，分析了兩種假定下數(shù)值計(jì)算結(jié)果之間的差別及其產(chǎn)生的原因。在上下層水深相同的前提下，盡管上下層流體密度差不同，但此時(shí)總體上數(shù)值解與理論解相吻合；盡管在自由表面假定下，靜水面處存在垂向速度，但其數(shù)值很小。在上下兩層流體密度差固定的前提下，上下兩層流體深度差值的改變會(huì)較為明顯地影響數(shù)值計(jì)算結(jié)果，以密躍層接近海洋上表面時(shí)表現(xiàn)得更為明顯，此時(shí)數(shù)值解與理論解差別較大，這體現(xiàn)了非線性的影響；在自由表面假定下出現(xiàn)了較為明顯的垂向速度，這體現(xiàn)了海洋上表面邊界條件的影響；在給定入射條件下，上層流體水深越小，上表面處速度越大。

鑒于在實(shí)際海洋中，一般密躍層以上水深遠(yuǎn)小于下層水深。當(dāng)計(jì)算運(yùn)動(dòng)幅值更大的內(nèi)孤立波時(shí)，可以設(shè)想此時(shí)邊界條件的影響會(huì)更加明顯，這便是本研究意義所在，也是下一步將要研究的內(nèi)容。本文工作也為內(nèi)波與表面波之間相互作用的研究提供了參考。