河流数学模型理论与实践案例

1、河流数学模型理论与实践案例一、案例案例性质：一维非恒定水流模型案例背景：（1）以长江为基本背景，构建一条河长100km的无支流入汇河段，并给出断面几何形态和其他物理参数。（2）给定入流条件Q符合e指数分布律，平均值取Q0=4000m3/s；并根据断面形态和河段坡度给出下游控制水位。结果要求：（1）依据Preissmann格式导出Holl-Cunge水动力学流量演算方法；（2）给出Holl-Cunge水动力学流量演算数值方法的三性分析，并给出该方法的收敛特性分析R1、R2图；（3）利用马斯京根法和Holl-Cunge方法，在上述案例背景上进行流量演算，给出数值演算结果并对两种方法的计算结果进行对

2、比分析；（4）考虑基本方程试证明其为对流扩散方程，理论阐明采用Holl-Cunge流量演算方法模拟该对流扩散方程的技术方法。二、Holl-Cunge水动力学流量演算方法 一维水流数学模型基于以下圣维南方程组：水流连续方程：（1）水流运动方程：（2）式中，为过水断面面积，为流量，为水位，为重力加速度，为能坡，为流程方向，为时间。 假定在某一确定断面上，流量和水位成一一对应的单值关系（3）这里为过水断面面积，为水深。若假定式（3）成立，则反函数存在，且也为单值函数，因此有（4）将式（4）代入连续方程式（1）有（5）这是一阶的双曲型方程式，因为变量Q的函数，故是拟线性的。在特征线上，有（6）（7）结

3、合式（5）有（8）由此可写出方程（5）的另一种形式为（9）在满足式（6）的条件下，方程（9）的另一等价形式及相应的解为 或 （10）这个解表明流量沿特征线族上为一常数，并不允许流量过程在任意时刻任何断面上发生任何变形输移。从这个结果可得出一个结论，流量水位成单值关系的假设，实质上描述的是流量的非线性对流输移，入流过程线自始至终将保持原来的形状。 由于方程（9）或式（5）是非线性的，故很难得到解析解，可考虑求其数值解。采用四点偏心Preissmann格式逼近方程（9）中关于时空偏微商项，有（11）这里和为权重因子，和分别为节点函数值，式中的暂取为常数。其网格布置如下图1所示。图1 差分网格示意图

4、令，并取，则从方程（11）解出为（12）式中（13）式（12）即为Holl-Cunge水动力学流量演算方法。 在点（j，n）上按泰勒级数展开方程（11）中各项节点函数值，并经整理后得（14）略去的高阶项得（15）式中（16）很明显，当和取任意值时，数值格式（12）实质上以二阶近似逼近于一个对流扩散方程（15），而不是纯对流方程（9）；只有当时，才有，差分格式（12）才完全逼近微分方程（9）。 考虑连续方程式（1）和忽略时空惯性项的运动方程式（2），有（17）（18）式中，为底坡，为流量系数，为谢才系数。 当河宽沿程不变时，有，则（19）为便于分析，对式（17）关于x求偏微商（20）再次对方程（

5、18）关于t求偏微商得（21）将方程（20）和方程（21）相加得（22）对K关于t求偏微商并使用方程（17）有（23）将式（23）代入式（22）化成一维线性抛物型方程（24） 要使差分格式（12）逼近物理扩散方程（24），则必须满足（25）和（26） 由此看出，忽略运动方程时空惯性项的圣维南方程组，在假定流量系数K是水深的单值函数的条件下，描述了沿程坦化的洪水扩散波，而不是运动波，且按式（26）选取数值计算参数，和，时可由差分方程（12）逼近。三、Holl-Cunge方法的三性分析 在用差分方程代替微分方程问题，将求解微分方程化为求代数方程的离散解时，必须进一步对差分方程的收敛性、稳定性和相容

6、性进行研究。（1）相容性 将一个微分方程用差分格式化为相应差分方程，当步长与趋近于零时，这个微分方程应当收敛于原微分方程，也就是说，相应的差分方程和微分方程之间的截断误差在任一时刻任一网格点上均应趋近于零，这样的差分方程和微分方程才是相容的。 对于上节得到的差分方程（11），假定的结点值是微分方程（9）的解（27）将，和在结点（j，n）处展开为泰勒级数有（28）（29）（30）并注意到（31）将式（28）（30）代入差分方程（11）有（32）式中，因在点（j，n）处满足式（27），即（33）则有截断误差（34）由此可看出，用差分方程（11）逼近微分方程（9）的截断误差的阶是：当和时，差分格式具

7、有一阶精度；当和时，具有二阶精度；当和时，格式具有高阶精度。当与趋近于零时有（35）即当步长与趋近于零时，截断误差趋近于零。可见，差分格式（11）与微分方程（9）是相容的。（2）稳定性 差分法计算中所产生的误差（舍入误差，参数误差等）随时间衰减或不增大，则称离散格式是稳定的，反之，则是不稳定的。分析差分方程稳定性有不同的方法，如矩阵方法，谐波分析法等。 由于偏微分方程（9）是线性的，其解可用Fourier级数表示，解的特性可由特解表征。不妨取Fourier级数解的第m阶分量为（36）这样，差分方程（11）的特解在网格结点上的值为：（37）（38）（39）（40）上几式中，为幅常数，与x和t无关

8、，L为波长，T为周期，将式（36）（39）代入差分方程（11）得（41）即（42）式中，为未知量结点值，令（43）则有（44）因，有（45）得（46）记（47）式（46）可写成（48）由此可解出（49）亦即（50）记（51）因（52）故式（50）变为（53）设取，和分别表示波频率的实部和虚部，可写成（54）比较式（53）与式（54）的实部和虚部得：（55）（56）上两式两边平方相加得（57）式中， 数值稳定要求幅因子，则必有，亦即（58）显然如果同时满足以下不等式（59）数值计算是稳定的。由不等式得（60）解得（61）同样由不等式可解得式（61），亦即只要按式（61）选择时空权重因子和，按普列

9、斯曼格式计算是数值稳定的。当时，稳定条件是，而当时，则要求。对于和取任意值，稳定性与水流流动方向有关，稳定域如图2所示。显然只要取和，条件（61）是得到满足的，数值计算无条件稳定。1.0无条件稳定条件稳定1.0条件稳定无条件稳定0.5条件稳定无条件不稳定0.5无条件不稳定条件稳定00.51.000.51.0 （a） （b）图2 稳定域图（3）收敛性 若求得的差分方程精确解在所考虑区域内的任意点（j，n）上，当离散的网格步长趋于零时，差分方程精确解是否趋近于微分方程的解，这就是收敛性问题。解析解的幅和相 将解式（36）代入微分方程（9）得（62）因此微分方程解的波速为（63）另一方面，微分方程（

10、9）的特征速度可通过特征线理论求得，它取决于系数，而系数是实数，且波长L和周期T也均为实数。因此，与的关系为实数范围的单一关系，这样可得到解式（36）的波幅因子的模为（64）它表明由微分方程（9）决定的解将不存在增减幅现象。数值解的幅和相由式（55）和式（56）相除得（65）可解得数值解波频率为（66）进而数值波速可表示成（67）解式（57）表示数值解幅因子，而解式（67）则表示数值解波速，根据莱恩收敛准则，得两个收敛因子为（68）（69）根据以上两式，以，和为参数，可计算收敛因子和与相对波长的关系曲线（图3图3），从图中可以看出：对于，当时，格式不存在减幅和相错位；当时，格式存在相错位。当的

11、取值小于0.5，取值大于0.5时，对于短波模拟，取值不宜大于10。而对于长波模拟，值可允许取大一些也能保证数值计算精度。对于确定的波长，和值愈接近于1，结果愈不理想，数值减幅将较严重。对于河床变形模拟，由于河床变形速度较小，从物理过程上讲，可选得较大。对，如果取得过小，致使时，数值计算将出现数值扰动现象。图3收敛因子与相对波长的关系曲线（，）图4收敛因子与相对波长的关系曲线（，）四、流量演算（1）案例背景 以长江为基本背景，构建一条河长100km的无支流入汇河段，并给出断面几何形态和其他物理参数。给定入流条件Q符合e指数分布律，平均值取Q0=4000m3/s；并根据断面形态和河段坡度给出下游控

12、制水位。断面布置及断面形态 河段坡度取0.001，全长100km，全河段共布置100个断面，断面间距为1km。假设断面为矩形断面，宽1000m，高50m。进口断面为x=0km，出口断面为x=100km。进口流量 进口流量过程为（70）其中，流量单位为m3/s，均值为4000 m3/s。进口流量过程如图5所示。图5 进口流量过程线下游控制水位 根据断面形态及河段坡度，取下游控制水位为常数20.0m。初始流量 给定初始时刻（t=0h时）各断面流量相等，为t=0h时的进口流量。初始时刻各断面流量分布如图6所示。图6 初始时刻流量沿程分布（2）马斯京根法 将连续方程（1）写成（71）式中，为微小河段内

13、河槽的槽蓄量，为微小河段内入流量，为出流量。 为了求解方程（71），在马斯京根法中不是直接求解运动方程（2），而是引进新的假设条件来使问题有解。 假定槽蓄量仅是入流与出流量的线性关系，即（72）式中，和为经验系数，可以通过实测资料由试错法确定。将方程（71）写成增量形式（73）式中（74）按式（72）可写出（75）（76）将式（74）、（75）、（76）代入式（73）可得（77）式中（78）式中，和分别为计算时段始末河段的入流量，和分别为相应的出流量。这样，出流断面上n+1时刻上的流量就可按式（77）计算。由于在使用式（77）中必需先已知入流断面上的流量过程线（和为已知）时才能演算出流断面与同

14、时刻的流量。这种演算方法在水文学中称为无预见期的演算法。 根据马斯京根方法，取时间步长，经验系数，给定初始时刻各断面的流量与进口流量一致，进行流量演算，结果如下图7和图8所示。 图7给出了沿程各断面流量随时间变化图，从图中可以看出，随着沿程距离的加大，出现波动的时间相应越迟，x=10km断面约在t=2h时就出现波动，x=20km处大约在t=3h时出现波动，x=30km处大约在t=4h时出现波动，x=40km处大约在t=6h时出现波动，x=50km处大约在t=7h时出现波动，x=60km处大约在t=8h时出现波动，x=70km处大约在t=10h时出现波动，x=80km处大约在t=11h时出现波动

15、，x=90km处大约在t=12h时出现波动，x=100km处大约在t=14h时出现波动。相应地，其流量达到峰值的时间也随着沿程距离的增大而延迟。图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左右，但洪水波的形状则随着时间的推进逐渐变缓，说明马斯金根法流量演算方法计算流量时存在洪水波的坦化现象。图8 各时刻流量沿程分布（3）Holl-Cunge方法 根据Holl-Cunge方法，取时间步长，经验系数，给定初始时刻各断面的流量与进口流量一致，进行流量演算，结果如下图9和图10所示。 图9给出了沿程各断面流量随时间变化图，从图中可以看出，随着沿程距离的加大，出现波动的时间相应越迟，x=10km断面约在t=3h时出现波动，x=20km处大约在t=5h时出现波动，x=30km处大约在t=8h时出现波动，x=40km处大约在t=10h时出现波动，x=50

17、km处大约在t=12h时出现波动，x=60km处大约在t=15h时出现波动，x=70km处大约在t=18h时出现波动，x=80km处大约在t=20h时出现波动，x=90km处大约在t=22h时出现波动，x=100km处在024h内几乎未出现波动，说明水波还没传播到此处。相应地，其流量达到峰值的时间也随着沿程距离的增大而延迟。图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）计算结果对比 根据流量演算结果，对马斯京根法和Holl-Cunge方法的计算结果进行对比，对比结果如下图11及图12所示。 图11给出了各断面流量随时间变化对比图，从图中可以看出，各断面，马斯京根法计算结果出现流量波动的时间较Holl-Cunge方法的早，说明马斯京根法计算的流量演进速度大于Holl-Cunge方法计算的流量演进速度，随着沿程距离的推进，两者计算的洪水波的演进时间差别越大。图11 各断面流量随时间变化对比 图12给出了各时刻流量沿程分布对比图，从图中可以看出，各时刻马斯京根法计算结果流量波传播的距离较Holly-Cunge方法的长，也说明了马斯京根法计算的流量演进速度大于Holl-Cunge方法计算的流量演进速度。另一方面，从图中可以看出，马斯京根法流量演算的洪水波较Holl-Cunge方法的演算结果更为平缓，说明马斯京根法计算的洪水波存在坦化现象，而Holl-Cunge方法则较完整地保持了原有的波形。图12 各时刻流量沿程分布对比五、结论 本文依据Preissmann格式理论推导了H