非淹沒雙直立圓柱系統(tǒng)波浪爬升的數(shù)值模擬

陳浩民，倪云林

(浙江海洋大學(xué) 海洋工程裝備學(xué)院，浙江 舟山 316022)

非淹沒直立圓柱系統(tǒng)是近海構(gòu)筑物和海洋平臺常見的結(jié)構(gòu)，當(dāng)波浪在傳播過程中遇到這些柱體時，就會產(chǎn)生反射、繞射和爬升。隨著波浪沿柱體表面的迅速攀升，會對海洋構(gòu)筑物產(chǎn)生撞擊作用，從而影響其承載力和穩(wěn)定性。關(guān)于波浪爬升問題，文獻(xiàn)[1]基于線性勢流理論，給出了規(guī)則波繞射單個直立圓柱的解析解；文獻(xiàn)[2]通過物理模型試驗(yàn)，研究了大波陡波浪作用下直立圓柱的波浪爬升問題；文獻(xiàn)[3]對規(guī)則波和隨機(jī)波浪中的圓柱爬升問題進(jìn)行了模型試驗(yàn)研究；Morris-Thomas[4]采用物理模型試驗(yàn)研究了波陡和散射參數(shù)對波浪爬升的影響。近年來，隨著計算機(jī)技術(shù)的快速發(fā)展，基于三維數(shù)值水池的CFD方法被應(yīng)用于波浪爬升問題的研究[5-6]。Morgan等[7]采用Open FOAM求解器，探討了網(wǎng)格大小、離散格式、時間步長等因素對非線性波與圓柱相互作用的影響；武昕竹等[8-10]通過求解Navier-Stokes方程，結(jié)合VOF模型捕捉自由液面，計算了圓柱周圍的波浪爬升效應(yīng)；吳昊等[11]基于雷諾時均Navier-Stokes方程（RANS）和連續(xù)方程，結(jié)合k-湍流模型，建立了波浪與非淹沒豎直圓柱相互作用的三維數(shù)值模型。

在研究波浪與海上建筑物相互作用時，緩坡方程憑借其精確的計算性和較小的計算量而被廣泛應(yīng)用。緩坡方程即聯(lián)合折射繞射方程，最初由Berkhoff[12]基于線性波浪理論，在緩坡假定的基礎(chǔ)上，采用小參數(shù)展開的方法推導(dǎo)得到。Booij[13]就該方程的適用條件，證明在地形坡度小于1∶3時，緩坡方程具有足夠的精度。考慮到大部分近海構(gòu)筑物或海洋平臺由兩個以上的直立圓柱構(gòu)成，所以本文將采用有限元方法求解緩坡方程，研究波浪正向入射條件下，上、下游非淹沒雙直立圓柱周圍的波浪爬升特性。

1 數(shù)值模型

1.1 控制方程

Berkhoff[12]推導(dǎo)的緩坡方程為：

式中：g為重力加速度。

當(dāng)水深h為常數(shù)時，緩坡方程轉(zhuǎn)化為Helmholtz方程：

圖1 計算域劃分示意Fig. 1 Division of computational domain

1.2 求解方法

本文采用Houston等[14-15]列出的有限元法求解緩坡方程。如圖1所示，求解時將計算域分為水深變化的內(nèi)域Ω1和水深恒定的外域Ω2，BI為內(nèi)域中的物面邊界，滿足全反射條件，BC為外域和內(nèi)域之間的公共邊界，滿足速度連續(xù)和速度勢連續(xù)邊界條件，B∞為無窮遠(yuǎn)邊界，滿足輻射邊界條件。

外域不劃分網(wǎng)格，利用Helmholtz方程的級數(shù)解析解和輻射邊界條件，以Hankel函數(shù)Gi(i=1, 2, ···, M,M為截斷項(xiàng))作為權(quán)函數(shù)，通過格林公式將外域面積分轉(zhuǎn)為邊界上的線積分，建立單元有限元方程：

各單元有限元方程經(jīng)總裝形成總矩陣方程組后可采用列主元高斯消去法求解。

1.3 模型驗(yàn)證

為了驗(yàn)證本文數(shù)值模型的正確與否，首先計算非淹沒單個直立圓柱周圍的波高分布。如圖2所示，圓柱直徑d=0.522 m，淹沒深度h=0.165 m；波浪沿x軸正方向入射，波高H0=0.079 m，周期T=0.94 s，波長L=1.046 m。圖2中圓環(huán)部分為內(nèi)域，大小為用Gambit軟件劃分網(wǎng)格。在兼顧計算精度和計算速率的前提下，1個波長范圍內(nèi)設(shè)置15個三角形六節(jié)點(diǎn)網(wǎng)格，內(nèi)域總計9 098個單元，18 444個節(jié)點(diǎn)。

計算結(jié)果繪于圖3。可以看出，本文的數(shù)值解波浪爬高Hi與文獻(xiàn)[1]提出的線性解析解吻合很好，最大波高位于圓柱迎浪點(diǎn)（θ=0°）位置處，其值約為入射波高的1.8倍，最小波高位于圓柱背浪點(diǎn)（θ=180°）左右兩側(cè)35°位置處，即θ=145°，其值約為入射波高的50%。這證明了本文數(shù)值模型的正確性，其計算結(jié)果準(zhǔn)確可靠。

圖2 非淹沒單個直立圓柱地形（單位：m）Fig. 2 Sketch of single non-submerged vertical circular cylinder (unit: m)

圖3 單個直立圓柱周圍波高分布比較Fig. 3 Comparison of calculated wave run-up with analytical solutions along circumference

2 非淹沒雙直立圓柱系統(tǒng)波浪爬高的數(shù)值模擬

為了進(jìn)一步研究非淹沒雙直立圓柱系統(tǒng)中圓柱間距對波浪爬升的影響，如圖4所示，在沿波浪傳播方向上布置兩個直立圓柱，分別將其命名為上游圓柱和下游圓柱，圓柱尺寸和波浪要素與上文模型驗(yàn)證中的算例參數(shù)相同，兩圓柱間距為S。針對不同的工況，S是唯一變量，其變化范圍是1 . 25L≤S≤3.00L，變化步長為0.25L，計算工況1～8對應(yīng)的圓柱間距分別為1.25L、1.50L、1.75L、2.00L、2.25L、2.50L、2.75L和3.00L。

圖4 非淹沒雙直立圓柱地形（單位：m）Fig. 4 Sketch of dual non-submerged vertical circular cylinder(unit: m)

2.1 上游圓柱波浪爬升特性

不同間距條件下，上游圓柱周圍波高分布隨圓柱間距 S 的變化情況及其與單個直立圓柱的對比見圖5。

相較于單個圓柱的情況，波高分布曲線存在較大的波動現(xiàn)象，具體表現(xiàn)為：

（1）當(dāng)圓柱間距為1/4波長的整數(shù)倍時，即 S /L=n/4(n=5,6,···,12)時，最大波高的出現(xiàn)位置和最大波高值與單個圓柱的情況基本相同，均位于圓柱迎浪點(diǎn)（θ=0°）處，最大值約為入射波高的1.8倍。

圖5 圓柱間距改變對上游圓柱波浪爬高的影響Fig. 5 Influence of spacing between two cylinders on wave run-up of upstream cylinder

2.2 下游圓柱波浪爬升特性

圖6 是下游圓柱周圍波高分布隨圓柱間距S的變化情況及其與單個直立圓柱的對比?？梢钥闯?，在不同圓柱間距的工況下，下游圓柱周圍的波高分布曲線與單個直立圓柱的情況相似，最大波高均位于圓柱迎浪點(diǎn)（θ=0°），最小波高均位于圓柱背浪點(diǎn)左右兩側(cè) 35°處（θ≈145°），最大值和最小值計算結(jié)果詳見表1。

分析表1可見：

（1）受上游圓柱的掩護(hù)作用，下游圓柱周圍相對波高的最大值和最小值均小于單個圓柱的情況。

圖6 圓柱間距改變對下游圓柱波浪爬升的影響Fig. 6 Influence of spacing between two cylinders on wave run-up of downstream cylinder

表1 下游圓柱波浪爬高最大值和最小值Tab. 1 The maximum and minimum values of wave run-up around downstream cylinder

3 結(jié) 語

（1）本文采用有限元方法建立了緩坡方程數(shù)值模型，利用該模型計算了非淹沒單個直立圓柱周圍的波高分布，計算結(jié)果與解析解吻合很好，證明了本文數(shù)值模型的正確性以及用其求解圓柱周圍波浪爬升問題的適用性。

（2）相較于單個圓柱，上游圓柱周圍波高分布曲線的波動較大，但它們的最大相對波高和最小相對波高發(fā)生位置基本一致。同時，當(dāng)圓柱間距為1/4波長的奇數(shù)倍時，即 S /L=(2n+1)/4(n=2,3,4,5)，最小相對波高明顯小于單個圓柱的波高，并在圓柱肩部出現(xiàn)第二峰值；而當(dāng)圓柱間距為1/4波長的偶數(shù)倍時，即S/L=2n/4(n=3,4,5,6)，最小相對波高則明顯大于單個圓柱的波高，并在圓柱肩部出現(xiàn)第二谷值。這是由上游圓柱的繞射和下游圓柱的反射所引起。

（3）下游圓柱周圍的波高分布曲線與單個圓柱的情況相似，但受上游圓柱的掩護(hù)作用，其最大相對波高和最小相對波高均小于單個圓柱的情況。

（4）吳昊等[11]基于黏性流理論，研究了非淹沒雙直立圓柱系統(tǒng)的波浪爬高問題，結(jié)果表明，當(dāng)圓柱間距較大時，下游圓柱壅水效果降低，波浪爬升減弱；當(dāng)圓柱間距較小時，圓柱之間壅水效果增強(qiáng)，上游圓柱背浪面的波浪爬高逐漸上升。上述結(jié)論與本文采用勢流理論得到的結(jié)果基本一致。