墩柱壅水數值模擬研究

董耀華，楊元平

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

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

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

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

1 物理水槽試驗

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

圖1 水槽平面圖

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

圖2 水槽裝置圖

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

表1 物模試驗水流條件

圖3 測點位置圖

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

圖4 壅水高度圖

2 平面二維水流數學模型

2.1 控制方程與求解

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

(1)

(2)

(3)

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

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

2.2 網格劃分

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

圖5 網格劃分情況

2.3 模型參數

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

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

(4)

(5)

(6)

(7)

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

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

2.3.1 直接分配

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

2.3.2 按Peclet數分配

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

(8)

2.3.3 按Smagorinsky方程分配

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

(9)

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

2.4 邊界條件

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

表2 數模試驗水流條件

2.5 模型驗證與建立

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

表3 紊動黏性系數取值

圖6 各組壅水曲線圖

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

3 紊動黏性系數敏感性分析

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

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

表4 敏感性分析方案

3.1 直接分配E值

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

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

ΔZx=α1E+β1

(10)

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

3.2 按Peclet數分配

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

圖8 按Peclet數分配結果分析

(11)

式中：α2、β2均為待定參數。

3.3 按Smagorinsky方程分配

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

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

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

4 結論與展望

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

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