基于ANSYS CFX的泥石流數值模擬

李戰國

(新疆昌吉方匯水電設計有限公司，新疆 昌吉 831100)

泥石流是由山坡和溝渠產生的小石塊和大石頭構成的水和泥沙固液兩相流體。當滿足大量松散的固體物質、地形條件和一定強度的降水條件時，則可能會發生滑坡[1]。當山體滑坡發生時，泥石流裹挾的壓力會對山谷地形以及生命財產造成一定程度的損害。山體滑坡引起的巨大土石堆積會使肥沃土地變為荒涼戈壁，其對自然的破壞力不言而喻。目前，對泥石流的研究主要包括泥石流流體模型的研究、泥石流的預測預警、泥石流狀況評估、泥石流特性的數值模擬、各種泥石流預防結構和系統的研究等[2]。

1 研究區概況

洞峽子溝泥石流位于大川河右岸，南天門路通過溝槽(圖1)。

圖1 洞峽子溝溝谷范圍圖

從2011和2012年累計的滑坡總數來看，本次滑坡為中型滑坡。洞峽子溝地區地形通常為高腐蝕性的中型山地地形，其海拔高差在1 100～2 260 m之間，地勢陡峭，溝渠形狀通常為不規則梯形。洞峽子溝泥石流主要威脅溝陽凱村居民18戶68人和移民安置區45戶183人住宅共45棟(即將建造)，潛在威脅財產近1 300萬元，同時威脅溝渠南天門道路交通安全。

2 模擬方法

2.1 計算流體動力學數值模擬方法

目前在流體計算中使用3種數值方法，即有限差分法、有限元法和有限體積法(也稱為控制體積法)。其中，有限體積法具有計算效率高、收斂性好的特點，因此在CFD軟件的開發中得到廣泛的應用[3]。

有限體積法將分析對象分割成多個彼此不重復的控制體積，然后采用控制方程以積分插值的形式用于每個控制體積，體積界面的物理尺寸近似于物理節點的數量。區域離散化的本質是用有限的離散點代替原來的連續空間。離散區域過程首先將計算出的區域劃分為若干個不重疊的細分(即網格)，以確定每個細分節點的位置以及這些點所指示的控制點。這個區域管制流程完成后，4個最基本的要素是：①節點，需要解決的物理數量的幾何位置；②應用控制體積、控制方程或規律的最小幾何因素；③接口，接口規定了與每個節點對應的控制點的接口位置；④網格線，連接相鄰節點形成的曲線簇。在單個進程中，使用節點作為控制點的代表，在該節點上定義和存儲每個控制點的物理數量。圖2、圖3顯示了A、B和P是節點編號，φA、φB為節點邊界條件值，N、E、S、W是北、東、南、西，n、e、w、s是相應節點單元編號，Δx和Δy是一維、二維有限體積法的P節點的X和Y方向的長度。本文采用ANSYS CFX軟件進行相關模擬。

圖2 一維問題有限體積法計算網格

圖3 二維問題有限體積法計算網格

2.2 泥石流流變模型的選擇

研究泥石流流變的特性對于研究其流動規律、了解其性質和數值模擬滑坡具有重要意義。泥石流的流變特性揭示了從泥石流獲得的剪切應力與相應的液體流變速度梯度之間的關系。選擇流變模型有助于更深入地了解泥石流流變的內在規律。目前，研究人員已經根據實驗數據的理論推導和統計分析開發了多種流變模型。最常用的流變模型有牛頓流體模型、空流體模型、假塑性流體模型、膨脹流體模型[4]等。實驗表明，典型的泥石流流體特性的數學模型公式如下:

式中：τ為剪切應力；τb為屈服應力；η為黏性系數;du/dy為剪應變。

剪切應力主要是由泥石流中黏性粒子的凝集結構和沉積物粒子之間的摩擦形成的，黏性粒子的凝集結構是賓漢極限剪切應力的主要組成部分，其影響因素與粒度和微粒物質的含量有關，粒子越細，微粒含量越高，產生賓漢極限剪切應力的臨界濃度就越低。實驗證實，賓漢流體模型更適合一般泥石流流體的特性，因此本次模擬使用賓漢流體模型。

2.3 幾何建模

幾何建模選擇洞峽子溝主渠道大壩建模。結構形式為普通重力塊壩。

將選定的研究流域、大壩尺寸和材質與現有研究模擬計算器相結合，根據賓厄姆流體和湍流模型選擇單向流固耦合方法。

1) 壩體大，不易變形，變形對整個流場的流體分布影響很小。

2) 若使用雙向流固耦合需要設置移動網格，難以確定移動網格計算的準確性，過于依賴手動經驗，以及在研究體變形較小的前提下生成移動網格的計算誤差可能會更大。

3) 雙向流固耦合計算對計算機性能要求很高，計算速度較慢。此模型最好使用單向流固耦合，因為在采用單向-固態耦合的100 s、步長設置為0.1 s的英特爾i7-7800k 32gb ddr4平面下，累積時間超過44 h。

4) 根據經驗，單向流固耦合在計算詳細應力應變時比雙向流固耦合更加精確，前提是變形較小。因此，使用單向流固耦合模擬此模型。流固耦合模擬是兩個條件：第一個是大規模泥漿-流池溢流，第二個是中小土流逐漸沉入整個倉庫時的情況。流場和壩寬5 m，耦合分析的初始流場數據來自整個流域流場分析的相應部分。使用CFX執行流場計算，然后導入ANSYS執行流固耦合計算。流固耦合模型見圖4。如果將流場網格“軀干大小”設置為0.5m并加密FSI面，則“面大小”設置為0.25 m。工作條件2，進口面小，進口面也加密，面大小設置為0.2 m，最終結果節點(nodes)為54 544 607，單位(Element)為511 280(圖5)。

圖4 流固耦合模型

圖5 流體網格FSI細節加密圖

需要特別指出的是，FSI面大小必須與壩體網格的大小相匹配。否則，將出現overflow(溢出)錯誤。最終獲得壩網的節點為2 614 341，單位為623 280(圖6)。

圖6 流體網格細節展示圖

“基于SYM面的屬性”(Boundary type)是symmetry屬性。模擬壩體左右兩側的工作狀態與中間的工作狀態相匹配，因此選擇20 m的壩長進行模擬(圖7)。

圖7 SYM流場邊界圖

2.4 流固耦合計算參數設置

見表1。

表1 流固耦合參數設置一覽表

3 模擬結果

根據流場數據計算出的流體逐漸滲透到整個水箱中時，自由水平參考圖8，圖8顯示了每個時間速度向量，總液體壓力參考圖8。

圖8 逐漸淤積至滿庫過流工況泥流流體與攔擋壩相互作用時不同時刻的流速矢量圖

t=1 s時，二流流體開始流入二流塊水壩，傾斜將使二流流體“水龍頭”部分的速度逐漸加快。t=4 s時，泥漿流體“水龍頭”部分速度為12.3 m/s。t=8 s，當泥漿流體到達壩底時，泥漿流體“水龍頭”速度為9.9 m/s，此時泥漿流開始沉入壩內。t=38.4 s時，上面的泥漿流體沉入水平持續30.4 s到壩頂，當泥漿流體的液體水平足以通過泥漿水壩，第一次流向大流體時，流體速度為t=8 s的“水龍頭”速度明顯下降，從而逐漸降低液體水平。如果t=54.2 s，則土流水平會持續提高，在整個計算過程中達到最高值，下游流體可能會增加，而大壩前下游速度為4.5 m/s，則此過程可以更準確地說明大壩位置對控制土流的重要性。

將結果流場數據導入耦合計算后，獲得內部實體的應力應變和位移等主要參數。當泥漿流再次沉入整個水箱時，通道和水壩的液體最大壓力位置在水壩入口底部，壓力值為1.891 0×105Pa，水壩前面的最大壓力位置是負流沖擊區，沖擊壓力約為1.225 2×105Pa。

大壩的動態響應分析結果顯示，最大大壩應力強度位于大壩進水口表面的頂部，最大應力值為4.304 2×105Pa，大壩入口表面的應力強度也較大。塊壩的彈性變形強度分布基本上與壩的應力強度分布相匹配，最大應變值為2.224 1×105Pa。

在整個水庫過流條件下，塊壩是整體拱形層分布。隨著壩高的增加，位移值也會增加，最大位移值為4.790 0×10-5m。

4 結 論

泥石流災害的防治大部分采用經驗和規范，很難對所有泥石流流體的管理采用相同的加固標準。因此,要對每條泥石流溝進行單獨的數值模擬計算，以便對后防措施提供參考。在此次泥石流的數值模擬中，選擇合理的泥石流模擬參數，以便軟件模擬泥石流的起始、流動和整個流域流場情況。本文研究成果可為泥石流流場和大壩流固耦合計算、泥石流流場的動態特性和攔截工程的動態響應分析、泥石流流場動力學研究以及泥石流防治工程的理論研究提供參考和借鑒。