悬移质泥沙有限元模型初探（柳海涛，吉祖稳，胡春宏）
(中国水利水电科学研究院 泥沙研究所)摘要：本文针对传统有限元方法在求解以对流作用为主的泥沙运动方程时所遇到的困难，分别采用迎风有限元和特征有限元法进行耦合求解。通过数值实验证实两种方法可以很好地处理对流项，其中尤以特征有限元法为最优。尽管目前该方法计算耗时稍多，但随着硬件技术的发展，将得到广泛的应用。关键词：对流作用为主的泥沙运动；迎风有限元法；特征有限元法；泥沙运动方程
基金项目：国家重点基础研究发展规划项目(G19990436)。
作者简介：柳海涛(1971-)，男，中国水科院工程师。
1 前言
目前，国内外泥沙数模的研制已有了相当大的进展，二维数模的理论研究已经成熟，以前鉴于计算速度较慢和理论上的复杂性，有限元法在水流和泥沙耦合计算中的应用在国内一直处于次要的地位，因此，发展适用于强对流、高含沙水流的高精度有限元泥沙模型，具有一定的现实意义。我们知道，对于大多数水流模型，在缓流情况下均能得到非常精确的解，当水流出现局部急流的情况下，只需在该处稍加一些粘滞性，即可将其抑制，并且它对整体的求解精度并无太大的影响[1]；而在以扩散为主的传质模型中，我们也可以通过增大扩散系数和滞后因子来避免出现浓度分布的不合理现象。但是，泥沙在河流中的运动机理非常复杂，一般情况下，无论是急流还是缓流，悬移质的对流作用始终处处占优，采用传统的有限元方法或者再通过增大泥沙扩散系数来保证计算结果合理是行不通的。同时，由于它还同床面发生交换从而改变水流条件，因而泥沙输运问题在理论上讲始终是一个非恒定问题。泥沙模型的计算除了借鉴水质模型中的数学方法外，还要对床面进行不断的修正来重新计算流场，通常的外部或人工的耦合求解方法只能适用于流速和含沙量的变化率均较小的情况。为此，本文主要研究迎风及特征有限元法在泥沙输移计算中的应用，编制相应的程序并进行算例的对比，从而得出合理的结论。
2 基本方程
水流控制方程如下
| (1) |
| (2) |
| (3) |
其中u，v为流速；h，zb为水深和床面高程；E为紊动粘滞系数；C为谢才系数，表达式为：C=1/nh1/6。
泥沙输移控制方程如下
| (4) |
河床变形方程
| (5) |
其中，εs为泥沙扩散系数(m2/s),γ′为淤积物干容重(kg/m3)，ω为颗粒沉速(m/s)，s和s*分别为垂线平均含沙量和水流挟沙力(kg/m3)，α为恢复饱和系数。为简化计算，本文采用如下假设：
(1)悬移质为均质沙并且床面有足够数量的床沙质可与其交换。
(2)泥沙扩散系数取一常数，挟沙力采用一维水流挟沙力公式：s*=k（U3/ghw)m。
3 本文数模的思想方法及实施
3.1 传统及迎风有限元数值方法
对于上述悬移质扩散方程，本文原先采用传统的有限元方法来编制计算程式，通过初步的数值实验证实，该方法对顺直无回流的区段计算结果较为合理，但在回流区以及含沙量和流速梯度均较大的区域（如回流区与主流区的过渡区域）均出现数值耗散，以至在回流区内含沙量出现负值。其中原因主要是，在突扩流场中回流区的泥沙主要靠回流来搬运，而不是依靠扩散作用，因此不能单纯依靠增大扩散系数来‘抹平’计算结果。传统的有限元法认为，任一点的含沙量为周围各点对其贡献的总和，夸大了下游节点对它的影响，其值就有可能变为负值。在已有的各类离散格式中，只有迎风有限元法和特征有限元法能够解决此类问题，它们均认为流场中任一点的含沙量只与其上游点的含沙量有关。据此我们拟对传统的有限元方法进行改进。
本文首先采用一种特殊的迎风格式[2]对传统的有限元方法作了一些改进，其中，对于时间导数项我们仍采用隐式差分方法，而对于对流项我们采用迎风有限元离散格式，对于扩散项和源（汇）项仍然使用传统的有限元方法。其基本原理是：对于空间任一节点，它的影响区域只能是其上游区域（如图1中阴影部分），在进行伽辽金积分时，对流项的积分区域只能取A1＋A2＋A3的区域面积。其中，上游垂线的交点M和N的值，采用插值方法来确定，即有下式
| SM=aSA+(1-a)SFSN=bSct(1-b)SD | (6) |
| a=LMF/LAFb=LND/LCD | (7) |
这样在对每个单元进行伽辽金积分时，应当先判断各个节点的流速方向及其影响域，如图1中(c)、(d)、(e)。在求解单元系数矩阵时，每行系数的积分域不同，如图1(c)中节点1的影响域不包含该单元，因此在系数矩阵中第1行的系数均为0；在图1(d)中节点2的影响域完全包含该单元，因此在系数矩阵中第2行系数不作改变；而在图1(e)中，节点3的影响域部分包含该单元，因此第3行系数应作调整，先插值求出M点的坐标，代替第2点的坐标，然后再进行积分。系数矩阵可以用(8)式表示
f′、e′、c′表达式与f、e、c相同，但积分域不同，α=L2M/L12 |
图1 单元系数矩阵元素的差别示意图 |
这样系数矩阵就成为不对称矩阵，并且由于引入系数a，使其具有“迎风”的性质。
由于对流项和其它各项的积分区域不同，所以在求出每个单元的系数矩阵后还需要存贮每个节点在该单元内的影响域面积，然后在总体合成时将其与其它单元内该节点相应的影响域进行累加，形成如图1(a)中的0点的总体有效影响域。最后将总体方程中各节点的系数除以相应的有效积分域面积，其中，对于对流项系数应除以其有效影响域面积，而对其它各项系数则除以总体积分域面积，即该节点相邻单元的总面积，这样方程的等式才能成立。采用上述迎风方法后，程序变得较为复杂，在计算中需要累加并存贮各节点的有效和总体两种积分域面积，使得程序运行速度比后来的特征有限元法要慢。而且我们也注意到系数矩阵中仍然引用了影响域之外的节点，如图1(e)中的节点2。因而在随后的突扩流场悬沙输移计算中，局部流速变化异常的地区仍会出现个别的含沙量浓度负值。然而该方法也具有它的优点，一是计算中时间步长可以取得较长，从而适合于求解长系列泥沙冲淤问题。二是作为一种选择，我们在以后的改进中可以通过修正迎风系数α来消除其数值耗散，保持计算的稳定性。
3.2 特征有限元方法[3]
前面提到过，在能够正确反映水流中泥沙输移的离散方法中，特征线方法是最精确的，只要我们能找到特征线，就可以在特征线上运用全导数算子对输运方程进行简化，这样在悬移质泥沙运动方程中，可以直接用下式的差分格式来代替时间导数项和对流项
DS/Dt|0=S0-SK/Δt | (9) |
其中SK为K点处的含沙量，K点的坐标为
| (10) |
若时间步长取得很短，则上式可以写作
| (11) |
这当然是一种简化表达，但其精度仍然要比迎风方法高。
这样，我们可以对泥沙输运方程作如下离散
| (12) |
上式中，基函数 φ=[ε,η,ξ]
单元节点的含沙量 SI=[S1S2S3]
单元节点的特征点处的含沙量 SKI=[SK1SK2SK3]
单元节点的水流挟沙力S*I=[S*1S*2S*3]
L为单元e落在求解域边界上的一边，对于内部单元，该项线积分为0。
这样，对于每一个单元e，我们均可得到其系数矩阵A和右端向量B，然后将其合成总体矩阵，即可求得tn时刻的含沙量S。在求得含沙量后，需要对河床高程进行修正，程序采用的仍然是河床变形方程
| (13) |
在求解Z时，先应用下列的导数差分格式[1]
| (14) |
将(14)式代入(13)式而得到tn时刻床面高程Zn的表达式
Zn=Zn-1+Δt/K[αω/γ(Sn-S*)]-1-K/KΔt· |  | (15) |
| 我们利用公式(15)对河床地形进行修正，然后再使用(14)式得出该时刻的地形变化率 |
| ,以便在下一时段修正地形时使用。 |
从上述理论可以看出，特征有限元模型中，最关键的是特征点K的确定及其含沙量SK的求解。如果我们知道tn时刻的流速场，则可以利用式(11)求出时刻的特征点K的坐标，然后再插值求出SK值。但在具体的程序实现中，就有几个问题需要解决：
(1)在求得xk、yk后，必须设法找到K点所在的单元，再插值求解该点的Sk值。本文的程序中，在每次时间步长完成后，将存贮K点所在的单元号，在下次时间步长开始时将首先从预先存贮的单元编号开始查找，若K点不在此单元内，则继续向该单元编号的上、下游进行双向查找，这样可以大大节省程序的运行时间，一旦找到该单元，我们可以根据坐标值xk、yk，运用空间线性插值求出Sk值。判别K点是否落在某单元内，只需求出其在该单元内的局部坐标值，若三者之和为1，则K点在该单元内。
(2)在求解边界节点特征点K的坐标时，由于边界形态复杂，常常会发现K点位于求解域之外，在本文的程序中边界节点做如下处理：我们假定边界上水流始终是沿边界流动的，因而边界节点的特征点也一定在边界上，具体说明见图2。对于D点，其特征点的坐标若仍按式(11)求解，则特征点的位置将会是M点，两点的距离为Ldm=V·Δt,但M点位于求解域以外，这显然是不对的。正确的位置应当是N点，它位于D点水流上游的边界上，并且它到D点的距离亦为V·Δt。若特征点的位置位于更远的N′处，则还需要沿E点的上游的边界继续查找。一旦找到点所在的单元边界即可利用边界节点的线性插值求解Sk值。因此程序中除了需要判别边界节点外，还必须记录各个边界节点的相互关系，以供在边界上搜索特征点时使用。
图2 边界上特征点的求解示意图 |
图3 上游流量曲线 |
4 模型合理性检验及结论
| 本文中针对上述改进方法，在已有的水流程式的基础上，均编制了水沙耦合计算程式，并对一个典型的单边突扩流场中悬移质泥沙输移进行模拟。计算部分网格如图4所示。具体计算条件如下：时间步长取90s，计算时段为2h，上游来流Q如图3所示，为一非恒定流，下游为恒定水位4m，糙率取0.025，挟沙力系数k=0.025，指数m=0.75[4]。初始时刻床面高程均为0.0m，一般地，随着上游来水量的变化，来沙量也是变化的，本文为简单起见，令上游进口的含沙量始终为5kg/m3。这里简单介绍一下特征有限元模型的计算结果，如图5所示。 |
图4 数值计算区域网格图 |
 |
| 图5 计算区域的地形及含沙量变化 Fig.5 Variation process of concentration and topography in the computational domain |
限于篇幅，这里仅截取上述三个时刻的计算结果。从床面高程的变化来看，计算区域经历了一个先冲后淤的过程，在时段初期，由于没有上游来沙补给，床面基本以冲刷为主，随着上游来沙的进入，计算区域自上而下逐渐转为淤积。通过计算结果可知，特征有限元法完全消除了含沙量浓度的数值损耗，特别是在主流和回流交界的地区，得出的含沙量分布更加符合理论分析结论。同时，我们也看到，有限元模型的求解精度同网格的剖分有很大的关系，如含沙量梯度较大的区域(参见图5)适当加密网格，效果会更好。
5 小结
(1)本文从传统的有限元方法出发，探讨了有限元数值解法在泥沙计算中的可行性，实验结果表明特征线法和带有迎风格式的有限元法是完全可以胜任河道泥沙计算的。鉴于计算机硬件技术的发展，该类方法将会有很大发展前景。
(2)在两种计算方法中，本文认为，迎风方法是一种折中方案，不但计算量大，再加上对边界节点上迎风格式的处理以及迎风指数的采用，使得程序变的庞杂，就目前来讲其计算效益仍较低，而特征线法相对较为简单快捷，计算效益较高，因此具有一定实际意义。
(3)尽管本文模型取得了一些进展，但由于在理论上采用了一些假设，因此程式还不能直接应用于生产实际，需要进一步完善。
参 考 文 献
[1] 柳海涛。河口水流数值模拟。河海大学硕士论文，1995.6.[2] 吴江航。二维非定常对流扩散与Navier-Stokes方程在不规则网格上的一种差分解。第一次偏微分方程数值解学术会议论文，1987.3.
[3] 金忠青等。水力学中的反问题。河海大学出版社，1995.7.
[4]韦直林。河道水流泥沙问题的一种有限元解法。武汉水利电力学院学报，1990.12.
