CN110705188B - 一维冰水耦合运动的高精度格式模拟方法 - Google Patents

Abstract

本发明涉及一维冰水耦合运动的高精度格式模拟方法，属于江河冰凌灾害预报领域。包括获取边界条件，控制方程，冰水耦合运动统一离散，面通量计算，冰水耦合条件下浅水波波速计算和高精度格式。本发明将冰与水动力学方程统一离散成高精度Godunov格式，采用HLL近似Riemann解来计算界面通量，采用等效水深代替断面平均水深，空间上对每个单元基本变量进行分段线性重构，时间上采用预测—校正方法，使数值解的精度整体提高到二阶。本发明提出的冰水耦合运动的高精度格式模拟方法对流冰期河道水位和冰厚有较高的预报精度，特别是采用等效水深法代替平均水深法计算浅水波波速，可明显提高水位预报精度，研究成果对于河道冰凌灾害防治意义重大。

Description

一维冰水耦合运动的高精度格式模拟方法

技术领域

本发明涉及江河冰凌灾害预报领域，特别涉及江河冰塞、冰坝发生位置以及冰塞、冰坝水位预报方法。

背景技术

我国高纬度严寒地区河道冰情严重，冰坝凌汛灾害频繁发生。河道内形成冰坝后，河道水位迅速升高，各别年份甚至超过历史最高洪水位，严重威胁河流沿岸居民的生命和财产安全。另外，流冰动量大，可使河道内水工建筑物或涉河建筑物遭到严重破坏。深入研究河冰动力输移的基本规律，对冰坝凌汛灾害防治具有十分重要的意义。

在封河和开河时段，水面冰塞过程模拟是河冰模拟的重要组成部分。冰塞模型分为静态冰塞模型和冰动力输移模型。Pariset和Hausser将水面冰塞划分为并置型、窄河阻塞型、宽河阻塞型三种类型，发展了静态冰塞模型理论。许多研究机构相继开发了多款一维静态非平衡冰塞模型并应用到工程分析中，如ICEJAM、RIVJAM、HEC-RAS等。因没有考虑冰的动力特性，静态冰塞模型不能确定何时、何地发生冰塞。此外，冰和水流的运动以及非恒定作用对冰塞演变和冰厚影响显著。Shen开发了一维河冰动力输移模型，将流冰内部应力计算划分为急流和缓流两种类型。Lal和Shen开发了基于MacCormack方法的一维河冰动力输移数值计算模型，由于冰层中存在压缩波，冰塞前缘流冰密度和冰厚变化剧烈，模型需要采用高精度格式数值计算方法进行模拟。Zufelt和Ettema对一维冰水耦合运动方程进行求解，冰应力计算采用静态莫尔—库伦准则，离散方法采用四点隐格式有限差分法。

现有的一维冰水耦合动力学数值模拟方法中，将冰和水分别采用不同方法进行离散，同时采用不同时间步长的两个方案进行求解，例如水动力方程通常采用四点隐格式有限差分法进行离散，冰动力学方程通常采用MacCormack方法进行离散，为了抑制数值震荡，需加入人工黏性项。这样的处理方式存在几个问题：①四点隐格式有限差分法具有一阶精度，对冰塞溃决等复杂水流条件模拟精度较差；②人工黏性项缺乏明确的物理意义，人为因素影响较大；③浅水波波速决定水动力界面通量，而冰的运动特性又直接影响浅水波波速，冰与水之间存在紧密耦合关系，而以往的研究成果忽略了冰对浅水波波速的影响，冰与水之间被认为是松散的耦合关系。

发明内容

本发明的提供一种一维冰水耦合运动的高精度格式模拟方法，用高精度Godunov格式对冰水耦合动力学方程进行统一离散并模拟，对流冰条件下的浅水波波速进行推导，实现了冰与水运动的紧密耦合和高精度模拟，以期为防凌减灾工作提供技术支撑。

本发明采取的技术方案是，包括下列步骤：

(1)边界条件

获取计算河段河道横断面数据，上边界流量Z_{_boundary}和水位Q_{_boundary}观测数据，上边界冰的密集度N_{_boundary}、冰厚t_{i_boundary}以及冰流速u_{_boundary}观测数据；

(2)控制方程

水动力方程包含的未知数为Z和Q，冰动力方程求解的未知数包括N、t_i、u，水动力与冰动力控制方程用守恒矢量型式来表示：

其中：

式中，A为水流断面面积；Q为流量；Z为水位；g为重力加速度；p_b、p_i分别为河床和岸冰湿周；τ_b、τ_w分别为在河底和冰—水界面上的剪切应力；τ_i为岸冰的剪切应力；τ_a′_x为沿河道风对水面的拖曳力；ρ_i为冰的质量密度；x为空间坐标；N为冰的密集度；t_i为冰厚；u为冰的流速；B为水面冰层宽度；E_m为由热力过程所引起的表层冰质量的变化率、包括降雪以及浮冰在水中的交换；τ_w为沿河水的应力；τ_ax为沿河风的应力；τ_B为河岸剪切力；Z为水位；σ_x为沿水平方向上冰的平均应力；

为水流方向上重力分量；E_a为由于动力过程和悬浮冰花交换所引起的冰面积密度的时间变率；R_a为由机械再分配引起的冰面积密度的时间变率；

(3)冰水耦合运动统一离散方法

采用中心格式的有限体积法对控制方程进行离散，设第i个单元的长度为△x_i，对水动力和冰动力统一型式的守恒方程进行长度积分可得：

式中，D_i表示i单元内D的平均值，利用Green定理可得离散方程如下：

式中，U_i为单元中心点的向量，为整个单元的平均值；F_i+1/2和F_i-1/2分别为i+1/2和i-1/2界面处的通量；

(4)界面通量计算方法

根据Godunov方法，每个单元中的变量近似为恒定值，界面通量可通过Riemann问题进行计算；采用HLL近似Riemann解来计算界面通量，该方法稳定并且计算简便，界面通量定义如下：

式中，S_L和S_R分别为左侧和右侧波速；

当S_L<0<S_R时，水和冰的界面通量统一表达式，

具体如下：

浅水水流波速及冰压缩波波速分别通过下式计算：

式中，V_L和V_R分别为界面左侧和右侧水流流速；

和

分别为界面左侧和右侧水力半径；u_mL和u_mR分别为界面左侧和右侧修正冰流速；s_L和s_R分别为界面左侧和右侧冰压缩波波速；

(5)冰水耦合条件下浅水波波速计算方法

浅水波波速是界面通量计算的关键变量，在有无冰体存在以及不同冰速条件下，水流运动所受阻力完全不同，可直接影响浅水波波速；在固定冰盖条件下，根据Einstein阻力划分原理，按零应力面将冰盖下过水断面划分为冰盖区A_i和床面区A_b；提出等效水深的定义：为了使冰盖下水流运动与畅流期水流运动具有相同的阻力形式，将冰盖区按照零应力面展开，此时的断面平均水深称为等效水深；宽浅型河道，等效水深h近似等于断面综合水力半径R_c，表达式如下：

式中，B_i为A_i对应的宽度；B_b为A_b对应的宽度；P_i和P_b分别为冰盖区湿周和床面区湿周；

浅水波波速计算表达式如下：

综合水力半径R_c计算表达式如下：

式中，α＝P_b/P_i，在流冰条件下，A_i、R_i以及P_i均会随冰速u变化而变化，冰盖区A_i和床面区A_b之间仍然存在零应力面，因此等效水深的概念仍然适用，当V_w<＝u时，A_i＝0，P_i＝0；当V_w>u时，A_i/A_b的表达式如下：

式中，

R′_i＝A_i/B，进一步推导P_i，表达式如下：

Q_i＝Q-Q_b  (24)

当u→V_w时，

是比

高阶的无穷小，因此P_i→0，并且α→∞，进而可得R_c＝R_b，这表明，当冰速接近断面平均流速时，A_i、P_i都趋于0，流冰对水流的阻力也趋于0，计算表达式与事实相符，在求得A_i和P_i后，进而可求得R_i及α，根据式(22)计算综合水力半径R_c，根据式(21)计算浅水波波速；

(6)高精度格式

HLL近似Riemann解可获得一阶精度数值解，对每个单元基本变量进行分段线性重构后，可获得空间二阶精度的数值解；

为了避免数值震荡，在估计坡度

时必须使用限制器，采用minmod限制器，2个变量，当其符号相同时取其最小值，否则取0，因此i单元的坡度表达式如下：

时间步上引入Hancock预测—校正方法，使数值解的精度整体提高到二阶，表达式如下：

通过式(30)的迭代可获得各变量的各个时间步和空间步的计算结果。

源项中，水面坡度的二阶精度表达式如下：

式中，

源项中的其它项采用显格式处理。

本发明优点是：将冰与水动力学方程统一离散成高精度Godunov格式，采用HLL近似Riemann解来计算界面通量，提出了水和冰的界面通量统一表达式；在浅水波波速计算中，采用等效水深代替断面平均水深，推导得出了流冰条件下浅水波波速的计算方法，实现了冰水运动的紧密耦合。空间上对每个单元基本变量进行分段线性重构，时间上采用预测—校正方法，使数值解的精度整体提高到二阶，取代了为抑制数值震荡而添加的人工黏性项，消除了人为影响因素。实例研究表明，本发明提出的冰水耦合运动的高精度格式模拟方法对流冰期河道水位和冰厚有较高的预报精度，特别是采用等效水深法代替平均水深法计算浅水波波速，可明显提高水位预报精度，研究成果对于河道冰凌灾害防治意义重大。

附图说明

图1是等效水深示意图；

图2是计算断面布置；

附图3是洛古河水文站断面水位计算值与实测值对比；

附图4是水面宽度的沿程变化；

附图5是冰厚的沿程变化；

附图6是冰速的沿程变化；

附图7是水位的沿程变化；

附图8是不同浅水波波速计算方法下沿程水位对比图。

具体实施方式

包括下列步骤：

(1)边界条件

获取计算河段河道横断面数据，上边界流量Z_{_boundary}和水位Q_{_boundary}观测数据，上边界冰的密集度N_{_boundary}、冰厚t_{i_boundary}以及冰流速u_{_boundary}观测数据；

(2)控制方程

水动力方程包含的未知数为Z和Q，冰动力方程求解的未知数包括N、t_i、u，水动力与冰动力控制方程用守恒矢量型式来表示：

其中：

式中，A为水流断面面积；Q为流量；Z为水位；g为重力加速度；p_b、p_i分别为河床和岸冰湿周；τ_b、τ_w分别为在河底和冰—水界面上的剪切应力；τ_i为岸冰的剪切应力；τ′_ax为沿河道风对水面的拖曳力；ρ_i为冰的质量密度；x为空间坐标；N为冰的密集度；t_i为冰厚；u为冰的流速；B为水面冰层宽度；E_m为由热力过程所引起的表层冰质量的变化率、包括降雪以及浮冰在水中的交换；τ_w为沿河水的应力；τ_ax为沿河风的应力；τ_B为河岸剪切力；Z为水位；σ_x为沿水平方向上冰的平均应力；

为水流方向上重力分量；E_a为由于动力过程和悬浮冰花交换所引起的冰面积密度的时间变率；R_a为由机械再分配引起的冰面积密度的时间变率；

(3)冰水耦合运动统一离散方法

采用中心格式的有限体积法对控制方程进行离散，设第i个单元的长度为△x_i，对水动力和冰动力统一型式的守恒方程进行长度积分可得：

式中，D_i表示i单元内D的平均值，利用Green定理可得离散方程如下：

式中，U_i为单元中心点的向量，为整个单元的平均值；F_i+1/2和F_i-1/2分别为i+1/2和i-1/2界面处的通量；

(4)界面通量计算方法

根据Godunov方法，每个单元中的变量近似为恒定值，界面通量可通过Riemann问题进行计算；采用HLL近似Riemann解来计算界面通量，该方法稳定并且计算简便，界面通量定义如下：

式中，S_L和S_R分别为左侧和右侧波速；

当S_L<0<S_R时，水和冰的界面通量统一表达式，

具体如下：

浅水水流波速及冰压缩波波速分别通过下式计算：

式中，V_L和V_R分别为界面左侧和右侧水流流速；

和

分别为界面左侧和右侧水力半径；u_mL和u_mR分别为界面左侧和右侧修正冰流速；s_L和s_R分别为界面左侧和右侧冰压缩波波速；

(5)冰水耦合条件下浅水波波速计算方法

浅水波波速是界面通量计算的关键变量，在以往的研究中，认为浅水波波速与断面平均水深h^1/2呈线性关系，浅水波波速与冰体中压缩波波速二者相互独立；实际上，在有无冰体存在以及不同冰速条件下，水流运动所受阻力完全不同，可直接影响浅水波波速；在固定冰盖条件下，根据Einstein阻力划分原理，按零应力面将冰盖下过水断面划分为冰盖区A_i和床面区A_b；本发明提出等效水深的定义：为了使冰盖下水流运动与畅流期水流运动具有相同的阻力形式，将冰盖区按照零应力面展开(见图1)，此时的断面平均水深称为等效水深；宽浅型河道，等效水深

近似等于断面综合水力半径R_c，表达式如下：

式中，B_i为A_i对应的宽度；B_b为A_b对应的宽度；P_i和P_b分别为冰盖区湿周和床面区湿周；

浅水波波速计算表达式如下：

综合水力半径R_c计算表达式如下：

式中，α＝P_b/P_i，在流冰条件下，A_i、R_i以及P_i均会随冰速u变化而变化，冰盖区A_i和床面区A_b之间仍然存在零应力面，因此等效水深的概念仍然适用，当V_w<＝u时，A_i＝0，P_i＝0；当V_w>u时，A_i/A_b的表达式如下：

式中，

R′_i＝A_i/B，进一步推导P_i，表达式如下：

Q_i＝Q-Q_b  (24)

当u→V_w时，

是比

高阶的无穷小，因此P_i→0，并且α→∞，进而可得R_c＝R_b，这表明，当冰速接近断面平均流速时，A_i、P_i都趋于0，流冰对水流的阻力也趋于0，计算表达式与事实相符，在求得A_i和P_i后，进而可求得R_i及α，根据式(22)计算综合水力半径R_c，根据式(21)计算浅水波波速；

(6)高精度格式

HLL近似Riemann解可获得一阶精度数值解，对每个单元基本变量进行分段线性重构后，可获得空间二阶精度的数值解；

为了避免数值震荡，在估计坡度

时必须使用限制器，这里采用以稳健性著称的minmod限制器，2个变量，当其符号相同时取其最小值，否则取0，因此i单元的坡度表达式如下：

时间步上引入Hancock预测—校正方法，使数值解的精度整体提高到二阶，表达式如下：

通过式(30)的迭代可获得各变量的各个时间步和空间步的计算结果。

源项中，水面坡度的二阶精度表达式如下：

式中，

源项中的其它项采用显格式处理。

下边通过具体实例来对本发明作进一步说明。

洛古河水文站位于黑龙江干流额尔古纳河与石勒喀河汇合口下游约10km处，因额尔古纳河和石勒喀河均为南北走向，春季多形成“倒开江”，大量流冰进入黑龙江干流，受河道形态等因素的影响，在洛古河水文站下游约6.6km(航标880)附近易形成河床阻塞型冰坝。选取了2009年比较典型的冰凌洪水过程作为实例，对本发明开发的冰水耦合运动的高精度格式模拟方法进行验证。

2009年4月，黑龙江上游出现了1895年以来的最早开江和比较突出的倒开江。4月13日15时洛古河水文站下游6km处形成冰坝，致使水位迅速上涨，冰坝不断发展，13日21:07水位上涨到98.07m(黄海高程为307.85m)，冰坝溃毁。

本模型选取洛古河水文站断面至下游约6.6km河段作为计算范围，布设18个计算断面，具体见图2。上游边界条件采用洛古河水文站实测流量、冰密集度、冰厚、冰速成果，见表1；在1#断面附近形成河床阻塞型冰坝，因此下游边界条件冰速取0，冰密集度取1.0。

表1洛古河水文站2009年实测冰流量成果

一维冰水耦合运动的高精度格式模拟方法基本步骤如下：

①按照式(27)～(29)对每个单元的基本变量进行线性重构，得到界面处各变量值；

②按照式(12)～(19)对浅水波波速和冰压缩波波速进行计算，其中基于等效水深法的浅水波波速按照式(21)～(26)进行计算；

③按照式(6)～(11)对界面通量进行计算；

④按照式(30)计算下个时段的Z、Q、N、t_i、u。

采用洛古河水文站实测水位值与计算值进行对比，结果见图3。从图3中可见，洛古河水文站断面水位观测结果与计算结果趋势基本一致，4月13日21:00时，实测水位为307.74m，而计算值为307.63m，计算值与实测值相差较小。从实测流量来看，计算时段内流量变化幅度比较小，最大、最小流量仅相差81m³/s，而水位变幅相对比较大，并且升降频繁，主要受来冰量大小以及下游冰塞情况的影响。水面宽度的沿程变化见图4，冰厚的沿程变化见图5，冰速的沿程变化见图6，水位的沿程变化见图7。

从水面宽度的沿程变化上看，18#断面至14#断面水面宽度逐渐变小，14#断面至11#断面水面宽度逐渐变大，11#断面至9#断面水面宽度又逐渐变小，9#断面以下，水面宽度变化较小。水面宽度的沿程变化对冰厚的沿程变化影响比较明显，如图5所示，14#断面水面宽度较小，其上游断面冰厚较大，特别是18:00时，15#断面冰厚为该时刻沿程冰厚的最大值，达到了2.24m。另外，受下游河床阻塞冰坝边界条件影响，1#断面上游冰厚逐渐增加。随着时间推移，上游来冰量不断进入计算河段，14#断面卡塞冰量逐渐向下游运动，冰量逐渐向1#断面推进，21:00时3#断面处冰厚达到了2.80m。冰速方面，随着冰体向下游推移的过程中，冰速也在逐渐加大，21:00时17#断面冰速达到最大值，为2.26m/s。

沿程水位计算方面，当来冰卡塞堆积到一定厚度以后，其上游断面水面比降变小，其下游断面水面比降变大。如18:00时，14#断面以上水面比降较小，而14#断面以下水面比降较大，7#断面以上水面比降小，而7#断面以下水面比降较大，相比于12:00时，局部水面比降发生了比较明显的变化。冰对水流运动的影响主要体现在2个方面：①冰—水界面的切应力，与冰速和水流流速差值正相关，当冰速降低，断面流量下降，间接影响断面水位；②流冰条件下浅水波波速，当冰速降低，断面等效水深和水力半径减小、浅水波波速降低，对流量和水位均产生直接影响，特别对水位，使其界面通量减小，致使卡塞处上游壅水，水面比降减小。

进一步对比平均水深法和等效水深法计算浅水波波速下的河道沿程水位(图8)计算结果可见，采用等效水深法计算浅水波波速，3#断面上游河道水面比降明显小于平均水深法的计算结果，洛古河水文站断面水位计算结果较平均水深法高8cm，等效水深法水位计算结果更接近实测值。这表明采用等效水深法计算浅水波波速可明显提高水位预报精度，同时也表明了冰与水之间存在紧密耦合的关系，而非松散耦合。

经本实例验证，本发明提出的冰水耦合运动高精度格式模拟方法对流冰期河道水位和冰厚有较高的预报精度。

Claims (1)

1.一维冰水耦合运动的高精度格式模拟方法，其特征在于：包括下列步骤：

(1)边界条件

获取计算河段河道横断面数据，上边界流量Z_{_boundary}和水位Q_{_boundary}观测数据，上边界冰的密集度N_{_boundary}、冰厚t_{i_boundary}以及冰流速u_{_boundary}观测数据；

(2)控制方程

水动力方程包含的未知数为Z和Q，冰动力方程求解的未知数包括N、t_i、u，水动力与冰动力控制方程用守恒矢量型式来表示：

其中：

式中，A为水流断面面积；Q为流量；Z为水位；g为重力加速度；p_b、p_i分别为河床和岸冰湿周；τ_b、τ_w分别为在河底和冰—水界面上的剪切应力；τ_i为岸冰的剪切应力；τ′_ax为沿河道风对水面的拖曳力；ρ_i为冰的质量密度；x为空间坐标；N为冰的密集度；t_i为冰厚；u为冰的流速；B为水面冰层宽度；E_m为由热力过程所引起的表层冰质量的变化率、包括降雪以及浮冰在水中的交换；τ_w为沿河水的应力；τ_ax为沿河风的应力；τ_B为河岸剪切力；Z为水位；σ_x为沿水平方向上冰的平均应力；为水流方向上重力分量；E_a为由于动力过程和悬浮冰花交换所引起的冰面积密度的时间变率；R_a为由机械再分配引起的冰面积密度的时间变率；

(3)冰水耦合运动统一离散方法

采用中心格式的有限体积法对控制方程进行离散，设第i个单元的长度为△x_i，对水动力和冰动力统一型式的守恒方程进行长度积分可得：

式中，D_i表示i单元内D的平均值，利用Green定理可得离散方程如下：

式中，U_i为单元中心点的向量，为整个单元的平均值；F_i+1/2和F_i-1/2分别为i+1/2和i-1/2界面处的通量；

(4)界面通量计算方法

根据Godunov方法，每个单元中的变量近似为恒定值，界面通量可通过Riemann问题进行计算；采用HLL近似Riemann解来计算界面通量，该方法稳定并且计算简便，界面通量定义如下：

式中，S_L和S_R分别为左侧和右侧波速；

当S_L<0<S_R时，水和冰的界面通量统一表达式，具体如下：

浅水水流波速及冰压缩波波速分别通过下式计算：

式中，V_L和V_R分别为界面左侧和右侧水流流速；和分别为界面左侧和右侧水力半径；u_mL和u_mR分别为界面左侧和右侧修正冰流速；s_L和s_R分别为界面左侧和右侧冰压缩波波速；

(5)冰水耦合条件下浅水波波速计算方法

浅水波波速是界面通量计算的关键变量，在有无冰体存在以及不同冰速条件下，水流运动所受阻力完全不同，可直接影响浅水波波速；在固定冰盖条件下，根据Einstein阻力划分原理，按零应力面将冰盖下过水断面划分为冰盖区A_i和床面区A_b；提出等效水深的定义：为了使冰盖下水流运动与畅流期水流运动具有相同的阻力形式，将冰盖区按照零应力面展开，此时的断面平均水深称为等效水深；宽浅型河道，等效水深近似等于断面综合水力半径R_c，表达式如下：

式中，B_i为A_i对应的宽度；B_b为A_b对应的宽度；P_i和P_b分别为冰盖区湿周和床面区湿周；

浅水波波速计算表达式如下：

综合水力半径R_c计算表达式如下：

式中，α＝P_b/P_i，在流冰条件下，A_i、R_i以及P_i均会随冰速u变化而变化，冰盖区A_i和床面区A_b之间仍然存在零应力面，因此等效水深的概念仍然适用，当V_w<＝u时，A_i＝0，P_i＝0；当V_w>u时，A_i/A_b的表达式如下：

式中，R′_i＝A_i/B，进一步推导P_i，表达式如下：

Q_i＝Q-Q_b           (24)

当u→V_w时，是比高阶的无穷小，因此P_i→0，并且α→∞，进而可得R_c＝R_b，这表明，当冰速接近断面平均流速时，A_i、P_i都趋于0，流冰对水流的阻力也趋于0，计算表达式与事实相符，在求得A_i和P_i后，进而可求得R_i及α，根据式(22)计算综合水力半径R_c，根据式(21)计算浅水波波速；

(6)高精度格式

HLL近似Riemann解可获得一阶精度数值解，对每个单元基本变量进行分段线性重构后，可获得空间二阶精度的数值解；

为了避免数值震荡，在估计坡度时必须使用限制器，采用minmod限制器，2个变量，当其符号相同时取其最小值，否则取0，因此i单元的坡度表达式如下：

时间步上引入Hancock预测—校正方法，使数值解的精度整体提高到二阶，表达式如下：

通过式(30)的迭代可获得各变量的各个时间步和空间步的计算结果；

源项中，水面坡度的二阶精度表达式如下：

式中，源项中的其它项采用显格式处理。

Images