悬沙运动的快速特征有限元模型柳海涛1 王小松1 夏庆福1 (1. 中国水利水电科学研究院，北京，100038，htliou@)摘要：本文在已有的研究基础上，发展了一种快捷的水沙特征有限元计算模型，来模拟以对流作用为主的悬沙运动及河床变形。
该模型采用一种预测-校正方法来对特征点进行搜索，提高了特征含沙量的计算精度，对粘土质河床的处理引入粘性土冲刷平衡流速作为判别指标。
另外对长系列计算下时间步长的取值、含沙量的数值耗散、大型矩阵的存储方法以及动边界处理也给出了具体方法。
通过典型算例表明，模型具有较好的适应性和可靠性。
关键词：特征有限元模型；特征点；对流为主的悬沙运动；冲淤平衡流速1 前言当前，二维泥沙数学模型在基本理论和计算方法上已相当成熟，国内外出现了许多模块化和商业化的计算软件，但这些软件在使用时一是没有源代码，用户无法根据实际需要进行修改；二是其采用的数值方法在一些特殊的情况下并不适用。
因此，为解决实际问题，在现有的模型的基础上进行二次开发是一种必要的手段，在文献[1]中，作者就提出了一种简洁的求解悬移质泥沙运动的特征有限元方法，该方法在求解精度和速度上基本可以胜任河道泥沙计算，其理论精度较以往的集中质量法以及许多常用的隐格式都要高。
然而，该方法在泥沙特征点的搜索、单元干湿判别、粘性土质的冲刷、大型稀疏矩阵的存贮和求解方法等方面的研究还不够深入。
本文在已有水流模型的基础上，发展了一种基于特征有限元方法的水沙计算模型，并对含沙量特征点的搜索方法、长系列水沙计算时间步长的取值、粘性土质和局部刚性河床的处理等方面进行研究，同时对模型中含沙量的数值耗散、单元的干湿变化的处理、大型矩阵的存储技术和计算方法等问题进行改进，使其更加实用。
最后通过典型算例来验证该模型的适应性和可靠性。
2 基本方程本文模型的水流控制方程如下0=∂∂+∂∂+∂∂yhv x hu t h (1) 2/12222222)()(v u Cgu y u E x u E h x z h gh y u hv x u hu t u h xy xx b +-⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+∂+∂-=∂∂+∂∂+∂∂ρ (2) 2/12222222)()(v u Cgv y v E x v E h y z h gh y v hv x v hu t v h xy xx b +-⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+∂+∂-=∂∂+∂∂+∂∂ρ (3) 其中u ，v 为流速（m/s ）；h 和b z 为水深和床面高程(m)；E 为综合紊动粘滞系数(N ·s/m 2)；C 为谢才系数，6/11h n C =。
泥沙输移控制方程0)()(*2222=-+∂∂+∂∂-∂∂+∂∂+∂∂s s h ys x s y s v x s u t s c s αωβε (4)河床变形方程)(*s s tz c -=∂∂αωβγ (5) 其中s ε为泥沙扩散系数(m 2/s)，γ为淤积物干容重(kg/m 3)，ω为颗粒沉速(m/s)，s 和*s 分别为垂线平均含沙量和水流挟沙力(kg/m 3)，α为恢复饱和系数，c β为冲刷判别系数，其作用将在后面介绍。
根据两相流理论，对于泥沙输移方程，因其忽略惯性力、升力、泥沙应力梯度等的作用，当泥沙运动以扩散作用为主时，仅适用于含细颗粒的低浓度挟沙水流，而对于包含粗颗粒的挟沙水流则应采用动理学模型来计算[2]；而当泥沙运动以对流作用为主，特别是当泥沙粒径足够小，颗粒驰豫时间远小于水流微元自身的拉格朗日积分时间尺度时，颗粒能充分跟随水流紊动，泥沙运动方程的精度就很高。
关于非均匀沙问题，除了考虑来沙条件外，对于床沙，文献[3]和[4]认为，应当将河床概化为表层、次表层和深层，但计算中将表层的厚度设定为2m ，这样整个河段冲淤基本都限于此层。
若进一步对该层采用分层记忆模式，就需要知道计算河段的床沙级配的垂线分布资料。
本文认为，对表层和次表层的划分要客观，采用分层记忆模式的厚度应依据河段的历史冲淤资料来确定，不应只包含表层。
其次，泥沙粒径级配对淤积的影响主要体现在较短距离内存在淤积相位差，如水库的淤积等，但河道内流速和水深变化复杂，泥沙淤积模式与此并不完全相同，河床各分层的泥沙粒径与其深度也没有直接的关系。
另外，泥沙浓度和分层固结对计算结果影响也不能忽略。
在粘土质河床中，当冲淤幅度较小时，作为简化，对表层泥沙可以直接采用一个特征粒径，同时该粒径值可以是沿程变化的。
限于作者的水平，上述问题还有待下一步研究，本文模型暂时采用如下假定：(1) 悬移质为较均质细沙，床沙质为粘性土或相同的均质细沙。
(2) 挟沙力采用一维水流挟沙力公式：m gh U k s )(3*ω=，其中k 和m 取常数。
(3) 泥沙扩散系数s ε取常数；恢复饱和系数α在冲刷时为1.0，淤积时为0.25。
3 模型的数值方法3.1 特征有限元方法对于对流占优问题，文献[4]认为，传统的Galerkin 法等价于中心差分格式，因缺乏足够的耗散，往往导致数值振荡，并为此采用了高分辨率格式对对流项进行重构，即通过引入几乎相等的扩散与反扩散以保证格式的高精度，同时利用限制因子保证影响系数的非负性和解的保单调性；而文献[3]中则采用了质量集中方法辅以预报-校正-迭代的算法，但其对于泥沙计算的效果没有作详细介绍。
本文认为，上述格式对于流速和水深的计算是有效的，但对于含沙量、传质的计算来说并不完全适用，因为该类方法只是以保证计算稳定收敛为主要目的，没有任何的物理意义，而且更重要的是，无法保证计算结果的非负性。
事实上，流场中任一点的含沙量的变化除了自身的紊动扩散外，主要与其上游的含沙量输运有关，例如回流区内的泥沙主要是靠回流本身从高流速区挟带来的，而不是周围节点泥沙的扩散作用。
因此，本文采用了文献[1]中的特征有限元方法，该方法只需沿着特征线找到特征点，就可以对输运方程进行简化。
具体方法如下：在悬移质泥沙运动方程中，时间导数项和对流项组成的全导数算子可以用下面的隐格式代替：tS S y S v x S u t S Dt DS K T T ∆-=∂∂+∂∂+∂∂=)( （6） 其中K S 为t T ∆-时刻的特征点K 处的含沙量值，K 点的坐标为⎪⎩⎪⎨⎧-=-=⎰⎰∆∆t T K t T K vdtY Y udt X X (7) 这样泥沙输运方程作如下离散 0)(*2222=⎥⎥⎦⎤⎢⎢⎣⎡-+⎪⎪⎭⎫ ⎝⎛∂∂+∂∂-⎰dA S S h S y x Dt DS e T T T ϕϕαωϕϕεϕϕϕ (8) 将式(6)代入上式，经整理得到如下有限元离散方程B S A I =∙ (9) 系数矩阵dA h dx y dy x dA y y x x t dAA T L e c T T e T T e T ϕϕαωβϕϕϕϕεϕϕϕϕεϕϕ⎰⎰⎰⎰+∂∂-∂∂-⎪⎪⎭⎫ ⎝⎛∂∂∂∂+∂∂∂∂+∆= (10) 右端向量⎰⎰∙+∙∆=e IT c K I eT S dA h S t dA B *ϕϕαωβϕϕ (11) 基函数],,[ζηεϕ=单元节点含沙量],,[321S S S S I =单元节点的特征点处的含沙量],,[321K K K K I S S S S =单元节点的水流挟沙力],,[*3*2*1*S S S S I =这样，对每一个单元e 求出单元系数矩阵A 和右端向量B ，然后将其合成总体矩阵，即可求得n T 时刻的含沙量S 。
对于河床高程的修正，仍采用式(5)中的河床变形方程，在求解Z 之前，先应用如下的隐式差分格式： 11)1(--∂∂-+∆-=∂∂n n n n t Z K t Z Z Kt Z(12)其中，K 值取1.0~2.0。
将上式代入式(5)中，得到求出n T 时刻的床面高程表达式:⎥⎦⎤⎢⎣⎡⎪⎭⎫ ⎝⎛∂∂∆---∆+=--1*11)]([n n C n n t Z t K K S S K t Z Z γαωβ (13)因此，在模型中，先用式(13)对河床进行修正，再使用式(12)求得该时刻的地形变化率，以便下一循环使用。
3.2 特征点的确定关于特征点的搜索方法，文献[1]对其落在计算域内部和计算节点位于第二类边界上（侧向边界）的情况做了较详细的研究，但是对于特征点位于计算域之外的情况，还未提出相应的处理方法。
同时，从式(7)可以看出，由于包含时间积分，在流线变化剧烈的地区，特征点的确定有一定困难，在文献[1]偏差，如图1中O /点流速矢0U 和半时间步长t ∆5.0点的流速矢U '，然后用矢量旋转的方法，矢1U ，1U 满足01U U =，且U U '||1全时间步长t ∆求出特征点K 二类边界之外的情况大大减少。
梯度，因此该方法不失为一种解决办法。
3.3 长系列冲刷问题的处理对于长系列泥沙计算问题，一般取较大的时间步长。
然而，当计算河段以冲刷为主，流速很大时，在靠近边界的特征线长度范围内，节点对应的特征点均位于边界之外，若时间步长过大，则用边界值代替特征值会有较大的误差。
因此，这里就产生一个矛盾，在采用隐格式后,时间步长可以取得较大，但为了保证泥沙计算精度,又要求较小的步长。
对此本文的处理方法是，对于水流和泥沙的计算分别采用不同的时间步长，在一个较长的水流计算步长T 内，进行N 次泥沙及床面冲淤计算，每次泥沙计算的步长为N T /，期间水流条件保持不变。
3.4 粘土质河床或刚性河床的处理方法对于粘土河床，参照铁道部第四勘测设计院的实测经验公式可知，当冲刷达到平衡时，流速和水深满足下面的关系式：5/3)1(33.0p lS h I V = (14) 式中S V 是冲刷平衡时的流速，p h 冲刷后的最大水深,l I 为粘土的液性指数。
本文模型以式(14)作为判断河道是否开始冲刷的标准，为此在泥沙运动方程和河床变形方程中，源（汇）项的系数除了恢复饱和系数α外，还有一个冲刷判断系数c β。
当垂线平均流速S V V ≥时，河道开始冲刷，冲刷的速度和强度可以参照一维挟沙力公式，此时有0.1=α，0.1=c β，因此含沙量允许暂时出现超饱和状态；当垂线平均流速s V V <时，有可能开始淤积，但如果水流含沙量处于不饱和状态，河道也不会发生淤积，此时有25.0=α，0.0=c β。
对于诸如衬砌河段或人工建筑物，则只允许淤积而不会冲刷，因此当垂线平均流速S V V ≥（粘土河床）或含沙量处于不饱和状态（散粒体河床）时，仍有0.0=α，0.0=c β，当发生淤积时，有25.0=α0.1=c β。
3.5 含沙量的数值耗散的处理从离散式(10)中可以看到，左端第三项为沿单元边界L 上的泥沙扩散通量，理论上讲，含沙量的空间导数处处连续，因此在总体系数矩阵中，对内部节点，该项的贡献为0，只对于边界节点，该项不为0，在许多泥沙数值模型中，为了节省计算量，将其略去，即假定边界上没有泥沙通过浓度扩散进入水体。
