， ，

(天津大學 a.水利工程仿真與安全國家重點實驗室， b.建筑工程學院，天津 300072)

0 引 言

為保證海洋結構物的安全性，須研究海浪與結構物的相互作用問題。圓柱形結構物為最簡單也最常見的海上結構物，國內外學者對其波浪爬坡現象進行試驗研究。KRIEBEL[1]進行圓柱形結構波浪爬坡的基準試驗。NIEDZWECKI等[2]采用小比例模型試驗，對固定圓柱結構的波浪爬坡作用進行測量，結果發(fā)現：除了波陡很小的情況之外，線性繞射理論計算出的波浪爬坡結果小于實際值。MARTIN等[3]進行關于圓柱型結構波浪爬坡的試驗，并將試驗結果與不同理論的計算結果進行對比，結果發(fā)現：試驗得到的爬坡結果大于理論計算值。De VOS等[4]進行關于規(guī)則波與不規(guī)則波爬坡的試驗，并通過試驗結果擬合出關于不規(guī)則波爬坡高度的計算公式。除此之外，國內外學者也對波浪爬坡高度進行理論計算研究。KRIEBEL[5]、ISAACSON等[6]、BUCHMANN等[7]使用二階頻域程序計算波浪爬坡高度。BUCHMANN等[8]、ISAACSON等[9]使用二階時域計算方法計算波浪爬坡高度。但這些研究方法計算出的爬坡結果與KRIEBEL的試驗結果相比吻合度不高。

近年來，隨著計算流體力學(Computational Fluid Dynamics，CFD)技術的高速發(fā)展，數值波浪水池技術的研究成為當今的熱點。數值波浪模擬可分為仿物理造波法(推板/搖板造波)和純數值造波法(源項造波)。在仿物理造波方面，顧挺鋒[10]以哈爾濱工程大學的海洋工程水池為原型，采用搖板造波原理模擬穩(wěn)定的線性規(guī)則波、不規(guī)則波，成功模擬出消波灘上的波浪爬高以及破碎現象。封星等[11]采用推板造波原理模擬二維線性規(guī)則波。在數值造波方法方面，李宏偉[12]將搖板造波方法與源函數造波方法的二維造波進行對比，并對三維模型的數值造波進行初步探討。季紅葉[13]構建三維數值波浪水池，模擬單色波及雙色波，并對MARINTEK波浪爬坡試驗進行模擬。

目前，國內外學者對波浪爬坡現象數以及數值波浪水池均進行了相關研究，但對波浪模擬精度、雙色波模擬等方面還需要更深入的研究。本文采用CFD方法，將數值造波理論應用于波浪爬坡的模擬中，使用數值波浪水池對MARINTEK的波浪爬坡試驗進行模擬，并與試驗結果以及其他方法進行對比。

1 數值波浪水池模擬基本理論

1.1 Navier-Stokes方程

數值水池模擬的基本方程為Navier-Stokes方程(簡稱N-S方程)。N-S方程包括連續(xù)性方程和動量方程2個部分。連續(xù)性方程為

(1)

式中：ux、uy、uz分別為x、y、z3個方向的速度。

動量方程為

(2)

式中：t為時間；ρ為流體密度；g為重力加速度；υ為動力黏性系數；p為壓力。

1.2 數值造波方法

采用質量源法進行波浪模擬。質量源法是在連續(xù)性方程中引入質量源函數。源項范圍取水池左側的第1層網格，源項頂端在靜水面下距靜水面1個波高處。對于僅沿x方向傳播的三維波浪問題，使用質量源函數進行造波，連續(xù)性方程應修改為

(3)

式中：q(x,z,t)為質量源函數，計算公式為

(4)

式中：xs為造波源在x方向的位置；qs(z,t)為造波強度。

式(4)表明：在造波源區(qū)域x=xs處造波強度為qs(z,t)，在其他位置造波強度為0。加入質量源后動量方程為

(5)

式中：q為質量源。

造波強度為

(6)

式中：u(xs,z,t)為在造波源處z高度的水質點在t時刻的水平速度；Δx為網格在x方向的長度。

1.3 數值消波方法

采用人工黏性消波法對水池進行消波處理。人工黏性消波法是仿照現實波浪水池中的消波層進行的，其主要思想是在消波區(qū)域中加入阻尼項消除反射波浪。在消波區(qū)域加入阻尼后，此處的動量方程為

(7)

式中：μ(x)為消波系數，是起點為0的單調遞增函數。

本文采用線性消波系數消波，即在消波區(qū)域內，消波系數沿水池水平方向呈線性分布，如圖1所示。

圖1 阻尼消波示意圖

消波系數消波的表達式為

(8)

式中：x0為消波區(qū)域前端起始位置；x1為消波區(qū)末端結束位置；x1-x0為消波區(qū)長度，一般取1個波長；α為消波經驗系數，選取合適的值可使水池獲得最好的消波效果。

1.4 數值水池模擬方法

本文采用有限體積法(Finite Volume Method，FVM)對偏微分方程進行離散，離散格式采用二階迎風格式，壓力速度耦合方法選用壓力隱式分離算法(Pressure Implicit Split Operator，PISO)。對于自由液面的捕捉采用較為常用的流體體積法(Volume of Fluid Method，VOF)。

2 MARINTEK波浪爬坡試驗模型

本文采用MARINTEK波浪爬坡模型試驗測量圓柱周圍各個點的波浪爬坡結果，并將數值水池模擬結果與其他公開發(fā)表的模擬結果進行對比分析。

模型試驗的原型為半潛式平臺的大尺度圓柱，服役海域水深為489 m，圓柱吃水24 m，圓柱直徑為16 m。模型試驗在挪威Trondheim的MARINTEK拖曳水池中進行，水池寬10 m，深10 m。根據原型和試驗水深比值作為模型與實物的比例，取1∶48.93，試驗圓柱直徑為0.327 m，圓柱吃水0.49 m，圓柱中心在水池中距造波區(qū)10 m的位置。波面監(jiān)測點位置如圖2所示，圖中：A1、A2、A3、A4方向分別為270°、225°、202.5°、180°；點1、2、3、4距圓柱中心的距離分別為0.17 m、0.19 m、0.26 m、0.33 m。試驗模擬的波浪條件包括3組線性單色波和3組雙色波，本文分別選取其中1組線性規(guī)則波和1組雙色波進行模擬，模型試驗的波浪狀況見表1。

圖2 波面監(jiān)測點位置示意圖

表1 模型試驗波浪狀況

3 波浪生成模擬

3.1 模型網格設置

根據表1的波浪參數模擬該深水波浪，在模擬時縮小原來拖曳水池的尺寸，選擇數值水池的尺寸為長26 m，寬5 m，深6 m，水深為5 m。消波范圍為沿水池長度20 ～26 m的范圍，造波源范圍為頂端距離靜水面1個波高，x方向寬度為1個網格。波浪載荷計算模型如圖3所示。

圖3 波浪載荷計算模型

為保證模擬的準確性，本文先進行波浪模擬，然后進行波浪與結構物作用的模擬，圓柱體尺寸和位置與模型試驗相同。水池左側即造波源一端采用Symmetry邊界條件，頂端采用Pressure-outlet邊界條件，其余2個面采用Wall邊界條件。由于采用源項造波法，模型網格中沒有動網格存在，對時間步長要求較低，時間步長取0.01 s，計算3 000步，即模擬時間為30 s。

在波浪模擬時，為保證波面的捕捉，在靜水面上、下分別取1個波高的范圍進行網格加密，z方向其他部分采用逐漸變稀疏的方式進行劃分；x方向在造波區(qū)以及波浪工作范圍內采用0.06 m的網格，在消波區(qū)域采用逐漸變稀疏的網格進行劃分；y方向平均分為50份，即采用0.1 m的網格進行劃分。

3.2 波面升高模擬結果

選取10 m處的結果進行對比，在M1單色波工況下的模擬結果如圖4所示。使用濾波程序進行統(tǒng)計，可以得到模擬波浪的周期為1.286 2 s，與理論值誤差為0.031%；模擬計算得到的波高為0.086 8 m，與理論值誤差為0.696%，M1工況的波浪模擬結果較準確。在B1雙色波工況下的模擬結果如圖5所示。將計算得到的波浪進行統(tǒng)計分析，與理論對比的結果見表2。

圖4 M1工況x=10 m處波面升高結果對比 圖5 B1工況x=10 m處波面升高結果對比

表2 B1工況10 m處雙色波統(tǒng)計結果對比

由表2可以看出：跨零周期計算結果準確度較高，但是波高結果存在一定的差異?？赡苡幸韵略颍?1)B1工況網格的劃分密度，對于B11單色波分量來說較疏，導致B11分量模擬結果與理論值相比有一定誤差，而在雙色波中這種誤差被擴大了；(2)源項區(qū)域頂端距離靜水面的高度應該選取1倍波高的距離，此距離偏大時會使得波高偏小，本文中雙色波模擬時采用的是1/3波高統(tǒng)計值，可能會有一定影響。

4 波浪爬坡模擬

4.1 模型網格設置

為了更精確地計算結構物周圍的壓力和波面，在結構物周圍進行網格劃分加密，自由液面處加密情況與3.1節(jié)相同。網格劃分如圖6所示。

圖6 波浪載荷模型結構加密網格俯視圖

4.2 模擬結果對比

爬坡模擬時，分別監(jiān)測A1、A2、A3、A4方向上對應的1、2、3、4點，獲得4個點的液面高度以判斷4個方向上的波浪爬坡情況。將模擬結果與公開發(fā)表的論文的試驗值以及A、B、C、D 4種不同數值模擬方法的計算結果進行對比[14]。方法A～D 的模擬參數見表3，表中：BEF為邊界元技術，自由表面使用格林方程；BER為邊界元技術，使用蘭金-格林方程；FD表示使用頻域計算；TD表示使用時域計算；2nd表示波浪的擾動采用二階算法計算，1st表示波高采用線性算法計算。

表3 不同模擬方法參數對比

線性規(guī)則波的爬坡采用M1波浪狀況模擬計算，對于M1工況本文的數值模擬計算值與方法 A，方法 B以及試驗值的對比結果如圖7所示，其中橫坐標表示測量點距圓柱中的距離與圓柱半徑的比值r/R，縱坐標表示爬坡高度與1/2波高的比值A/(H/2)。

圖7 M1工況不同方向的波浪爬高計算結果對比

由圖7可以看出：4個方向中方法A與方法B的計算結果較吻合；在A1、A2方向上，本文的計算結果與方法A、方法B和試驗值的結果均吻合較好；在A3方向上，試驗值比數值計算值大出很多，而本文計算結果與其他數值模擬結果較為接近；在A4方向上本文的計算結果比試驗值稍大，而其他數值模擬結果比試驗值稍小。

觀察每個方向上幾個點的平均高度并分別將每種計算結果進行比較，對比結果如圖8所示，可以看出：本文計算值和方法B的波面增高按照A1～A4的順序逐漸增大；方法A中A4方向平均波面增高比A3方向接近但稍微小一些；試驗值A3方向的波面增高明顯大于其他幾個方向，其余3個方向的波面增高按照A1、A2、A4的順序逐漸增大，根據實際經驗，A3方向這種突然增大是不合理的，因而A3方向的試驗數據忽略不使用。根據檢測位置的布置情況，A1與波浪傳播方向相互垂直，A4與波浪傳播方向相互平行，A2、A3都與波浪傳播方向呈一定角度。波浪速度在4個方向的分量按照A1～A4的順序逐漸增大，根據動量定理，波浪與結構物的相互作用效果按照A1～A4的順序也逐漸明顯，因此可以推測，波面升高平均值應該按照A1～A4的順序逐漸增大?？梢姡瑔紊ㄓ嬎憬Y果與試驗及其他模擬結果相比基本吻合。

圖8 M1工況不同方向平均波面高度對比結果

雙色波選取B1波浪狀況進行模擬計算，對于B1工況本文的數值模擬計算值與方法B、C、D和試驗值的對比結果如圖9所示，其中橫坐標表示測量點距圓柱中心的距離與圓柱半徑的比值r/R，縱坐標表示爬坡高度A。

圖9 B1工況不同方向的波浪爬高結果對比

由圖9可以看出：在A1方向上，本文的計算結果與其他模擬和試驗的結果基本吻合；在A2和A4方向上，本文計算值大于試驗值，而其他數值模擬結果要小于試驗值，但本文計算結果更接近試驗值；在A3方向上，沒有給出試驗值，本文計算結果與其他模擬結果相差較大，因此A3方向上結果應舍棄，不予考慮。

觀察每個方向上幾個點的平均高度并分別將每種計算結果進行比較，結果如圖10所示，可以看出：除A3方向上試驗值未給出外，其余幾種計算結果的波面升高均按照A1～A4的順序逐漸增大，與實際情況相符。由此可見，雙色波平均波面升高計算結果略大于試驗值，而B、C、D 3種模擬方法的計算結果小于試驗值，但整體變化趨勢基本一致。

圖10 B1工況不同方向平均波面高度對比結果

5 結 論

根據MARINTEK的波浪爬坡模型試驗結果，本文展開數值模擬對比工作。采用數值波浪水池技術，針對單色波和雙色波2種工況，模擬了波浪和圓柱周圍不同位置的波浪爬坡結果。波浪結果與理論計算值進行對比，而波浪爬坡結果與試驗值和其他數值方法進行對比。得出以下結論：

(1) 本文計算得到的線性規(guī)則波的模擬周期及波高結果與理論值相比的誤差均小于1%；雙色波模擬的周期與理論值基本相符，波高存在一定的誤差，但基本滿足計算要求。

(2) 與其他已發(fā)表文獻的數值模擬結果相比，本文計算得到的波浪爬坡結果略大于試驗值，除了A3方向外，其他方向的計算結果比其他幾種模擬方法更接近試驗值。這證明本文所述數值方法的計算精度要高于上述進行對比的其他幾種數值方法。

[ 1 ] KRIEBEL D. Nonlinear Wave Diffraction by Vertical Circular Cylinder Part 1:Diffraction Theory[J].Ocean Engineering,1990, 17(04): 345-377.

[ 2 ] NIEDZWECKI J M, DUGGAL S D. Wave Run-up and Forces on Cylinders in Regular and Random Waves[J]. Journal of Waterway Port Coastal and Ocean Engineering,2004, 118(06): 615-634.

[ 3 ] MARTIN A J,EASSON W J,BRUCE T. Run-up on Columns in Steep, Deep Water Regular Waves [J]. Journal of Waterway Port Coastal and Ocean Engineering, 2001, 127(01): 26-32.

[ 4 ] De VOS L, FRIGAARD P, De ROUCK J. Wave Run-up on Cylindrical and Cone Shaped Foundations for Offshore Wind Turbines[J]. Coastal Engineering, 2007, 54(01): 17-29.

[ 5 ] KRIEBEL D. Non-Linear Wave Interaction with a Vertical Circular Cylinder Part II: Wave Run-up[J]. Ocean Engineering, 1992, 19 (01): 75-99.

[ 6 ] ISAACSON M D S Q, CHEUNG K. Correction Factors for Non-Linear Run-up and Wave Force on a Large Cylinder[J]. Canadian Journal of Civil Engineering, 1994, 21(05): 762-769.

[ 7 ] BUCHMANN B, FERRANT P, SKOURUP J. Run-up on a Body in Waves and Current, Fully Non-Linear and Finite Order Calculations[C]//13th International Workshop on Water Waves and Floating Bodies, 1998.

[ 8 ] BUCHMANN B, SKOURUP J,CHEUNG K. Run-up on a Structure Due to Waves and Current[C]//Proceeding of the Seventh International Offshore and Polar Engineering Conference, 1997.

[ 9 ] ISAACSON M D S Q, CHEUNG K. Time Domain Solution for Wave Current Interaction with a Two Dimensional Body[J].Applied Ocean Research,1993, 15 (01): 39-52.

[10] 顧挺鋒. 海洋工程水池波浪生成的數值模擬[D]. 哈爾濱: 哈爾濱工程大學,2009.

[11] 封星, 吳宛青, 吳文峰, 等.二維數值波浪水槽在FLUENT中的實現[J]. 大連海事大學學報, 2010, 36(03):94-101.

[12] 李宏偉. 數值水池造波方法研究[D]. 哈爾濱: 哈爾濱工程大學, 2009.

[13] 季紅葉. 波浪對大尺度固定柱狀結構物作用的數值水池模擬[D]. 天津: 天津大學, 2016.

[14] NIELSEN F G.Comparative Study on Airgap Under Floating Platforms and Run-up Along Platform Columns[J]. Journal of Marine Structures, 2003(16): 97-134.