河流數(shù)學(xué)模型理論與實(shí)踐案例

1、河流數(shù)學(xué)模型理論與實(shí)踐案例一、案例案例性質(zhì)：一維非恒定水流模型案例背景：（1）以長(zhǎng)江為基本背景，構(gòu)建一條河長(zhǎng)100km的無(wú)支流入?yún)R河段，并給出斷面幾何形態(tài)和其他物理參數(shù)。（2）給定入流條件Q符合e指數(shù)分布律，平均值取Q0=4000m3/s；并根據(jù)斷面形態(tài)和河段坡度給出下游控制水位。結(jié)果要求：（1）依據(jù)Preissmann格式導(dǎo)出Holl-Cunge水動(dòng)力學(xué)流量演算方法；（2）給出Holl-Cunge水動(dòng)力學(xué)流量演算數(shù)值方法的三性分析，并給出該方法的收斂特性分析R1、R2圖；（3）利用馬斯京根法和Holl-Cunge方法，在上述案例背景上進(jìn)行流量演算，給出數(shù)值演算結(jié)果并對(duì)兩種方法的計(jì)算結(jié)果進(jìn)行對(duì)

2、比分析；（4）考慮基本方程試證明其為對(duì)流擴(kuò)散方程，理論闡明采用Holl-Cunge流量演算方法模擬該對(duì)流擴(kuò)散方程的技術(shù)方法。二、Holl-Cunge水動(dòng)力學(xué)流量演算方法 一維水流數(shù)學(xué)模型基于以下圣維南方程組：水流連續(xù)方程：（1）水流運(yùn)動(dòng)方程：（2）式中，為過(guò)水?dāng)嗝婷娣e，為流量，為水位，為重力加速度，為能坡，為流程方向，為時(shí)間。 假定在某一確定斷面上，流量和水位成一一對(duì)應(yīng)的單值關(guān)系（3）這里為過(guò)水?dāng)嗝婷娣e，為水深。若假定式（3）成立，則反函數(shù)存在，且也為單值函數(shù)，因此有（4）將式（4）代入連續(xù)方程式（1）有（5）這是一階的雙曲型方程式，因?yàn)樽兞縌的函數(shù)，故是擬線性的。在特征線上，有（6）（7）結(jié)

3、合式（5）有（8）由此可寫(xiě)出方程（5）的另一種形式為（9）在滿足式（6）的條件下，方程（9）的另一等價(jià)形式及相應(yīng)的解為 或 （10）這個(gè)解表明流量沿特征線族上為一常數(shù)，并不允許流量過(guò)程在任意時(shí)刻任何斷面上發(fā)生任何變形輸移。從這個(gè)結(jié)果可得出一個(gè)結(jié)論，流量水位成單值關(guān)系的假設(shè)，實(shí)質(zhì)上描述的是流量的非線性對(duì)流輸移，入流過(guò)程線自始至終將保持原來(lái)的形狀。 由于方程（9）或式（5）是非線性的，故很難得到解析解，可考慮求其數(shù)值解。采用四點(diǎn)偏心Preissmann格式逼近方程（9）中關(guān)于時(shí)空偏微商項(xiàng)，有（11）這里和為權(quán)重因子，和分別為節(jié)點(diǎn)函數(shù)值，式中的暫取為常數(shù)。其網(wǎng)格布置如下圖1所示。圖1 差分網(wǎng)格示意圖

4、令，并取，則從方程（11）解出為（12）式中（13）式（12）即為Holl-Cunge水動(dòng)力學(xué)流量演算方法。 在點(diǎn)（j，n）上按泰勒級(jí)數(shù)展開(kāi)方程（11）中各項(xiàng)節(jié)點(diǎn)函數(shù)值，并經(jīng)整理后得（14）略去的高階項(xiàng)得（15）式中（16）很明顯，當(dāng)和取任意值時(shí)，數(shù)值格式（12）實(shí)質(zhì)上以二階近似逼近于一個(gè)對(duì)流擴(kuò)散方程（15），而不是純對(duì)流方程（9）；只有當(dāng)時(shí)，才有，差分格式（12）才完全逼近微分方程（9）。 考慮連續(xù)方程式（1）和忽略時(shí)空慣性項(xiàng)的運(yùn)動(dòng)方程式（2），有（17）（18）式中，為底坡，為流量系數(shù)，為謝才系數(shù)。 當(dāng)河寬沿程不變時(shí)，有，則（19）為便于分析，對(duì)式（17）關(guān)于x求偏微商（20）再次對(duì)方程（

5、18）關(guān)于t求偏微商得（21）將方程（20）和方程（21）相加得（22）對(duì)K關(guān)于t求偏微商并使用方程（17）有（23）將式（23）代入式（22）化成一維線性拋物型方程（24） 要使差分格式（12）逼近物理擴(kuò)散方程（24），則必須滿足（25）和（26） 由此看出，忽略運(yùn)動(dòng)方程時(shí)空慣性項(xiàng)的圣維南方程組，在假定流量系數(shù)K是水深的單值函數(shù)的條件下，描述了沿程坦化的洪水?dāng)U散波，而不是運(yùn)動(dòng)波，且按式（26）選取數(shù)值計(jì)算參數(shù)，和，時(shí)可由差分方程（12）逼近。三、Holl-Cunge方法的三性分析 在用差分方程代替微分方程問(wèn)題，將求解微分方程化為求代數(shù)方程的離散解時(shí)，必須進(jìn)一步對(duì)差分方程的收斂性、穩(wěn)定性和相容

6、性進(jìn)行研究。（1）相容性 將一個(gè)微分方程用差分格式化為相應(yīng)差分方程，當(dāng)步長(zhǎng)與趨近于零時(shí)，這個(gè)微分方程應(yīng)當(dāng)收斂于原微分方程，也就是說(shuō)，相應(yīng)的差分方程和微分方程之間的截?cái)嗾`差在任一時(shí)刻任一網(wǎng)格點(diǎn)上均應(yīng)趨近于零，這樣的差分方程和微分方程才是相容的。 對(duì)于上節(jié)得到的差分方程（11），假定的結(jié)點(diǎn)值是微分方程（9）的解（27）將，和在結(jié)點(diǎn)（j，n）處展開(kāi)為泰勒級(jí)數(shù)有（28）（29）（30）并注意到（31）將式（28）（30）代入差分方程（11）有（32）式中，因在點(diǎn)（j，n）處滿足式（27），即（33）則有截?cái)嗾`差（34）由此可看出，用差分方程（11）逼近微分方程（9）的截?cái)嗾`差的階是：當(dāng)和時(shí)，差分格式具

7、有一階精度；當(dāng)和時(shí)，具有二階精度；當(dāng)和時(shí)，格式具有高階精度。當(dāng)與趨近于零時(shí)有（35）即當(dāng)步長(zhǎng)與趨近于零時(shí)，截?cái)嗾`差趨近于零。可見(jiàn)，差分格式（11）與微分方程（9）是相容的。（2）穩(wěn)定性 差分法計(jì)算中所產(chǎn)生的誤差（舍入誤差，參數(shù)誤差等）隨時(shí)間衰減或不增大，則稱(chēng)離散格式是穩(wěn)定的，反之，則是不穩(wěn)定的。分析差分方程穩(wěn)定性有不同的方法，如矩陣方法，諧波分析法等。 由于偏微分方程（9）是線性的，其解可用Fourier級(jí)數(shù)表示，解的特性可由特解表征。不妨取Fourier級(jí)數(shù)解的第m階分量為（36）這樣，差分方程（11）的特解在網(wǎng)格結(jié)點(diǎn)上的值為：（37）（38）（39）（40）上幾式中，為幅常數(shù)，與x和t無(wú)關(guān)

8、，L為波長(zhǎng)，T為周期，將式（36）（39）代入差分方程（11）得（41）即（42）式中，為未知量結(jié)點(diǎn)值，令（43）則有（44）因，有（45）得（46）記（47）式（46）可寫(xiě)成（48）由此可解出（49）亦即（50）記（51）因（52）故式（50）變?yōu)椋?3）設(shè)取，和分別表示波頻率的實(shí)部和虛部，可寫(xiě)成（54）比較式（53）與式（54）的實(shí)部和虛部得：（55）（56）上兩式兩邊平方相加得（57）式中， 數(shù)值穩(wěn)定要求幅因子，則必有，亦即（58）顯然如果同時(shí)滿足以下不等式（59）數(shù)值計(jì)算是穩(wěn)定的。由不等式得（60）解得（61）同樣由不等式可解得式（61），亦即只要按式（61）選擇時(shí)空權(quán)重因子和，按普列

9、斯曼格式計(jì)算是數(shù)值穩(wěn)定的。當(dāng)時(shí)，穩(wěn)定條件是，而當(dāng)時(shí)，則要求。對(duì)于和取任意值，穩(wěn)定性與水流流動(dòng)方向有關(guān)，穩(wěn)定域如圖2所示。顯然只要取和，條件（61）是得到滿足的，數(shù)值計(jì)算無(wú)條件穩(wěn)定。1.0無(wú)條件穩(wěn)定條件穩(wěn)定1.0條件穩(wěn)定無(wú)條件穩(wěn)定0.5條件穩(wěn)定無(wú)條件不穩(wěn)定0.5無(wú)條件不穩(wěn)定條件穩(wěn)定00.51.000.51.0 （a） （b）圖2 穩(wěn)定域圖（3）收斂性 若求得的差分方程精確解在所考慮區(qū)域內(nèi)的任意點(diǎn)（j，n）上，當(dāng)離散的網(wǎng)格步長(zhǎng)趨于零時(shí)，差分方程精確解是否趨近于微分方程的解，這就是收斂性問(wèn)題。解析解的幅和相 將解式（36）代入微分方程（9）得（62）因此微分方程解的波速為（63）另一方面，微分方程（

10、9）的特征速度可通過(guò)特征線理論求得，它取決于系數(shù)，而系數(shù)是實(shí)數(shù)，且波長(zhǎng)L和周期T也均為實(shí)數(shù)。因此，與的關(guān)系為實(shí)數(shù)范圍的單一關(guān)系，這樣可得到解式（36）的波幅因子的模為（64）它表明由微分方程（9）決定的解將不存在增減幅現(xiàn)象。數(shù)值解的幅和相由式（55）和式（56）相除得（65）可解得數(shù)值解波頻率為（66）進(jìn)而數(shù)值波速可表示成（67）解式（57）表示數(shù)值解幅因子，而解式（67）則表示數(shù)值解波速，根據(jù)萊恩收斂準(zhǔn)則，得兩個(gè)收斂因子為（68）（69）根據(jù)以上兩式，以，和為參數(shù)，可計(jì)算收斂因子和與相對(duì)波長(zhǎng)的關(guān)系曲線（圖3圖3），從圖中可以看出：對(duì)于，當(dāng)時(shí)，格式不存在減幅和相錯(cuò)位；當(dāng)時(shí)，格式存在相錯(cuò)位。當(dāng)?shù)?/p>

11、取值小于0.5，取值大于0.5時(shí)，對(duì)于短波模擬，取值不宜大于10。而對(duì)于長(zhǎng)波模擬，值可允許取大一些也能保證數(shù)值計(jì)算精度。對(duì)于確定的波長(zhǎng)，和值愈接近于1，結(jié)果愈不理想，數(shù)值減幅將較嚴(yán)重。對(duì)于河床變形模擬，由于河床變形速度較小，從物理過(guò)程上講，可選得較大。對(duì)，如果取得過(guò)小，致使時(shí)，數(shù)值計(jì)算將出現(xiàn)數(shù)值擾動(dòng)現(xiàn)象。圖3收斂因子與相對(duì)波長(zhǎng)的關(guān)系曲線（，）圖4收斂因子與相對(duì)波長(zhǎng)的關(guān)系曲線（，）四、流量演算（1）案例背景 以長(zhǎng)江為基本背景，構(gòu)建一條河長(zhǎng)100km的無(wú)支流入?yún)R河段，并給出斷面幾何形態(tài)和其他物理參數(shù)。給定入流條件Q符合e指數(shù)分布律，平均值取Q0=4000m3/s；并根據(jù)斷面形態(tài)和河段坡度給出下游控

12、制水位。斷面布置及斷面形態(tài) 河段坡度取0.001，全長(zhǎng)100km，全河段共布置100個(gè)斷面，斷面間距為1km。假設(shè)斷面為矩形斷面，寬1000m，高50m。進(jìn)口斷面為x=0km，出口斷面為x=100km。進(jìn)口流量 進(jìn)口流量過(guò)程為（70）其中，流量單位為m3/s，均值為4000 m3/s。進(jìn)口流量過(guò)程如圖5所示。圖5 進(jìn)口流量過(guò)程線下游控制水位 根據(jù)斷面形態(tài)及河段坡度，取下游控制水位為常數(shù)20.0m。初始流量 給定初始時(shí)刻（t=0h時(shí)）各斷面流量相等，為t=0h時(shí)的進(jìn)口流量。初始時(shí)刻各斷面流量分布如圖6所示。圖6 初始時(shí)刻流量沿程分布（2）馬斯京根法 將連續(xù)方程（1）寫(xiě)成（71）式中，為微小河段內(nèi)

13、河槽的槽蓄量，為微小河段內(nèi)入流量，為出流量。 為了求解方程（71），在馬斯京根法中不是直接求解運(yùn)動(dòng)方程（2），而是引進(jìn)新的假設(shè)條件來(lái)使問(wèn)題有解。 假定槽蓄量?jī)H是入流與出流量的線性關(guān)系，即（72）式中，和為經(jīng)驗(yàn)系數(shù)，可以通過(guò)實(shí)測(cè)資料由試錯(cuò)法確定。將方程（71）寫(xiě)成增量形式（73）式中（74）按式（72）可寫(xiě)出（75）（76）將式（74）、（75）、（76）代入式（73）可得（77）式中（78）式中，和分別為計(jì)算時(shí)段始末河段的入流量，和分別為相應(yīng)的出流量。這樣，出流斷面上n+1時(shí)刻上的流量就可按式（77）計(jì)算。由于在使用式（77）中必需先已知入流斷面上的流量過(guò)程線（和為已知）時(shí)才能演算出流斷面與同

14、時(shí)刻的流量。這種演算方法在水文學(xué)中稱(chēng)為無(wú)預(yù)見(jiàn)期的演算法。 根據(jù)馬斯京根方法，取時(shí)間步長(zhǎng)，經(jīng)驗(yàn)系數(shù)，給定初始時(shí)刻各斷面的流量與進(jìn)口流量一致，進(jìn)行流量演算，結(jié)果如下圖7和圖8所示。 圖7給出了沿程各斷面流量隨時(shí)間變化圖，從圖中可以看出，隨著沿程距離的加大，出現(xiàn)波動(dòng)的時(shí)間相應(yīng)越遲，x=10km斷面約在t=2h時(shí)就出現(xiàn)波動(dòng)，x=20km處大約在t=3h時(shí)出現(xiàn)波動(dòng)，x=30km處大約在t=4h時(shí)出現(xiàn)波動(dòng)，x=40km處大約在t=6h時(shí)出現(xiàn)波動(dòng)，x=50km處大約在t=7h時(shí)出現(xiàn)波動(dòng)，x=60km處大約在t=8h時(shí)出現(xiàn)波動(dòng)，x=70km處大約在t=10h時(shí)出現(xiàn)波動(dòng)，x=80km處大約在t=11h時(shí)出現(xiàn)波動(dòng)

15、，x=90km處大約在t=12h時(shí)出現(xiàn)波動(dòng)，x=100km處大約在t=14h時(shí)出現(xiàn)波動(dòng)。相應(yīng)地，其流量達(dá)到峰值的時(shí)間也隨著沿程距離的增大而延遲。圖7 各斷面流量隨時(shí)間變化 圖8給出了各時(shí)刻流量沿程分布情況，從圖中可以看出，隨著時(shí)間推進(jìn)，進(jìn)口水波逐漸向前推進(jìn)，在t=4h時(shí)水波大約傳播至x=40km處；在t=8h時(shí)水波大約傳播至x=60km處；在t=12h時(shí)水波大約傳播至x=90km處，此時(shí)，進(jìn)口流量達(dá)到峰值；在t=16h時(shí)，水波已傳播至x=100km處，流量峰值傳播至x=30km附近；在t=20h時(shí)，流量峰值傳播至x=60km附近；在t=24h時(shí)，流量峰值傳播至x=80km附近。洪水波的流量峰值

16、基本保持在入口流量峰值4000m3/s左右，但洪水波的形狀則隨著時(shí)間的推進(jìn)逐漸變緩，說(shuō)明馬斯金根法流量演算方法計(jì)算流量時(shí)存在洪水波的坦化現(xiàn)象。圖8 各時(shí)刻流量沿程分布（3）Holl-Cunge方法 根據(jù)Holl-Cunge方法，取時(shí)間步長(zhǎng)，經(jīng)驗(yàn)系數(shù)，給定初始時(shí)刻各斷面的流量與進(jìn)口流量一致，進(jìn)行流量演算，結(jié)果如下圖9和圖10所示。 圖9給出了沿程各斷面流量隨時(shí)間變化圖，從圖中可以看出，隨著沿程距離的加大，出現(xiàn)波動(dòng)的時(shí)間相應(yīng)越遲，x=10km斷面約在t=3h時(shí)出現(xiàn)波動(dòng)，x=20km處大約在t=5h時(shí)出現(xiàn)波動(dòng)，x=30km處大約在t=8h時(shí)出現(xiàn)波動(dòng)，x=40km處大約在t=10h時(shí)出現(xiàn)波動(dòng)，x=50

17、km處大約在t=12h時(shí)出現(xiàn)波動(dòng)，x=60km處大約在t=15h時(shí)出現(xiàn)波動(dòng)，x=70km處大約在t=18h時(shí)出現(xiàn)波動(dòng)，x=80km處大約在t=20h時(shí)出現(xiàn)波動(dòng)，x=90km處大約在t=22h時(shí)出現(xiàn)波動(dòng)，x=100km處在024h內(nèi)幾乎未出現(xiàn)波動(dòng)，說(shuō)明水波還沒(méi)傳播到此處。相應(yīng)地，其流量達(dá)到峰值的時(shí)間也隨著沿程距離的增大而延遲。圖9 各斷面流量隨時(shí)間變化 圖10給出了各時(shí)刻流量沿程分布情況，從圖中可以看出，隨著時(shí)間推進(jìn)，進(jìn)口水波逐漸向前推進(jìn)，在t=4h時(shí)水波大約傳播至x=20km處；在t=8h時(shí)水波大約傳播至x=30km處；在t=12h時(shí)水波大約傳播至x=50km處，此時(shí)，進(jìn)口流量達(dá)到峰值；在t=

18、16h時(shí)，水波已傳播至x=60km處，流量峰值傳播至x=20km附近；在t=20h時(shí)，水波已傳播至x=80km處，流量峰值傳播至x=30km附近；在t=24h時(shí)，水波已傳播至x=90km處，流量峰值傳播至x=50km附近。圖10 各時(shí)刻流量沿程分布（4）計(jì)算結(jié)果對(duì)比 根據(jù)流量演算結(jié)果，對(duì)馬斯京根法和Holl-Cunge方法的計(jì)算結(jié)果進(jìn)行對(duì)比，對(duì)比結(jié)果如下圖11及圖12所示。 圖11給出了各斷面流量隨時(shí)間變化對(duì)比圖，從圖中可以看出，各斷面，馬斯京根法計(jì)算結(jié)果出現(xiàn)流量波動(dòng)的時(shí)間較Holl-Cunge方法的早，說(shuō)明馬斯京根法計(jì)算的流量演進(jìn)速度大于Holl-Cunge方法計(jì)算的流量演進(jìn)速度，隨著沿程距離的推進(jìn)，兩者計(jì)算的洪水波的演進(jìn)時(shí)間