二維間斷有限元水動力模型與波浪模型實時耦合研究

郭江華, 孔令臣，張慶河，王容基，冉國全

(1.天津大學 水利工程仿真與安全國家重點實驗室，天津 300072；2.中交天津港灣工程設計院有限公司，天津 300461；3.中冶賽迪重慶信息技術有限公司，重慶 400013)

間斷有限元(Discontinuous Galerkin，DG)方法具有高階精度和局部守恒性，對復雜地形有很好的適應性和易于并行[1]，近年來得到了越來越廣泛的應用[2-7]。就二維水動力模擬而言，Aizinger和Dawson[8]較早建立了間斷有限元模型，Kubatko等[9]將hp自適應間斷有限元模型應用于以對流為主的水體運動模擬，Brus等[10]將高階間斷有限元應用于海岸水動力學問題。針對河口海岸淺水流動，比利時的魯汶大學利用DG開發(fā)了SLIM非結構化網(wǎng)格二維和三維水動力模型，并應用于Fly河污染淤泥入侵Torres海峽的模擬[11]、Titicaca湖藻華模擬預測[12]、湄公河下游河流來沙對河口環(huán)境特性影響[13]等一系列研究，其中的二維模塊應用于珊瑚礁海域大范圍潮流模擬[14-16]。天津大學團隊提出了一種任意非結構化網(wǎng)格無積分DG方法[17]，并相繼開發(fā)了水動力、泥沙等數(shù)值模型[18-19]。值得指出的是，在近岸海洋環(huán)境條件下，波浪與水流耦合作用往往對波浪傳播和水流運動有著重要影響，因此目前基于有限體積、有限元或有限差分法建立的共享二維和三維水動力模型，如FVCOM、Delft3D、TELEMAC、ROMS、ADCIRC等均通過不同方法實現(xiàn)了實時波流耦合模擬[20-24]，而現(xiàn)有DG水動力模型尚缺乏針對波流耦合模擬的研究。

為此，本文將借鑒已有波流耦合模型的實現(xiàn)方法，在李文俊等[18]建立的二維非結構化網(wǎng)格無積分間斷有限元水動力模型的基礎上，采用我國清華大學自主開發(fā)的C-Coupler2耦合器與第三代波浪模型SWAN進行實時雙向耦合，建立波浪與節(jié)點間斷有限元二維水動力的實時波流耦合模型。

1 二維波流耦合模型

1.1 水動力模型及其DG離散

水流控制方程為如下淺水方程

(1)

(2)

式中：ρw為流體密度；z為底坡高程；τbx和τby分別為水底床面沿x和y方向的摩阻應力；DFx和DFy為波流共同作用下的水平擴散項。Sxx、Sxy、Syx、Syy為沿水深積分的波浪輻射應力分量，如式(3)所示

(3)

(4)

式中：Er為水滾動能。

方程的離散采用DG方法，首先將計算域劃分為Ne個不重疊單元，對于第k個單元Γk，滿足式(1)，在Γk上最高不超過p階的局部多項式空間xp(Γk)上選擇一組基函數(shù)Φk(x)，用以對精確解近似，使得殘差最小，殘差方程如下

(5)

(6)

對式(6)應用兩次分部積分及格林公式可得到如下方程形式

(7)

將解用基函數(shù)表示，經(jīng)過推導可得空間半離散形式

(8)

式中：Mk為單元總體質量矩陣；Dk,*為在*方向的微分矩陣；Me,k為單元內(nèi)每條邊對應的質量矩陣。最后采用顯式二階Runge-Kutta方法對時間進行離散。

1.2 波浪模型

波浪模型采用第三代海浪模型SWAN(Simulation WAves Nearshore)，可用于模擬海洋、河口、海岸和湖泊等水域中波浪的生成與傳播[26-27]。SWAN采用二維波作用譜密度來描述隨機波浪場，在笛卡爾坐標下模型控制方程如下

(9)

式中：左端第一項表示波作用量N隨時間變化，第二項和第三項表示波作用量以速度cx、cy在x和y方向上的傳播，第四項表示波作用量N在相對頻率σ空間上的變化，第五項表示水流及地形變化引起的波浪折射，右端項表示源項，包括風能輸入和波浪相互作用波浪破碎、底摩阻等引起的能量損耗。關于SWAN模型的數(shù)值離散，詳見文獻[27]。

1.3 波流耦合的實現(xiàn)

波流耦合利用C-Coupler2耦合器實現(xiàn)，耦合框架如圖1所示。C-Coupler2是由我國清華大學獨立自主研發(fā)的地球模式耦合器，C-Coupler2耦合器具有并行軟件架構，可生成適用不同耦合模式的耦合器實例[28-29]。在本耦合系統(tǒng)中，水動力模型和波浪模型均采用相同的非結構化網(wǎng)格，二者同時獨立運行，根據(jù)指定的耦合頻率通過C-Coupler2實現(xiàn)流場與波浪場物理量交換，其中，二維水動力模型向SWAN模型提供水位η、x和y方向上的垂向平均流速u和v，同時SWAN提供有效波高Hs、波向θ和破波率Qb等波浪要素。

圖1 耦合框架圖

2 模型驗證

波浪在近岸區(qū)域傳播時，隨著地形變化，會發(fā)生淺水變形、折射、繞射以及波浪破碎等波浪變形現(xiàn)象，在此過程中波高發(fā)生變化，并通過波浪輻射應力引起近岸波浪流現(xiàn)象。同時，近岸區(qū)水位變化與水流流速也會明顯引起波浪場的變化，波流耦合模型應能夠合理描述受水流影響的波浪近岸傳播變形及波浪影響下的近岸水流。為此下面通過系列已有實驗驗證波流耦合模型的合理性。

2.1 裂流通道波生流實驗

Hamm在水池中進行了裂流通道地形波生流實驗[30]，實驗在水池中進行，地形如圖2所示，中間有下凹的裂流通道，坡度為1:30，入射波浪采用方向譜，其表達式為S(f,θ)=Sγ(f)cosnθ，Sγ采用JONSWAP譜，γ取為3.3，方向分布參數(shù)n取6，入射波高0.07 m，譜峰周期1.25 s，破波指標γ=0.71。

圖2 裂流通道地形及波浪場

波浪在向岸傳播過程中，由于地形變化產(chǎn)生折射，并且產(chǎn)生波能輻聚和輻散，波浪傳播矢量場如圖2所示，與不耦合模型相比，耦合模型受水流作用明顯，在裂流通道附近波高增加。在破波帶內(nèi)，由于下凹地形影響，波浪增水在沿岸形成水位差異產(chǎn)生的靜壓力差形成沿岸流，沿岸流在此區(qū)域匯合，以裂流形式流向外海，如圖3所示。圖4顯示了通過裂流中軸線處(Y=15 m)的波高分布，其中實線和虛線分別是通過耦合和非耦合模型計算的波高沿程變化圖，可以看出，如果不考慮水流的耦合作用，單純波浪模型不能反映出裂流通道內(nèi)波浪的變化。實際上，由于裂流通道處波浪受到和波浪傳播方向相反的逆向水流影響，波高明顯增加，只有利用波流耦合模型才能獲得與實驗結果一致的波高分布。

圖4 裂流通道軸線波高沿程變化

2.2 T型堤波生流實驗

Gravens和Wang[31]在波浪水池中進行了T型防波堤掩護下的動床實驗，以研究斜向不規(guī)則波產(chǎn)生沿岸流作用下泥沙輸移。實驗布置如圖5所示，防波堤位于Y22與Y26之間區(qū)域，由四臺造波機進行造波，有效波高為0.22 m，周期1.5 s，波向角為6.5，破波指標γ=0.73。實驗過程中測量了斷面處的沿岸流流速與波高，這里采用靠近防波堤的Y22、Y26測點處的沿岸流流速與波高進行驗證。

模型計算結果與測量結果吻合良好，如圖6、圖7所示，說明模型能夠反映出正確的水動力規(guī)律和T型堤附近波浪變形規(guī)律。在T型堤掩護區(qū)，波浪以繞射形式傳入，掩護區(qū)內(nèi)外輻射應力梯度較大，于是在防波堤右側形成順時針環(huán)流，而左側由于防波堤阻擋沿岸流形成了逆時針環(huán)流，如圖8。因此本文波流耦合模型可以應用于復雜地形，計算結果合理且與實驗結果相符合。

圖6 Y22斷面沿岸流流速分布及斷面波高沿程分布

圖8 T型堤局部流場

3 模型應用

將模型應用于黃驊港大風過程中水動力特性研究。根據(jù)實測資料顯示，在2006年3月10日到3月14日期間，渤海區(qū)域內(nèi)出現(xiàn)較大風速，本文將對此時間段內(nèi)黃驊港附近海域波流動力特性進行研究。波浪和水動力計算均采用為大、小模型嵌套的方式，采用相同的三角形網(wǎng)格，大小模型網(wǎng)格劃分如圖9，其中大模型節(jié)點數(shù)為27 221，網(wǎng)格單元數(shù)為52 522；小模型節(jié)點數(shù)為17 652，網(wǎng)格單元數(shù)34 320，水動力計算模型采用2階精度。圖9～圖11中A、B分別為波浪測站與潮位測站，P0、P1、P2點分別為0 m、-2 m、-4 m等深線(以平均水面為基準)上的計算測點。波浪和水動力計算均考慮風場作用，采用ERA5風場數(shù)據(jù)，空間分辨率為0.25°，時間分辨率為1 h。

圖9 大、小模型網(wǎng)格劃分與各測點位置

根據(jù)潮位與波高驗證結果(圖10、圖11)，計算結果與實測潮位和波高吻合良好，說明模型能夠合理地模擬大風過程中黃驊港區(qū)域潮流與波浪的運動。

圖10 潮位驗證(以平均水面為基準，0時刻為2006-03-11 12:00)

由于波浪破碎過程中將伴隨能量損失，導致波高等波要素發(fā)生變化，從而在近岸區(qū)域出現(xiàn)明顯的增減水現(xiàn)象。圖12為P0、P1和P2點的水位歷時曲線，其中實線和虛線分別為耦合模型和未耦合模型計算結果，可以看出由于波浪作用，近岸地區(qū)的水位相對于只考慮風增水的水動力模型會明顯增加，即產(chǎn)生波浪增水，P0點(0 m等深線)增水將近0.2 m，P1點(-2 m等深線)最高增水0.1 m，并且隨著水深增加增水減小。

圖12 各點水位變化

圖13為點P0、P1和P2的波浪歷時曲線，其中實線和虛線分別為耦合模型和未耦合模型計算結果，結果顯示潮位變化和潮流對波浪的作用顯著，其中P1點(-2 m等深線)區(qū)域波高增大可達0.8 m左右，并且當水深越小時作用越明顯。

圖13 各點波高變化

總體而言，波流相互作用越近岸越明顯，隨著水深的增加而減弱。相較于非耦合模型，耦合模型的計算結果更為合理。

4 結論與展望

本文在李文俊等[18]建立的無積分節(jié)點間斷有限元二維水動力模型的基礎上，利用C-Coupler2耦合器與波浪模型SWAN進行實時耦合，建立了二維波流耦合間斷有限元模型。水動力計算中考慮了輻射應力和波面水滾等源項，波浪模型通過實時更新水位、流速等條件，實現(xiàn)了模型的實時雙向耦合。模型通過Hamm裂流通道實驗、Gravens和Wang的T型防波堤實驗進行驗證，數(shù)值模擬結果與實驗結果吻合良好，模型能夠反映出由于輻射應力梯度產(chǎn)生的波生沿岸流、由于沿岸不均勻地形產(chǎn)生的裂流和堤后環(huán)流、以及水流對波浪傳播的影響。將模型應用于黃驊港海域模擬計算，較好地描述了現(xiàn)場波浪和水動力變化過程。

波流作用下海岸泥沙運動和地形演變是海岸動力學關心的重要問題，今后將進一步在波流耦合模型的基礎上發(fā)展基于間斷有限元方法的波流耦合作用下泥沙運動和地形演變模型。