王飛　張彬　齊劍峰

摘要：基于FLUENT軟件的動網(wǎng)格更新技術(shù)和用戶自定義函數(shù)功能實(shí)現(xiàn)了橋墩局部沖刷過程的三維動態(tài)模擬。以Melville經(jīng)典沖刷試驗(yàn)為原型，建立數(shù)值模型。將河床面設(shè)置為主要的動邊界，當(dāng)床面結(jié)點(diǎn)瞬時(shí)剪應(yīng)力大于臨界剪應(yīng)力時(shí)，結(jié)點(diǎn)位置下移，表現(xiàn)為沖刷，引入Van Rijn提出的沉積輸運(yùn)函數(shù)來控制河床面各結(jié)點(diǎn)的運(yùn)動速度。數(shù)值模擬結(jié)果在流場形態(tài)，沖坑發(fā)生發(fā)展過程及沖坑形態(tài)均與試驗(yàn)結(jié)果較為吻合，模擬的沖坑深度略小于試驗(yàn)結(jié)果，誤差約13%。誤差產(chǎn)生的主要原因?yàn)榛诶字Z平均N-S的湍流模型不能有效地反應(yīng)鈍形橋墩前端湍流脈動的影響。

關(guān)鍵詞：局部沖刷；動網(wǎng)格更新；數(shù)值模擬；橋墩；Melville沖刷試驗(yàn)

中圖分類號：U442 文獻(xiàn)標(biāo)識碼：A 文章編號：1672-1683（2017）02-0132-06

沖積河道在遇到阻水構(gòu)筑物（橋墩、丁壩等）時(shí)，構(gòu)筑物周圍河床的局部沖刷對阻水構(gòu)筑物的穩(wěn)定有很大的影響。對于橋墩而言，水流在遇到橋墩后，由于橋墩的阻水使得過水面積減小，墩周流速增大，河床剪應(yīng)力增加，墩周河床沉積物被水流搬運(yùn)，墩周床面高程逐步降低，并產(chǎn)生沖坑，導(dǎo)致橋墩基礎(chǔ)的埋深減小，進(jìn)而會導(dǎo)致橋梁的倒塌，甚至生命和財(cái)產(chǎn)的損失。

橋墩的局部沖刷是一個動態(tài)的發(fā)展過程，影響因素眾多，空間分布具有很強(qiáng)的三維特性，這就使得沖刷模型試驗(yàn)成為以往研究沖刷問題的主要手段。但模型實(shí)驗(yàn)存在費(fèi)用高，無法普遍應(yīng)用，條件單一，存在模型尺寸效應(yīng)等不確定因素，數(shù)值模擬方法的不斷改進(jìn)使得其作為一種研究手段越來越顯示出其不可替代的作用。

近年來，國內(nèi)外學(xué)者針對橋墩沖刷三維性態(tài)發(fā)展開展了一系列數(shù)值模擬研究。Ehteram運(yùn)用SSIIM軟件對橋臺的沖刷過程進(jìn)行了三維模擬，得到了沖刷坑深度和形狀并與試驗(yàn)結(jié)果進(jìn)行了比較。Khosronejad對不同橫截面形狀的橋墩進(jìn)行了三維動床模擬，采用了流固耦合曲線浸入邊界的技術(shù)。Kim采用大渦模擬的方法對相鄰的兩個圓柱形墩的局部沖刷坑進(jìn)行了模擬，得到的最大沖深位置與試驗(yàn)結(jié)果較為一致。韋雁機(jī)等基于Open-FOAM開源軟件的動網(wǎng)格技術(shù)，用輸沙率計(jì)算床面地形隨時(shí)間的變化，構(gòu)建起樁周局部沖刷的動態(tài)三維數(shù)學(xué)模型。祝志文等根據(jù)床底泥沙的單寬體積輸沙率得到河床高程坐標(biāo)的瞬時(shí)變化，采用邊界自適應(yīng)網(wǎng)格技術(shù)修改動邊界計(jì)算域網(wǎng)格，得到圓柱形橋墩周圍局部沖刷坑的演化過程。以上研究對局部沖刷的數(shù)值模擬起到了很好的推動作用，在實(shí)際應(yīng)用中多少都存在一些不足的地方，如采用虛擬的浸入邊界、地形函數(shù)等很難與實(shí)際條件一致，大渦模擬或分離渦的模擬計(jì)算消耗極大，基于單寬體積輸沙率來計(jì)算河床的變形，計(jì)算過程復(fù)雜，涉及到梯度的計(jì)算，會使誤差增加等。

本文基于CFD計(jì)算軟件FLUENT的動網(wǎng)格技術(shù)和用戶白定義函數(shù)（UDF）功能，使用基于雷諾平均N-S模型的Realizable k-e湍流模型，將床面瞬時(shí)剪應(yīng)力和臨界剪應(yīng)力帶入Van Rijn沉積輸運(yùn)函數(shù)，得到床面坐標(biāo)的變化，通過網(wǎng)格重構(gòu)和彈性光順結(jié)合的方法來不斷修正變形較大的網(wǎng)格，實(shí)現(xiàn)了局部沖刷過程的動態(tài)模擬。

1數(shù)值模型

1.1物理模型的選取

在已有的沖刷試驗(yàn)中，Melville和Raudkivi（MR）對局部沖刷坑發(fā)展的三個不同的階段進(jìn)行了相對較詳細(xì)的定量描述。所以本次研究選取MR的經(jīng)典沖刷試驗(yàn)資料建立數(shù)值模型，并進(jìn)行對比分析。MR試驗(yàn)水槽長19 m，寬0.456 m，水深0.15m，平均來流流速為0.25 m/s。模型布置見圖1。MR分別選取初始定床階段（測定了河床面附近的流速），中間發(fā)展階段（沖刷30 min時(shí)，沖坑深度達(dá)到0.04 m）和沖刷終止的平衡階段進(jìn)行分析，給出了詳細(xì)的試驗(yàn)結(jié)果。試驗(yàn)中前30 min發(fā)展較為劇烈，而之后沖刷發(fā)展開始緩慢。30 min時(shí)的沖刷深度達(dá)到總沖刷深度的75%。數(shù)值模擬取前30 min進(jìn)行研究。

1.2計(jì)算域及網(wǎng)格劃分

試驗(yàn)研究表明，圓柱形橋墩繞流流態(tài)以x軸基本呈對稱分布，因此取其中一半作為本次模擬的計(jì)算域，以縮短計(jì)算時(shí)間。根據(jù)試驗(yàn)布置和模擬要求，將計(jì)算域高度設(shè)置為15 cm，寬度設(shè)置為3d即15.24cm，圓柱上游距橋墩中心為3d，出口處要求尾流充分發(fā)展，所以設(shè)置下游距橋墩中心為10d^[s]即50.8 cm。取圓柱豎向?yàn)閦軸方向，床面為x-y面水流方向?yàn)閤軸正方向。具體見圖2。

本次模擬為動床模擬，橋墩周圍由于沖刷作用使得局部變形較大，且變形不規(guī)則，所以選用適應(yīng)變形能力較強(qiáng)的四面體非結(jié)構(gòu)化網(wǎng)格。為了提高床面剪應(yīng)力的模擬精度，在床面設(shè)置了0.2 cm（約0.5d50）邊界層，且邊界層會隨著床面結(jié)點(diǎn)的移動而跟隨移動，這樣更大程度上保證了床面剪應(yīng)力獲取的精度。因此本次模擬在床面邊界層內(nèi)為三棱柱體網(wǎng)格，其余部分為四面體網(wǎng)格。由于墩周及靠近床面的部位各物理量梯度較大，采用尺寸函數(shù)功能（最小網(wǎng)格尺寸0.3 cm，比率1.2，最大網(wǎng)格尺寸1.5cm）對局部網(wǎng)格進(jìn)行加密。動床面靠近墩周的部位考慮到湍流邊界層及后續(xù)發(fā)生較大的局部變形，采用尺寸函數(shù)（最小網(wǎng)格尺寸0.2 mm，比率1.05，最大網(wǎng)格尺寸1.5 cm）進(jìn)行加密。整個模型共劃分網(wǎng)格單元數(shù)127 224個，見圖2。

1.3湍流模型

針對湍流求解，大渦模擬（LES）和分離渦模擬（DES）方法對通過橋墩的大尺度渦的動力特性能夠精確的預(yù)測，但是由于其計(jì)算消耗過大，應(yīng)用到工程中有很大的挑戰(zhàn)。沖刷達(dá)到平衡的時(shí)間尺度（小時(shí)或天）比湍流脈動的時(shí)間尺度（秒或更?。┮蟮亩?，如此大的懸殊使得用LES和DES方法進(jìn)行沖刷的水動力耦合模擬不太實(shí)際。本次模擬采用的湍流模型為更經(jīng)濟(jì)實(shí)用的雷諾平均N-S模型。以往的研究表明，對于圓柱型墩，雷諾平均N-S模型的缺陷在于其不能有效捕獲上游面橋墩與河床相接處的高能湍流渦，而這樣的湍流脈動對于沖刷的發(fā)展是有影響的。因此可以預(yù)見，采用基于雷諾平均N-S模型的局部沖刷模擬在圓柱型橋墩前緣的河床沖刷深度會低于實(shí)際值。因此，本次采用動床模擬沖坑深度，可以比較得出基于雷諾平均N-S模型的數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果的誤差大小并分析誤差的來源。另外，橋墩前緣的湍流脈動，對于橋墩前緣為非圓柱形狀的情況，如尖角型，可能會產(chǎn)生不同的流動類型和沖刷動力。所以，橋墩形狀，尤其是前緣的形狀對數(shù)值模擬結(jié)果的精度在文末進(jìn)行了探討。

1.4邊界條件

由于計(jì)算域選取的流場入口段距離較小，所以要經(jīng)過計(jì)算給定一個穩(wěn)定的，邊界層充分發(fā)展的流速剖面，作為速度入口邊界條件。因此，在三維計(jì)算之前，首先建立二維無圓柱流場并給定速度入口條件讓其充分發(fā)展，模擬結(jié)果和Melville試驗(yàn)結(jié)果均顯示出充分發(fā)展的速度剖面分布基本符合最廣泛使用的karman-Prandtl對數(shù)流速分布公式，即：

（1）

底部河床指定為粗糙壁面，根據(jù)Melville試驗(yàn)結(jié)果，其有效粗糙高度取為2d50。橋墩面，水槽側(cè)壁均設(shè)置為光滑壁面。頂面設(shè)置為對稱邊界來模擬自由水面。由于取一半流場進(jìn)行數(shù)值模擬，沿x軸剖分出來的面均設(shè)置為對稱邊界。具體邊界設(shè)置見圖3。

1.5動態(tài)網(wǎng)格更新

CFD模擬中，流場的計(jì)算采用單相，瞬態(tài)求解。在FLUENT中激活動網(wǎng)格，河床面設(shè)置為動邊界，用DEFINE_GRID_MOTION宏命令來控制邊界各個節(jié)點(diǎn)的運(yùn)動。在每個時(shí)間步開始計(jì)算時(shí)，比較床面（動邊界）各結(jié)點(diǎn)實(shí)時(shí)剪應(yīng)力T與床沙起動臨界剪應(yīng)力值τ^σ，若存在超臨界剪應(yīng)力（τ-τ^cr>0），該點(diǎn)表現(xiàn)為沖刷，結(jié)點(diǎn)下移，否則表現(xiàn)為靜止，結(jié)點(diǎn)位置不變。局部沖刷問題，局部網(wǎng)格變形較大且不太規(guī)則，隨著沖刷坑逐漸發(fā)展，局部網(wǎng)格必然變大或扭曲造成數(shù)值發(fā)散，所以為了保證各區(qū)域變形后網(wǎng)格尺寸不至于過大或者過于扭曲而使數(shù)值發(fā)散，選用局部網(wǎng)格重構(gòu)與彈性光順相結(jié)合的方式進(jìn)行動態(tài)網(wǎng)格更新，網(wǎng)格變形前后對比圖見圖4。

在FLUENT軟件中，通過用戶白定義函數(shù)（UDF）獲取床面實(shí)時(shí)剪應(yīng)力值并存儲。平床下床沙起動的臨界剪應(yīng)力值通過由希爾茲公式推導(dǎo)得出的臨界剪應(yīng)力的計(jì)算方法：

（2）

（3）

隨著沖刷的發(fā)展，床面開始出現(xiàn)坡度并逐漸增加，由于床沙重力在水流方向產(chǎn)生分力，床沙的起動臨界剪應(yīng)力將不再等于平床下的臨界剪應(yīng)力值。本文采用Dey提出的經(jīng)驗(yàn)方程來計(jì)算變臨界剪應(yīng)力：

（4）

除了通過剪應(yīng)力確定各結(jié)點(diǎn)在每個時(shí)間步是否沖刷產(chǎn)生向下位移外，還需要確定各結(jié)點(diǎn)向下位移的大小。在每個時(shí)間步內(nèi)，各結(jié)點(diǎn)向下的位移等于時(shí)間步長與網(wǎng)格移動速度的積。通過調(diào)整時(shí)間步長可以控制在每一時(shí)間步內(nèi)網(wǎng)格移動增量處于一個較小的值。底邊界網(wǎng)格移動的速度根據(jù)Van Rijn基于水槽試驗(yàn)提出的沉積輸運(yùn)函數(shù)來表示：

（5）

2模型驗(yàn)證

2.1流場對比

在進(jìn)行動床模擬前先進(jìn)行了定床條件下的流場模擬。從圖5可以看出，墩前流線的分離和墩后尾渦的形態(tài)，墩前垂直剖分面上的下降水流都基本與試驗(yàn)結(jié)果吻合。Melville試驗(yàn)定床沖刷時(shí)，觀測到在橋墩附近流速明顯變大，最大流速在迎流面中軸線兩側(cè)±100°的位置，數(shù)值約為1.5倍的平均來流速度即0.37 m/s。圖6是數(shù)值模擬得出的流速等值線分布圖和床面剪應(yīng)力分布圖，無論是最大流速的位置還是大小都與試驗(yàn)結(jié)果基本一致。最大河床剪應(yīng)力的位置與最大流速的位置基本一致，這也與試驗(yàn)結(jié)果基本吻合，說明流速的增加引起了床面剪應(yīng)力的增加，當(dāng)床面剪應(yīng)力超過床沙起動應(yīng)力時(shí)就表現(xiàn)為沖刷。

2.2沖刷過程對比

試驗(yàn)和數(shù)值模擬結(jié)果均表明橋墩的局部沖刷經(jīng)歷了一個先強(qiáng)后弱的過程（圖7）。在沖刷開始階段，試驗(yàn)和模擬結(jié)果高度吻合，沖刷速率很大。當(dāng)沖刷約7 min后，沖刷速率減緩，試驗(yàn)深度開始大于數(shù)值模擬的深度。在30 min時(shí)，試驗(yàn)深度達(dá)到4 cm，數(shù)值模擬深度為3.5 cm，數(shù)值模擬最大沖深比試驗(yàn)結(jié)果小了13%。

從Melville試驗(yàn)及其他研究者的研究成果可以看出，沖刷首先出現(xiàn)在橋墩兩側(cè)（迎流面中軸線兩側(cè)45°～100°之間），沖刷的產(chǎn)生是由于橋墩兩側(cè)水流的加速引起了床面剪應(yīng)力超過床沙起動剪應(yīng)力，從而產(chǎn)生了沖刷。從圖8看出數(shù)值模擬沖坑也首先出現(xiàn)在橋墩兩側(cè)45°～100°之間，與試驗(yàn)結(jié)果吻合。沖刷幾分鐘后，橋墩前端下降水流引起的湍流渦系開始顯現(xiàn)并逐漸發(fā)揮對床沙沖刷的加速作用，沖坑開始由橋墩兩側(cè)往橋墩前端貫通，最終沖坑最深的位置處于橋墩迎流面正前方到兩側(cè)75°范圍內(nèi)，如圖9試驗(yàn)沖坑等值線所示。數(shù)值模擬沖坑等值線整體形態(tài)與試驗(yàn)結(jié)果相近，主要區(qū)別在于橋墩前端沖深比試驗(yàn)結(jié)果小。正如1.3節(jié)所述，數(shù)值模擬采用的雷諾平均N-S模型對圓柱型墩正前端湍流脈動對沖刷的貢獻(xiàn)不能有效的模擬出來，使得橋墩正前端的沖坑深度低于試驗(yàn)結(jié)果。

從以上分析可以看出，本文所采用的動態(tài)網(wǎng)格模擬沖坑的方法能夠較好地反應(yīng)沖刷發(fā)展的動態(tài)過程，模擬結(jié)論可信。

3橋墩橫截面形狀影響的探討

從前面的分析可以看出，對于圓柱型的橋墩，產(chǎn)生沖刷的原因主要來源于以下兩個方面：（1）由于橋墩的存在，局部流速加大，引起局部河床剪應(yīng)力超過了臨界剪應(yīng)力，引起墩周局部沖刷；（2）湍流脈動速度引起的瞬時(shí)的河床剪應(yīng)力增大引起的局部沖刷。本文所采用的數(shù)值模擬方法是基于雷諾平均N-S模型的，因此對于湍流脈動速度估計(jì)不足，從而使得墩頭處模擬深度與實(shí)際相比偏小。以往的研究表明，墩頭鈍度不同，湍流脈動強(qiáng)度也不同。比如方形橋墩，鈍度最大，圓端型次之，尖角型墩鈍度因子（BF）為0。Olcmen和Simpson經(jīng)過研究發(fā)現(xiàn)，湍流脈動和馬蹄渦強(qiáng)度隨障礙物前緣鈍度的減小而逐漸降低，當(dāng)鈍度因子為O時(shí)，湍流脈動微弱到無法測得。因此本文所采用的數(shù)值模擬方法，對于尖角型墩可以達(dá)到準(zhǔn)確預(yù)測，隨著鈍度的增大，此方法在墩頭處的預(yù)測深度會比實(shí)際的偏小，已有的試驗(yàn)研究成果也已觀測到，對于尖角型墩，其最大沖刷主要位于兩側(cè)角部而墩前端的沖刷深度很小。而圓端型和方型的墩最大沖深均處于墩前端及附近一個較小的范圍內(nèi)，這也間接地證明了以上分析的合理性。具體更詳細(xì)的驗(yàn)證有待后續(xù)進(jìn)一步的研究。

4結(jié)論

床面的沖刷類似于柔性變形，因此河床底面坐標(biāo)的改變要對河床面的各個結(jié)點(diǎn)進(jìn)行控制。各個結(jié)點(diǎn)在每個時(shí)間步是否下移和下移多少主要取決于兩個關(guān)鍵參數(shù)：床面剪應(yīng)力和臨界剪應(yīng)力。通過對模型底部網(wǎng)格設(shè)置邊界層來提高床面剪應(yīng)力的提取精度。臨界起動剪應(yīng)力會隨著床面坡度的改變而變化，所以采用變臨界剪應(yīng)力的方法更符合實(shí)際。各結(jié)點(diǎn)運(yùn)動的速度通過Van Rijn沉積輸運(yùn)函數(shù)來獲得運(yùn)動的速度，速度乘以時(shí)間步長就是在每個時(shí)間步內(nèi)坐標(biāo)移動的大小。與求解輸沙率來求得床面坐標(biāo)變化的方法相比，本方法不需求解標(biāo)量的梯度和微分方程，計(jì)算簡單，避免了梯度求解產(chǎn)生的誤差，沖刷的速率更接近試驗(yàn)結(jié)果。本方法在應(yīng)用時(shí)，鈍形前端橋墩的最大沖深比實(shí)際情況會小約13%，本方法是在FLUENT軟件基礎(chǔ)上的二次開發(fā)，可塑性較高，在后續(xù)的研究中可通過嘗試增強(qiáng)橋墩前端床面剪應(yīng)力的方法來減少鈍形橋墩前端沖深預(yù)測不足的缺點(diǎn)。