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

(天津大學(xué) 水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室，天津 300072)

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

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

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

1 二維水動(dòng)力數(shù)學(xué)模型

1.1 控制方程

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

(1)

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

源項(xiàng)S考慮了底坡、科氏力、底摩阻和表面風(fēng)應(yīng)力影響，最終表達(dá)式為

(2)

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

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

f=2ωsinφ

(3)

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

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

(4)

(5)

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

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

(6)

為了推導(dǎo)方便，下面的數(shù)值離散仍以矢量形式為主。

1.2 節(jié)點(diǎn)間斷有限元數(shù)值離散

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

(7)

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

(8)

對(duì)上式用分部積分和格林公式可得

(9)

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

(10)

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

(11)

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

(12)

(13)

其中，

(14)

(15)

(16)

對(duì)式(9)左端第二、第三項(xiàng)再運(yùn)用一次分部積分，得到方程組

(17)

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

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

1.3 無積分格式

傳統(tǒng)間斷有限元方法在求解式(17)時(shí)采用數(shù)值積分，這種方法要求被積函數(shù)在積分節(jié)點(diǎn)處的值已知，因此需要先計(jì)算出積分節(jié)點(diǎn)處的守恒變量值與通量值，從而造成計(jì)算量較大而影響計(jì)算效率。為了提高計(jì)算效率，Atkins和Shu[17]曾提出一種無積分格式，通過假定單元內(nèi)雅克比行列式為常數(shù)，將數(shù)值積分運(yùn)算轉(zhuǎn)化為標(biāo)準(zhǔn)單元內(nèi)矩陣計(jì)算。這一方法省去了數(shù)值積分過程，但其只適用于結(jié)構(gòu)化平行四邊形網(wǎng)格，對(duì)于任意四邊形單元并不適用。李龍翔和張慶河[10]提出的無積分格式無需假定單元內(nèi)雅克比行列式為常數(shù)，可適用于任意四邊形網(wǎng)格。該格式將積分節(jié)點(diǎn)和插值節(jié)點(diǎn)均選定為同一套邊界點(diǎn)的Legendre-Gauss-Lobatto(LGL)節(jié)點(diǎn)，在計(jì)算式(17)時(shí)將任意四邊形單元通過雅克比變換映射到標(biāo)準(zhǔn)四邊形單元，并將體積分中的單元?jiǎng)偠染仃嚤硎緸橘|(zhì)量矩陣與導(dǎo)數(shù)矩陣乘積形式，最終可得到標(biāo)準(zhǔn)單元坐標(biāo)系下包含NP個(gè)方程的方程組(18)

(18)

式中：NP為單元插值節(jié)點(diǎn)個(gè)數(shù)，1≤j≤NP，r、s為標(biāo)準(zhǔn)單元坐標(biāo)變量，J為雅克比變換矩陣行列式，由式(19)給出

(19)

無積分格式將數(shù)值積分轉(zhuǎn)化為上述矩陣的運(yùn)算，避免了求解積分節(jié)點(diǎn)處的函數(shù)值。同時(shí)，確定標(biāo)準(zhǔn)單元后，其導(dǎo)數(shù)矩陣、體積分質(zhì)量矩陣和面積分質(zhì)量矩陣都為定值，可以預(yù)先調(diào)用矩陣運(yùn)算庫計(jì)算，從而有效減少計(jì)算量并節(jié)省大量?jī)?chǔ)存空間。特別需要指出的是，當(dāng)采用高階精度格式時(shí)，無積分運(yùn)算求解數(shù)值積分產(chǎn)生的誤差與采用高階基函數(shù)近似產(chǎn)生的截?cái)嗾`差同階[18]。因此，無積分方法同樣對(duì)于格式計(jì)算誤差影響不大。

1.4 干濕處理及和諧特性

1.4.1 干濕處理

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

1.4.2 和諧性

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

(20)

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

(21)

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

2 模型驗(yàn)證

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

2.1 赤道孤立Rossby波

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

(22)

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

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

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

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

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

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

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

2.2 平底坡Stommel問題

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

(23)

(24)

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

(25)

2-a 模擬結(jié)果水位等值線2-b 模擬結(jié)果流速等值線圖2 穩(wěn)定狀態(tài)水深等值線圖及流速等值線圖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，各參數(shù)選取均與文獻(xiàn)[29]相同。

在科氏力、風(fēng)應(yīng)力和底摩阻的共同作用下，計(jì)算域內(nèi)初始會(huì)在中心形成一個(gè)凸起的波，之后波峰慢慢向左側(cè)偏移，當(dāng)偏移一定距離后達(dá)到穩(wěn)定狀態(tài)，波的形狀和位置均不再發(fā)生變化。結(jié)果達(dá)到穩(wěn)定時(shí)計(jì)算域內(nèi)水位等值線和流速等值線如圖2所示，與文獻(xiàn)[29]中的結(jié)果相近。

2.3 平緩岸灘波浪爬坡

為了說明本模型的干濕處理能力，考慮一個(gè)由海底滑坡引起的波浪在平緩岸灘上爬坡的算例。本算例長(zhǎng)度為[-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時(shí)爬坡和回落過程波面圖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。計(jì)算域劃分為10 000個(gè)網(wǎng)格，計(jì)算時(shí)間為280 s，多項(xiàng)式階數(shù)P=2。圖3-b給出了t=160 s、175 s和220 s三個(gè)時(shí)刻的波面模擬結(jié)果。由圖可知，本模型模擬的波浪爬坡(t=220 s)、回落(t=160 s、t=175 s)波面和Carrier等[30]給出的半解析解吻合良好，說明模型具有較為可靠的干濕處理能力。

2.4 斜坡靜水算例

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

(28)

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

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

3 模型應(yīng)用

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

為了驗(yàn)證加入源項(xiàng)后的二維水動(dòng)力模型模擬實(shí)際潮流的合理性，選取三亞市紅塘灣附近海域，模擬了2016年4月26日～27日一個(gè)完整大潮期間的潮流運(yùn)動(dòng)。本算例計(jì)算域共劃分為2 448個(gè)任意四邊形網(wǎng)格，如圖5-a所示，網(wǎng)格劃分由外海至近岸逐漸加密。由于該海域受北部灣復(fù)雜潮流影響，為了充分反映模擬海域的潮流運(yùn)動(dòng)，本文首先選取北部灣區(qū)域建立大模型，運(yùn)用Mike21的FM模塊模擬了大潮期間北部灣大范圍的潮流運(yùn)動(dòng)，本算例開邊界使用的潮位數(shù)據(jù)由大模型計(jì)算結(jié)果插值得到。本算例中科氏力系數(shù)根據(jù)網(wǎng)格點(diǎn)的實(shí)際緯度得到，糙率n根據(jù)當(dāng)?shù)啬嗌沉椒植技暗匦翁攸c(diǎn)確定取值為0.017，計(jì)算階數(shù)P=1。

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

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

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

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

4 結(jié)論與展望

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