王藝之，韓新宇，董 勝

(中國海洋大學(xué) 工程學(xué)院，山東 青島 266100)

引 言

在海岸工程中，防波堤發(fā)揮著維護港池穩(wěn)定、保障港內(nèi)生產(chǎn)工作平穩(wěn)進行、保護生態(tài)岸線等重要作用，因此對防波堤迎浪面所受的波壓力也得到了廣泛的研究。防波堤一般分為混合堤、直立堤、特殊型式三種，其中直立堤是建立最早、應(yīng)用最廣泛的防波堤型式。

1963年侯穆堂和李玉成[1]通過物理試驗對直立堤上的波壓力進行了研究，他們將直立堤前的波浪分為立波、破碎波、過渡波三種形態(tài)，并對其中破碎型波壓力的計算進行了研究。隨后他們又研究了堤前破碎波的特征及其對混合堤和直立堤的作用[2]，提出破碎波壓力分布呈上大下小，且在靜水面處產(chǎn)生最大波壓力。Kirkg?z[3-6]從1992年到2005年期間對直立堤和斜坡堤上波浪破碎的荷載進行了一系列的研究，分析了堤前相對水深、堤前波浪形態(tài)、等因素對波浪作用力的影響，并提出了一種計算直立堤上破碎壓強的理論方法。他還提出，當垂直波面沖擊到混合堤上時，將造成最大波浪荷載。Yu等[7]研究了直立堤在斜向波作用下的水力學(xué)特性以及堤上波浪力的分布規(guī)律，并分析了波浪方向、波陡以及波浪的不規(guī)則性對直立堤所受波浪力的影響。

隨著科技發(fā)展，數(shù)值模擬由于具有節(jié)省人力物力、計算效率高、后處理方便等優(yōu)勢，得到了迅速發(fā)展，產(chǎn)生了多種多樣的求解方法，并在波浪與海工建筑物相互作用方面也得到了廣泛的應(yīng)用：黃蕙等[8]利用加源項造波的方法探究了波浪與涵洞式削角直立堤的越浪、透浪情況；李東洋等[9]基于OpenFoam開源軟件模擬了不規(guī)則波作用下扭王字塊護面斜坡堤的越浪過程；李昌良等[10]利用Ansys-Fluent軟件模擬了規(guī)則波與透空防波堤的相互作用，并討論了防波堤和波浪要素的變化對透空防波堤透射系數(shù)的影響；謝海清等[11]利用FLOW-3D軟件模擬了狹窄庫區(qū)河道滑坡涌浪產(chǎn)生及傳播過程，并對最大涌浪高度及其沿程傳播衰減規(guī)律進行了公式擬合；李晶晶等[12]利用格子Boltzmann方法對土石混合體的滲流場進行了模擬，并對其滲流特性進行了分析。

在利用數(shù)值模擬方法求解CFD問題時，常用的方法可以分為網(wǎng)格法和無網(wǎng)格法兩大類，其中，網(wǎng)格法主要包括有限差分法、有限體積法、有限元以及邊界元等，無網(wǎng)格法主要包括光滑粒子法(SPH)以及格子Boltzmann方法(LBM)等。本文使用的模擬方法，分別選取利用網(wǎng)格法中有限差分法求解為基礎(chǔ)的FLOW-3D，以及利用無網(wǎng)格法中LBM為基礎(chǔ)的XFLOW軟件進行數(shù)值模擬，探討這兩種軟件在模擬波浪對直立堤波壓力問題上的準確性，并選擇其中較為準確的方法進一步模擬，討論直立堤上波浪力隨入射波周期增大時的變化。在進行數(shù)值建模時，兩種模型均采用推波板造波、斜坡消波，并在防波堤附近進行網(wǎng)格加密或粒子細化。本文分別對兩種水位下、波峰和波谷分別作用時的波壓力進行了分析，并與于龍基論文中提供的實驗數(shù)據(jù)[13]進行了對比，所得結(jié)論對防波堤設(shè)計有一定參考價值。

1 數(shù)值模型

基于無網(wǎng)格法中LBM的數(shù)值模擬軟件XFLOW，其控制方程為BGK逼近下的格子Boltzmann方程：

(1)

式中，f=f(x,ξ,t)是單粒子分布函數(shù);ξ是微觀速度;λ是與碰撞有關(guān)的松弛時間;feq是平衡分布函數(shù)。

基于網(wǎng)格法中有限差分法的數(shù)值模擬軟件FLOW-3D，使用FAVOR網(wǎng)格處理技術(shù)以及VOF自由表面捕捉方法，達到精確定義幾何外形以及自由液面模型的效果。其控制方程為如下：

連續(xù)性方程：

(2)

N-S方程：

(3)

(4)

(5)

式中，u、v、w分別是x、y、z方向的速度；Ax、Ay、Az分別是三個方向的流動面積分數(shù)；Gx、Gy、Gz分別是三個方向的重力加速度分量；fx、fy、fz分別是三個方向?qū)?yīng)的粘滯力加速度；VF為體積分數(shù)；ρ為流體密度。具體原理可以分別參考文獻[14,15]，此處不再贅述。

2 數(shù)值模擬

2.1 數(shù)值模型設(shè)定

2.1.1 XFLOW數(shù)值模型設(shè)定

在XFLOW中建立二維水槽，長度為20 m，高度為1.2 m，防波堤與推板中心之間的水平距離為5.0L(L代表入射波波長)，在水槽末端設(shè)置消浪斜坡以達到消浪的效果。為提高計算效率，本文選取粒子間距為2.0 cm，僅在自由表面上下總高度約一個波高的矩形范圍內(nèi)以及整個防波堤附近進行自適應(yīng)粒子細化至1.0 cm。模擬時長為30 s。

2.1.2 FLOW-3D數(shù)值模型設(shè)定

在FLOW-3D中建立長度為20 m、寬度為0.5 m、高度為1.2 m的三維水槽，并將其劃分為三段：推板至防波堤迎浪面前約1倍入射波長處(記為A點)為造波段；從A點至防波堤后1 m處(記為B點)為試驗段；從B點至水槽末端為消波段，消波段放置消浪斜坡以達到消浪的效果。為保證試驗結(jié)果的準確性，同時提高計算效率，對造波段和試驗段中波高范圍內(nèi)的網(wǎng)格，以及試驗段中壓力測點高度范圍內(nèi)的網(wǎng)格進行加密，由于縮尺后波高為182.5 mm，且波高范圍內(nèi)包含不少于10個網(wǎng)格，故加密范圍內(nèi)網(wǎng)格邊長不大于18.25 mm，其它部分網(wǎng)格邊長約為25 mm。防波堤與推波板的相對位置根據(jù)率波情況確定，放置點率波結(jié)果應(yīng)與目標波高基本吻合。模擬時長同樣為30 s。

2.2 工況設(shè)定

共設(shè)兩種水位，分別為極端高水位(水深22.62 m)和設(shè)計高水位(水深21.88 m)，波高為7.3 m，波周期為8.9 s。在進行數(shù)值模擬時，對原模型進行縮放，取長度比尺為λl=40，時間比尺為λT=401/2，壓強比尺為λp=40，壓力比尺為λF=403。斷面尺寸如圖1所示，(斷面尺寸按照原型尺度標注，單位：mm)。在防波堤迎浪面布置壓強測點，其中A1—A6位于防浪墻上(A2位于設(shè)計高水位，A3位于極端高水位)，C1—C5位于沉箱段。以直立堤迎浪面為y軸，防波堤底面為x軸，測點坐標如表1所示(坐標尺度為縮尺后的尺度，單位：mm)。對應(yīng)的物理模型在長50 m、寬1.2 m、深1.2 m的水槽中進行，模型按照重力相似準則進行設(shè)計，模型縮放比例與數(shù)值模擬所用比尺相同，用于對比的波壓力試驗數(shù)據(jù)來自于文獻[13]。

圖1 防波堤斷面及測點示意圖

表1 測點坐標

在研究波壓力隨入射波周期變化情況時，增設(shè)110%T、120%T、130%T和140%T四種周期，與兩種水位共組合出十種工況進行計算。

2.3 數(shù)值造波及驗證

根據(jù)文獻[16]中的波浪理論的適用范圍，本文波浪采用斯托克斯二階波進行模擬比較合適。Madsen于1971年提出斯托克斯波的推板運動公式[17]，對于二階波，波高、波長和水深應(yīng)滿足HL3/d3<8π2/3的條件，推板的運動方程可以表示為：

(6)

相應(yīng)的推板速度方程為：

(7)

以極端高水位為例分別對數(shù)值水槽進行率波，驗證造波效果。得到的對比結(jié)果如圖2所示。兩種方法得到的波形和波高都能較好地與理論波形吻合。

圖2 數(shù)值造波與理論波形對比

3 基于數(shù)值結(jié)果的分析和討論

3.1 迎浪面壓力與試驗結(jié)果對比

圖3給出了分別在波峰和波谷作用下兩種不同的數(shù)值模擬方法得到的各測點波壓力與試驗結(jié)果[13]的對比。橫坐標為波壓力，縱坐標為測點坐標，橫、縱坐標均還原為原型尺度。從圖中可以看出，波峰作用時，最大波壓力發(fā)生在靜水位附近，在防波堤設(shè)計時應(yīng)考慮對靜水位附近進行加固。總體來說，兩種方法得到的壓力曲線均與試驗結(jié)果趨勢一致，吻合較好。

圖3 兩種數(shù)值模擬結(jié)果與試驗結(jié)果對比

從計算過程來看，兩種模擬方法各有利弊?；跓o網(wǎng)格法的XFLOW的優(yōu)勢是可以并行計算，且測點位置可以在后處理時設(shè)定，只需計算一次即可得到該參數(shù)下任意測點的結(jié)果，但其計算時間較長，無法較快地得到反饋并對模型進行調(diào)整；基于網(wǎng)格法的FLOW-3D的優(yōu)勢是計算速度快，能夠比較快地得到反饋，及時調(diào)整模型，但它需要在計算之前對測點地位置進行設(shè)定，只有設(shè)定好的測點能夠得到結(jié)果，一旦該測點需要調(diào)整位置就需要對整個模型重新計算。兩種方法共同的優(yōu)勢是可以對特定區(qū)域的網(wǎng)格進行加密或者自適應(yīng)細化粒子，顯著提高了計算效率。

從計算結(jié)果來看，分別對四種情況下防波堤整個迎浪面的波壓力進行積分，并計算其與試驗結(jié)果的誤差，具體誤差如表2。從表中可以看出，針對本文的模型，F(xiàn)LOW-3D數(shù)值模型的誤差除設(shè)計高水位下波峰作用時略大于X-FLOW數(shù)值模型以外，其它工況下的誤差均小于X-FLOW數(shù)值模型，且整體的平均誤差也相對較小。因此對于本文模型，F(xiàn)LOW-3D能夠得到更好的模擬結(jié)果，故在研究波壓力隨入射波周期變化情況時，選用FLOW-3D數(shù)值模型進行計算。

表2 兩種模擬方法所得波壓力的誤差

3.2 迎浪面壓力隨周期變化情況

圖4給出了沿高程分布的測點的波壓力分別在波峰、波谷作用下，隨周期變化的情況。從圖中可以看出，隨著入射波周期的增加，波峰作用時壓力最大點出現(xiàn)在靜水位附近，各測點波壓力均有增長，增幅均勻減小，極端高水位下當周期達到120%T后以及設(shè)計高水位下周期達到130%T后，各測點波壓力隨周期變化增幅極小；波谷作用時，C2測點以上各測點露出水面，故壓力基本一致，未發(fā)生變化，C2測點以下各測點壓力均有所減小，且從110%T開始各測點波壓力較為穩(wěn)定，隨周期變化的幅度極小。

圖4 迎浪面波壓力隨波周期的變化

4 結(jié)論

本文分別基于網(wǎng)格法和無網(wǎng)格法建立了數(shù)值水槽，模擬了規(guī)則波與直立堤的相互作用，以及波壓力隨入射波周期的變化，結(jié)論如下：

(1)兩種模擬方法都可以較好地模擬二階波對直立堤的作用力，沿高程的波壓分布曲線趨勢與實驗結(jié)果吻合良好。

(2)針對本文的計算模型，從計算過程來說，兩種方法存在差異；綜合幾種工況的平均誤差來看，F(xiàn)LOW-3D模擬結(jié)果的誤差更小。

(3)在其他波要素不變的情況下，隨入射波周期增大，波峰作用時，波壓力逐漸增大，且增幅均勻減??；在極端高水位時波浪周期T達到120%，或者設(shè)計高水位時波浪周期T達到130%的情況下，波壓力基本穩(wěn)定；波谷作用時，波壓力先減小，當周期T達到110%后，波壓力基本穩(wěn)定。