弧形防浪墻上波浪荷載的數值模擬*

雷雪蓮，那鑫宇，孫 亮，周銀軍

(1.武漢理工大學 船海與能源動力工程學院，湖北 武漢 430063；2.中信建筑設計研究總院有限公司，湖北 武漢 430010；3.長江科學院 河流研究所，湖北 武漢 430010)

波浪運動到近岸時會對沿海的堤岸、碼頭、海塘等產生較大的沖擊，造成各類海岸建筑物不同程度的損壞。護岸結構及其頂部的防浪墻作為一種常見的海岸防護措施，可以防止海浪對陸地侵蝕、保障沿岸建筑物的安全。弧形防浪墻是對傳統(tǒng)直立式防浪墻的改進，能夠有效減少胸墻上方的越浪量、降低堤頂高程和工程造價。作為一種新型的結構形式，弧形防浪墻可以利用弧形面的引導將大部分的波浪水體導回海洋，對其波浪沖擊過程的研究可以通過物理模型試驗[1]和數值模擬[2]兩種方法進行。近年來越來越多的研究人員使用開源計算流體力學工具OpenFOAM[3]建立數值波浪水槽分析結構物所受到的波浪荷載。

OpenFOAM(Open Field Operation and Manipulation)是一個完全基于C++編寫的計算連續(xù)介質力學的自由開源程序庫。1993年，以Weller和Jasak為代表的帝國理工大學學者開發(fā)了OpenFOAM的前身FOAM。2004年Weller將FOAM的所有源代碼開放，并更名為OpenFOAM。OpenFOAM中的多相流求解器基于有限體積法對計算區(qū)域進行離散，使用流體體積法(volume of fluid)捕捉自由液面，同時可以根據水體流態(tài)的變化采用層流模型或者適當的湍流模型。

當雷諾數較小時，水的運動表現為層流，數值分析時可以采用基于層流的數值模型。Chen等[4]使用OpenFOAM分析了固定圓柱與波浪的非線性相互作用，將基于層流模型的數值計算結果與在丹麥水動力研究所進行的物理模型試驗進行比較，驗證了數值結果的精度。Chen等[5]基于層流模型分析了漂浮方箱與波浪的非線性相互作用，分析表明在共振頻率附近基于勢流理論計算所得到的方箱的轉動角度會大大高于試驗的測量值，而基于黏性流理論的數值結果與試驗的測量值吻合良好。孫亮等[6]基于OpenFOAM建立了數值波浪水槽，使用層流模型模擬規(guī)則波在潛堤地形上的傳播變形，將數值結果分別與代爾夫特理工大學的試驗數據、基于Boussinesq方程的數值結果對比，證明使用OpenFOAM可以精確地模擬波浪在潛堤地形上的傳播變形。

當雷諾數超過臨界范圍時，水的流態(tài)會從層流轉變?yōu)橥牧鳎谑褂肙penFOAM進行數值分析時需要引入適當的湍流模型。Jacobsen等[7]在原有OpenFOAM平臺的框架內加入了造波和消波功能建立數值波浪水槽，同時對原有的k-ω湍流模型進行了修正，通過分析波浪在斜坡上的傳播變形并與試驗測量結果比較，驗證了所開發(fā)數值波浪水槽的可靠性。Devolder等[8]提出在使用OpenFOAM中原有的k-ωSST湍流模型進行數值分析時，所得到的水與空氣交界面處的湍流黏度會遠大于實際值，從而導致波浪在數值水槽的傳播過程中發(fā)生顯著衰減；在湍流模型中加入了浮力修正項使得數值水槽中的波高保持穩(wěn)定，并通過分析單樁式風機基礎附近的波浪爬高驗證了修正的湍流模型的有效性。Devolder等[9]在所建立的數值波浪水槽中采用浮力修正的k-ω和k-ωSST湍流模型，模擬了波浪在斜坡上的破碎過程，并與Ting等[10]的試驗測量數據進行對比，再一次證明了在OpenFOAM的原有湍流模型引入浮力修正項的必要性。

為了分析圓弧半徑的大小對弧形防浪墻迎浪面所受波浪載荷的影響，李雪艷[11]在大連理工大學海岸和近海工程國家重點實驗室的溢流水槽對直立堤弧形防浪墻所受到的波浪作用進行了模型試驗研究。本文在OpenFOAM框架內采用Jacobsen等開發(fā)的造波、消波方法建立數值波浪水槽，分析波浪對弧形防浪墻的沖擊，并與李雪艷得到的試驗數據進行比較，驗證數值結果的精度。數值模擬中分別基于層流和湍流兩種假設，湍流模型選用了Devolder等開發(fā)的浮力修正的k-ωSST湍流模型。

1 數學模型

1.1 控制方程

數值分析中采用有限體積法求解雷諾平均N-S方程(RANS方程)。對于不可壓縮流體，RANS方程由連續(xù)性方程和動量守恒方程組成：

(1)

(2)

式中：t為時間；ui、uj(i，j=1，2，3指代x，y，z向)為流體的速度；ρ為流體密度；μeff為有效動態(tài)黏度，與計算中流體的體積分數α和湍流動力黏性ρνt相關；p*為超過流體靜力學的壓力；Fb為外部體積力(包括重力)；fσ為表面張力(在實際計算中可以忽略)。

基于RANS方程分析波浪與結構物相互作用時，將水和空氣看作一種可變密度的單一流體，在整個計算域內實現速度場、壓力場的共享。水和空氣的界面由基于體積分數α的流體體積法計算獲得，α=0代表單元完全被空氣充滿，α=1代表單元完全被水充滿，α在0～1表示單元中同時存在水和空氣。

1.2 湍流模型

OpenFOAM有多種湍流模型可以直接調用，主要是通過輸運方程求解湍流運動黏性νt。已有的文獻表明k-ωSST湍流模型在分析水流與結構物相互作用時具有較好的適用性，Devolder等以原有的k-ωSST湍流模型為基礎，考慮水和空氣的密度差異同時增加了浮力修正項Gb。Devolder等最終所建立的浮力修正的k-ωSST湍流模型為：

ρPk-ρβ*ωk+Gb

(3)

(4)

Pk=min(G，10β*kω)

(5)

(6)

(7)

(8)

式中：k為湍動能；ω為耗散率；v為流體的運動黏度；S為流體的平均應變率；F1和F2為計算σω、σk和γ所用到的函數；β*=0.09；α1=0.31；其他變量的介紹與計算方法見文獻[9]。

2 數值計算結果

2.1 數值水槽的建立

物理模型試驗中水槽長23.0 m、寬0.8 m、高0.8 m，試驗水深為0.4 m。直立堤弧形防浪墻試驗中設計了3種模型(圖1)，圓弧半徑R分別為0.45、0.67、0.98 m，靜水面以上弧形高度為0.32 m。試驗中在弧形防浪墻迎浪面中心線處安裝了壓力傳感器，3個模型測點布置相同。圓弧半徑R=0.67 m時圓弧上1#～4#號測點的位置見圖2。對表1所列的4種工況下弧形防浪墻受到的波浪沖擊過程進行模擬。工況2的波陡最大(H/L=0.082)，可以認為波浪的非線性最強、湍流效應的影響更為明顯。

圖1 物理模型試驗中的直立堤弧形防浪墻(單位：m)

圖2 R=0.67 m時弧形防浪墻上壓力測點布置

表1 波浪參數

根據直立堤弧形防浪墻的物理模型(圖1)和表1不同工況組合建立相應的二維數值波浪水槽。如圖3所示，對于工況1，波浪水槽長度設置為7.5 m、高度為0.72 m、水槽寬度為0.01 m。水槽左端邊界使用Jacobsen等所提出的方法產生入射波，水槽底部和右端邊界設置為壁面邊界，水槽頂部邊界設置為空氣自由出流邊界。由于目前所分析的問題為二維問題，水槽的前后邊界設置為空。為了避免入射波與右端壁面邊界相遇后產生的反射波作用到左端邊界，在水槽左邊界之前設置了4.0 m長的松弛區(qū)。對于表1中的其他工況，水槽高度保持不變，水槽長度取4倍波長，松弛區(qū)為2倍波長。由于波浪沖擊到弧形防浪墻時流體的運動是非定常的，數值模擬時需要設定初始條件。設置初始水位為0.4 m，計算域中初始速度為0。在計算過程中開啟可調節(jié)時間步功能，并且控制整個計算域中最大庫朗數小于0.25。

圖3 弧形防浪墻數值波浪水槽(單位：m)

2.2 網格收斂性分析及測點壓力計算

對于目前所建立的二維數值模型，沿水槽寬度方向只需要劃分1個網格。根據Chen等提出的建議，水平方向上每個波長范圍內至少需要劃分70個網格，垂直方向上在靜水面附近每個波高范圍內至少需要劃分8個網格。對于圖3所示的數值水槽，在水平方向上劃分了750個網格，平均尺寸為0.01 m。為了更準確地模擬波浪，沿著圖4所示方向由上、下兩邊向靜水面位置進行了漸變加密，上、下邊界的網格尺寸在垂直方向的高度是水面附近的2倍。計算中保持水平網格數量不變，在水深方向上采用了3套網格。網格1在垂直方向有36個單元，計算域內共有2.7萬個網格；網格2在垂直方向劃分了72個單元，總計5.4萬個網格；網格3在垂直方向劃分了144單元，總計10.8萬個網格。

圖4 弧形防浪墻附近的網格劃分

選用時間過程線完整的工況1的物理模型試驗數據進行網格收斂性分析?；趯恿骷僭O計算所得到的工況1中弧形防浪墻(R=0.45 m)上2#和3#測點處的壓力隨時間的變化見圖5，可以發(fā)現基于3種網格的數值計算結果基本一致。在保證計算精度的前提下，為了節(jié)約計算時間，在后面的分析中均采用與網格1相同的單元劃分。

圖5 工況1中弧形防浪墻(R=0.45 m)2#和3#測點處的壓力變化過程

分析時所采用的網格確定后，進一步將基于層流假設和湍流假設的計算結果與李雪艷在物理模型試驗所測量的數據進行對比。工況1中弧形防浪墻(R=0.45 m)上2#和3#測點處的壓力變化見圖6。圖6中點為試驗數據、實線為基于層流模型的數值模擬結果、虛線為采用浮力修正的k-ωSST湍流模型的數值結果。試驗中2#測點處壓力峰值在0.42～0.62 kPa波動，平均值約為0.55 kPa?；趯恿髂Ｐ湍M的壓力峰值在0.51～0.69 kPa波動，平均值約0.57 kPa，與試驗數據誤差約為3.6%?；诟×π拚膋-ωSST湍流模型的壓力峰值在0.45 ～ 0.58 kPa波動，平均值為0.52 kPa，與試驗數據誤差約為5.4%。試驗中3#處壓力峰值在0.19～0.34 kPa波動，平均值約為0.28 kPa?；趯恿髂Ｐ湍M的壓力峰值在0.29～0.45 kPa波動，平均值約0.34 kPa，與試驗數據誤差約21%。基于浮力修正的k-ωSST湍流模型模擬的壓力峰值在0.24～0.33 kPa，平均值約為0.28 kPa，與試驗數據吻合較好。

圖6 基于層流和湍流模型的壓力計算結果

2.3 弧形防浪墻上的壓力分布

圖7～9分別為工況2～4下不同圓弧半徑時弧形防浪墻上1#～4#測點處壓力峰值的平均值。由于在物理模型試驗中沒有指定計算峰值平均值時所選擇的周期數，也在一定程度上影響了數值結果與試驗數據之間的對比。

圖7 工況2下不同測點處壓力峰值的平均值

圖8 工況3下不同測點處壓力峰值的平均值

圖9 工況4下不同測點處壓力峰值的平均值

由圖7～9可以觀察到，不同圓弧半徑時防浪墻上的壓力都隨測點高度的增加非線性減小。表2為相鄰兩測點的壓力差值與高度差的比值，可以看到隨著測點位置的升高比值減小。

表2 不同測點壓力峰值的非線性減小

在工況2中基于浮力修正的k-ωSST湍流模型計算所得到的壓力值與基于層流模型的計算結果有顯著差異，同一測點的計算值層流模型結果比湍流模型高0.044～0.126 kPa；在工況3和工況4中，浮力修正的k-ωSST湍流模型模擬的結果與層流模型相比差別并不顯著，僅相差0.021～ 0.055 kPa。表3為計算所得到的壓力平均峰值的相對誤差，即數值結果與試驗數據的差值與試驗數據之比。湍流模型的誤差大部分集中在60%內，2組計算值大于100%，最大誤差152.86%；而層流模型則有7組計算值誤差大于100%，其中3組大于200%，最大誤差309.15%?？傮w來說，基于目前湍流模型的計算結果更接近試驗中的測量值。

表3 不同工況壓力峰值計算結果的相對誤差

2.4 弧形防浪墻所受到的波浪荷載

由于防浪墻半徑R=0.45 m在工況1中的試驗數據較完整，因此將此組合下計算所得波浪荷載的歷時曲線與試驗數據進行比較。圖10a)中模型試驗測得的水平波浪荷載峰值在50～65 N波動，平均值為55.76 N?；诟×π拚膋-ωSST湍流模型得到的水平荷載峰值在50～70 N波動，平均值為57.07 N，與基于層流的計算結果相比波動較小，更接近于試驗中的測量值。圖10b)中基于層流模型所得到的垂向波荷載峰值在15～27 N波動，平均值為19.82 N，與試驗數據的誤差為32.78%?；诟×π拚膋-ωSST湍流模型所得到的垂向波載荷峰值在13～18 N，平均值約為15.04 N，與試驗數據的誤差約為2.54%。波浪荷載計算結果的對比再次表明，基于浮力修正的k-ωSST湍流模型能夠更準確地模擬波浪對弧形防浪墻的沖擊過程。

圖10 R=0.45 m在工況1中弧形防浪墻的波浪荷載

3 結論

1)基于開源軟件OpenFOAM建立數值波浪水槽模擬了波浪對弧形防浪墻的沖擊過程?；趯恿髂Ｐ团c湍流模型的數值結果存在一定差異，尤其對于波陡較大的工況；基于湍流模型的數值結果更接近試驗數據，可為分析波浪與海岸結構物相互作用的強非線性問題提供可靠的依據。

2)在不同波浪要素條件下，直立堤弧形防浪墻迎浪面所受到的波浪力隨著圓弧半徑的增大而減小，不同圓弧半徑防浪墻上的壓力隨高度的增加非線性減小，實際工程中可以根據這一特性對防浪墻的構造形式進行優(yōu)化。