楊俊彥，杜永成，楊 立

（海軍工程大學(xué) 動(dòng)力工程學(xué)院,武漢 430033）

0 引 言

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

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

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

1 非線性內(nèi)波的KdV方程模型

1.1 KdV方程內(nèi)波模型

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

式中：u(x,z,t)、w(x,z,t)分別代表水平與垂向流速；t為時(shí)間；g代表重力加速度，正向?yàn)榇怪毕蛏戏较颍沪汛頍o(wú)量綱的密度，相應(yīng)的有量綱密度為ρ0(1+ρ)；ρ0表示為參考密度；p為無(wú)量綱壓強(qiáng)，對(duì)應(yīng)的有量綱的壓強(qiáng)為ρ0g(H-z) +ρ0p，H代表水的深度。在海底，取z= 0；在z=H處，認(rèn)為是剛性海面。

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

式中，

密度ρ可表示為背景的平均密度與擾動(dòng)密度ρ'之和，N代表環(huán)境的浮頻率。對(duì)于小振幅長(zhǎng)內(nèi)波，流函數(shù)Ψ(x,z,t)可用小參數(shù)漸進(jìn)展開(kāi)[10]：

式中，μ=H2/L2，ε代表內(nèi)波波動(dòng)振幅與深度之比，η(x,t)為內(nèi)波的幅值，Φ(z)為內(nèi)波垂向位移的模態(tài)函數(shù)，上標(biāo)（i,j）代表對(duì)應(yīng)于O(εi μj)項(xiàng)的垂直結(jié)構(gòu)函數(shù)。將式（8）代入式（5）和式（6），可得求解流函數(shù)Ψ(x,z,t)幅值的方程，即EKdV方程[11]：

式中，c為長(zhǎng)內(nèi)波相速度（long wave phase speed），α表示平方非線性系數(shù)，α1代表立方非線性系數(shù)或稱高階非線性系數(shù)，β表示色散系數(shù)。對(duì)潛艇激發(fā)的內(nèi)波，由于波幅較小，可忽略高階非線性項(xiàng)的影響，當(dāng)α1=0時(shí)，方程（9）變?yōu)镵dV方程：

內(nèi)波垂向位移的模態(tài)函數(shù)Φ(z)可由下式計(jì)算[10]：

1.2 非線性和頻散系數(shù)計(jì)算

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

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

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

2 KdV方程數(shù)值求解方法

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

本文采用全隱有限差分方法數(shù)值求解KdV方程，對(duì)式（10）的時(shí)間項(xiàng)取

式中，j為時(shí)間序號(hào)，i為單元序號(hào)。式（10）中的一、三階導(dǎo)數(shù)項(xiàng)取

式（10）中的η項(xiàng)取為

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

式中，Δt為時(shí)間步長(zhǎng)，Δx為空間步長(zhǎng)。

3 潛艇激發(fā)內(nèi)波的傳播

3.1 內(nèi)波傳播模型驗(yàn)證

Shishkina（1996）[7]在熱分層環(huán)境中進(jìn)行了孤立波的共振成生實(shí)驗(yàn)，實(shí)驗(yàn)水槽中的溫度、密度和浮頻率分布如圖1所示。當(dāng)模型速度與第二模態(tài)內(nèi)波相速度相等時(shí)，將產(chǎn)生第二模態(tài)孤立波，實(shí)驗(yàn)得到的激發(fā)內(nèi)波頻率為0.100 6 rad/s。實(shí)驗(yàn)測(cè)量得到的孤立波波速為V2=0.013 6 m/s，內(nèi)孤波幅值為η0=0.006 m,內(nèi)波的特征半倍波寬為L(zhǎng)2=0.6 m。

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

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

圖2 由TH方法計(jì)算得到的內(nèi)波水平相速與波數(shù)的關(guān)系Fig.2 Relationship between horizontal phase velocity of internal wave and wave number calculated by TH method

圖3 由TH方法計(jì)算得到的內(nèi)波波函數(shù)Fig.3 Wave function of internal wave calculated by TH method

3.2 內(nèi)波的傳播計(jì)算與討論

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

圖4 不同時(shí)刻內(nèi)孤波的空間分布Fig.4 Spatial distribution of solitary waves at different times

圖5 內(nèi)孤波隨時(shí)間空間的演化Fig.5 Evolution of solitary waves with time and space

對(duì)激發(fā)的第一模態(tài)內(nèi)波,由計(jì)算得到非線性系數(shù)α=0.025,色散系數(shù)β=0.017 6，利用KdV 方程計(jì)算可得到內(nèi)波隨時(shí)間和空間的演化。圖6和圖7為第一模態(tài)內(nèi)波隨時(shí)間與空間分布。由圖可見(jiàn),對(duì)初始幅值為一定值的內(nèi)孤波，隨著時(shí)間的演化，內(nèi)孤波的幅值逐漸減小，并產(chǎn)生色散效應(yīng)，與實(shí)驗(yàn)觀測(cè)結(jié)果一致，如圖8 所示[14]。對(duì)第三模態(tài)內(nèi)波，由于非線性效應(yīng)較弱，色散效應(yīng)相對(duì)較強(qiáng)，內(nèi)波的時(shí)間空間演化規(guī)律與第一模態(tài)內(nèi)波類似。

圖6 內(nèi)波隨時(shí)間空間的演化Fig.6 Evolution of internal waves with time and space

圖7 不同時(shí)刻內(nèi)波的空間分布Fig.7 Spatial distribution of internal waves at different times

圖8 Koop等實(shí)驗(yàn)得到的色散內(nèi)波波列[14]Fig.8 Dispersiver wave train obtained by Koop et al[14]

4 結(jié) 語(yǔ)

當(dāng)潛艇航速與內(nèi)波各模態(tài)的水平波速相等時(shí)，潛艇運(yùn)動(dòng)將共振激發(fā)內(nèi)孤波。本文考慮色散效應(yīng)、非線性效應(yīng)的影響，建立了潛艇內(nèi)波傳播的非線性KdV 數(shù)學(xué)模型，采用TH 數(shù)值方法求解內(nèi)波各模態(tài)的波函數(shù)，由已知的當(dāng)?shù)馗☆l率計(jì)算出內(nèi)波非線性KdV 數(shù)學(xué)模型中的色散系數(shù)和非線性系數(shù)。采用有限差分方法求解KdV 方程，獲得了潛艇內(nèi)波的傳播規(guī)律。利用前人研究得到的實(shí)驗(yàn)數(shù)據(jù)對(duì)潛艇激發(fā)的內(nèi)孤波的傳播進(jìn)行了驗(yàn)證，計(jì)算得到的內(nèi)波傳播相速度和內(nèi)波幅值與實(shí)驗(yàn)測(cè)量結(jié)果的誤差均小于20%，說(shuō)明KdV模型能有效描述潛艇激發(fā)的內(nèi)孤波的傳播。

本文采用了無(wú)粘不可壓縮流體二維控制方程研究潛艇激發(fā)的內(nèi)孤波的傳播。由于潛艇激發(fā)的內(nèi)波通常是三維的，下一步將進(jìn)一步開(kāi)展?jié)撏Ъぐl(fā)三維內(nèi)波的傳播規(guī)律研究。