波浪與帶有窄縫結構作用的非線性數值模擬

寧德志，蘇曉杰，滕 斌

（大連理工大學 海岸和近海工程國家重點實驗室，遼寧 大連 116024）

波浪與帶有窄縫結構作用的非線性數值模擬

寧德志，蘇曉杰，滕 斌

（大連理工大學 海岸和近海工程國家重點實驗室，遼寧 大連 116024）

針對波浪與帶有窄縫結構作用產生的流體共振問題，采用基于域內源造波技術的時域高階邊界元方法并在窄縫內流體引入人工阻尼，建立了自由水面滿足完全非線性邊界條件的二維時域數值波浪水槽模型。求解中采用混合歐拉-拉格朗日方法追蹤流體瞬時水面，運用四階龍格庫塔方法更新下一時間步的波面和速度勢，利用加速度勢的方法來求得作用結構上的瞬時波浪荷載。通過與已發表文獻的數值與實驗結果對比，驗證了模型的準確性。同時通過大量的數值計算研究了共振條件入射波浪非線性對反射波高、透射波高和窄縫內波高、及結構波浪荷載和壓力分布的影響。

窄縫；數值波浪水槽；流體共振；源造波技術；高階邊界元

0 引 言

在海洋工程結構作業過程中，經常會出現多浮體結構并排布置聯合作業的情況，在浮體間會出現相對結構尺度小很多的窄縫。如由很多獨立模塊組成的超大型浮體，這些模塊間存在窄縫；浮式生產儲卸油系統（FPSO）與穿梭油輪并聯進行油氣外輸時和大型船舶在碼頭前系泊時也都存在這種窄縫問題。窄縫內的水動力相互作用是十分復雜的水動力干涉問題，窄縫內水體在某些頻率波浪作用下會發生共振現象，誘發很大的波浪爬高和荷載，進而對結構的作業安全帶來很大的負面影響。

關于多結構間窄縫內水體共振問題，國內外學者已開展了很多研究探索其產生機理及相關的水動力相互作用問題等。譬如，Miao等[1]采用漸近匹配法研究了帶狹縫二維雙箱的共振現象，指出了共振頻率與方箱的吃水深度和狹縫寬度的關系，并給出狹縫很小時雙箱的理論共振頻率。滕斌等[2]和何廣華等[3]采用比例邊界有限元方法分別研究了波浪與兩箱和三箱結構作用下窄縫內水體共振現象。Saitoh等（2006）[4]和 Iwataet等（2007）[5]分別對不同入射波浪作用下兩個方箱和三個方箱之間窄縫的波高變化進行了試驗研究，發現窄縫內最大共振波高可以達到入射波高的五倍，共振頻率與共振波高與箱體的吃水深度、窄縫寬度和箱體個數成一定的函數關系。Kristiansen和Faltinsen（2009）[6]對波浪與具有窄縫的二維箱體和固定式岸壁結構之間的相互作用開展了數值和試驗模擬，發現線性和非線性勢流數值結果都比共振條件下的窄縫內試驗波高要大，尤其是線性結果更大。Kristiansen和Faltinsen（2010）[7]對波浪與具有窄縫的二維浮式箱體和固定岸壁式結構間的相互作用進行了數值和試驗模擬，研究了共振條件下浮體三個自由度的運動響應與水動力系數。Lu等（2010，2011）[8-9]采用粘性流模型和改進的線性勢流模型對結構間窄縫內水體共振問題進行了研究，發現在選取一定的人工阻尼系數條件下，線性勢流模型也可以得到與試驗和粘性流模型一樣的結果。Li和Zhang（2016）[10]基于完全非線性勢流理論采用邊界元方法對兩個做升沉運動方箱間窄縫內水體運動進行了模擬分析。Feng和Bai（2015）[11]采用基于入散射分離技術的完全非線性高階邊界元方法對兩個三維方箱間窄縫內容水體共振運動進行了模擬。

在上述研究工作基礎上，本文通過采用高階邊界元方法建立自由水面滿足完全非線性邊界條件的時域數值波浪水槽，采用域內源造波方法造波并布置前置阻尼層消除二次反射影響，進而可以在較短計算域內實現長時間模擬；在窄縫水體內引入常人工阻尼系數來近似共振條件下粘性耗散，進而求解波浪與具有窄縫兩結構的相互作用問題。通過與Saitoh等[4]實驗結果及Lu等[9]的CFD數值結果進行對比證明了本模型的準確性，并進一步通過大量數值計算研究了入射波浪非線性對反射波高、透射波高、窄縫內共振波高和波浪荷載等的影響，以及共振條件下結構各側面的壓力分布等。

1 數學模型

本文研究的波浪與具有窄縫的兩方箱作用布置如圖1所示，建立的笛卡爾坐標系Oxz，z=0位于靜水面上，且z軸向上為正，x軸右方向為正。計算域包含自由水面Гf和固體邊界ГN（包括水底Гd和箱體Гb）。波浪由控制垂直源造波面（Гs）的流量密度產生。圖中h為水槽靜水深，W為箱體寬度，D為箱體吃水深度，Wg為兩箱體間窄縫的寬度，點I、II、III和IV分別代表箱體A和箱體B迎浪點和背浪點，而點Ⅴ是窄縫中間位置水面點。b、c、d、b0、c0、d0分別代表箱體A和箱體B的迎浪面、背浪面和底面。在流體無粘、不可壓縮和流動無旋的假定下，整個流域的速度可用速度勢的梯度來描述，上述問題的控制方程為由速度勢滿足的泊松（Poisson）方程[12-13]，即

圖1 水槽示意圖Fig.1 Definition sketch of the wave flume

式中：q*（xs,z, t）=2vδ（x- xs）為造波源強度，造波位置x=xs（本文均取xs=0），v為流體質點水平速度，本文給定二階Stokes速度解析解。

在自由水面上，滿足完全非線性動力學和運動學邊界條件，本文中采用混合歐拉-拉格朗日方法更新自由水面，利用物質導數d/d t=?/?t+v·▽，并在計算域的上游和下游區域的自由水面分別布置人工阻尼層來吸收從結構物反射回來的波浪與出流波浪[14]，在窄縫內布置一常參數人工阻尼來近似由于渦和分流引起的粘性耗散，Kim（2003）[15]用此方法模擬了自振頻率下容器內液體晃蕩的粘性耗散。自由水面邊界條件可以寫成以下形式：

式中：η代表自由水面的鉛垂位移，g是重力加速度，X0=（x0, ）0 是指水質點靜止時的位置，造波源位置為x=0，阻尼系數μ1用于計算域邊界的兩個阻尼層，阻尼系數μ2用于窄縫內自由水面，其數值根據試驗粘性耗散來確定，L為阻尼層長度，本文取1.5倍波長，即1.5λ，k是波數，ω是波浪角頻率，滿足如下線性色散方程關系

在固定的結構表面和底面邊界上，流體法向速度為0，滿足固壁不可滲透邊界條件，即：

假定自由水面初始時是靜止的，即：

在整個流域內對速度勢應用格林第二定理，可得到如下邊界積分方程[14]：

式中：p=（x0,z0）為源點，q=（x, z）為場點，C為固角系數，Ω代表整個流域，G是簡單格林函數，考慮到水底鏡像，可以表示為如下形式：

式中：r1為p和q兩點距離，r2為p和q關于水底鏡像之間距離。

本文用三節點高階邊界元離散計算域成一些曲線單元，單元內任一點的幾何坐標和速度勢等物理量可以用二次形狀函數插值得到。積分方程經高階邊界元離散后，可通過求解線性方程組得到未知量。計算中認為當前時刻物面上的速度勢法向導數和自由水面上的速度勢是已知的，根據積分方程計算當前時刻物面上的速度勢和自由水面上的速度勢法向導數，然后應用四階Runga-Kutta法，根據自由水面條件式（2）計算下一時刻的水質點位置和自由水面上的速度勢，再用二次形狀函數在舊單元上插值求得新節點上的物理量來對自由水面網格重新劃分，重新應用積分方程計算下一時刻物面上的速度勢和自由水面上的速度勢法向導數。這樣計算周而復始，直到計算結束[16-17]。

求解作用在結構上的波浪力F= ｛fx,fy｝可通過在瞬時物體濕表面Гb上做壓強積分得到

速度勢的時間導數φt也滿足Laplace方程

在自由水面Гf上，φt由Bernoulli方程給出：

在固定邊界上滿足

進而通過求解積分方程

可以求得φt，其中系數矩陣與（6）式中相同，不用重新建立，qt*為源強qt的時間導數，可解析求得。最后通過（8）式求得作用在物體上的波浪力。

2 數值計算及討論

2.1 模型準確性和穩定性

作為算例，本文以Saitoh等（2006）[6]的實驗來進行數值模擬波浪與具有窄縫的兩箱體的相互作用問題，驗證本文數學模型的準確性。這里選用的實驗參數為水槽靜水深h=0.5m，箱體A和B的寬度W=0.5m，吃水深度D=0.252m，入射波高H0=0.024m，箱體間窄縫寬度Wg=0.03m、0.05m和0.07m。在數值模型中，計算域長度取7.5倍波長，在水槽的左右兩端各布置1.5倍波長的阻尼層，造波源位于x=0，箱體A側面邊界b位于距離造波源2.5倍波長的位置，然后依次按Wg調整箱體B的位置。通過開展數值收斂性實驗，自由水面上每個波長布置15個單元，窄縫內布置2個單元，箱體側面邊界b、d、b0、d0均布置6個單元，底面邊界c和c0均布置12個單元；時間步長Δt=T/60 s，每個算例模擬30個周期。

圖3給出了窄縫寬度Wg分別為0.03m、0.05m和0.07m情況下窄縫中心位置V點處無量綱波高Hg/H0隨入射波波數kh的變化關系，及本文數值結果與實驗數據[4]、粘性流模型數值結果[9]的比較。通過數值實驗測算，窄縫中間水面的人工粘性系數分別取為μ2=0.04，0.03，0.03；數值模擬中，波高Hg以測試點的波面時間序列中波峰值與波谷值和的平均來計算。從圖中可以看出，三種窄縫寬度情況下，窄縫內數值模擬波面高度分別在kh=1.7，1.58和1.47時達到最大，也即窄縫內流體發生共振，波高分別為入射波高的4.3倍、5.2倍和5.0倍。整體上兩種數值結果與實驗數據均符合得很好，在個別峰值處甚至本文結果比粘性流模型結果與實驗數據吻合得更好，說明所建立模型在取得合適的人工粘性系數情況下可以準確模擬窄縫內流體共振問題。

圖2 窄縫中間位置V處無因次波高Hg/H0與波數kh間的關系Fig.2 Distribution of non-dimensionalwave heightwith respect to incidentwave frequency

圖3給出了窄縫寬度Wg為0.05m時上述波況下作用在箱體A和B上的x和z向無量綱化波浪荷載幅值隨波數kh的變化，及本文結果與粘性流模型[9]和線性勢流模型[9]結果的對比，圖中波浪荷載幅值以波浪力時間穩定段序列中峰值與相鄰谷值和的一半的平均來計算。從圖中可以看出，作用在兩個箱體上的波浪荷載也均是在共振頻率處達到最大，同時三種數值結果均吻合得很好，也進一步說明本模型可以準確模擬窄縫內共振流體作用在結構上的波浪荷載。

圖3 窄縫寬度Wg為0.05m時作用在兩箱體波浪荷載隨波數的分布Fig.3 Distribution of dimensionlesswave forces on twin boxes against kh at Wg=0.05m

圖4給出了波數kh=1.58，窄縫寬度Wg=0.05 m，入射波高H0=0.024 m情況下t=26T和30T時的整個計算域波面分布，圖中波面空白處為兩個箱體所在的位置。從圖中可以看出，兩個時刻的波面曲線已經完全重合，包括箱體A前的反射波，箱體B后的透射波和兩個箱體之間的窄縫內波面；在水槽兩端的阻尼層吸收波浪的效果都很理想，兩端的波面基本趨于0，箱體A前的波浪形成了穩定的立波，說明其反射回去的波浪透過造波源被前端阻尼層完全吸收，而對入射波浪沒有產生影響。以上現象說明本模型的模擬結果已達到穩定。圖5給出了箱體A和B的水平荷載和垂向荷載的時間歷程，從圖中可以看出，水平力和垂向力在數值上均達到穩定，在共振條件下，窄縫內水體的作用占主導，其作用的側面d和b0法向量相反，箱體A和B的水平力反相位；而作用在箱體A底面c上的反射浪非線性要強于作用在箱體B底面c0上的透射浪，導致了作用在箱體A和B的垂向力在相位和幅值上的差異。

圖4 t=26T和30T兩個時刻的水槽波面分布Fig.4 Snapshotofwave elevation along the wave flume at t=26T and 30T

圖5 箱體的水平和垂向荷載時間歷程Fig.5 Time series ofwave loads on boxes in horizontal and vertical directions

2.2 數值結果

下面以上述算例中窄縫寬度Wg=0.05m為模擬對象，進一步分析入射波浪非線性對波浪爬高和荷載的影響，以及共振條件下壓力分布等問題。圖6給出了共振條件下（即，kh=1.58）入射波高H0分別為0.04m、0.06m和0.08m時窄縫內的波面分布情況，圖中波面通過除以入射波幅（H0/2）進行無量綱化處理。計算中對窄縫內三點II、V、III的波面時間歷程進行記錄，通過對比發現三點波面歷程完全相同，說明窄縫內流體沒有發生x方向波動，而是沿著z向做整體波動，故而可以用窄縫內任一點代表整個流體波動情況。從圖6可以看出，入射波高H0=0.08m對應的曲線無量綱化波高最大，傳播速度最快，H0=0.04m對應的曲線無量綱化波高最小，傳播速度最慢，這說明了大波高入射波浪所體現的強非線性作用。

圖7為入射波高H0=0.024 m時結構迎浪側（I點）、背浪側（IV點）和窄縫內（V點）的無量綱化波高Hg/H0隨入射波數kh的變化關系。從圖中可以看出，在共振頻率處kh=1.58，窄縫內V點和背浪側IV點波高均達到最大，而迎浪側I點則為最小，這是因為當窄縫內共振時能量達到最大，同時會分別有一部分能量向迎浪側和背浪側傳遞，與迎浪側波浪能量作用起到減弱效果，而與背浪側波浪能量作用起到增強的效果。在高頻區域，窄縫內波浪和背浪側波浪都趨于很小值，而迎浪側波高Hg/H0趨于常值2，這是因為對應短波大多被結構迎浪面全反射，進入到窄縫和背浪側的波浪成分很少的緣故。

圖6 共振條件不同入射波幅對應的窄縫內波面時間歷程Fig.6 Time series ofwave elevation atnarrow gap for different inputwave amplitude at resonance

圖7 不同位置無因次波高與入射波頻率的關系Fig.7 Dimensionlesswave height against kh at different positions

圖8 共振時不同位置無因次波高與入射波波高的關系Fig.8 Dimensionlesswave heightagainst H0at different positions at resonance

圖8為共振條件下kh=1.58，迎浪側（I點）、背浪側（IV點）和窄縫內（V點）的無量綱化波高Hg/H0隨入射波高H0的變化關系。從圖中可以看出，在H0≤0.03 m時，三條曲線均近似為水平線，表現為波浪的線性特性；反之，波高Hg/H0均隨著入射波高H0的增加而增大，呈現為波浪的非線性特性，在H0= 0.08 m時，Hg/H0在I點、IV點和V點的值分別達到2.7、5.6和0.46，而在H0=0.01m時，三點處波高僅為入射波高的1.85倍、5.1倍和0.34倍。

圖9 共振條件下箱體A和B的無因次波浪荷載隨入射波高的變化關系Fig.9 Dimensionlesswave loads on boxes A and B against H0at resonance

圖9是共振條件下kh=1.58，箱體A和B的無因次化水平荷載與垂向荷載隨入射波高H0的變化關系。從圖中可以看出，當H0≥0.03m，作用在兩個箱體上的波浪荷載如波高變化一樣呈非線性特性，且總體上作用在迎浪側箱體A的荷載要大于背浪側的箱體B，這與透射到箱體后的波浪貢獻很小相關。從圖9（a）中還可以發現，當H0≥0.07m，作用在箱體A上的水平荷載開始減小，而作用在箱體B上的水平荷載一直處于增大趨勢，分析其原因是入射波高增大，波浪非線性增強，高階波浪貢獻相應增大，而根據立波理論，作用在迎浪側b上的偶數階（尤其是二階差頻項）波浪荷載的相位與線性部分呈反相位，使得箱體A整體波浪荷載下移。

圖10是共振條件下kh=1.58，箱體A和B各個側面上的最大波壓力空間分布，其中圖10（a）中的橫坐標x=0分別代表兩個箱體在窄縫內的下端點。從圖10（a）中可以看出，由于受窄縫內共振波浪的影響，兩箱體內側壓力最大且相當，隨著向箱體外側延伸，底面壓力開始減小，由于透射波浪最小，所以箱體B上的壓力減小最快，而在箱體A的外側邊緣位置，由于其立波的形成，使得該位置壓力又開始增大。從圖10（b）中可以看出，由于窄縫內發生共振，產生很大的水體垂向運動，所以作用在側面d和b0上的壓力是相同的，同時也是各個側面中最占優的，而由于透射浪成分很小，作用在箱體B背浪面d0上的壓力最小。

圖10 窄幅共振條件下箱體底面和側面最大波壓力分布Fig.10 Distribution ofmaximum pressure along bottom and lateral sides of boxes at resonance

3 結 論

本文基于域內源造波的時域高階邊界元方法建立波浪與具有窄縫的兩箱體結構相互作用的完全非線性數值水槽模型，對窄縫內流體共振條件下反射波高、透射波高、窄縫內波高、作用在箱體上的波浪荷載和各側面上的最大壓力分布等進行了模擬研究。通過與已發表實驗數據和數值結果進行對比驗證，表明本文所建立數學模型可以準確模擬波浪與具有窄縫結構相互作用過程，且在較小計算域內可長時間模擬得到穩定的結果，沒有在入射邊界發生二次反射現象。研究發現，窄縫內水體發生共振時，透射波高和窄縫內波高均達到最大，而反射波高則為最小；入射波高H0=0.08m時的反射波高、窄縫內共振波高和透射波高分別比H0=0.01m時增大0.41倍、0.1倍和0.35倍，說明箱體迎浪側非線性最強，窄縫內流體非線性最弱；作用在箱體上的波浪荷載非線性影響明顯，整體上為作用在迎浪側結構上的荷載大于背浪側結構，但當入射波浪非線性增大到一定程度時，由于立波作用引起高階波浪力的相位與線性部分相反，會使得作用在迎浪側上的整體水平力會小于作用在背浪側的結構；共振時，作用在窄縫內側面上的波壓力最大，而作用在迎浪側結構底面上的波壓力則呈現兩端大中間小的特點。本研究還可以進一步拓展到波浪與具有窄縫的多結構相互作用、波浪與具有窄縫的船舶和碼頭之間的相互作用等問題，進而為海洋工程作業安全和設計提供參考和依據。

參 考 文 獻：

[1]Miao G P,Ishida H,Saitoh T.Influence of gaps between multiple floating bodies on wave forces[J].China Ocean Engineering,2000,14(4):407-422.

[2]滕 斌,何廣華,李博寧,程 亮.應用比例邊界有限元法求解狹縫對雙箱水動力的影響[J].海洋工程,2006,24(2)：29-37. Teng B,He G H,Li B N,Cheng L.Research on the hydrodynamic influence from the gap between twin caissons by a scaled boundary finite elementmethod[J].Ocean Engineering,2006,24(2):29-37.

[3]何廣華,滕 斌,李博寧,程 亮.應用比例邊界有限元法研究波浪與帶狹縫三箱作用的共振現象[J].水動力學研究與進展A輯,2006,21(3):418-424. He G H,Teng B,Li BN,Cheng L.Research on the hydrodynamic influence from the gaps between three identical boxes by a scaled boundary finite elementmethod[J].Journal of Hydrodynamics,Ser.A,2006,21(3):418-424.

[4]Saitoh T,Miao G P,Ishida H.Theoretical analysis on appearance condition of fluid resonance in a narrow gap between twomodules of very large floating structure[C]//Proceedings of the Third Asia-Pacific Workshop on Marine Hydrodynamics.Shanghai,China,2006,170-175.

[5]Iwata H,Saitoh T and Miao G P.Fluid resonance in narrow gaps of very large floating structure composed of rectangular modules[C]//Proceedings of the Fourth International Conference on Asian and Pacific Coasts.Nanjing,China,2007,815-826.

[6]Kristiansen T,Faltinsen O M.Studies on resonantwatermotion between a ship and a fixed terminal in shallow water[J]. Journal of Offshore Mechanics and Arctic Engineering,2009,131:021102

[7]Kristiansen T,Faltinsen OM.A two-dimensional numerical and experimental study of resonant coupled ship and pistonmodemotion[J].Applied Ocean Research,2010,32:158-176.

[8]Lu L,Cheng L,Teng B,Zhao M.Numerical investigation of fluid resonance in two narrow gaps of three identical rectangular structures[J].Applied Ocean Research,2010,32(2):177-190.

[9]Lu L,Teng B,Cheng L,Sun L,Chen X B.Modelling ofmulti-bodies in close proximity under water waves-fluid resonance in narrow gaps[J].Science China Physics,Mechanics&Astronomy,2011,54(1):16-25.

[10]Li L,Zhang C.Analysis of wave resonance in gap between two heaving barges[J].Ocean Engineering,2016,117:210-220.

[11]Feng X,BaiW.Wave resonances in a narrow gap between two barges using fully nonlinear numerical simulation[J].Applied Ocean Research,2015,50:119-129.

[12]Brorsen M,Larsen J.Source generation of nonlinear gravity waveswith the boundary integral equationmethod[J].Coastal Engineering,1987,11:93-113.

[13]Ning D Z,Teng B,Eatoack Taylor R,Zang J.Numerical simulation of non-linear regular and focused waves in an infinite water-depth[J].Ocean Engineering,2008,35(8-9):887-899.

[14]Tanizawa K.Long time fully nonlinear simulation of floating bodymotionswith artificial damping zone[J].The Society of Naval Architects of Japan,1996,180:311-319.

[15]Kim YW.Artificial damping in water wave problems I:Constant damping[J].International Journal of Offshore and Polar Engineering,2003,13(2):88-93.

[16]周斌珍,寧德志,滕 斌.造波板運動造波實時模擬[J].水動力學研究與進展A輯,2009,24(4):1-12. Zhou B Z,Ning D Z,Teng B.Real-time simulation of waves generated by a wavemaker[J].Journal of Hydrodynamics, Ser.A,2009,24(4):1-12.

[17]陳麗芬,寧德志,滕 斌,宋偉華.潛堤上波流傳播的完全非線性數值模擬[J].力學學報,2011,43(5):834-843. Chen L F,Ning D Z,Teng B,Song W H.Fully nonlinear numerical simulation for wave-current propagation over a submerged bar[J].Chinese Journal of Theoretical and Applied Mechanics,2011,43(5):834-843.

Nonlinear numerical simulation of wave action on objectsw ith narrow gap

NING De-zhi,SU Xiao-jie,TENG Bin
(State Key Laboratory of Coastal and Offshore Engineering,Dalian University of Technology,Dalian 116024,China)

Based on a time-domain higher-order boundary elementmethod(HOBEM)with wave generation by the inner-domain source,a two-dimensional time-domain numericalwave flume is developed to investigate the fluid resonance due to the interaction between wave and objectswith a narrow gap.In the numericalmodel,the artificial damping is introduced into the fluid at gap and the fully nonlinear boundary conditions are satisfied on the instantaneous free surface.In the solving process,themixed Eulerian-Lagrangianmethod is adopted to track the transientwater surface and the 4th Runga-Kutta technique is used to refresh the velocity potential and free surface at the next time step.The acceleration potential technique is adopted to calculate the transientwave loads along the wetted object surface.By comparison with the published experimental and numerical data,the proposed model is validated.Numerical experiments are performed to study the effects of the incidentwave nonlinearity on reflection wave height,transmission wave height,wave height at narrow gap,wave loads and pressure distribution at resonance.

narrow gap;numericalwave flume;fluid resonance;source generation technique; higher-order boundary element

O353.2

：Adoi:10.3969/j.issn.1007-7294.2017.02.003

2016-07-20

國家自然科學基金資助項目（51679036，51490672）；教育部新世紀優秀人才支持計劃（NCET-13-0076）；中央高校基本科研業務費專項資金（DUT14YQ206）

寧德志（1975-），男，教授，博士生導師，E-mail:dzning@dlut.edu.cn；

蘇曉杰（1989-），男，博士研究生，E-mail:676413719@qq.com。

1007-7294（2017）02-0143-09