基于無積分節點間斷有限元的二維水沙模型：(1)水動力

李文俊，張慶河，李龍翔，冉國全

(天津大學 水利工程仿真與安全國家重點實驗室，天津 300072)

數值模擬方法在二維河流與海岸自由水面流動模擬中得到了廣泛應用，它可以合理預測水動力變化規律，并為模擬二維泥沙運動和地形演變、預測污染物輸移、指導工程建設[1]等工作提供基礎。過去幾十年間，二維水動力模型取得了相當多的成果，有限差分法、有限體積法和有限元法等都在求解二維淺水方程中獲得了大量應用。最近若干年來，間斷有限元方法也逐漸開始應用于二維淺水方程的數值模擬，一系列研究者的工作使得這一模型逐步具有了計算干濕、保證和諧特性、保證水深為正等特性[2-4]。這一模型也逐漸被應用到立波、潰壩以及大范圍潮流、海嘯、風暴潮等一系列復雜水動力模擬中[5-8]。

間斷有限元方法融合了有限元法和有限體積法的諸多優點，不僅可以通過提高單元插值基函數的階數來獲得高精度數值結果，而且可以對數值通量采用有限體積法中的Riemann算子精確計算，確保間斷捕捉的準確性和質量的守恒性。模型可采用非結構化網格，能夠很好地模擬復雜岸線，同時可應用多種自適應網格算法。另外，在求解給定單元內部的變量時，間斷有限元只需要相鄰單元的變量，使得處理器之間的信息傳遞量保持最小，易于實現算法的并行化[9]。當然，間斷有限元方法也存在未知數個數較多、計算量較大的缺點。針對這一問題，李龍翔和張慶河[10]建立了任意四邊形網格的無積分節點間斷有限元模型。本模型相比于已有間斷有限元淺水模型，將數值積分轉換為矩陣運算，省略傳統數值計算中求解積分節點處函數值的過程，減少計算量的同時減少了單元的內存儲量大小。同時本模型可采用任意四邊形網格計算，相比于三角形網格可以有效地減少未知量個數，從而減少計算量。

為了開發高精度高效間斷有限元二維水沙模型，本文將在李龍翔和張慶河[10]基于任意四邊形網格無積分節點間斷有限元方法求解對流擴散方程的模型基礎上，將控制方程擴展為二維淺水方程，并加入相關源項，建立間斷有限元二維水動力模型。該模型將考慮底坡源項、科氏力、風應力、底摩阻和地形干濕變化，以充分反映在各種實際條件影響下淺水流動的二維水動力特性。

1 二維水動力數學模型

1.1 控制方程

非恒定二維淺水方程的守恒形式為

(1)

式中：Q=[h,hu,hv]T為守恒變量，其中h為水深，u、v分別為沿x和y方向的水平流速。通量項E和G分別為x，y方向的通量，表達式為

源項S考慮了底坡、科氏力、底摩阻和表面風應力影響，最終表達式為

(2)

式中：η=h+z為水位，z為底坡高程，ρ為水體密度。

科氏力系數f可由式(3)給出[11]

f=2ωsinφ

(3)

式中：ω為地球繞地軸旋轉的角速度，大小為ω=7.29×10-5rad/s，φ為地理緯度。

τbx和τby分別為水底床面沿x和y方向的摩阻應力，τwx和τwy分別為自由表面沿x和y方向的風場切應力,分別由經驗公式(4)和(5)給出[12]

(4)

(5)

式中：n為糙率，Cw為無因次風應力系數，ρa為空氣密度，ωw為水體自由表面以上10 m處的風速，β為風向角，一般取風向與x方向的夾角。

上述方程寫成以分量形式表示的方程即

(6)

為了推導方便，下面的數值離散仍以矢量形式為主。

1.2 節點間斷有限元數值離散

利用在單元上連續可微、單元間允許存在跳躍的分段光滑函數Qh逼近Q[13]，Qh表達式為

(7)

在節點間斷有限元方法中，插值基函數與試驗函數ji(x)均采用Lagrange插值函數定義，定義分段光滑函數Qh所在的空間為Vh。利用Qh替換控制方程中的Q，并用試驗函數v∈Vh乘以替換后的公式，在單元上積分，可得

(8)

對上式用分部積分和格林公式可得

(9)

式中：F表示通過邊界的法向通量，n為邊界?Ωk處的單位外法向量，式(9)便是弱形式的空間離散方程組。由于間斷有限元方法中物理通量F在邊界兩側的值不相同，故將通過邊界的法向通量用數值通量F*代替[14]。利用二維淺水方程的旋轉不變性，將整體坐標系繞原點逆時針旋轉角度φ(φ為邊界外法線向量與x軸正方向的夾角)，使坐標軸x向與邊界外法線向量平行，并將原點平移至對應邊界中點處，從而建立局部坐標系。由此將數值通量的求解轉化為局部坐標系下關于守恒變量向量Q的一維Riemann問題的求解。本文采用HLL格式[15]的數值通量，其形式為

(10)

其中，QL和QR、EL和ER分別為局部坐標系下交界面左右的守恒變量、通量向量值

(11)

SL和SR為波速，分別由式(12)和式(13)計算得到[16]

(12)

(13)

其中，

(14)

(15)

(16)

對式(9)左端第二、第三項再運用一次分部積分，得到方程組

(17)

式(17)為強形式的空間離散方程組。相對于弱形式而言，強形式無積分離散格式的存儲量更小。

本文采用4階5步顯式Runge-Kutta(RK)格式對半離散方程(17)進行時間積分，最終可得到不同時刻具有高階時間與空間精度的數值解。

1.3 無積分格式

傳統間斷有限元方法在求解式(17)時采用數值積分，這種方法要求被積函數在積分節點處的值已知，因此需要先計算出積分節點處的守恒變量值與通量值，從而造成計算量較大而影響計算效率。為了提高計算效率，Atkins和Shu[17]曾提出一種無積分格式，通過假定單元內雅克比行列式為常數，將數值積分運算轉化為標準單元內矩陣計算。這一方法省去了數值積分過程，但其只適用于結構化平行四邊形網格，對于任意四邊形單元并不適用。李龍翔和張慶河[10]提出的無積分格式無需假定單元內雅克比行列式為常數，可適用于任意四邊形網格。該格式將積分節點和插值節點均選定為同一套邊界點的Legendre-Gauss-Lobatto(LGL)節點，在計算式(17)時將任意四邊形單元通過雅克比變換映射到標準四邊形單元，并將體積分中的單元剛度矩陣表示為質量矩陣與導數矩陣乘積形式，最終可得到標準單元坐標系下包含NP個方程的方程組(18)

(18)

式中：NP為單元插值節點個數，1≤j≤NP，r、s為標準單元坐標變量，J為雅克比變換矩陣行列式，由式(19)給出

(19)

無積分格式將數值積分轉化為上述矩陣的運算，避免了求解積分節點處的函數值。同時，確定標準單元后，其導數矩陣、體積分質量矩陣和面積分質量矩陣都為定值，可以預先調用矩陣運算庫計算，從而有效減少計算量并節省大量儲存空間。特別需要指出的是，當采用高階精度格式時，無積分運算求解數值積分產生的誤差與采用高階基函數近似產生的截斷誤差同階[18]。因此，無積分方法同樣對于格式計算誤差影響不大。

1.4 干濕處理及和諧特性

1.4.1 干濕處理

近岸二維水動力模型模擬實際問題的難點在于處理復雜的干濕邊界條件[19]。對于水位變化引起的岸線變化的處理，最主要的兩種方法是動邊界法和干濕算法。其中動邊界往往需要持續地對網格進行重新劃分，所需計算量偏大。因此，二維水動力模型大多采用干濕算法。干濕處理需要解決三個主要問題：(1)保持水深為正值；(2)消除非物理壓力梯度造成的虛假振蕩；(3)保持干濕交界處的數值穩定[20]。

1.4.2 和諧性

在守恒型淺水方程模型中，水位梯度項被分為壓力項和底坡源項之和

(20)

為了使控制方程離散后該等式始終成立，即滿足靜水條件下水位為常數、流速為0的和諧(Well-balanced)特性[24]，需要引入一些特殊處理手段，如：改變單元水力模型，消去底坡項；或者對底坡源項進行改造，使底坡源項與壓力項在靜水條件下能夠平衡[25]。本模型采用改造底坡源項的方法，將底坡源項作如下拆分

(21)

并采用Bollermann[26]的方法，將底坡方程用一個連續分段線性函數近似，以保證模型的和諧特性。

2 模型驗證

按照上述數值模式建立了二維水動力模型，下面將分別利用不同算例對模型進行驗證，其中包括：(1)利用赤道孤立Rossby波驗證模型描述淺水長波和科氏力的能力；(2)利用平底坡Stommel問題驗證模型描述風應力作用以及科氏力與底摩阻聯合作用的能力；(3)利用平緩岸灘波浪爬坡算例驗證模型的干濕處理能力；(4)利用斜坡靜水算例驗證模型的和諧特性。

2.1 赤道孤立Rossby波

本算例模擬的是赤道處孤立Rossby波的運動情況，算例的計算域取為[0,48]×[0,16]的矩形區域，初始的孤立Rossby波位于計算域中心，即坐標[24,8]處，其0階初始波形由式(22)給出[27]

(22)

式中：Γ(x)=0.771a2sech2(ax)，a為確定孤立波振幅的參數，取為0.395。

算例求解的是無量綱化的二維淺水方程，方程中源項除了底坡源項僅包括科氏力項。由于弱非線性和弱色散之間的平衡，該問題的攝動解預測赤道孤立Rossby波會保持其形態向西行進[28]。采用計算網格為4 181個任意四邊形單元。

1-a t=0 s1-b t=10 s

1-c t=20 s1-d t=40 s圖1 赤道孤立Rossby波算例模擬結果示意圖Fig.1 Simulation result of equatorial solitary Rossby wave

算例共計算了40個單位時長，多項式階數P=2，四周邊界均視為墻面，采用可滑移邊界條件，模擬結果如圖1所示。由圖可知，在孤立Rossby波向西傳播過程中，有一部分水體分離出來以孤立Kelvin波的形式向東傳播，這一現象是由于給定的初始條件不是精確解造成的。攝動解預測孤立Rossby波波峰相位速度為-0.77 m/s，最終會運動到x=-15.68處，波峰高度為0.153 4；Giraldo和Warburton運用間斷有限元方法模擬的孤立Rossby波波峰相位速度為-0.78 m/s，最終會運動到x=-15.56處，波峰高度為0.148 4；本模型計算結果孤立Rossby波波峰相位速度為-0.78 m/s，最終會運動到x=-15.64處，波峰高度為0.158 4。本文模型模擬結果與攝動解和已有數值模型的結果均接近。

為了說明提高多項式階數對計算精度的作用，分別計算了孤立Rossby波在多項式階數P分別為1、2、3時的結果，三個算例的計算網格數均為1 600個。由于本算例沒有精確解，本文選取25 600個網格下P=3的計算結果作為精確解，分別計算了不同階數下模擬所得水位的二范數誤差，如表1所示。由表可知，在相同的網格劃分條件下，隨著多項式階數的提高，模擬的誤差逐漸減小，模型的計算精度逐漸提高。

表1 不同階數下孤立Rossby波模擬結果誤差比較Tab.1 L2 errors of simulation results in equatorial solitary Rossby wave problem

2.2 平底坡Stommel問題

Stommel問題是科氏力、風應力和底摩阻達到平衡穩定狀態時的解，它涉及到淺水方程實際應用中會遇到的多種物理過程，是海洋環流問題的一個很好的原型[29]。參考文獻[29]的算例設置，算例中計算域大小為106×106m的矩形區域，計算域內水深均為1 000 m，四周邊界視為壁面，采用可滑移邊界條件。計算域共劃分為502個任意四邊形網格，算例選取多項式階數P=2。本算例中科氏力系數f由式(23)給出，風場切應力由式(24)給出，底摩阻應力由式(25)給出

(23)

(24)

τbx=ργuh，τby=ργvh

(25)

2-a 模擬結果水位等值線2-b 模擬結果流速等值線圖2 穩定狀態水深等值線圖及流速等值線圖Fig.2 The free surface elevation and the velocity magnitude at the steady state of the nonlinear Stommel problem

式中:f0=10-4，β=10-111/m，L=106m；τ0=0.2 N/m2；ρ=1 000 kg/m3，γ=2×10-6，各參數選取均與文獻[29]相同。

在科氏力、風應力和底摩阻的共同作用下，計算域內初始會在中心形成一個凸起的波，之后波峰慢慢向左側偏移，當偏移一定距離后達到穩定狀態，波的形狀和位置均不再發生變化。結果達到穩定時計算域內水位等值線和流速等值線如圖2所示，與文獻[29]中的結果相近。

2.3 平緩岸灘波浪爬坡

為了說明本模型的干濕處理能力，考慮一個由海底滑坡引起的波浪在平緩岸灘上爬坡的算例。本算例長度為[-500 m,50 000 m]，底坡坡度為1/10，初始水位由式(26)給出

η=a1e-k1(x-x1)2-a2e-k2(x-x2)2
u=0,t=0

(26)

(27)

3-a 初始波面3-b t=160 s、175 s、220 s時爬坡和回落過程波面圖3 波浪初始波面及模擬波面圖Fig.3 The initial wave profile and the simulation result of wave profile at t= 160 s, 175 s and 220 s

初始波形如圖3-a所示，初始速度為0。計算域劃分為10 000個網格，計算時間為280 s，多項式階數P=2。圖3-b給出了t=160 s、175 s和220 s三個時刻的波面模擬結果。由圖可知，本模型模擬的波浪爬坡(t=220 s)、回落(t=160 s、t=175 s)波面和Carrier等[30]給出的半解析解吻合良好，說明模型具有較為可靠的干濕處理能力。

2.4 斜坡靜水算例

為了說明模型的和諧特性，本文選取了Bonev等[20]采用的斜坡靜水算例。計算域為[-1.5 m,1.5 m]×[-1.5 m,1.5 m]的矩形區域，底坡b僅為x的函數，表達式如式(28)所示

(28)

4-a t=1 s時斜坡靜水水位4-b t=1 s時斜坡靜水x向通量圖4 斜坡靜水算例t=1 s時計算結果Fig.4 The free surface elevation and the flux magnitude at t=1 s of the water at rest on a sloping beach

靜水面設置在z=0 m處，g=9.8，最小水深H0=1×10-4m，四周邊界均設置為可滑移固壁邊界，多項式階數P=2。這里將未對底坡源項進行拆分的間斷有限元淺水模型稱為傳統方法，圖4-a分別繪制了模擬至1 s時本模型和傳統方法計算的水位，由圖可知，本模型計算的水位基本保持為靜水位而傳統方法計算的水位則在干濕交接附近有微弱波動。圖4-b分別繪制了模擬至1 s時本模型和傳統二維淺水方程計算的x向通量，此圖可以更直觀地看出傳統方法存在虛假流動的問題，而本模型通過對底坡源項進行拆分，保證了各項階數的平衡，從而很好地保持了斜坡上的靜水特性。

3 模型應用

5-a 網格劃分5-b 測點位置圖5 三亞紅塘灣算例網格劃分及測點布置示意圖Fig.5 Computing grid and locations of measuring stations of Hongtang Bay in Sanya

為了驗證加入源項后的二維水動力模型模擬實際潮流的合理性，選取三亞市紅塘灣附近海域，模擬了2016年4月26日～27日一個完整大潮期間的潮流運動。本算例計算域共劃分為2 448個任意四邊形網格，如圖5-a所示，網格劃分由外海至近岸逐漸加密。由于該海域受北部灣復雜潮流影響，為了充分反映模擬海域的潮流運動，本文首先選取北部灣區域建立大模型，運用Mike21的FM模塊模擬了大潮期間北部灣大范圍的潮流運動，本算例開邊界使用的潮位數據由大模型計算結果插值得到。本算例中科氏力系數根據網格點的實際緯度得到，糙率n根據當地泥沙粒徑分布及地形特點確定取值為0.017，計算階數P=1。

圖6 2016年4月26日～27日大潮期潮位模擬值與實測值對比Fig.6 Comparison of surface elevation between simulation results and observed data

從2016年4月26日～27日的水文實測資料與計算結果進行驗證，選取大潮T1、T2兩個測站的潮位資料，V1～V6共6個測站的流速流向測站資料，潮位及流速測站位置如圖5-b所示。潮位驗證結果如圖6所示，流速驗證結果如圖7所示，流向驗證結果如圖8所示。由圖可知，模型模擬的潮流運動與實際觀測的情況吻合良好，模型計算結果較為可靠。

圖7 2016年4月26日～27日大潮期流速模擬值與實測值對比Fig.7 Comparison of current speed between simulation results and observed data

圖8 2016年4月26日～27日大潮期流向模擬值與實測值對比Fig.8 Comparison of current direction between simulation results and observed data

4 結論與展望

基于節點無積分間斷有限元和任意四邊形網格建立二維淺水方程模型，在初始控制方程中加入了科氏力、風應力和底摩阻項，使得模型可以模擬風吹流及實際潮流運動。采用赤道孤立Rossby波算例驗證了科氏力項的合理性，通過求解Stommel問題的穩定解對科氏力、風應力和底摩阻共同作用下的模型進行了驗證。在源項驗證完畢后，利用平緩岸灘波浪爬坡算例和斜坡靜水算例分別說明了模型的干濕處理能力和靜水模擬的和諧特性。模型最后應用于三亞市紅塘灣海域，對實際大潮過程進行了模擬，模擬結果與全潮水文觀測的水位和流速、流向結果吻合良好。本文建立的模型已經具備了模擬河口海岸二維水動力的基本功能，今后將進一步在模型中加入粘性應力項及輻射應力項，完善波流耦合功能進而實現實際情況下河口海岸地區風、浪、潮相互作用的風暴潮過程模擬。