非結(jié)構(gòu)網(wǎng)格數(shù)值模擬映射的河流能量公式推導(dǎo)及探析

李彬, 孫東坡

(1.中水北方勘測(cè)設(shè)計(jì)研究有限責(zé)任公司,天津 300222; 2.華北水利水電大學(xué) 水利學(xué)院,河南 鄭州 450046)

河流水動(dòng)力要素的千變?nèi)f化及河床演變趨勢(shì),本質(zhì)上都是由河流能量的分布及調(diào)整轉(zhuǎn)換決定的。從伯努利方程到河流熵理論,河流能量研究始終是一個(gè)被關(guān)注的話題。目前,國(guó)內(nèi)外學(xué)者在河流能量方面的研究主要集中在應(yīng)用最小能耗率原理或最小熵產(chǎn)生原理分析、預(yù)測(cè)自然條件下或是修建攔河壩、丁壩等一系列河勢(shì)控導(dǎo)工程后,河流某一河段的河床整體演變趨勢(shì),而沒(méi)有從能量分布及轉(zhuǎn)換角度研究河流內(nèi)部的調(diào)整規(guī)律[1-6]。周筑寶與徐國(guó)賓、楊志達(dá)(YANG C T)曾就最小能耗率原理或最小熵產(chǎn)生原理可否應(yīng)用于河流系統(tǒng)研究進(jìn)行了深入的討論[7-10],說(shuō)明將這些涉及河流能量分布及耗散的原理應(yīng)用到河流系統(tǒng)的理論分析還需進(jìn)一步完善,應(yīng)用條件還有待于進(jìn)一步研究。

由于河流的河床演變過(guò)程是一個(gè)河段內(nèi)部(微觀)調(diào)整、全河段(宏觀)再平衡的過(guò)程,河床演變必然伴隨著河流能量的變化;同樣,河流能量的改變也會(huì)有與之相匹配的河床變化。因此,可以通過(guò)分析河流全河段各區(qū)域的能量分布及變化情況來(lái)研究河流系統(tǒng)能量的變化趨勢(shì),進(jìn)而預(yù)測(cè)河流的演變趨勢(shì)。本文基于計(jì)算流體力學(xué)理論,從河流網(wǎng)格單元體系的能量分布特征出發(fā),分析了河流能量(機(jī)械能)的組成;針對(duì)河流網(wǎng)格單元體系中的非結(jié)構(gòu)三角網(wǎng)格單元,分析了三角網(wǎng)格單元內(nèi)水體的能量(機(jī)械能)形式,推導(dǎo)出了計(jì)算單元河流機(jī)械能各部分的表達(dá)式,引入泥沙起動(dòng)公式對(duì)河流機(jī)械能的變化性質(zhì)及河床調(diào)整的內(nèi)部轉(zhuǎn)換特征進(jìn)行了討論,并以概化后的順直河段為例,進(jìn)行了初步驗(yàn)證和應(yīng)用研究。

1 基于離散三角單元的河流機(jī)械能

在平面二維數(shù)值模擬中,流場(chǎng)被離散成三角網(wǎng)格單元體系統(tǒng)。根據(jù)歐拉的分析理念,流場(chǎng)的物理屬性便可由原來(lái)的微元體物理屬性轉(zhuǎn)變?yōu)槿蔷W(wǎng)格單元體物理屬性;水流的總機(jī)械能便可由微元體機(jī)械能集合轉(zhuǎn)而用三角網(wǎng)格單元體機(jī)械能集合表征。每一個(gè)離散三角網(wǎng)格單元體的水流機(jī)械能包括三角網(wǎng)格單元體的水流重力勢(shì)能和動(dòng)能,此部分能量亦可稱為動(dòng)水總流能量,而其中的水流重力勢(shì)能稱為動(dòng)水重力勢(shì)能。本文定義三角網(wǎng)格單元體底面(河床床面)高程與重力勢(shì)能基準(zhǔn)之間的流體重力勢(shì)能為靜水重力勢(shì)能,分別對(duì)三角網(wǎng)格單元體內(nèi)的動(dòng)水總流能量公式、河床床面所具有的靜水重力勢(shì)能公式以及總機(jī)械能公式進(jìn)行推導(dǎo),同時(shí)引入泥沙起動(dòng)公式對(duì)總機(jī)械能的變化性質(zhì)及內(nèi)部轉(zhuǎn)換特征進(jìn)行討論。

根據(jù)MIKE21軟件的控制方程組(沿水深平均的二維淺水流動(dòng)質(zhì)量和動(dòng)量守恒控制方程組)及其離散格式,模型計(jì)算出的流速V、水位z和水深h等水力要素保存在三角網(wǎng)格單元體的幾何重心點(diǎn)上[11]。圖1給出了過(guò)流斷面和水力要素與三角網(wǎng)格單元體的幾何關(guān)系。其中,圖1(b)是圖1(a)在水深hn/2處的橫截面(三角網(wǎng)格單元)剖分圖。圖1中:L為過(guò)流斷面;V為三角網(wǎng)格單元的平均流速；u為三角網(wǎng)格單元沿x軸方向的流速；v為三角網(wǎng)格單元沿y軸方向的流速。

圖1 過(guò)流斷面和水力要素與三角網(wǎng)格單元體之間的幾何關(guān)系示意

1.1 流場(chǎng)三角網(wǎng)格單元的水流能量推導(dǎo)

1.1.1 標(biāo)量法網(wǎng)格單元邊流速計(jì)算

因?yàn)閡2+v2≠0,所以過(guò)流斷面L(與流速V方向垂直)與xoy平面的交線L交和x軸的夾角γ分為4種計(jì)算情況:

(1)

因?yàn)?xj-xi)2+(yj-yi)2≠0,所以三角網(wǎng)格單元各邊與x軸的夾角βk分為4種計(jì)算情況:

(2)

式中:i=1、2、3;j=1、2、3;k=1、2、3且i≠j≠k,遵循右手定則。例如:A(x1,y1)B(x2,y2)邊為邊l3,l3與x軸的夾角為β3,l3與交線L交的夾角為α3,通過(guò)l3的水流流速為V3,通過(guò)l3的流量為Q3。其余兩邊依次循環(huán)類推。

根據(jù)三角形外角性質(zhì),三角網(wǎng)格單元各邊與交線L交的夾角αk=|γ-βk|。根據(jù)過(guò)流斷面交線L交與三角網(wǎng)格單元各邊的幾何關(guān)系,三角網(wǎng)格單元各邊的平均流速為:

Vk=V·|cosαk|, 0≤αk<π。

(3)

式中:Vk為通過(guò)三角網(wǎng)格單元各邊的平均流速,m/s;V為通過(guò)三角網(wǎng)格單元的平均流速,m/s。

1.1.2 矢量法網(wǎng)格單元邊流速計(jì)算

因?yàn)?xj-xi)2+(yj-yi)2≠0,所以三角網(wǎng)格單元的各邊垂面通過(guò)重心點(diǎn)(Xn,Yn)的法向量J在各邊上的坐標(biāo)(Xk,Yk)可由以下方程組計(jì)算得到:

(4)

式中:i=1、2、3;j=1、2、3;k=1、2、3且i≠j≠k,遵循右手定則。

根據(jù)流速V與三角網(wǎng)格單元各邊垂面法向量J的矢量關(guān)系及矢量乘法法則,三角網(wǎng)格單元各邊的平均流速Vk為:

(5)

式中Vk為正值時(shí)表示流入網(wǎng)格,為負(fù)值時(shí)表示流出網(wǎng)格。

1.1.3 網(wǎng)格單元通量推導(dǎo)

根據(jù)三角網(wǎng)格單元各邊的平均流速Vk,三角網(wǎng)格單元各邊通過(guò)的流量可用下式計(jì)算:

Qk=Vkhklk。

(6)

式中:Qk為通過(guò)三角網(wǎng)格單元各邊的流量,m3/s;hk為三角網(wǎng)格單元各邊的水深,其等于三角網(wǎng)格單元的平均水深hn,m;lk為三角網(wǎng)格單元各邊的邊長(zhǎng),m。

因?yàn)槿蔷W(wǎng)格單元流入的流量和流出的流量相等,只能是三角網(wǎng)格單元的一個(gè)邊流入、另兩個(gè)邊流出或是兩個(gè)邊流入、另一個(gè)邊流出。所以,三角網(wǎng)格單元的流量值等于三角網(wǎng)格單元3條邊流量值和的1/2或流量最大邊的值,即:

(7)

以質(zhì)量計(jì)的通量可通過(guò)類比法應(yīng)用公式(7)確定。

1.1.4 網(wǎng)格單元能量方程推導(dǎo)

根據(jù)伯努利能量方程[12],每個(gè)單元格的動(dòng)水總流能量方程為:

(8)

對(duì)于天然河流,Pn取自由面相對(duì)壓強(qiáng),Pn=0;取動(dòng)能修正系數(shù)αn等于1則:

(9)

式中:ρ為水的密度,ρ=1 000 kg/m3;g為重力加速度,g=9.81 m2/s;Qn為通過(guò)三角網(wǎng)格單元的流量,m3/s;Hn為三角網(wǎng)格單元的總水頭,m;zn為三角網(wǎng)格單元的平均水位,m;z0為計(jì)算重力勢(shì)能的基準(zhǔn)水位,m;Vn為三角網(wǎng)格單元的平均流速,m/s;n為網(wǎng)格單元編號(hào)。

1.2 靜水三角網(wǎng)格單元的重力勢(shì)能推導(dǎo)

結(jié)合前述的靜水重力勢(shì)能定義,由圖1可知靜水重力勢(shì)能的計(jì)算公式為:

(10)

其中，

式中:e1靜n為靜水重力勢(shì)能,kJ;sn為三角網(wǎng)格單元的面積,m2。

1.3 三角網(wǎng)格單元總機(jī)械能

假設(shè)河床泥沙為散粒體均勻沙,泥沙起動(dòng)流速可用沙莫夫公式確定[13]，如下:

(11)

式中:Ue為泥沙起動(dòng)流速,m3/s;ρs為泥沙的密度,一般取2 650 kg/m3;d為泥沙粒徑,m。

根據(jù)河流動(dòng)力學(xué)原理,當(dāng)泥沙起動(dòng)流速Ue等于三角網(wǎng)格單元的水流平均流速Vn時(shí),可以得到河床沖刷穩(wěn)定條件下的三角網(wǎng)格單元極限平均水深hn，

(12)

將式(12)代入式(10)可得靜水重力勢(shì)能e1靜n的計(jì)算公式為:

(13)

由公式(9)與公式(13)可知:Vn增加,E動(dòng)n增加,e1靜n減小;反之,E動(dòng)n減小,e1靜n增加。水流是泥沙輸移的原動(dòng)力,河道水流運(yùn)動(dòng)的天然不穩(wěn)定性體現(xiàn)在流速的波動(dòng)上,這種波動(dòng)必然導(dǎo)致動(dòng)水總流能量與靜水重力勢(shì)能的相互轉(zhuǎn)換;由于河床泥沙運(yùn)動(dòng)相較于水流運(yùn)動(dòng)具有一定的滯后性(兩者波動(dòng)有相位差),這種滯后性又會(huì)進(jìn)一步影響能量的轉(zhuǎn)換特征。根據(jù)河流最小能耗原理[3],總機(jī)械能En總是趨向最小化,其計(jì)算公式如下:

En=E動(dòng)n+e1靜n=e1動(dòng)n+e2動(dòng)n+e1靜n→min。

(14)

式中:En為總機(jī)械能;符號(hào)“→”表示趨向性;“min”表示最小值。

應(yīng)用公式(9)和公式(10)或公式(9)和公式(13)對(duì)模擬河段每個(gè)三角網(wǎng)格單元不同形式的能量分別進(jìn)行計(jì)算并求和,便可得到模擬河段處于沖刷平衡狀態(tài)時(shí)不同形式能量的數(shù)值,還可得到其分布狀況。通過(guò)流場(chǎng)能量的分布特征及其調(diào)整轉(zhuǎn)換情況,分析河流工程區(qū)域三角網(wǎng)格單元的能量變化情況,就可以對(duì)河流水動(dòng)力的變化及演變趨勢(shì)進(jìn)行預(yù)測(cè)。

2 用于能量分析的流場(chǎng)質(zhì)量守恒驗(yàn)證

2.1 模擬條件

構(gòu)造長(zhǎng)為1 000 m，寬為100 m，渠底高程為-1.500 m的平坡順直河道作為模擬河段,對(duì)推導(dǎo)的公式進(jìn)行流場(chǎng)質(zhì)量守恒驗(yàn)證。模型進(jìn)口流量為100 m3/s,出口水位為0.500 m,模型網(wǎng)格單元節(jié)點(diǎn)個(gè)數(shù)為1 390,三角形網(wǎng)格單元總數(shù)為2 558。模型四邊固定邊界的網(wǎng)格單元邊長(zhǎng)設(shè)置為10.0 m,區(qū)域內(nèi)網(wǎng)格單元面積設(shè)置為不超過(guò)60.0 m2,并剖分生成網(wǎng)格單元。模型局部網(wǎng)格剖分及出口斷面網(wǎng)格單元編號(hào)如圖2所示。

圖2 模型局部網(wǎng)格剖分及出口斷面網(wǎng)格單元編號(hào)示意

2.2 結(jié)果及分析

應(yīng)用MIKE21軟件進(jìn)行平面二維水動(dòng)力數(shù)值模擬,得到模型出口斷面網(wǎng)格單元的水力要素(表1),據(jù)此采用三角網(wǎng)格單元流量公式計(jì)算得到模型出口斷面的網(wǎng)格單元流量(表2和表3)。應(yīng)用MIKE21軟件自帶的斷面流量監(jiān)測(cè)功能監(jiān)測(cè)到的模型出口斷面流量為99.909 7 m3/s,比給定流量偏小0.090 3 m3/s,誤差非常小,僅為0.903‰,可以認(rèn)為MIKE21軟件計(jì)算出的流量守恒。

由表2和表3可知:①采用兩種三角網(wǎng)格單元流量計(jì)算方法得到的流量模擬計(jì)算值在小數(shù)點(diǎn)后保留4位有效數(shù)字時(shí)沒(méi)有差別,兩種方法的計(jì)算精度相當(dāng)。②應(yīng)用本文流量公式得到的模型出口斷面三角網(wǎng)格單元的流量和為100.678 4 m3/s,比給定流量稍大。分析其原因是:在模型出口邊界三角網(wǎng)格單元系統(tǒng)中存在鈍角三角形網(wǎng)格單元,導(dǎo)致相鄰網(wǎng)格單元的流量有重復(fù)計(jì)量。從圖2可以看出:33#網(wǎng)格單元為鈍角三角形網(wǎng)格單元,造成它與相鄰的1 552#網(wǎng)格單元之間有通量傳遞,部分流量計(jì)重。若將Q網(wǎng)33替換為33#網(wǎng)格單元出口斷面的單邊流量Q3,扣除其與1 552#網(wǎng)格單元之間的流量重計(jì)量,則模型出口斷面流量為99.790 7 m3/s,僅比給定流量偏小0.209 3 m3/s,誤差為2.093‰。因此,在對(duì)模型下開邊界附近的網(wǎng)格進(jìn)行剖分時(shí),應(yīng)注意避免出現(xiàn)鈍角三角形網(wǎng)格單元,這樣采用文中所給的兩種計(jì)算三角網(wǎng)格單元流量的方法都可以獲得滿意的計(jì)算精度。

通過(guò)以上分析基本可以認(rèn)為:應(yīng)用推導(dǎo)出的兩種三角網(wǎng)格單元流量計(jì)算方法得到的河流流量均守恒;模擬結(jié)果真實(shí)可信,兩種方法均可用于三角網(wǎng)格單元的通量及能量計(jì)算。

表1 模型出口斷面網(wǎng)格單元的水力要素統(tǒng)計(jì)

表2 由公式計(jì)算出的模型出口斷面網(wǎng)格單元流量(流速標(biāo)量法)

表3 由公式計(jì)算出的模型出口斷面網(wǎng)格單元流量(流速矢量法)

3 模擬河段的能量分布及分析

對(duì)上節(jié)概化順直河道的水動(dòng)力模擬計(jì)算結(jié)果應(yīng)用公式(9)和公式(10)計(jì)算模擬河段的能量,其中z0取模擬河段河底高程最低值z(mì)河底min(-1.500 m),以保證所有計(jì)算三角網(wǎng)格單元的勢(shì)能均不小于零。計(jì)算模擬河段每個(gè)三角網(wǎng)格單元的動(dòng)水重力勢(shì)能e1動(dòng)n(包含壓能)、動(dòng)水動(dòng)能e2動(dòng)n、動(dòng)水機(jī)械能E動(dòng)n、靜水重力勢(shì)能e1靜n和總機(jī)械能En,累加求和后得到模擬河段不同形式能量的值,結(jié)果見表4。圖3給出了對(duì)基于流速標(biāo)量法的三角網(wǎng)格單元流量結(jié)果應(yīng)用公式(9)和公式(10)計(jì)算得到的模擬河段不同形式能量的分布,圖4給出了模擬河段不同形式能量的面密度(單位面積上的能量)分布。

表4 概化模擬河段的單位時(shí)間能量 kJ

圖3 概化模擬河段的單位時(shí)間能量分布

圖4 概化模擬河段的單位時(shí)間能量面密度分布

因?yàn)槟M河段是基于河床平坡的概化順直河道,且選取的勢(shì)能計(jì)算基準(zhǔn)是z河底min,所以e1靜n等于零,因此En與E動(dòng)n的數(shù)值相等、分布相同。分析表4可知:①基于流速標(biāo)量法、流速矢量法的流量結(jié)果,應(yīng)用公式(9)和公式(10)計(jì)算得到的模擬河段不同形式能量的值幾乎無(wú)差別,計(jì)算精度相當(dāng)。②選取z河底min作為勢(shì)能的計(jì)算基準(zhǔn),在本次模擬條件(常水位)下,模擬河段總體上e1動(dòng)占E動(dòng)的比例為994.109‰,e2動(dòng)占E動(dòng)的比例為5.891‰。這表明:當(dāng)選取z河底min作為勢(shì)能的計(jì)算基準(zhǔn)時(shí),在常水位條件下模擬河段所具有的能量主要是動(dòng)水重力勢(shì)能e1動(dòng)。這也從理論上給出了洪水期的河流特別是地上懸河的潰堤水流為什么會(huì)具有驚人破壞力的一個(gè)佐證。

圖3反映了模擬河段不同形式河流能量的平面分布特征。圖3表明:河流能量在空間上不是均勻分布的,沿縱橫向(網(wǎng)格單元體系間)都存在能量梯度,亦即都存在能量的傳遞與轉(zhuǎn)換。

圖4給出的不同形式河流能量的面密度分布特征,顯示了比圖3更為連續(xù)的分布變化特點(diǎn),同時(shí)顯示出能量縱向分布有波動(dòng)變化特征,也間接體現(xiàn)了沿程的能耗特點(diǎn);動(dòng)能和總機(jī)械能的縱橫向變化也表征了各網(wǎng)格單元體系間的能量梯度與傳遞特性。圖3和圖4表明:河流能量的損失一方面與邊界條件有關(guān);另一方面也遵循最小能耗原理,即總是使能量損失趨向最小化,它決定了能量耗散形式,同時(shí)也制約著流場(chǎng)中不同形式能量間的分配關(guān)系。

基于河流網(wǎng)格單元體系的能量分布特征,必然帶有網(wǎng)格的特征并會(huì)受其影響:模擬河段網(wǎng)格剖分不均勻,并非均一等面積網(wǎng)格;網(wǎng)格單元能量特征值受網(wǎng)格精度影響。這是需要在前期的網(wǎng)格剖分以及在模擬計(jì)算結(jié)果分析中必須關(guān)注的問(wèn)題。本文旨在提出一種基于計(jì)算流體力學(xué)的河流能量研究思路和方法,對(duì)于網(wǎng)格布置與精度(單元尺寸及其均勻化程度)對(duì)模擬結(jié)果的影響,以及通過(guò)分析三角網(wǎng)格單元能量的變化預(yù)測(cè)河流演變趨勢(shì),將在后續(xù)的研究中進(jìn)行。

4 結(jié)語(yǔ)

為了研究基于模擬流場(chǎng)的河流能量分布及轉(zhuǎn)換情況,通過(guò)分析河流能量(機(jī)械能)的組成,利用非結(jié)構(gòu)三角網(wǎng)格單元推導(dǎo)出了數(shù)值模擬區(qū)河流能量(機(jī)械能)的表達(dá)式,采用水動(dòng)力數(shù)值模擬軟件構(gòu)造了概化順直河道模型,并對(duì)推導(dǎo)的公式進(jìn)行了不同計(jì)算方法的守恒性驗(yàn)證和模擬計(jì)算,得到以下結(jié)論。

1)從能量角度分析,河流能量(機(jī)械能)包括由流動(dòng)水體產(chǎn)生的動(dòng)水總流能量和河床床面高程具備的靜水重力勢(shì)能。其中,動(dòng)水總流能量包括動(dòng)水重力勢(shì)能(包含壓能)和動(dòng)水動(dòng)能。河道水流的波動(dòng)導(dǎo)致動(dòng)水總流能量與靜水重力勢(shì)能頻繁轉(zhuǎn)換,河床泥沙運(yùn)動(dòng)相較于水流運(yùn)動(dòng)的滯后性又進(jìn)一步影響能量的轉(zhuǎn)換特征。根據(jù)河流最小能耗原理,河流能量(機(jī)械能)的消耗總是趨向最小化,它決定了能量耗散形式和不同形式能量間的分配關(guān)系。

2)河流能量在空間上不是均勻分布的,沿縱橫向(網(wǎng)格單元體系間)都存在能量梯度,亦即都存在能量的傳遞與轉(zhuǎn)換。河流能量面密度的分布特征顯示出能量縱向分布有波動(dòng)變化特征,也間接體現(xiàn)了沿程的能耗特點(diǎn);動(dòng)能和總機(jī)械能的縱橫向變化也表明各網(wǎng)格單元體系間存在能量梯度,河流能量存在波動(dòng)性分布不均與空間的相互傳遞。

3)基于非結(jié)構(gòu)三角網(wǎng)格單元推導(dǎo)出的兩種三角網(wǎng)格單元流量計(jì)算方法均可用于三角網(wǎng)格單元的流量計(jì)算,且計(jì)算精度相當(dāng)。應(yīng)用流速矢量法計(jì)算三角網(wǎng)格單元流量可使三角網(wǎng)格單元各邊的過(guò)流情況(大小及流向)一目了然。

4)基于非結(jié)構(gòu)三角網(wǎng)格單元推導(dǎo)出的河流能量(機(jī)械能)公式可用于河流能量計(jì)算。在基于非結(jié)構(gòu)三角網(wǎng)格單元的河流平面二維數(shù)值模擬中,可以通過(guò)該公式得到河段能量(機(jī)械能)的數(shù)值及其分布態(tài)勢(shì)。

5)本文研究對(duì)象為概化的順直河道模型,對(duì)于河道形態(tài)不規(guī)則的實(shí)際河流,三角網(wǎng)格單元尺寸及其均勻化程度對(duì)河流全河段能量計(jì)算結(jié)果的影響，以及沖刷平衡的極限狀態(tài)的討論有待于進(jìn)一步研究。