非恒定流泥沙数学模型原理及其应用
摘要：本文介绍了一维非恒定泥沙数学模型的原理和方法，并且对于泥沙数学模型中的一些关键技术等问题进行了较详细的讨论。利用已建立的一维非恒定泥沙数学模型，对黄河小北干流1981年洪水资料进行了模拟计算，计算结果表明：黄河干流、渭河、北洛河计算水位、流量过程线与实测值符合良好，可用于黄河的洪水演进计算。
关键词：非恒定流 泥沙数学模型 河床变形
1 泥沙数学模型基本方程
明渠或天然河流常被考虑作为一维流动，根据洪水波运动的圣维南方程、泥沙连续方程和泥沙扩散方程，可以简化推导出一维非恒定流泥沙数学模型的基本方程，其形式为水流连续方程
| (1) |
水流动量方程 |
|
| (2) |
泥沙连续方程 |
|
| (3) |
不平衡输沙方程 |
|
| (4) |
挟沙力方程 |
|
S*=f(Q,A,B,ω,S) | (5) |
其中 Q为流量，A为断面面积，B为断面宽度，Z为水位，K为流量模数，S为断面平均含沙量，S*为断面平均挟沙力，g为重力加速度，α为泥沙非平衡恢复饱和系数，ω为泥沙颗粒沉速，Ad为断面冲淤面积，γ′s为泥沙干容重。
2 计算方法简介
一维非恒定泥沙数学模型的计算采用非耦合方法，首先求解水流连续方程和动量方程，然后求解水流挟沙力、泥沙不平衡输沙方程和泥沙连续方程，具体求解过程如下。
2.1 水流方程的求解
首先利用Preissmann隐式差分格式将水流连续方程和动量方程转化为差分方程，然后对差分方程进行线性化，在线性化过程中，略去增量的乘积项，最后得到以下线性方程组
A1jΔQj+B1jΔZj+C1jΔQj+1+D1jΔZj+1=E1j | (6) |
A2jΔQj+B2jΔZj+C2jΔQj+1+D2jΔZj+1=E2j | (7) |
其中 Aij、Bij、Cij、Dij、Eij（i=1,2）为第j单元河段差分方程的系数（j=1,2,......,N-1,其中N为断面个数）。
给定边界条件
ΔQ1=Q1 n+1-Qn1=Q1(tn+1)-Q1n | (8) |
ΔZN=ZNn+1-ZNn=ZN(tn+1)-ZNn | (9) |
方程(6)、(7)式及边界条件，共有2N个未知数，2N个方程，可以求解。由于差分方程中的系数包含有未知数，方程求解不能直接求出未知变量，因此方程求解时必须进行迭代处理。下面给出用追赶法求解的步骤，追赶方程为
ΔQj=FjΔZj+Gj | (10) |
ΔZj=HjΔQj+1+IjΔZj+1+Jj | (11) |
其中Hj、Ij、Jj、，Fj、Gj为追赶系数。
2.2 水流挟沙力、动床阻力、河宽变化、床沙级配调整采用[2]的方法
2.3 不平衡输沙方程求解
利用迎风格式，将(3)式离散为差分方程，整理后得
| (12) |
当Q≥0时，利用上边界条件，自上而下计算各断面含沙量；当Q<0时，利用下边界条件由下至上计算各断面含沙量。
2.4 河床变形及淤积量计算
由式(3)与式(4)相减可以得到河床变形方程为 |
|
| (13) |
将上式写成差分方程，各断面淤积面积为 |
|
| (14) |
第j河段的淤积量为 |
|
ΔWj=(ΔAj+ΔAj+1)Δxj/2 | (15) |
3 水流内边界的处理
水流内边界是指河道的几何形状的不连续或水力特性的不连续点。例如，河流的汇合点、河流分流、局部河段内生产堤决口等等。在这些内部边界处，圣维南方程组和单一河道泥沙不平衡输沙方程等都不再适用，必须根据其水力特性作特殊处理。内边界条件通常包含两个相容条件：即流量的连续条件和能量守恒条件（或动量守恒条件）。本模型主要考虑了以下几个类型的内边界处理。
3.1 水沙的汇入或汇出
如图1所示，假设汇入或汇出点上下断面满足以下条件
|
|
|
3.2 支流从干流分流
如图2所示，干流和支流上断面之间应满足连续方程和能量方程
|
|
|
|
|
|
图1 水沙的汇入或汇出示意图 | 图2 支流从干流分流示意图 | 图3 支流汇入干流示意图 |
支流从干流分流时，将干流分流断面按干支流流量比分为两部分，忽略时变项，将方程(4)式直接写成差分形式，求得干流和支流下游断面的含沙量为
|
|
其中φ为分流系数，由实测资料确定。
3.3 支流汇入干流
当支流汇入干流时，与支流从干流分流类似，干流和支流上断面之间也满足连续方程和能量方程
|
|
|
将汇流断面按干支流流量比分为两部分，忽略时变项，将方程(4)式直接写成差分形式，求得汇流断面的含沙量为
|
|
|
其中 φ为分流系数，由实测资料确定。
4 泥沙数学模型的应用
4.1 区域和时段
计算区域为龙门、华县、河津、状头至潼关河段，黄河干流上有渭河和汾河汇入，在渭河上有北洛河汇入。本模型同时模拟黄河干流、渭河、北洛河三条河流的洪水演进过程，各汇入点作为内边界处理，汾河仅作为已知水、沙过程线汇入黄河干流。上边界条件为龙门、华县、河津、?NFDA1?头各站，水沙为进口控制的已知条件，下边界条件为潼关站出口控制水位。
计算时段为1981年汛初第一场洪水，洪水时间为7月3号至7月14号，洪水持续时间为12天，龙门最大流量为6400m3/s、最大含沙量为298.0kg/m3；华县最大流量为970m3/s、最大含沙量为117.0kg/m3。
原始大断面资料采用1981年汛前实测大断面资料。
4.2 计算结果与实测值的比较
图4为潼关站计算流量与实测值的比较、计算含沙量与实测值的比较。由图4可以看出，计算的潼关出口流量过程线与实测过程线比较符合，计算洪峰最大值和相应洪峰传播时间与实测值比较接近;计算潼关出口含沙量过程线与实测值也是比较一致的。由图4分析得出，本模型能够比较好地模拟出三条河流的洪水传播过程和泥沙冲淤调整过程。
|
图4 潼关站流量、含沙量计算值与实测值的比较 |
Comparison of calculated and measured discharge and concentrations at Tongguan station |
|
图5 华阴站流量、水位计算值与实测值的比较 |
Comparison of calculated discharge and water levels with measured at Huayin station |
图5为华阴站计算水位与实测值的比较、计算流量与实测值的比较。由图中可以看出，计算水位、流量过程线与相应实测过程线比较符合。同时也看出，在洪水初期，华阴站出现了倒灌现象，流量出现了负值，本模型也比较好地模拟出了黄河干流倒灌渭河的现象。
图6为朝邑站计算水位与实测值的比较、计算流量与实测值的比较。可以看出，计算水位、流量过程线与相应实测过程线比较一致。由于计算北洛河各断面的水位和流量，是通过一系列内边界条件计算得到的，因此也可以说明本模型所采用的内边界处理技术是比较合理的。
5 结语 本文较详细地介绍了一维非恒定泥沙数学模型的原理和计算方法，并对泥沙数学模型中的一些关键技术等问题进行了讨论。利用一维非恒定泥沙数学模型，对黄河小北干流1981年洪水资料进行了模拟计算，结果表明：模型比较好地模拟了黄河干流、渭河、北洛河洪水传播特性，洪峰传播时间和洪峰最大值与实测值比较一致；模型对于汇流区水流和泥沙计算所采用的处理技术也是比较合理的，对于干流倒灌支流或支流倒灌干流的现象，模型也是比较能够适应的，因此，模型可用于多条河流的洪水演进计算。 |
|
参 考 文 献
[1] Lyn,D.A. Goodwin,P. Stabling of a eneral Preissmann scheme. J. of Hydra. Eng. Vol. 113,No.1,1987.
[2] 梁国亭，张仁。黄河小北干流一维泥沙冲淤数学模型，人民黄河，1996年9月。