MCC-FS內(nèi)孤立波模型算法及速度場研究

趙彬彬， 張?zhí)焘暎?王戰(zhàn),2， 段文洋， 王澤航 (.哈爾濱工程大學(xué) 船舶工程學(xué)院，黑龍江 哈爾濱 5000； 2.青島哈爾濱工程大學(xué)創(chuàng)新發(fā)展中心，山東 青島 266555)

內(nèi)孤立波是一種最大振幅發(fā)生在密度穩(wěn)定層化海洋內(nèi)部的波動，通常由洋流通過不平坦的海底而產(chǎn)生，是海洋中頻繁發(fā)生的現(xiàn)象[1-2]。我國南海是內(nèi)孤立波的頻發(fā)海域,被評價為研究大幅內(nèi)波的“天然實驗室”，Huang等[3]在南海觀測到波幅240 m，流速2.55 m/s的大幅內(nèi)孤立波。大幅內(nèi)孤立波在傳播中攜帶著巨大能量，對海洋結(jié)構(gòu)物存在嚴重威脅。因此，內(nèi)孤立波的研究對于南海海洋環(huán)境的認識和開發(fā)具有極大的現(xiàn)實意義。在內(nèi)波問題的研究上，當密度躍層的厚度較小時，可以簡化為兩層流體間的內(nèi)波問題，即每一層流體的密度均為常數(shù)，兩層流體界面處存在密度跳躍。在兩層內(nèi)孤立波的水池實驗研究上，文獻[4-7]采用重力塌陷法對兩層流體間內(nèi)孤立波的波形、波速以及速度場等特性進行研究。黃文昊等[8]采用雙推板法生成內(nèi)孤立波，對實驗中內(nèi)孤立波的振幅實現(xiàn)了自主控制。國際上流行的內(nèi)孤立波理論模型主要有Korteweg-deVries (KdV)模型、MCC模型以及High-Level Green-Naghdi (HLGN)模型等。KdV模型是當前應(yīng)用最廣泛的內(nèi)波模型，但是對大幅內(nèi)孤立波進行研究時，KdV模型給出的結(jié)果往往不盡如人意[4]。

考慮到大幅內(nèi)孤立波的強非線性，文獻[9-10]假設(shè)特征波長較總水深更長，振幅與總水深相當，從完全非線性歐拉方程出發(fā)，采用層平均水平速度方法，建立了適用于剛蓋近似下大幅內(nèi)孤立波的MCC-RL模型。在MCC-RL模型中，由于沒有引入小幅波假定，同時其形式簡單，因此被廣泛地應(yīng)用到淺水模式下的大幅內(nèi)孤立波問題的研究中，并展示了良好的性能[7-8,10-11]。

當上、下層流體密度相差較大時，自由面對內(nèi)孤立波特性的影響不可忽略，引入剛蓋近似假設(shè)是不合理的。文獻[5,12]發(fā)現(xiàn)考慮自由面效應(yīng)，MCC-FS模型得到的波形結(jié)果相對MCC-RL模型要更窄，同樣更接近實驗值。然而，Kodaira等[5]求解MCC-FS模型時，在MCC-FS模型的4個方程中消去上、下層平均速度，推導(dǎo)得到只保留自由面以及內(nèi)波界面2個未知數(shù)的2個非線性微分方程，求解出自由面和內(nèi)波波面之后，可計算上、下層的層平均速度。本文中則是對MCC-FS模型4個方程進行直接求解。另外，MCC-FS模型的速度場結(jié)果從未被展示以及驗證。

本文通過與實驗值以及他人數(shù)值結(jié)果的對比，對本文算法以及速度場結(jié)果進行了驗證。推導(dǎo)出MCC-FS模型的算法以及基于MCC-FS模型速度場結(jié)果，并研究了不同密度比對淺水模式下內(nèi)孤立波特性的影響。

1 保留自由面的MCC模型

考慮自由面的兩層流體間內(nèi)孤立波傳播示意圖如圖1所示。假定兩層流體均無粘、不可壓縮。在未受擾動時，上層流體厚度為h1，密度為ρ1；下層流體厚度為h2，密度為ρ2。建立大地直角坐標系oxz，其中，ox軸與未擾動的內(nèi)界面重合，指向內(nèi)孤立波傳播方向。η1、η2分別為上層流體與下層流體的厚度。

圖1 兩層流體間內(nèi)孤立波Fig.1 The sketch of internal solitary wave in a two-layer system

保留自由面的MCC-FS模型方程為[5]：

(1)

(2)

(3)

(4)

式中：ζ1表示自由面位移；ζ2表示內(nèi)波界面位移。

η1=h1+ζ1-ζ2

(5)

η2=h2+ζ2

(6)

(7)

(8)

Gi作為一個運算規(guī)則，其定義為：

(9)

通過消去速度項，得到了簡化后的含二階導(dǎo)數(shù)項的非線性微分方程[5]：

(10)

(11)

式中：i,j=1,2；式子中的上標撇號代表對水平方向上的X的求導(dǎo)，需給出波速c求解上式。該方程組中包含η1、η22個未知數(shù)，封閉可解。

更多關(guān)于MCC-FS模型原始算法的內(nèi)容參考文獻[5]。

圖2 原始算法與本文算法對比Fig.2 Comparison between original and the algorithm in this paper

2 MCC-FS模型的數(shù)值算法

在MCC-FS模型的本文算法中，將大地坐標系轉(zhuǎn)換到隨著內(nèi)波面相速度平動的穩(wěn)態(tài)坐標系oxz：

(12)

式(12)代入MCC-FS模型的方程(1)～(4)，可得到穩(wěn)態(tài)坐標系下的MCC-FS模型方程。采用差商代替導(dǎo)數(shù)的方法進行空間離散，在這里以內(nèi)波界面位移ζ2為例，其他各量的表示方法類似，各階導(dǎo)數(shù)可表示為：

(13)

(14)

2ζ2(X+ΔX)+ζ2(X+2ΔX)]

(15)

將內(nèi)孤立波的波峰設(shè)在X=0處，無論是內(nèi)波界面位移、自由面位移還是上、下層平均速度均滿足對稱性。同時，在無窮遠處，內(nèi)波界面位移、自由面位移以及上、下層平均速度均為0，同樣以內(nèi)波界面位移ζ2為例，邊界條件可以表示為：

ζ2(-ΔX)=ζ2(ΔX)

(16)

ζ2(-2ΔX)=ζ2(2ΔX)

(17)

ζ2(X)=0 (X→∞)

(18)

數(shù)值求解采用牛頓迭代法，當相鄰迭代步得到變量的絕對值之差小于10-7時，迭代停止。

Camassa等[11]給出了MCC-RL模型中計算流體質(zhì)點水平速度分布的近似方法，該方法可以應(yīng)用于MCC-FS模型上。

1)對于上層流體：

(19)

2)對于下層流體：

(20)

根據(jù)式(19)、(20)以及MCC-FS模型得到的自由面以及內(nèi)波波面結(jié)果，可以得到MCC-FS模型的水平速度分布。在垂向速度的計算上，Choi等[10]給出了通過層平均水平速度得到垂向速度為：

(21)

式中：i=1代表上層流體；i=2代表下層流體。本文算法可以直接得到上、下層流體的層平均速度，可以更方便地求解MCC-FS模型的垂向速度。

在數(shù)值計算中，MCC-FS模型本文算法更為直接方便，在MATHMATICA中，得到數(shù)值結(jié)果的時間要小于1 min，計算耗時很短。

3 MCC-FS模型速度場算法結(jié)果與分析

本文基于Kodaira等[5]的算例，比較了本文算法、原始算法、實驗值得到的內(nèi)孤立波的波形以及波速結(jié)果，以驗證本文算法的準確性。并基于鄒麗等[7]的PIV實驗結(jié)果，對本文算法以及給出的MCC-FS模型速度場求解方法進行進一步驗證。最后，研究了不同密度比對內(nèi)孤立波特性的影響。

3.1 內(nèi)孤立波速度場研究

Kodaira等[5]的算例中，上層流體的密度為ρ1=856 kg/m3，厚度為h1=0.05 m；下層流體的密度為ρ2=996 kg/m3，厚度為h2=0.25 m，對幅值為a/h1為-0.50、-0.77、-0.99以及-1.21的內(nèi)孤立波進行計算。內(nèi)孤立波波形結(jié)果如圖3所示。圖3中，本文算法與Kodaira等[5]的原始算法得到的波形結(jié)果完全重合，與實驗值相比，波形結(jié)果也很準確。

圖3 內(nèi)孤立波波形(ρ1/ρ2=0.859,h1/h2=1/5)Fig.3 Profiles of internal solitary waves (ρ1/ρ2=0.859,h1/h2=1/5)

圖4為不同波幅內(nèi)孤立波的波速曲線，其中，速度無量綱化參數(shù)c0為：

圖4 內(nèi)孤立波波速(ρ1/ρ2=0.859,h1/h2=1/5)Fig.4 Speed of internal solitary waves (ρ1/ρ2=0.859,h1/h2=1/5)

(22)

可以看出，本文算法與Kodaira等[5]得到的內(nèi)孤立波波速結(jié)果也完全重合。需要指出的是，對于該算例，Kodaira等[5]并未給出速度場的數(shù)值結(jié)果與實驗結(jié)果。

鄒麗等[7]的算例中，上層流體的密度為ρ1=941 kg/m3，厚度為h1=0.05 m；下層流體的密度為ρ2=1 003 kg/m3，厚度為h2=0.25 m；幅值為a/h1=-0.875,-1.323的內(nèi)孤立波波形的實驗結(jié)果，MCC-FS模型本文算法與實驗結(jié)果對比如圖5所示。在該實驗布置下，MCC-FS模型本文算法得到的波形結(jié)果同實驗值基本一致。

圖5 內(nèi)孤立波波形(ρ1/ρ2=0.938,h1/h2=1/5)Fig.5 Profiles of internal solitary waves (ρ1/ρ2=0.938,h1/h2=1/5)

鄒麗等[7]基于PIV技術(shù)得到了內(nèi)孤立波的速度場結(jié)果。其中，實驗中內(nèi)孤立波幅值分別為a/h1=-0.24,-0.59,-0.70,-0.96。需要說明的是，鄒麗等[7]并未直接給出內(nèi)孤立波的波幅，本文中內(nèi)孤立波的波幅是對文章中給出的內(nèi)孤立波波面描點后得到的。圖6中，基于MCC-FS模型本文算法得到的速度場與鄒麗等[7]得到的實驗速度場基本一致。

針對幅值為a/h1=-0.59的內(nèi)孤立波，鄒麗等[7]給出了不同深度下水平流速時歷曲線以及波峰豎直剖面水平速度分布的實驗結(jié)果。

圖7中，ρ1/ρ2=0.938,h1/h2=1/5,a/h1=-0.59時，MCC-FS模型本文算法預(yù)報的不同深度的流體質(zhì)點水平速度時歷曲線與鄒麗等[7]的實驗結(jié)果很接近。

圖8中，MCC-FS模型本文算法得到的波峰所在豎直剖面水平速度分布的數(shù)值結(jié)果同實驗值吻合的很好，尤其是下層流體速度的預(yù)報上，基本與實驗值重合。

3.2 不同密度比的淺水模式內(nèi)孤立波特性研究

在不同密度比內(nèi)孤立波特性的研究上，算例的布置為上層流體的厚度h1=0.05 m，下層流體的厚度h2=0.25 m，下層流體的密度ρ2=1 000 kg/m3，改變上層流體的密度ρ1研究不同密度比的影響，上層密度ρ1=800,900,999 kg/m3，內(nèi)孤立波幅值a=-0.025 m。這里需要說明的是，由于速度的無量綱化參數(shù)c0在不同的上、下層密度比是不同的，因此，接下來展示的結(jié)果均有量綱。

圖9中，h1/h2=1/5,a/h1=-0.5時，上、下層密度越接近，得到的內(nèi)孤立波波形越寬，上、下層密度相差越大，得到的內(nèi)孤立波波形越窄。圖10中，h1/h2=1/5,a/h1=-0.5時，在上、下層密度相差不大時，自由面波形幅值很小，隨著上、下層密度相差的越大，得到的自由面幅值越大。上、下層密度比為0.999、0.9以及0.8時，自由面最大位移ζ1/h1分別為0.003、0.033 2、0.071 9。這一結(jié)論同樣也解釋了對于上、下層流體密度接近的情況，引入剛蓋假設(shè)是合理的，而在上、下層密度相差較大時需要考慮自由面的影響。

圖9 不同密度比內(nèi)孤立波波形Fig.9 Profiles of internal solitary wave with different density ratios

圖10 不同密度比自由面波形Fig.10 Profiles on the free surface with different density ratios

圖11中，h1/h2=1/5,a/h1=-0.5時隨著密度相差越大，上、下層流體質(zhì)點的水平速度變得更大，內(nèi)波界面上的速度跳躍值也越大。圖12中，h1/h2=1/5,a/h1=-0.5時，隨著上、下層密度比的增大，內(nèi)孤立波的波速變小，上下層密度比為0.8到0.95時，波速隨密度比近似成線性變化，密度比為0.95～1時，隨著密度比增大，波速減小的較快。

圖12 不同密度比內(nèi)孤立波波速Fig.12 Wave speed of internal solitary waves with different density ratios

4 結(jié)論

1)通過與原算法的結(jié)果、文獻[5，7]實驗的波形結(jié)果進行比較，表明了本文提出的本文算法準確可靠，本文算法可同時求得內(nèi)波波形、自由表面波波形以及上、下層平均速度，更為直接、方便。

2)通過與文獻[7]實驗的速度場結(jié)果對比，證明了MCC-FS模型可以準確的預(yù)報淺水模式內(nèi)孤立波的速度場。

3)基于MCC-FS模型的速度場結(jié)果，發(fā)現(xiàn)內(nèi)孤立波波峰位置的豎直剖面上，隨著上、下層密度差增大，上、下層流體質(zhì)點的水平速度變得更大，內(nèi)波界面上的速度跳躍值也增大。