橋墩局部沖刷發展過程的三維動網格模擬

王飛　張彬　齊劍峰

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

關鍵詞：局部沖刷；動網格更新；數值模擬；橋墩；Melville沖刷試驗

中圖分類號：U442 文獻標識碼：A 文章編號：1672-1683（2017）02-0132-06

沖積河道在遇到阻水構筑物（橋墩、丁壩等）時，構筑物周圍河床的局部沖刷對阻水構筑物的穩定有很大的影響。對于橋墩而言，水流在遇到橋墩后，由于橋墩的阻水使得過水面積減小，墩周流速增大，河床剪應力增加，墩周河床沉積物被水流搬運，墩周床面高程逐步降低，并產生沖坑，導致橋墩基礎的埋深減小，進而會導致橋梁的倒塌，甚至生命和財產的損失。

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

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

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

1數值模型

1.1物理模型的選取

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

1.2計算域及網格劃分

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

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

1.3湍流模型

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

1.4邊界條件

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

（1）

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

1.5動態網格更新

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

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

（2）

（3）

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

（4）

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

（5）

2模型驗證

2.1流場對比

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

2.2沖刷過程對比

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

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

從以上分析可以看出，本文所采用的動態網格模擬沖坑的方法能夠較好地反應沖刷發展的動態過程，模擬結論可信。

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

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

4結論

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