胡偉，紀師明，任效忠，王琪，馬玉祥

（1.大連海洋大學基建管理處，遼寧大連116023；2.中交一航局第三工程有限公司，遼寧大連116083；3.大連九成市政設計有限公司，遼寧大連116000；4.大連理工大學運載工程與力學學部船舶工程學院，遼寧大連116024；5.大連理工大學海岸和近海工程國家重點實驗室，遼寧大連116024）

?

作用于橢圓沉箱上的非線性波浪力研究

胡偉1，紀師明2，任效忠3，王琪4，馬玉祥5

（1.大連海洋大學基建管理處，遼寧大連116023；2.中交一航局第三工程有限公司，遼寧大連116083；3.大連九成市政設計有限公司，遼寧大連116000；4.大連理工大學運載工程與力學學部船舶工程學院，遼寧大連116024；5.大連理工大學海岸和近海工程國家重點實驗室，遼寧大連116024）

摘要：為深入研究橢圓沉箱的水動力特性，建立了非線性波浪與結構物相互作用的三維時域數(shù)值計算模型，該模型中N-S方程的求解基于有限差分法，采用VOF方法追蹤流體自由表面，使用部分單元體方法處理橢圓沉箱結構物的曲面邊界，用計算圓沉箱波浪力對該模型進行了驗證，并將該模型應用于相對波高（H/d）、相對波長（kD）和相對長寬比（B/D）對橢圓沉箱波浪力的影響研究。結果表明：相對波高對橢圓沉箱波浪力有重要影響，其值越大波浪非線性的影響越明顯；與沉箱無量綱正向波浪力相比，相對長寬比對無量綱橫向波浪力影響更大。研究表明，本研究中建立的數(shù)值模型可應用于橢圓沉箱波浪力的計算，但實際設計中需依據(jù)建設地點的風浪特征與掩護條件，確定合適的沉箱相對長寬比，以保證結構足夠的安全系數(shù)。

關鍵詞：橢圓沉箱；波浪力；VOF方法；部分單元體

橢圓沉箱是一種應用于開敞式深水碼頭的新型結構型式，該結構由三部分組成，沉箱的前后兩部分是半圓形斷面，中間部分是矩形斷面。為了克服兩排圓沉箱的不均勻沉降，專業(yè)人員設計了橢圓沉箱結構，避免了采用大型圓沉箱對應強浪方向斷面大、波浪壅高的不利影響。橢圓沉箱重力墩首次在大連礦石碼頭二期工程中被實施應用［1］，解決了大型圓沉箱方案碼頭面高程較高的難題，降低了建設成本，并方便運營作業(yè)。之后，橢圓沉箱在青島港董家口港區(qū)40萬t礦石碼頭［2］建設中作為主體結構再次得到應用。橢圓沉箱與傳統(tǒng)小型圓沉箱結構方案相比，減少了50%沉箱數(shù)量，且沉箱上部不搭設聯(lián)結前后沉箱間的聯(lián)系梁，上部預制安裝構件數(shù)量不僅大大減少，而且碼頭整體性更好。橢圓沉箱在使用性能、安全性能和經(jīng)濟效益方面均有突出優(yōu)勢，具有廣泛的推廣應用前景。

橢圓沉箱結構作為新型深水重力碼頭結構型式，需深入研究其水動力特性。任效忠等［3-4］進行了橢圓沉箱波浪力的試驗研究工作，但試驗僅考慮了沉箱相對長寬比B/D為0.63的模型。任效忠等［5］還利用二維源分布法研究了波浪對橢圓沉箱作用的線性波浪力。然而，對于較大的波浪，由線性理論計算的橢圓沉箱波浪力比試驗結果小很多，因此，需要建立非線性時域數(shù)值模型。本研究中，建立了非線性波浪與結構物相互作用的三維時域數(shù)值計算模型，并將該數(shù)值模型應用于相對波高、相對波長和相對長寬比對橢圓沉箱波浪力的影響研究中。模型中N-S方程的求解基于有限差分法，采用VOF方法追蹤流體自由表面［6-8］，采用部分單元體法處理橢圓沉箱的邊界，較好地滿足了該模型的計算效率與計算精度要求。

1　數(shù)值模型

1.1控制方程

基于流體為粘性不可壓縮的基礎假設，流體控制方程采用連續(xù)性方程和N-S方程。流體體積函數(shù)F用于追蹤流體自由面，單元體內流體占有的體積比率用F表示，流體控制方程為

其中：u、v、w分別為x、y、z軸的速度分量；t為時間參數(shù)；g為重力加速度；ρ為流體密度；p為壓力；νk為運動粘性系數(shù)；θ為部分單元體參數(shù)，θ∈［0，1］。計算參數(shù)：ρ=1000 kg/m3，g=9.80 m/s2，νk=1.002×10-6m2/s。

單元體中的F值定義為單元體被流體占有的體積比率，即F=1表示充滿流體，F(xiàn)=0表示無流體，0＜F＜1表示該單元體未充滿流體，一種情形是單元體與自由表面相交的自由表面單元體，另一種情形是含有比單元體尺度小的空氣泡。

1.2數(shù)值計算方法

1.2.1有限差分方法 將交錯網(wǎng)格應用于有限差分求解控制方程，如圖1所示。參數(shù)p、F位于單元體的中心。速度分量u定義在單元體的右側面中心，記為ui+1/2，j，k；v定義在單元體頂面的中心，記為vi，j+1/2，k；w定義在單元體前面的中心，記為wi，j，k+1/2。ARi，j，k、ATi，j，k和AFi，j，k表示3個面非邊界部分面積的比率 （分別對應單元體右側面、頂面和前側面）；而VCi，j，k表示單元體非邊界部分體積的比率。N-S方程離散中時間項的處理采用時間向前差分格式，粘性項基于二階中心差分格式；用基于迎風差分與二階中心差分的混合差分格式處理對流項。

本研究中，應用SOLA-VOF算法讓每個控制體滿足差分方程。每一計算時間步長，求解出各個控制體的速度場和壓力。用前一時刻計算的流場結果代入顯式格式的動量方程，得到當前時刻的流場近似值。為滿足連續(xù)方程，充滿流體的每個計算單元體的壓力和速度必須進行修正。通過反復迭代法修正每一個單元體的壓力，內部流體單元體修正后使連續(xù)方程達到預設精度要求。表面單元體首先滿足自由表面的動力邊界條件，其壓力通過自由表面處壓力和內部流體單元壓力線性插值求解得到。

圖1　單元體結構示意圖Fig.1 Sketch of a staggered cell

每個計算流體單元經(jīng)過上述插值過程得到速度場與壓力，然后通過流體體積函數(shù)F計算得到新時刻的流體分布。流體體積函數(shù)F是階梯函數(shù)，自由面變化由單元體流體通量的計算確定，單元體邊界面上的流體通量采用施主與受主單元體模型計算，最后計算得到F函數(shù)。

1.2.2邊界條件 本研究中，三維數(shù)值波浪水池如圖2所示，L為水池長度 （m），B為水池寬度（m），d為水深 （m），L1為水池有效區(qū)域長度（m），L2為速度衰減區(qū)域長度 （m）。

圖2　數(shù)值波浪水池示意圖Fig.2 Sketch of a numerical wave tank

在水池的左側 （S1=SA1A2D2D1）設置造波推板。依據(jù)可吸收造波理論，造波中產(chǎn)生行進波的同時還產(chǎn)生附加波動將水池反射抵消，以消除反射波對計算的影響，附加波動與反射波幅值相等，相位相反［9-10］?；跈E余波理論，造波推板速度表達式為

其中：η（k，t）為第kth塊造波板前的實際波面；η0（k，t）為依據(jù)橢余波理論需要的波面；d為水深；c為橢余波波速。

采用數(shù)值衰減與開邊界兩種處理方法，減小波浪反射對數(shù)值計算的干擾。與波浪水池右端（S2= SB1B2C2C1）相鄰的L2區(qū)域是數(shù)值衰減區(qū)域，L2區(qū)域中垂向速度的衰減系數(shù)為

其中：r0為數(shù)值衰減的起點；L2為數(shù)值衰減區(qū)的長度。數(shù)值衰減區(qū)中的垂向速度為由控制方程得到的垂向速度乘以衰減系數(shù)。為進一步減小反射干擾，依據(jù)Sommerfeld輻射條件，右側開邊界的出流速度可以表示為

其中：k為波數(shù)；ω為角頻率；η為開邊界處的自由面。

橢圓沉箱邊界由部分單元體方法處理。波浪水池的側面邊界（S5=SA1B1C1D1、S6=SA2B2C2D2）和底面邊界（S3=SA1A2B2B1）設置為自由滑移直墻邊界。在數(shù)值計算中，自由滑移直墻邊界的法向速度和切向速度的法向導數(shù)均設為0。

2　數(shù)值模型的驗證

2.1數(shù)值波浪水池的穩(wěn)定性

本研究中，對數(shù)值波浪水池的主要指標數(shù)值耗散、計算穩(wěn)定性和造波重復性均進行了驗證。波高H=1.0～4.0 m，周期T=8.0～13.0 s，水深d=20 m。水池總長L=800 m，其中水池有效區(qū)域長L1= 600 m，水池寬度B=260 m。單元體網(wǎng)格尺度為Δx=2 m、Δy=3 m和Δz=2 m。表1列出了沿水池距造波板一倍波長 （λ）、兩倍波長 （2λ）、三倍波長 （3λ）、四倍波長 （4λ）位置處的5個波的平均波高。圖3為位于兩倍波長處的典型波面歷時曲線。圖4為數(shù)值波浪水池的瞬時波面。從數(shù)值計算結果可知，數(shù)值波浪水池的波浪均勻性和重復性均比較好，故數(shù)值波浪水池能滿足計算所需。

2.2波浪力計算

表1　數(shù)值水池中不同位置的波高值Tab.1 Wave height at different positions in a tank

圖3　數(shù)值波浪水池波面歷時曲線 （x=2λ）Fig.3 Time history of wave surface at location（x=2λ）

圖4　數(shù)值波浪水池波面分布示意圖Fig.4 Instantaneous wave surface in a tank

在弱非線性條件下，用本研究中建立的數(shù)值模型計算直立圓柱水平波浪力，并將數(shù)值計算結果與波浪繞射理論［11-14］的解析結果進行對比。計算時，設水深d=20 m，波高H=1.9 m，波浪周期T= 8.0～13.0 s，圓柱半徑=10 m。直立圓柱的曲面邊界采用部分單元體處理，邊界周圍的單元節(jié)點平面坐標如圖5所示。直立圓柱水平波浪力的數(shù)值計算結果及比較見表2，數(shù)值計算結果比解析結果大7%～11%。

圖5　直立圓柱網(wǎng)格示意圖Fig.5 Sketch of cells on a vertical circular cylinder

表2　直立圓柱水平波浪力與解析解的比較Tab.2 Comparison of the horizontal wave force on a vertical circular cylinder with analytical solution

3　橢圓沉箱結構波浪力的計算

橢圓沉箱結構示意圖見圖6，由3部分組成。沉箱的前部分和后部分為半圓形斷面，半圓直徑為D；沉箱的中間部分為矩形斷面，矩形長度為B。

計算參數(shù)中：設H=1.0～4.0 m，T=8.0～13.0 s，d=20 m，相對波長kD=0.7～1.5；橢圓沉箱半圓形斷面直徑為20 m，相對長寬比B/D= 0～0.8；入射波角度α=0°和α=45°；迭代過程收斂精度ε設為0.005。采用部分單元體法處理橢圓沉箱曲面邊界，橢圓沉箱邊界周圍的單元節(jié)點平面坐標如圖7所示 （單位為m）。計算設置初始時間步長Δt=0.04 s，計算過程中時間步長能實現(xiàn)自動調整，以符合Courant條件和擴散穩(wěn)定條件。

通過本研究中建立的數(shù)值模型計算得到橢圓沉箱的橢余波波浪力，并進行無量綱化處理。無量綱正向波浪力和無量綱橫向波浪力的表達式為

圖6　橢圓沉箱結構示意圖Fig.6 Sketch of an ellipse caisson

圖7　圓沉箱計算網(wǎng)格示意圖 （B/D=0.6）Fig.7 Sketch of the cells on ellipse caisson（B/D=0.6）

將B/D=0.6及不同相對波高下橢圓沉箱波浪力的數(shù)值計算結果和基于二維源分布法的線性波浪力結果進行比較，結果見圖8-A、B。圖8中，實心圖例和實線代表數(shù)值計算結果，空心圖例和虛線代表線性結果。比較結果顯示：數(shù)值計算的無量綱正向波浪力和無量綱橫向波浪力均比線性結果大；隨著相對波高H/d的增大，非線性特征的影響在計算結果中體現(xiàn)的越明顯，說明H/d對橢圓沉箱波浪力有重要影響。H/d越大，用數(shù)值模型計算的無量綱波浪力越大，越不容忽視波浪非線性的影響。

圖8　橢圓沉箱非線性和線性波浪力計算結果的對比（B/D=0.6）Fig.8 Comparison of nonlinear and linear wave forces on the quasi-ellipse caisson（B/D=0.6）

圖9-A、B分別為波浪入射角度為0°和45°時，不同相對長寬比下橢圓沉箱波浪力的數(shù)值計算結果，其中，H/d=0.15，B/D=0～0.8。隨著相對波長kD的增大，無量綱波浪力和呈現(xiàn)增大趨勢；隨著B/D的增大，減小，而呈增大趨勢，相對，B/D對有更大影響。

圖9　不同相對長寬比與無量綱橢圓沉箱波浪力的關系（H/d=0.150）Fig.9 Nondimensional wave force of the quasi-ellipse caisson versus kD at different relative lengthwidth ratio（H/d=0.150）

因此，實際設計中選擇沉箱B/D時，應依據(jù)建設地點的風浪特征與掩護條件，確定合適的沉箱B/D，在保證結構足夠安全的同時，充分發(fā)揮橢圓沉箱的優(yōu)勢。在同樣碼頭建設寬度條件下，減小沉箱正向迎浪面積，可增大沉箱間的間距，削弱沉箱間的相互干擾，同時也縮減了建設工程量與工程投資。

4　結論

（1）本研究中基于有限差分法，并結合VOF方法和部分單元體方法，建立了非線性波浪與結構物相互作用的三維時域數(shù)值計算模型，用該數(shù)值模型模擬了橢圓沉箱的非線性波浪力，并得到橢圓沉箱波浪力的非線性特征。對于大型結構橢圓沉箱，波浪力的非線性特征隨相對波高的逐漸增大而更加明顯，波浪非線性的影響不容忽視。

（2）隨著沉箱相對長寬比的增大，無量綱正向波浪力逐漸減小，而無量綱橫向波浪力逐漸增大。與無量綱正向波浪力相比，沉箱相對長寬比對無量綱橫向波浪力影響更大。波浪斜向入射（45°）結果顯示，當沉箱相對長寬比為0.6和0.8時，無量綱橫向波浪力明顯大于無量綱正向波浪力。

（3）選用合適的沉箱相對長寬比是設計中必須重視的問題。對掩護條件好、大風向和大浪向常年集中的海域，適合采用相對長寬比較大的橢圓沉箱，充分發(fā)揮橢圓沉箱的優(yōu)勢；而掩護條件差、復雜多變的海域宜選擇相對長寬比較小的橢圓沉箱，以增強沉箱的抗傾覆能力，保證結構使用安全。

參考文獻：

［1］白景濤，胡家順.橢圓沉箱墩式碼頭新結構的研究與設計［J］.水運工程，2006（11）：25-30.

［2］張志明，胡家順，鄭小楠，等.青島港董家口港區(qū)40萬t礦石碼頭設計關鍵技術［J］.水運工程，2011，457（9）：61-67，127.

［3］Ren Xiaozhong，Wang Yongxue，Wang Guoyu.The wave force on the multiple quasi-ellipse caissons［C］//Proceedings of the Fourth International Conference on Asian and Pacific Coasts.Nanjing China：APAC，2007：201-208.

［4］任效忠，王永學，王國玉.準橢圓沉箱群墩結構波浪力試驗研究［J］.大連理工大學學報，2009，49（6）：944-950.

［5］任效忠，王永學，王國玉.不規(guī)則波作用下準橢圓沉箱單墩的波浪荷載［J］.船海工程，2009，38（1）：85-89.

［6］Hirt C W，Nichols B D.Volume of Fluid（VOF）method for the dynamics of free boundaries［J］.Journal of Computational Physics，1981，39（1）：201-225.

［7］萬德成，繆國平.數(shù)值模擬波浪翻越直立方柱［J］.水動力學研究與進展，1998，13（3）：363-370.

［8］齊鵬.近岸波浪對結構物作用耦合數(shù)值模型及其工程應用［D］.大連：大連理工大學，2001.

［9］李曉磊，欒曙光，陳勇，等.立方體人工魚礁背渦流的三維渦結構［J］.大連海洋大學學報，2012，27（6）：572-577.

［10］Wang Yongxue，Zang Jun，Qiu Dahong.Numerical model of cnoidal wave flume［J］.China Ocean Engineering，1999，13（4）：391-398.

［11］王永學，Su T C.圓柱容器液體晃動問題的數(shù)值計算［J］.空氣動力學學報，1991，9（1）：112-119.

［12］Maccamy R C，F(xiàn)uchs R A.Wave forces on piles：a diffraction theory［R］.Technical Memorandum No.69，California：US Beach E-rosion Board，1954.

［13］桂勁松，溫志超，畢恩凱，等.漁港建設標準中碼頭岸線長度的確定［J］.大連海洋大學學報，2015，30（5）：558-562.

［14］李玉成，滕斌.波浪對海上建筑物的作用［M］.2版.北京：海洋出版社，2002.

中圖分類號：U656.22

文獻標志碼：A

DOI：10.16535/j.cnki.dlhyxb.2016.03.019

文章編號：2095-1388（2016）03-0338-06

收稿日期：2016-03-25

基金項目：國家自然科學基金資助項目 （51422901）

作者簡介：胡偉 （1968—），男，工程師。E-mail：190045408@qq.com

Numerical study of nonlinear wave force on ellipse caisson

HU Wei1，JI Shi-ming2，REN Xiao-zhong3，WANG Qi4，MA Yu-xiang5
（1.Construction Management Department，Dalian Ocean University，Dalian 116023，China；2.No.3 Engineering Co.Ltd.，F(xiàn)irst Harbor Engineering Co.Ltd.，China Communications Construction Co.，Dalian 116083，China；3.Dalian Jiucheng Municipal Design Co.Ltd.，Dalian 116000，China；4.School of Naval Architecture and Ocean Engineering，F(xiàn)aculty of Vehicle Engineering and Mechanics，Dalian University of Technology，Dalian 116024，China；5.State Key Laboratory of Coastal and Offshore Engineering，Dalian University of Technology，Dalian 116024，China）

Abstract：A three dimensional numerical model of interaction of nonlinear wave with structures was developed in this paper to investigate the wave force on ellipse caisson.Navier-Stokes equations were established by a finite difference method，the problem of structure surface boundary was dealt by a partial cell method and the volume of fluid（VOF）method was employed to trace the free surface.The numerical model was verified by calculation of the wave force on circular caisson，and is used to investigate the effects of the relative wave height（H/d），relative caisson width（kD），and relative length-width ratio（B/D）on the wave forces of the ellipse caisson.It was shown that the relative wave height had a significant effect on the wave forces of the caisson，with the increase in the relative wave height，with more obvious nonlinear characteristics.Compared with the non-dimensional inline wave force，the relative length-width ratio was found to have significant influence on the non-dimensional transverse wave force.Since it is very important to consider the characteristics of wind wave and sheltering conditions for selecting the appropriate relative length-width ratio，ensuring the structure having enough safety coefficient，the numerical model in this paper can be used to calculate the wave force on ellipse caisson.

Key words：ellipse caisson；wave force；VOF method；partial cell method