水利工程论文-河道平面二维水沙数学模型的有限元方法.doc

水利工程论文-河道平面二维水沙数学模型的有限元方法摘要：采用有限元方法建立起一套河道平面二维水流泥沙数学模型。在前人研究的基础上，采用了质量集中的处理方法，提出了压缩存储的方法，从而大大减少了计算存储量。针对有限元法时间步长需取得较短问题，采用了“预报-校正-迭代”的算法，提出了“非恒定-恒定-非恒定流”的算法，既能解决工程实际问题，又大大减少了计算量。作者以下荆江监利河段为例进行泥沙冲淤计算，计算结果与实测值符合较好，从而证明了模型的可靠性。关键词：水流泥沙有限元模型验证三峡工程建成后，水库将拦蓄大量泥沙，下泄水流含沙量减小，对三峡工程坝下游河道将产生以冲刷为主的影响，包括对荆江河段的河势及荆江大堤带来影响。为研究坝下游重点河段的河床冲淤分布、河势变化、近岸流速变化等问题，一维模型显得无能为力，但可采用平面二维模型来解决。有限元方法可采用无结构化网格，能很好地模拟不规则的几何形状，因此很适合于对天然河道的模拟。然而，正如其它方法一样，有限元法也有它的缺点，主要是计算存储量和运算量较大。为扬长避短，使有限元方法能运用到对天然河道的模拟上来，本模型运用质量集中4的方法将系数矩阵转化为三对角矩阵，并提出了紧凑的分块压缩存储方法，从而大大减少了计算存储量，使得计算能在一般微机上进行。采用质量集中方法的不足之处是时间步长需取得较短，且在河道模拟中尤为突出(因河道比较窄长，网格需划分很细，而该法的稳定性要求时间步长与网格尺度成正比)。针对该问题，笔者采用了“预报-校正-迭代5”的算法，该法可加大时间步长，同时有效避免了数值震荡。针对长系列水沙条件下计算量较大问题，作者又提出了“非恒定-恒定-非恒定流”的算法，该算法既能解决工程实际问题，又大大减少了计算量，使有限元方法能够很好地运用于河道水流泥沙问题的实际计算。1基本方程平面二维水流方程(1)(2)(3)(4)推移质不平衡输移方程6(5)河床变形方程由悬移质引起的河床变形方程为(6)由推移质引起的河床变形方程为(7)以上各式中U，V分别为垂线平均流速在x,y方向上的分量；Zs、Zb和H分别为水位、河底高程和水深；g为重力加速度；vt为水流紊动粘性系数；为水的密度；x、y、舄瓂分别为底部切应力在x和y方(x、y)=，向上的分量：C为谢才系数，常用曼宁公式计算：C=H1/6/n；S和S*分别为垂线平均含沙量和挟沙力；N和N*分别为推移质输沙量和推移质输沙能力折算成全水深的泥沙浓度；s为泥沙紊动扩散系数；为泥沙沉速；为床沙干容重；为悬移质泥沙恢复饱和系数，淤积时取0.25，冲刷时取1.0；为推移质泥沙恢复饱和系数，取0.25。悬移质泥沙共分成8组，水流挟沙力和分组挟沙力级配采用李义天方法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)由此，经整理得到如下有限元方程(13)(14)(15)(16)(17)其中i,j,k=1,2,3,，NP.3数值解法为书写方便，采用通用变量P来代替方程(13)(17)中的变量Zs、U、V、S和N，用FP表示方程中等号右边项，并用对角矩阵来代替Aij，则方程可化成统一的形式：dPi/dt=FPi。对该方程的求解，模型采用了“预报-校正-迭代5”的计算方法。其方法为：用二阶显式Adams公式作为预测公式，梯形公式(隐式)作校正公式，可构造Adams二阶PC公式。离散后的方程不需联接，计算过程稳定性好，时间步长可取得较长，同时还有效避免了数值振荡。具体计算过程如下预报(18)校正(19)迭代给定误差洌杂谒械1iNP，若，则令：，否则令：=；转到校正，继续迭代，直到满足精度要求为止。当求出Hi、Ui、Vi、Si和Ni后，代入河床变形方程便可求得冲淤变形后的河床高程Zn+1bi，这样便完成了一个时段的计算。以上角标n表时段，*表预报值，*表校正值，无*标记的表示该时段的计算结果。在水位求解过程中还采用了Kawahara“选择系数集中4”的经验处理方法，引入了选择性集中系数：=0.81.0，其作用相当于对计算结果进行“光滑”处理。由于系数矩阵为NPNP阶矩阵，若直接进行存储会占用大量的内存，从而导致计算无法进行。笔者发现系数矩阵中绝大部分为零元素，为此提出了压缩存储的方法，即只对系数矩阵中的非零元素进行存储和计算，同时用辅助数组指明相应非零元素所在的位置。对于长系列水沙条件的计算，提出了“非恒定-定非-定流”的算法，即将整个非恒定流量过程概化为梯级恒定流，对每个梯级恒定流采用非恒定流方程，以时间步长为迭代参数进行迭代，直到得到恒定的流场。3数模计算有关问题处理4.1初始条件及边界条件初始条件：可由初始刻实测资料给出。计算所需的初始条件在实际计算中一般难以全部获得，只能通过估算加以补足，当然估算得越接近实际越好。边界条件：进口给定流量、含沙量及其级配，出口给定水位；对于不滑动岸边界，取U=0,V=0；对于滑动岸边界，取边界法线方向流速分量为0。对于活动边界，采用了动边界模拟技术：对每个膖时段，采用计算的水位及水深值判别和区分水域和陆域计算节点。对陆域计算节点，使其保持一较小富余水深(Hmin=0.001m)，并取其糙率为一个接近于无穷大(如1010)的正数。4.2三角形网格划分基于网格生成的基本思想，提出了四边形法。对于单一河道，首先根据河势将河道剖分为若干大的四边形，再分别将这些大的四边形剖分为若干小的四边形，最后将四边形的较短对角线相连，即形成三角形网格。对于分汊河道，将每条汊道当作单一河道处理后再进行拼接，便可形成整个计算域的网格。根据该方法，我们编制了通用的计算程序，只需输入少量的信息，计算机便能自动生成网格，并给出节点坐标和单元关联信息，最后配以屏幕显示及图形自动绘制。4.3紊动粘性系数根据零方程紊流模型，紊动粘性系数由vtU*H公式确定。其中U*为摩阻流速；为常数，经调试确定。对于泥沙扩散系数s,可近似取=vt。4.3糙率通过实测资料反求，并根据河道中不同的部位分块调试糙率。4.4床沙级配计算方法在泥沙冲淤频繁的河段，由于水流与泥沙的相互作用，使某组泥