黄河下游平面二维水沙运动模拟的有限元方法
摘要：本文以水流和泥沙运动的研究成果为基础，建立了平面二维水沙模型。用有限元方法导出了本模型的离散方程式，用质量集中和预估校正法处理、迭代求解方程组。以黄河济南河段1976年汛期的洪水演进和河床演变为例从水位过程、流量过程、断面平均流速和最大流速、全域流速场和河床断面冲淤形态等方面，对模型进行了验证模拟。数值模拟计算结果与实测或物理模型试验结果符合较好。从而证明了本模型可靠性。
关键词：黄河下游 有限元方法 验证
1 前言
对于河道上修建桥涵等跨河工程、引程等建筑物后，人们所关心的河势变化、流速分布、河床局部冲淤形态和壅水等问题，一维模型是无能为力的，但可以用二维模型解决。在进行黄河水沙运动的模拟过程中，一些模型往往需要对方程组中的许多参数进行经验处理，而且在处理复杂的河道边界以及方程的离散和求解等方面还有许多问题值得研究。本模型以黄河科研究院在泥沙运动规律方面的研究成果为基础，建立模型框架。并选用适用于浑水的群体沉速公式计算泥沙沉速，引入适应于从清水到高含沙的水流挟沙能力计算公式和动床阻力计算公式，克服数学模型参数过多而不能通用的缺陷。
对模型研究区域的离散和方程的求解是数值模拟的关键。在离散方法方面，各种离散方法都有其优缺点，而且对每种方法又可分为多种形式，根据各种离散方法的特点和河道形态，为使对区域的离散既能很好地拟合长宽比很大，而且弯曲复杂的河道边界，又能根据河势、主流和水深的变化对局部区域加密细划，本模型选用了有限元法。由于离散后形成的方程组的计算量巨大，用一般计算方法对线性方程组求解不能满足要求，本模型应用质量集中的方法划系数矩阵为三对角矩阵，并用预估校正法处理、迭代求解方程组。大大减少了计算量。
2 基本方程和定解条件
2.1 基本方程
水流运动的基本方程为
| (1) |
| (2) |
式中Ui为垂线平均流速；H为水深；Z为水位；C为谢才系数。 C=1/nR1/6(水力半径R＝H)，g，ρ，ε分别为重力加速度，水密度和粘滞系数。f为科氏力系数(f=2ωsinφ，ω为地球自转角速度，φ为地理纬度)。 |
为系数矩阵泥沙运动方程为[1]
| (3) |
| (4) |
其中U为流速；ωi为泥沙沉速；S和S* 分别为水流含沙量和挟沙力，f1为泥沙非饱和系数；K1为考虑紊流脉动在水平方向产生的扩散作用及泥沙存在产生的附加影响而引入的修正系数，简称为附加系数； α*为平衡含沙量分布系数，详见[1]。
水流挟沙力计算，采用[2]中的公式计算水流挟沙力。
河床糙率计算，应用[3]中的糙率计算公式，可以描述水力泥沙因子的变化对摩阻特性的影响。
2.2 定解条件?
边界条件。对入流边界，给出水流流速(或单宽流量)和含沙量过程；对出口边界，给出水位过程线或流速(单宽流量)过程线；对固壁边界，法向流速为0，水流沿切线方向流速非0。
初始条件。给出在计算的初始时刻地形、流速、水位和含沙量等物理量的初始值。
3 有限元离散模式的建立和方程求解
3.1 有限元离散模式的建立
有限元网格的选择。根据需要，本次选用三角形常应变单元类型离散研究区域。
??设整个区域共划分为NE个单元，I个节点，单元节点总体编号为i, i=1, 2,…I。
对水沙方程(1) ～(4)，用Galerkin加权余量法逼近；对研究区域进行剖分；对单元节点和整体节点分别编号，建立局部节点编号系统和整体节点编号系统，并确定两个编号系统的关系；在离散区域的基础上，求出未知量在每个单元上的形函数；把形函数代入Galerkin积分表达式进行单元分析，建立局部有限元方程式；对所有局部有限元求和、总体合成建立总体有限元方程式，加上本质边界条件，即可得到本问题的有限元方程式。
对问题进行有限元分析
|
|
|
|
式δH，δU，δV，δS表示求解变量的变分，Ω表示域。
对以上二阶导数项利用Green公式分部积分公式，并设任一单元第i节点的平均流速，水位和含沙量及对应的加权函数的形函数值分别为Ui，Vi，Hi，Si和U*i,V*i,H*i,S*i。则有
|
|
|
|
|
|
|
|
将上述表达式代入方程可得有限元方程式
MijdZsi/dt=-Pij(HU)j-Pzij(HV)j |
MijdUj/dt=-NijUj-gP1ijZsj-Mij/ρ-ε(Qij+Rij)Uj |
MijdVj/dt=-NijVj-gP2ijZsj-Mijτyi/ρ-ε(Qij+Rij)Vj |
MijdSj/dt=-NijSj-εs(Qij+Rij)Sj-K1α*ωMij/[(f1S-S*)/H]j i, j=1, 2, 3,…I |
其中 |
|
|
|
|
|
|
|
| i,j,k=1,2,…I |
其中 |
对以上各式，当在同一项中含有两个相同的下标时，就意味着该项表示在单元体内全体同节点编号项的叠加；根据以上积分式先对各个单元分析，再叠加全域各单元方程式并使之满足边界条件即得到整体有限元方程式。
对上述有限元方程在时间上对变量Z、U、V和S用向前差分格式离散
|
|
|
|
可得到一个I×I常微分方程组。用常规方法对I×I阶方程组求解，计算量大，本文用质量集中方法[4]化方程组的系数矩阵成对角矩阵，方程组就可以容易解出。
3.2 方程的求解
用于预估校正法迭代求解方程组
第一时段的计算采用欧拉格式，第二时段采用预估校正格式，即对函数值fn进行预估、校正。用 ｜fn-f*n｜≤ε进行判断，若上式成立，则fn=f*n，否则令fn |
(fn+fn-1)/2，再由上式求f*n进行继续迭代，直到满足精度要求。当求出Ui, Vi, Hi和Si后代入有关方程，即
求出冲淤变形的河床高程
，完成方程的求解。4 模型的验证
河段基本概况 黄河下游北展滞洪区以南的北店子至后张庄河段，长约30km，是受工程控制的弯曲性河道，河道纵比降约为1?，其中北店子至泺口铁桥河段长20km，堤距一般0.7km至1.5km。泺口铁桥至后张庄长约10km，堤距1.5km～2.5km，汛期平均悬移质中值粒径d50=0.02mm～0.027mm，河床质中值粒径D50=0.07mm～0.11mm，河床糙率n=0.013～0.016。
模型区域的选择和网格划分。选萨口断面为距泺口水文站上游约20km的北店子处，进口断面距其下游的曹家圈断面、郑家店断面分别约3km和9km；出口断面选择在距泺口断面下游约10km的后张庄断面附近。计算区域内共有曹家圈、郑家店、洛口和后张庄四个大断面，两岸边界选择在险工、护滩控导工程和大堤等较稳固的工程所连接的边线上。整个计算区域划分成128个小断面，3728个三角网格单元体，共有2000个节点，三角形沿河宽方向最小长度为30m。一般沿河宽方向网格边长为主槽40m～50m左右，滩地100m左右。网格图如图1所示。
验证时段及水沙条件 选用1976年汛期8月2日至10月15日及对应的水沙过程，共75天。 初始条件 假设初始地形、流场、水位和含沙量为一定的常数。 边界条件 给出入流断面各结点的流速(单宽流量)和含沙量过程和出口断断面的水位过程线；对两岸边界，根据需要分别设定滑动和不滑动边界条件。 |
|
验证结果及其分析。
图2给出了郑家店断面和泺口断面的水位随时间的变化过程。从图中可以看出，数模计算结果与实测值符合较好。
图3给出了泺口断面流量过程与实测流量过程的比较图，可以看出，计算值与实测值基本一致。
|
图2 水位比较图 |
Water stages comparisons |
|
图3 流量比较图 |
Comparison of discharge |
表1 计算流速与实测流速比较表
Compasion of calculated and field velocities
日期 | 8.11 | 8.23 | 9.4 | 9.6 | 9.14 | 9.28 | 10.2 | |
流量(m3/s) | 1550 | 3510 | 6610 | 7410 | 5240 | 3850 | 2120 | |
平均流速 | 实测 | 2.06 | 2.62 | 2.64 | 2.83 | 2.66 | 2.41 | 2.08 |
(m/s) | 计算 | 1.86 | 2.50 | 2.50 | 2.69 | 2.53 | 2.33 | 2.00 |
最大流速 | 实测 | 3.04 | 3.61 | 3.84 | 3.86 | 3.65 | 3.26 | 2.86 |
(m/s) | 计算 | 2.76 | 3.40 | 3.56 | 3.72 | 3.60 | 3.19 | 2.70 |
|
图4不同流量流场图 |
Flow fields with different dischsrge |
表1给出了一些流量的泺口断面最大流速和最小流速的计算与实测值的比较，从图4给出了流量为2600m3/s和7800m3/s下的流场的计算结果; 图4给出了汛前小水、洪峰期大水全域流速分布，从这些图中可以看出，不同流量下河势、主流线的变化、水流漫滩、流速在主槽和滩地的不同分布情况。
|
图5 横断面冲淤变形比较图 |
Cross-sectional erosion and deposition comparisons |
图5给出了郑家店、泺口重要断面的1976年汛前6月10日实测地形、汛后10月15日数学模型计算地形与实测地形及物理模型试验的对照结果。从图中可以看出，数模计算结果在其冲淤部位和冲淤面积等方面基本上与实测数据或物理模型试验数据相符合。
5 结语
1. 模型采用了黄河泥沙研究的新成果，理论基础可靠。?
2. 采用有限元法对区域进行离散，容易处理和更好地拟合不规则河道边界；网格划分灵活，可以对局部区域任意加密；有限元程序模块可移植性强，很容易进行不同网格形状单元离散模式之间转换。
?3. 验证结果表明，在水位、流量、流速和河床冲淤等方面，数模计算结果与实测数据或物理模型试验结果符合较好，从而证明了模型的可靠性。
参 考 文 献
[1] 张红武，吕昕。 弯道水力学。 电力出版社，1993.?
[2] 张红武，张清。黄河水流挟沙力计算公式。人民黄河,1992，(1).?
[3] 赵连军，张红武。紊流阻力系数的计算。′95全国水动力学研讨会文集。海洋出版社，1995.
[4] 李东风，赖国璋。 串列双圆柱绕流的有限元模拟。力学在实践中的应用。林业出版社, 1992.5.
[5] 谢贻权，何福保。弹性和塑性力学中的有限元方法。机械出版社，1981.