基于數(shù)值模擬的矩形明渠床面切應力分析

黃宇航，周曉泉，周文桐，姚 暢，周芳齡

(四川大學 水力學與山區(qū)河流開發(fā)保護國家重點實驗室，四川 成都 610065)

人工修建或自然形成的渠道被稱為明渠，在該渠道流動并擁有與大氣接觸的自由表面的流動稱為明渠流，明渠流表面上各點的壓強皆為大氣壓強。對于明渠水流，無論是在工程實踐或是理論研究中，邊界切應力都是一個重要參數(shù)，水流阻力、泥沙運動、污染物擴散等水力學問題都需要關(guān)注邊界切應力。

國內(nèi)外許多學者對矩形明渠邊界切應力進行了試驗研究。惠遇甲等通過矩形水槽試驗分析了光滑明渠流和粗糙明渠流阻力系數(shù)的變化、分割，以及剪切應力沿床面和邊壁的分布規(guī)律。Knight等根據(jù)矩形光滑水槽試驗得出的邊界剪應力大小及分布情況，提出了墻體所受總剪力的百分比與寬深比的函數(shù)關(guān)系，并給出了平均壁面、床層和床層中心線剪應力隨長寬比的函數(shù)表達式。胡旭躍等通過對矩形明渠均勻流運動的力學分析，探討了平均剪應力沿邊界的分布問題，對Knight等的公式進行了改進；此外，胡旭躍等還在分析Rhodes、Ghosh等的資料基礎上發(fā)現(xiàn)愛因斯坦阻力公式可適用邊壁和床面阻力計算，若將連續(xù)方程、幾何關(guān)系及愛因斯坦阻力公式聯(lián)合使用，便可解決阻力分割計算中的方程組不封閉問題。張書農(nóng)、唐存本對水槽底部和側(cè)壁的關(guān)系進行了理論推導，得出了消除側(cè)壁影響的槽底阻力計算式，該式在任何粗糙情況下均適用。陳小芳、何建京在分析Preston管量測試驗資料的基礎上，測量了不同水力條件下均勻流與非均勻流的底壁切應力，試驗表明：均勻流時，切應力與坡度及流量呈正相關(guān)；非均勻流時，在相對水深相同時，其相對底壁切應力大致相同。黃紀中進行的水槽試驗結(jié)果表明，明渠槽底平均切應力與弗勞德數(shù)具有函數(shù)關(guān)系，可供矩形水槽試驗進行參考。Pedro等將兩種算法遞歸圖（RP）和特征面識別算法合并使用，并與一個標準回歸模型相結(jié)合，來預測明渠邊界切應力分布，并且將該方法應用于梯形、圓形渠道邊界切應力分布預測，證明了該方法可以準確估算明渠邊界切應力分布。Kamalini等對不同糙率下的復合明渠邊界切應力做出了研究，提出了一種利用邊界剪切分布和動量傳遞系數(shù)計算復合通道交界處流動分布的方法。Kamalin等又提出了一種改進的直線型復合河道深度平均流速和床層剪切應力預測方法，得到了N-S方程深度積分湍流形式的解析解，并將該解析解成功地應用于多種復合通道的試驗和現(xiàn)場實例。Fatemeh等建立了復合明渠彎道邊界剪應力計算的解析模型，利用達西公式計算邊界剪應力分布，并與梯形和矩形渠道、直道和彎曲道試驗數(shù)據(jù)比較，證明建立的解析模型可用于簡單明渠和復合明渠的邊界剪應力預測。Behzad等提出了一種新的邊界剪應力模型，利用明渠的速度分布估計邊界剪應力的分布。

迄今為止，學者們對邊界切應力進行了大量的試驗與研究，并提出了許多方法解讀和預測邊界切應力的大小和分布。但是，即便在簡單的光滑均勻流棱柱明渠里，由于斷面形狀、河道坡度、水流結(jié)構(gòu)的復雜性，使得測量和計算邊界切應力仍然具有挑戰(zhàn)性。大多數(shù)的邊界切應力求解都是建立在寬深比之上，或者是進行直接試驗得出經(jīng)驗公式，本文通過建立直道矩形明渠模型，進行均勻流和非均勻流受力分析，進一步歸納出床面切應力通用公式。

1 數(shù)學模型

1.1 控制方程

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

水流連續(xù)方程：

水流動量方程：

式（1）～（3）中：

U

、

V

為垂線平均流速在

X

、

Y

方向的分量；

z

為水位；

H

為水深；

C

為謝才系數(shù)；

f

為柯氏力系數(shù)，

f

=2ωsin

φ

（ ω為 地球自轉(zhuǎn)角速度，

φ

為當?shù)鼐暥龋?p>g

為重力加速度；

v

為湍動黏性系數(shù)。

1.2 模型建立

河流的流動問題極為復雜，除了重力這個絕對的驅(qū)動因素外，還受到河底與河岸對水流的阻礙作用、風對水體的拖拽作用等因素的影響。若去掉各種不確定因素后，便成為矩形明渠，只有渠道寬

B

、坡降

J

、和糙率

n

這3個因素。本文將直道明渠設計為寬度

W

=0.5 m，長度

L

=20 m，比降

J

=0.1%，河床糙率

n

=0.015，將其定義為R2（Rectangular）模型；放大10倍、100倍的模型（按相似理論中的重力相似準則進行縮尺，其中糙率系數(shù)的縮尺為 λ=λ）分別定義為R3、R4模型，比降皆為0.1%；模型R2J尺寸同R2，比降為0.5%；模型R4J尺寸同R4，比降為0.05%。各模型的基本參數(shù)見表1。

表1 各模型參數(shù)
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 網(wǎng)格及邊界條件設置

各比尺的直道模型網(wǎng)格均按寬度方向10等分，長度方向100等分進行劃分，共劃分了1 000個8節(jié)點的四邊形網(wǎng)格。

在上游進口給定流量條件，在出口給定水位條件，不同比尺的水位條件嚴格按阻力相似準則規(guī)定。先固定下游水位，通過試算判定流場是否為均勻流，如果不是，則調(diào)整流量，直至流場為均勻流。

在數(shù)值模擬計算中，直道明渠水深從上游往下游總是深-淺-再深-再淺，如果上、下游同相位的“深”水位的差值較小，則認為流動是均勻流。設置出口邊界條件后，若上游水深大，則流量過大，必須降低流量；反之，加大流量。本文將明渠上下游水位差(

H

-

H

)/

H

小于1×10，作為均勻流的判別標準。

各模型典型水深見圖1。R2的水深變化有一個多周期；R3有十多個周期；R4應有百多個周期，但只有100個網(wǎng)格，已經(jīng)沒法分辨水深周期，因此按預定位置水深進行評估來判斷是否為均勻流。

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

2 計算結(jié)果分析

一般的明渠均勻流公式（謝才公式）如下：

式中，

A

為斷面面積，

R

為渠道斷面水力半徑，

J

為河道比降或坡降，

C

為謝才系數(shù)。

曼寧公式定義如下：

式中，

n

為曼寧系數(shù)或糙率，通常認為是常數(shù)。當渠道寬無窮大，或者不考慮渠道邊壁影響，同時均勻流水深處處相等時，水力半徑

R

即為渠道水深

H

，這時的曼寧公式可以寫為：

將平面2維流動模型中的水力半徑

R

用水深

H

替代，便有：

2.1 直道明渠均勻流計算結(jié)果分析

將R2模型定義為標準算例，當流量為0.1 m/s時，水位經(jīng)試算確定為0.243 5 m，其他模型均按與標準算例對應的比尺模型計算。

按不同比尺、不同比降的模型進行穩(wěn)定流數(shù)值模擬計算，各算例統(tǒng)計如表2所示。

表2 直道明渠各比尺模型數(shù)模結(jié)果
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

統(tǒng)計，則所有直道明渠均勻流數(shù)據(jù)可以歸納出統(tǒng)一的明渠流均勻流的公式：

式（9）中的

Q

/

B

即 為式（6）中的單寬流量

H

·

u

，因此，可以認為該模型模擬出來的參數(shù)是可靠的。

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

如果對直道明渠應用動量定理，將進口斷面定義為1-1斷面（水深

h

、流速

u

、水體質(zhì)量

m

），出口定義為2-2斷面（水深

h

、流速

u

、水體質(zhì)量

m

），則有：

因為是均勻流，有

h

=

h

=

H

=

H

， 且

m

=

m

，

u

=

u

，重力分量與阻力匹配，簡化后有均勻流河床切應力公式：

式中，

H

為均勻流水深， τ為均勻流河床切應力。將式（11）和式（9）相結(jié)合，可以得到其他的切應力表達方式：當消除坡降因子

J

時，均勻流河床切應力公式表示為：

式中，

u

為均勻流速度。當消除水深因子

H

時 ，均勻流河床切應力公式表示為：

也可用單寬流量

H

·

u

表達均勻流河床切應力公式：

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

天然河流的流動遠沒有達到均勻流。如果當水深

H

=

H

時為均勻流， 那么當同樣流量下

H

>

H

時，或下游與均勻流同樣水位的條件下，降低進口流量，一定有一個恒定流的解沒有達到均勻流，需提高流量才能達到，所以命名為欠均勻流；同理，當

HH

為超均勻流，或均勻流條件下，保持下游水位不變，將上游流量增加，這樣的情況下也有一恒定流的解，因為流量高于均勻流，所以稱為超均勻流。以表1為模板，計算9個非均勻流的算例，下游水位保持不變，欠均勻流流量調(diào)至均勻流的1/4至1/2；超均勻流流量調(diào)至均勻流的2至5倍，統(tǒng)計結(jié)果見表3，并將其恒定流的解（進口水深

h

、出口水深

h

）一并列出。

表3 非均勻流各算例統(tǒng)計結(jié)果
Tab. 3 Statistical results of each example of non-uniform flow

流態(tài) 模型 算例 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個網(wǎng)格單元的流動數(shù)據(jù)進行積分，即應用如下表達式：

式中： ρ為水體密度，

g

為重力加速度，

h

、

h

分別為進、出口水流水深，

u

、

u

分別為進、出口的水流流速，

J

為河道比降，

l

為河道長度；等號左側(cè)前兩項分別為進口水壓力

P

、出口水壓力

P

，第3項為整個水體重力的水平分量

G

，第4項為河床阻力

F

，等號右側(cè)為動量改變量Δ

M

。除了切應力τ 之外，所有的量均可以準確統(tǒng)計出來，各算例具體受力見表4。

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

流態(tài) 模型 算例 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），計算得出的指數(shù)

b

的結(jié)果見表3，顯然除了第7算例外，其余8個算例均一致，取其平均數(shù)

b

≈2.335 5。

則非均勻流的河床切應力有如下的表達式：

當消除坡降因子

J

時，非均勻流河床切應力公式表達為：

式中，

u

為非均勻流流速。當用單寬流量

H

u

表達時，非均勻流河床切應力的計算公式為：

當用流速

u

表達時，非均勻流河床切應力的計算公式為：

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

h

h

，即“欠”均勻流狀態(tài)；第7算例為超均勻流邊界條件下的“超”均勻流，然而，現(xiàn)實中是否存在這樣的流動還需用試驗驗證。由表4還可以看出，非均勻流的重力分量

G

和阻力

F

差異較大，這也是導致壅水或跌水的根本原因。因此，式（17）～（19）對所有的直道河流流動均是普適的，尤其是式（17），只要知道流速

u

、水深

h

和糙率

n

，無論其坡降如何，便能求出其床面切應力。將式（12）和（16）代入式（17），確定

b

=7/3，消去

J

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

式（20）同均勻流床面切應力式（12）結(jié)構(gòu)相同，因此，無論流動是否為均勻流，式（20）均成立。通用的切應力公式可以將復雜的流動問題簡單化，將流體力學問題（求解偏微分方程）變?yōu)楹唵蔚拇鷶?shù)問題。

3 公式對比分析

許多學者都對非均勻流床面切應力進行過研究，也給出了一些經(jīng)驗公式。陳小芳等由試驗結(jié)果推導出相對河床切應力 τ/τ與相對水深

h

/

h

存在大致相同的經(jīng)驗公式：

金中武等也通過水槽試驗擬合出雍水條件下的經(jīng)驗公式：

張小峰等通過數(shù)模擬合出了經(jīng)驗公式：

為方便比較，將本文公式（16）改寫為如下形式：

將式（21）～（24）應用到表3的9個非均勻流的數(shù)據(jù)算例中，對比5個公式，驗證本文公式的合理性。對各公式切應力沿流動方向的100個單元進行積分，得出河床阻力見表5。由表4、5可以看出：當 τ >τ時，式（21）～（24）計算的結(jié)果比動量定理推導出的河床阻力偏??；反之，計算結(jié)果偏大。顯然，式（16）（即式（25））既符合動量定理又有物理意義。

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

N流態(tài) 模型 算例 式（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個水位，數(shù)據(jù)較少，并且由于非均勻流的水位波動較大，測量出來的水位對擬合結(jié)果也會有影響。張小峰等擬合出來的公式具有明確的物理意義，但并未用動量定理驗證。本文每個算例的每個流量都取了劃分出來的1 000個8節(jié)點四邊形網(wǎng)格中的直道中心處間隔相同的101個點位進行受力分析，每個點位都嚴格遵守動量方程，切應力同進出口水壓力、水體重力的水平分量及動量的改變?nèi)咧拖嗥ヅ?；每兩個點之間取1次切應力的平均值，最后再將所有點位的阻力求和，以近似于沿程積分出河床阻力，這將大大減少水位波動等不確定因素帶來的誤差影響，最終8個算例都得出相同的結(jié)論：相對河床切應力 τ/τ與相對水位（

h

/

h

）遵循式（25），因此，可以認為本文推導公式更合理且有物理意義。

4 結(jié) 論

本文對5種比尺的均勻流直道模型以及4種比尺不同流量下的非均勻流直道模型分別在均勻流和非均勻流情況下進行了數(shù)值模擬，對分析模擬出來的數(shù)據(jù)得出以下結(jié)論：

1）通過直道明渠均勻流的模擬數(shù)據(jù)結(jié)果，根據(jù)所有直道明渠的水位流量關(guān)系圖可以歸納出明渠均勻流公式：

Q

·

n

/

B

/

J

=

H

，并在此基礎上進行受力分析得出均勻流切應力表達式（消除坡降

J

）：τ=ρ

g

(

nu

)/

H

。2）假定明渠非均勻流切應力與非均勻流有同樣的形式 τ=(

u

/

u

)τ， 對直道明渠非均勻流的模擬數(shù)據(jù)結(jié)果進行分析，比例因子 (

u

/

u

)中指數(shù)

b

=3/7，最終得到消去

J

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

g

(

nu

)/

h

，該公式也同樣適用于垂直立面的切應力計算。

3）本文公式是基于動量定理及運用受力分析嚴格推導出來的，繪制出來的曲線與陳小芳、金中武等用試驗數(shù)據(jù)擬合出來的曲線相似，但當 τ>τ時，式（21）～（24）計算的結(jié)果比動量定理推導出的河床阻力偏??；反之，計算結(jié)果偏大。

4）由本文推導出的床面切應力公式可知，無論矩形明渠的流動是否均勻，均可通過水深

h

、流速

u

、糙率

n

這3個參數(shù)求出床面切應力的代數(shù)解，相比于已有非均勻流的切應力公式，本文推導出的通用切應力計算公式更加直接，且已在平面2維流動模型中應用。