一种二维非结构干滩条件下流量边界的处理方法与流程

1.本发明属于水利工程技术领域，特别涉及一种二维非结构干滩条件下流量边界的处理方法。


背景技术：

2.二维洪水演进计算是防洪减灾领域的一项重要业务，能够进行准确的模拟洪水演进过程的水动力模型是该领域的一项重要工具。进入二十一世纪以来，以求解riemann近似解为基础的godunov格式由于其强大的激波捕捉和大梯度水面模拟能力，成为洪水演进分析计算领域的主要计算格式。为了适应复杂的计算区域边界，该类模型的求解通常是采用具有良好守恒性能的有限体积法对二维浅水控制方程进行非结构离散。
3.在工程应用中，二维洪水演进模型常用于模拟蓄滞洪区泄洪或大坝溃坝后下游洪水的演进情况。这种情景下，计算区域在初始状态通常是没有水流的，也就是属于干滩的状态，在蓄滞洪区分洪口位置或大坝溃口位置通常会给定一个流量边界条件，然后通过二维洪水演进模型把洪水演进的过程模拟出来。在干滩条件下，如何让水流准确的进入入口各单元是一个难点问题。现在常用的做法是在入口单元内事先给定一个很小的水深值，根据入口单元对应的边占入口边总长度的权重将入流量数值进行分配，使得每个入口单元都有水进入。这种做法的缺陷有两点，首先入口单元给定初始水深值会破坏真实的入流情景，因为入口处各单元的底高程通常是有很大差异的，在计算的初始状态，不是所有的入口单元一开始都会有水流进入的；另一点缺陷是，在进行入口流量分配的时候，进入到每个单元的流量跟单元自身已有的水深是相关的，同样入口边长度下，对应单元水深越大，进入的流量也会越大。如何更精准的处理二维非结构干滩条件下流量边界，依然是当前领域的一个热点和难点问题。


技术实现要素：

4.本发明的目的在于提供一种二维非结构干滩条件下流量边界的处理方法，可以准确的将入口的流量值分配到计算区域的各入口单元内，提高二维非结构干滩条件下流量边界的处理精度，进而提升洪水演进过程整体的计算精度。
5.本发明是通过以下技术方案实现的：
6.一种二维非结构干滩条件下流量边界的处理方法，在计算的初始阶段，将水流入口单元对应的各入口边连接起来组成一维入口控制断面，根据临界流条件计算初始时刻的入流流量值在该断面下对应的水位值，该水位值减去入口单元的底高程即为各入口单元的水深值，根据入口单元的水深值与入口边法线方向流速的对应关系分配各入口单元的入流量；
7.具体步骤如下：
8.1)获取计算所需的基础数据，主要包括地形数据、土地利用类型数据以及入流流量过程数据；
9.2)采用四边形或三角形非结构网格剖分计算区域，剖分完成后，采用平底模型模式对网格高程进行赋值，即采用该区域的地形数据为非结构网格单元的中心点赋值底高程值，在每个计算单元内，地形是平的；根据计算区域的土地利用类型数据为网格单元赋值糙率值；计算的初始水深和流速值均设置为0，满足干滩计算条件；
10.3)选取计算区域的上边界为流量边界，该边界对应的单元称为入口单元，各入口单元对应的边界边为入口边；选取计算区域的下边界为自由出流边界，该边界对应的单元称为出口单元；除入口单元和出口单元外，计算区域内部的单元称为常规单元；
11.4)获取t时刻对应的入口处流量值q
t
，根据cfl条件获取计算时间步长dt；
12.5)不失一般性，在计算的起始步t＝0时，假定q0的值不为0，将入口单元对应的入口边连接起来组成一维入口控制断面，根据临界流条件迭代求解q0对应的控制断面水位z0；该水位值z0如果高于第i个入口单元的底高程zi，则用z0减去zi可以获得该入口单元的水深值，记为hi，反之，该水位值z0如果低于等于第i个入口单元的底高程zi，则该单元的水深hi为0；
13.6)据manning公式可知，计算区域中各处的流速与水深的2/3次幂存在着正比的相关关系，依据该相关关系，根据各入口单元的水深值和入口处流量值计算入口单元对应的入口边法线方向流速；
14.7)根据入口单元的水深值和入口单元对应的入口边法线方向流速，计算通过各入口单元对应的入口边处的法向数值通量值；
15.8)采用完整二维浅水方程组描述洪水运动过程，采用roe格式的近似riemann解来计算区域内部常规单元间的数值通量；进行下边界自由出流边界边数值通量的计算；对通过各单元的数值通量进行求和计算，更新各单元中心水力要素值到t+dt时刻；
16.9)令t＝t+dt，重复步骤4)-8)，直到计算结束；
17.进一步的，步骤5)中的临界流条件如下式所示：
[0018][0019]
式中，a(z0)为入口控制断面水位为z0时对应的过水面积，h(z0)为入口控制断面水位为z0时对应的控制断面平均水深，g为重力加速度。
[0020]
进一步的，步骤6)中各入口单元对应的入口边法线方向流速如下式所示：
[0021][0022]
式中，ui为第i条入口边法线方向流速，li为第i条入口边的长度，n为流量边界处入口单元的总个数。
[0023]
进一步的，步骤7)中通过流量边界第i条入口边的数值通量如下式所示：
[0024][0025]
式中，f
i＊
为第i条入口边法线方向的数值通量，n
x
，ny分别为第i条入口边法线方向的单位向量在x方向和y方向的分量，g为重力加速度。
[0026]
本发明的优点和有益效果是：
[0027]
本发明方法可以准确的确定干滩条件下真实的入口单元个数，可以准确有效的将二维非结构干滩条件下的入口流量值分配到各入流单元内，在保证入口处入流量守恒的条件下，更合理的考虑了干滩入口处水流的运动特点，为二维非结构有限体积法数值框架下干滩情形的入流边界处理提供了一种新的解决途径。
附图说明
[0028]
图1为本发明的一种二维非结构干滩条件下流量边界的处理方法流程图；
[0029]
图2为入口单元对应的入口边相互连接后形成的入口控制断面示意图；
[0030]
图中z0为起始流量q0经临界流条件迭代计算后获得的水位值；
[0031]
图3为toce河算例计算结果对比图；
[0032]
图中(a)为传统流量边界处理方法的水深分布计算结果，图中(b)为本发明所提出的流量边界处理方法的水深分布计算结果。
具体实施方式
[0033]
实施例一：
[0034]
下面结合附图1和实施例对本发明作进一步说明。
[0035]
本发明为一种二维非结构干滩条件下流量边界的处理方法。考虑初始计算时入口单元无水分布的特点，为正确确定入口单元的水流分布，将水流入口单元对应的各入口边连接起来组成一维入口控制断面，根据临界流条件计算初始时刻的入流流量值在该断面下对应的水位值，该水位值减去入口单元的底高程即为各入口单元的水深值。根据入口单元水深值与入口边法线方向流速的对应关系分配各入口单元的入流量。该方法的具体过程如下：
[0036]
1)获取计算所需的基础数据，主要包括比例尺不小于1:1万的dlg地形数据，该数据主要用来提取计算区域的几何边界以及对剖分好的二维网格进行高程插值；分辨率不小于10m的土地利用类型dom数据，该数据主要用来设置计算区域所需的糙率值；以及入流流量过程数据。
[0037]
2)采用四边形或三角形非结构网格剖分计算区域，剖分完成后，采用平底模型模式对网格高程进行赋值，即采用该区域的地形数据为非结构网格单元的中心点赋值底高程值，在每个计算单元内，地形是平的；根据计算区域的土地利用类型数据为网格单元赋值糙率值；计算的水力要素值均设置在网格单元的中心节点，初始水深和流速值均设置为0，满足初始为干滩的计算条件；
[0038]
3)选取计算区域的上边界为流量边界，该边界对应的单元称为入口单元，各入口单元对应的边界边为入口边；选取计算区域的下边界为自由出流边界，该边界对应的单元称为出口单元；除入口单元和出口单元外，计算区域内部的单元称为常规单元；
[0039]
4)获取t时刻对应的入口处流量值q
t
，根据cfl条件获取计算时间步长dt；
[0040][0041]
其中，其中u，v为网格单元中心x，y方向的流速分量，h为网格单元中心水深，g为重
力加速度，n
cfl
为cfl数，dt为计算时间步长，l
l,lr
为网格单元中心到对应边中点之间的距离。
[0042]
5)不失一般性，在计算的起始步t＝0时，假定q0的值不为0，将入口单元对应的入口边连接起来组成一维入口控制断面，控制断面的示意图如附图2所示，根据临界流条件迭代求解q0对应的控制断面水位z0，临界流条件如下式所示：
[0043][0044]
式中，q0为入口处初始时刻流量值，z0为临界流条件下q0对应的水位值，a(z0)为入口控制断面水位为z0时对应的过水面积，h(z0)为入口控制断面水位为z0时对应的入口控制断面平均水深，g为重力加速度。
[0045]
该水位值z0如果高于第i个入口单元的底高程zi，则用z0减去zi可以获得该入口单元的水深值，记为hi，反之，该水位值z0如果低于等于第i个入口单元的底高程zi，则该单元的水深hi为0。
[0046]
6)据manning公式可知，计算区域中各处的流速与水深的2/3次幂存在着正比的相关关系，依据该相关关系，根据入口各单元的水深值hi和入流量q
t
可以获取该入口单元对应的入口边处的法线方向流速ui，具体公式如下式所示：
[0047][0048]
式中，q
t
为t时刻入口处流量值，ui为第i条入口边法线方向的流速，hi为第i个入口单元内的水深值，li为第i条入口边的长度，n为流量边界处入口单元的总个数。
[0049]
7)根据入口处的hi和ui值，可以计算通过各入流单元对应入流边处的法向数值通量值，具体公式形式如下：
[0050][0051]
式中，f
i＊
为第i条入口边法线方向的数值通量，n
x
，ny分别为第i条入口边法线方向的单位向量在x方向和y方向的分量，g为重力加速度。
[0052]
8)采用完整二维浅水方程组描述洪水运动过程，采用roe格式的近似riemann解来计算区域内部常规单元间的数值通量；进行下边界自由出流边界边数值通量的计算；对通过各单元的数值通量进行求和计算，更新各单元中心水力要素值到t+dt时刻；该部分具体的数值求解过程可参见如下著作(张大伟.基于godunov格式的堤坝溃决水流数值模拟[m],2014,中国水利水电出版社)。
[0053]
9)令t＝t+dt，重复步骤4)-8)，直到计算结束。
[0054]
图3为toce河算例计算结果的对比图。该算例的基本情况如下：enel-hydro在意大利米兰市按照1:100的比尺建立了米兰市toce河上游约5km的物理模型，该物理模型作为欧盟支持的cadam项目的一个标准模型，常被用来检验干滩条件下洪水模拟的效果。该物理模型的尺寸大约是50m
×
11m。模型按照5cm的空间步长给出了河道高程点分布，比较精确的描绘了toce河的真实地形。该物理模型的试验材质为混凝土，所以计算糙率采用cadam推荐的
0.0162。该物理模型上端为流量入流边界，下端为一自由出流边界，初始条件下计算区域内的水深为0，满足干滩计算条件。具体的入流边界数值和其它详细描述可参考如下文献(张大伟,堤坝溃决水流数学模型及其应用研究[d]，清华大学,2008.)。计算时，采用0.2m的空间步长均匀离散计算区域，共生成21176个三角形网格单元。图3(a)为采用传统入流边界处理方法获得的水深分布计算结果，图3(b)为采用本发明提出的方法获得的水深分布计算结果，通过对比可以看出，传统处理方法下，入口水流分布不合理，在入口两端的地形较高处也分配了水流，而本发明提出的方法获得的入口水流运动情况与真实情况吻合较好，成功的避免了在初始阶段直接给所有入口单元赋值小水深进行入口流量分配这种传统处理方法所带来的缺陷。综上说明本发明提出的二维非结构干滩条件下流量边界的处理方法是成功的。
[0055]
上述的实施例仅是本发明的部分体现，并不能涵盖本发明的全部，在上述实施例以及附图的基础上，本领域技术人员在不付出创造性劳动的前提下可获得更多的实施方式，因此这些不付出创造性劳动的前提下获得的实施方式均应包含在本发明的保护范围内。