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

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

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

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

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

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

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

1 數(shù)學(xué)模型

1.1 控制方程

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

(1)

(2)

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

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

1.2 湍流模型

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

ρPk-ρβ*ωk+Gb

(3)

(4)

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

(5)

(6)

(7)

(8)

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

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

2.1 數(shù)值水槽的建立

物理模型試驗中水槽長23.0 m、寬0.8 m、高0.8 m，試驗水深為0.4 m。直立堤弧形防浪墻試驗中設(shè)計了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)，可以認為波浪的非線性最強、湍流效應(yīng)的影響更為明顯。

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

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

表1 波浪參數(shù)

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

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

2.2 網(wǎng)格收斂性分析及測點壓力計算

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

圖4 弧形防浪墻附近的網(wǎng)格劃分

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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