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

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

泥沙運(yùn)動(dòng)對(duì)于河口海岸演變、航道淤積以及河口治理等問題有著重要的影響，一直是河口海岸地形演變過程的研究重點(diǎn)。近年來，隨著計(jì)算機(jī)技術(shù)的發(fā)展，數(shù)值模型越來越成為研究泥沙運(yùn)動(dòng)的重要手段。目前，二維水沙模型在解決工程泥沙問題中的應(yīng)用已十分普及。過去幾十年來，有限差分法、有限元法和有限體積法等經(jīng)典數(shù)值方法均在二維水沙模型中得到大量應(yīng)用。有限差分法建立在經(jīng)典數(shù)學(xué)逼近理論的基礎(chǔ)上，簡(jiǎn)潔直觀，但差分法一般不易保證計(jì)算的守恒性[1-2]；有限單元法計(jì)算精度相對(duì)較高，但不適于計(jì)算間斷問題[3-4]；有限體積法在實(shí)際應(yīng)用中，可以保證不同網(wǎng)格下的積分守恒，能夠計(jì)算間斷問題，但很難在一般非規(guī)則網(wǎng)格上取得高精度[5-6]。與上述三種經(jīng)典方法相比，間斷有限元法(Discontinuous Galerkin Method)因在數(shù)值精度、間斷捕捉能力、網(wǎng)格自適應(yīng)特性、緊致性、高度可并行性以及適用于非結(jié)構(gòu)化網(wǎng)格等方面的優(yōu)勢(shì)，近年來逐漸開始在二維水沙模型的研究中得到應(yīng)用[7]。

Gourgue[8]較早采用間斷有限元方法建立了二維水沙數(shù)值模型，并應(yīng)用于Scheldt河口泥沙運(yùn)動(dòng)研究。石寶海和陳煥貞[9]基于迎風(fēng)間斷有限元方法建立了平面二維水沙輸運(yùn)模型，泥沙運(yùn)動(dòng)采用挾沙力模型描述，并給出了一個(gè)理想算例結(jié)果。趙張益和張慶河[10](2014)采用間斷有限元法建立了具有二階精度的非結(jié)構(gòu)化網(wǎng)格二維懸沙數(shù)值模型，模型能夠有效限制間斷處的數(shù)值振蕩，保證穩(wěn)定性，但是該模型未考慮推移質(zhì)輸沙和地形演變的影響。

上述研究表明，間斷有限元可以用于二維水沙運(yùn)動(dòng)的模擬，然而，這些模型中包含的正交項(xiàng)是基于全階積分的，這需要大量的計(jì)算資源[11]。為了減少這種計(jì)算量，Atkins和Shu[12]提出了無積分節(jié)點(diǎn)間斷有限元格式，該格式避免了離散方程時(shí)的積分過程，能夠有效減少計(jì)算時(shí)間以及存儲(chǔ)空間，并可以保持間斷有限元固有的緊致性和穩(wěn)健性，具有良好的應(yīng)用前景，但是該方法要求單元的雅克比系數(shù)為常數(shù)，不能適用于任意四邊形網(wǎng)格。李龍翔和張慶河[13]在此基礎(chǔ)上建立了適用于任意四邊形網(wǎng)格的無積分節(jié)點(diǎn)間斷有限元模型，用矩陣運(yùn)算進(jìn)行數(shù)值積分，有效地減小了四邊形網(wǎng)格上模型計(jì)算所需的存儲(chǔ)空間與計(jì)算量。李文俊等[14]在二維淺水方程中加入了科氏力、風(fēng)應(yīng)力、底摩阻項(xiàng)，基于無積分節(jié)點(diǎn)間斷有限元建立了充分反映各種實(shí)際條件影響的二維水動(dòng)力模型。

本文將在李文俊等[14]二維水動(dòng)力模型基礎(chǔ)上，進(jìn)一步基于無積分節(jié)點(diǎn)間斷有限元建立二維泥沙運(yùn)動(dòng)與地形演變模型，同時(shí)，本模型采用了李龍翔和張慶河[15]提出的采用權(quán)重重構(gòu)方法的節(jié)點(diǎn)型斜率限制器來消除間斷處產(chǎn)生的振蕩，為二維高效高精度無積分節(jié)點(diǎn)間斷有限元二維水沙模型的建立奠定基礎(chǔ)。

1 二維泥沙運(yùn)動(dòng)及地形演變數(shù)學(xué)模型

1.1 水動(dòng)力方程

模型的水動(dòng)力計(jì)算采用二維淺水方程作為控制方程，該方程是由三維N-S方程在采用靜壓假定，忽略粘性以及垂向加速度的情況下，沿水深積分得到的

(1)

式中：U= [h,hu,hv]T，F(xiàn)(U)=[E(U),G(U)]T，S(U)分別代表守恒變量、通量項(xiàng)以及源項(xiàng)，表達(dá)式寫為

(2)

(3)

式中：h為總水深；g為重力加速度；u和v分別為x方向和y方向的垂線平均流速；z為底部高程；η=h+z為水位；Sfx和Sfy分別為x方向和y方向底摩阻項(xiàng)，可表示為

Sfx=n2u(u2+v2)1/2h-4/3，Sfy=n2v(u2+v2)1/2h-4/3

(4)

式中：n為曼寧系數(shù)。

1.2 泥沙運(yùn)動(dòng)控制方程

泥沙運(yùn)動(dòng)暫時(shí)只考慮非粘性泥沙，按照其運(yùn)動(dòng)型式可分為推移質(zhì)和懸移質(zhì)，下面分別給出描述泥沙運(yùn)動(dòng)的推移質(zhì)運(yùn)動(dòng)公式和懸移質(zhì)運(yùn)動(dòng)控制方程。

1.2.1 推移質(zhì)運(yùn)動(dòng)控制方程

推移質(zhì)運(yùn)動(dòng)控制方程一般用經(jīng)驗(yàn)公式表示，本文采用van Rijn[16](2007)公式

qbx=0.015uh(d50/h)1.2Me1.5qby=0.015vh(d50/h)1.2Me1.5

(5)

式中：qbx和qby分別為x和y方向的推移質(zhì)泥沙濃度；d50為泥沙中值粒徑；Me表達(dá)式為

Me=max(0，|uc|-ucr,c)/[(s-1)gd50]0.5

(6)

式中：uc為水平方向的合速度；s為泥沙密度與液體密度的比值；ucr,c為泥沙運(yùn)動(dòng)的臨界流速，表達(dá)式為

(7)

式中：d90表示一個(gè)樣品的累計(jì)粒度分布數(shù)達(dá)到90%時(shí)所對(duì)應(yīng)的粒徑。

1.2.2 懸移質(zhì)運(yùn)動(dòng)控制方程

懸移質(zhì)泥沙運(yùn)動(dòng)控制方程采用沿水深平均的二維對(duì)流擴(kuò)散方程

(8)

式中：c為沿垂向平均后的泥沙濃度；εx和εy為泥沙的水平擴(kuò)散系數(shù)；Fs為源項(xiàng)。

源項(xiàng)可以采用挾沙力公式或參考濃度法描述，本文采用參考濃度法[17]進(jìn)行計(jì)算

Fs=P-D=(ca-c0)wf

(9)

式中：P為泥沙起懸通量；D為泥沙沉積通量；wf為泥沙沉速。ca和c0是定義在參考高度處a的濃度，ca與底床切應(yīng)力有關(guān)，表示懸沙平衡時(shí)參考高度處的泥沙濃度，c0與對(duì)流擴(kuò)散方程的垂向平均濃度有關(guān)，表示在水深平均濃度對(duì)應(yīng)于參考高速度處的濃度。當(dāng)ca>c0時(shí)，泥沙起懸，當(dāng)ca

(10)

式中：τs,max為最大底部應(yīng)力；τcr為泥沙顆粒臨界啟動(dòng)切應(yīng)力；βd為垂向濃度分布轉(zhuǎn)換系數(shù)，由下式計(jì)算

(11)

(12)

式中：κ為馮卡門系數(shù)，一般取0.4；u*c為底摩阻流速。

1.3 地形演變方程

推移質(zhì)和懸移質(zhì)的泥沙濃度會(huì)引起底床的沖刷和淤積變形，底床的地形演變可通過以下方程描述

(13)

式中：p為泥沙的孔隙率。

2 控制方程數(shù)值離散

水動(dòng)力模型中淺水方程的數(shù)值離散過程已在文獻(xiàn)[14]中給出，下面重點(diǎn)給出懸沙輸移對(duì)流擴(kuò)散方程與地形演變方程的離散過程。

2.1 對(duì)流擴(kuò)散方程與地形演變方程的離散

式(8)中含有2階偏導(dǎo)項(xiàng)，為將方程改為一階偏導(dǎo)方程，引入輔助變量[19]

Q=εhc

(14)

此時(shí)二維對(duì)流擴(kuò)散方程可改寫為

(15)

為方便離散，將(13)、(15)式寫為

(16)

式中：W為守恒通量，W=[hc,z]T，P(W)=[K(W),L(W)]和R(W)分別代表守恒變量通量項(xiàng)以及源項(xiàng)，表達(dá)式寫為

(17)

(18)

與水動(dòng)力模型的離散過程類似[14]，將求解域劃分為Ne個(gè)非重疊單元Ωe[6]，定義Ωe上至多為p階的多項(xiàng)式空間Vp(Ωe)，取Lagrange插值函數(shù)作為試驗(yàn)函數(shù)φVP(Ωe)，將試驗(yàn)函數(shù)乘以式 (14)和(16)，并在單元Ωe上積分，可得

(19)

(20)

對(duì)式(18)進(jìn)行分部積分，得到離散方程

(21)

式中：n為本單元邊界處的外法向向量，由于相鄰單元在交界面兩側(cè)的值不同，為保證兩側(cè)單元流入與流出的通量守恒，引入數(shù)值通量(hc)*，采用中心通量計(jì)算[20]

(hc)*=0.5(hc-+hc+)

(22)

(23)

2.2 時(shí)間遞進(jìn)

對(duì)控制方程離散后的強(qiáng)形式進(jìn)行進(jìn)一步離散，得到如下的半離散格式

(24)

式中：uh為未知變量的矢量。

式(28)的時(shí)間遞進(jìn)可以采取歐拉法、龍格庫塔法等求解。為了盡可能減小時(shí)間離散誤差，本文采用顯式四階龍格庫塔方法[22]用于半離散方程的時(shí)間遞進(jìn)。

(25)

式中：系數(shù)ai，bi，ci的取值參照文獻(xiàn)[22]。

3 模型驗(yàn)證與應(yīng)用

下面分別利用前人進(jìn)行的泥沙懸揚(yáng)、沙波演變和潰壩波沖淤水槽試驗(yàn)對(duì)建立的泥沙運(yùn)動(dòng)與地形演變模型進(jìn)行驗(yàn)證。最后，還將模型應(yīng)用于海南省紅塘灣附近海域的模擬，與現(xiàn)場(chǎng)實(shí)測(cè)懸沙資料進(jìn)行比較。

3.1 泥沙懸揚(yáng)算例

該算例為van Rijn[23]在圖1所示的水槽中所做的泥沙懸揚(yáng)實(shí)驗(yàn)。水槽中分為定床段和鋪沙段，鋪沙段長30 m、寬0.5 m、高0.7 m。水深為0.25 m，平均流速0.67 m/s，底床泥沙中值粒徑為0.23 mm。Lin 和 Falconer[24]曾用數(shù)學(xué)模型對(duì)此實(shí)驗(yàn)進(jìn)行了模擬，本文模擬時(shí)的相關(guān)系數(shù)均參考文獻(xiàn)[24]進(jìn)行取值。懸沙沉降速度取為0.022 m/s，水平擴(kuò)散系數(shù)取為0.002 m2/s。

圖2為模型模擬結(jié)果與實(shí)測(cè)垂向平均含沙量結(jié)果的對(duì)比，結(jié)果表明，本文模型的數(shù)值解與實(shí)測(cè)值基本接近，泥沙輸移趨勢(shì)相同，模型可以有效地模擬泥沙懸揚(yáng)以及輸移過程。

圖1 實(shí)驗(yàn)布置圖Fig.1 Layout of experiment圖2 垂向平均含沙量模擬值與實(shí)測(cè)值對(duì)比Fig.2 Comparison of simulated and measured vertical average sediment concentration

3.2 沙波演變算例

該算例在不考慮懸移質(zhì)泥沙條件下模擬沙波地形演變過程[25]，用以檢驗(yàn)?zāi)Ｐ湍M地形演變數(shù)值模型的合理性。初始沙波為正弦形，源項(xiàng)取為0，推移質(zhì)表達(dá)式取為

(26)

式中：a和b為常數(shù)；u=(ux，0)為沿水深平均的流速；D=(Dx，0)為流量；Δy=1.2 m為寬度；a=0.001 s2·m-1；b=3；Dx=1 m3·s-1，初始地形的函數(shù)可參考一維算例[26]。

圖3為y=0.75 m處剖面在t=500 s時(shí)模型模擬結(jié)果與文獻(xiàn)[25]給出的解析解的對(duì)比，結(jié)果表明，本文模型的模擬值與解析解吻合較好，模型可以精確地模擬地形演變過程。圖4為本文模型計(jì)算結(jié)果與文獻(xiàn)[25]基于有限體積方法給出考慮人工擴(kuò)散項(xiàng)的計(jì)算結(jié)果與解析解誤差的對(duì)比圖，可以看到，本文的模擬值誤差更小，精度更高。

圖3 沙波演變模型模擬值與解析解的對(duì)比Fig.3 Comparisons of simulated values and analytical solutions圖4 間斷有限元與文獻(xiàn)[25]考慮擴(kuò)散項(xiàng)有限體積法誤差對(duì)比Fig.4 Error comparisons of discontinuous Galerkin method and finite volume method with additional diffusion in reference[25]

3.3 潰壩動(dòng)床算例

圖5 潰壩算例的實(shí)驗(yàn)布置圖 (m)Fig.5 Layout of dam break experiment

該算例出自Leal的潰壩實(shí)驗(yàn)[27]，實(shí)驗(yàn)在一個(gè)長19.2 m、寬0.5 m、高0.7 m的水槽中進(jìn)行，圖5為初始時(shí)實(shí)驗(yàn)布置的示意圖，水槽分為左右兩端，左段在固定床面之上鋪沙高度為0.19 m，水深0.4 m，右側(cè)鋪沙高度為0.071 m，水深0.075 m。泥沙中值粒徑為1.22 mm，泥沙沉速為9.92 cm/s。

圖6和圖7分別顯示了t=1 s和t=4 s時(shí)模型模擬與實(shí)驗(yàn)測(cè)量水位值和地形值的比較情況。雖然水位、地形模擬值與實(shí)測(cè)值在局部有一定差異，二者變化趨勢(shì)總體上吻合較好，可以認(rèn)為模型能夠合理模擬復(fù)雜水流及其引起的泥沙沖淤導(dǎo)致的快速地形演變問題。

圖6 t=1 s時(shí)模擬值與實(shí)測(cè)值的對(duì)比Fig.6 Comparisons of simulated and measured values at t=1 s圖7 t=4 s時(shí)模擬值與實(shí)測(cè)值的對(duì)比Fig.7 Comparisons of simulated and measured values at t=4 s

3.4 模型在三亞紅塘灣泥沙運(yùn)動(dòng)模擬中的應(yīng)用

圖8 三亞紅塘灣泥沙測(cè)站布置圖Fig.8 Layout of sediment station in Hongtang Bay of Sanya

為了進(jìn)一步考察二維水沙模型模擬實(shí)際海域水沙運(yùn)動(dòng)的合理性，本節(jié)選取三亞市紅塘灣附近海域，模擬了2016年4月26日～27日一個(gè)完整大潮期間的泥沙運(yùn)動(dòng)。本算例計(jì)算域與文獻(xiàn)[9]相同，潮流模擬結(jié)果與實(shí)測(cè)資料的比較見文獻(xiàn)[14]，這里主要給出懸沙運(yùn)動(dòng)模擬結(jié)果。根據(jù)現(xiàn)場(chǎng)實(shí)測(cè)資料分析，懸沙中值粒徑取為0.3 mm，懸沙沉降速度取為0.035 m/s。

現(xiàn)場(chǎng)全潮水文觀測(cè)包含潮流與懸沙測(cè)量，這里選取S1～S6共六個(gè)測(cè)站的懸浮泥沙資料與模擬結(jié)果進(jìn)行比較。泥沙測(cè)站位置如圖8所示。模擬結(jié)果與實(shí)測(cè)結(jié)果比較見圖9。由圖可知，模型模擬的泥沙運(yùn)動(dòng)與實(shí)際觀測(cè)情況吻合良好?？梢哉J(rèn)為，本文建立的水沙模型能夠合理模擬比較復(fù)雜的現(xiàn)場(chǎng)實(shí)際情況。

圖9 2016年4月26日～27日大潮期懸沙濃度模擬值與實(shí)測(cè)值對(duì)比Fig.9 Comparison of simulated and measured suspended sediment concentration during high tide period (2016-04-26～27)

4 結(jié)論

基于無積分節(jié)點(diǎn)間斷有限元方法建立了泥沙輸移二維數(shù)學(xué)模型，懸沙運(yùn)動(dòng)通過對(duì)流擴(kuò)散方程進(jìn)行離散求解，方程的源項(xiàng)采用參考濃度法，推移質(zhì)輸沙采用經(jīng)驗(yàn)公式進(jìn)行計(jì)算，地形演變過程則通過推移質(zhì)和懸移質(zhì)輸沙控制的床面變形方程進(jìn)行求解。建立的泥沙輸移與地形演變數(shù)值模型獲得了已有水槽沖刷試驗(yàn)、沙波演變解析解和潰壩動(dòng)床實(shí)驗(yàn)的驗(yàn)證。模型還應(yīng)用于三亞紅塘灣大潮過程懸沙運(yùn)動(dòng)的模擬，模擬結(jié)果與全潮水文觀測(cè)多點(diǎn)含沙量過程吻合較好。模型驗(yàn)證和應(yīng)用結(jié)果表明，基于無積分節(jié)點(diǎn)間斷有限元建立的水沙模型可以合理模擬泥沙運(yùn)動(dòng)及其導(dǎo)致的地形演變問題。

為了更全面描述河口海岸泥沙二維運(yùn)動(dòng)問題，模型將進(jìn)一步考慮波流共同作用下二維水動(dòng)力及泥沙運(yùn)動(dòng)。