密度連續(xù)變化水體內(nèi)孤立波數(shù)值模擬研究

李景遠(yuǎn)，張慶河，陳同慶

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

內(nèi)孤立波是發(fā)生在密度穩(wěn)定層化海水內(nèi)部的一種波動(dòng)，是實(shí)際海洋環(huán)境中常見(jiàn)的海洋動(dòng)力學(xué)現(xiàn)象．內(nèi)孤立波在傳播過(guò)程中所攜帶的巨大能量可誘發(fā)水體強(qiáng)烈輻聚輻散和突發(fā)性強(qiáng)流[1-2]，對(duì)海域的生態(tài)環(huán)境、油氣資源開(kāi)發(fā)、水下潛航器安全等具有重要影響[3-5]，其生成、傳播、演變機(jī)制研究一直受到國(guó)內(nèi)外學(xué)者廣泛關(guān)注．

目前關(guān)于內(nèi)孤立波的研究除了一部分現(xiàn)場(chǎng)觀測(cè)外[6]，主要集中于實(shí)驗(yàn)室物理模擬和數(shù)值模擬工作．谷夢(mèng)夢(mèng)等[7]對(duì)過(guò)山脊內(nèi)孤立波進(jìn)行了實(shí)驗(yàn)研究，給出了過(guò)山脊內(nèi)孤立波的演變規(guī)律．黃文昊等[8]基于內(nèi)孤立波理論解發(fā)展了一種雙推板內(nèi)孤立波實(shí)驗(yàn)室造波方法，根據(jù)實(shí)驗(yàn)結(jié)果分析了KdV、eKdV、MCC和mKdV 理論的適用條件．黃鵬起等[9]在二維內(nèi)波水槽中研究了內(nèi)孤立波在斜坡上的演變過(guò)程，分析了內(nèi)孤立波破碎時(shí)的湍動(dòng)能、湍流耗散率等變化規(guī)律．Kao 等[10]配置了具有有限厚度密度躍層的內(nèi)孤立波水槽，結(jié)果表明，相較于兩層模型，密度躍層的存在會(huì)削減內(nèi)孤立波的傳播速度．Segur 等[11]的實(shí)驗(yàn)也有類似結(jié)論．杜輝等[12]、武軍林等[13]、黃文昊等[14]、王玲玲等[15]利用重力式分層流水槽研究了內(nèi)孤立波對(duì)水中結(jié)構(gòu)物的作用．這些實(shí)驗(yàn)研究對(duì)認(rèn)識(shí)內(nèi)孤立波生成和演變以及內(nèi)孤立波對(duì)結(jié)構(gòu)物的作用具有重要意義，但實(shí)驗(yàn)研究受到實(shí)驗(yàn)室尺度等的制約，數(shù)值模型可以彌補(bǔ)這方面的不足，已成為內(nèi)孤立波研究的重要途徑．付東明等[16]、黃文昊等[17]基于兩層流體模型，結(jié)合內(nèi)孤立波理論，通過(guò)雙推板、速度入口、質(zhì)量源造波等方法，建立了內(nèi)波數(shù)值水槽，并獲得了與理論解吻合良好的模擬結(jié)果．關(guān)暉等[18]基于Navier-Stokes 方程，使用開(kāi)源軟件Gerris 模擬建立了密度強(qiáng)分層內(nèi)孤立波數(shù)值水槽，并研究了有航速潛艇遭遇內(nèi)孤立波時(shí)的流場(chǎng)形態(tài)和受力載荷特性.徐鑫哲[19]基于勢(shì)流理論，利用Fluent 模擬平板拍擊造波過(guò)程建立了密度強(qiáng)分層二維數(shù)值內(nèi)波水槽．王飛[20]基于Fluent 建立了密度強(qiáng)分層的三維內(nèi)波數(shù)值水槽，研究了內(nèi)孤立波波致流場(chǎng)及其對(duì)豎直圓柱體的荷載特性．

總體來(lái)看，現(xiàn)有大多數(shù)數(shù)值模型研究采用忽略密度躍層厚度的兩層流體模型來(lái)研究?jī)?nèi)孤立波．亞洲海洋國(guó)際聲學(xué)實(shí)驗(yàn)(ASIAEX)的定點(diǎn)溫度鏈現(xiàn)場(chǎng)觀測(cè)數(shù)據(jù)表明，我國(guó)南海沿水深方向存在一定厚度的密度躍層，這是南海內(nèi)孤立波頻發(fā)的先決條件[21-22]．從海洋的實(shí)際物理環(huán)境看，兩層流體的密度強(qiáng)分層模型不能完全表征實(shí)際海洋環(huán)境下內(nèi)孤立波的真實(shí)特性.

為此，為了克服強(qiáng)分層模型對(duì)于內(nèi)孤立波數(shù)值模擬的不足，本文通過(guò)求解Navier-Stokes 方程和密度輸運(yùn)方程，建立密度連續(xù)變化內(nèi)孤立波數(shù)值模型，并對(duì)其進(jìn)行合理性驗(yàn)證．本文的內(nèi)容安排如下：第1 節(jié)對(duì)本文模型控制方程和離散過(guò)程進(jìn)行描述；第2 節(jié)給出密度連續(xù)變化內(nèi)孤立波數(shù)值水池的實(shí)現(xiàn)過(guò)程；第3節(jié)對(duì)數(shù)值模擬網(wǎng)格敏感性進(jìn)行分析，推薦合適的網(wǎng)格尺度；第4 節(jié)對(duì)模型的合理性進(jìn)行驗(yàn)證，分析和研究密度連續(xù)變化流體中內(nèi)孤立波的傳播、演變特性，并與傳統(tǒng)兩層流體模型進(jìn)行比較．

1 控制方程與數(shù)值算法

1.1 控制方程

本文基于開(kāi)源軟件OpenFOAM，通過(guò)添加密度輸運(yùn)方程，創(chuàng)建了密度連續(xù)變化的不可壓縮流體新求解器iswFoam．由于實(shí)際海洋環(huán)境中穩(wěn)定密度躍層的密度變化較小，故引入Boussinesq 近似是合適的．本文模型流體的控制方程為考慮Boussinesq 近似的雷諾時(shí)均N-S 方程(RANS)[23]，即

1.2 控制方程離散及數(shù)值求解

采用有限體積法對(duì)第1.1 節(jié)中的控制方程進(jìn)行求解，變量均定義在單元中心．在(t ,t+Δ t)時(shí)段內(nèi)，對(duì)控制體單元P 做時(shí)間空間積分，將控制方程式(1)和(6)微分方程轉(zhuǎn)化為如下積分方程：

有限體積法離散過(guò)程中應(yīng)用高斯定理將體積分轉(zhuǎn)化為面積分，再將所得面積分形式轉(zhuǎn)化為單元體各面累加求和的形式．下面以積分形式的密度擴(kuò)散項(xiàng)(式(10))為例說(shuō)明這一過(guò)程，圖1 為離散過(guò)程中所需參量的示意．

圖1 三維笛卡爾坐標(biāo)系下密度擴(kuò)散項(xiàng)的離散示意Fig.1 Schematic of the density diffusion term in a discrete elements in a three-dimensional Cartesian coordinate system

按照上述離散過(guò)程，可獲得動(dòng)量方程式(8)和密度對(duì)流擴(kuò)散方程式(9)的離散形式為

動(dòng)量方程對(duì)流項(xiàng)中非線性項(xiàng)采取線性化處理，即對(duì)流項(xiàng)通量采用當(dāng)前時(shí)刻已知的通量．為了防止同位網(wǎng)格系統(tǒng)下求解出現(xiàn)壓力棋盤效應(yīng)，面上的速度計(jì)算采用Rhie-Chow 差值[25]．

離散后的速度壓力耦合方程采用pimple 算法求解，pimple 的基本思想是將每個(gè)時(shí)間步長(zhǎng)內(nèi)用simple穩(wěn)態(tài)算法求解(也就是將每個(gè)時(shí)間步內(nèi)看成穩(wěn)態(tài)流動(dòng))，時(shí)間步長(zhǎng)的步進(jìn)用piso 算法來(lái)完成[26]．每個(gè)時(shí)間步長(zhǎng)內(nèi)，速度壓力求解完成后，將速度代入離散后的密度對(duì)流擴(kuò)散方程式(12)便可求出下一時(shí)刻的密度場(chǎng)ρn1+，更新速度場(chǎng)用于下一時(shí)刻求解．

2 內(nèi)孤立波的生成

目前，描述內(nèi)孤立波的理論主要有KdV、eKdV、mKdV、MCC 理論等[27]，這些理論模型的解析解均基于密度強(qiáng)分層的兩層流體模型給出．對(duì)于更接近實(shí)際的密度連續(xù)變化情況下的海洋內(nèi)孤立波，目前尚未見(jiàn)到解析解，因此內(nèi)孤立波的生成仍按兩層流體解析解給定初值．本文選用可模擬大振幅內(nèi)孤立波的eKdV 方程理論解作為初值．

2.1 eKdV方程及其穩(wěn)態(tài)解

在圖2 所示的坐標(biāo)系中，設(shè)未受擾動(dòng)的兩層流體界面位于Oxy 平面，坐標(biāo)原點(diǎn)在兩層界面處，內(nèi)孤立波沿Ox 軸正向傳播．設(shè)各層均為不可壓縮流體，上層流體厚度和密度分別為 h1和 ρ1，下層流體厚度和密度分別為 h2和 ρ2．ζ ( x, t)表示內(nèi)孤立波的波面函數(shù)，則eKdV 理論的內(nèi)孤立波控制方程及其穩(wěn)態(tài)解[8]如下：

式中：1u 和 u2分別為上、下層流體水質(zhì)點(diǎn)平均水平速度；a 為內(nèi)孤立波振幅；λeKdV為內(nèi)孤立波的特征波長(zhǎng)；ceKdV為內(nèi)孤立波理論解的相速度．

當(dāng)c3＝0 時(shí)，eKdV 方程即可退化為KdV 方程，eKdV 理論解退化為KdV 理論解．

圖2 兩層流體內(nèi)孤立波示意Fig.2 Schematic of ISWs in a two-fluid system

2.2 內(nèi)孤立波的生成

利用第2.1 節(jié)中內(nèi)孤立波波面理論解ζ 以及上、下層平均水平速度1u 和 u2對(duì)計(jì)算域進(jìn)行初始化.Yang 等[28]在中國(guó)南海密度分布的研究發(fā)現(xiàn)，冬季南海的密度垂向分布呈現(xiàn)雙曲正切型式，大多數(shù)學(xué)者在進(jìn)行南海潮流、內(nèi)波等研究時(shí)，亦主要采用雙曲正切型式，如Chen 等[29]和Lamb[30]．基于南海的實(shí)際情況，本文密度的初始垂向分布采用Aghsaee 等[31]提出的雙曲正切形式，即

式中：ρ ( z )為計(jì)算域內(nèi)連續(xù)變化的密度場(chǎng)；ρ1為計(jì)算域頂部的密度值；2ρ 為計(jì)算域底部的密度值；zpyc為密度躍層中心位置；dpyc為密度躍層厚度的1/2；z為垂向坐標(biāo)，當(dāng)dpyc趨近于0 時(shí)，密度躍層厚度趨近于0，此時(shí)密度連續(xù)變化模型可轉(zhuǎn)化為兩層流體模型.

3 內(nèi)孤立波模擬網(wǎng)格敏感性分析

在對(duì)內(nèi)孤立波進(jìn)行模擬時(shí)，需要選擇合理的網(wǎng)格尺度和時(shí)間步長(zhǎng)，以在獲得所需模擬精度的條件下提高計(jì)算效率．由于內(nèi)孤立波理論波長(zhǎng)為無(wú)限長(zhǎng)，因此需要對(duì)內(nèi)孤立波波長(zhǎng)進(jìn)行合理簡(jiǎn)化，確定適合分析的有效波長(zhǎng)．下面通過(guò)引入有效勢(shì)能(available potential energy，APE)概念，確定合理的有效波長(zhǎng)，并在此基礎(chǔ)上進(jìn)一步討論內(nèi)孤立波模擬的網(wǎng)格尺寸．

3.1 內(nèi)孤立波“有效波長(zhǎng)”的確定

考慮一個(gè)二維區(qū)域( x, z)∈[－l , l] × [－H,0]，內(nèi)孤立波有效勢(shì)能[32-33]可定義為

設(shè)波長(zhǎng)λ所對(duì)應(yīng)的截取有效勢(shì)能為APEλ，利用總有效勢(shì)能對(duì)其進(jìn)行無(wú)量綱化得APEnon＝APEλ/APEt，無(wú)量綱截取有效勢(shì)能APEnon亦可理解為截取勢(shì)能占總勢(shì)能的百分比．

利用Michallet 等[34]提出的內(nèi)孤立波特征波長(zhǎng)L對(duì)波長(zhǎng)λ進(jìn)行無(wú)量綱化，得

特征波長(zhǎng)L 的計(jì)算式為

基于 eKdV 理論，在上、下層水深比分別為h1/h2＝1/9、2/8 和3/7 的條件下，設(shè)置系列波高，通過(guò)計(jì)算一系列波長(zhǎng)λ所對(duì)應(yīng)的截取有效勢(shì)能APEλ，獲得APEnon隨λnon的變化規(guī)律，如圖3 所示．由圖3 可知，因?yàn)閮?nèi)孤立波能量主要集中于波峰附近區(qū)域，APEnon隨λnon增大而迅速增大趨近于1．在相同水深比情況下，較大波高對(duì)應(yīng)的APEnon隨λnon更快趨近于1，主要原因?yàn)榧性趦?nèi)孤立波中間區(qū)域的有效勢(shì)能隨波高的增加所占總勢(shì)能比重增大．總體來(lái)看，當(dāng)APEnon＞99.8%時(shí)，所截取的有效勢(shì)能隨著波長(zhǎng)的增加幾乎不產(chǎn)生變化，因此，取APEnon＝99.8%時(shí)所對(duì)應(yīng)的波長(zhǎng)作為有效波長(zhǎng)λA是合理的，能充分反映內(nèi)孤立波的特性．下面的有關(guān)討論主要結(jié)合λA進(jìn)行．

圖3 無(wú)量綱有效勢(shì)能隨波長(zhǎng)變化關(guān)系Fig.3 Relationship of dimensionless available potential energy with wavelength

3.2 空間步長(zhǎng)測(cè)試

以h1/h2＝2/8、波高H＝16.5 cm 為例分析計(jì)算空間步長(zhǎng)對(duì)模擬結(jié)果的影響，測(cè)點(diǎn)設(shè)置在距離造波中心8λA處．表1 和表2 列出了內(nèi)孤立波在不同網(wǎng)格系統(tǒng)下傳播至8λA位置測(cè)點(diǎn)的波高結(jié)果及其與理論解的對(duì)比．

表1 垂直方向波高范圍內(nèi)空間步長(zhǎng)測(cè)試取值Tab.1 Values of the space step test in the vertical wave height range

表2 沿波浪傳播方向空間步長(zhǎng)測(cè)試取值Tab.2 Values of the space step test along the wave propagation direction

由表1 可知，隨著垂向步長(zhǎng)的減小，測(cè)點(diǎn)位置處的波高值逐漸增加，誤差逐漸減小，且在垂向步長(zhǎng)減小至H/15 時(shí)趨于穩(wěn)定．由表2 可知，水平向步長(zhǎng)Aλ/20 誤差率最大，水平向步長(zhǎng)達(dá)到Aλ/80 后，誤差隨著步長(zhǎng)的減小基本不再減小，同時(shí)波高基本穩(wěn)定．綜合考慮波形分辨率和計(jì)算效率，本文模型水平向步長(zhǎng)取Aλ/80、垂向步長(zhǎng)取H/15 是合適的．

4 數(shù)值模擬結(jié)果與分析

4.1 算例設(shè)置

本文結(jié)合黃文昊等密度強(qiáng)分層水槽實(shí)驗(yàn)，設(shè)置一系列數(shù)值算例(表3)．內(nèi)孤立波數(shù)值水槽長(zhǎng)150 m，總水深為1 m，密度分布按照式(18)設(shè)定，其中ρ1為998 kg/m3、ρ2為1 025 kg/m3．計(jì)算域內(nèi)網(wǎng)格尺寸：水平向步長(zhǎng)取λA/80；垂向步長(zhǎng)取H/15．

表3 算例列表Tab.3 Detailed settings of the cases

4.2 數(shù)值模擬結(jié)果與理論解的對(duì)比及分析

為驗(yàn)證本文密度連續(xù)變化內(nèi)孤立波模型(iswFoam)的造波能力，選取不同時(shí)刻的波面場(chǎng)圖與理論解進(jìn)行對(duì)比．圖4 以工況1 為例給出了密度等值線分布，圖中黑色實(shí)線為浮頻率最大時(shí)對(duì)應(yīng)的等密度線，將該等密度線作為本文分析的內(nèi)孤立波波面．

圖5 給出了表3 中4 種工況連續(xù)分層內(nèi)孤立波波面數(shù)值造波結(jié)果，并與相同條件下的eKdV 理論解進(jìn)行了比較．結(jié)果表明，數(shù)值造波結(jié)果與理論解基本一致，波高的誤差百分比小于2%，本文模型對(duì)于不同上下層水深比，不同波高條件，均能造出符合理論解的內(nèi)孤立波．圖6 給出了工況4 不同時(shí)刻的密度場(chǎng)圖，由圖可知內(nèi)孤立波在傳播過(guò)程中，波面穩(wěn)定傳播，且清晰無(wú)破碎．

圖4 密度等值線分布Fig.4 Schematic diagram of density contours

圖6 不同時(shí)刻密度場(chǎng)圖(工況4：h1/h2＝3/7，a＝16.5 cm)Fig.6 Density contours of different time(Case 4：h1/h2＝3/7，a＝16.5 cm)

圖7 密度連續(xù)變化模型與兩層流體模型的密度剖面及浮頻率剖面Fig.7 Density profile and buoyancy frequency profile in continuously stratified and two-layer fluid model

圖8 密度連續(xù)模型與兩層流體模型內(nèi)孤立波傳播數(shù)值模擬結(jié)果Fig.8 Comparison of numerical simulation results of internal solitary wave propagation in continuously stratified and two-layer fluid model

4.3 密度連續(xù)變化與兩層流體模型的對(duì)比

為了比較密度連續(xù)變化和強(qiáng)分層模型的不同，利用OpenFAOM 多相流模型對(duì)密度強(qiáng)分層內(nèi)孤立波進(jìn)行模擬．圖7 給出了密度連續(xù)變化與兩層流體模型的密度剖面及浮頻率剖面的對(duì)比情況，其中流體浮頻率為1 000 kg/m3．對(duì)于強(qiáng)分層模型，除了密度分布不同，算例其他設(shè)置均與本文模型相同，圖8 以工況2 和工況4 為例，給出了兩種模型的模擬結(jié)果．圖8 表明，受密度躍層的影響，密度連續(xù)變化模型的內(nèi)孤立波波速(cpyc)小于強(qiáng)分層模型中的內(nèi)孤立波波速(ctwo)，即密度躍層的存在使內(nèi)孤立波波速相對(duì)于密度強(qiáng)分層條件有所減小．兩層流體的密度強(qiáng)分層模型本質(zhì)上忽略了實(shí)際海洋環(huán)境下的密度躍層厚度，可能過(guò)高估計(jì)了實(shí)際海洋中內(nèi)孤立波傳播速度．該結(jié)論與Kao等[10]和Segur 等[11]的實(shí)驗(yàn)結(jié)果一致，即相較于強(qiáng)分層模型，有限厚度密度躍層的存在會(huì)削減內(nèi)孤立波的傳播速度．Thiem 等[35]也得到類似的結(jié)論，指出兩層流體模型的內(nèi)孤立波波速結(jié)果高于實(shí)測(cè)結(jié)果．值得注意的是，雖然本文密度垂向分布選用雙曲正切型式，但是通過(guò)Segur 等[11]的實(shí)驗(yàn)分析可知只要存在密度躍層，兩層流體中的內(nèi)孤立波波速就會(huì)被削減，所以采用其他型式的密度垂向分布也會(huì)得到與本文相同的結(jié)果．

對(duì)于同一種工況，除了密度方面的不同，兩種模型的初始條件、算例設(shè)置及網(wǎng)格均相同，所以相同工況下初始時(shí)刻的總能量和傳播過(guò)程中的數(shù)值耗散是相同的．在總能量和數(shù)值耗散相同的情況下，兩層流體模型的內(nèi)孤立波衰減大于本文密度連續(xù)變化模型(圖8)，表明兩層流體模型由于忽略密度躍層厚度，可能導(dǎo)致其過(guò)高估計(jì)實(shí)際海洋環(huán)境下內(nèi)孤立波傳播過(guò)程中的能量耗散及波高衰減．

4.4 紊流影響分析

為了分析前文各個(gè)工況下紊流對(duì)內(nèi)孤立波傳播演變過(guò)程的影響，圖9 給出了工況2 所對(duì)應(yīng)的紊流動(dòng)能耗散率．由計(jì)算結(jié)果可知，內(nèi)孤立波穩(wěn)定傳播時(shí)，最大紊流動(dòng)能耗散率出現(xiàn)在波谷附近的躍層內(nèi)，量級(jí)為10-6，紊流對(duì)內(nèi)孤立波的穩(wěn)定傳播以及密度場(chǎng)的變化幾乎無(wú)影響，這一結(jié)論與Small 等[24]研究?jī)?nèi)孤立波時(shí)給出的結(jié)論一致．需要指出的是，本文主要對(duì)內(nèi)孤立波的穩(wěn)定傳播過(guò)程進(jìn)行了模擬，在研究淺水區(qū)內(nèi)波破碎等過(guò)程時(shí)紊流的影響是不可忽略的．

圖9 紊流動(dòng)能耗散率結(jié)果(工況2，t＝100 s，藍(lán)色實(shí)線為等密度線)Fig.9 Graph of results of turbulent energy dissipation rate(Case 2，t＝100 s，the solid blue line represents the iso-density line)

5 結(jié) 論

本文基于雷諾時(shí)均Navier-Stokes 方程，加入密度輸運(yùn)方程，將內(nèi)孤立波理論解作為初始場(chǎng)，建立了密度連續(xù)變化內(nèi)孤立波數(shù)值模型，并對(duì)其合理性進(jìn)行了驗(yàn)證．通過(guò)引入有效勢(shì)能概念，將內(nèi)孤立波的無(wú)限波長(zhǎng)合理簡(jiǎn)化為有效波長(zhǎng)，并對(duì)數(shù)值模擬時(shí)在保證模擬精度前提下的合理網(wǎng)格尺寸進(jìn)行了分析．在此基礎(chǔ)上，探討了密度連續(xù)變化模型與兩層流體模型的差異性．主要結(jié)論如下．

(1) 本文所建模型在不同上、下層水深比(h1/h2)、不同波高條件下均能造出符合理論解及符合傳播演變規(guī)律的內(nèi)孤立波．

(2) 對(duì)于 eKdV 理論，將無(wú)量綱化有效勢(shì)能APEnon＝99.8%時(shí)所對(duì)應(yīng)的波長(zhǎng)取為有效波長(zhǎng)Aλ是合理的，能充分反映內(nèi)孤立波特性．

(3) 綜合考慮波形分辨率和計(jì)算效率，推薦內(nèi)孤立波模擬水平向空間步長(zhǎng)取1/80 有效波長(zhǎng)，垂向空間步長(zhǎng)取1/15 內(nèi)孤立波波高．

(4) 兩層流體模型由于忽略密度躍層厚度，導(dǎo)致其可能過(guò)高估計(jì)實(shí)際海洋環(huán)境下內(nèi)孤立波的傳播速度和波高衰減．

本文所建模型考慮了密度的連續(xù)變化和密度躍層的厚度，更接近實(shí)際海洋環(huán)境情況，能更好地反映實(shí)際海洋環(huán)境下內(nèi)孤立波的各項(xiàng)特性．