有限深兩層流體內(nèi)孤立波非線性數(shù)值模擬入口邊界條件研究

王 旭，尤云祥，黃文昊

(上海交通大學海洋工程國家重點實驗室，上海 200240)

內(nèi)孤立波是最大振幅發(fā)生在密度分層海水內(nèi)部的一種常見海洋動力學現(xiàn)象，不僅是海洋能量級串中的一個重要環(huán)節(jié)，也是對海洋工程結構物安全性產(chǎn)生重要影響的一類海洋環(huán)境因素。在南海流花油田的早期延長測試期間，就曾發(fā)生過因內(nèi)孤立波產(chǎn)生的突發(fā)性強流而導致纜繩拉斷、船體碰撞，甚至拉斷和擠破漂浮軟管等事故［1］。在南海陸豐油田的早期延長測試期間，也曾發(fā)生過內(nèi)孤立波產(chǎn)生的突發(fā)性強流使半潛式鉆井船與錨定油輪在連接輸油管道時發(fā)生困難等問題［2］。

非線性和頻散效應是密度分層海洋中影響內(nèi)孤立波傳播演化特性的兩個基本因素，非線性效應趨向于使波形變陡，而頻散效應則趨向于使波形變得平坦，內(nèi)孤立波的穩(wěn)定傳播是非線性與頻散效應在一定尺度上達到平衡的結果，可以用KdV、eKdV和MCC等理論模型來描述［3］。設ε=a/h為非線性參數(shù)，μ=( h /λ)2為色散參數(shù)(其中a 和λ分別為內(nèi)孤立波振幅和特征寬度，h為垂向特征尺度)，KdV理論是建立在弱非線性、弱色散且兩者平衡(即ε=O(μ)?1)條件下的一類內(nèi)孤立波理論［4］。研究表明，在小振幅情況下KdV理論與實驗結果相符，但在大振幅情況KdV理論明顯與實驗結果不符［5-8］。

隨著內(nèi)孤立波振幅的增大，高階非線性項的影響會逐漸增大，而且非線性和色散性不再保持平衡［9］，這正是KdV理論不能適用于大振幅內(nèi)孤立波的主要原因之一。這個缺陷可以通過在KdV方程中加入立方非線性項等途徑解決，所得修正方程稱為eKdV方程［10］。但eKdV理論依然是建立在弱非線性和弱色散條件下的，為克服這兩類理論中需要弱非線性這個限制性條件的缺陷，從完全非線性歐拉方程出發(fā)，Choi和Camassa［11］建立了一個強非線性和弱色散的兩層Green-Naghdi模型，在定態(tài)條件下該理論可以簡化為一個常微分方程形式的非線性方程，稱為MCC理論。這是建立在強非線性和弱色散條件下的一類定態(tài)內(nèi)孤立波理論，研究表明在大振幅內(nèi)孤立波的情況，MCC理論與實驗結果相符［12］。

上述KdV、eKdV和MCC理論均是在弱色散條件下建立的，而對KdV和eKdV理論還需要弱非線性這個限制性條件，但關于弱非線性和弱色散這兩個限制性條件迄今仍是定性的。最近，黃文昊等［13］利用大型密度分層水槽，采用雙推板內(nèi)孤立波造波方法，研究了弱非線性和弱色散這兩個條件的定量表征問題，結果表明在給定密度分層及內(nèi)孤立波振幅條件下，KdV、eKdV和MCC理論的適用性與非線性和色散參數(shù)的某種組合有關。

各種海洋環(huán)境條件下深海浮式平臺載荷特性及其運動響應的合理確定，是其設計與應用中的一項關鍵性指標。程友良等［14］和蔡樹群等［15-16］將Morison公式與KdV理論結合，而XIE等［17］則將Morison公式與MCC理論結合，研究了內(nèi)孤立波作用下直立貫底圓柱體的載荷特性問題。最近，尤云祥等［18-19］將Morison公式與eKdV理論結合，研究了內(nèi)孤立波作用下張力腿和半潛式平臺的載荷與動力響應問題，而宋志軍等［20］將Morison公式與KdV理論結合，研究了內(nèi)孤立波作用下Spar平臺的載荷與動力響應問題。在這些研究中，Morison公式中的慣性力和拖曳力系數(shù)是參照表面波的方法選取的，這種選取方法尚缺乏理論和實驗依據(jù)。

隨著計算流體力學(CFD)技術的發(fā)展，采用CFD方法來研究內(nèi)孤立波的生成傳播及其與海上結構物的相互作用問題，已成為目前國際上的研究熱點之一。由于采用CFD方法可以直接獲得內(nèi)孤立波與海洋結構物相互作用過程中的速度場及壓力場的變化特性，因此可以直接獲得內(nèi)孤立波作用下海上結構物的載荷等水動力特性。在這種方法中，通常在入口邊界條件中以定態(tài)內(nèi)孤立波的解析解(包括KdV、eKdV和MCC等)來驅動完全非線性模型，進而模擬內(nèi)孤立波的生成傳播及其與海上結構物的作用特性。

陳鈺等［21］、關輝等［22］和李水娟等［23］利用 KdV 理論解作為入口邊界條件，而高原雪等［24］利用 MCC 理論解作為入口邊界條件，采用CFD方法研究了內(nèi)孤立波的生成傳播問題。付東明等［25］利用KdV理論解作為入口邊界條件，而陳杰等［26］和關輝等［27］利用eKdV理論解作為入口邊界條件，采用CFD方法研究了內(nèi)孤立波作用下水下潛體的載荷特性問題。劉碧濤等［28］利用eKdV理論解作為入口邊界條件，采用CFD方法研究了內(nèi)孤立波作用下深海立管的載荷特性問題。

根據(jù)文獻［13］的實驗結果，在給定密度分層及內(nèi)孤立波振幅條件下，由于KdV、eKdV和MCC理論的適用性與非線性和色散參數(shù)的某種組合有關，因此如果選擇不合適的內(nèi)孤立波理論解作為CFD數(shù)值模擬的入口邊界條件，則有可能會導致振幅及波形都不可控的數(shù)值模擬結果。因此，在給定密度分層及內(nèi)孤立波振幅條件下，如何選擇合適的內(nèi)孤立波理論解作為CFD數(shù)值模擬的入口邊界條件是利用CFD方法來研究內(nèi)孤立波的生成傳播及其與海上結構物相互作用問題時需要研究和解決的關鍵之一。

有鑒于此，以有限深兩層流體中內(nèi)孤立波誘導上下層深度平均水平速度為入口邊界條件，以理想流體完全非線性歐拉方程為控制方程，結合VOF動態(tài)界面跟蹤方法，對內(nèi)孤立波的生成特性進行完全非線性數(shù)值模擬。在此基礎上，研究在給定密度分層及內(nèi)孤立波振幅條件下如何選擇合適的內(nèi)孤立波理論解作為完全非線性數(shù)值模擬入口邊界條件等問題。

1 數(shù)值模擬方法

設兩層流體均為理想不可壓縮流體，當流體處于靜平衡狀態(tài)時，上層流體深度與密度分別為h1和ρ1，而下層流體深度與密度分別為h2和ρ2，水深為h=h1+h2。建立直角坐標系oxz，其中ox軸位于流體靜止時兩層流體的界面上，oz軸垂直向上為正。對水表面作剛蓋假定，假定水底為不可滲透平坦的剛性固壁。設ζ為兩層流體界面位移，ui、wi為流場中的速度矢量，其中i=1、2分別表示上下層流體。

在上下層流體中，流場控制方程如下:

其中，g為重力加速度，各物理量的下標t、x和z分別表示關于時間和空間坐標的導數(shù)，pi為壓力。

在兩層流體界面上，要求滿足如下連續(xù)性條件:

在水面和水底，要求滿足如下固壁條件:

流場計算的控制區(qū)域如圖1所示，包括內(nèi)孤立波生成傳播區(qū)和消波區(qū)兩個區(qū)域。采用速度入口方法生成內(nèi)孤立波，當形成穩(wěn)定的內(nèi)孤立波后，對所生成內(nèi)孤立波的傳播特性進行監(jiān)測分析。

圖1 流場計算控制區(qū)域Fig.1 Governing domain for flow computation

速度入口方法為:在上層流體入口處，入口速度取為u1;在下層流體入口處，入口速度取為。內(nèi)孤立波在生成與傳播過程中，兩層流體的界面會發(fā)生變化，采用VOF(volume of fluid)方法追蹤兩層流體界面的變化［29］。在VOF方法中，用體積分數(shù)aq表示單元內(nèi)第q相流體所占體積與該單元的體積之比，其中q=1表示上層流體，q=2表示下層流體。當aq=1時，表示單元內(nèi)全部為第q相流體，而當aq=0時，表示單元不存在第q相流體，對含界面的單元則有0＜aq＜1。在每一個流體單元中各相的體積分數(shù)之和為1，即=1，而且對于任一相的體積分數(shù)aq，要求滿足如下連續(xù)方程:

采用海綿層方法對水槽末端的內(nèi)孤立波進行消波處理，該方法通過添加人工粘性項實現(xiàn)阻尼消波之目的。為此，對消波區(qū)，在動量方程的右端添加源項，將其修改為

其中，海綿層衰減系數(shù)μ依據(jù)文獻［30］中的方法確定。

采用KdV、eKdV和MCC理論計算內(nèi)孤立波誘導上下層流體中的層深度平均水平速度。為此，設a和λ分別為內(nèi)孤立波振幅與特征寬度，定義非線性參數(shù)ε=a/h和色散參數(shù)μ=(h/λ)2，那么在弱色散和弱非線性(即ε?1而且μ?1)條件下，界面位移ζ滿足如下的eKdV方程［3］

其中

當內(nèi)孤立波為定態(tài)平面前進波時，eKdV方程(10)有如下形式的解［3］:

其中

其中，cekdv為內(nèi)孤立波eKdV理論解的相速度。

當ε=O(μ)?1，即當內(nèi)孤立波是弱非線性、弱色散且兩者平衡時，在eKdV方程中的立方非線性項可以忽略，即可取c3=0，這時eKdV理論解(13)即退化為如下的KdV理論解［3］

其中，ckdv為內(nèi)孤立波KdV理論解的相速度。

在上述理論中，均要求內(nèi)孤立波是弱非線性和弱色散的，對KdV理論還要求兩者是平衡的。因此，這兩類內(nèi)孤立波理論一般只適用于振幅較小的情況。為克服這個缺陷，Choi和Camassa［11］建立了一個強非線性和弱色散的內(nèi)孤立波理論模型，稱為MCC理論模型，其定態(tài)解可表示為

其中，cMCC為內(nèi)孤立波MCC理論解的相速度，

a-和a+(a-＜a+)分別為下面方程的兩個根

式中:

對于KdV理論解來說，只要h1/h2≠ ( h1/h2)c，則不存在極限振幅，即無論振幅有多大，該理論均能給出對應的內(nèi)孤立波解。Camassa等［12］指出，KdV理論無極限振幅的這個缺陷，正是其不能適用于大振幅內(nèi)孤立波的原因之一。與KdV理論的一個重要不同之處是，其它兩類內(nèi)孤立波理論均存在極限振幅，其中eKdV和MCC理論解的極限振幅分別為［3］

2 結果與分析

文中將結合該文獻中的相關實驗結果進行數(shù)值模擬與分析。為此，數(shù)值水槽主尺度、上下層流體密度及其深度比均與該文一致。在數(shù)值模擬中，內(nèi)孤立波生成傳播區(qū)的長度為24 m，而消波區(qū)的長度為6 m。生成傳播區(qū)內(nèi)網(wǎng)格的縱向間距為0.03 m，垂向網(wǎng)格劃分以底部向上0.4 m位置為界，上部區(qū)域取網(wǎng)格垂向間距為0.005 m，下部區(qū)域首層垂向間距取為0.005 m，為保證網(wǎng)格平滑過渡，后續(xù)網(wǎng)格間距按照1.02的比例逐漸加寬;消波區(qū)內(nèi)垂向網(wǎng)格間距同生成傳播區(qū)的情況一致，而縱向首層網(wǎng)格間距取0.03 m，后續(xù)網(wǎng)格間距按1.04的比例逐漸加寬;整個計算域的網(wǎng)格單元共計60 800個。

利用商業(yè)軟件Fluent進行數(shù)值模擬，初始時水槽中流場物理量均取為零，時間步長為0.01 s。計算區(qū)域采用結構化網(wǎng)格進行離散，在兩層流體界面附近，網(wǎng)格需進行加密處理。控制方程采用有限體積法(FVM)進行離散，對流項采用二階迎風格式離散，時間步推進選用二階隱式時間匹配格式，兩相界面的構造方法選用幾何重構法，應用非定常壓力的隱式算子分割法(PISO)進行壓力-速度的耦合計算。

文獻［13］的實驗結果表明，存在一個臨界色散參數(shù)μ0=0.1，以此為基準對弱色散條件μ?1進行定量化，可將其分為弱色散(μ＜μ0)和強色散(μ≥μ0)兩種情況;在弱色散(μ＜μ0)條件下，對非線性條件也進行定量化，可將其分為弱非線性(ε≤μ)、中等非線性(μ＜ε≤)和強非線性(ε＞)三種情況;在系列實驗基礎上，獲得了KdV、eKdV和MCC理論的如下適用性條件:KdV理論適用于弱非線性和弱色散情況，eKdV理論適用于中等非線性和弱色散情況，而MCC理論適用于強非線性和弱色散或強色散情況。

在三種上下層流體深度比情況下，文獻［13］中系列實驗所得內(nèi)孤立波振幅如表1所示。在本文數(shù)值模擬中，初始內(nèi)孤立波的設計振幅ad也取為表中所列相應的值。由式(21)可知，當h1∶h2=3∶7和1∶9時，h=-0.182 6和 -0.211 8，而/h=-0.196 7和 -0.396 7;當 h1∶h2=2∶8時，/h =-0.234 1，而/h=-0.296 7。由此可知，當h1∶h2=3∶7和1∶9時，各設計振幅均在eKdV和MCC理論解的極限振幅內(nèi);而當h1∶h2=2∶8時，各設計振幅雖均在MCC理論解的極限振幅內(nèi)，但最后三個振幅已超出eKdV理論解的極限振幅。

表1 數(shù)值模擬中的無因次設計振幅ad/hTab.1 Non-dimensional design amplitudes for numerical simulations ad/h

在數(shù)值模擬中，對給定的設計振幅ad，利用KdV、eKdV和MCC理論解分別計算其誘導的上下層流體中的層深度平均水平速度(見式(6))，將其作為入口邊界條件來驅動完全非線性數(shù)值模型(見式(1)-(5))。在上述初始內(nèi)孤立波設計振幅ad下，這里重點研究如何選擇合適的定態(tài)內(nèi)孤立波理論解來驅動完全非線性數(shù)值模型的問題。

為此，設εm和μm分別為數(shù)值模擬所得波的非線性和色散參數(shù)，而εd和μd分別為在設計振幅ad下利用內(nèi)孤立波理論解計算得到的非線性和色散參數(shù)。在圖2中，給出了當h1/h2=3∶7時，(εm，μm)和(εd，μd)的結果比較。其中:“▲”表示(εm，μm)，“2”表示(εd，μd)。在給定設計振幅ad下，點▲和2之間的橫向距離表示μm和μd之間的絕對誤差，而點▲和2之間的垂向距離表示εm和εd之間的絕對誤差。由圖可知，在給定設計振幅ad下，利用eKdV理論解作為入口邊界條件，所得μm和μd以及εm和εd之間的絕對誤差最小。因此，在表1所示各設計振幅下，利用eKdV理論解作為入口邊界條件，數(shù)值模擬所得內(nèi)孤立波振幅及設計振幅誤差最小，而且兩者波形最為吻合。

圖2 當 h1∶h2=3∶7 時 (εm，μm)和 (εd，μd)結果比較Fig.2 Comparisons between(εm，μm)and(εd，μd)when h1∶h2=3∶7

由圖2(a)可知，在各設計振幅ad下，KdV理論解和數(shù)值模擬所得波的非線性和色散參數(shù)均屬于弱色散和中等非線性區(qū)域，而文獻［13］的實驗表明，KdV理論適用于弱色散和弱非線性的情況，因此采用KdV理論解作為入口邊界條件本身就是不合適的。由圖2(b)可知，在各設計振幅ad下，eKdV理論解和數(shù)值模擬所得波的非線性和色散參數(shù)也均屬于弱色散和中等非線性區(qū)域，而文獻［13］的實驗表明，eKdV理論適用于弱色散和中等非線性的情況，兩者正好是一致的，即采用eKdV理論解作為入口邊界條件是合適的。由圖2(c)可知，MCC理論解和數(shù)值模擬所得波的非線性和色散參數(shù)均屬于弱色散和中等非線性區(qū)域，而文獻［13］的實驗表明，MCC理論適用于弱色散和強非線性或強色散的情況，因此采用MCC理論解作為入口邊界條件本身就是不合適的。

在圖3中，給出了當h1/h2=3∶7時，數(shù)值模擬結果與初始內(nèi)孤立波及相關實驗結果的比較。由圖中第一列可知，當利用KdV理論解作為入口邊界條件時，數(shù)值模擬所得波的振幅及其波形與初始輸入波的振幅及其波形差異都很大，而且隨著設計振幅的增大，數(shù)值模擬結果還會出現(xiàn)顯著的尾波現(xiàn)象。由圖中第二列可知，當利用eKdV理論解作為入口邊界條件時，數(shù)值模擬所得波的振幅及其波形與初始輸入波和文獻［13］中實驗所得波振幅及其波形之間的差異都很小，而且數(shù)值模擬結果沒有明顯的尾波現(xiàn)象。由圖中第三列可知，與利用KdV理論解作為入口邊界條件的結果相比，當利用MCC理論解作為入口邊界條件時，數(shù)值模擬所得波的振幅及其波形與初始輸入波的振幅及其波形之間的差異明顯減小，但這種差異要比利用eKdV理論解作為入口邊界條件的結果大，因此采用eKdV理論解作為入口邊界條件更為合適。

由圖2可知，當h1/h2=3∶7時，在各設計振幅ad下，利用KdV理論解作為入口邊界條件，數(shù)值模擬所得波均屬于弱色散和中等非線性的。因此，根據(jù)文獻［13］中的實驗結果，如果該波為內(nèi)孤立波，則應該屬于eKdV型的。為此，在圖4中，給出了當ad/h=0.034 0、0.105 5和0.158 0時，利用KdV理論解作為入口邊界條件所得數(shù)值模擬結果與等振幅eKdV理論解的比較。由圖4(a)可知，在較小無因次振幅的情況，數(shù)值模擬結果的波形與等振幅eKdV理論解波形一致，這表明數(shù)值模擬所得波為eKdV型內(nèi)孤立波。由圖4(b)和4(c)可知，在較大無因次振幅的情況，數(shù)值模擬所得波的波形與等振幅eKdV理論解的波形差異很大，這表明在這種情況下利用KdV理論解作為入口邊界條件是不合適的。

圖3 當h1∶h2=3∶7時數(shù)值模擬結果與初始內(nèi)孤立波及相關實驗結果比較Fig.3 Comparisons of numerical results with theoretical and experimental ones when h1∶h2=3∶7

圖4 當h1∶h2=3∶7時利用KdV理論解作為入口邊界條件所得數(shù)值模擬結果與eKdV理論解比較Fig.4 Comparisons of numerical results with eKdV theoretical solutions based on the KdV input conditions when h1∶h2=3∶7

圖5給出了當h1/h2=2∶8時，(εm，μm)和(εd，μd)的結果比較。由圖可知，對表1中前3個設計振幅ad，利用KdV理論解作為入口邊界條件，所得μm和μd以及εm和εd之間的誤差最小，而且利用KdV理論解和數(shù)值模擬所得波的非線性和色散參數(shù)均屬于弱色散和弱非線性區(qū)，與文獻［13］實驗結果一致，即這時采用KdV理論解作為入口邊界條件最為合適。對表1中第4到9個設計振幅ad，利用eKdV理論解作為入口邊界條件，所得μm和μd以及εm和εd之間的誤差最小，而且利用eKdV理論解和數(shù)值模擬所得波的非線性和色散參數(shù)均屬于弱色散和中等非線性區(qū)，與文獻［13］實驗結果一致，即這時采用eKdV理論解作為入口邊界條件最為合適。對表1中最后4個設計振幅ad，利用MCC理論解作為入口邊界條件，所得μm和μd以及εm和εd之間的誤差最小，而且利用MCC理論解和數(shù)值模擬所得波的非線性和色散參數(shù)均屬于弱色散和強非線性區(qū)域，與文獻［13］實驗結果一致，即這時采用MCC理論解作為入口邊界條件最為合適。

圖5 當 h1∶h2=2∶8 時 (εm，μm)和 (εd，μd)結果比較Fig.5 Comparisons between(εm，μm)and(εd，μd)when h1∶h2=2∶8

在圖6中，給出了當h1/h2=2∶8時，數(shù)值模擬結果與初始內(nèi)孤立波及相關實驗結果的比較。由圖可知，當ad/h=0.062 1時，與利用eKdV和MCC理論解作為入口邊界條件的結果相比，利用KdV理論解作為入口邊界條件所得數(shù)值模擬波的振幅及其波形與初始輸入波最為吻合，而且與文獻［13］中的實驗一致。當ad/h=0.130 4時，與利用KdV和MCC理論解作為入口邊界條件的結果相比，利用eKdV理論解作為入口邊界條件所得數(shù)值模擬波的振幅及其波形與初始輸入波最為吻合，而且與文獻［13］中的實驗一致。當ad/h=0.243 3時，與利用KdV和eKdV理論解作為入口邊界條件的結果相比，利用MCC理論解作為入口邊界條件所得數(shù)值模擬波的振幅及其波形與初始輸入波最為吻合，而且與文獻［13］中的實驗一致。

圖6 當h1∶h2=2∶8時數(shù)值模擬結果與初始內(nèi)孤立波及相關實驗結果比較Fig.6 Comparisons of numerical results with theoretical and experimental ones when h1∶h2=2∶8

由圖5可知，當h1/h2=2∶8時，對表1中最后10個設計振幅ad，利用KdV理論解作為入口邊界條件，數(shù)值模擬所得波均屬于弱色散和中等非線性的。因此，根據(jù)文獻［13］中實驗結果，如果該波為內(nèi)孤立波，則應該屬于eKdV型的。為此，在圖7中，給出了當ad/h=0.130 4和0.243 3時，利用KdV理論解作為入口邊界條件所得數(shù)值模擬結果與等振幅eKdV理論解的比較。由圖可知，數(shù)值模擬結果的波形與等振幅eKdV理論解波形一致，這表明數(shù)值模擬所得波為eKdV型內(nèi)孤立波。

圖7 當h1∶h2=2∶8時，利用KdV理論解作為入口邊界條件所得數(shù)值模擬結果與eKdV理論解比較Fig.7 Comparisons of numerical results with eKdV theoretical solutions based on the KdV input conditions when h1∶h2=2∶8

圖8 當 h1∶h2=1∶9 時 (εm，μm)和 (εd，μd)結果比較Fig.8 Comparisons between(εm，μm)and(εd，μd)when h1∶h2=1∶9

在圖8中，給出了當h1∶h2=1∶9時，(εm，μm)和 (εd，μd)的結果比較。由圖可知，在給定設計振幅 ad下，利用MCC理論解作為入口邊界條件，所得μm和μd以及εm和εd之間的誤差最小。因此，在表1所示各設計振幅下，利用MCC理論解作為入口邊界條件，數(shù)值模擬所得內(nèi)孤立波振幅與設計振幅誤差最小，而且兩者波形最為吻合。由圖8(a)和8(b)可知，在各設計振幅ad下，KdV和eKdV理論解和數(shù)值模擬所得波的非線性和色散參數(shù)均屬于強弱色散區(qū)域。文獻［13］的實驗表明，KdV理論適用于弱色散和弱非線性的情況，而eKdV理論適用于弱色散和中等非線性的情況，因此采用KdV和eKdV理論解作為入口邊界條件本身就是不合適的。由圖8(c)可知，在各設計振幅ad下，利用MCC理論解和數(shù)值模擬所得波的非線性和色散參數(shù)均屬于強弱色散區(qū)域，而文獻［13］的實驗表明，MCC理論適用于強弱色散的情況，兩者正好是一致的，即采用MCC理論解作為入口邊界條件是合適的。

在圖9中，給出了當h1∶h2=1∶9時，數(shù)值模擬結果與初始內(nèi)孤立波及相關實驗結果的比較。由圖中第一列和第二列可知，當利用KdV和eKdV理論解作為入口邊界條件時，隨著設計振幅的增大，數(shù)值模擬所得波的振幅及其波形與初始輸入波的振幅及其波形之間的差異也隨著增大。由圖中第三列可知，當利用MCC理論解作為入口邊界條件時，數(shù)值模擬所得波的振幅及其波形與初始輸入波和文獻［13］中實驗所得波振幅及其波形之間的差異都很小。

圖9 當h1∶h2=1∶9時數(shù)值模擬結果與初始內(nèi)孤立波及相關實驗結果比較Fig.9 Comparisons of numerical results with theoretical and experimental ones when h1∶h2=1∶9

圖10 當h1∶h2=1∶9時利用KdV和eKdV理論解作為入口邊界條件所得數(shù)值模擬結果與MCC理論解比較Fig.10 Comparisons of numerical results with MCC theoretical solutions based on the KdV and eKdV input conditions when h1∶h2=1∶9

由圖8可知，當h1/h2=1∶9時，在各設計振幅ad下，利用KdV和eKdV理論解作為入口邊界條件，數(shù)值模擬所得波均屬于強色散的。因此，根據(jù)文獻［13］中的實驗結果，如果該波為內(nèi)孤立波，則應該屬于MCC型的。為此，在圖10中，給出了當ad/h=0.041 8，0.137 3和0.151 3時，利用KdV和eKdV理論解作為入口邊界條件所得數(shù)值模擬結果與等振幅MCC理論解的比較。由圖可知，數(shù)值模擬結果的波形與等振幅MCC理論解波形一致，這表明數(shù)值模擬所得波為MCC型內(nèi)孤立波。

3 結語

這里為有限深水內(nèi)孤立波非線性數(shù)值模擬中如何選擇合適的理論解作為入口邊界條件提供了理論依據(jù)。所得入口邊界條件選擇方法僅限于有限深水的情況，針對無限深水條件下應如何選擇合適的內(nèi)孤立波理論解作為入口邊界條件等問題，尚需做進一步研究。

［1］ 陳景輝.南海流花11-1深海油田開發(fā)工程［J］.中國海洋平臺，1996，11(1):43-45.(CHEN Jing-hui.Liuhua 11-1 deep sea oil-field development project in South China Sea［J］.China Offshore Platform，1996，11(1):43-45.(in Chinese))

［2］ Bole J B，Ebbesmeyer C C，Romea R D.Soliton currents in the South China Sea:measurements and theoretical modeling［C］//The 16thOffshore Technology Conference.1994:367-376.

［3］ Helfrich K R，Melville W K.Long nonlinear internal waves［J］.Ann.Rev.Fluid Mech.，2006，38:395-425.

［4］ Ostrovsky L A，Stepanyants Yu A.Do internal solitons exist in the ocean［J］.Rev.Geophys.，1989，27(3):293-310.

［5］ Michallet H，Barthelemy E.Experimental study of interfacial solitary waves［J］.J.Fluid Mech.，1998，366:159-177.

［6］ Grue J，Jensen A，Rusas P E，et al.Properties of large-amplitude internal waves［J］.J.Fluid Mech.，1999，380:257-278.

［7］ Sveen J K，Guo Y，Davies P A，et al.On the breaking of internal solitary waves at a ridge［J］.J.Fluid Mech.，2002，469:161-188.

［8］ Walker S A，Martin A J，Easson W J，et al.Comparison of laboratory and theoretical internal solitary wave kinematics［J］.Journal of Waterway，Port，Coastal and Ocean Engineering，2003，129(5):210-218.

［9］ Brandt P，Rubino A，Alpers W，et al.Internal waves in a strait of messina studied by a numerical model and synthetic aperture radar form the ERS1/2 satellites［J］.J.Phys.Oceanography.，1997，27(5):648-663.

［10］Djordjevic V D.The fission and disintegration of internal solitary waves moving over two-dimensional topography［J］.J.Phys.Oceanography，1978，8(6):1016-1024.

［11］Choi W，Camassa R.Fully nonlinear internal waves in a two-fluid system［J］.J.Fluid Mech.，1999，396:1-36.

［12］Camassa R，Choi W，Michallet H，et al.On the realm of validity of strongly nonlinear asymptotic approximations for internal waves［J］.J.Fluid Mech.，2006，549:1-23.

［13］黃文昊，尤云祥，王旭，等.有限深兩層流體中內(nèi)孤立波造波實驗及其理論模型［J］.物理學報，2013，62(8):1-14.(HUANG Wen-hao，YOU Yun-xiang，WANG Xu，et al.Wave-making experiments and theoretical models for internal solitary waves in a two-layer fluid of finite depth［J］.Acta Phys.Sin.，2013，62(8):1-14.(in Chinese))

［14］CHENG Y L，LI J C，LIU Y F.The induced flow field by internal solitary wave and its action on cylindrical piples in the stratified ocean［C］//The 4th International Conference on Fluid Mechanics.2004:296-299.

［15］CAI S Q，LONG X M，GAN Z J.A method to estimate the forces exerted by internal solitons on cylinder piles［J］.Ocean Engineering，2003，30(5):673-689.

［16］CAI S Q，LONG X M，WANG S G.Force and torques exerted by internal solitions in shear flows on cylinder piles［J］.Applied Ocean Research，2008，30(1):72-77.

［17］XIE J S，JIANG Y J，YANG L G.Strongly nonlinear internal solution load on a small vertical circular cylinder in two-layer fluids［J］.Applied Mathematical Modelling，2010，34(8):2089-2101.

［18］尤云祥，李巍，何景異，等.海洋內(nèi)孤立波中張力腿平臺的水動力特性［J］.上海交通大學學報，2010，44:12-17.(YOU Yun-xiang，LI Wei，HE Jing-yi，et al.Hydrodynamic characteristics of tension leg platforms in ocean internal solitary waves［J］.Journal of Shanghai Jiao Tong University，2010，44:12-17.(in Chinese))

［19］尤云祥，李巍，胡天群，等.內(nèi)孤立波中半潛平臺動力響應特性［J］.海洋工程，2012，30(2):1-7.(YOU Yun-xiang，LI Wei，HU Tian-qun，et al.Dynamic responses of a semi-submersible platform in internal solitary waves［J］.The Ocean Engineering，2012，30(2):1-7.(in Chinese))

［20］宋志軍，勾瑩，滕斌，等.內(nèi)孤立波作用下Spar平臺的運動響應［J］.海洋學報，2010，32:12-19.(SONG Zhi-jun，GOU Ying，TENG Bin，et al.The motion responses of a Spar platform under internal solitary wave［J］.Acta Oceanologica Sinica，2010，32(2):12-19.(in Chinese))

［21］陳鈺，朱良生.基于 FLUENT的海洋內(nèi)孤立波數(shù)值水槽模擬［J］.海洋技術，2009，28(4):72-75.(CHEN Yu，ZHU Liang-sheng.Numerical wave tank simulation of oceanic internal solitary waves based on FLUENT［J］.Ocean Technology，2009，28(4):72-75.(in Chinese))

［22］關暉，蘇曉冰，田俊杰.分層流體中內(nèi)孤立波數(shù)值造波方法及其應用［J］.力學季刊，2011，32(2):218-223.(GUAN Hui，SU Xiao-bing，TIAN Jun-jie.Numerical method of solitary wave generation in a stratified flow and its applications［J］.Chinese Quarterly Of Mechanics，2011，32(2):218-223.(in Chinese))

［23］李水娟，余真真，羅龍洪，等.三維數(shù)值波浪水槽內(nèi)孤立波特性分析［J］.水電能源科學，2012，30(7):96-99.(LI Shui-juan，YU Zhen-zhen，LUO Long-hong，et al.Analysis of internal solitary waves characteristics of 3d numerical wave flume［J］.Water Resources and Power，2012，30(7):96-99.(in Chinese))

［24］高原雪，尤云祥，王旭，等.基于 MCC理論的內(nèi)孤立波數(shù)值模擬［J］.海洋工程，2012，30(4):29-36.(GAO Yuan-xue，YOU Yun-xiang，WANG Xu，et al.Numerical simulation for the internal solitary wave based on MCC theory［J］.The Ocean Engineering，2012，30(4):29-36.(in Chinese))

［25］付東明，尤云祥，李巍.兩層流體中內(nèi)孤立波與潛體相互作用數(shù)值模擬［J］.海洋工程，2009，27(3):38-44.(FU Dong-ming，YOU Yun-xiang，LI Wei.Numerical simulation of internal solitary waves with submerged body in a two-layer fluid［J］.The Ocean Engineering，2009，27(3):38-44.(in Chinese))

［26］陳杰，尤云祥，劉曉東，等.內(nèi)孤立波與有航速潛體相互作用數(shù)值模擬［J］.水動力學研究與進展，2010，25(3):344-350.(CHEN Jie，YOU Yun-xiang，LIU Xiao-dong，et al.Numerical simulation of interaction of internal solitary waves with a moving submarine［J］.Chinese Journal of Hydrodynamics，2010，25(3):344-350.(in Chinese))

［27］關暉，魏崗，杜輝.內(nèi)孤立波與潛艇相互作用的水動力學特性［J］.解放軍理工大學學報:自然科學版，2012，13(5):577-582.(GUAN Hui，WEI Gang，DU Hui，et al.Hydrodynamic properties of interactions of three-dimensional internal solitary waves with submarine［J］.Journal of PLA University of Science and Technology，Natural Science Edition，2012，13(5):577-582.(in Chinese))

［28］劉碧濤，李巍，尤云祥，等.內(nèi)孤立波與深海立管相互作用數(shù)值模擬［J］.海洋工程，2011，29(4):1-7.(LIU Bi-tao，LI Wei，YOU Yun-xiang，et al.Numerical simulation of interaction of internal solitary waves with deep-sea risers［J］.The Ocean Engineering，2011，29(4):1-7.(in Chinese))

［29］Hirt C W，Nichols B D.Volume of fluid(VOF)method for the dynamics of free boundaries［J］.J.Computational Physics，1981，39(1):201-225.

［30］韓鵬.基于VOF方法的不規(guī)則波阻尼消波研究［D］.大連:大連理工大學，2008.(HAN Peng.The study of damping absorber for irregular waves based on VOF method［D］.Dalian:Dalian University of Technology，2008.(in Chinese))