唐　鵬， 于定勇， Bai Wei， 田　艷

(1.中國海洋大學(xué)工程學(xué)院，山東 青島 266100; 2.新加坡國立大學(xué)，新加坡 119077)

?

海洋工程中直立圓柱波浪爬升問題的數(shù)值研究?

唐鵬1， 于定勇1， Bai Wei2， 田艷1

(1.中國海洋大學(xué)工程學(xué)院，山東 青島 266100; 2.新加坡國立大學(xué)，新加坡 119077)

本文基于Navier-Stokes方程和連續(xù)性方程，利用VOF方法捕捉自由液面和松弛造波法實行造波，建立基于黏性流體理論的三維數(shù)值波浪水槽，并在該數(shù)值波浪水槽內(nèi)對圓柱周圍的波浪爬升效應(yīng)進(jìn)行數(shù)值模擬。選用兩種不同的波浪參數(shù)，探討波陡參數(shù)對柱體周圍的波浪爬升效應(yīng)及其所受載荷的影響。與勢流模型相比，本文采用的黏性流模型可以較好的捕捉實驗中出現(xiàn)的二次波峰現(xiàn)象。計算表明，本文建立的數(shù)值水槽可以較好用于直立圓柱波浪爬升問題的計算。

波浪爬升；圓柱；數(shù)值水槽；OpenFOAM

引用格式：唐鵬， 于定勇， Bai Wei， 等. 海洋工程中直立圓柱波浪爬升問題的數(shù)值研究[J]. 中國海洋大學(xué)學(xué)報(自然科學(xué)版), 2016, 46(10): 116-122.

TANG Peng, YU Ding-Yong, Bai Wei， et al. Numerical simulation of wave run-up on cylindrical offshore structures[J].Periodical of Ocean University of China, 2016, 46(10): 116-122.

波浪在傳播的過程中遇到結(jié)構(gòu)物時，會產(chǎn)生繞射效應(yīng)，波浪爬升可描述為入射波浪在傳播過程中遇部分浸入水中的結(jié)構(gòu)物時波浪沿其表面向上爬升的非線性現(xiàn)象。波浪爬升過程的物理解釋為: 當(dāng)入射波浪遇部分浸入水中的柱體時，其周圍的自由波面將發(fā)生劇烈的變化，其中部分水體繞過立柱繼續(xù)朝前傳播，而受立柱阻礙的另一部分波浪開始實現(xiàn)能量的轉(zhuǎn)化，即由動能轉(zhuǎn)化成勢能，導(dǎo)致波浪沿立柱表面迅速向上攀升[1]。對海洋平臺，隨著波浪爬升往往會造成平臺甲板上浪和對平臺底板抨擊作用，從而引起海洋平臺的破壞。為避免波浪爬升對浮式平臺的破壞，在設(shè)計中，應(yīng)該充分考慮波浪爬升，對波浪爬升幅度有準(zhǔn)確的估計，保證平臺甲板有足夠的標(biāo)高。同時，平臺甲板也不能過高，甲板過高往往會產(chǎn)生工程成本增加和對平臺穩(wěn)定性方面產(chǎn)生影響，因此對結(jié)構(gòu)物爬升問題研究具有非常重要的現(xiàn)實意義。

Galvin和Hallermeier等[2]首次采用物理模型試驗方法對圓柱的波浪爬升進(jìn)行了研究。其通過采用安裝于圓柱壁面附近的多個浪高儀來獲取其圓柱周圍自由水面的波形分布。Nielsen[3]采用規(guī)則波對圓形和帶導(dǎo)角的方柱進(jìn)行了物理模型試驗，研究了兩種柱體形狀對波浪爬升效應(yīng)和柱體周圍波面分布的影響，其結(jié)果表明，與圓柱相比，帶導(dǎo)角的方柱周圍的局部波浪放大率比圓柱更大，非線性效應(yīng)更強(qiáng)。Contento G[4]等針對波陡較大、深水規(guī)則波浪下圓形直立柱表面的波浪爬升效應(yīng)進(jìn)行試驗研究。當(dāng)A/D?1(A為入射波波幅，D為圓柱直徑)時候，必須考慮黏性效應(yīng)的影響，黏性效應(yīng)對波浪爬升的貢獻(xiàn)較大，不宜采用勢流理論求解。Donald G. Danmeier[5]基于黏性流模型，利用ComFLOW軟件數(shù)值模擬了規(guī)則波作用下重力式平臺的波浪爬升效應(yīng)。結(jié)果表明，入射波陡較小時，線性波浪爬升占主導(dǎo)，隨著波陡的不斷增大，平臺周圍的波峰逐漸尖瘦，波谷趨于平坦，波浪爬升的高階非線性特征越來越明顯，且在平臺的迎浪面和側(cè)面位置，一個波浪周期內(nèi)出現(xiàn)二次波峰現(xiàn)象。

OpenFOAM(Open Source Field Operation and Manipulation)是一個在Linux下運(yùn)行的計算流體力學(xué)(CFD)類庫，基于有限體積法FVM(Finite Volume Method)求解偏微分方程，是一個針對不同的的流動編寫不同的C++程序集合。本文利用有限體積法求解N-S方程、VOF方法捕捉自由液面和松弛造波法實行造波，建立三維數(shù)值波浪水槽。通過修改OpenFoam中已有兩相流求解器interFoam中的interFoam.C文件編譯適合本文的求解器waveFoam，進(jìn)行波浪與直立圓柱爬升問題的數(shù)值研究。

1　數(shù)值方法

1.1 控制方程

本文中，固液兩相均為不可壓縮流體，不可壓縮黏性流的連續(xù)方程和N-S方程表示為：

連續(xù)方程：

(1)

N-S方程：

(2)

式中：ui為速度，i,j=1,2,3代表三個不同維度，即u1=u，u2=v，u3=w；p為壓力；ρ為流體密度；μ為分子黏性系數(shù)；t為時間；gi為重力。

VOF方法[6]是目前常用捕捉自由液面的方法之一，通過求解體積分?jǐn)?shù)函數(shù)來捕捉自由液面。

(3)

γ為流體單元的體積分?jǐn)?shù)，0≤γ≤1。在空氣中γ=0，在水中γ=1，自由表面位置0<γ<1。

1.2 數(shù)值水槽

本文采用Jacobsen N G[7]提出的松弛造波法實現(xiàn)數(shù)值造波和消波，建立三維數(shù)值波浪水槽。松弛造波法首先是由Mayer等[8]在利用海綿層(Sponge layer)消波方法中提出的。松弛造波法的基本原理是在造波區(qū)和消波區(qū)內(nèi)每一時刻對速度以如下形式進(jìn)行更新。

u=(1-αR)utarget+αRumodel，

(4)

式中；umodel是通過求解N-S方程得到的速度值；utarget是期望得到的目標(biāo)波浪速度值；αR是與空間位置有關(guān)的加權(quán)函數(shù), 在邊界位置αR=0。可以通過下式來確定：

(5)

2　數(shù)值計算和結(jié)果分析

2.1 三維數(shù)值波浪水槽的建立和驗證

為了驗證本文數(shù)值計算方法是否正確，將利用新建立的三維數(shù)值波浪水槽模型來實現(xiàn)線性波的模擬。該數(shù)值波浪水池長20m，寬度為3m，高度為1m，水深h=0.65m。x軸正向與波浪傳播方向一致。反射波往往會影響到實驗的結(jié)果，因此消波效果的對于數(shù)值水槽模擬非常重要，本算例在出口處3m設(shè)置為消波區(qū)，模型相關(guān)參數(shù)如圖2所示。邊界條件設(shè)置為：入口邊界(inlet)patch；出口邊界(outlet)patch；底部邊界(bottom)wall；頂部邊界(atmosphere)patch；前后壁面(front和back)邊界wall。采用的波浪參數(shù)為：波高H=0.08m，周期T=1.5s，波長λ=3.058m。

圖1　入口和出口區(qū)域松弛層中αR和χR變化圖Fig.1　Sketch of the variation of αR(χR) for both inlet and outlet relaxation zones

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

合理的網(wǎng)格分布對生成理想的波浪場是非常重要的，本文的數(shù)值計算中，在水平方向上采用均勻網(wǎng)格,在豎直方向上采用非均勻網(wǎng)格，自由液面附近對網(wǎng)格進(jìn)行加密處理從而更加準(zhǔn)確的獲得自由表面位置。表1給出了數(shù)值計算中采用的3種網(wǎng)格參數(shù)。

表1　數(shù)值計算網(wǎng)格參數(shù)Table 1　Mesh information for simulation

建立數(shù)值波浪水槽是為了實現(xiàn)對波浪條件的模擬，評價和驗證一個數(shù)值波浪水槽是否可行的一個重要參數(shù)為波高，圖3給出了x=7.25m處位移歷程曲線與理論解的對比圖，可以看出生成波形穩(wěn)定，數(shù)值計算波高接近理論波高值，說明本模型可以很好的實現(xiàn)對波浪條件的模擬。圖4給出了t=15s時刻下的3種網(wǎng)格生成的波形與理論解的對比圖，從圖中可以看出，當(dāng)波浪傳遞到水槽的最右邊界時，波高已經(jīng)趨于零，說明本模型設(shè)置的消波區(qū)有很好的消波作用，消波區(qū)能夠很好地吸收波能。隨著網(wǎng)格的加密，三維數(shù)值水槽生成的波形圖更加符合理論解。同時，由于波浪的非線性產(chǎn)生的紊動耗散產(chǎn)生的原因，產(chǎn)生波形與理論解也存在一定差異。

圖3　x=7.25m處位移歷程曲線與理論解的對比圖Fig.3　Comparison of the elevation between simulation and analytical result at x=7.25m

圖4　t=15s時刻下的3種網(wǎng)格的波形與 理論結(jié)果的對比圖Fig.4　Comparison of the wave profile between simulation and analytical results at t=15s

因此，本文建立的三維數(shù)值波浪水槽可以較好地模擬波浪水槽的基本特性，可以用于模擬波浪水槽中進(jìn)行的斷面及三維實驗。

2.2 波浪爬升數(shù)值計算

本算例模擬的是規(guī)則波繞過直立圓柱的爬升現(xiàn)象。幾何模型采用類似上節(jié)數(shù)值水槽的幾何模型。長度x方向為20m，其中x=17～20m設(shè)置為消波區(qū)，同時在x=0～4m設(shè)置入口消波區(qū)，用來消除波浪與結(jié)構(gòu)物相互作用產(chǎn)生的反射波，消除反射波對入口造波的影響。寬度z方向為7m，高度y方向為1m，水深h=0.65m，圓柱直徑d=2r=0.5m，距離入口邊界為7.5m。采用兩種波浪參數(shù)，算例一：波高H=0.08m，周期T=1.5s，波長λ=3.058m，波陡H/λ=0.026；算例二：波高H=0.19m，周期T=1.5s，波長λ=3.058m，波陡H/λ=0.062。

數(shù)值計算中，在水平方向上采用非均勻網(wǎng)格,靠近圓柱處采用加密網(wǎng)格，Δxmin=Δzmin=0.005m，在豎直方向上采用非均勻網(wǎng)格，自由液面附近對網(wǎng)格進(jìn)行加密處理從而更加準(zhǔn)確的獲得自由表面位置，Δymin=0.01m。為考慮波浪的紊流對結(jié)構(gòu)物的影響，本文采用基于k-ε兩方程湍流模型，其中k為湍動能，ε為湍動能耗散率。圖5給出了圓柱周圍網(wǎng)格細(xì)部圖。數(shù)值計算中采用的最大柯朗數(shù)為0.25，采用可變時間步長保證柯朗數(shù)滿足要求。

圖6給出了算例一和算例二圓柱在波浪作用下的無量綱化的水平波浪力歷程曲線。其中：ρ為流體密度；g為重力加速度；r為圓柱半徑；a為入射波的波幅?？梢钥吹接捎诓ɡ说闹芷谛赃\(yùn)動，圓柱也受到水平的周期振蕩力。

圖6　水平波浪力歷程曲線Fig.6　Time history of the horizontal force on the cylinder

圖7給出了算例一中圓柱表面在x=7.25m處位移與圖3中水槽處位移歷程對比圖，即有結(jié)構(gòu)物和無結(jié)構(gòu)物相同位置處的自由表面處位移隨時間變化對比曲線。從圖中明顯可以看出，有結(jié)構(gòu)物的自由表面處位移增大，在t=9.5、11、12.5、14和15.5s時波峰作用該位置波幅達(dá)到最大，最大波幅為a=0.0477m，意味著波浪爬升位移達(dá)到0.0077m，爬升幅度與波幅的比值為19.25%。圖8給出了算例二中圓柱表面在x=7.25m圓柱自由表面處位移圖，其中最大波幅為a=0.231m，意味著波浪爬升幅度達(dá)到0.136m，爬升幅度與波幅的比值為143.15%，接近1.5倍的入射波波幅。Isaacson[9]和Kriebel[10]模型試驗表明在極端情況下，波浪作用在圓柱的爬升幅度最大可以超過1.5倍的波幅。隨著波陡H/λ的增大，圓柱周圍的波峰逐漸尖瘦，有更多的動能轉(zhuǎn)化為勢能，導(dǎo)致波浪沿立柱的爬升幅度增大，尤其在圓柱的迎浪面，因此在波浪作用下形成的波浪爬升對于結(jié)構(gòu)物的影響是不能忽視。同時發(fā)現(xiàn)當(dāng)波面下降到最低位置時，波面出現(xiàn)了小的波動，也就是二次波峰(Secondary crest)現(xiàn)象。

圖7　自由表面處位移歷程對比圖Fig.7　Comparison of the surface elevation time history

圖8　算例二 圓柱表面在x=7.25m處位移圖Fig.8　Time history of the surface elevation at x=7.25m in case 2

圖9　算例二 圓柱表面在x=7.5m處位移圖Fig.9　Time history of the surface elevation at x=7.5m in case 2

圖9給出了算例二中圓柱表面在x=7.5m處自由表面位移圖，相比圖8可以更加明顯發(fā)現(xiàn)，當(dāng)波面下降到最低位置時，波面出現(xiàn)了較大波動，波浪爬升的高階非線性特征越來越明顯，且在圓柱的迎浪面和側(cè)面位置，一個波浪周期內(nèi)出現(xiàn)二次波峰現(xiàn)象。Morris Thomas[11]指出二次波峰現(xiàn)象是由于入射波的波前在繞過柱體兩側(cè)之后在柱子背后相遇并相互干擾，由此產(chǎn)生的自由波會向上游方向回流并與上游過來的波相互疊加，從而導(dǎo)致二次波峰的出現(xiàn)。Trulsen[12]基于完全非線性勢流模型數(shù)值計算了圓柱爬升問題，將數(shù)值計算結(jié)果與實驗結(jié)果對比中，實驗結(jié)果中也出現(xiàn)了二次波峰現(xiàn)象，而基于完全非線性勢流模型和二階勢流結(jié)果沒有捕捉到二次波峰現(xiàn)象，說明二次波峰現(xiàn)象確實存在，勢流模型無法捕捉到該現(xiàn)象，而本文基于黏性流模型可以很好的捕捉這一現(xiàn)象。

圖10和11分別給出了算例一和算例二圓柱不同點(diǎn)在一個周期上的位移圖。從圖中可以看到，圖形以α=180°呈對稱變化，在t=(7+5/15)T時，圓柱迎浪面處(α=180°)表面位移達(dá)到最大。同時，可以明顯看出，算例二(圖11)圓周上位移變化幅度比算例一(圖10)大，說明在波陡較大的情況下，圓柱的自由表面變化更加劇烈。

圖10　算例一圓周不同點(diǎn)在一個周期上的位移圖Fig.10　Wave profiles on the waterline of the cylinder at a whole wave period in case 1

圖11　算例二 圓周不同點(diǎn)在一個周期上的位移圖Fig.11　Wave profiles on the waterline of the cylinder at a whole wave period in case2

圖12給出了算例二在一個完整周期圓柱附近的自由表面位移等值量圖。從圖中可以清楚的看到在一個完整波浪周期作用下圓柱上的位移變化。與圖11對應(yīng)，在t=(7+5/15)T時，圓柱迎浪面(α=180°)處自由表面位移達(dá)到最大，即在此處波浪爬升幅度最大。

圖12　算例二 在一個完整周期圓柱附近的 自由表面位移等值量圖Fig.12　Wave profiles contour of the cylinder at a whole wave period in case2

3　結(jié)論

(1)基于Navier-Stokes方程和連續(xù)性方程，利用VOF方法捕捉自由液面和松弛造波法實行造波，建立基于黏性流體理論的三維數(shù)值波浪水槽。通過數(shù)值計算結(jié)果與理論結(jié)果對比，對比結(jié)果符合性較好，可以證明本文建立的三維數(shù)值波浪水槽可以很好的實現(xiàn)數(shù)值造波和消波，較好地模擬了波浪水槽的基本特性。

(2)利用建立的數(shù)值水槽對圓柱周圍的波浪爬升效應(yīng)進(jìn)行數(shù)值模擬。選用兩種典型工況條件，其中當(dāng)入射波陡較小時，波浪爬升以線性特征為主，波浪爬升波幅較小。隨著波陡的不斷增大，圓柱周圍生成的波峰逐漸尖瘦，波谷開始出現(xiàn)波動，出現(xiàn)了與實驗相似的二次波峰現(xiàn)象，波浪爬升波幅增大，接近1.5倍的入射波波幅，此時波浪爬升的高階非線性特征占主導(dǎo)。與勢流模型相比，本文采用的黏性流模型可以較好的捕捉實驗中出現(xiàn)的二次波峰現(xiàn)象。

(3)在入射波波陡較大的情況下，通過比較圓周不同點(diǎn)在一個周期上的位移可以發(fā)現(xiàn)，圓柱附近的自由表面變化幅度也較大。在圓柱迎浪面(α=180°)處自由表面位移達(dá)到最大。

[1]單鐵兵, 楊建民, 李欣, 等. 深海平臺立柱周圍波浪非線性爬升研究進(jìn)展[J]. 海洋工程, 2012,30(1):151-160.

Shan T, Yang J, Li X, et al. Review of the research on non-linear wave run-up around columns of deepwater platform[J]. The Ocean Engineering,2012, 30(1):151-160.

[2]Galvin C J,Hallermeier R J. Wave runup on vertical cylinders[J]. Coastal Engineering Proceedings, 1972,1(13):65-71.

[3]Nielsen F G. Comparative study on airgap under floating platforms and run-up along platform columns[J]. Marine structures, 2003,16(2):97-134.

[4]Contento G, Francescutto A, Lalli F. Nonlinear wave loads on single vertical cylinders:Pressure and wave field measurements and theoretical predictions[C]. San Francisco： International Society of Offshore and Polar Engineers, 1998.24-32.

[5]Danmeier D G, Seah R K, Finnigan T, et al. Validation of wave run-up calculation methods for a gravity based structure[C]. San Francisco： International Conference on Offshore Mechanics and Arctic Engineering, 2008.265-274.

[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]Jacobsen N G, Fuhrman D R, Freds E J. A wave generation toolbox for the open‐source CFD library: OpenFoam[J]. International Journal for Numerical Methods in Fluids, 2012,70(9):1073-1088.

[8]Mayer S, Garapon A, Soerensen L S. A fractional step method for unsteady free-surface flow with applications to non-linear wave dynamics[J]. International Journal for Numerical Methods in Fluids, 1998,28(2):293-315.

[9]Isaacson MdSQ. Wave runup around large circular cylinder[J]. Journal of the Waterway Port Coastal and Ocean Division, 1978, 104(1): 69-79.

[10]Kriebel D L. Nonlinear wave interaction with a vertical circular cylinder. Part II: Wave run-up[J]. Ocean Engineering, 1992, 19(1): 75-99.

[11]Morris-Thomas M T. Wave run up on vertical columns of an offshore structure[D]. Perth :The University of Western Australia, 2003.

[12]Trulsen K, Teigen P. Wave scattering around a vertical cylinder: Fully nonlinear potential flow calculations compared with low order perturbation results and experiment[C]. San Francisco： International Conference on Offshore Mechanics and Arctic Engineering,2002.359-367.

責(zé)任編輯陳呈超

Numerical Simulation of Wave Run-up on Cylindrical Offshore Structures

TANG Peng1, YU Ding-Yong1, Bai Wei2, TIAN Yan1

(1.College of Engineering, Ocean University of China, Qingdao 266100, China； 2.National University of Singapore, 119077, Singapore)

Wave run-up and wave impact can cause unexpected damage to offshore structures. Therefore, the design of offshore structures requires accurate predictions of the maximum wave elevation to maintain sufficient airgap below the platform deck. Accurate prediction of wave run-up can both help reducing building costs and avoid the risk of wave impact and damage to the platform. For the increasing number of offshore platforms built for ocean oil and gas exploration, the investigation of wave run-up becomes more and more significant for the design of fixed offshore structures. In this study, the interaction between wave and a surface-piercing vertical cylinder is investigated numerically based on waveFoam solver, which is developed under the framework of the OpenFOAM, an open source code library writing with fully objectoriented C++ language. The governing equations are discretized by using the finite volume method. Two fluids are considered and the volume of fluids (VOF) technique is employed to capture water-air interface. A 3D numerical wave tank of the viscous fluid theory was created based on Navier-Stokes equations and continuity equations, used by VOF method to trace free surface and relaxation wave method to generate wave. Meanwhile make use of the 3D wave tank to simulate wave run-up on a vertical cylinder. The model is validated by compare the extreme wave profile with analytical results. Then the simulation of extreme wave-cylinder interactions is carried out. Two different wave parameters was adopted, wave steepness parameters were discussed for the wave effect and run-up around the vertical cylinder. Compared with the potential flow model, viscous flow model used in this paper can better capture the secondary crest phenomenon appeared in the experiments. Calculations show that the numerical tank established in this paper can be better used to solve the problem of wave run-up on cylindrical offshore structures.

wave run-up; vertical cylinder; numerical wave tank; OpenFOAM

國家留學(xué)基金委項目(CSC 201306330022)資助

2014-11-12；

2015-02-15

唐鵬(1985-)，男，博士生。tangpeng851008@163.com

TV 139

A

1672-5174(2016)10-116-07

10.16441/j.cnki.hdxb.20140342

Suppoted by China Scholarship Council(CSC 201306330022)