摘要： 采用有限元方法 建立起一套河道平面二维水流泥 沙数学模型。在前人研究 的基础上，采用了质量集中的处理方法，提出了压缩存储的方法， 从而大大减少了计算 存储量。针对有限元法时间步长需取得较短问题 ，采用了“预报-校正-迭代”的算法，提出了“非恒定-恒定-非恒定流”的算法，既能解决工程实际问题，又大大减少了计算量。作者以下荆江监利河段为例进行泥沙冲淤计算，计算结果与实测值符合较好，从而证明了模型的可靠性。
关键词： 水流泥沙 有限元 模型验证
三峡工程建成后，水库将拦蓄大量泥沙，下泄水流含沙量减小，对三峡工程坝下游河道将产生以冲刷为主的影响 ，包括对荆江河段的河势及荆江大堤带来影响。为研究坝下游重点河段的河床冲淤分布、河势变化、近岸流速变化等问题，一维模型显得无能为力，但可采用平面二维模型来解决。有限元方法可采用无结构 化网格，能很好地模拟不规则的几何形状，因此很适合于对天然河道的模拟。
然而，正如其它方法一样，有限元法也有它的缺点，主要是计算存储量和运算量较大。为扬长避短，使有限元方法能运用到对天然河道的模拟上来，本模型运用质量集中［4 ］ 的方法将系数矩阵转化为三对角矩阵，并提出了紧凑的分块压缩存储方法，从而大大减少了计算存储量，使得计算能在一般微机上进行。采用质量集中方法的不足之处是时间步长需取得较短，且在河道模拟中尤为突出( 因河道比较窄长，网格需划分很细，而该法的稳定性要 求时间步长与网格尺度成正比) 。针对该问题，笔者采用了“ 预报- 校正- 迭代［5 ］ ” 的算法，该法可加大时间步长，同时有效避免了数值震荡。针对长系列水沙条件下计算量较大问题，作者又提出了“ 非恒定- 恒定- 非恒定流” 的算法，该算法既能解决工程实际问 题，又大大减少了计算量，使有限元方法能够很好地运用于河道水流泥沙问题的实际计算。
1 基本方程
平面二维水流方程
(1)
(2)
(3)
悬移质泥沙扩散方程
(4)
推移质不平衡输移方程［6］
(5)
河床变形方程
由悬移质引起的河床变形方程为
<![endif]> (6)
由推移质引起的河床变形方程为
(7)
以上各式中U ，V 分别为垂线平均流速在x,y 方向上的分量；Zs 、Zb 和H 分别为水位、河底高程和水深；g 为重力加速度；vt 为水流紊动粘性系数；ρ 为水的密度；τx 、τy 、舄瓂分别为底部切应力在x 和y 方(τx 、τy )=，向上的分量：C 为谢才系数，常用曼宁公式计算：C=H1/6 /n ；S 和S* 分别为垂线平均含沙量和挟沙力；N 和N* 分别为推移质输沙量和推移质输沙能力折算成全水深的泥沙浓度；εs 为泥沙紊动扩散系数；ω 为泥沙沉速；γ′ 为床沙干容重；α 为悬移质泥沙恢复饱和系数，淤积时取0.25 ，冲刷时取1.0 ；β 为推移质泥沙恢复饱和系数，取0.25 。
悬移质泥沙共分成8组，水流挟沙力和分组挟沙力级配采用李义天方法［7］ 进行计算。
2 有限元方程
将整个计算域剖分成一个三角形网格系统。每个三角形为一个单元，其编号为e,e=1,2,3 ， … ，NE,NE 为单元总数。单元三个顶点为节点，其局部编号为j=1,2,3( 以逆时针为序) 。节点的整体编号为i,i=1,2,3,… ，NP ，NP 为节点总数。节点的整体编号与局部编号计算前一定要规定好。
引入插值函数：f=fi (t)φi(x,y),φ 为形函数。对方程(1) ～(5) 中的变量用插值函数近似表示，并使用伽辽金法［1 ］ 对方程进行整理变形，可得到积分方程
(8)
(9)
（10）
(11)
(12)
由此，经整理得到如下有限元方程
(13)
(14)
(15)
(16)
(17)
其中
i,j,k=1,2,3,… ，NP.
3 数值解法
为书写方便，采用通用变量P 来代替方程(13) ～(17) 中的变量Zs 、U 、V 、S 和N ，用FP 表示方程中等号右边项，并用对角矩阵 Aij ，则方程可化成统一的形式： dPi /dt=FPi 。
对该方程的求解，模型采用了“ 预报- 校正- 迭代［5 ］ ” 的计算 方法 。其方法为：用二阶显式Adams 公式作为预测公式，梯形公式( 隐式) 作校正公式，可构造Adams 二阶PC 公式。离散后的方程不需联接，计算过程稳定性好，时间步长可取得较长，同时还有效避免了数 值振荡。具体计算过程如下
预报 (18)
校正 (19)
迭代给定误差洌杂谒械1≤i≤NP ，若 ，则令： ，否则令： ；转到校正，继续迭代，直到满足精度要求为止。当求出Hi 、Ui 、Vi 、Si 和Ni 后，代入河床变形方程便可求得冲淤变形后的河床高程Zn 1 bi ，这样便完成了一个时段的计算。以上角标n 表时段，* 表预报值，** 表校正值，无* 标记的表示该时段的计算结果。
在水位求解过程中还采用了Kawahara“ 选择系数集中［4 ］ ” 的经验处理方法，引入了选择性集中系数： α=0.8 ～1.0 ， 其作用相当于对计算结果进行“ 光滑” 处理。
由于系数矩阵为NP×NP 阶矩阵，若直接进行存储会占用大量的内存，从而导致计算无法进行。笔者发现系数矩阵中绝大部分为零元素，为此提出了压缩存储的方法，即只对系数矩阵中的非零元素进行存储和计算，同时用辅助数组指明相应非零元素所在的位置。
对于长系列水沙条件的计算，提出了“ 非恒定- 定非- 定流” 的算法，即将整个非恒定流量过程概化为梯级恒定流，对每个梯级恒定流采用非恒定流方程，以时间步长为迭代参数进行迭代，直到得到恒定的流场。
3 数模计算有关问题 处理
4.1 初始条件及边界条件
初始条件：可由初始刻实测资料给出。计算所需的初始条件在实际计算中一般难以全部获得，只能通过估算加以补足，当然估算得越接近实际越好。
边界条件：进口给定流量、含沙量及其级配，出口给定水位；对于不滑动岸边界，取U=0,V= 0 ；对于滑动岸边界，取边界法线方向流速分量为0 。
对于活动边界，采用了动边界模拟技术：对每个膖时段，采用计算的水位及水深值判别和区分水域和陆域计算节点。对陆域计算节点，使其保持一较小富余水深(Hmin =0.001m) ，并取其糙率为一个接近于无穷大( 如1010 ) 的正数。
4.2 三角形网格划分
基于网格生成的基本思想，提出了四边形法。对于单一河道，首先根据河势将河道剖分为若干大的四边形，再分别将这些大的四边形剖分为若干小的四边形，最后将四边形的较短对角线相连，即形成三角形网格。对于分汊河道，将每条汊道当作单一河道处理后再进行拼接，便可形成整个计算域的网格。根据该方法，我们编制了通用的计算程序，只需输入少量的信息，计算机便能自动生成网格，并给出节点坐标和单元关联信息，最后配以屏幕显示及图形自动绘制。
4.3 紊动粘性系数
根据零方程紊流模型，紊动粘性系数由vt αU* H 公式确定。其中U* 为摩阻流速；α 为常数，经调试确定。对于泥沙扩散系数εs , 可近似取ε=vt 。
4.3 糙率
通过实测资料反求，并根据河道中不同的部位分块调试糙率。
4.4 床沙级配计算方法
在泥沙冲淤频繁的河段，由于水流与泥沙的相互作用，使某组泥沙发生冲刷时，另一组泥沙可能发生淤积。因此，床沙级配的计算应能反映泥沙冲淤交替过程。本模型将计算河段内河床泥沙概化成三层，即表层(交换层)、次表层(扰动层)和深层。当表层某组泥沙发生冲刷时，处在次表层的该组泥沙将受到扰动，根据扰动强度的大小，确定该层泥沙进入表层参加交换的量。
5 计算实例
河段概况：监利河弯段位于下荆江，上起塔市驿，下止于沙夹边，长约20km 。该河段平滩河宽约1400m ，最宽处约3200m ( 乌龟洲) ，最窄处约1000m ，平滩水位下的平均水深约11m 。河段内有一高程为30m ( 黄海) 左右的乌龟洲( 洲长约7km ，宽约2km ) 将河道分成左右两汊。本河段内有姚圻脑水文站，该站多年平均流量为11.370m3 ／s ，多年平均含沙量为1.121kg /m3 。河床主要由粉质粘土、砂粘土和细砂组成，河床表层床沙中值粒径一般在0.00 ～0.22mm 之间。
图1 计算河段河势及网格布置图 Fig.1 The calculated river reach and the grids
计算区域的选取及网格划分：选取塔市驿至沙夹边长约20km 的河段作为计算河段。采用三角形网格，将整个计算河段划分为100 个断面，6732 个三角形网格单元，共有3500 个节点，见图1 。
地形条件：以1993年10月实测河道地形为计算起始地形，并用1996年10月实测地形进行验证。
计算时段及水沙条件：选用1993年10月～1996年10月水沙资料，按流量共划分为120个计算时段。进口按姚圻脑水文站流量和分组含沙量给出，出口水位根据沙夹边与姚圻脑水位相关关系给出。
计算结果分析 ：
图2 为姚圻脑计算水位与实测水位对比图，由图中可见计算水位与实测水位吻合较好。
图3为不同流量下流场示意图，由图中可见大水漫滩，小水归槽，水流平顺，流场合理。
图4为1996年10月计算和实测的河道地形对比图，由图中可见通过计算得到的地形等高线与实测河道地形等高线的位置和范围基本吻合。
图2 计算与实测水位对比图 Fig.2 Comparison between calculated and
图3 不同流量下流场示意图 Fig.3 Flow fields for different discharges
图4 计算与实测河道地形比较图 Fig.4 Comparison between calculated and measured channel
从验证时段的冲淤总量来看，河段内实际淤积总量为2820万m3 ，计算淤积总量为2852万m3 ，相对误差仅为1.1%。由此可见，本模型能较好反映本河段的河床变形情况。
6 结论
(1) 对于有限单元法在河道平面二维水流泥沙计算中的应用 作了一些成功的探讨，建立了一套行之有效的解法。有限元法具有能够很好地模拟复杂边界和河道地形等优点，但若不经过处理，就需要大量的存储单元和很长的运算时间。为克服这些缺点，作者采用了质量集中的方法，并提出了压缩存储的方法，使得计算存储量大为减少。另外，还采用了“ 预报- 校正- 迭代” 的时间推进算法，并提出了“ 非恒定- 恒定- 非恒定流” 的算法，大大减少了计算量。
(2) 针对有限元前期单元剖分工作量巨大的问题，提出并运用四边形法开发了三角形网格自动生成系统。该系统不仅能适用于单一河道，而且能适用于分汊和支流出汇、入汇等各种复杂情况。
(3) 为检验模型的可靠性，文中以下荆江的监利河段为例进行了验证计算。计算结果表明该模型能较好模拟本河段的水流泥沙运动及河床变形情况。
(4) 由于天然河道的复杂性，目前 关于泥沙运动理论 还不够成熟，且验证资料尚不充分，因此数模计算精度有待作进一步的提高。
参考 文献 计算 水力学［M］.清华大学出版社，1995.9. 问题 的一种有限元解法［J］.武汉水利电力学院学报，1990年6月. 理论 ［J］.水利水运科学 研究 ，1983. 工业 出版社，1993.9 .
减小字体
增大字体
= 
