河道平面二维水沙数学模型的有限元方法（精简版）

河道平面二维水沙数学模型的有限元方法
河道平面二维水沙数学模型的有限元方法
摘要：采用有限元方法建立起一套河道平面二维水流泥沙数学模型。

在前人研究的基础上，采用了质量集中的处理方法，提出了压缩存储的方法，从而大大减少了计算存储量。

针对有限元法时间步长需取得较短问题，采用了“预报-校正-迭代”的算法，提出了“非恒定-
恒定-非恒定流”的算法，既能解决工程实际问题，又大大减少了计算量。

作者以下荆江监利河段为例进行泥沙冲淤计算，计算结果与实测值符合较好，从而证明了模型的可
摘要：采用有限元方法建立起一套河道平面二维水流泥沙数学模型。

在前人研究的基础上，采用了质量集中的处理方法，提出了压缩存储的方法，从而大大减少了计算存储量。

针对有限元法时间步长需取得较短问题，采用了“预报-校正-迭代”的算法，提出了“非恒定-
恒定-非恒定流”的算法，既能解决工程实际问题，又大大减少了计算量。

作者以下荆江监利河段为例进行泥沙冲淤计算，计算结果与实测值符合较好，从而证明了模型的可靠性。

关键词：水流泥沙有限元模型验证三峡工程建成后，水库将拦蓄大量泥沙，下泄水流含沙量减小，对三峡工程坝下游河道将产生以冲刷为主的影响，包括对荆江河段的河势及荆江大堤带来影响。

为研究坝下游重点河段的河床冲淤分布、河势变化、近岸流速变化等问题，一维模型显得无能为力，但可采用平面二维模型来解决。

有限元方法可采用无结构化网格，能很好地模拟不规则的几何形状，因此很适合于对天然河道的模拟。

然而，正如其它方法一样，有限元法也有它的缺点，主要是计算存储量和运算量较大。

为扬长避短，使有限元方法能运用到对天然河道的模拟上来，本模型运用质量集中［4］的方法将系数矩阵转化为三对角矩阵，并提出了紧凑的分块压缩存储方法，从而大大减少了计算存储量，使得计算能在一般微机上进行。

采用质量集中方法的不足之处是时间步长需取得较短，且在河道模拟中尤为突出（10）迭代给定误差洌杂谒械1≤i≤NP，若，则令：，否则令：= ；转到校正，
继续迭代，直到满足精度要求为止。

当求出Hi、Ui、Vi、Si和Ni后，
代入河床变形方程便可求得冲淤变形后的河床高程Zn 1bi，这样便
完成了一个时段的计算。

以上角标n表时段，*表预报值，**表校正值，无*标记的表示该时段的计算结果。

在水位求解过程中还采用了Kaahara“选择系数集中［4］”的经验处理方法，引入了选择性集中系数：α=0 8～1 0，其作用相当于对计算结果进行“光滑”处理。

由于系数矩阵为NP×NP阶矩阵，若直接进行存储会占用大量的内存，从而导致计算无法进行。

笔者发现系数矩阵中绝大部分为零元素，为此提出了压缩存储的方法，即只对系数矩阵中的非零元素进行存储和计算，同时用辅助数组指明相应非零元素所在的位置。

对于长系列水沙条件的计算，提出了“非恒定-定非-定流”的算法，即将整个非恒定流量过程概化为梯级恒定流，对每个梯级恒定流采用非恒定流方程，以时间步长为迭代参数进行迭代，直到得到恒定的流场。

3 数模计算有关问题处理4 1 初始条件及边界条件初始条件：可由初始刻
实测资料给出。

计算所需的初始条件在实际计算中一般难以全部获得，只能通过估算加以补足，当然估算得越接近实际越好。

边界条件：进口给定流量、含沙量及其级配，出口给定水位；对于不滑动岸边界，取U=0,V= 0；对于滑动岸边界，取边界法线方向流速分量为0。

对于活动边界，采用了动边界模拟技术：对每个膖时段，采用计算的水位及水深值判别和区分水域和陆域计算节点。

对陆域计算节点，使其保持一较小富余水深，并取其糙率为一个接近于无穷大的正数。

4 2 三角形网格划分基于网格生成的基本思想，提出了四边形法。

对于单一河道，首先根据河势将河道剖分为若干大的四边形，再分别将这些大的四边形剖分为若干小的四边形，最后将四边形的较短对角线相连，即形成三角形网格。

对于分汊河道，将每条汊道当作单一河道处理后再进行拼接，便可形成整个计算域的网格。

根据该方法，我们编制了通用的计算程序，只需输入少量的信息，计算机便能自动生成网格，并给出节点坐标和单元关联信息，最后配以屏幕显示及图形自动绘制。

4 3 紊动粘性系数根据零方程紊流模型，紊动粘性系数由vtαU*H公式确定。

其中U*为摩阻流速；α为常数，经调试确定。

对于泥沙扩散系数εs ,可近似取ε=vt。

4 3 糙率通过实测资料反求，并根据河道中不同的部位分块调试糙率。

4 4 床沙级配计算方法在泥沙冲淤
频繁的河段，由于水流与泥沙的相互作用，使某组泥沙发生冲刷时，另一组泥沙可能发生淤积。

因此，床沙级配的计算应能反映泥沙冲淤交替过程。

本模型将计算河段内河床泥沙概化成三层，即表层、次表层和深层。

当表层某组泥沙发生冲刷时，处在次表层的该组泥沙将受到扰动，根据扰动强度的大小，确定该层泥沙进入表层参加交换的量。

5 计算实例河段概况：监利河弯段位于下荆江，上起塔市驿，下止于沙夹边，长约20km。

该河段平滩河宽约1400m，最宽处约3200m，最窄处约1000m，平滩水位下的平均水深约11m。

河段内有一高程为30m 左右的乌龟洲将河道分成左右两汊。

本河段内有姚圻脑水文站，该站多年平均流量为11 370m3／s，多年平均含沙量为1 1kg/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月水沙资料，按流量共划分为0个计算时段。

进口按姚圻脑水文站流量和分组含沙量给出，出口水位根据沙夹边与姚圻脑水位相关关系给出。

计算结果分析：图2为姚圻脑计算水位与实测水位对比图，由图中可见计算水位与实测水位吻合较好。

图3为不同流量下流场示意图，由图中可见大水漫滩，小水归槽，水流平顺，流场合理。

图4为1996年10月计算和实测的河道地形对比图，由图中可见通过计算得到的地形等高线与实测河道地形等高线的位置和范围基本吻合。

图2 计算与实测水位对比图Fig 2 Comparison beteen calculated andmeasured ater level 图3 不同流量下流场示意图Fig 3 Flo fields for different discharges图4 计算与实测河道地形比较图Fig 4 Comparison beteen calculated and measured channel 从验证时段的冲淤总量来看，河段内实际淤积总量为2820万m3，计算淤积总量为2852万m3，相对误差仅为1 1%。

由此可见，本模型能较
好反映本河段的河床变形情况。

6 结论对于有限单元法在河道平面二维水流泥沙计算中的应用作了一些成功的探讨，建立了一套行之有效的解法。

有限元法具有能够很好地模拟复杂边界和河道地形等优点，但若不经过处理，就需要大量的存储单元和很长的运算时间。

为克服这些缺点，作者采用了质量集中的方法，并提出了压缩存储的方法，使得计算存储量大为减少。

另外，还采用了“预报-校正-迭代”的时间推进算法，并提出了“非恒定-恒定-非恒定流”的算法，大大减少了计算量。

针对有限元前期单元剖分工作量巨大的问题，提出并运用四边形法开发了三角形网格自动生成系统。

该系统不仅能适用于单一河道，而且能适用于分汊和支流出汇、入汇等各种复杂情况。

为检验模型的可靠性，文中以下荆江的监利河段为例进行了验证计算。

计算结果表明该模型能较好模拟本河段的水流泥沙运动及河床变形
情况。

由于天然河道的复杂性，目前关于泥沙运动理论还不够成熟，且验证资料尚不充分，因此数模计算精度有待作进一步的提高。

参考文献［1］周雪漪计算水力学［M］清华大学出版社，1995 9 ［2］谢鉴衡主编河流模拟［M］水利电力出版社，1988 6 ［3］张瑞瑾，谢鉴衡等河流泥沙动力学［M］水利电力出版社，1989 6 ［4］Kaahara, M , Hirano, H , Tsubota, K , and Inayuki, K Selective Lumping Finite Element Method for Shallo Water Flo, Internation Journal for Nu merical Method in Fluids, 1980,vol 2:89-1 ［5］韦直林河道水流泥沙问题的一种有限元解法［J］武汉水利电力学院学报，1990年6月［6］窦国仁河口海岸全沙模型相似理论［J］水利水运科学研究，1983 ［7］李义天等三峡水库变动回水区泥沙数学模型研究及初步运用长江三峡工程泥沙与航运关键技术研究专题研究报告集武汉工业出版社，1993 9 ［8］陈国祥三峡以下重点浅滩河段二维水流泥沙数值模拟长江三峡工程泥沙与航运关键技术研究
专题研究报告集武汉工业出版社，1993 9 ［9］长江水利委员会三峡工程泥沙研究湖北科学技术出版社，1997 10。