分層流中潛艇共振激發內孤波的傳播

楊俊彥，杜永成，楊 立

（海軍工程大學 動力工程學院,武漢 430033）

0 引 言

1834年8月，John Scott Russell首先發現了在自由表面上傳播的孤立波，這是當時傳統Airy理論所無法描述的。1895年，Korteweg與Devries發現了孤立波傳播滿足的非線性模型，描述了內孤波的演化過程。孤立波理論在內波研究中的應用始于20世紀60年代，1966年，Benjamin[1]首先導出了淺水界面波單向傳播的KdV 方程，這一方程對小振幅內波的傳播誤差較??；針對在深海中傳播的內波，Benjamin（1967）[2]和Ono（1975）[3]給出了在無限深流體中界面波傳播的Benjamin-Ono 方程，簡稱BO 方程；Joseph（1975）[4]推導了一個有限深度的內波傳播方程。以上模型均屬于一維內孤波傳播模型。Kubota（1978）[5]和Lee(1974)[6]等人提出了中長波內波的傳播方程，即ILW 方程。當水深趨于淺水時，其解趨于KdV方程；當水深趨于深水時，其解趨于BO方程。

對海洋中內孤波的傳播和演化，人們常用KdV 方程或淺水有限深度方程來研究。由于KdV 方程能很好地處理非線性效應和色散效應之間的平衡，且方程簡單，許多水下孤立波傳播的數值模型都是基于KdV 理論。Koop 與Butler、Segur 和Hammack 檢驗了KdV 方程和有限水深方程的有效性，發現即使在KdV 方程的假設條件僅有少量滿足的條件下，它仍具有很好的適用性。Shishkina(1994)[7]在熱分層水槽中，實驗研究了細長拖曳模型共振激發的內孤波，發現激發的內孤波與KdV方程描述的內孤波波形是一致的。

在實際海洋中航行的潛艇，由于潛艇激發的內波幅值較小[8]，大約為10 米量級，高階非線性影響較小，本文采用KdV方程來研究潛艇激發內波的傳播特性。

1 非線性內波的KdV方程模型

1.1 KdV方程內波模型

在Boussinesq近似下，1996年，Lamb和Yan[9]給出了無粘不可壓縮流體的二維控制方程組：

式中：u(x,z,t)、w(x,z,t)分別代表水平與垂向流速；t為時間；g代表重力加速度，正向為垂直向上方向；ρ代表無量綱的密度，相應的有量綱密度為ρ0(1+ρ)；ρ0表示為參考密度；p為無量綱壓強，對應的有量綱的壓強為ρ0g(H-z) +ρ0p，H代表水的深度。在海底，取z= 0；在z=H處，認為是剛性海面。

在公式（1）與（2）中消除壓強p，導出的流函數Ψ(x,z,t)和密度擾動控制方程為[10]

式中，

密度ρ可表示為背景的平均密度與擾動密度ρ'之和，N代表環境的浮頻率。對于小振幅長內波，流函數Ψ(x,z,t)可用小參數漸進展開[10]：

式中，μ=H2/L2，ε代表內波波動振幅與深度之比，η(x,t)為內波的幅值，Φ(z)為內波垂向位移的模態函數，上標（i,j）代表對應于O(εi μj)項的垂直結構函數。將式（8）代入式（5）和式（6），可得求解流函數Ψ(x,z,t)幅值的方程，即EKdV方程[11]：

式中，c為長內波相速度（long wave phase speed），α表示平方非線性系數，α1代表立方非線性系數或稱高階非線性系數，β表示色散系數。對潛艇激發的內波，由于波幅較小，可忽略高階非線性項的影響，當α1=0時，方程（9）變為KdV方程：

內波垂向位移的模態函數Φ(z)可由下式計算[10]：

1.2 非線性和頻散系數計算

2004年，Grimshaw等[11]給出了公式（10）中非線性系數和色散系數的計算公式：

根據蔡樹群和甘子鈞（1995）[12]在無Baussinesq 近似下的Thomson-Haskell 方法（簡稱T-S 法）可求解得到Φ(z)的數值解和相速度c。由Φ(z)和方程（12）～（13），可計算得到非線性系數和色散系數。對常系數KdV方程，其理論解為

式中，η0為孤立波的振幅。

2 KdV方程數值求解方法

求解非線性KdV 方程的數值解法有多種，如有限差分法和偽譜法（pseudospectral method）等。1965年，Zabusky與Kruskal[15]采用有限差分方法給出了KdV方程的數值解。

本文采用全隱有限差分方法數值求解KdV方程，對式（10）的時間項取

式中，j為時間序號，i為單元序號。式（10）中的一、三階導數項取

式（10）中的η項取為

將公式（16）～（19）代入式（10），可得KdV方程的全隱式差分格式為

式中，Δt為時間步長，Δx為空間步長。

3 潛艇激發內波的傳播

3.1 內波傳播模型驗證

Shishkina（1996）[7]在熱分層環境中進行了孤立波的共振成生實驗，實驗水槽中的溫度、密度和浮頻率分布如圖1所示。當模型速度與第二模態內波相速度相等時，將產生第二模態孤立波，實驗得到的激發內波頻率為0.100 6 rad/s。實驗測量得到的孤立波波速為V2=0.013 6 m/s，內孤波幅值為η0=0.006 m,內波的特征半倍波寬為L2=0.6 m。

圖1 Shishkina（1996）實驗的溫度、密度和浮頻率分布Fig.1 Distribution of temperature,density and buoyancy frequency in Shishkina's（1996）experiment

根據內波理論可計算出在該實驗條件下的前3模態的水平波數kh為1.282、7.367 7、10.268 rad/m,相速度ch為0.078 4、0.013 6、0.009 8 m/s，如圖2 所示。計算得到的各模態內波波函數見圖3。根據內波KdV 方程非線性系數和色散系數計算公式（12）和（13），可計算得到內波第一至第三模態非線性系數α分別為0.025、1.68和0.071 6,色散系數β分別為0.017 6、0.000 232和0.000 155。

圖2 由TH方法計算得到的內波水平相速與波數的關系Fig.2 Relationship between horizontal phase velocity of internal wave and wave number calculated by TH method

圖3 由TH方法計算得到的內波波函數Fig.3 Wave function of internal wave calculated by TH method

3.2 內波的傳播計算與討論

對實驗激發的第二模態內孤波，計算得到非線性系數α=1.68,色散系數β=0.000 232。由KdV 方程的孤波解，可給出初始時刻的孤波波形，利用KdV 方程計算可得到內波隨時間和空間的演化。圖4和圖5為不同時刻激發的第二模態內孤波的空間分布。由圖可見，在一定初始幅值條件下，內孤波隨時間幅值逐漸增大，在140 s 時內波幅值為0.006 5 m,與實驗在同時刻測量得到的內波幅值0.006 m 相比，基本一致，誤差為8.33%，說明利用KdV方程可以較為準確地模擬內波幅值隨時間的演化。

圖4 不同時刻內孤波的空間分布Fig.4 Spatial distribution of solitary waves at different times

圖5 內孤波隨時間空間的演化Fig.5 Evolution of solitary waves with time and space

對激發的第一模態內波,由計算得到非線性系數α=0.025,色散系數β=0.017 6，利用KdV 方程計算可得到內波隨時間和空間的演化。圖6和圖7為第一模態內波隨時間與空間分布。由圖可見,對初始幅值為一定值的內孤波，隨著時間的演化，內孤波的幅值逐漸減小，并產生色散效應，與實驗觀測結果一致，如圖8 所示[14]。對第三模態內波，由于非線性效應較弱，色散效應相對較強，內波的時間空間演化規律與第一模態內波類似。

圖6 內波隨時間空間的演化Fig.6 Evolution of internal waves with time and space

圖7 不同時刻內波的空間分布Fig.7 Spatial distribution of internal waves at different times

圖8 Koop等實驗得到的色散內波波列[14]Fig.8 Dispersiver wave train obtained by Koop et al[14]

4 結 語

當潛艇航速與內波各模態的水平波速相等時，潛艇運動將共振激發內孤波。本文考慮色散效應、非線性效應的影響，建立了潛艇內波傳播的非線性KdV 數學模型，采用TH 數值方法求解內波各模態的波函數，由已知的當地浮頻率計算出內波非線性KdV 數學模型中的色散系數和非線性系數。采用有限差分方法求解KdV 方程，獲得了潛艇內波的傳播規律。利用前人研究得到的實驗數據對潛艇激發的內孤波的傳播進行了驗證，計算得到的內波傳播相速度和內波幅值與實驗測量結果的誤差均小于20%，說明KdV模型能有效描述潛艇激發的內孤波的傳播。

本文采用了無粘不可壓縮流體二維控制方程研究潛艇激發的內孤波的傳播。由于潛艇激發的內波通常是三維的，下一步將進一步開展潛艇激發三維內波的傳播規律研究。