管道內泥沙粗顆粒沖淤過程的數值模擬

謝潤瑜 韓迅*

（四川大學水力學與山區河流開發保護國家重點實驗室，四川 成都 610065）

1 概述

在城市排水系統中，由于污水所挾帶的固體顆粒比重比水大，同時管道中的水流特性存在間歇性，使得排水管道內經常出現了沉積物淤積的現象，大量的淤積會導致管道損壞[1]。當排水系統出現淤積甚至堵塞時，可以使用人工清淤或者水力沖洗的方式進行疏通清理。國內大多數城市依舊使用人工清理的方式，而歐洲許多國家已經針對管道水力自清理方法進行了比較多的研究，也已經獲得了許多專利，如法國的水力平衡閥、德國和瑞士的水力自凈系統[2]等等。因此，對于排水管道內沉積物的沖淤研究是十分有意義的，本文將使用CFD方法對管道內泥沙沖淤進行三維數值研究。

由于真實泥沙粒徑小，實驗測量技術和設備要求高，數值模擬逐漸成為當前世界研究該問題的有效方式。數值模擬逐漸成為當前世界研究該問題的有效方式，模擬方式可分為數學模型和形態學模型。數學模型通常是建立在對流擴散方程基礎上，描述泥沙的質量輸移。形態學模型則可描述泥沙顆粒的具體運動形態，有無耦合、半耦合或完全耦合三種形式。伴隨著高性能計算機技術的發展，研究者們逐漸擺脫了“利用濃度表示顆粒分布，通過求解泥沙連續性方程研究泥沙輸移”的數值模擬方法，轉而以單個泥沙顆粒為研究對象，利用顆粒離散元法求解數以萬計的顆粒運動，得到的結果也更為真實[3]。

本文采用三維流固耦合CFD代碼模型（CgLES_Ycode）對圓管內泥沙粗顆粒運動進行數值模擬研究，采用浸沒邊界法和顆粒離散元法來模擬管道內水流和泥沙顆粒之間的耦合作用。首先，本文驗證了靜水條件下單個顆粒自由沉降運動以及明渠紊流中泥沙顆粒的輸移；接著，基于浸沒邊界法進行了圓管內水流流動特性研究；最后，在不同流速條件下對管段內均勻顆粒沙床沖刷進行數值研究。

2 數值方法

2.1 浸沒邊界法。本文粘性不可壓流體的控制方程為Navier-Stoke 方程。采用具有二階時間精度的Adams—Bashfort格式對其進行時間離散，可得控制方程的守恒形式為：

式中，Un+1為將歐拉點插值到拉格朗日點上的流體速度，表示為Un+1=I（un+1）；而Vn+1為拉格朗日點上的固體速度。這里引入插值函數I（φ）和分布函數D（φ）,具體描述參考文獻[4]。將式(1)代入式(3),經過整理，可得到定義在歐拉點上的附加體積力為：

式中，Xp和Vp分別為顆粒的位移矢量和速度矢量，Fp為作用在顆粒上的合力矢量，m 為顆粒的質量。作用在顆粒上的力包括有效重力FG，水流拖曳力FD和升力FL，法向接觸力Fn和切向摩擦力Fs。由于模擬泥沙粒徑較小，模型不考慮顆粒轉動。Fp的表達式如下：

式中，r 為顆粒與顆粒（壁面）之間間距，（xi，yi，zi）和（xj，yj，xj）分別為第i、j 個顆粒（或壁面）的坐標。

3 驗證算例

3.1 單個顆粒的自由沉降驗證。首先，本文進行單個顆粒在靜水中沉降運動的模擬以驗證單個顆粒與流體之間的耦合作用正確性，在重力方向上顆粒受到的流體拖曳力與有效重力平衡時，顆粒達到最終沉降速度。在驗證算例中，顆粒直徑為0.1cm，長方體流域高度和寬度分別為6.4cm和0.4cm，顆粒初始高度為6.3cm，網格尺寸為粒徑的1/4，重力加速度為980cm/s2，顆粒與流體密度分別為2.6g/cm3和1.0g/cm3，流體運動粘度為0.01cm2/s。經過對模擬數據的處理，可得到如圖1 的數據結果。

圖1 給出了使用該代碼程序模擬得到的顆粒沉降速度V和位移X歷時曲線，并且與其他學者所模擬的數值結果進行了比較，結果擬合比較良好。此外，本研究得到的最終沉降速度為14.6963cm/s，與2021 年陸程等人報導的數值結果(14.7031cm/s)以及2017 年宋等人的試驗結論(15.1392cm/s)差距不大。

圖1 顆粒自由沉降剖面速度云圖、三維示意圖以及速度和位移歷時曲線

3.2 明渠紊流中泥沙顆粒的輸移驗證。本文通過模擬明渠紊流中泥沙顆粒運動來考察顆粒輸移情況，計算域分為順流向、壁面法向和展向，在三個方向上的大小為6.0×4.0×1.0，采用64×64×64 的均勻網格進行空間離散。數值模擬時，順流向和展向均采用周期性邊界條件，底面和頂面分別采用不可滑移和自由滑移邊界條件，基于高度和摩阻流速的雷諾數Res=1000。顆粒直徑為1.0，顆粒與流體密度比為2.65。

本文通過調整顆粒重力加速度以獲得一定范圍的希爾茲參數，同時水流條件保持不變，這樣能夠保證每個工況的顆粒雷諾數相同，使得模擬結果更具有可比性。四種不同希爾茲參數的沙床輸移示意圖和無量綱輸沙率隨希爾茲數變化的散點圖如圖2 所示。

圖2 沙床輸移示意圖、無量綱輸沙率隨希爾茲數變化的散點圖

圖2 給出的模擬散點圖與經典的和修正的Meyer-Peter 和Muller 公式進行了比較，可以看到隨著希爾茲參數的增大，無量綱輸沙率逐漸遞增，并且結果均介于兩個經驗曲線之間，吻合十分良好。

3.3 圓管內水流流動特性研究。為了能夠進一步實現圓管內水流條件下的顆粒運動，本文采用浸沒邊界法進行圓管層流的數值模擬，確定不同雷諾數時沿程阻力系數的數值結果與理論計算結果間誤差，對拋物線速度入口邊界條件下的流體速度分布進行與理論曲線的比較，以驗證內部流場隨雷諾數變化的情況是否正確。

當雷諾數小于2000 時，流體在圓管中的流動被視為層流流動，因此選取雷諾數Re 介于50 到2000 之間的8 個數值進行層流模擬。通過三維建模獲得內徑為D=1.2cm、長度為L=12.0cm 的圓形直管有限元模型，將其浸沒于一個長方體計算域(12.0×1.4×1.4cm3)中，分辨率為192×48×48，如圖3 所示。入口邊界條件采用拋物線入流以減少流向發展距離，出口邊界設置為自由出流，管壁有限元結點作為浸沒邊界點通過耦合計算獲得不可滑移墻邊界。

圖3 圓管模型示意圖、沿程阻力系數與誤差隨Re 變化散點圖以及不同雷諾數下速度分布曲線

當管內流體充分發展后，數值結果與理論計算的沿程阻力系數吻合良好，相對誤差在2.0%左右，是可以接受的。此外，模擬得到的速度分布曲線與層流理論計算的速度分布曲線擬合的很好，這說明使用該代碼模型進行管內流模擬的數值結果是可靠的。

4 結果分析及討論

顆粒在流場中的運動被認為主要受拖曳力的影響，因此流速是顆粒運動的主導因素，本研究接下來將針對不同流速對管道的水力沖淤能力進行模擬和討論。為了節省計算資源，考慮使用較小管徑以模擬更多的顆粒。本文模擬的圓管內徑為D=1.0cm、長度為L=6.0cm，管道底部空間鋪設了均勻沙顆粒床，以表示間歇流所導致的沉積泥沙，顆粒質心間距等于顆粒直徑大小d=0.025cm，顆粒床覆蓋了管道一半的空間。

在管道滿流條件下，首先進行了三個低速層流工況的模擬，入口邊界分別采用平均流速為10cm/s、15cm/s、20cm/s 三種拋物線速度入流。由于主流速度驅動，顆粒上游床面先受到侵蝕，形成向下游擴散的沙坡；隨著泥沙顆粒逐漸向下游堆積，沙峰也逐漸向下游轉移；隨著顆粒數減少，沙床高度逐漸下降且趨于平坦；由于低速水流動力提供給顆粒運動的動力有限以及層流的速度場分布特點，當沙床低于某一高度，由水流提供的拖曳力不能能夠支持顆粒運動，最終大量的殘余顆粒形成靜止沙床留在管道內。從圖4 結果中可以看出，隨著層流流速遞增，沙床受侵蝕程度越大，沙峰越靠近下游，達到穩定的時間更早。

圖4 三種低速工況的顆粒床侵蝕情況、無量綱輸沙率和殘余顆粒數歷時曲線以及殘余顆粒百分比隨速度變化散點圖

在考慮管道泥沙淤積問題上，現行的排水管道設計標準一般以自清流速作為管道設計的最小流速。因此，為了能夠找到一個合適的最小設計流速范圍且能進一步跟實際工程問題結合起來，本研究繼續進行了六個湍流工況的模擬，模擬工況的流速區間從30～120cm/s。值得說明的是，這六個工況將采用湍流對數分布速度入流作為入口邊界條件，在邊壁處采用基于浸沒邊界法和大渦模擬相結合的湍流壁面模型[4]進行剪應力的修正，以得到一個正確的管內湍流。圖4 給出了湍流工況下最終殘余顆粒百分散點圖，并與三個低速層流工況進行了對比，可以看到從層流過度到湍流后，沖淤效果顯著提升，從20cm/s 到30cm/產生了約50%的跨度，這里可以認為是由于湍流與層流速度場特性不一致導致的，與層流場相比較，湍流場具有更大范圍的高速區，同時湍流邊界層幾乎貼近壁面。當流速大于30cm/s，清淤效果均達到90%以上；而當流速達到70cm/s，清淤效果更是超過了95%；流速超過100cm/s 后，清淤能力提升不再明顯。據查找相關手冊，滿流條件下，美國[5]和英國[6]對于城市排水管最小設計流速給出的標準分別為：90cm/s 和100cm/s，中國GBJ47-87 手冊和給排水工程快速設計手冊[7]中給出的是75cm/s。因此，可以看出當前模擬結果是能夠與實際的排水管設計標準相吻合的。

5 結論

本研究應用顆粒離散元和點球浸沒邊界法研究了水平管道內泥沙沖淤問題，并分析了不同流速下泥沙沉積床受到侵蝕的情況以及最終清淤能力。研究結果表明：a.在層流條件下，隨著流速增大，水流對顆粒沉積床的侵蝕能力雖然提升，但最終會殘留較多顆粒，清淤能力有限；b.當入流條件從層流過度到湍流后，沖淤效果顯著提升，這是由于湍流與層流速度場特性不一致導致的，湍流邊界層幾乎貼近壁面，與層流場相比較，其具有更大范圍的高速區；c. 當流速達到70cm/s，清淤效果更是超過了95%；流速超過100cm/s 后，清淤能力提升不再明顯。這與現行的國內外排水管最小設計流速標準相吻合的。