基于勢流理論的內孤立波追趕數值模擬

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

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

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

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

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

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

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

1.1 多域邊界元法

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

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

▽2φk=0

(1)

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

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

(2)

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

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

(3)

運動學邊界條件為:

(4)

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

(5)

(6)

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

界面S2處，根據內孤立波波形給定速度勢的初始法向導數為：

(7)

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

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

1.2 移動計算域

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

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

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

1.3 網格和時間步無關性驗證

以本文最大計算波幅a1/h1=0.4為例，進行網格和時間步無關性驗證。

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

圖3 BEM-FS模型網格無關性驗證Fig.3 The convergence verification of profile in regard to mesh sizes for BEM-FS model

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

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

2.1 理論推導

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

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

(8)

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

(9)

(10)

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

(11)

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

(12)

2.2 網格和時間步無關性驗證

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

表1 離散設置情況

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

2.3 初始條件

根據KdV方程理論解，設定波面初始條件為：

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

sech2((X+X1)/λ2)

(13)

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

3 波形演化及追趕分類

3.1 波形演化

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

圖6 內孤立波波速對比Fig.6 The comparison of wave speed of ISWs

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

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

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

圖7 BEM-FS數值模型下內孤立波波幅分別為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內孤立波追趕過程(由于內孤立波追趕作用時間長計算域大，為能在一張圖內完整清晰展現作用過程，本文所有時空分布下的波形演化過程圖均為相對作用圖，即固定大波幅內孤立波的水平位置)。此類別內孤立波追趕在作用中心時刻，會有明顯的相互作用過程。

當2個內孤立波波幅比2.941

圖8 BEM-FS數值模型下內孤立波波幅分別為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

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

圖9 BEM-FS數值模型下內孤立波波幅分別為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譜方法離散求解下內孤立波波幅分別為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譜方法離散求解下內孤立波波幅分別為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 幅值變化曲線

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

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

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

對圖12做進一步分析，顯示2種方法對于模擬內孤立波追趕存在諸多差異，具體參數對比如表2所示。在KdV方程Fourier譜方法離散求解下，內孤立波經過追趕相互作用后，波幅幾乎恢復到原始波幅，但邊界元法波幅存在一定的衰減。此外，表1顯示2種方法所得作用中心波幅a3/h1以及波幅演化時間也存在分歧。在圖12中，取內孤立波初始波幅的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時間內，2種方法所得的作用前時間t3-t2和作用后時間t4-t3分別近似相等，這說明此段時間內2種方法皆顯示追趕過程幾乎是對稱的。但在t5-t1時間內，對比作用前時間t3-t1和作用后時間t5-t3顯示，KdV方程Fourier譜方法離散求解顯示追趕過程幾乎對稱的，但BEM-FS模型顯示追趕過程呈現出明顯的不對稱性。

表2 2種方法作用參數對比

3.3 作用中心波幅變化曲線

圖13 BEM-FS數值模型下作用中心波幅a3/h1隨內孤立波初始波幅比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隨初始內孤立波波幅比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 自由表面波表現形式

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

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

4 結論

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

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

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

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

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