河道平面二维水沙数学模型的有限元方法(1)程(2)

3　数值解法

　　为书写方便，采用通用变量p来代替方程(13)～(17)中的变量zs、u、v、s和n，用fp表示方程中等号右边项，并用对角矩阵 来代替a_ij，则方程可化成统一的形式： dp_i/dt=fp_i。

对该方程的求解，模型采用了“预报-校正-迭代^［5］”的计算方法。其方法为：用二阶显式adams公式作为预测公式，梯形公式(隐式)作校正公式，可构造adams二阶pc公式。离散后的方程不需联接，计算过程稳定性好，时间步长可取得较长，同时还有效避免了数 值振荡。具体计算过程如下

预报　　　　　　　　　　　　　　　　(18)

校正　　　　　　　　　　　　　(19)

　迭代给定误差洌杂谒械1≤i≤np，若　　　，则令： _{，否则令：= ；转到校正，继续迭代，直到满足精度要求为止。当求出hi、ui、vi、si和ni后，代入河床变形方程便可求得冲淤变形后的河床高程zⁿ⁺¹bi，这样便完成了一个时段的计算。以上角标n表时段，*表预报值，**表校正值，无*标记的表示该时段的计算结果。}

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

由于系数矩阵为np×np阶矩阵，若直接进行存储会占用大量的内存，从而导致计算无法进行。笔者发现系数矩阵中绝大部分为零元素，为此提出了压缩存储的方法，即只对系数矩阵中的非零元素进行存储和计算，同时用辅助数组指明相应非零元素所在的位置。

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

3　数模计算有关问题处理

4.1　初始条件及边界条件

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

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

对于活动边界，采用了动边界模拟技术：对每个膖时段，采用计算的水位及水深值判别和区分水域和陆域计算节点。对陆域计算节点，使其保持一较小富余水深(h_min=0.001m)，并取其糙率为一个接近于无穷大(如10¹⁰)的正数。

4.2　三角形网格划分

基于网格生成的基本思想，提出了四边形法。对于单一河道，首先根据河势将河道剖分为若干大的四边形，再分别将这些大的四边形剖分为若干小的四边形，最后将四边形的较短对角线相连，即形成三角形网格。对于分汊河道，将每条汊道当作单一河道处理后再进行拼接，便可形成整个计算域的网格。根据该方法，我们编制了通用的计算程序，只需输入少量的信息，计算机便能自动生成网格，并给出节点坐标和单元关联信息，最后配以屏幕显示及图形自动绘制。

4.3　紊动粘性系数

根据零方程紊流模型，紊动粘性系数由v_tαu_*h公式确定。其中u_*为摩阻流速；α为常数，经调试确定。对于泥沙扩散系数ε_s ,可近似取ε=v_t。

4.3　糙率

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

4.4　床沙级配计算方法

在泥沙冲淤频繁的河段，由于水流与泥沙的相互作用，使某组泥沙发生冲刷时，另一组泥沙可能发生淤积。因此，床沙级配的计算应能反映泥沙冲淤交替过程。本模型将计算河段内河床泥沙概化成三层，即表层(交换层)、次表层(扰动层)和深层。当表层某组泥沙发生冲刷时，处在次表层的该组泥沙将受到扰动，根据扰动强度的大小，确定该层泥沙进入表层参加交换的量。

5　计算实例

河段概况：监利河弯段位于下荆江，上起塔市驿，下止于沙夹边，长约20km。该河段平滩河宽约1400m，最宽处约3200m(乌龟洲)，最窄处约1000m，平滩水位下的平均水深约11m。河段内有一高程为30m(黄海)左右的乌龟洲(洲长约7km，宽约2km)将河道分成左右两汊。本河段内有姚圻脑水文站，该站多年平均流量为11.370m³／s，多年平均含沙量为1.121kg/m³。河床主要由粉质粘土、砂粘土和细砂组成，河床表层床沙中值粒径一般在0.00～0.22mm之间。

图1 计算河段河势及网格布置图 （科教论文网 lw.NsEac.com编辑整理）
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
measured water level

图3　不同流量下流场示意图 （科教论文网 lw.nSeAc.com编辑发布）
fig.3　flow fields for different discharges

图4　计算与实测河道地形比较图
fig.4　comparison between calculated and measured channel 

 　　从验证时段的冲淤总量来看，河段内实际淤积总量为2820万m³，计算淤积总量为2852万m³，相对误差仅为1.1%。由此可见，本模型能较好反映本河段的河床变形情况。

6　结论

　　(1)对于有限单元法在河道平面二维水流泥沙计算中的应用作了一些成功的探讨，建立了一套行之有效的解法。有限元法具有能够很好地模拟复杂边界和河道地形等优点，但若不经过处理，就需要大量的存储单元和很长的运算时间。为克服这些缺点，作者采用了质量集中的方法，并提出了压缩存储的方法，使得计算存储量大为减少。另外，还采用了“预报-校正-迭代”的时间推进算法，并提出了“非恒定-恒定-非恒定流”的算法，大大减少了计算量。

　　(2)针对有限元前期单元剖分工作量巨大的问题，提出并运用四边形法开发了三角形网格自动生成系统。该系统不仅能适用于单一河道，而且能适用于分汊和支流出汇、入汇等各种复杂情况。

　　(3)为检验模型的可靠性，文中以下荆江的监利河段为例进行了验证计算。计算结果表明该模型能较好模拟本河段的水流泥沙运动及河床变形情况。

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

共3页: 2