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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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