河道平面二维水沙数学模型的有限元方法

1、河道平面二维水沙数学模型的有限元方法    然而，正如其它方法一样，有限元法也有它的缺点，主要是计算存储量和运算量较大。为扬长避短，使有限元方法能运用到对天然河道的模拟上来，本模型运用质量集中 4 的方法将系数矩阵转化为三对角矩阵，并提出了紧凑的分块压缩存储方法，从而大大减少了计算存储量，使得计算能在一般微机上进行。采用质量集中方法的不足之处是时间步长需取得较短，且在河道模拟中尤为突出 ( 因河道比较窄长，网格需划分很细，而该法的稳定性要 求时间步长与网格尺度成正比 ) 。针对该问题，笔者采用了 “ 预报 - 校正 - 迭代 5 ” 的算法，该法可加大时

2、间步长，同时有效避免了数值震荡。针对长系列水沙条件下计算量较大问题，作者又提出了 “ 非恒定 - 恒定 - 非恒定流 ” 的算法，该算法既能解决工程实际问 题，又大大减少了计算量，使有限元方法能够很好地运用于河道水流泥沙问题的实际计算。 1 基本方程 平面二维水流方程 (1) (2) (3) 悬移质泥沙扩散方程 (4) 推移质不平衡输移方程6 (5) 河床变形方程 由悬移质引起的河床变形方程为 <!endif> (6) 由推移质引起的河床变形方程为 (7) 以上各式中 U ， V 分别为垂线平均流速在 x,y 方向上的分量； Zs 、 Zb 和 H 分别为水位、河底高程和水深； g

3、 为重力加速度； vt 为水流紊动粘性系数； 为水的密度； x 、 y 、舄瓂分别为底部切应力在 x 和 y 方 (x 、 y)= ，向上的分量： C 为谢才系数，常用曼宁公式计算： C=H1/6/n ； S 和 S* 分别为垂线平均 含沙量和挟沙力； N 和 N* 分别为推移质输沙量和推移质输沙能力折算成全水深的泥沙浓度； s 为泥沙紊动扩散系数； 为泥沙沉速； 为床沙干容重； 为悬移质泥沙恢复饱和系数，淤积时取 0.25 ，冲刷时取 1.0 ； 为推移质泥沙恢复饱和系数，取 0.25 。      悬移质泥沙共分成8组，水流挟沙力和分组挟沙

4、力级配采用李义天方法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) 由此，经整理得到如下有限元方

5、程 (13) (14) (15) (16) (17) 其中 i,j,k=1,2,3, ， NP. 对该方程的求解，模型采用了 “ 预报 - 校正 - 迭代 5 ” 的 计算 方法 。其方法为：用二阶显式 Adams 公式作为预测公式，梯形公式 ( 隐式 ) 作校正公式，可构造 Adams 二阶 PC 公式。离散后的方程不需联接，计算过程稳定性好，时间步长可取得较长，同时还有效避免了数 值振荡。具体计算过程如下 预报 (18) 校正 (19) 迭代给定误差洌杂谒械 1iNP ，若 ，则令：      ，否则令： = ；转到校正，继续迭代，直到满足

6、精度要求为止。当求出 Hi 、 Ui 、 Vi 、 Si 和 Ni 后，代入 河床变形方程便可求得冲淤变形后的河床高程 Zn+1bi ，这样便完成了一个时段的计算。以上角标 n 表时段， * 表预报值， * 表校正值，无 * 标记的表示该时段的计算结果。     在水位求解过程中还采用了 Kawahara“ 选择系数集中 4 ” 的 经验 处理方法，引入了选择性集中系数： =0.8 1.0 ，    其作用相当于对计算结果进行 “ 光滑 ” 处理。    由于系数矩阵为 NP

7、×NP 阶矩阵，若直接进行存储会占用大量的内存，从而导致计算无法进行。笔者发现系数矩阵中绝大部分为零元素，为此提出了压缩存储的方法，即只对系数矩阵中的非零元素进行存储和计算，同时用辅助数组指明相应非零元素所在的位置。 对于长系列水沙条件的计算，提出了 “ 非恒定 - 定非 - 定流 ” 的算法，即将整个非恒定流量过程概化为梯级恒定流，对每个梯级恒定流采用非恒定流方程，以时间步长为迭代参数进行迭代，直到得到恒定的流场。 3 数模计算有关 问题 处理 4.1 初始条件及边界条件 初始条件：可由初始刻实测资料给出。计算所需的初始条件在实际计算中一般难以全部获得，只能通过估算加以补足，当然估

8、算得越接近实际越好。 边界条件：进口给定流量、含沙量及其级配，出口给定水位；对于不滑动岸边界，取 U=0,V= 0 ；对于滑动岸边界，取边界法线方向流速分量为 0 。 对于活动边界，采用了动边界模拟技术：对每个膖时段，采用计算的水位及水深值判别和区分水域和陆域计算节点。对陆域计算节点，使其保持一较小富余水深 (Hmin=0.001m) ，并取其糙率为一个接近于无穷大 ( 如 1010) 的正数。 4.2 三角形网格划分 基于网格生成的基本思想，提出了四边形法。对于单一河道，首先根据河势将河道剖分为若干大的四边形，再分别将这些大的四边形剖分为若干小的四边形，最后将四边形的较短对角线相连，即形成三

9、角形网格。对于分汊河道，将每条汊道当作单一河道处理后再进行拼接，便可形成整个计算域的网格。根据该方法，我们编制了通用的计算程序，只需输入少量的信息， 计算机 便能自动生成网格，并给出节点坐标和单元关联信息，最后配以屏幕显示及图形自动绘制。 4.3 紊动粘性系数 根据零方程紊流模型，紊动粘性系数由 vtU*H 公式确定。其中 U* 为摩阻流速； 为常数，经调试确定。对于泥沙扩散系数 s , 可近似取 =vt 。 4.3 糙率 通过实测资料反求，并根据河道中不同的部位分块调试糙率。 4.4 床沙级配计算方法 在泥沙冲淤频繁的河段，由于水流与泥沙的相互作用，使某组泥沙发生冲刷时，另一组泥沙可能发生淤

10、积。因此，床沙级配的计算应能反映泥沙冲淤交替过程。本模型将计算河段内河床泥沙概化成三层，即表层(交换层)、次表层(扰动层)和深层。当表层某组泥沙发生冲刷时，处在次表层的该组泥沙将受到扰动，根据扰动强度的大小，确定该层泥沙进入表层参加交换的量。 5 计算实例 河段概况：监利河弯段位于下荆江，上起塔市驿，下止于沙夹边，长约 20km 。该河段平滩河宽约 1400m ，最宽处约 3200m ( 乌龟洲 ) ，最窄处约 1000m ，平滩水位下的平均水深约 11m 。河段内有一高程为 30m ( 黄海 ) 左右的乌龟洲 ( 洲长约 7km ，宽约 2km ) 将河道分成左右两汊。本河段内有姚圻脑水文站

11、，该站多年平均流量为 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月实测地形

12、进行验证。 计算时段及水沙条件：选用1993年10月1996年10月水沙资料，按流量共划分为120个计算时段。进口按姚圻脑水文站流量和分组含沙量给出，出口水位根据沙夹边与姚圻脑水位相关关系给出。 计算结果 分析 ： 图 2 为姚圻脑计算水位与实测水位对比图，由图中可见计算水位与实测水位吻合较好。 图3为不同流量下流场示意图，由图中可见大水漫滩，小水归槽，水流平顺，流场合理。 图4为1996年10月计算和实测的河道地形对比图，由图中可见通过计算得到的地形等高线与实测河道地形等高线的位置和范围基本吻合。 图2计算与实测水位对比图 Fig.2 Comparison between calculate

13、d andmeasured water level     图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) 针对有限元前期单元剖分