MCC-FS內孤立波模型算法及速度場研究

趙彬彬， 張天鈺， 王戰,2， 段文洋， 王澤航 (.哈爾濱工程大學 船舶工程學院，黑龍江 哈爾濱 5000； 2.青島哈爾濱工程大學創新發展中心，山東 青島 266555)

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

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

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

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

1 保留自由面的MCC模型

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

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

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

(1)

(2)

(3)

(4)

式中：ζ1表示自由面位移；ζ2表示內波界面位移。

η1=h1+ζ1-ζ2

(5)

η2=h2+ζ2

(6)

(7)

(8)

Gi作為一個運算規則，其定義為：

(9)

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

(10)

(11)

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

更多關于MCC-FS模型原始算法的內容參考文獻[5]。

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

2 MCC-FS模型的數值算法

在MCC-FS模型的本文算法中，將大地坐標系轉換到隨著內波面相速度平動的穩態坐標系oxz：

(12)

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

(13)

(14)

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

(15)

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

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

(16)

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

(17)

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

(18)

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

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

1)對于上層流體：

(19)

2)對于下層流體：

(20)

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

(21)

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

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

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

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

3.1 內孤立波速度場研究

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的內孤立波進行計算。內孤立波波形結果如圖3所示。圖3中，本文算法與Kodaira等[5]的原始算法得到的波形結果完全重合，與實驗值相比，波形結果也很準確。

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

圖4為不同波幅內孤立波的波速曲線，其中，速度無量綱化參數c0為：

圖4 內孤立波波速(ρ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]得到的內孤立波波速結果也完全重合。需要指出的是，對于該算例，Kodaira等[5]并未給出速度場的數值結果與實驗結果。

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

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

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

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

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

3.2 不同密度比的淺水模式內孤立波特性研究

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

圖9中，h1/h2=1/5,a/h1=-0.5時，上、下層密度越接近，得到的內孤立波波形越寬，上、下層密度相差越大，得到的內孤立波波形越窄。圖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。這一結論同樣也解釋了對于上、下層流體密度接近的情況，引入剛蓋假設是合理的，而在上、下層密度相差較大時需要考慮自由面的影響。

圖9 不同密度比內孤立波波形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時隨著密度相差越大，上、下層流體質點的水平速度變得更大，內波界面上的速度跳躍值也越大。圖12中，h1/h2=1/5,a/h1=-0.5時，隨著上、下層密度比的增大，內孤立波的波速變小，上下層密度比為0.8到0.95時，波速隨密度比近似成線性變化，密度比為0.95～1時，隨著密度比增大，波速減小的較快。

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

4 結論

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

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

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