一种基于河网水沙二维有限控制体积的预测方法与流程

本发明涉及一种基于二维有限控制体积的预测方法，尤其涉及一种基于河网水沙二维有限控制体积的预测方法，属于河流动力学计算领域。

背景技术：

近年来，泥沙数学模型研究已取得了长足的进展，广泛应用于水库、河道、湖泊以及河口的河床冲淤变形预测。为了提高预测精度，除了致力于河床变形机理和计算模式的研究外，水沙运动方程组的数值离散算法的改进和创新亦是不可缺少的重要环节。河网一维水沙有限差分法因其简单实用，已成为河网水沙数值模拟的主要算法，但在实际模拟工作中存在缺陷，主要表现在：无法反映河道主槽与滩地的水沙交换量；河网汊点的分沙计算仍采用带有经验系数的计算模式，汊点分沙计算精度不够高。河道二维水沙差分算法也同样存在不够完善的方面，如有限差分法除非矩形网格很细，否则难以准确描述计算域边界和地形，不仅增加计算工作量，而且还导致局部流场和泥沙浓度场模拟精度欠佳；由于有限差分法属于节点算法，难以准确实现计算单元以及整个计算域上水量、动量以及泥沙量的平衡，影响了水沙运动数值模拟的计算精度和数值稳定性；边界条件处理繁复或不准确，内点与边界点格式不一致。数值实验表明，适用于湖泊等平底网格单元的二维水沙有限体积法的算法难以适用于河道。

技术实现要素：

发明目的：本发明旨在提供一种方法简单、易于移植、稳定性好、运行效率提高的基于河网水沙二维有限控制体积的预测方法。

技术方案：本发明的基于河网水沙二维有限控制体积的预测方法，包括如下步骤：

(1)网格划分方式使用四边型斜低网格单元，按照cfl小于等于1.0时计算特征时间步长其中，i为第i号单元，δx^j、δy^j表示四边形第j边的边长(j＝1～4)；ui、vi、hi分别为第i单元的x方向流速、y方向流速和水深，g为重力加速度常量，相同面积的四边形的时间步长几乎是三角形的时间步长的1倍；

(2)设流速在控制体单元内为常数分布，单元界面处左右两边的x方向垂线平均的水平流速分量分别为ul、ur，单元界面处左右两边的水深分别为hl、hr，则单元界面左右两边的单宽流量分别为ql＝hlul、qr＝hrur，选用osher格式来计算斜底跨单元界面的法向水流数值通量向量fn；

(3)将斜底跨单元界面的法向水流数值通量向量fn与单元界面的含沙量相乘，得到单元界面处的泥沙数值通量，代入有限控制体积水沙离散方程组中。

对于平底网格单元需要加密2～3倍单元数才能逼近滩地斜坡，所以采用斜底网格单元有利于减少网格单元数量，且利于提高计算效率。

进一步地，步骤(2)中斜底跨单元界面的法向水流数值通量向量fn采用二维非恒定浅水方程组的守恒形式计算：

其中，t为时间，守恒物理量w、x向通量向量f、y向通量向量g以及源项向量d分别为：

式中，h为水深、u是x方向垂线平均的水平流速分量、v是y方向垂线平均的水平流速分量，g为重力加速度；和分别为x和y方向的水底底坡；和分别为x和y方向的摩阻坡度，采用二维形式的曼宁阻力公式计算；q为单元旁侧入流。

采用二维不恒定浅水方程组的守恒形式计算是为了保证格式的守恒性，以及适用于含间断或陡梯度的流动；采用有限体积法对水沙微分方程组进行数值离散，其实质是逐单元进行水量、动量和沙量平衡，物理意义清晰，准确满足积分形式的守恒律，成果无守恒性误差，且能够处理含间断或陡梯度的流动。

进一步地，记ωi为第i单元的区域，利用格林公式，可得二维不恒定浅水方程组的守恒形式的有限体积为：

式中，ai为ωi的面积；为ωi的外法向单元向量；dl为线积分微元；为非齐次项在ωi上的算术平均，wi为ωi的守恒物理量。

进一步地，记为斜底跨单元界面的法向水流数值通量，时间积分采用显向前格式，二维不恒定浅水方程组的守恒形式的有限体积为

式中，求和号下的指标j表示某一单元i的第j边，lj为单元界面边长，wiⁿ表示第n个单位时间的守恒物理量，wiⁿ⁺¹表示第n+1个单位时间的守恒物理量。

进一步地，步骤(2)中，在计算fn前，记单元界面处左右两边的y方向垂线平均的水平流速分量分别为vl、vr，计算单元i的左状态向量(hl,hlul,hlvl)^t，计算相邻单元j的右状态向量(hr,hrur,hrvr)^t。

在控制体单元的动量平衡中，按照斜底面坡的水深变化计算静水压力项，避免底床不平引起的反力作用。模拟试验表明，采用斜底单元的二维水沙计算模式，明显减少网格单元数目，提高了计算效率，改善了河网河道流场计算稳定性和提高了模拟精度。

有益效果：与现有技术相比，本发明具有如下显著优点：

本发明的算法简单，易于移植；而且稳定性好、运行效率大幅提高。

附图说明

图1为本发明的斜底单元；

图2为斜底单元和平底单元的对比图；

图3松虎水系网格图。

具体实施方式

下面结合实施例对本发明的技术方案作进一步说明。

以长江中游洞庭湖区松虎河网水系为例，松虎二维水沙计算模块，见图3，以松滋口、太平口和澧水津市为入流进沙条件，目平湖水位为下边界的二维河道泥沙冲淤计算模块；松虎水系共划分为3032个四边形网格单元，其中，主槽部分沿流向方向网格单元最小边长250m、最大边长800m，滩地部分最小边长400m、最大边长1200m。

按照上述单元划分，得到特征时间步长如表1。

表1单元网格特征时间步长单位：秒

根据上述特征时间步长，选择等于特征时间步长为模型计算时间步长，考虑水沙长系列计算的效率需求，故针对特征时间步长较小的单元进行再一次的调整，即适当加大小单元边长、减小大单元边长如表2，使得网格单元的特征时间步长尽可能均衡，经过调整后计算时间步长加大1.5倍，且计算稳定性得到了提高。

表2调整后的单元网格特征时间步长单位：秒

本实施例的基于河网水沙二维有限控制体积的预测方法，包括如下步骤：

(1)网格划分方式使用四边型斜低网格单元，按照cfl为小于等于1.0时计算特征时间步长其中，i为第i号单元，δx^j、δy^j表示四边形第j边的边长(j＝1～4)；ui、vi、hi分别为第i单元的x方向流速、y方向流速和水深，g为重力加速度常量。经过按照特征时间步长调整网格单元边长见表2，使得特征时间步长较为均衡，能够较好提高计算稳定性和计算效率。

(2)设流速在控制体单元内为常数分布，单元界面处左右两边的x方向垂线平均的水平流速分量分别为ul、ur，单元界面处左右两边的水深分别为hl、hr，则单元界面左右两边的单宽流量分别为ql＝hlul、qr＝hrur，选用osher格式来计算斜底跨单元界面的法向水流数值通量向量fn；

(3)将斜底跨单元界面的法向水流数值通量向量fn与单元界面的含沙量相乘，得到单元界面处的泥沙数值通量，代入有限控制体积水沙离散方程组中。

进一步地，步骤(2)中斜底跨单元界面的法向水流数值通量向量fn采用二维非恒定浅水方程组的守恒形式计算：

其中，t为时间，守恒物理量w、x向通量向量f、y向通量向量g以及源项向量d分别为：

式中，h为水深、u是x方向垂线平均的水平流速分量、v是y方向垂线平均的水平流速分量，g为重力加速度；和分别为x和y方向的水底底坡；和分别为x和y方向的摩阻坡度，采用二维形式的曼宁阻力公式计算；q为单元旁侧入流。

进一步地，记ωi为第i单元的区域，利用格林公式，可得二维不恒定浅水方程组的守恒形式的有限体积为：

式中，ai为ωi的面积；为ωi的外法向单元向量；dl为线积分微元；为非齐次项在ωi上的算术平均，wi为ωi的守恒物理量。

进一步地，记为斜底跨单元界面的法向水流数值通量，时间积分采用显向前格式，二维不恒定浅水方程组的守恒形式的有限体积为

式中，求和号下的指标j表示某一单元i的第j边，lj为单元界面边长，wiⁿ表示第n个单位时间的守恒物理量，wiⁿ⁺¹表示第n+1个单位时间的守恒物理量。

进一步地，步骤(2)中，在计算fn前，记单元界面处左右两边的y方向垂线平均的水平流速分量分别为vl、vr，计算单元i的左状态向量(hl,hlul,hlvl)^t，计算相邻单元j的右状态向量(hr,hrur,hrvr)^t。

如图1-2所示，至少需要3倍数量的平底单元，才能实现与斜底单元相近的效果。