第2章　淮河中游河道一維、二維耦合水動(dòng)力數(shù)學(xué)模型

2.1　洪水演進(jìn)數(shù)學(xué)模型

隨著計(jì)算機(jī)和計(jì)算技術(shù)的不斷發(fā)展，水動(dòng)力數(shù)學(xué)模型發(fā)展迅速，在工程上得到了廣泛應(yīng)用，能較好地模擬復(fù)雜邊界條件下的水流運(yùn)動(dòng)，并與實(shí)體模型耦合運(yùn)用相互驗(yàn)證，相互提供邊界，可以較好地解決工程實(shí)際問(wèn)題。

淮河中游河勢(shì)復(fù)雜，有單一河道，有分汊河道，還有行洪區(qū)和湖泊。低水期時(shí)，水流主要在河道中運(yùn)動(dòng)。洪水期，隨著水位的上升，上游來(lái)流超過(guò)河道泄流能力時(shí)，通過(guò)分洪、滯洪等措施，使水流通過(guò)閘門(mén)、口門(mén)、漫堤等方式進(jìn)入到行洪區(qū)內(nèi)。為了正確地模擬上述水流的運(yùn)動(dòng)情況，需要根據(jù)不同的水流特征采用不同的模擬方法，以達(dá)到提高精度和效率的目的:①蓄洪洼地由于流速較小，水面比降小，主要考慮其水量調(diào)蓄的作用，采用水庫(kù)型水量平衡法;②干支流河道水流運(yùn)動(dòng)，需要計(jì)算水流的過(guò)流能力及水位，采用一維非恒定河網(wǎng)、堰閘聯(lián)合模型求解;③行洪區(qū)及湖泊由于流速大小及流向隨空間和時(shí)間變化，采用平面二維淺水方程求解。

2.1.1　一維水流數(shù)學(xué)模型

2.1.1.1　控制方程組及離散

河網(wǎng)一維水動(dòng)力模型的控制方程為Saint Venant方程組。

連續(xù)方程:

動(dòng)量方程:

式中:t為時(shí)間坐標(biāo);x為河道沿程坐標(biāo);Q為流量;Z為水位;A為過(guò)水?dāng)嗝娴拿娣e;B為水面寬度;K為流量模數(shù);g為重力加速度;q為旁側(cè)入流流量;n為河道糙率系數(shù);R為水力半徑。

利用Abbott六點(diǎn)隱式差分格式離散上述控制方程組，該離散格式在每個(gè)網(wǎng)格點(diǎn)并不同時(shí)計(jì)算水位和流量，而是按順序交替計(jì)算水位和流量，分別稱為h點(diǎn)和Q點(diǎn)，如圖2.1-1所示。該格式無(wú)條件穩(wěn)定，可以在相當(dāng)大的Courant數(shù)下保持穩(wěn)定，可以取較長(zhǎng)的時(shí)間步長(zhǎng)以節(jié)約計(jì)算時(shí)間[1-7]。

采用如圖2.1-2所示的Abbott六點(diǎn)隱式差分格式，連續(xù)性方程中的各項(xiàng)可以寫(xiě)為

于是連續(xù)性方程可以寫(xiě)為

同樣動(dòng)量方程中各項(xiàng)可以寫(xiě)為

于是動(dòng)量方程在流量點(diǎn)上的差分格式為

式(2.1-3)整理后可記為

式(2.1-4)整理后可記為

以上式中:αj、βj、γj、δj分別為離散方程的系數(shù)。

2.1.1.2　定解條件

初始條件:給定初始時(shí)刻t=0時(shí)，所有計(jì)算節(jié)點(diǎn)Q和h值。非恒定流數(shù)值計(jì)算表明，初始條件對(duì)于計(jì)算的初期階段會(huì)顯示影響，但這種影響將隨著計(jì)算時(shí)間的延伸逐步消失。

邊界條件分三類:

(1)水位邊界條件:邊界處水位h=h(t)。

(2)流量邊界條件:邊界處流量Q=Q(t)。

(3)水位—流量關(guān)系邊界條件:邊界處Q=f(h)給定。

2.1.1.3　求解方法

如前所述，河道內(nèi)任一點(diǎn)的水力參數(shù)Z(水位h或流量Q)與相鄰的網(wǎng)格點(diǎn)的水力參數(shù)的關(guān)系可以表示為統(tǒng)一的線性方程:

上式中的系數(shù)可分別由式(2.1-5)或式(2.1-6)計(jì)算。

假設(shè)一河道有n個(gè)網(wǎng)格點(diǎn)，因?yàn)楹拥赖氖啄┚W(wǎng)格點(diǎn)總是水位點(diǎn)，所以n是奇數(shù)。對(duì)于河網(wǎng)的所有網(wǎng)格點(diǎn)寫(xiě)出式(2.1-7)，可以得到n個(gè)線性方程:

其中第一個(gè)方程的Hus和最后一個(gè)方程中Hds分別是上、下游汊點(diǎn)的水位。某一河道第一個(gè)網(wǎng)格點(diǎn)的水位等于與之相連河段上游汊點(diǎn)的水位:α1=-1，β1=1，γ1=0，δ1=0。同樣，α1=0，β1=1，γ1=-1，δ1=0。

對(duì)于單一河道，只要給出上下游水位邊界，即Hus和Hds為已知，就可用消元求解方程組式(2.1-8)。對(duì)于河網(wǎng)問(wèn)題，由方程組式(2.1-8)，通過(guò)消元法可以將河道內(nèi)任意點(diǎn)的水力參數(shù)(水位或流量)表示為上下游汊點(diǎn)水位的函數(shù):

只要先求河網(wǎng)各汊點(diǎn)的水位，就可用式(2.1-9)求解河段任意網(wǎng)格點(diǎn)的水力參數(shù)。

圖2.1-3為河網(wǎng)汊點(diǎn)方程示意圖，圍繞汊點(diǎn)的控制體連續(xù)方程為

將式(2.1-10)右邊第二式的三項(xiàng)分別以式(2.1-9)替代，可以得到

式中:H為該汊點(diǎn)的水位;HA，us、HB，us分別為支流A、B上游端汊點(diǎn)水位;HC，ds為支流C下游端汊點(diǎn)水位。

在式(2.1-11)中，將某個(gè)汊點(diǎn)水位表示為與之直接相連的河道汊點(diǎn)水位的線性函數(shù)。同樣，對(duì)于河網(wǎng)所有汊點(diǎn)(假設(shè)為N個(gè))，可以得到N個(gè)類似的方程(汊點(diǎn)方程組)。在邊界水位或流量為已知的情況下，可以利用高斯消元法直接求解汊點(diǎn)方程組，得到各個(gè)汊點(diǎn)的水位，進(jìn)而回代式(2.1-9)求解河道任意網(wǎng)格點(diǎn)的水位和流量。

大型稀疏矩陣求解計(jì)算時(shí)間主要取決于矩陣主對(duì)角線非零元素的寬度，可通過(guò)對(duì)河網(wǎng)節(jié)點(diǎn)進(jìn)行優(yōu)化編碼的方法來(lái)降低汊點(diǎn)方程組系數(shù)矩陣的帶寬，使之稱為主對(duì)角元素占優(yōu)的矩陣，從而方便了方程組的求解，并大大減少了計(jì)算時(shí)間。

若在河道邊界節(jié)點(diǎn)上給出水位的時(shí)間變化過(guò)程，即h=h(t)，此時(shí)，假設(shè)邊界所在河道編號(hào)為j，則邊界上的汊點(diǎn)方程為

若在河道邊界節(jié)點(diǎn)上給出流量的時(shí)間變化過(guò)程，即Q=Q(t)。

對(duì)如圖2.1-4所示的控制體，應(yīng)用連續(xù)方程可以得到:

將以式(2.1-9)代入式(2.1-13)，可以得到

若在河道邊界節(jié)點(diǎn)上給出的是水位流量關(guān)系Q=Q(h)，其處理方法同流量邊界，得到與式(2.1-14)類似的方程，只是方程中的由水位流量關(guān)系計(jì)算得到。

平原河網(wǎng)大多有堰、閘等水工建筑物，此時(shí)Saint Venant方程已經(jīng)不再適用，必須根據(jù)堰、閘的水力學(xué)特性作特殊處理。在模型中堰、閘通常作為流量點(diǎn)處理，根據(jù)相鄰水位點(diǎn)水位關(guān)系采用寬頂堰或孔口出流計(jì)算過(guò)閘流量，得到式(2.1-6)類似的方程[1-4]。

2.1.2　平面二維水流數(shù)學(xué)模型

2.1.2.1　基本方程

大部分河流湖泊都是淺水水流，滿足以下假定:①具有自由表面;②以重力為主要的驅(qū)動(dòng)力，以水流和固體邊界之間及水流內(nèi)部摩阻力為主要的耗散力;③水平流速沿水深近似成均勻分布;④垂向流速和垂向加速度可忽略，水壓接近靜壓分布。基于以上假定，對(duì)三維水流的運(yùn)動(dòng)方程沿水深進(jìn)行積分，可得二維淺水水流控制方程。

連續(xù)性方程:

動(dòng)量方程:

式中:ξ為水位;h為水深，h=ξ-Zb，Zb為河床底高程;u、v分別為x，y方向的垂向平均流速;f為科氏力系數(shù)，f=2ωsinφ，ω為地球自轉(zhuǎn)的角速度，φ為計(jì)算水域的地理緯度;τbx、τby分別為x方向和y方向的底部摩阻，=n為曼寧系數(shù);τsx、τsy分別為風(fēng)對(duì)自由表面x方向和y方向的剪切力，cd為風(fēng)阻力系數(shù)，ρa為空氣密度，uw，vw為x方向和y方向的風(fēng)速;ρ為水密度;Ex、Ey分別為x方向和y方向的水流紊動(dòng)黏性系數(shù);q為源或匯單位面積的流量，源取正，匯取負(fù);u*、v*分別為源或匯輸入輸出時(shí)在x和y方向的流速。

上述二維淺水方程可以寫(xiě)成如下統(tǒng)一的形式:

式中:

q=(h，hu，hv)

2.1.2.2　定解條件

初始條件:給定初始時(shí)刻t=0時(shí)，計(jì)算域內(nèi)所有計(jì)算變量(u，v，ζ)的初始值。

u(x，y，t)t=0=u0(x，y)

v(x，y，t)t=0=v0(x，y)

ξ(x，y，t)t=0=ξ0(x，y)

邊界條件:包括固壁邊界和開(kāi)邊界。

(1)固壁邊界:近壁區(qū)因?yàn)榉肿羽ば暂^大，雷諾數(shù)較低，為避免在壁面處加密網(wǎng)格，通常采用壁面函數(shù)法和滑移邊界來(lái)處理近壁處的流速。壁面函數(shù)法通過(guò)一組半經(jīng)驗(yàn)公式直接將壁面上的物理量與湍流核心區(qū)內(nèi)待求的未知量直接聯(lián)系起來(lái);滑移邊界條件，即=0。本文計(jì)算采用滑移邊界條件處理固壁邊界。

(2)開(kāi)邊界:可以給定水位、流量或流速的過(guò)程。

給定水位邊界:ξ=ξ(t);

給定流量邊界:Q=Q(t);

給定法向流速過(guò)程:un=un(t)。

2.1.2.3　離散方法

Mike21 FM采用非結(jié)構(gòu)有限體積法離散控制方程[1-4]。有限體積法中使用的非結(jié)構(gòu)網(wǎng)格通常由三角形或四邊形網(wǎng)格構(gòu)成，為了準(zhǔn)確的逼近水下地形，這里僅采用三角形網(wǎng)格。

采用網(wǎng)格中心式布置計(jì)算節(jié)點(diǎn)，見(jiàn)圖2.1-5。將式(2.1-18)在控制體上進(jìn)行積分得:

式中:Ω為控制體平面域;dA為面積分微元。采用高斯格林公式將面積分化為沿其周界的線積分，得

式中:?Ω為控制體周界(逆時(shí)針?lè)较?;d S為線積分微元;nx、ny為周界上外法向單位向量;分別為守恒物理量、源匯項(xiàng)在控制體內(nèi)平均，當(dāng)精度小于或等于二階時(shí)，它們分別等于變量和源、匯項(xiàng)在三角形形心處的值。

記跨控制體邊界的法向通量Fn(q)=F(q)nx+G(q)ny，考慮到本次計(jì)算采用三角網(wǎng)格，故式(2.1-19)又可寫(xiě)為

式中:Lj為第j邊的邊長(zhǎng);Fn(q)j為第j邊平均法向數(shù)值通量。

由于Fn(q)沿控制體邊界一般為非線性分布，用數(shù)值求積公式計(jì)算可取得較高的精度。因此式(2.1-20)可寫(xiě)為

式中:為數(shù)值求積公式計(jì)算的第j邊的平均法向數(shù)值通量Fn(qj，k);ωk為積分權(quán);n為積分權(quán)總數(shù);qj，k為守恒物理量在j邊，第k積分點(diǎn)的值。

對(duì)于空間離散后得到的常微分方程，本文采用以下兩種常微分方程解算器進(jìn)行時(shí)間離散。記式(2.1-21)等號(hào)右端項(xiàng)為L(zhǎng)(q)，則式(2.1-21)可寫(xiě)為

(1)一階顯格式:

式中:為控制體內(nèi)守恒物理量的平均值，q由重構(gòu)的控制體內(nèi)守恒物理量，為(x，y)的函數(shù)。

(2)二階Runge-kutta法:

至此，所有問(wèn)題歸結(jié)為如何求解跨界面處的數(shù)值通量。目前常用的途徑是基于間斷的思想，認(rèn)為在控制體界面兩側(cè)的物理量存在間斷，這樣在每個(gè)控制界面處可以形成一個(gè)黎曼問(wèn)題。可將F(qn)取為圖2.1-6所示的一維黎曼問(wèn)題的解。這樣一維黎曼問(wèn)題可以表示為

初始條件為

其中qL和qR分別為間斷左右的物理量值。

Mike21 FM采用Roe格式求解以上黎曼問(wèn)題，為避免數(shù)值振蕩，使用了二階TVD限制器。

2.1.3　一維、二維耦合水流數(shù)學(xué)模型

一維模型計(jì)算簡(jiǎn)便，適合模擬河道內(nèi)水流運(yùn)動(dòng)情況;二維模型需要占用更多的計(jì)算資源，適合模擬行洪區(qū)和湖泊內(nèi)水流運(yùn)動(dòng)情況。通過(guò)建立一維、二維耦合水流數(shù)學(xué)模型，將上述兩種模型連成一體，聯(lián)合求解，既可以發(fā)揮出一維模型快速方便的特點(diǎn)，同時(shí)又能獲得局部范圍的細(xì)部信息。Mikeflood包括以下三種類型的連接方式:

(1)標(biāo)準(zhǔn)連接。標(biāo)準(zhǔn)連接主要用于將一個(gè)或多個(gè)二維網(wǎng)格單元連接到一個(gè)一維模型的末端。這種連接方式適用于河口復(fù)雜水流的計(jì)算中，具體應(yīng)用如圖2.1-7所示。

(2)側(cè)向連接。側(cè)向連接主要用于一串二維網(wǎng)格側(cè)向連接到一維河道中(一部分河道或整條河道)。這種連接方式適用于模擬水流從河道進(jìn)入行洪區(qū)的過(guò)程，具體應(yīng)用如圖2.1-8所示。

(3)結(jié)構(gòu)物連接。結(jié)構(gòu)物連接主要應(yīng)用于將兩個(gè)分開(kāi)的二維網(wǎng)格通過(guò)一維結(jié)構(gòu)物連接成一個(gè)整體。這種連接方式適用于在二維網(wǎng)格中模擬結(jié)構(gòu)物對(duì)水流的影響，具體應(yīng)用如圖2.1-9所示。

對(duì)于上述連接方式，由于連接點(diǎn)不屬于邊界點(diǎn)，因此不論對(duì)于一維問(wèn)題、還是二維問(wèn)題都缺少邊界條件，無(wú)法獨(dú)立求解。需要通過(guò)在連接點(diǎn)處補(bǔ)充物理量之間的關(guān)系(水位、流量相等)，從而實(shí)現(xiàn)一維、二維模型的耦合。

按照耦合求解方式的不同，將其分為兩大類。標(biāo)準(zhǔn)連接和側(cè)向連接采用交替計(jì)算的方法，該方法是在一個(gè)時(shí)步內(nèi)把一個(gè)模型計(jì)算出來(lái)的物理量(水位或流量)作為另一個(gè)模型的邊界條件，兩模型交替計(jì)算，互賦邊界的方法完成整個(gè)計(jì)算過(guò)程，這種方式在時(shí)間步長(zhǎng)不大時(shí)，省去了迭代的過(guò)程，節(jié)省了計(jì)算的時(shí)間;結(jié)構(gòu)物連接是采用直接求解的方法，該方法將連接點(diǎn)處的連接條件作為定解條件引入方程組中，通過(guò)補(bǔ)充迭代關(guān)系式，從而實(shí)現(xiàn)一維、二維的耦合計(jì)算。這種方式每一時(shí)步都需完成迭代計(jì)算，計(jì)算量大但穩(wěn)定性較高，可以適用較大的時(shí)間步長(zhǎng)。