非结构网格数值模拟映射的河流能量公式推导及探析

李彬, 孙东坡

(1.中水北方勘测设计研究有限责任公司,天津 300222; 2.华北水利水电大学 水利学院,河南 郑州 450046)

河流水动力要素的千变万化及河床演变趋势,本质上都是由河流能量的分布及调整转换决定的。从伯努利方程到河流熵理论,河流能量研究始终是一个被关注的话题。目前,国内外学者在河流能量方面的研究主要集中在应用最小能耗率原理或最小熵产生原理分析、预测自然条件下或是修建拦河坝、丁坝等一系列河势控导工程后,河流某一河段的河床整体演变趋势,而没有从能量分布及转换角度研究河流内部的调整规律[1-6]。周筑宝与徐国宾、杨志达(YANG C T)曾就最小能耗率原理或最小熵产生原理可否应用于河流系统研究进行了深入的讨论[7-10],说明将这些涉及河流能量分布及耗散的原理应用到河流系统的理论分析还需进一步完善,应用条件还有待于进一步研究。

由于河流的河床演变过程是一个河段内部(微观)调整、全河段(宏观)再平衡的过程,河床演变必然伴随着河流能量的变化;同样,河流能量的改变也会有与之相匹配的河床变化。因此,可以通过分析河流全河段各区域的能量分布及变化情况来研究河流系统能量的变化趋势,进而预测河流的演变趋势。本文基于计算流体力学理论,从河流网格单元体系的能量分布特征出发,分析了河流能量(机械能)的组成;针对河流网格单元体系中的非结构三角网格单元,分析了三角网格单元内水体的能量(机械能)形式,推导出了计算单元河流机械能各部分的表达式,引入泥沙起动公式对河流机械能的变化性质及河床调整的内部转换特征进行了讨论,并以概化后的顺直河段为例,进行了初步验证和应用研究。

1 基于离散三角单元的河流机械能

在平面二维数值模拟中,流场被离散成三角网格单元体系统。根据欧拉的分析理念,流场的物理属性便可由原来的微元体物理属性转变为三角网格单元体物理属性;水流的总机械能便可由微元体机械能集合转而用三角网格单元体机械能集合表征。每一个离散三角网格单元体的水流机械能包括三角网格单元体的水流重力势能和动能,此部分能量亦可称为动水总流能量,而其中的水流重力势能称为动水重力势能。本文定义三角网格单元体底面(河床床面)高程与重力势能基准之间的流体重力势能为静水重力势能,分别对三角网格单元体内的动水总流能量公式、河床床面所具有的静水重力势能公式以及总机械能公式进行推导,同时引入泥沙起动公式对总机械能的变化性质及内部转换特征进行讨论。

根据MIKE21软件的控制方程组(沿水深平均的二维浅水流动质量和动量守恒控制方程组)及其离散格式,模型计算出的流速V、水位z和水深h等水力要素保存在三角网格单元体的几何重心点上[11]。图1给出了过流断面和水力要素与三角网格单元体的几何关系。其中,图1(b)是图1(a)在水深hn/2处的横截面(三角网格单元)剖分图。图1中:L为过流断面;V为三角网格单元的平均流速；u为三角网格单元沿x轴方向的流速；v为三角网格单元沿y轴方向的流速。

图1 过流断面和水力要素与三角网格单元体之间的几何关系示意

1.1 流场三角网格单元的水流能量推导

1.1.1 标量法网格单元边流速计算

因为u2+v2≠0,所以过流断面L(与流速V方向垂直)与xoy平面的交线L交和x轴的夹角γ分为4种计算情况:

(1)

因为(xj-xi)2+(yj-yi)2≠0,所以三角网格单元各边与x轴的夹角βk分为4种计算情况:

(2)

式中:i=1、2、3;j=1、2、3;k=1、2、3且i≠j≠k,遵循右手定则。例如:A(x1,y1)B(x2,y2)边为边l3,l3与x轴的夹角为β3,l3与交线L交的夹角为α3,通过l3的水流流速为V3,通过l3的流量为Q3。其余两边依次循环类推。

根据三角形外角性质,三角网格单元各边与交线L交的夹角αk=|γ-βk|。根据过流断面交线L交与三角网格单元各边的几何关系,三角网格单元各边的平均流速为:

Vk=V·|cosαk|, 0≤αk<π。

(3)

式中:Vk为通过三角网格单元各边的平均流速,m/s;V为通过三角网格单元的平均流速,m/s。

1.1.2 矢量法网格单元边流速计算

因为(xj-xi)2+(yj-yi)2≠0,所以三角网格单元的各边垂面通过重心点(Xn,Yn)的法向量J在各边上的坐标(Xk,Yk)可由以下方程组计算得到:

(4)

式中:i=1、2、3;j=1、2、3;k=1、2、3且i≠j≠k,遵循右手定则。

根据流速V与三角网格单元各边垂面法向量J的矢量关系及矢量乘法法则,三角网格单元各边的平均流速Vk为:

(5)

式中Vk为正值时表示流入网格,为负值时表示流出网格。

1.1.3 网格单元通量推导

根据三角网格单元各边的平均流速Vk,三角网格单元各边通过的流量可用下式计算:

Qk=Vkhklk。

(6)

式中:Qk为通过三角网格单元各边的流量,m3/s;hk为三角网格单元各边的水深,其等于三角网格单元的平均水深hn,m;lk为三角网格单元各边的边长,m。

因为三角网格单元流入的流量和流出的流量相等,只能是三角网格单元的一个边流入、另两个边流出或是两个边流入、另一个边流出。所以,三角网格单元的流量值等于三角网格单元3条边流量值和的1/2或流量最大边的值,即:

(7)

以质量计的通量可通过类比法应用公式(7)确定。

1.1.4 网格单元能量方程推导

根据伯努利能量方程[12],每个单元格的动水总流能量方程为:

(8)

对于天然河流,Pn取自由面相对压强,Pn=0;取动能修正系数αn等于1则:

(9)

式中:ρ为水的密度,ρ=1 000 kg/m3;g为重力加速度,g=9.81 m2/s;Qn为通过三角网格单元的流量,m3/s;Hn为三角网格单元的总水头,m;zn为三角网格单元的平均水位,m;z0为计算重力势能的基准水位,m;Vn为三角网格单元的平均流速,m/s;n为网格单元编号。

1.2 静水三角网格单元的重力势能推导

结合前述的静水重力势能定义,由图1可知静水重力势能的计算公式为:

(10)

其中，

式中:e1静n为静水重力势能,kJ;sn为三角网格单元的面积,m2。

1.3 三角网格单元总机械能

假设河床泥沙为散粒体均匀沙,泥沙起动流速可用沙莫夫公式确定[13]，如下:

(11)

式中:Ue为泥沙起动流速,m3/s;ρs为泥沙的密度,一般取2 650 kg/m3;d为泥沙粒径,m。

根据河流动力学原理,当泥沙起动流速Ue等于三角网格单元的水流平均流速Vn时,可以得到河床冲刷稳定条件下的三角网格单元极限平均水深hn，

(12)

将式(12)代入式(10)可得静水重力势能e1静n的计算公式为:

(13)

由公式(9)与公式(13)可知:Vn增加,E动n增加,e1静n减小;反之,E动n减小,e1静n增加。水流是泥沙输移的原动力,河道水流运动的天然不稳定性体现在流速的波动上,这种波动必然导致动水总流能量与静水重力势能的相互转换;由于河床泥沙运动相较于水流运动具有一定的滞后性(两者波动有相位差),这种滞后性又会进一步影响能量的转换特征。根据河流最小能耗原理[3],总机械能En总是趋向最小化,其计算公式如下:

En=E动n+e1静n=e1动n+e2动n+e1静n→min。

(14)

式中:En为总机械能;符号“→”表示趋向性;“min”表示最小值。

应用公式(9)和公式(10)或公式(9)和公式(13)对模拟河段每个三角网格单元不同形式的能量分别进行计算并求和,便可得到模拟河段处于冲刷平衡状态时不同形式能量的数值,还可得到其分布状况。通过流场能量的分布特征及其调整转换情况,分析河流工程区域三角网格单元的能量变化情况,就可以对河流水动力的变化及演变趋势进行预测。

2 用于能量分析的流场质量守恒验证

2.1 模拟条件

构造长为1 000 m，宽为100 m，渠底高程为-1.500 m的平坡顺直河道作为模拟河段,对推导的公式进行流场质量守恒验证。模型进口流量为100 m3/s,出口水位为0.500 m,模型网格单元节点个数为1 390,三角形网格单元总数为2 558。模型四边固定边界的网格单元边长设置为10.0 m,区域内网格单元面积设置为不超过60.0 m2,并剖分生成网格单元。模型局部网格剖分及出口断面网格单元编号如图2所示。

图2 模型局部网格剖分及出口断面网格单元编号示意

2.2 结果及分析

应用MIKE21软件进行平面二维水动力数值模拟,得到模型出口断面网格单元的水力要素(表1),据此采用三角网格单元流量公式计算得到模型出口断面的网格单元流量(表2和表3)。应用MIKE21软件自带的断面流量监测功能监测到的模型出口断面流量为99.909 7 m3/s,比给定流量偏小0.090 3 m3/s,误差非常小,仅为0.903‰,可以认为MIKE21软件计算出的流量守恒。

由表2和表3可知:①采用两种三角网格单元流量计算方法得到的流量模拟计算值在小数点后保留4位有效数字时没有差别,两种方法的计算精度相当。②应用本文流量公式得到的模型出口断面三角网格单元的流量和为100.678 4 m3/s,比给定流量稍大。分析其原因是:在模型出口边界三角网格单元系统中存在钝角三角形网格单元,导致相邻网格单元的流量有重复计量。从图2可以看出:33#网格单元为钝角三角形网格单元,造成它与相邻的1 552#网格单元之间有通量传递,部分流量计重。若将Q网33替换为33#网格单元出口断面的单边流量Q3,扣除其与1 552#网格单元之间的流量重计量,则模型出口断面流量为99.790 7 m3/s,仅比给定流量偏小0.209 3 m3/s,误差为2.093‰。因此,在对模型下开边界附近的网格进行剖分时,应注意避免出现钝角三角形网格单元,这样采用文中所给的两种计算三角网格单元流量的方法都可以获得满意的计算精度。

通过以上分析基本可以认为:应用推导出的两种三角网格单元流量计算方法得到的河流流量均守恒;模拟结果真实可信,两种方法均可用于三角网格单元的通量及能量计算。

表1 模型出口断面网格单元的水力要素统计

表2 由公式计算出的模型出口断面网格单元流量(流速标量法)

表3 由公式计算出的模型出口断面网格单元流量(流速矢量法)

3 模拟河段的能量分布及分析

对上节概化顺直河道的水动力模拟计算结果应用公式(9)和公式(10)计算模拟河段的能量,其中z0取模拟河段河底高程最低值z河底min(-1.500 m),以保证所有计算三角网格单元的势能均不小于零。计算模拟河段每个三角网格单元的动水重力势能e1动n(包含压能)、动水动能e2动n、动水机械能E动n、静水重力势能e1静n和总机械能En,累加求和后得到模拟河段不同形式能量的值,结果见表4。图3给出了对基于流速标量法的三角网格单元流量结果应用公式(9)和公式(10)计算得到的模拟河段不同形式能量的分布,图4给出了模拟河段不同形式能量的面密度(单位面积上的能量)分布。

表4 概化模拟河段的单位时间能量 kJ

图3 概化模拟河段的单位时间能量分布

图4 概化模拟河段的单位时间能量面密度分布

因为模拟河段是基于河床平坡的概化顺直河道,且选取的势能计算基准是z河底min,所以e1静n等于零,因此En与E动n的数值相等、分布相同。分析表4可知:①基于流速标量法、流速矢量法的流量结果,应用公式(9)和公式(10)计算得到的模拟河段不同形式能量的值几乎无差别,计算精度相当。②选取z河底min作为势能的计算基准,在本次模拟条件(常水位)下,模拟河段总体上e1动占E动的比例为994.109‰,e2动占E动的比例为5.891‰。这表明:当选取z河底min作为势能的计算基准时,在常水位条件下模拟河段所具有的能量主要是动水重力势能e1动。这也从理论上给出了洪水期的河流特别是地上悬河的溃堤水流为什么会具有惊人破坏力的一个佐证。

图3反映了模拟河段不同形式河流能量的平面分布特征。图3表明:河流能量在空间上不是均匀分布的,沿纵横向(网格单元体系间)都存在能量梯度,亦即都存在能量的传递与转换。

图4给出的不同形式河流能量的面密度分布特征,显示了比图3更为连续的分布变化特点,同时显示出能量纵向分布有波动变化特征,也间接体现了沿程的能耗特点;动能和总机械能的纵横向变化也表征了各网格单元体系间的能量梯度与传递特性。图3和图4表明:河流能量的损失一方面与边界条件有关;另一方面也遵循最小能耗原理,即总是使能量损失趋向最小化,它决定了能量耗散形式,同时也制约着流场中不同形式能量间的分配关系。

基于河流网格单元体系的能量分布特征,必然带有网格的特征并会受其影响:模拟河段网格剖分不均匀,并非均一等面积网格;网格单元能量特征值受网格精度影响。这是需要在前期的网格剖分以及在模拟计算结果分析中必须关注的问题。本文旨在提出一种基于计算流体力学的河流能量研究思路和方法,对于网格布置与精度(单元尺寸及其均匀化程度)对模拟结果的影响,以及通过分析三角网格单元能量的变化预测河流演变趋势,将在后续的研究中进行。

4 结语

为了研究基于模拟流场的河流能量分布及转换情况,通过分析河流能量(机械能)的组成,利用非结构三角网格单元推导出了数值模拟区河流能量(机械能)的表达式,采用水动力数值模拟软件构造了概化顺直河道模型,并对推导的公式进行了不同计算方法的守恒性验证和模拟计算,得到以下结论。

1)从能量角度分析,河流能量(机械能)包括由流动水体产生的动水总流能量和河床床面高程具备的静水重力势能。其中,动水总流能量包括动水重力势能(包含压能)和动水动能。河道水流的波动导致动水总流能量与静水重力势能频繁转换,河床泥沙运动相较于水流运动的滞后性又进一步影响能量的转换特征。根据河流最小能耗原理,河流能量(机械能)的消耗总是趋向最小化,它决定了能量耗散形式和不同形式能量间的分配关系。

2)河流能量在空间上不是均匀分布的,沿纵横向(网格单元体系间)都存在能量梯度,亦即都存在能量的传递与转换。河流能量面密度的分布特征显示出能量纵向分布有波动变化特征,也间接体现了沿程的能耗特点;动能和总机械能的纵横向变化也表明各网格单元体系间存在能量梯度,河流能量存在波动性分布不均与空间的相互传递。

3)基于非结构三角网格单元推导出的两种三角网格单元流量计算方法均可用于三角网格单元的流量计算,且计算精度相当。应用流速矢量法计算三角网格单元流量可使三角网格单元各边的过流情况(大小及流向)一目了然。

4)基于非结构三角网格单元推导出的河流能量(机械能)公式可用于河流能量计算。在基于非结构三角网格单元的河流平面二维数值模拟中,可以通过该公式得到河段能量(机械能)的数值及其分布态势。

5)本文研究对象为概化的顺直河道模型,对于河道形态不规则的实际河流,三角网格单元尺寸及其均匀化程度对河流全河段能量计算结果的影响，以及冲刷平衡的极限状态的讨论有待于进一步研究。