天然河道一维水动力高精度数值模拟方法
阅读说明：本技术 天然河道一维水动力高精度数值模拟方法 (High-precision numerical simulation method for one-dimensional water power of natural river ) 是由 孙万光 马军 栾宇东 范宝山 杨海滔 唐振华 姜彪 蒋攀 赵浩 于 2021-08-26 设计创作，主要内容包括：本发明提供了天然河道一维水动力高精度数值模拟方法,涉及计算流体动力学技术领域,包括守恒型圣维南方程、采用Godunov格式的有限体积法对控制方程离散、HLLC近似Riemann求解器、变量空间重构、采用二阶龙格—库塔离散的时间分裂方法来处理源项,本发明基于Godunov格式,提出了HLLC近似Riemann求解器的通量计算方法,并拓展至守恒型圣维南方程；提出了针对天然河道复杂断面几何形状下的变量空间重构方法：依据过流断面面积和静力矩等效原则将河道断面概化成矩形,通过线性插值构造单元界面处断面几何形状,根据水位重构结果计算界面两侧过流断面面积和静力矩的重构值,保证计算格式守恒。本发明为天然河道水动力及环境水力学高精度数值模拟提供了新的方法。(The invention provides a high-precision numerical simulation method for one-dimensional hydrodynamic force of a natural river channel, which relates to the technical field of computational fluid dynamics and comprises a conservation type Saint-Venn equation, a finite volume method adopting a Godunov format to discretize a control equation, an HLLC approximate Riemann solver, a variable space reconstruction and a time splitting method adopting second-order Runge-Kutta discretization to process source items, wherein the invention provides a flux calculation method of the HLLC approximate Riemann solver based on the Godunov format and expands the flux calculation method to the conservation type Saint-Venn equation; a variable space reconstruction method aiming at the complex cross section geometry of a natural river channel is provided: the river channel section is generalized into a rectangle according to the equivalent principle of the flow cross-section area and the static moment, the cross-section geometry at the unit interface is constructed through linear interpolation, the reconstruction values of the flow cross-section area and the static moment on the two sides of the interface are calculated according to the water level reconstruction result, and the conservation of the calculation format is ensured. The invention provides a new method for high-precision numerical simulation of natural river hydrodynamic force and environmental hydraulics.)
技术领域
本发明涉及计算流体动力学技术领域，具体为天然河道一维水动力高精度数值模拟方法。
背景技术
天然河道一维非恒定流数值模拟在水利规划设计以及防洪减灾等领域应用十分广泛；天然河道特别是山区性河道断面几何形状快速变化、河底比降变化大，水流流态可能出现急流、缓流交替现象；另外，高纬度严寒地区江河冰塞及冰坝现象频繁发生，冰坝溃决后坝下将出现激波间断，这些都属于浅水间断流动，因间断处水力要素突变，具有较强非线性，传统数值计算方法(如有限差分)往往失效。
对于间断问题，控制方程上，守恒数值格式的解可收敛至守恒控制方程的弱解，采用非守恒型控制方程求解间断问题将会得到错误结果，这些方程仅对光滑流动有效，对于间断问题，根据跳跃条件计算出的激波波速是错误的；数值计算方法上，Godunov格式因具备较强的不连续问题处理能力而广受欢迎。在Godunov格式下，将齐次浅水方程组的Riemann问题定义为具有分段变量恒定的初值问题，采用精确Riemann求解器可求得单元界面处的通量值，精度高，但需要对非线性代数恒等式进行迭代求解，增加了计算资源的消耗。为了简化计算，多种近似Riemann求解器得到快速发展，如Roe法和HLL法等， Roe线性化Riemann求解器缺点如下：①在跨临界流或激波模拟中，为避免非物理解必须进行熵修正；②在强稀疏波作用下(水流非常浅的区域)，线性化Riemann求解器计算水深出现负值；③在强波相互作用下，线性化Riemann解一般缺乏鲁棒性。HLL近似Riemann求解器将波族简化为2个，只适用于具有2个方程的一维系统，对接触间断和剪切波模拟精度较差。HLLC求解器为3波族模型，可准确求解接触间断和剪切波问题；对比HLL求解器，HLLC求解器可扩充浓度组分方程，为天然河道环境水力学数值模拟创造了有利条件，推广应用前景十分广阔，以往的研究中，多采用HLLC求解器对浅水方程进行求解，而天然河道一维水动力通常以圣维南方程作为控制方程，采用 HLLC求解器对圣维南方程进行求解至今少见报道，与浅水方程相比，守恒型圣维南方程需考虑河宽及其沿程变化产生的影响，实际应用中，需重新推导基于守恒型圣维南方程的HLLC求解器中通量计算公式。
Godunov格式下，为了获得高精度数值解，一般需要进行变量空间重构，以MUSCL(Monotonic Upstream-centered Scheme for Conservation Laws)为代表的变量空间重构方法应用比较广泛，但该方法仅适用于断面渐变条件下的变量空间重构，天然河道断面几何形状复杂，河宽和水深沿程快速变化，直接采用MUSCL法进行变量空间重构会产生较大误差，基于此，本发明设计了天然河道一维水动力高精度数值模拟方法，以解决上述问题。
发明内容
本发明的目的在于提供天然河道一维水动力高精度数值模拟方法，采用守恒型圣维南方程作为天然河道一维非恒定流控制方程，基于 Godunov格式，提出了天然河道断面几何形状快速变化条件下的变量空间重构方法，推导了基于守恒型圣维南方程的HLLC求解器通量计算公式，为天然河道复杂水动力数值模拟提供一种高精度、简便的方法。
为实现上述目的，本发明提供如下技术方案：天然河道一维水动力高精度数值模拟方法，包括如下步骤：
S1:守恒型圣维南方程；
S2:采用Godunov格式的有限体积法对控制方程离散；
S3:基于守恒型圣维南方程的HLLC近似Riemann求解器；
S4:变量空间重构；
S5:采用二阶龙格—库塔离散的时间分裂方法来处理源项。
优选的，采用守恒型圣维南方程组，表达式如下：
式中：U为变量；F为通量；S为源项；t为时间；x为空间坐标；A＝ A(x,t)为过流断面面积；Q为流量；g为重力加速度；S0为床面比降；Sf为摩阻比降；I1和I2分别为静力矩和侧压力，表达式如下：
式中：h为水深，b(x,η)为断面宽度。
优选的，提出了基于守恒型圣维南方程的HLLC近似Riemann求解器的通量计算方法，将该求解器由浅水方程拓展至守恒型圣维南方程：
式中：UL和UR分别为界面左侧和右侧变量；FL和FR分别为界面左侧和右侧通量；U*L和U*R分别为中间波左侧和右侧变量，为待求变量；F*L和F*R分别为中间波左侧和右侧通量；SL和SR分别为界面左侧和右侧波速。
为了计算界面通量，还需U*L和U*R的值；这里引入如下假设：
A*L＝A*R＝A*,Q*L＝Q*R＝Q*,S*＝Q*/A* (11)
事实上，上述假设对于精确Riemann求解器也是正确的。推导得出基于守恒型圣维南方程的中间波变量U*＝[A*,Q*]T以及中间波波速S*的表达式：
对跨越波速SL、SR和S*的状态分别运用Rankine-Hugoniot条件，可得：
将式(2)中第一个分量(即质量守恒方程)带入式(12)中前两项，并将式(11)带入，可得：
上式是关于Q*和A*的方程组，经求解可得：
优选的，提出了针对天然河道复杂断面几何形状下的变量空间重构方法：依据过流断面面积和静力矩等效原则将河道断面概化成矩形，通过线性插值构造单元界面处断面几何形状，根据水位重构结果计算界面两侧过流断面面积和静力矩的重构值，保证计算格式守恒。具体步骤如下：
a)天然河道断面概化方法
给定河道水位，首先计算过流断面面积A和静力矩I1，依据过流断面面积和静力矩等效原则将天然河道断面概化成矩形，即：A′＝A、 I′1＝I1，其中A′、I′1分别为等效过流断面面积和等效静力矩。由此可推导出过流断面等效水深h′、等效河宽B′和等效河底高程Z′b，表达式如下：
水深是波速计算的关键参数，天然河道通常采用h＝A/B来计算断面平均水深，其中B为水面宽度,然而对于复式断面，B和h均存在突变，对变量空间重构产生明显影响，而等效水深和等效河宽则不存在突变，此概化方法保证单元水体受力条件不变，物理概念清晰，等效河宽和等效水深的概化结果唯一；
b)单元界面处断面几何形状构造
根据概化后相邻河道断面的宽度和底高程，通过线性插值构造单元界面处断面几何形状；以i+1/2界面处断面几何形状构造为例，表达式如下：
c)变量空间重构
需空间重构的变量包括A、Q和I1；先采用MUSCL法构造界面处的水位和流量然后基于构造界面处过流断面面积和静力矩，表达式如下：
与现有技术相比，本发明的有益效果是：本发明基于Godunov格式，提出了基于守恒型圣维南方程的HLLC近似Riemann求解器的通量计算方法，将该求解器由浅水方程拓展至守恒型圣维南方程；提出了针对天然河道复杂断面几何形状下的变量空间重构方法：依据过流断面面积和静力矩等效原则将河道断面概化成矩形，通过线性插值构造单元界面处断面几何形状，根据水位重构结果计算界面两侧过流断面面积和静力矩的重构值，保证计算格式守恒，本发明为天然河道水动力及环境水力学高精度数值模拟提供了新的方法。
当然，实施本发明的任一产品并不一定需要同时达到以上所述的所有优点。
附图说明
为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为中心格式有限体积法单元离散；
图2为HLLC近似Riemann求解器波结构；
图3为复式断面示意图；
图4为断面平均水深和等效水深对比；
图5为河道平面示意图及断面布置；
图6为河道典型实测断面及概化断面图；
图7为断面最高水位计算值与实测值对比
附图中，各标号所代表的部件列表如下：
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。
请参阅图1-7，本发明提供一种技术方案：天然河道一维水动力高精度数值模拟方法，包括如下步骤：
S1:守恒型圣维南方程；
S2:采用Godunov格式的有限体积法对控制方程离散；
S3:基于守恒型圣维南方程的HLLC近似Riemann求解器；
S4:变量空间重构；
S5:采用二阶龙格—库塔离散的时间分裂方法来处理源项。
其中，采用守恒型圣维南方程组，表达式如下：
式中：U为变量；F为通量；S为源项；t为时间；x为空间坐标；A＝ A(x,t)为过流断面面积；Q为流量；g为重力加速度；S0为床面比降；Sf为摩阻比降；I1和I2分别为静力矩和侧压力，表达式如下：
式中：h为水深，b(x,η)为断面宽度。
其中，提出了基于守恒型圣维南方程的HLLC近似Riemann求解器的通量计算方法，将该求解器由浅水方程拓展至守恒型圣维南方程：
式中：UL和UR分别为界面左侧和右侧变量；FL和FR分别为界面左侧和右侧通量；U*L和U*R分别为中间波左侧和右侧变量，为待求变量；F*L和F*R分别为中间波左侧和右侧通量；SL和SR分别为界面左侧和右侧波速；
为了计算界面通量，还需U*L和U*R的值；这里引入如下假设：
A*L＝A*R＝A*,Q*L＝Q*R＝Q*,S*＝Q*/A* (11)
事实上，上述假设对于精确Riemann求解器也是正确的。推导得出基于守恒型圣维南方程的中间波变量U*＝[A*,Q*]T以及中间波波速S*的表达式：
对跨越波速SL、SR和S*的状态分别运用Rankine-Hugoniot条件，可得：
将式(2)中第一个分量(即质量守恒方程)带入式(12)中前两项，并将式(11)带入，可得：
上式是关于Q*和A*的方程组，经求解可得：
其中，提出了针对天然河道复杂断面几何形状下的变量空间重构方法：依据过流断面面积和静力矩等效原则将河道断面概化成矩形，通过线性插值构造单元界面处断面几何形状，根据水位重构结果计算界面两侧过流断面面积和静力矩的重构值，保证计算格式守恒。具体步骤如下：
a)天然河道断面概化方法
给定河道水位，首先计算过流断面面积A和静力矩I1，依据过流断面面积和静力矩等效原则将天然河道断面概化成矩形，即：A′＝A、 I′1＝I1，其中A′、I′1分别为等效过流断面面积和等效静力矩，由此可推导出过流断面等效水深h′、等效河宽B′和等效河底高程Z′b，表达式如下：
水深是波速计算的关键参数，天然河道通常采用h＝A/B来计算断面平均水深，其中B为水面宽度,然而对于复式断面，B和h均存在突变，对变量空间重构产生明显影响，而等效水深和等效河宽则不存在突变，此概化方法保证单元水体受力条件不变，物理概念清晰，等效河宽和等效水深的概化结果唯一；
b)单元界面处断面几何形状构造
根据概化后相邻河道断面的宽度和底高程，通过线性插值构造单元界面处断面几何形状；以i+1/2界面处断面几何形状构造为例，表达式如下：
c)变量空间重构
需空间重构的变量包括A、Q和I1；先采用MUSCL法构造界面处的水位和流量然后基于构造界面处过流断面面积和静力矩，表达式如下：
本实施例的一个具体应用为：
某河流经盆地和峡谷，河道平面形态变化较大，常水位下河宽变化范围为200m～1000m，河道平面示意图见图5，此河段长度 59.44km，布置实测大断面18条(0#～17#)。
河道典型横断面见图6；盆地河段为复式断面，河宽较宽，而峡谷段为“V”字形断面，河宽较小，河床纵断面起伏剧烈，深泓最低点高程-48.67m，深泓最高点高程为33.3m。该河段1994年发生了频率P＝4％洪水，洪峰流量为44400m3/s，洪水过后对各断面洪痕进行了准确测量，该场次洪水水面线数据完备、准确，可以此作为河道糙率率定以及计算结果验证的依据。
天然河道一维水动力高精度数值模拟方法步骤如下：
①按照式(19)～(21)对每个单元的基本变量进行线性重构，得到界面处各变量值；
②对浅水波波速进行计算；
③按照式(7)、式(14)对界面通量进行计算；
④计算下个时段的A、Q。
将断面计算最高水位与实测洪痕进行对比，结果见图7(本发明方法简称“HLLC”)，为了便于对比，基于HEC-RAS软件，采用相同的断面数据和边界条件，建立一维非恒定流数值计算模型，对1994年典型洪水过程进行计算(简称“HEC-RAS”)；HEC-RAS软件非恒定流控制方程为非守恒型，模型采用有限差分进行求解，文献中非恒定流控制方程同为非守恒型，基于Godunov格式，采用HLL近似 Riemann求解器求解，采用该文方法同步进行对比(简称“HLL”)。
从图7中可见，基于HLLC算法的计算结果与实测值吻合良好，而HEC-RAS和HLL法模拟结果与实测值有较大偏差，偏差的起始点位于9#断面，9#～8#断面河道水面比降明显高于实测值，9#断面上游水位计算结果均高于实测值，其中HEC-RAS法偏差最大可达1.69m(位于11#断面)，HLL法偏差最大可达1.42m(位于10#断面)。
天然河道断面几何形状快速变化，断面间水深和过流断面面积关系差别较大，以8#断面和9#断面为例(见图6)，最高水位条件下8# 断面过流断面面积为9595m2，等效水深为47.13m；而9#断面过流断面面积为17133m2，为8#断面的1.79倍，等效水深为33.9m，为8#断面的0.72倍，若直接进行变量空间重构则会产生较大误差，通量计算结果严重失真，同时也会产生格式不守恒问题。人为构造单元界面处断面几何形状，采用界面处基于MUSCL法的水位重构值计算与之对应的过流断面面积和静力矩，以此作为变量空间重构结果，成功避开了断面几何形状快速变化对变量空间重构的影响，同时保证了计算格式的守恒。
在本说明书的描述中，参考术语“一个实施例”、“示例”、“具体示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。
以上公开的本发明优选实施例只是用于帮助阐述本发明。优选实施例并没有详尽叙述所有的细节，也不限制该发明仅为所述的具体实施方式。显然，根据本说明书的内容，可作很多的修改和变化。本说明书选取并具体描述这些实施例，是为了更好地解释本发明的原理和实际应用，从而使所属技术领域技术人员能很好地理解和利用本发明。本发明仅受权利要求书及其全部范围和等效物的限制。