CN111368410B - 一种基于河网水沙二维有限控制体积的预测方法 - Google Patents
一种基于河网水沙二维有限控制体积的预测方法 Download PDFInfo
- Publication number
- CN111368410B CN111368410B CN202010123743.3A CN202010123743A CN111368410B CN 111368410 B CN111368410 B CN 111368410B CN 202010123743 A CN202010123743 A CN 202010123743A CN 111368410 B CN111368410 B CN 111368410B
- Authority
- CN
- China
- Prior art keywords
- water
- unit
- dimensional
- vector
- interface
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 title claims abstract description 62
- 239000004576 sand Substances 0.000 title claims abstract description 25
- 238000000034 method Methods 0.000 title claims abstract description 23
- 238000004364 calculation method Methods 0.000 claims abstract description 27
- 230000004907 flux Effects 0.000 claims abstract description 22
- 239000013049 sediment Substances 0.000 claims description 22
- 230000001133 acceleration Effects 0.000 claims description 6
- 230000006872 improvement Effects 0.000 abstract description 2
- 238000004422 calculation algorithm Methods 0.000 description 6
- 238000004088 simulation Methods 0.000 description 6
- 238000011160 research Methods 0.000 description 3
- 241000282376 Panthera tigris Species 0.000 description 2
- 235000008331 Pinus X rigitaeda Nutrition 0.000 description 2
- 235000011613 Pinus brutia Nutrition 0.000 description 2
- 241000018646 Pinus brutia Species 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000011010 flushing procedure Methods 0.000 description 2
- 238000013178 mathematical model Methods 0.000 description 2
- 238000013459 approach Methods 0.000 description 1
- 230000032823 cell division Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 230000002706 hydrostatic effect Effects 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000012821 model calculation Methods 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000012876 topography Methods 0.000 description 1
Landscapes
- Complex Calculations (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种基于河网水沙二维有限控制体积的预测方法,本发明的计算步骤如下:优化网格单元设置;选用Osher格式来计算斜底跨单元界面的法向水流数值通量向量;将斜底跨单元界面的法向水流数值通量向量与单元界面的含沙量相乘,得到单元界面处的泥沙数值通量,代入有限控制体积水沙离散方程组中,本发明的有限控制体积计算的预测方法具有方法简单,易于移植,稳定性好、运行效率提高等显著优势。
Description
技术领域
本发明涉及一种基于二维有限控制体积的预测方法,尤其涉及一种基于河网水沙二维有限控制体积的预测方法,属于河流动力学计算领域。
背景技术
近年来,泥沙数学模型研究已取得了长足的进展,广泛应用于水库、河道、湖泊以及河口的河床冲淤变形预测。为了提高预测精度,除了致力于河床变形机理和计算模式的研究外,水沙运动方程组的数值离散算法的改进和创新亦是不可缺少的重要环节。河网一维水沙有限差分法因其简单实用,已成为河网水沙数值模拟的主要算法,但在实际模拟工作中存在缺陷,主要表现在:无法反映河道主槽与滩地的水沙交换量;河网汊点的分沙计算仍采用带有经验系数的计算模式,汊点分沙计算精度不够高。河道二维水沙差分算法也同样存在不够完善的方面,如有限差分法除非矩形网格很细,否则难以准确描述计算域边界和地形,不仅增加计算工作量,而且还导致局部流场和泥沙浓度场模拟精度欠佳;由于有限差分法属于节点算法,难以准确实现计算单元以及整个计算域上水量、动量以及泥沙量的平衡,影响了水沙运动数值模拟的计算精度和数值稳定性;边界条件处理繁复或不准确,内点与边界点格式不一致。数值实验表明,适用于湖泊等平底网格单元的二维水沙有限体积法的算法难以适用于河道。
发明内容
发明目的:本发明旨在提供一种方法简单、易于移植、稳定性好、运行效率提高的基于河网水沙二维有限控制体积的预测方法。
技术方案:本发明的基于河网水沙二维有限控制体积的预测方法,包括如下步骤:
(1)网格划分方式使用四边型斜低网格单元,按照CFL小于等于1.0时计算特征时间步长其中,i为第i号单元,Δxj、Δyj表示四边形第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表示第n个单位时间的守恒物理量,Wi n+1表示第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号单元,Δxj、Δyj表示四边形第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表示第n个单位时间的守恒物理量,Wi n+1表示第n+1个单位时间的守恒物理量。
进一步地,步骤(2)中,在计算Fn前,记单元界面处左右两边的y方向垂线平均的水平流速分量分别为vL、vR,计算单元i的左状态向量(hL,hLuL,hLvL)T,计算相邻单元j的右状态向量(hR,hRuR,hRvR)T。
如图1-2所示,至少需要3倍数量的平底单元,才能实现与斜底单元相近的效果。
Claims (4)
1.一种基于河网水沙二维有限控制体积的预测方法,其特征在于,包括如下步骤:
(1)网格划分方式使用四边形斜底网格单元,按照CFL小于等于1.0时计算特征时间步长其中,i为第i号单元,Δxj、Δyj表示四边形第j边的边长;ui、vi、hi分别为第i单元的x方向流速、y方向流速和水深,g为重力加速度常量;
(2)设流速在控制体单元内为常数分布,单元界面处左右两边的x方向垂线平均的水平流速分量分别为uL、uR,单元界面处左右两边的水深分别为hL、hR,则所述单元界面左右两边的单宽流量分别为qL=hLuL、qR=hRuR,选用Osher格式来计算斜底跨单元界面的法向水流数值通量向量Fn;
步骤(2)中,在计算所述Fn前,记单元界面处左右两边的y方向垂线平均的水平流速分量分别为vL、vR,计算单元i的左状态向量(hL,hLuL,hLvL)T,计算相邻单元j的右状态向量(hR,hRuR,hRvR)T;
(3)将所述斜底跨单元界面的法向水流数值通量向量Fn与单元界面的含沙量相乘,得到单元界面处的泥沙数值通量,代入有限控制体积水沙离散方程组中。
2.根据权利要求1所述的基于河网水沙二维有限控制体积的预测方法,其特征在于:步骤(2)中所述斜底跨单元界面的法向水流数值通量向量Fn采用二维不恒定浅水方程组的守恒形式计算:
其中,t为时间,守恒物理量W、x向通量向量F、y向通量向量G以及源项向量D分别为:
式中,h为水深、u是x方向垂线平均的水平流速分量、v是y方向垂线平均的水平流速分量,g为重力加速度;和/>分别为x和y方向的水底底坡;/>和/>分别为x和y方向的摩阻坡度,采用二维形式的曼宁阻力公式计算;q为单元旁侧入流。
3.根据权利要求2所述的基于河网水沙二维有限控制体积的预测方法,其特征在于:记Ωi为第i单元的区域,利用格林公式,可得所述二维不恒定浅水方程组的守恒形式的有限体积为:
式中,Ai为Ωi的面积;为Ωi的外法向单元向量;dl为线积分微元;/>为非齐次项在Ωi上的算术平均,wi为Ωi的守恒物理量。
4.根据权利要求3所述的基于河网水沙二维有限控制体积的预测方法,其特征在于:记为斜底跨单元界面的法向水流数值通量向量,时间积分采用显向前格式,所述二维不恒定浅水方程组的守恒形式的有限体积为
式中,求和号下的指标j表示某一单元i的第j边,lj为单元界面边长,Wi n表示第n个单位时间的守恒物理量,Wi n+1表示第n+1个单位时间的守恒物理量。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010123743.3A CN111368410B (zh) | 2020-02-27 | 2020-02-27 | 一种基于河网水沙二维有限控制体积的预测方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010123743.3A CN111368410B (zh) | 2020-02-27 | 2020-02-27 | 一种基于河网水沙二维有限控制体积的预测方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN111368410A CN111368410A (zh) | 2020-07-03 |
CN111368410B true CN111368410B (zh) | 2024-04-02 |
Family
ID=71208240
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010123743.3A Active CN111368410B (zh) | 2020-02-27 | 2020-02-27 | 一种基于河网水沙二维有限控制体积的预测方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111368410B (zh) |
Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107784141A (zh) * | 2016-08-31 | 2018-03-09 | 施勇 | 一种二维有限控制体积计算的加速方法 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
ES2792357T3 (es) * | 2009-04-20 | 2020-11-11 | Exxonmobil Upstream Res Co | Procedimiento para predecir el flujo de fluido |
-
2020
- 2020-02-27 CN CN202010123743.3A patent/CN111368410B/zh active Active
Patent Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107784141A (zh) * | 2016-08-31 | 2018-03-09 | 施勇 | 一种二维有限控制体积计算的加速方法 |
Non-Patent Citations (1)
Title |
---|
施勇."长江中下游水沙输运及其调控数学模型研究".《中国博士学位论文全文数据库》.2007,第60-95页. * |
Also Published As
Publication number | Publication date |
---|---|
CN111368410A (zh) | 2020-07-03 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN110543731B (zh) | 一种高含沙河流水库水沙分置调度模拟方法及*** | |
CN107045568B (zh) | 基于动态规划逐次逼近法的河道糙率反演方法 | |
CN106599374B (zh) | 一种适用于径流潮汐河口段的河相关系及其推导方法 | |
CN109271672B (zh) | 一种河-湖-泵站相互影响作用下的河道水面线计算方法 | |
CN110210109B (zh) | 一种河网中堰闸工程反向水流的数值模拟方法及*** | |
CN104933268B (zh) | 一种基于一维非恒定流数值模型的洪水分析方法 | |
CN103886141A (zh) | 一种自动批量计算河道断面信息的方法 | |
CN111400974B (zh) | 一种估算矩形河渠壁面和床面切应力的方法 | |
CN111597732B (zh) | 一种使用汊点影响区水面梯度的河网水流数值模拟方法 | |
CN107256312B (zh) | 一种基于潮流环境下海湾纳潮变化量计算方法 | |
CN111368410B (zh) | 一种基于河网水沙二维有限控制体积的预测方法 | |
CN106320255B (zh) | 一种漫滩水流滩槽过流量的计算方法 | |
CN109145396B (zh) | 一种基于植被分布的河道糙率分区率定方法 | |
CN106383997A (zh) | 一种反推三峡水库区间入流过程的计算方法 | |
CN114330156A (zh) | 一种基于中心圆柱矩形槽比尺效应的渠道测流法 | |
CN110147646A (zh) | 一种二维非结构洪水数值模拟框架下线性挡水构筑物的过流处理方法 | |
CN107784141B (zh) | 一种二维有限控制体积计算的加速方法 | |
CN114357903B (zh) | 基于河床形态参数的冲积河流阻力计算、模拟方法及装置 | |
CN116933685A (zh) | 一种河道型水库调控方法、装置、电子设备及存储介质 | |
CN114595596B (zh) | 一种水库淤积物密实沉降速率的模拟方法 | |
Duan et al. | Improved 2D shallow water model able to capture the effects of complex bathymetric features through their subgrid modeling | |
CN115564238A (zh) | 一种基于水动力机制的小水库洪水预报调洪分析方法 | |
CN114547869B (zh) | 一种二维非结构干滩条件下流量边界的处理方法 | |
Han et al. | Numerical experiments on stagnation points influenced by the Three Gorges Dam in the Yangtze Estuary | |
CN117010727B (zh) | 一种基于水平年的径流序列一致性改正方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |