董耀華，楊元平

(1.中水東北勘測設計研究有限責任公司，長春 130021； 2.浙江省水利河口研究院，杭州 310020)

墩柱壅水[1]計算是橋梁設計的重要內(nèi)容，也是防洪影響評價的核心任務[2-3]。在橋梁水力計算中，除天然水流與自然河床外，總是以墩柱壅水計算作為主要內(nèi)容[4]，因此墩柱壅水計算有著重要的實際意義。

近些年來，數(shù)值模擬逐漸成為墩柱壅水問題研究的重要手段[5-10]。而在數(shù)值模擬中，計算參數(shù)對于模擬準確性有重要影響，特別是紊動黏性系數(shù)(eddy viscosity)與壅水高度密切相關[11]，但對紊動黏性系數(shù)的作用規(guī)律研究仍然不足。目前，在利用數(shù)值模擬分析橋墩壅水的研究中，通常選取率定的方式確定紊動黏性系數(shù)的取值[5,12-13]。該方法依賴實測壅水資料，當缺乏壅水資料時，模型的準確性也就難以保證。張瑋[14]等探討了紊動黏性系數(shù)3種取值方式對壅水計算的影響，結果表明3種方法均可以成功地進行壅水計算。

本文擬設計物理模型試驗，并利用試驗數(shù)據(jù)建立平面二維水動力數(shù)學模型，然后應用模型分析紊動黏性系數(shù)3種取值方法與壅水高度的敏感性關系，研究3種紊動黏性系數(shù)取值方法與壅水高度的各自的作用規(guī)律，研究成果可為墩柱壅水數(shù)值計算提供一定參考。

1 物理水槽試驗

試驗在浙江省河口海岸重點實驗室中進行，見圖1。

圖1 水槽平面圖

水槽為長50 m、寬4 m、高0.5 m的矩形平底水槽，水槽兩端各有一組水泵，通過水流控制室的軟件控制水泵的頻率，可以調(diào)節(jié)水槽中水流的水位與流速；水槽中心位置放置直徑11 cm圓柱形墩柱；水位測量采用超聲波探頭，精度為0.1 mm，探頭固定在可向三維方向移動的高精度機床上，通過軟件輸入坐標，即可將探頭移動至所需位置測量，水位測量裝置及試驗水槽局部見圖2。此外，通過在墩柱上包裹細網(wǎng)格紙、利用相機拍攝墩柱淹沒位置來測量墩柱前后水位差的方式，輔助觀測壅水高度。

圖2 水槽裝置圖

本試驗共設計3組工況，見表1，均為恒定流。表1中，流速、水深為樁柱位置處的平均流速與水深。為獲取壅水數(shù)據(jù)，每組工況均進行有墩柱、無墩柱放置兩組試驗。在圖3中，以墩柱中心為坐標原點，對X軸上8個測點進行數(shù)據(jù)采集，通過計算有無墩柱試驗數(shù)值差得到壅水數(shù)值；探頭數(shù)據(jù)采集頻率為20 Hz，考慮到水流的波動性，探頭采集時間為1 min，取1 min數(shù)據(jù)平均值為采集結果值。

表1 物模試驗水流條件

圖3 測點位置圖

試驗結果見圖4。考慮到測量誤差及水位波動，壅水值有一定隨機性，特別是遠離墩柱位置處壅水值偏小的測點，波動較為明顯，但整體變化規(guī)律仍十分清晰，可以作為數(shù)學模型驗證的依據(jù)。

圖4 壅水高度圖

2 平面二維水流數(shù)學模型

2.1 控制方程與求解

應用SMS軟件RMA2模塊建立二維水流數(shù)學模型，RMA2模塊的控制方程是將三維流動基本方程沿水深積分后平均得到如下方程組[6]。其中，式(1)為水流連續(xù)性方程，式(2)、式(3)分別為X、Y方向上的動量守恒方程。

(1)

(2)

(3)

式中：h為水深；u、v為X、Y方向的流速；x、y、t為直線坐標系和時間；ρ為流體密度；E為紊動黏性系數(shù)，xx為X軸表明上的法線方向；g為重力加速度；α為底部標高；n為曼寧系數(shù)；ζ為風應力系數(shù)；vα為風速；Ψ為風向；w為地球角速度；φ為當?shù)鼐暥取?/p>

模型的求解過程為：①通過導入地形信息文件或手動輸入散點，生成計算網(wǎng)格；②根據(jù)地形復雜程度和流速梯度的大小，對網(wǎng)格的疏密性及部分節(jié)點進行調(diào)整；③將模型計算網(wǎng)格和控制條件導入RMA2計算模塊，對控制方程進行時間與空間上的離散，時間離散采用差分法，空間離散采用有限元法，并采用伽遼金加權余量法把控制方程從偏微分方程組轉變成代數(shù)方程組進行求解計算；④計算結果后處理，輸出圖文成果[15]。

2.2 網(wǎng)格劃分

參照物模試驗建立數(shù)學模型，水槽長50 m、寬4 m，圓柱形墩柱直徑為0.11 m。墩柱位于水槽中心，墩柱及附近區(qū)域用三角形網(wǎng)格劃分，其它區(qū)域按四邊形網(wǎng)格劃分，將墩柱作為不透水區(qū)域處理，網(wǎng)格單元數(shù)量為18 212，單元節(jié)點數(shù)量為40 973，最大網(wǎng)格尺度0.35 m，最小網(wǎng)格尺度為0.01 m。網(wǎng)格劃分見圖5。

圖5 網(wǎng)格劃分情況

2.3 模型參數(shù)

曼寧糙率n與紊動黏性系數(shù)E取值對計算結果有重要影響[11]。糙率與水力坡度相關，根據(jù)實測比降率定得到糙率n為0.013。

紊動黏性系數(shù)E與壅水高度值有關。紊動黏性系數(shù)E定義為影響湍流交換的控制參數(shù)，而湍流(Turbulence)通?？梢远x為速度隨時間變化的影響，以及與其空間梯度相關的動量交換，是流體作不規(guī)則運動的一種流動狀態(tài)。紊動黏性系數(shù)E的作用方程為：

(4)

(5)

(6)

(7)

式中：μ為分子黏度(Molecular viscosity)；u′、v′分別為X、Y方向的湍流速度波動。

紊動黏性系數(shù)E包括分子黏性力項和湍流速度波動項，但后者比前者大幾個數(shù)量級，因此分子黏度通常被忽略。有3種基本方法可控制紊動黏性系數(shù)E，具體如下：

2.3.1 直接分配

按對應計算水流條件，直接為計算域分配E值。

2.3.2 按Peclet數(shù)分配

該方法允許模型在每次迭代后根據(jù)提供的Peclet數(shù)即P值自動調(diào)整E值。P值是基于每個元素內(nèi)的唯一大小和計算出的速度，P值定義了平均元素速度值、元素長度、流體密度和E之間的關系，計算公式如下：

(8)

2.3.3 按Smagorinsky方程分配

該方法可以根據(jù)計算出的速度實時調(diào)整紊動黏性系數(shù)E值。相對于按Peclet數(shù)分配，它考慮了速度梯度，以確定合適的湍流系數(shù)來滿足流體動力學模擬中的條件。Smagorinsky方程如下：

(9)

式中：TBFACT為方程系數(shù)；A為單元的面積；?u/?x和?v/?y為速度分量的變化率；E為紊動黏性系數(shù)。

2.4 邊界條件

RMA2模型通過給定入口流量和出口水位設定邊界條件。流量邊界條件由物模試驗墩柱位置流速乘以過水斷面面積(水深乘以水槽寬度)計算得到；水位邊界條件通過擬合物模試驗墩柱位置的水深、流速率定得到。邊界條件取值見表2。

表2 數(shù)模試驗水流條件

2.5 模型驗證與建立

確定數(shù)學模型的網(wǎng)格樣式、邊界條件設置、參數(shù)糙率值后，將紊動黏性系數(shù)以默認值設置，建立數(shù)學模型，并通過物理模型試驗的壅水數(shù)據(jù)率定修正紊動黏性系數(shù)取值，最終數(shù)值模擬計算結果與對應紊動黏性系數(shù)取值見表3與圖6。

表3 紊動黏性系數(shù)取值

圖6 各組壅水曲線圖

由于物模試驗墩后測點出現(xiàn)局部跌水,導致在流態(tài)上與數(shù)學模型存在差異[13],數(shù)值模擬時墩后水面下降模擬效果并不理想，故主要以墩前壅水驗證紊動黏性系數(shù)取值。分析模擬結果，工況C1-工況C3的3種取值方法中，數(shù)學模型模擬壅水曲線與實測物模試驗數(shù)據(jù)一致，數(shù)值模擬結果較為準確，可以認為表3的參數(shù)選取較為合理。此外，在靠近墩柱迎流測，按Smagorinsky方程分配的數(shù)模壅水曲線較另外兩種方法結果低，水面線形態(tài)也與物模試驗觀察到的水面線形態(tài)接近。

3 紊動黏性系數(shù)敏感性分析

墩柱壅水數(shù)值模擬中，紊動黏性系數(shù)取值與壅水高度值密切相關。本節(jié)分析紊動黏性系數(shù)3種取值方法與壅水高度的敏感性關系，為紊動黏性系數(shù)的取值提供參考。

取工況C2進行數(shù)值模擬，在第2節(jié)數(shù)模結果基礎上設計試驗方案，分別研究3種紊動黏性系數(shù)方法取值與壅水高度的影響關系，方案設計見表4。其中，M8組方案為第二節(jié)工況C2數(shù)模紊動黏性系數(shù)取值。

表4 敏感性分析方案

3.1 直接分配E值

圖7(a)為數(shù)模計算結果生成壅水高度曲線圖，圖7(b)為數(shù)模計算生成的紊動黏性系數(shù)取不同值時的墩前各位置的壅水高度數(shù)據(jù)。由圖7(a)可知，直接分配法中，壅水高度隨E值增大而增大。由圖7(b)中數(shù)據(jù)擬合可以得到E值與壅水高度的關系式(10)，并據(jù)此推出墩前各位置的壅水高度與E值之間具有明顯的線性關系，即E值增加，壅水高度也按固定比例增加。

圖7 E值直接分配結果分析

ΔZx=α1E+β1

(10)

式中：ΔZx為橫坐標x處的壅水高度；α1、β1均為待定參數(shù)。

3.2 按Peclet數(shù)分配

圖8(a)為數(shù)模計算結果生成壅水高度曲線圖，圖8(b)為數(shù)模計算生成的紊動黏性系數(shù)取不同值時的墩前各位置的壅水高度數(shù)據(jù)。由圖8(a)可知，按Peclet數(shù)分配法中，壅水高度隨P值增大而減小。由圖8(b)中數(shù)據(jù)擬合可以得到P值與壅水高度的關系式(11)，并據(jù)此推出墩前各位置的壅水高度與P值間具有明顯的乘冪關系，壅水高度隨P值的增加而減小。當P趨近于0時，壅水高度大幅增加；隨著P值增大，壅水高度的減小幅度逐漸趨緩。

圖8 按Peclet數(shù)分配結果分析

(11)

式中：α2、β2均為待定參數(shù)。

3.3 按Smagorinsky方程分配

圖9為數(shù)模計算結果生成壅水高度曲線圖，表5為數(shù)模計算生成的紊動黏性系數(shù)取不同值時墩前各位置的壅水高度數(shù)據(jù)。由圖9可知，按Smagorinsky方程分配中，壅水高度隨Smagorinsky系數(shù)增大而增大。為進一步分析表5，計算表5中相鄰系數(shù)的壅水高度差值可知，壅水高度隨Smagorinsky系數(shù)的增大而增加，但壅水高度的增加幅度隨Smagorinsky系數(shù)增大逐漸減小。

圖9 按Smagorinsky方程分配壅水曲線圖

表5 Smagorinsky系數(shù)取不同值時中軸線各處的壅水高度 /cm

4 結論與展望

本文設計了墩柱壅水物理模型試驗實測樁柱壅水，建立了描述墩柱壅水的平面二維水流模型，并應用該模型分析紊動黏性系數(shù)3種取值方法取值對壅水高度的影響。主要結論如下：紊動黏性系數(shù)的3種取值方式中，直接分配法中E值與壅水高度成線性關系，E值增加，壅水值以固定比例增加；按Peclet數(shù)分配法中P值與壅水高度成乘冪關系，壅水高度隨P值的增加而減小，壅水高度的增加幅度隨Smagorinsky系數(shù)增大逐漸減?。话碨magorinsky系數(shù)分配法中系數(shù)TBFACT與壅水高度無明顯的函數(shù)關系，壅水高度隨Smagorinsky系數(shù)TBFACT的增大而增加，但壅水高度的增加幅度隨Smagorinsky系數(shù)TBFACT增大逐漸減小。

本文給出了紊動黏性系數(shù)3種取值方法與壅水高度的影響關系，在實際工作中，可以為紊動黏性系數(shù)的取值提供一定參考。