墩柱壅水?dāng)?shù)值模擬研究

董耀華，楊元平

(1.中水東北勘測(cè)設(shè)計(jì)研究有限責(zé)任公司，長(zhǎng)春 130021； 2.浙江省水利河口研究院，杭州 310020)

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

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

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

1 物理水槽試驗(yàn)

試驗(yàn)在浙江省河口海岸重點(diǎn)實(shí)驗(yàn)室中進(jìn)行，見圖1。

圖1 水槽平面圖

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

圖2 水槽裝置圖

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

表1 物模試驗(yàn)水流條件

圖3 測(cè)點(diǎn)位置圖

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

圖4 壅水高度圖

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

2.1 控制方程與求解

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

(1)

(2)

(3)

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

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

2.2 網(wǎng)格劃分

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

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

2.3 模型參數(shù)

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

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

(4)

(5)

(6)

(7)

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

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

2.3.1 直接分配

按對(duì)應(yīng)計(jì)算水流條件，直接為計(jì)算域分配E值。

2.3.2 按Peclet數(shù)分配

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

(8)

2.3.3 按Smagorinsky方程分配

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

(9)

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

2.4 邊界條件

RMA2模型通過(guò)給定入口流量和出口水位設(shè)定邊界條件。流量邊界條件由物模試驗(yàn)墩柱位置流速乘以過(guò)水?dāng)嗝婷娣e(水深乘以水槽寬度)計(jì)算得到；水位邊界條件通過(guò)擬合物模試驗(yàn)墩柱位置的水深、流速率定得到。邊界條件取值見表2。

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

2.5 模型驗(yàn)證與建立

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

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

圖6 各組壅水曲線圖

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

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

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

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

表4 敏感性分析方案

3.1 直接分配E值

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

圖7 E值直接分配結(jié)果分析

ΔZx=α1E+β1

(10)

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

3.2 按Peclet數(shù)分配

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

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

(11)

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

3.3 按Smagorinsky方程分配

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

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

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

4 結(jié)論與展望

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

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