基于TELEMAC與SISYPHE耦合模型的圍堰附近砂床局部沖刷的三維數值模擬

譚成宇

(東臺市水利建設有限公司，江蘇 鹽城 224200)

圍堰是指為執行水利工程而修建的臨時性的圍護水工建筑物。它可以防止水土流入工程區域，從而使得受保護區域在施工期內的工程建設得以正常運行[1]。但是，圍堰的修建一般會縮小河道的斷面尺寸，提高局部流速，從而增大河床剪切應力。當河床的剪切應力超過河床的臨界切應力時，就會引發河床的局部沖刷，可能造成施工圍堰乃至整個工程的破壞[2]。因此，分析與預測施工圍堰造成的局部沖刷可以為圍堰的優化設計與安全評價提供必要的信息[3]，具有非常重要的現實意義。河床沖刷的預測分析一般可以采用經驗法、試驗法和數值模擬法。經驗法操作簡單，但是精確度往往難以滿足要求。試驗法較為精確，是當前最可靠的方法，但是，大尺度的試驗的時間與物力成本過高，而小尺度的物理模擬則可能會因為比例縮放問題而引發預測精確度的下降。隨著近年來計算機技術與資源的不斷發展，數值模擬法在水利工程領域得到了越來越廣泛的應用[4- 6]。TELEMAC是由法國電力集團開發的一款計算流體力學平臺[7]。它的一項重要優勢是免費且開源，因此，采用它建立模型可以降低成本、掌握源代碼以及為以后增減模型功能提供基礎。但是，標準的TELEMAC無法進行泥沙沖刷淤積的模擬。因此，本文提出將SISYPHE泥沙沖刷淤積模塊與TELEMAC進行耦合，并采用該耦合模塊完成圍堰周圍河床沖刷淤積三維數值模擬的方法。據筆者所知，基于TELEMAC與SISYPHE的圍堰局部沖刷三維數值模擬尚未被報道過，而其是否有能力準確模擬圍堰的局部沖刷亟需驗證。

因此，本文的主要目的是建立基于TELEMAC與SISYPHE耦合模型的圍堰附近泥沙沖刷模型并對其進性能進行評價。首先，進行圍堰局部沖刷的物理試驗，得到圍堰附近的沖刷實測數據，該數據不僅可用于本文數值模擬的驗證中，也可以為其它模型的驗證提供資料。其次，應用TELEMAC與SISYPHE耦合模型建立數值模型并根據試驗工況條件進行模擬。最后，將模擬結果與試驗值進行比較并對其精確度進行分析。結果表明，基于TELEMAC與SISYPHE耦合模型的圍堰附近泥沙沖刷模擬結果接近實測值，滿足工程要求，因此具有較好的推廣價值。

1 方法

1.1 試驗方法

不同的數值模型所得結果往往相差較大，因此，必須要對其進行驗證。本文首先進行物理試驗，以便于為數值模型的驗證與評價提供實測數據。試驗的基本設置如圖1所示。試驗水槽長度2.5m、寬度0.35m。圍堰沿x軸方向的長度為5cm、沿y軸方向的長度為15cm。水槽利用水泵進行水流循環，上游(圖1左側)為水流入口，流量為0.005m3/s，下游(圖1右側)為水流自由出口。水槽鋪有石英砂，其D50粒徑為0.001m。試驗主要關注超臨界流的沖刷過程，因此，當水流達到出口后則停止水循環，之后采用高度測量儀測量砂床的沖刷情況。根據初步試驗結果可知，本試驗中的主要泥沙沖刷淤積位于距離圍堰的y軸正方向邊緣(以下簡稱為“前邊緣”)約5cm的范圍以內，其它區域的砂床高程變化幾乎忽略不計，因此選取該范圍為研究區域。因此，研究區域的x坐標介于-5～10cm之間，y坐標介于-7.5～5cm之間。測量點分布間距為1cm，因此在研究范圍內共設置72個測量點，其分布如圖1中圈狀標記所示。

圖1 試驗設置與測量點分布圖

1.2 模擬方法

數值模擬在開源平臺TELEMAC中進行。將之與SISYPHE模塊進行耦合，從而可以進行泥沙沖刷淤積的計算。

水動力模擬可采用標準的Navier-Stokes方程，其表達式為：

(1)

(2)

(3)

(4)

式中，x、y、z—順水流方向、橫向和縱向的坐標；U、V、W—在x，y和z這3個方向上的速度分量；p—壓強；patm—大氣壓強；ν—黏性系數；ρ0—參考密度；Δρ—密度的變化值；Fx、Fy—其它的外力。

河床的高程變化、即泥沙的沖刷與淤積采用泥沙守恒公式[8- 9]，其表達式為：

(5)

式中，γ—砂床的孔隙率；zb—砂床的高程；t—時間；qB—近底泥沙通量；D—沉降率；E—與夾帶率。其中泥沙通量的計算采用經典的Van Rijn方程[10]計算。

水動力與泥沙沖刷淤積的耦合過程如圖2所示。首先，讀取初始的河床高程信息，并采用TELEMAC的標準水動力求解器計算流速、壓強等水力參數。之后，將求得的水力參數提供給SISYPHE的泥沙沖刷淤積模塊，該模塊首先也是讀取原始的河床信息，之后，求解泥沙守恒方程求得河床高程的變化。最后，新的河床高程信息被傳送回水動力模塊作為新的原始河床高程，并進行下一時間點的水動力模擬。

圖2 數值模型流程示意圖

2 結果與討論

2.1 試驗結果

本次試驗的結果呈現于圖3中，其中橫軸為x坐標，縱軸為砂床高程的變化值。比較圖3中各圖的河床高程變化量級，沿y軸負方向觀察，當斷面距離圍堰較遠時(如y=5cm)，河床沖刷深度約為0.5cm，當較為接近時(如y=2.5cm和y=0時)，河床沖刷深度接近1.5cm，沖刷深度逐漸增加。當斷面繼續沿y軸負方向漸遠時(如y=-2.5cm)，最大沖刷深度逐漸變小至1cm左右，當斷面遠離圍堰前邊緣時(如y=-2.5cm和y=-7.5cm時)，河床沖刷深度繼續變小，最終接近于0。沖刷深度在x軸方向的沖刷深度分布可以通過觀察各圖中河床高程變化值隨x值的變化情況進行分析，但y=-2.5cm和y=-7.5cm這兩個斷面的沖刷深度較小，且沿x軸的分布比較均勻，因此不作比較。對于其它各斷面，泥沙的沖刷深度大體上呈現由x=-5cm至x=0逐漸增大的趨勢。在x=5cm斷面的沖刷深度一般大于x=0斷面處的

圖3 實測與模擬河床高程變化值

深度。綜合上述分析，在y方向上，泥沙沖刷在圍堰前邊緣較為嚴重；在x方向上，泥沙沖刷在圍堰的上游邊緣較為嚴重。因此，圖1中橫坐標為x=0、縱坐標為y=0的位置可能出現最大的沖刷深度、穩定風險最高，之后沖刷深度沿各個方向隨距離的增大而逐漸變小。

2.2 模擬結果

圖3中實線代表模擬結果。對模擬曲線進行觀察可知，在x小于0的范圍內(即圍堰上游)，泥沙沖刷深度沿著橫軸正方向而逐漸增大，在x大于5cm的范圍內(即圍堰下游)，泥沙沖刷深度沿著橫軸正方向而逐漸減小。在y大于0的范圍內(即圍堰前邊緣側)，泥沙沖刷深度沿著y軸負方向而逐漸增大，在y小于0的范圍內(即圍堰后邊緣側)，泥沙沖刷深度沿著y軸負方向而逐漸減小。以上結論與實驗結果基本一致，因此可以認為該數字模型可以非常準確地模擬出泥沙沖刷深度的分布規律。

圖4 實測值與模擬值的散點比較圖

為進一步分析定量地判斷數值模型的模擬性能，做散點圖于圖4中，其中橫坐標代表試驗值，而縱坐標代表數值模型所提供的結果。圖中實線為等值線，當各散點距離其較遠時，則認為模擬結果偏差較大，反之則較小。當各散點坐落于等值線的上部時，則認為模擬值偏向于高估了淤積值(即低估沖刷值)，反之則有低估了淤積值(即高估沖刷值)的趨勢。根據圖4發現，各散點非常接近于等值線，且較為均勻地分布于等值線的兩側，因此可以認為其可以非常準確地預測圍堰的局部沖刷。圖4種的虛線為15%精度線，大部分散點位于兩條15%的范圍以內，即大部分模擬值的精確度均保持在15%以內。個別散點超過15%的范圍，但這幾個散點較接近于0，即泥沙沖刷深度較小，因此，誤差的絕對值也較小，而且，超出15%精度線的點主要分布于等值線下方，即高估了泥沙沖刷深度，偏向保守。因此，模擬結果證明該模型在施工圍堰的局部沖刷模擬中效果良好，滿足一般的工程要求。

2.3 誤差分析

為進一步了解該模型的模擬性能，對數值模擬的結果進行量化的誤差分析，所選用的指標包括平均絕對百分比誤差(MAPE)、平均相對誤差(MRE)、平均絕對誤差(MAE)、均方根誤差(RMSE)、決定系數(R2)，平均偏差(MBE)和標準化的均方根誤差(NRMSE)。各誤差指標的計算公式分別為：

(6)

(7)

(8)

(9)

(10)

(11)

(12)

式中，xcalc—模擬值，；xobs—實測值。

結果匯總于表1中，在y=-2.5cm，y=-5cm，和y=-7.5cm斷面均有沖刷深度值為0的數據，因此這幾個界面的MAPE、MRE和NRMSE值無意義、不列于表1中。根據表1中的MBE和MRE值可知，各斷面模擬值的平均偏差和平均相對誤差非常接近于0，因此整體上不偏向于高估或者低估。在絕對誤差方面，MAE和RMSE值的平均值為0.05cm和0.06cm，均較小，取其中較大值，可以認為該模擬的絕對誤差為0.06cm。在相對誤差方面，MAPE的平均值為7.34%，而NRMSE的平均值為0.1，即10%，均滿足誤差小于15%的一般工程要求，取其中較大值，可以認為該模擬的相對誤差在10%以內。

3 結語

通過進行圍堰附近局部河床沖刷的物理試驗，得到圍堰周邊沖刷點的高程變化數據，該數據不僅可用于本文數值模擬的驗證中，也可以為其它模型的驗證提供資料。基于TELEMAC與SISYPHE的耦合模型對該試驗進行了三維數值模擬，并且與試驗結果進行對比。模擬結果與實測的砂床變化高程數據具有較好的吻合度，模擬所得的沖刷深度分布規律與試驗基本一致，平均誤差值約為0.06cm、相對誤差約為10%。因此，證明了該方法在施工圍堰局部沖刷計算中的有效性，為提高當前的圍堰附近河床沖刷淤積的預測水平提供了有效的參考。本方法可以用于進行數值模擬試驗以分析不同的影響因素對砂床沖刷的影響，從而有助于設計合理的圍堰設計方案并進行有效的工程安全評，因此具有重要的現實意義。TELEMAC和SISYPHE均開源代碼，因此，在未來可以繼續對其進行改進，以進一步提高其模擬精確度。

表1 數值模擬誤差分析