CN113536643B - 一种基于数字孪生的长河段滩槽演变预测方法及系统 - Google Patents

Abstract

本发明涉及了一种基于数字孪生的长河段滩槽演变预测方法及系统。本发明将二维水沙数学模型和正态物理模型结合，利用二维水沙数学模型模拟整个长河段的演变情况，利用正态物理模型对局部桥区进行模拟，获得高精度的局部桥区的计算结果，替换二维水沙数学模型模拟的该部分的结果，以克服二维水沙数学模型对局部桥区模拟的精度低的技术缺陷，提高桥区河段泥沙输移规律的模拟精度。

Description

一种基于数字孪生的长河段滩槽演变预测方法及系统

技术领域

本发明涉及地质监测技术领域，特别是涉及一种基于数字孪生的长河段滩槽演变预测方法及系统。

背景技术

近几十年来，人类活动的强烈干扰改变了河流自然输沙过程、输沙量和径流的年内分配，滩槽原有的形态及及其演变趋势必然发生自适应调整，长江中下游原有的河势演变规律将被打破，部分河段深泓摆动频繁，河势演变更趋复杂，滩槽的冲淤特征与冲淤量、床面形态等发生改变。

随着长江干线过江通道建设数量逐渐增多，近10年建设的桥梁主桥墩均涉水，并处于活动的边心滩区域。随着桥梁建成时间的增长，桥墩周围形成局部冲刷坑，改变洲滩形态，而河道形态调整又会对水沙输移特性产生影响，使得未守护的洲滩演变规律和航道条件变化趋势更为复杂。

处于冲淤平衡或接近平衡的自然河流，受到了建筑物的约束和控制，改变了水流的方向或流速的大小，或限制了原有浅滩的变形，则原有的平衡状态被破坏，将引起有限范围内的局部冲刷。如闸墩附近，枢纽下游，丁坝坝头以及堤防工程的根部，护滩带周围等，都是容易发生局部冲刷的地方。长河段二维水沙数学模型在系统治理和上下游洲滩联动规律模拟方面具有优势，但桥墩或整治建筑物周围具有明显的三维水流特性，二维模型不能真实反映桥墩周围极限冲刷深度，由于桥墩周围床面形态模拟与实际情况不一致，进而影响桥区河段泥沙输移规律的模拟精度。

发明内容

本发明的目的是提供一种基于数字孪生的长河段滩槽演变预测方法及系统，以提高桥区河段泥沙输移规律的模拟精度。

为实现上述目的，本发明提供了如下方案：

本发明提供一种基于数字孪生的长河段滩槽演变预测方法，所述方法包括如下步骤：

构建长河段的二维水沙数学模型和长河段内的局部桥区的局部正态物理模型；

根据第i年的年初的水文数据，采用离散方法求解所述二维水沙数学模型，获得第i年年末的水流运动数据和河床冲淤厚度数据，所述水流运动数据包括x方向流速、y方向流速、悬移质含沙量和水位；

基于所述第i年年末的水流运动数据，采用局部正态物理模型进行局部模拟实验，获得第i年的局部桥区的模拟形变量；

利用第i年的局部桥区的模拟形变量替换河床冲淤厚度数据中的局部桥区的河床冲淤厚度，得到第i年年末的替换后的河床冲淤厚度数据；

将第i年年末的水流运动数据和第i年年末的替换后的河床冲淤厚度数据，作为第i+1年年初的水文数据，令i的数值增加1，返回步骤“根据第i年的年初的水文数据，采用离散方法求解所述二维水沙数学模型，获得第i年年末的水流运动数据和河床冲淤厚度数据”，直到到达预测年份。

可选的，所述离散方法包括有限体积法、有限差分法和/或有限元法。

可选的，所述二维水沙数学模型包括：

水流连续方程：

X方向和Y方向水流运动方程：

悬沙不平衡输运方程：

推移质不平衡输移方程：

其中，ζ为水位，t为时间，h为静水深，u、v为流速矢量

沿X、Y方向的分量，g为重力加速度，ν_e为水流紊动黏性系数，c为谢才系数，f为科氏参量，

ω为地球自转角速度，

为地理纬度，ρ为水密度；S为单位水体垂线平均含沙量，

s＝ρ_sc，c为单位水体体积浓度，ρ_s为含沙密度；ν_t和σ_S分别为第一Schmidt数和第二Schmidt数；ω_S为泥沙沉速，下标i表示第i组泥沙，Φ_s为泥沙扩散通量，N_i、

分别为推移质输沙量和推移质输沙能力折算成相应水深的泥沙浓度，β_i为推移质泥沙恢复饱和系数，ω_si为第i组泥沙的沉速，H为水深，H＝h+ζ。

可选的，所述二维水沙数学模型的边界条件包括：

上游进口条件：

下游出口条件：ζ_j2＝ζ_out；

固壁条件：U_j3＝V_j3＝0，

收敛控制条件：

其中，Q_in(t)表示上游来流量，ξ为水位，U_j1和V_j1分别表示上游进口控制点j1在x方向和y方向的速度，h_j1表示上游进口控制点j1处的静水深，dy_j1为离散网格间距，α_j1为动能修正系数，其值大小取决于过水断面上流速分布情况，流速分布越均匀，其值越接近于1，一般取1.05～1.1，ζ_j2为下手出口控制点j2的水位，ζ_out表示给定出口水位，U_j3和V_j3分别表示固壁的控制点j3在x方向和y方向的速度，S_j3为固壁的控制点j3处的含沙量，N_j3为固壁的控制点j3处的底沙，

为梯度值，b_max为最大质量源，Q_j表示控制点j的断面流量，

表示第n次迭代运算时控制点j和控制点j’之间的x方向的速度差，

表示第n次迭代运算时控制点j和控制点j’之间的水位差，

表示第n+1次迭代运算时控制点j和控制点j’之间的x方向的速度差，

表示第n+1次迭代运算时控制点j和控制点j’之间的水位差。

可选的，所述根据第i年的年初的水文数据，采用离散方法求解所述二维水沙数学模型，获得第i年年末的水流运动数据和河床冲淤厚度数据，具体包括：

根据第i年的年初的水文数据，采用并行运算的方式，采用离散方法求解所述二维水沙数学模型，获得第i年年末的水流运动数据和河床冲淤厚度数据。

一种基于数字孪生的长河段滩槽演变预测系统，所述系统包括：

模型构建模块，用于构建长河段的二维水沙数学模型和长河段内的局部桥区的局部正态物理模型；

数学模型计算模块，用于根据第i年的年初的水文数据，采用离散方法求解所述二维水沙数学模型，获得第i年年末的水流运动数据和河床冲淤厚度数据，所述水流运动数据包括x方向流速、y方向流速、悬移质含沙量和水位；

物理模型计算模型，用于基于所述第i年年末的水流运动数据，采用局部正态物理模型进行局部模拟实验，获得第i年的局部桥区的模拟形变量；

局部数据替换模块，用于利用第i年的局部桥区的模拟形变量替换河床冲淤厚度数据中的局部桥区的河床冲淤厚度，得到第i年年末的替换后的河床冲淤厚度数据；

返回模块，用于将第i年年末的水流运动数据和第i年年末的替换后的河床冲淤厚度数据，作为第i+1年年初的水文数据，令i的数值增加1，返回步骤“根据第i年的年初的水文数据，采用离散方法求解所述二维水沙数学模型，获得第i年年末的水流运动数据和河床冲淤厚度数据”，直到到达预测年份。

可选的，所述离散方法包括有限体积法、有限差分法和/或有限元法。

可选的，所述二维水沙数学模型包括：

水流连续方程：

X方向和Y方向水流运动方程：

悬沙不平衡输运方程：

推移质不平衡输移方程：

其中，ζ为水位，t为时间，h为静水深，u、v为流速矢量

沿X、Y方向的分量，g为重力加速度，ν_e为水流紊动黏性系数，c为谢才系数，f为科氏参量，

ω为地球自转角速度，

为地理纬度，ρ为水密度；S为单位水体垂线平均含沙量，

s＝ρ_sc，c为单位水体体积浓度，ρ_s为含沙密度；ν_t和σ_S分别为第一Schmidt数和第二Schmidt数；ω_S为泥沙沉速，下标i表示第i组泥沙，Φ_s为泥沙扩散通量，N_i、

分别为推移质输沙量和推移质输沙能力折算成相应水深的泥沙浓度，β_i为推移质泥沙恢复饱和系数，ω_si为第i组泥沙的沉速，H为水深，H＝h+ζ。

可选的，所述二维水沙数学模型的边界条件包括：

上游进口条件：

下游出口条件：ζ_j2＝ζ_out；

固壁条件：U_j3＝V_j3＝0，

收敛控制条件：

其中，Q_in(t)表示上游来流量，ξ为水位，U_j1和V_j1分别表示上游进口控制点j1在x方向和y方向的速度，h_j1表示上游进口控制点j1处的静水深，dy_j1为离散网格间距，α_j1为动能修正系数，其值大小取决于过水断面上流速分布情况，流速分布越均匀，其值越接近于1，一般取1.05～1.1，ζ_j2为下手出口控制点j2的水位，ζ_out表示给定出口水位，U_j3和V_j3分别表示固壁的控制点j3在x方向和y方向的速度，S_j3为固壁的控制点j3处的含沙量，N_j3为固壁的控制点j3处的底沙，

为梯度值，b_max为最大质量源，Q_j表示控制点j的断面流量，

表示第n次迭代运算时控制点j和控制点j’之间的x方向的速度差，

表示第n次迭代运算时控制点j和控制点j’之间的水位差，

表示第n+1次迭代运算时控制点j和控制点j’之间的x方向的速度差，

表示第n+1次迭代运算时控制点j和控制点j’之间的水位差。

可选的，所述数学模型计算模块，具体包括：

数学模型计算子模块，用于根据第i年的年初的水文数据，采用并行运算的方式，采用离散方法求解所述二维水沙数学模型，获得第i年年末的水流运动数据和河床冲淤厚度数据。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明公开了一种基于数字孪生的长河段滩槽演变预测方法，所述方法包括如下步骤：构建长河段的二维水沙数学模型和长河段内的局部桥区的局部正态物理模型；根据第i年的年初的水文数据，采用离散方法求解所述二维水沙数学模型，获得第i年年末的水流运动数据和河床冲淤厚度数据，基于所述第i年年末的水流运动数据，采用局部正态物理模型进行局部模拟实验，获得第i年的局部桥区的模拟形变量；利用第i年的局部桥区的模拟形变量替换河床冲淤厚度数据中的局部桥区的河床冲淤厚度，得到第i年年末的替换后的河床冲淤厚度数据。本发明将二维水沙数学模型和正态物理模型结合，利用二维水沙数学模型模拟整个长河段的演变情况，利用正态物理模型对局部桥区进行模拟，获得高精度的局部桥区的计算结果，替换二维水沙数学模型模拟的该部分的结果，以克服二维水沙数学模型对局部桥区模拟的精度低的技术缺陷，提高桥区河段泥沙输移规律的模拟精度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明提供的一种基于数字孪生的长河段滩槽演变预测方法的流程图；

图2为本发明提供的数学模型、物理模型耦合求解流程示意图；

图3为本发明提供的泥沙坍塌调整示意图；

图4为本发明提供的长江下游河段耦合模型模拟范围示意图；

图5为本发明提供的物理模型水槽平面布置图；

图6为本发明提供的物理模型水槽纵剖面布置示意图；

图7为本发明提供的局部正态物理模型桥墩局部冲刷后地形图，图7a为系列年第1年末的地形图，图7b为系列年第2年末的地形图，图7c为系列年第3年末的地形图，图7d为系列年第5年末的地形图，7e为系列年第10年末的地形图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的目的在于提出一种基于数字孪生的长河段滩槽演变预测方法及系统，以提供一种长河段复杂边界条件下基于数字孪生的数学模型和物理模型耦合运行方法，实现兼顾局部河床冲淤调整和上下游之间联动规律的模拟。可以用于内河长河段系统治理效果模拟，也可用于桥区河段“桥梁工程+航道工程+水沙调控”综合作用下滩槽演变规律预测。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

河工模型试验一般采用变态模型，其水流运动在垂直方向是不相似的。而局部冲刷坑附近的水流运动三维特性明显，局部的旋涡与环流等水流结构是决定冲刷坑形态和大小的最重要的因素之一，因而较少利用变态模型研究纵剖面的急剧变化，多采用正态模型研究局部冲刷问题。

长江中下游河床冲淤多变，桥梁工程和航道整治工程的实施，形成了新的河势约束边界，进一步影响桥区河段浅滩形态及航道尺度。桥梁工程与航道工程耦合作用下形成了复杂的控制边界条件，江中桥墩的存在一定程度上改变了泥沙输移路径，桥墩对于河床冲淤的影响长期存在，特别是周围局部冲刷坑形成后，桥区河段边心滩周期性变化规律将受到影响。

数字孪生“Digital Twin”是由美国NASA首先提出的，一般认为包含以下几层含义：

(1)在数字空间内，使用高精度的数学模型来描述和模拟现实世界中的事物。

(2)将现实世界中采集的真实信息反映到数学模型，使之随现实进行更新。

(3)在数字空间内，使用模型和信息进行预测性的仿真分析和可视化。

本发明采用数学模型和物理模型，利用数字孪生技术进行数据融合、数据衍生，进而实现长河段高精度的复杂边界条件下滩槽演变规律模拟预测。

本发明是通过以下技术方案实现的：

河工模型试验一般采用变态模型，其水流运动在垂直方向是不相似的。而局部冲刷坑附近的水流运动三维特性明显，局部的旋涡与环流等水流结构是决定冲刷坑形态和大小的最重要的因素之一，因而较少利用变态模型研究纵剖面的急剧变化，多采用正态模型研究局部冲刷问题。

长江中下游河床冲淤多变，桥梁工程和航道整治工程的实施，形成了新的河势约束边界，进一步影响桥区河段浅滩形态及航道尺度。桥梁工程与航道工程耦合作用下形成了复杂的控制边界条件，江中桥墩的存在一定程度上改变了泥沙输移路径，桥墩对于河床冲淤的影响长期存在，特别是周围局部冲刷坑形成后，桥区河段边心滩周期性变化规律将受到影响。

本发明提供一种基于数字孪生的数学模型和物理模型耦合运行方法。在传统二维水沙数学模型基础上，优化泥沙模型中河床坡度对推移质输沙的影响，基于OpenACC对二维水沙模型进行了GPU并行计算开发，提高二维水沙数学模型的计算效率和稳定性。通过局部河段正态物理模型试验模拟桥墩周围局部冲刷坑形态，如图7所示，将物理模型试验得出的桥墩周围床面形态赋予数学模型，对数学模型中桥墩周围床面形态进行调整，再通过长河段二维水沙数学模型模拟桥墩局部冲刷坑形成后不同水动力条件和泥沙条件下滩槽间水动力变化和底沙输移规律。

物理模型每个试验年末桥墩周围河床变形Δz_物数据均作为长河段二维水沙数学模型的边界输入条件。具体实现方法如下：当长河段水沙数学模型计算至系列年第一年末时，将局部正态物理模型试验第一年末对桥墩周围河床形态的试验结果Δz_{物(第一年末)}替换数学模型桥墩周围河床冲淤厚度η_{(第一年末)}，以此作为第二年计算的起始地形；当数学模型计算至系列年第二年末时，将局部正态物理模型试验第二年末对桥墩周围河床形态的试验结果Δz_{物(第二年末)}替换数学模型桥墩周围河床冲淤厚度η_{(第二年末)}，以此作为第三年计算的起始地形；以此类推，直至系列年末整个计算过程结束。

如图1和2所示，本发明提供本发明提供一种基于数字孪生的长河段滩槽演变预测方法，所述方法包括如下步骤：

步骤101，构建长河段的二维水沙数学模型和长河段内的局部桥区的局部正态物理模型。以长江下游某河段为例，如图4所示为长江下游某河段耦合模型模拟范围示意图。

步骤102，根据第i年的年初的水文数据，采用离散方法求解所述二维水沙数学模型，获得第i年年末的水流运动数据和河床冲淤厚度数据，所述水流运动数据包括x方向流速、y方向流速、悬移质含沙量和水位。

步骤102具体包括：

为了考虑河段内上下游洲滩联动变化，首先，建立长河段二维水沙数学模型，采用式(1)、(2)、(3)、(4)，分别离散后，求得U、V、S、ζ，离散方法可采用有限体积法、有限差分法和有限元法等，数学模型由压力校正法进行压力(水位)～流速耦合求解，得出正确的流速场和水位值、以及含沙量场。

二维水流泥沙数学模型控制方程和初始边界条件如下：

平面二维水沙数学模型控制方程包括水流连续方程、X方向水流运动方程和Y方向水流运动方程、悬移质不平衡输移方程、推移质不平衡输移方程和河床变形方程。通过这些方程进行联立求解，组成方程组，计算得到流速U、V，水位ζ、含沙量S，底沙N，以及河床变形η。

1)、水流连续方程：

2)、X方向和Y方向水流运动方程：

式中：ζ——水位(m)；t——时间(s)；h——静水深(m)；u、v——流速矢量

沿X、Y方向的分量(m/s)；g——重力加速度(m/s²)；ν_e——水流紊动黏性系数(m²/s)；c——谢才系数(m^1/2/s)；f——科氏参量(s^-1)，

ω为地球自转角速度，

为地理纬度；ρ——水密度(kg/m³)。

3)、悬沙不平衡输运方程

由

沿水深积分，并假定由流速和含沙量沿垂线分布不均匀在积分时产生的修正系数：

引入冲淤平衡时的挟沙能力S^*，得：

式中：S为单位水体垂线平均含沙量，

s＝ρ_sc，c为单位水体体积浓度；ν_t＝ν_mt；σ_S＝σ_c为Schmidt数；ω_S为泥沙沉速，下标i表示非均匀泥沙分组情况。

式(4)中，

对于水面z＝ζ泥沙扩散通量为零边界条件：

对于底部z＝-h泥沙扩散通量：

一般认为悬沙粒径很细时，不论泥沙沿水深分布是否处于平衡状态，含沙量沿水深变化不大，上式表示为：Φ_s＝αω_s(S^*-S)。

式中α＝α^*P_r为系数，此表达式在泥沙输移数模计算中得到广泛运用。关于表达式中的系数α，韩其为在研究悬沙二维扩散方程的边界条件时，定义为恢复饱和系数，并如下关系：

在数学模型计算中，垂线恢复饱和系数α取值范围为0.25～1.0，淤积状态取α＝0.25；冲刷状态取α＝1.0。

分组挟沙力S_n ^*计算按照窦国仁、赵士清的模式，将非均匀沙按其粒径大小分成N₀组，S_n表示n组粒径得含沙量，P_n表示此粒径在悬沙总含沙量S中所占的比值：

总挟沙力：

挟沙力级配：

分组挟沙力：S_n ^*＝P_n ^*S_*，

式中：0<α<1，ω_n为第n组粒径的沉速，ω_m为非均匀平均沉速。

4)、推移质不平衡输移方程

根据推移质不平衡非均匀输沙原理，通过推移质水深折算推导出底沙不平衡输沙输移方程：

对于非均匀沙，推移质不平衡输移方程采用如下形式：

式中，N_i、

分别为推移质输沙量和推移质输沙能力折算成相应水深的泥沙浓度，β_i为推移质泥沙恢复饱和系数，下标i表示第i组粒径泥沙对应的变量。对于非均匀沙第i组粒径泥沙输沙率。

推移质输沙率的计算公式众多，目前较常用的有：van Rijn公式、窦国仁公式、岗恰洛夫公式等，对于非均匀沙第i组粒径泥沙输沙率考虑了隐蔽系数ηi：

总的输沙能力；

P_bi为第i组粒径泥沙所占的百分比；

为第i组粒径泥沙推移质输沙率。

窦国仁推移质输沙率公式为：

式中：

V_ki为第i组粒径泥沙临界启动流速：

式中：d_i为第i组泥沙粒径，γ_o为稳定干容重，γ_o＝1650kg/m³；ε_k为粘结力参数(天然沙ε_k＝2.56cm³/s²)；σ为薄膜水厚度，σ＝0.21×10^-4cm，γ′_o为床面泥沙干容重，对于细沙γ′_o＝γ_o，C为谢才系数：

Δ为床面糙度：

K为系数，对于沙质推移质k取0.01，这样方程(6)式中

可写成：

数学模型在计算开始前，要预先给定初始条件和边界条件，具体方法如下：

①、初始条件

给定初始条件时刻t＝0时，计算域内所有计算变量(U、V、ζ、Si、Ni)的初值，给出悬沙级配和分段床沙级配。

②、上、下游控制边界条件

上游进口条件：给定上游来流Q_in、沙量S_in和进口推移质输沙率q_b，进口各点流速

V_i＝0。进口各控制点流速由下式迭代算出：

式中：U_j，h_j为进口计算网格点沿y方向流速和水深，d_yj为离散网格间距，V_j＝0。

下游出口条件：给定水位ζ_out。

③、固壁条件

流速采用非滑移边界条件，其边壁流速给定为零，即U＝V＝0；对于含沙量S_i，底沙N_i计算中采用法向梯度为零条件：

④、收敛控制条件

控制连续方程最大质量源b_max和通过各断面流量Q_j，

流速：

水位：

离散方法可采用有限体积法、有限差分法和有限元法等，得出正确的流速场和水位值，以及含沙量场。

下面以有限体积法离散为例，说明水流运动方程求解过程。具体离散求解过程如下，在任意一个给定的三角形控制体内对方程(1)进行积分，并利用格林公式，得

式中：A_i为单元i的面积。

写成离散形式得

或

类似地，对方程(2)和(3)进行积分，采用时间向前差分，得到

式中：

为第i个单元的第j边上水平数值通量，采用基于黎曼解的Osher格式计算；A_i为第i个单元的面积，l_i,j为第i个单元的第j条边的边长；D_hU、D_hV分别为x、y方向水平紊动扩散项，其积分形式可表示为

对以上方程进行整理可得

通过上式，可求得水位ζ、x方向流速U、y方向流速V，通过类似离散方法，离散方程(4)和(6)，可求得S、N。

在长河段二维水沙数学模型计算过程中，为了进一步提高桥墩周围泥沙输移的模拟精度，在局部地形替换的基础上，针对局部冲刷床面坡度较大的特殊情况，分析床面坡度对推移质输沙的影响，对已有的考虑坡度影响的推移质输沙率修正方法进行改进，在求解河床变形的过程中，除了考虑冲淤造成的床面变化外，还将对床面发生坍塌调整的过程进行处理。在斜坡床面上，除了作用在泥沙颗粒上的水流剪切力外，重力的切向分力也是推移质输运的一个动力因素。为了考虑坍塌的影响，在数值模型中，加入河床发生坍塌调整的模块。在根据泥沙模型计算完河床变形后，对床面网格进行扫描。每个时间步网格节点高程根据地貌模型调整后，都需要判断每个单元平面外法向方向与Z轴正向的夹角大小，当发现相邻角点间的倾角大于泥沙休止角时，要进行局部调整。

如图3所示，为了考虑坍塌的影响，在数值模型中，对公式(6)中推移质输沙率计算结果进行局部修正。具体操作方法如下：根据泥沙模型计算完河床变形后，对床面网格进行扫描。每个时间步网格节点高程根据地貌模型调整后，都需要判断每个单元平面外法向方向与Z轴正向的夹角大小，当发现相邻角点间的倾角大于泥沙休止角时，要进行局部调整，调整模式如图3所示，在扫描时发现角点A(x_A，y_A，z_A)和角点B(x_B，y_B，z_B)之间倾角大于休止角

需要将较高的B点降至B′点，A点升至A′点，使AB之间的倾角降到

表示如下：

经过对床面的一次扫描调整后可能会出现新的陡坡，因此需要对整个床面进行反复的扫描检查直至床面上所有相邻角点间的倾角均在休止角范围内。

在长河段二维水沙数学模型计算过程中，为了提高长河段二维数学模型的计算效率，对二维水沙模型进行GPU并行计算改造。具体方法如下：

将串行程序进行并行化改造前，首先需要对原串行程序进行分析，找出程序中的热点，也就是程序中计算最集中，花费时间最长的地方；然后通过对热点区域的数据相关性，算法的可并行性进行分析，将热点区域进行并行化，以达到最好的并行化效果。

从一个串行程序得到一个并行程序的工作一般由下面几个步骤构成：

①对串行程序进行分析，找出其内在的可并行性，提出适合的并行算法，将计算分解为多个子任务；

②将任务分配给各个进程，一般一个处理器启动一个进程；

③在各个进程之间协调数据的访问、通信及同步。

串行程序采用非耦合解的形式，水流计算过程和泥沙计算过程分别进行，即水流恒定后，再将u、v、h等水力要素代入泥沙方程中进行求解，因此，可以单独考虑水流和泥沙的并行实现过程。

通过IntelVtune工具对串行程序运行60s的采样数据表明，71.56％左右的CPU时间被求解代数方程组子程序占用，21.6％的CPU时间被求解方程组系数子程序占用，因此首先考虑这部分的并行化处理。

下面以水流子程序为例说明并行化设计过程。

通过OpenACC提供的内核并行化导语，将这些热点进行改造，结合NVIDIA公司提供的OpenACC Toolkit工具，可以对并行化改造后的程序计算效率进行测试。

本程序对绕节点、单元和边的执行循环进行内核并行化改造，使用的计算构件是parallel。当程序执行遇到一个加速器paralle构件时，程序就创建一个或多个gang来执行这个加速器的并行区域。gang的数量、每个gang中worker的数量和每个worker中vector的数量在变个并行区域的存续期内保持为常数，其数量设置在paralle导语的子语中进行给定。

程序实现的伪代码为：

计算效率测试以长河段二维水沙数学模型为例，GPU计算采用的显卡为AMDRadeon R54302GB，CPU为intel Core i7-8700@3.2GHz六核。对比了不同网格数量条件下CPU串行计算和GPU并行计算的时间和效率，从中可以看出，采用OpennACC加速后的程序计算效率明显高于CPU串行计算，在网格数量较少时加速比相对较低，约为8倍；但随着网格节点数量的增加，加速比明显增大，可以达到20倍以上。

步骤103，基于所述第i年年末的水流运动数据，采用局部正态物理模型进行局部模拟实验，获得第i年的局部桥区的模拟形变量；

针对长河段内桥区局部开展正态动床物理模型试验，其中试验边界条件，如模型进口U、V、S，模型出口ζ，均由上述长河段数学模型给定。通过试验，采集记录典型年或系列年试验过程中每个试验年末桥墩周围河床变形Δz_物。典型水文年和系列水文年的选取应根据研究问题的需要从不利的角度选取不同的水沙年组合作为系列水文年，系列水文年宜包含丰水丰沙、中水中沙和小水少沙等特征水文年。物理模型设计、制作和试验方法如下：

根据场地条件、实验室供水能力，拟选用宽度5m的水槽进行研究，初步选择水槽概化模型比尺为1：60的正态。采用正态模型研究河床局部冲刷，模型设计应满足水流运动、泥沙运动和河床变形等相似条件。

则各模型比尺如下：

平面比尺：α_L＝60

垂直比尺：α_H＝60

流速比尺：α_V＝α_H ^1/2＝7.746

糙率比尺：

重量相似比尺：

试验水槽为长40m、宽5m的矩形水槽，模型的进出口控制系统采用自动控制系统。其中进口流量控制选用电磁流量计控制模型入流；模型出口水位选用远程控制的自动翻板门进行调节。模型水位由自动水位仪读取，读取精度可达到0.1mm。断面流速观测采用小威龙测速仪或无线旋桨流速仪等。如图5和6所示，水槽的上游设格墙以调整平顺水流，尾部设置沉沙池，水槽中间20m范围为动床试验段。沿程设置12把水尺观测水位，其中12#水尺控制尾门水位，水尺均设置在水槽两侧，3#～10#用以观测桥墩附近水位变化。

动床模型水沙系统由独立的清水和加沙系统组成，清水系统由水泵、电机、管道、平水塔、量水堰组成，加沙系统主要由3m3的圆形搅拌池和1.5kw的三相变速电机及输沙管道等组成。输沙管上安装的闸阀控制加入模型中的含沙水流的流量，模型试验时根据长河段二维数学模型提供的沙量S和推移质输沙率q_b，自动开启一定开度加沙孔，同时要求含沙水流在模型前充分混合后进入模型。

如图5和6所示，物理模型采用研究河段最新实测地形制模，动床模型地形断面采用0.70mm厚的镀锌铁皮制作(可重复使用)。动床范围平均铺沙厚度约15cm。

尾门控制水位取长河段二维水沙模型计算各级流量时水槽尾门处的水位ζ，水槽上游试验控制流速U、V由长河段二维水沙模型计算得到，模型进口沙量S和推移质输沙率q_b均由长河段二维水沙数学模型提供。

在试验过程中主要监测以下内容：

①沿程水位变化：观测水槽两侧的水位变化。

②进口断面的流速分布测量：每级流量在动床上游的定床部分(1#断面)进行断面流速分布测量，以了解进口断面的流速分布与流量闭合情况。

③冲刷过程中断面表面流速监测：在桥墩周围表面附近(5#～6#水尺断面)布置有断面测量表面流速。

④冲刷过程中断面冲刷形态监测：在桥墩周围附近(5#～6#水尺断面上、下游及两侧0.5m处)布置有断面监测冲刷地形的变化，如图7(图7a-图7e)为不同年份的局部正态物理模型桥墩局部冲刷后地形。

⑤每个典型年末冲刷后地形测量，得到Δz_物。

局部正态物理模型试验结束后，进行长河段二维水沙数学模型计算，模型计算采用的典型水文年与系列水文年与物理模型一致。采用式(8)、(9)，分别离散后求得悬移质冲淤引起的河床变形η_si和推移质冲淤引起的河床变形η_bi，再通过式(9)，求得河床总的冲淤厚度η。

由悬移质冲淤引起的河床变形方程为：

式中：η_si为第i组粒径悬移质泥沙引起的冲淤厚度。γ₀为床面泥沙干容重。

由推移质冲淤引起的河床变形方程为：

式中：η_bi为第i组粒径推移质泥沙引起的冲淤厚度。

河床总的冲淤厚度：

步骤104，利用第i年的局部桥区的模拟形变量替换河床冲淤厚度数据中的局部桥区的河床冲淤厚度，得到第i年年末的替换后的河床冲淤厚度数据；

数学模型中桥墩周围局部地形修正。采用局部正态物理模型每个试验年末桥墩周围河床变形Δz_物数据替换长河段二维水沙数学模型计算结果η。具体实现方法如下：

给定长河段二维水沙数学模型初始和边界条件后，开始计算，在计算过程中，每隔7～10天对河床地形进行一次更新，即用新时刻地形(η_t)替换上一时刻地形(η_t-1)，当计算至系列年第一年末时，得到第一年末的河床冲淤厚度η_{(第一年末)}，将步骤2中局部正态物理模型试验第一年末对桥墩周围河床形态的试验结果Δz_{物(第一年末)}替换步骤3中数学模型计算的桥墩周围河床冲淤厚度η_{(第一年末)}，完成第一年末桥墩周围河床地形替换；

第二年开始计算时，数学模型以替换后的第一年末地形作为初始地形继续计算，当数学模型计算至系列年第二年末时，将局部正态物理模型试验第二年末对桥墩周围河床形态的试验结果Δz_{物(第二年末)}替换数学模型桥墩周围河床冲淤厚度η_{(第二年末)}，完成第二年末桥墩周围河床地形替换；

以此类推，直至系列年末整个计算过程结束。最终得到系列年末，数学模型整个计算区域内河床冲淤变形η_{(系列年末)}。

步骤105，将第i年年末的水流运动数据和第i年年末的替换后的河床冲淤厚度数据，作为第i+1年年初的水文数据，令i的数值增加1，返回步骤“根据第i年的年初的水文数据，采用离散方法求解所述二维水沙数学模型，获得第i年年末的水流运动数据和河床冲淤厚度数据”，直到到达预测年份。

一种基于数字孪生的长河段滩槽演变预测系统，所述系统包括：

模型构建模块，用于构建长河段的二维水沙数学模型和长河段内的局部桥区的局部正态物理模型；

数学模型计算模块，用于根据第i年的年初的水文数据，采用离散方法求解所述二维水沙数学模型，获得第i年年末的水流运动数据和河床冲淤厚度数据，所述水流运动数据包括x方向流速、y方向流速、悬移质含沙量和水位。

所述数学模型计算模块，具体包括：

数学模型计算子模块，用于根据第i年的年初的水文数据，采用并行运算的方式，采用离散方法求解所述二维水沙数学模型，获得第i年年末的水流运动数据和河床冲淤厚度数据。

所述离散方法包括有限体积法、有限差分法和/或有限元法。

可选的，所述二维水沙数学模型包括：

水流连续方程：

X方向和Y方向水流运动方程：

悬沙不平衡输运方程：

推移质不平衡输移方程：

其中，ζ为水位，t为时间，h为静水深，u、v为流速矢量

沿X、Y方向的分量，g为重力加速度，ν_e为水流紊动黏性系数，c为谢才系数，f为科氏参量，

ω为地球自转角速度，

为地理纬度，ρ为水密度；S为单位水体垂线平均含沙量，

s＝ρ_sc，c为单位水体体积浓度，ρ_s为含沙密度；ν_t和σ_S分别为第一Schmidt数和第二Schmidt数；ω_S为泥沙沉速，下标i表示第i组泥沙，Φ_s为泥沙扩散通量，N_i、

分别为推移质输沙量和推移质输沙能力折算成相应水深的泥沙浓度，β_i为推移质泥沙恢复饱和系数，ω_si为第i组泥沙的沉速，H为水深，H＝h+ζ。

所述二维水沙数学模型的边界条件包括：

上游进口条件：

下游出口条件：ζ_j2＝ζ_out；

固壁条件：U_j3＝V_j3＝0，

收敛控制条件：

其中，Q_in(t)表示上游来流量，ξ为水位，U_j1和V_j1分别表示上游进口控制点j1在x方向和y方向的速度，h_j1表示上游进口控制点j1处的静水深，dy_j1为离散网格间距，α_j1为动能修正系数，其值大小取决于过水断面上流速分布情况，流速分布越均匀，其值越接近于1，一般取1.05～1.1，ζ_j2为下手出口控制点j2的水位，ζ_out表示给定出口水位，U_j3和V_j3分别表示固壁的控制点j3在x方向和y方向的速度，S_j3为固壁的控制点j3处的含沙量，N_j3为固壁的控制点j3处的底沙，

为梯度值，b_max为最大质量源，Q_j表示控制点j的断面流量，

表示第n次迭代运算时控制点j和控制点j’之间的x方向的速度差，

表示第n次迭代运算时控制点j和控制点j’之间的水位差，

表示第n+1次迭代运算时控制点j和控制点j’之间的x方向的速度差，

表示第n+1次迭代运算时控制点j和控制点j’之间的水位差。

物理模型计算模型，用于基于所述第i年年末的水流运动数据，采用局部正态物理模型进行局部模拟实验，获得第i年的局部桥区的模拟形变量；

局部数据替换模块，用于利用第i年的局部桥区的模拟形变量替换河床冲淤厚度数据中的局部桥区的河床冲淤厚度，得到第i年年末的替换后的河床冲淤厚度数据；

返回模块，用于将第i年年末的水流运动数据和第i年年末的替换后的河床冲淤厚度数据，作为第i+1年年初的水文数据，令i的数值增加1，返回步骤“根据第i年的年初的水文数据，采用离散方法求解所述二维水沙数学模型，获得第i年年末的水流运动数据和河床冲淤厚度数据”，直到到达预测年份。

本发明利用长河段二维水沙数学模型和局部正态物理模型的优势，取长补短，弥补了二维数学模型在局部冲刷模拟中的不足，提高了二维数学模型在局部冲刷模拟精度，解决了长河段复杂边界条件下河床变形耦合模拟技术，实现了兼顾局部河床冲淤调整和上下游之间联动规律的模拟。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种基于数字孪生的长河段滩槽演变预测方法，其特征在于，所述方法包括如下步骤：

构建长河段的二维水沙数学模型和长河段内的局部桥区的局部正态物理模型；

根据第i年的年初的水文数据，采用离散方法求解所述二维水沙数学模型，获得第i年年末的水流运动数据和河床冲淤厚度数据，所述水流运动数据包括x方向流速、y方向流速、悬移质含沙量和水位；

基于所述第i年年末的水流运动数据，采用局部正态物理模型进行局部模拟实验，获得第i年的局部桥区的模拟形变量；

利用第i年的局部桥区的模拟形变量替换河床冲淤厚度数据中的局部桥区的河床冲淤厚度，得到第i年年末的替换后的河床冲淤厚度数据；

将第i年年末的水流运动数据和第i年年末的替换后的河床冲淤厚度数据，作为第i+1年年初的水文数据，令i的数值增加1，返回步骤“根据第i年的年初的水文数据，采用离散方法求解所述二维水沙数学模型，获得第i年年末的水流运动数据和河床冲淤厚度数据”，直到到达预测年份。

2.根据权利要求1所述的基于数字孪生的长河段滩槽演变预测方法，其特征在于，所述离散方法包括有限体积法、有限差分法和/或有限元法。

3.根据权利要求1所述的基于数字孪生的长河段滩槽演变预测方法，其特征在于，所述二维水沙数学模型包括：

水流连续方程：

X方向和Y方向水流运动方程：

悬沙不平衡输运方程：

推移质不平衡输移方程：

其中，ζ为水位，t为时间，h为静水深，u、v为流速矢量

沿X、Y方向的分量，g为重力加速度，ν_e为水流紊动黏性系数，c为谢才系数，f为科氏参量，

ω为地球自转角速度，

为地理纬度，ρ为水密度；S_i为单位水体垂线平均含沙量，

s＝ρ_sc，c为单位水体体积浓度，ρ_s为含沙密度；ν_t和σ_S分别为第一Schmidt数和第二Schmidt数；ω_S为泥沙沉速，下标i表示第i组泥沙，Φ_s为泥沙扩散通量，N_i、

分别为推移质输沙量和推移质输沙能力折算成相应水深的泥沙浓度，β_i为推移质泥沙恢复饱和系数，ω_si为第i组泥沙的沉速，H为水深，H＝h+ζ。

4.根据权利要求3所述的基于数字孪生的长河段滩槽演变预测方法，其特征在于，所述二维水沙数学模型的边界条件包括：

上游进口条件：

下游出口条件：ζ_j2＝ζ_out；

固壁条件：U_j3＝V_j3＝0，

收敛控制条件：

其中，Q_in(t)表示上游来流量，ξ为水位，U_j1和V_j1分别表示上游进口控制点j1在x方向和y方向的速度，h_j1表示上游进口控制点j1处的静水深，dy_j1为离散网格间距，α_j1为动能修正系数，ζ_j2为下手出口控制点j2的水位，ζ_out表示给定出口水位，U_j3和V_j3分别表示固壁的控制点j3在x方向和y方向的速度，S_j3为固壁的控制点j3处的含沙量，N_j3为固壁的控制点j3处的底沙，

为梯度值，b_max为最大质量源，Q_j表示控制点j的断面流量，

表示第n次迭代运算时控制点j和控制点j’之间的x方向的速度差，

表示第n次迭代运算时控制点j和控制点j’之间的水位差，

表示第n+1次迭代运算时控制点j和控制点j’之间的x方向的速度差，

表示第n+1次迭代运算时控制点j和控制点j’之间的水位差。

5.根据权利要求1所述的基于数字孪生的长河段滩槽演变预测方法，其特征在于，所述根据第i年的年初的水文数据，采用离散方法求解所述二维水沙数学模型，获得第i年年末的水流运动数据和河床冲淤厚度数据，具体包括：

根据第i年的年初的水文数据，采用并行运算的方式，采用离散方法求解所述二维水沙数学模型，获得第i年年末的水流运动数据和河床冲淤厚度数据。

6.一种基于数字孪生的长河段滩槽演变预测系统，其特征在于，所述系统包括：

模型构建模块，用于构建长河段的二维水沙数学模型和长河段内的局部桥区的局部正态物理模型；

数学模型计算模块，用于根据第i年的年初的水文数据，采用离散方法求解所述二维水沙数学模型，获得第i年年末的水流运动数据和河床冲淤厚度数据，所述水流运动数据包括x方向流速、y方向流速、悬移质含沙量和水位；

物理模型计算模型，用于基于所述第i年年末的水流运动数据，采用局部正态物理模型进行局部模拟实验，获得第i年的局部桥区的模拟形变量；

局部数据替换模块，用于利用第i年的局部桥区的模拟形变量替换河床冲淤厚度数据中的局部桥区的河床冲淤厚度，得到第i年年末的替换后的河床冲淤厚度数据；

返回模块，用于将第i年年末的水流运动数据和第i年年末的替换后的河床冲淤厚度数据，作为第i+1年年初的水文数据，令i的数值增加1，返回步骤“根据第i年的年初的水文数据，采用离散方法求解所述二维水沙数学模型，获得第i年年末的水流运动数据和河床冲淤厚度数据”，直到到达预测年份。

7.根据权利要求6所述的基于数字孪生的长河段滩槽演变预测系统，其特征在于，所述离散方法包括有限体积法、有限差分法和/或有限元法。

8.根据权利要求6所述的基于数字孪生的长河段滩槽演变预测系统，其特征在于，所述二维水沙数学模型包括：

水流连续方程：

X方向和Y方向水流运动方程：

悬沙不平衡输运方程：

推移质不平衡输移方程：

其中，ζ为水位，t为时间，h为静水深，u、v为流速矢量

沿X、Y方向的分量，g为重力加速度，ν_e为水流紊动黏性系数，c为谢才系数，f为科氏参量，

ω为地球自转角速度，

为地理纬度，ρ为水密度；S为单位水体垂线平均含沙量，

s＝ρ_sc，c为单位水体体积浓度，ρ_s为含沙密度；ν_t和σ_S分别为第一Schmidt数和第二Schmidt数；ω_S为泥沙沉速，下标i表示第i组泥沙，Φ_s为泥沙扩散通量，N_i、

分别为推移质输沙量和推移质输沙能力折算成相应水深的泥沙浓度，β_i为推移质泥沙恢复饱和系数，ω_si为第i组泥沙的沉速，H为水深，H＝h+ζ。

9.根据权利要求8所述的基于数字孪生的长河段滩槽演变预测系统，其特征在于，所述二维水沙数学模型的边界条件包括：

上游进口条件：

下游出口条件：ζ_j2＝ζ_out；

固壁条件：U_j3＝V_j3＝0，

收敛控制条件：

其中，Q_in(t)表示上游来流量，ξ为水位，U_j1和V_j1分别表示上游进口控制点j1在x方向和y方向的速度，h_j1表示上游进口控制点j1处的静水深，dy_j1为离散网格间距，α_j1为动能修正系数，ζ_j2为下手出口控制点j2的水位，ζ_out表示给定出口水位，U_j3和V_j3分别表示固壁的控制点j3在x方向和y方向的速度，S_j3为固壁的控制点j3处的含沙量，N_j3为固壁的控制点j3处的底沙，

为梯度值，b_max为最大质量源，Q_j表示控制点j的断面流量，

表示第n次迭代运算时控制点j和控制点j’之间的x方向的速度差，

表示第n次迭代运算时控制点j和控制点j’之间的水位差，

表示第n+1次迭代运算时控制点j和控制点j’之间的x方向的速度差，

表示第n+1次迭代运算时控制点j和控制点j’之间的水位差。

10.根据权利要求6所述的基于数字孪生的长河段滩槽演变预测系统，其特征在于，所述数学模型计算模块，具体包括：

数学模型计算子模块，用于根据第i年的年初的水文数据，采用并行运算的方式，采用离散方法求解所述二维水沙数学模型，获得第i年年末的水流运动数据和河床冲淤厚度数据。

Images