基于数值模拟的矩形明渠床面切应力分析

黄宇航，周晓泉，周文桐，姚 畅，周芳龄

(四川大学 水力学与山区河流开发保护国家重点实验室，四川 成都 610065)

人工修建或自然形成的渠道被称为明渠，在该渠道流动并拥有与大气接触的自由表面的流动称为明渠流，明渠流表面上各点的压强皆为大气压强。对于明渠水流，无论是在工程实践或是理论研究中，边界切应力都是一个重要参数，水流阻力、泥沙运动、污染物扩散等水力学问题都需要关注边界切应力。

国内外许多学者对矩形明渠边界切应力进行了试验研究。惠遇甲等通过矩形水槽试验分析了光滑明渠流和粗糙明渠流阻力系数的变化、分割，以及剪切应力沿床面和边壁的分布规律。Knight等根据矩形光滑水槽试验得出的边界剪应力大小及分布情况，提出了墙体所受总剪力的百分比与宽深比的函数关系，并给出了平均壁面、床层和床层中心线剪应力随长宽比的函数表达式。胡旭跃等通过对矩形明渠均匀流运动的力学分析，探讨了平均剪应力沿边界的分布问题，对Knight等的公式进行了改进；此外，胡旭跃等还在分析Rhodes、Ghosh等的资料基础上发现爱因斯坦阻力公式可适用边壁和床面阻力计算，若将连续方程、几何关系及爱因斯坦阻力公式联合使用，便可解决阻力分割计算中的方程组不封闭问题。张书农、唐存本对水槽底部和侧壁的关系进行了理论推导，得出了消除侧壁影响的槽底阻力计算式，该式在任何粗糙情况下均适用。陈小芳、何建京在分析Preston管量测试验资料的基础上，测量了不同水力条件下均匀流与非均匀流的底壁切应力，试验表明：均匀流时，切应力与坡度及流量呈正相关；非均匀流时，在相对水深相同时，其相对底壁切应力大致相同。黄纪中进行的水槽试验结果表明，明渠槽底平均切应力与弗劳德数具有函数关系，可供矩形水槽试验进行参考。Pedro等将两种算法递归图（RP）和特征面识别算法合并使用，并与一个标准回归模型相结合，来预测明渠边界切应力分布，并且将该方法应用于梯形、圆形渠道边界切应力分布预测，证明了该方法可以准确估算明渠边界切应力分布。Kamalini等对不同糙率下的复合明渠边界切应力做出了研究，提出了一种利用边界剪切分布和动量传递系数计算复合通道交界处流动分布的方法。Kamalin等又提出了一种改进的直线型复合河道深度平均流速和床层剪切应力预测方法，得到了N-S方程深度积分湍流形式的解析解，并将该解析解成功地应用于多种复合通道的试验和现场实例。Fatemeh等建立了复合明渠弯道边界剪应力计算的解析模型，利用达西公式计算边界剪应力分布，并与梯形和矩形渠道、直道和弯曲道试验数据比较，证明建立的解析模型可用于简单明渠和复合明渠的边界剪应力预测。Behzad等提出了一种新的边界剪应力模型，利用明渠的速度分布估计边界剪应力的分布。

迄今为止，学者们对边界切应力进行了大量的试验与研究，并提出了许多方法解读和预测边界切应力的大小和分布。但是，即便在简单的光滑均匀流棱柱明渠里，由于断面形状、河道坡度、水流结构的复杂性，使得测量和计算边界切应力仍然具有挑战性。大多数的边界切应力求解都是建立在宽深比之上，或者是进行直接试验得出经验公式，本文通过建立直道矩形明渠模型，进行均匀流和非均匀流受力分析，进一步归纳出床面切应力通用公式。

1 数学模型

1.1 控制方程

通常的平面2维流动模型的基本方程如下：

水流连续方程：

水流动量方程：

式（1）～（3）中：

U

、

V

为垂线平均流速在

X

、

Y

方向的分量；

z

为水位；

H

为水深；

C

为谢才系数；

f

为柯氏力系数，

f

=2ωsin

φ

（ ω为 地球自转角速度，

φ

为当地纬度）；

g

为重力加速度；

v

为湍动黏性系数。

1.2 模型建立

河流的流动问题极为复杂，除了重力这个绝对的驱动因素外，还受到河底与河岸对水流的阻碍作用、风对水体的拖拽作用等因素的影响。若去掉各种不确定因素后，便成为矩形明渠，只有渠道宽

B

、坡降

J

、和糙率

n

这3个因素。本文将直道明渠设计为宽度

W

=0.5 m，长度

L

=20 m，比降

J

=0.1%，河床糙率

n

=0.015，将其定义为R2（Rectangular）模型；放大10倍、100倍的模型（按相似理论中的重力相似准则进行缩尺，其中糙率系数的缩尺为 λ=λ）分别定义为R3、R4模型，比降皆为0.1%；模型R2J尺寸同R2，比降为0.5%；模型R4J尺寸同R4，比降为0.05%。各模型的基本参数见表1。

表1 各模型参数
Tab. 1 Properties of each model

模型编号 比尺 长度 L/m 宽度 B/m 坡降 J 糙率n R2 1.0 20.0 0.5 0.001 0 0.015 0 R3 10.0 200.0 5.0 0.001 0 0.022 0 R4 100.0 2 000.0 50.0 0.001 0 0.032 3 R2J 1.0 20.0 0.5 0.005 0 0.015 0 R4J 100.0 2 000.0 50.0 0.000 5 0.032 3

1.3 网格及边界条件设置

各比尺的直道模型网格均按宽度方向10等分，长度方向100等分进行划分，共划分了1 000个8节点的四边形网格。

在上游进口给定流量条件，在出口给定水位条件，不同比尺的水位条件严格按阻力相似准则规定。先固定下游水位，通过试算判定流场是否为均匀流，如果不是，则调整流量，直至流场为均匀流。

在数值模拟计算中，直道明渠水深从上游往下游总是深-浅-再深-再浅，如果上、下游同相位的“深”水位的差值较小，则认为流动是均匀流。设置出口边界条件后，若上游水深大，则流量过大，必须降低流量；反之，加大流量。本文将明渠上下游水位差(

H

-

H

)/

H

小于1×10，作为均匀流的判别标准。

各模型典型水深见图1。R2的水深变化有一个多周期；R3有十多个周期；R4应有百多个周期，但只有100个网格，已经没法分辨水深周期，因此按预定位置水深进行评估来判断是否为均匀流。

图1 各模型的典型水深Fig. 1 Typical water depth of each model

2 计算结果分析

一般的明渠均匀流公式（谢才公式）如下：

式中，

A

为断面面积，

R

为渠道断面水力半径，

J

为河道比降或坡降，

C

为谢才系数。

曼宁公式定义如下：

式中，

n

为曼宁系数或糙率，通常认为是常数。当渠道宽无穷大，或者不考虑渠道边壁影响，同时均匀流水深处处相等时，水力半径

R

即为渠道水深

H

，这时的曼宁公式可以写为：

将平面2维流动模型中的水力半径

R

用水深

H

替代，便有：

2.1 直道明渠均匀流计算结果分析

将R2模型定义为标准算例，当流量为0.1 m/s时，水位经试算确定为0.243 5 m，其他模型均按与标准算例对应的比尺模型计算。

按不同比尺、不同比降的模型进行稳定流数值模拟计算，各算例统计如表2所示。

表2 直道明渠各比尺模型数模结果
Tab. 2 Simulation results of each scale of straight open channel

R2 R3 R4 R2J R4J H/m Q/(m3·s-1) H/m Q/(m3·s-1) H/m Q/(m3·s-1) H/m Q/(m3·s-1) H/m Q/(m3·s-1)0.020 00.001 3— — — —0.153 60.103 1205 098.250 0.243 50.100 02.43531.670 524.35010 008.800——5023 475.000 0.200 00.072 02.00022.814 620.0007 210.120— —0.300 00.141 63.00044.841 530.00014 171.100— — — —0.400 00.228 94.00072.428 040.00022 888.500— — — —0.500 00.332 05.000105.052 050.00033 199.000— — — —2.000 03.346 520.0001 058.750 0— — — —5.000 015.408 050.0004 875.700 0— — — —

若各比尺模型按

Q

·

n

/

W

/

J

统计，则所有直道明渠均匀流数据可以归纳出统一的明渠流均匀流的公式：

式（9）中的

Q

/

B

即 为式（6）中的单宽流量

H

·

u

，因此，可以认为该模型模拟出来的参数是可靠的。

2.2 直道明渠均匀流受力分析

如果对直道明渠应用动量定理，将进口断面定义为1-1断面（水深

h

、流速

u

、水体质量

m

），出口定义为2-2断面（水深

h

、流速

u

、水体质量

m

），则有：

因为是均匀流，有

h

=

h

=

H

=

H

， 且

m

=

m

，

u

=

u

，重力分量与阻力匹配，简化后有均匀流河床切应力公式：

式中，

H

为均匀流水深， τ为均匀流河床切应力。将式（11）和式（9）相结合，可以得到其他的切应力表达方式：当消除坡降因子

J

时，均匀流河床切应力公式表示为：

式中，

u

为均匀流速度。当消除水深因子

H

时 ，均匀流河床切应力公式表示为：

也可用单宽流量

H

·

u

表达均匀流河床切应力公式：

2.3 直道明渠非均匀流受力分析

天然河流的流动远没有达到均匀流。如果当水深

H

=

H

时为均匀流， 那么当同样流量下

H

>

H

时，或下游与均匀流同样水位的条件下，降低进口流量，一定有一个恒定流的解没有达到均匀流，需提高流量才能达到，所以命名为欠均匀流；同理，当

HH

为超均匀流，或均匀流条件下，保持下游水位不变，将上游流量增加，这样的情况下也有一恒定流的解，因为流量高于均匀流，所以称为超均匀流。以表1为模板，计算9个非均匀流的算例，下游水位保持不变，欠均匀流流量调至均匀流的1/4至1/2；超均匀流流量调至均匀流的2至5倍，统计结果见表3，并将其恒定流的解（进口水深

h

、出口水深

h

）一并列出。

表3 非均匀流各算例统计结果
Tab. 3 Statistical results of each example of non-uniform flow

流态 模型 算例 H/m Q/(m3·s-1) h1/m h2/m b 2.00 10 1.835 2.000 2.334 2 5.00 25 4.808 5.000 2.333 3 20.00 200 19.803 20.000 2.339 4 50.00 1 000 49.802 50.000 2.344 1 R3欠均匀流R4 5 24.35 2 500 22.455 24.350 2.334 6 50.00 10 000 48.129 50.000 2.334超均匀流 R3 7 2.00 50 2.792 2.030 2.042 8 5.00 500 4.389 4.999 2.332 R4 9 24.35 40 000 15.960 24.350 2.333

在整个计算域的水体上应用动量定理，并对沿流动方向的100个网格单元的流动数据进行积分，即应用如下表达式：

式中： ρ为水体密度，

g

为重力加速度，

h

、

h

分别为进、出口水流水深，

u

、

u

分别为进、出口的水流流速，

J

为河道比降，

l

为河道长度；等号左侧前两项分别为进口水压力

P

、出口水压力

P

，第3项为整个水体重力的水平分量

G

，第4项为河床阻力

F

，等号右侧为动量改变量Δ

M

。除了切应力τ 之外，所有的量均可以准确统计出来，各算例具体受力见表4。

表4 各算例具体受力情况
Tab. 4 Specific force of each calculation example

流态 模型 算例 P1/N P2/N Gf/N ΔM/(N·s) Fb/N 1 16 485.018 19 583.960 3 753.390 -179.578 834.027 2 113 180.823 122 403.296 9 604.036 -199.375 580.939 3 1 920 101.797 1 958 464.483 38 976.548 -793.927 1 407.789 4 12 143 606.990 12 240 412.810 97 729.376 -3 178.663 4 102.217 R4 5 2 468 821.317 2 903 047.565 458 283.142 -8 649.921 32 706.815 6 11 341 342.830 12 240 437.290 960 880.390 -31 050.683 92 836.609超均匀流 R3 7 38 162.681 20 171.841 5 106.255 13 425.185 9 671.910 8 94 330.095 122 334.319 9 194.422 -277 191.780 258 381.979 R4 9 1 247 152.697 2 902 899.732 394 795.643 -13 791 541.800 12 530 590.400 R3欠均匀流

假定非均匀流的切应力公式与式（12）具有相同的形式，速度增加则切应力增加，反之则减少，设其与速度的

b

次方成正比，因此，切应力的表达式如下：

将式（16）代入式（15），计算得出的指数

b

的结果见表3，显然除了第7算例外，其余8个算例均一致，取其平均数

b

≈2.335 5。

则非均匀流的河床切应力有如下的表达式：

当消除坡降因子

J

时，非均匀流河床切应力公式表达为：

式中，

u

为非均匀流流速。当用单宽流量

H

u

表达时，非均匀流河床切应力的计算公式为：

当用流速

u

表达时，非均匀流河床切应力的计算公式为：

由表3可知：在欠均匀流条件下，式（17）～（19）均成立；在超均匀流条件下，除第7算例外，均成立。成立的算例的共同特点，为

h

h

，即“欠”均匀流状态；第7算例为超均匀流边界条件下的“超”均匀流，然而，现实中是否存在这样的流动还需用试验验证。由表4还可以看出，非均匀流的重力分量

G

和阻力

F

差异较大，这也是导致壅水或跌水的根本原因。因此，式（17）～（19）对所有的直道河流流动均是普适的，尤其是式（17），只要知道流速

u

、水深

h

和糙率

n

，无论其坡降如何，便能求出其床面切应力。将式（12）和（16）代入式（17），确定

b

=7/3，消去

J

求得通用的床面切应力公式：

式（20）同均匀流床面切应力式（12）结构相同，因此，无论流动是否为均匀流，式（20）均成立。通用的切应力公式可以将复杂的流动问题简单化，将流体力学问题（求解偏微分方程）变为简单的代数问题。

3 公式对比分析

许多学者都对非均匀流床面切应力进行过研究，也给出了一些经验公式。陈小芳等由试验结果推导出相对河床切应力 τ/τ与相对水深

h

/

h

存在大致相同的经验公式：

金中武等也通过水槽试验拟合出雍水条件下的经验公式：

张小峰等通过数模拟合出了经验公式：

为方便比较，将本文公式（16）改写为如下形式：

将式（21）～（24）应用到表3的9个非均匀流的数据算例中，对比5个公式，验证本文公式的合理性。对各公式切应力沿流动方向的100个单元进行积分，得出河床阻力见表5。由表4、5可以看出：当 τ >τ时，式（21）～（24）计算的结果比动量定理推导出的河床阻力偏小；反之，计算结果偏大。显然，式（16）（即式（25））既符合动量定理又有物理意义。

表5 各公式计算的阻力
Tab. 5 Resistance calculated by each formula

N流态 模型 算例 式（25） 式（21） 式（22） 式（23） 式（24）1 834.027 1 039.185 946.414 883.168 1 139.246 2 580.939 949.605 744.006 890.006 996.897 R3 3欠均匀流1 407.789 2 866.153 2 067.704 2 869.829 2 654.673 4 4 102.217 7 803.315 5 767.432 7 668.606 7 562.515 5 R4 32 706.810 50 634.873 40 741.571 46 414.716 54 500.440 6 92 836.610 132 496.511 111 465.942 117 059.312 146 780.100超均匀流7 R3 9 671.912 8 902.099 8 112.135 9 482.511 9 535.652 8 258 382.000 81 811.101 70 450.180 110 503.622 151 290.300 R4 9 12 530 590.000 3 651 916.326 3 139 404.646 4 964 744.890 7 100 289.000

陈小芳等的试验中1个流量只对应4个水位，数据较少，并且由于非均匀流的水位波动较大，测量出来的水位对拟合结果也会有影响。张小峰等拟合出来的公式具有明确的物理意义，但并未用动量定理验证。本文每个算例的每个流量都取了划分出来的1 000个8节点四边形网格中的直道中心处间隔相同的101个点位进行受力分析，每个点位都严格遵守动量方程，切应力同进出口水压力、水体重力的水平分量及动量的改变三者之和相匹配；每两个点之间取1次切应力的平均值，最后再将所有点位的阻力求和，以近似于沿程积分出河床阻力，这将大大减少水位波动等不确定因素带来的误差影响，最终8个算例都得出相同的结论：相对河床切应力 τ/τ与相对水位（

h

/

h

）遵循式（25），因此，可以认为本文推导公式更合理且有物理意义。

4 结 论

本文对5种比尺的均匀流直道模型以及4种比尺不同流量下的非均匀流直道模型分别在均匀流和非均匀流情况下进行了数值模拟，对分析模拟出来的数据得出以下结论：

1）通过直道明渠均匀流的模拟数据结果，根据所有直道明渠的水位流量关系图可以归纳出明渠均匀流公式：

Q

·

n

/

B

/

J

=

H

，并在此基础上进行受力分析得出均匀流切应力表达式（消除坡降

J

）：τ=ρ

g

(

nu

)/

H

。2）假定明渠非均匀流切应力与非均匀流有同样的形式 τ=(

u

/

u

)τ， 对直道明渠非均匀流的模拟数据结果进行分析，比例因子 (

u

/

u

)中指数

b

=3/7，最终得到消去

J

的通用的床面切应力公式： τ =ρ

g

(

nu

)/

h

，该公式也同样适用于垂直立面的切应力计算。

3）本文公式是基于动量定理及运用受力分析严格推导出来的，绘制出来的曲线与陈小芳、金中武等用试验数据拟合出来的曲线相似，但当 τ>τ时，式（21）～（24）计算的结果比动量定理推导出的河床阻力偏小；反之，计算结果偏大。

4）由本文推导出的床面切应力公式可知，无论矩形明渠的流动是否均匀，均可通过水深

h

、流速

u

、糙率

n

这3个参数求出床面切应力的代数解，相比于已有非均匀流的切应力公式，本文推导出的通用切应力计算公式更加直接，且已在平面2维流动模型中应用。