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

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

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

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

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

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

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

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

1.1 控制方程

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

(1)

(2)

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

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

1.2 湍流模型

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

ρPk-ρβ*ωk+Gb

(3)

(4)

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

(5)

(6)

(7)

(8)

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

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

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

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

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

圖2 R=0.67 m時(shí)弧形防浪墻上壓力測(cè)點(diǎn)布置

表1 波浪參數(shù)

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

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

2.2 網(wǎng)格收斂性分析及測(cè)點(diǎn)壓力計(jì)算

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

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

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

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

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

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

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

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

圖7 工況2下不同測(cè)點(diǎn)處壓力峰值的平均值

圖8 工況3下不同測(cè)點(diǎn)處壓力峰值的平均值

圖9 工況4下不同測(cè)點(diǎn)處壓力峰值的平均值

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

表2 不同測(cè)點(diǎn)壓力峰值的非線性減小

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

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

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

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

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

3 結(jié)論

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

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