李永剛， 鄒麗,2,3， 胡英杰， 張九鳴， 裴玉國

(1.大連理工大學(xué) 船舶工程學(xué)院，遼寧 大連 116024； 2.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，遼寧 大連 116024； 3.高技術(shù)船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心， 上海 200240)

內(nèi)波屬于重力波，是發(fā)生在密度穩(wěn)定層化海水內(nèi)部的一種波動[1]。內(nèi)孤立波是一種特殊內(nèi)波，只有一個波峰或波谷，在傳播過程中波形恒定，能量集中，是一種災(zāi)害性的海洋環(huán)境因素。

衛(wèi)星遙感和現(xiàn)場觀測表明，內(nèi)孤立波相互作用是海洋中的常見現(xiàn)象[2-4]。當(dāng)2個孤立波幾乎沿同一方向傳播，相互作用時間較長，為“強”相互作用過程[5]，孤立波追趕是“強”相互作用中的一種特殊情況。目前針對內(nèi)孤立波追趕的研究較少，僅部分學(xué)者對表面孤立波追趕進行了研究。Lax[6]基于KdV方程，依據(jù)2孤立波初始波幅之比將孤立波追趕分成了3類。Weidman等[7]在Lax[6]分類范圍內(nèi)進行了實驗研究。Zabusky等[8]發(fā)現(xiàn)孤立波追趕后波形保持其初始形狀，但存在相位變化。Mirie等[9]通過數(shù)值研究也發(fā)現(xiàn)孤立波追趕后大波存在相位超前，小波存在相位滯后，但2個孤立波波幅可以恢復(fù)到原來的99%，尾波要比孤立波直面碰撞時小得多。Fenton等[10]利用傅里葉級數(shù)對孤立波追趕進行了數(shù)值求解，發(fā)現(xiàn)Lax[3]的分類與其數(shù)值結(jié)果存在差別。Craig等[11]利用實驗和數(shù)值模擬也觀察到了這種類似的差異，并對Lax[6]的分類結(jié)果進行了修正。Tong等[12]基于HPC方法對孤立波追趕進行了數(shù)值模擬，并驗證了Craig等[11]的分類結(jié)果。

此外，內(nèi)波在傳播過程中對海洋表面的輻聚和輻散作用會產(chǎn)生特有的表面特征，合成孔徑雷達可以捕捉到此自由表面波，這是發(fā)現(xiàn)、探測內(nèi)波的重要依據(jù)，因此有必要揭示內(nèi)孤立波追趕時的自由表面特征[12-16]。

鑒于內(nèi)孤立波和表面孤立波有許多共同的特征，本文采用包含自由表面效應(yīng)的多域邊界元數(shù)值模型和KdV方程Fourier譜方法離散求解分別對內(nèi)孤立波追趕進行了數(shù)值模擬，探討了2種方法中不同初始波幅比下內(nèi)孤立波追趕波形演化過程的差異，并通過考慮自由表面效應(yīng)的BEM-FS模型，揭示了內(nèi)孤立波追趕時自由表面波的表現(xiàn)形態(tài)。

1 移動計算域多域邊界元法

1.1 多域邊界元法

內(nèi)孤立波的雷諾數(shù)通常較大[17]，粘性效應(yīng)往往可以忽略不計。Zou等[18]將邊界元法進行推廣，根據(jù)流體密度的不同，將計算域劃分為2個計算子域，在無粘無旋不可壓縮的2層強分層理想流體中，對內(nèi)孤立波傳播進行了模擬，建立了包含自由表面效應(yīng)的多域邊界元數(shù)值模型(以下簡稱為BEM-FS模型)。計算域如圖1所示，圖中包含上下2層計算域，其中h1=0.1、ρ1=998和φ1分別代表上層流體水深、密度和速度勢，h2=0.9，ρ2=1 025和φ2分別代表下層流體水深、密度和速度勢，ai和ηi為界面處內(nèi)孤立波的波幅和波面高度，asi和ηsi為上層流體自由表面處內(nèi)孤立波誘導(dǎo)的自由表面波的波幅和波面高度，其中i=1, 2分別代表大、小波幅內(nèi)孤立波，x軸位于界面處，z軸垂直向上。

上下2層流體速度勢滿足拉普拉斯方程：

▽2φk=0

(1)

式中k=1, 2分別代表上層計算域和下層計算域。

左右兩側(cè)和底部均為不可穿透邊界條件：

(2)

圖1 內(nèi)孤立波追趕計算域示意Fig.1 The schematic of computational domain for overtaking of ISWs

上層流體自由表面S1處的動力學(xué)邊界條件為：

(3)

運動學(xué)邊界條件為:

(4)

在界面S2處，上下兩層流體的法向速度勢和壓力保持連續(xù)，因此滿足：

(5)

(6)

通過解析理論給定內(nèi)孤立波初始波形。當(dāng)內(nèi)孤立波波幅與淺層水深之比約大于0.4時，弱非線性的KdV理論將與實驗結(jié)果存在較大偏差[19]。因此當(dāng)內(nèi)孤立波波幅與淺層水深之比a/h1<0.4時，本文使用KdV方程理論解作為初始波形，當(dāng)內(nèi)孤立波波幅與淺層水深之比a/h1≥0.4時，采用文獻[20]考慮自由表面效應(yīng)的MCC理論給定初始波形。初始時刻上層流體上表面S1設(shè)為水平，經(jīng)過時間步更新演化，隨著內(nèi)孤立波波形趨于穩(wěn)定，會在自由表面誘導(dǎo)出相應(yīng)的自由表面波。

界面S2處，根據(jù)內(nèi)孤立波波形給定速度勢的初始法向?qū)?shù)為：

(7)

如圖1所示，大波幅內(nèi)孤立波與小波幅內(nèi)孤立波在左右兩側(cè)，兩波初始法相速度勢設(shè)為正負相同，以實現(xiàn)兩波同向傳播。在初始時刻，x1和x2處的波面波幅均小于10-5h1，且x2-x1>0.6，以確保在追趕發(fā)生前，波形傳播已穩(wěn)定。

關(guān)于BEM-FS模型的詳細描述以及具體解法，可以參照文獻[18]。

1.2 移動計算域

內(nèi)孤立波追趕為“強”相互作用過程，2個波沿同一方向傳播，波與波相對移動速度很慢[5]。因此孤立波追趕需要較長的作用時間及巨大的計算區(qū)域，本文提出移動計算域的改進方法，以實現(xiàn)在較小的固定長度計算域內(nèi)實現(xiàn)內(nèi)孤立波追趕的數(shù)值模擬。

如圖1所示，在t1時刻，當(dāng)左側(cè)大波幅內(nèi)孤立波每向右傳播一個計算網(wǎng)格長度時，下一計算時刻t1+Δt，整體計算域?qū)⒏S內(nèi)孤立波向右移動一個計算網(wǎng)格長度。在t1+Δt時刻，S4～S7邊界處速度勢認為與t1時刻近似相等，S1～S3處最左邊舍棄一個節(jié)點速度勢和法向速度，最右邊增加一個節(jié)點，其速度勢和法向速度與左側(cè)相鄰節(jié)點的速度勢和法向速度近似相等；反之，當(dāng)大波向前傳播未滿一個計算網(wǎng)格長度時，計算域不移動。當(dāng)內(nèi)孤立波離固壁S4～S7足夠遠時，移動計算域與固定計算域相比計算誤差很小。圖2為內(nèi)孤立波波谷距固壁20個單位長度時，40 s計算時刻分別采用移動計算域與固定計算域所得的內(nèi)孤立波追趕波形，由圖可知，移動計算域可以實現(xiàn)波形演化的準(zhǔn)確模擬。

圖2 40 s時移動、固定計算域下內(nèi)孤立波追趕波形對比Fig.2 The comparison of interface profiles by mobile computational domain and fixed computational domain at 40 s

1.3 網(wǎng)格和時間步無關(guān)性驗證

以本文最大計算波幅a1/h1=0.4為例，進行網(wǎng)格和時間步無關(guān)性驗證。

如圖3所示，當(dāng)無量綱化網(wǎng)格尺度Δx/h1<0.5時，波形收斂，為保證計算精度，網(wǎng)格尺寸大小選取Δx/h1=0.4。如圖4顯示，選取的時間步皆收斂，為保證計算精度，選取Δt=0.01 s作為時間步大小。

圖3 BEM-FS模型網(wǎng)格無關(guān)性驗證Fig.3 The convergence verification of profile in regard to mesh sizes for BEM-FS model

圖4 時間步無關(guān)性驗證Fig.4 The convergence verification of profile in regard to time steps

2 KdV方程Fourier譜方法離散求解

2.1 理論推導(dǎo)

在弱非線性、弱色散條件下，KdV方程的形式為[21]：

ζt+(c0+c1ζ)ζx+c2ζxxx=0

(8)

利用Fourier譜方法確定該方程離散格式，對非線性方程(8)線性化得：

(9)

(10)

根據(jù)傅里葉變換的微分性質(zhì)，對時間項采用向前差分為[22]：

(11)

采用四階龍格-庫塔法對其進行迭代求解。則傅里葉變換系數(shù)的迭代形式可表示為：

(12)

2.2 網(wǎng)格和時間步無關(guān)性驗證

為使快速傅里葉變換有效，離散點數(shù)應(yīng)盡可能選擇為2的整數(shù)次冪。以本文最大計算無量綱化波高a1/h1=0.4為例，采用不同網(wǎng)格尺寸和時間步長，在10個單位長度內(nèi)，數(shù)值模擬內(nèi)孤立波傳播了15 s，對網(wǎng)格和時間步進行了驗證，離散設(shè)置情況如表1所示。在不同離散設(shè)置情況下對內(nèi)孤立波波形與KdV理論波形進行了比較，如圖5所示，設(shè)置2～4下波形均收斂，綜合考慮時間成本，采用設(shè)置3中的網(wǎng)格尺寸和時間步長用于內(nèi)孤立波追趕的模擬。

表1 離散設(shè)置情況

圖5 Fourier譜方法網(wǎng)格與時間步無關(guān)性驗證Fig.5 The convergence verification in regard to mesh sizes and time steps for FSM

2.3 初始條件

根據(jù)KdV方程理論解，設(shè)定波面初始條件為：

ζ(X)=a1·sech2(X/λ1)+a2·

sech2((X+X1)/λ2)

(13)

式中：a1和λ1分別為大波幅孤立波波幅和特征波長；a2和λ2分別為小波幅孤立波波幅和特征波長；X1為2個內(nèi)孤立波的距離。在初始時刻，邊界處內(nèi)孤立波波幅以及X1范圍內(nèi)內(nèi)孤立波波幅均小于10-5h1，以消除壁面效應(yīng)和確保在波形穩(wěn)定之前，2波沒有發(fā)生碰撞。

3 波形演化及追趕分類

3.1 波形演化

圖6為文獻[20]、本文BEM-FS模型解和KdV方程理論解下不同無量綱化波高a/h1下內(nèi)孤立波波速的對比圖。BEM-FS模型數(shù)值解與實驗測得波速吻合較好，但KdV理論解波速整體偏小。由于KdV方程Fourier譜方法離散求解和BEM-FS數(shù)值模擬2種方法在計算內(nèi)孤立波波速上存在差異，故本文僅在各自時間軸下分別對內(nèi)孤立波追趕波形演化過程及其分類進行研究。

圖6 內(nèi)孤立波波速對比Fig.6 The comparison of wave speed of ISWs

Lax[6]基于KdV方程，將非線性演化方程與線性算子關(guān)聯(lián)，依據(jù)兩表面孤立波初始波幅之比a1/a2，將孤立波追趕分成了3類：

Craig[11]在完全歐拉方程哈密頓系統(tǒng)中利用數(shù)值方法對文獻[6]的分類進行了修正：

本文一系列波幅比內(nèi)孤立波追趕波形演化顯示，BEM-FS數(shù)值模型模擬結(jié)果和Craig[11]分類保持一致。當(dāng)2個內(nèi)孤立波波幅相差較大，波幅比3.536

圖7 BEM-FS數(shù)值模型下內(nèi)孤立波波幅分別為a1/h=0.4，a2/h=0.113追趕波形演化Fig.7 Overtaking of two ISWs of height a1/h=0.4，a2/h=0.113 by BEM-FS

圖7(b)時空分布下的波形演化過程圖更為形象地展示了類別3內(nèi)孤立波追趕過程(由于內(nèi)孤立波追趕作用時間長計算域大，為能在一張圖內(nèi)完整清晰展現(xiàn)作用過程，本文所有時空分布下的波形演化過程圖均為相對作用圖，即固定大波幅內(nèi)孤立波的水平位置)。此類別內(nèi)孤立波追趕在作用中心時刻，會有明顯的相互作用過程。

當(dāng)2個內(nèi)孤立波波幅比2.941

圖8 BEM-FS數(shù)值模型下內(nèi)孤立波波幅分別為a1/h1=0.4，a2/h1=0.133 追趕波形演化Fig.8 Overtaking of two ISWs of height a1/h1=0.4，a2/h1=0.133 by BEM-FS

當(dāng)2個內(nèi)孤立波波幅接近，波幅比a1/a2<2.941時，如圖9(a)所示，雖然在2個內(nèi)孤立波追趕過程中，也會出現(xiàn)大波波幅逐漸減小，小波波幅逐漸增大的過程，但在每一計算時刻，都存在2個明確的波谷，在作用中心時刻470.14 s會呈現(xiàn)出比類別2更為明顯分離的2個等高波谷。與類別2和類別3的本質(zhì)不同在于，如圖9(b)包含時間軸的三維波形演化過程圖所示，在整個計算過程中，大波幅內(nèi)孤立波并沒有將小波幅內(nèi)孤立波吞沒，沒有發(fā)生大波幅內(nèi)孤立波追趕并超越前方小波幅內(nèi)孤立波的明顯過程，而是大波幅內(nèi)孤立波將波幅不斷地傳遞給前方小波幅內(nèi)孤立波。在作用中心時刻過后，后方內(nèi)孤立波的波幅繼續(xù)傳遞給前方內(nèi)孤立波，最終導(dǎo)致前方內(nèi)孤立波波幅增大，后方內(nèi)孤立波波幅減小，因此在類別1內(nèi)孤立波追趕整個作用過程中，都是以這種波幅傳遞的方式完成的。

圖9 BEM-FS數(shù)值模型下內(nèi)孤立波波幅分別為a1/h1=0.4，a2/h1=0.3追趕波形演化圖Fig.9 Overtaking of two ISWs of height a1/h1=0.4，a2/h1=0.3 by BEM-FS

圖10 Fourier譜方法離散求解下內(nèi)孤立波波幅分別為a1/h1=0.4，a2/h1=0.133…追趕波形演化Fig.10 Overtaking of two ISWs of height a1/h1=0.4,a2/h1=0.133… by FSM

圖11 Fourier譜方法離散求解下內(nèi)孤立波波幅分別為a1/h1=0.4，a2/h1=0.152 67追趕波形演化Fig.11 Overtaking of two ISWs of height a1/h1=0.4,a2/h1=0.152 67 by FSM

3.2 幅值變化曲線

與內(nèi)孤立波直面碰撞出現(xiàn)約為初始波幅2倍的最大碰撞波高不同[21-23]，2個內(nèi)孤立波追趕時，對于所有類別下的波形演化，都存在大波波幅逐漸減小，小波波幅逐漸增大，追趕過程中波形的波幅始終介于初始2個內(nèi)孤立波波幅之間，Craig[11]在對孤立波追趕的數(shù)值研究中也有相同的發(fā)現(xiàn)。

以波幅比a1/a2=4內(nèi)孤立波追趕為例，對初始波幅為a1/h1=0.4內(nèi)孤立波的波幅大小隨時間變化的情況進行了擬合，圖12分別為2種方法所得擬合變化曲線。由圖可知，初始波幅為a1/h1=0.4的內(nèi)孤立波在追趕過程中，波幅先逐漸減小再逐漸增大，在作用中心時刻t3，達到波幅的最小值，本文定義此波幅為作用中心波幅a3/h1。

圖12 2種方法下幅值隨時間變化曲線Fig.12 The curves of trough versus time by two methods

對圖12做進一步分析，顯示2種方法對于模擬內(nèi)孤立波追趕存在諸多差異，具體參數(shù)對比如表2所示。在KdV方程Fourier譜方法離散求解下，內(nèi)孤立波經(jīng)過追趕相互作用后，波幅幾乎恢復(fù)到原始波幅，但邊界元法波幅存在一定的衰減。此外，表1顯示2種方法所得作用中心波幅a3/h1以及波幅演化時間也存在分歧。在圖12中，取內(nèi)孤立波初始波幅的80%和90%做水平線，分別與曲線交于2點，定義t4-t2和t5-t1為波幅演化為80%和90%的作用總時間，以作用中心時刻t3為分界點，把作用總時間分別分為作用前時間t3-t2、t3-t1和作用后時間t4-t3、t5-t3。如表2所示，對于其波幅演化到初始波幅80%和90%的作用總時間t4-t2以及t5-t1，邊界元法所得作用總時間為KdV方程Fourier譜方法離散求解的2倍左右。在t4-t2時間內(nèi)，2種方法所得的作用前時間t3-t2和作用后時間t4-t3分別近似相等，這說明此段時間內(nèi)2種方法皆顯示追趕過程幾乎是對稱的。但在t5-t1時間內(nèi)，對比作用前時間t3-t1和作用后時間t5-t3顯示，KdV方程Fourier譜方法離散求解顯示追趕過程幾乎對稱的，但BEM-FS模型顯示追趕過程呈現(xiàn)出明顯的不對稱性。

表2 2種方法作用參數(shù)對比

3.3 作用中心波幅變化曲線

圖13 BEM-FS數(shù)值模型下作用中心波幅a3/h1隨內(nèi)孤立波初始波幅比a1/a2的變化曲線Fig.13 The amplitude of acting center a3/h1 versus the amplitude ratio a1/a2 of the initial ISWs by BEM-FS

圖14 Fourier譜方法離散求解下作用中心波幅a3/h1隨初始內(nèi)孤立波波幅比a1/a2變化曲線Fig.14 The amplitude of acting center a3/h1 versus the amplitude ratio a1/a2 of the initial ISWs by FSM

3.4 自由表面波表現(xiàn)形式

Zou等[23]在非線性自由表面條件下，對內(nèi)孤立波“弱”相互作特例——內(nèi)孤立波直面碰撞作用過程中自由表面波的表現(xiàn)形態(tài)進行了研究，其研究表明內(nèi)孤立波與自由表面波的波形演化在形態(tài)和時間上存在明顯差異，即當(dāng)2個內(nèi)孤立波發(fā)生碰撞，伴隨著波高增大時，2個自由表面波也會相互融合，但其波高并非一味增大，而是呈現(xiàn)出上下起伏的現(xiàn)象。以上研究結(jié)果表明，在內(nèi)孤立波直面碰撞和追趕中，內(nèi)孤立波對自由表面波的影響是不同的。當(dāng)內(nèi)孤立波傳播時，會在自由表面誘導(dǎo)一個與內(nèi)孤立波極性相反的、波幅小于內(nèi)孤立波的類孤立波形式的自由表面波，其傳播速度與內(nèi)孤立波傳播速度近似相等[18,20,24]。內(nèi)孤立波誘導(dǎo)的自由表面波會造成海面粗糙度的變化，合成孔徑雷達可以捕捉到內(nèi)孤立波所產(chǎn)生的自由表面波，從而進行內(nèi)孤立波預(yù)報與反演[25-26]。本文在BEM-FS模型中，將上層流體上表面設(shè)置為非線性自由表面條件，得到了內(nèi)孤立波追趕時的自由表面波波形演化過程。以波幅比為a1/a2=2內(nèi)孤立波追趕為例，圖15顯示了自由表面波的波形演化，為顯示清晰，省略了部分上層流體，且自由表面波擴大了3倍，內(nèi)孤立波縮小了10倍。如圖所示，當(dāng)2個內(nèi)孤立波未發(fā)生相互作用時，自由表面波的傳播速度近似于內(nèi)孤立波的傳播速度，并以類孤立波的形式進行傳播；當(dāng)內(nèi)孤立波發(fā)生相互作用時，自由表面波也開始相互作用；當(dāng)內(nèi)孤立波以初始形狀發(fā)生分離時，自由表面波也以類孤立波的形式發(fā)生分離。因此在內(nèi)孤立波經(jīng)歷長時間的追趕作用過程中，自由表面波也表現(xiàn)出2個孤立波之間的追趕，兩者的波形演化在形態(tài)和時間上幾乎保持一致。

圖15 內(nèi)孤立波追趕時自由表面波波形演化Fig.15 Waveform evolution of free surface waves for overtaking

4 結(jié)論

1)依據(jù)兩內(nèi)孤立波波形是否融合成一個類孤立波波形以及是否發(fā)生相互作用等演化特征，根據(jù)初始波幅之比，可以將內(nèi)孤立波追趕分為3類。BEM-FS模型數(shù)值模擬結(jié)果和Craig[11]表面孤立波追趕分類保持一致；KdV方程Fourier譜方法離散求解下內(nèi)孤立波追趕分類和Lax[6]模擬結(jié)果分類保持一致。

2)內(nèi)孤立波追趕所有類別下的波形演化，都存在大波幅內(nèi)孤立波波幅逐漸減小，小波幅內(nèi)孤立波波幅逐漸增大，追趕演化過程中的波形波幅始終介于初始2個內(nèi)孤立波波幅之間，BEM-FS模型和KdV方程Fourier譜方法離散求解所得的幅值變化曲線存在諸多差異。

3)作用中心波幅擬合曲線存在拐點，BEM-FS模型下曲線拐點位于Craig[11]分類中類別1)和類別2)分界點附近，而KdV方程Fourier譜方法離散求解下曲線拐點位于Lax[6]分類中類別1)和類別2)分界點附近。

4)與內(nèi)孤立波直面碰撞作用過程中自由表面波的表現(xiàn)形態(tài)不同，內(nèi)孤立波追趕過程中所誘導(dǎo)的自由表面孤立波與內(nèi)孤立波的波形演化在形態(tài)和時間上幾乎保持一致。

存在于真實的海洋中的內(nèi)孤立波追趕，是否也會因為長時間相互作用導(dǎo)致完成追趕過程后波幅存在較大衰減，而呈現(xiàn)出類似于本文邊界元法FS條件BEM模型所得的模擬結(jié)果，還需相關(guān)實驗以及現(xiàn)場觀測的進一步驗證。