李戰(zhàn)國

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

泥石流是由山坡和溝渠產(chǎn)生的小石塊和大石頭構(gòu)成的水和泥沙固液兩相流體。當(dāng)滿足大量松散的固體物質(zhì)、地形條件和一定強(qiáng)度的降水條件時(shí)，則可能會(huì)發(fā)生滑坡[1]。當(dāng)山體滑坡發(fā)生時(shí)，泥石流裹挾的壓力會(huì)對山谷地形以及生命財(cái)產(chǎn)造成一定程度的損害。山體滑坡引起的巨大土石堆積會(huì)使肥沃土地變?yōu)榛臎龈瓯?，其對自然的破壞力不言而喻。目前，對泥石流的研究主要包括泥石流流體模型的研究、泥石流的預(yù)測預(yù)警、泥石流狀況評估、泥石流特性的數(shù)值模擬、各種泥石流預(yù)防結(jié)構(gòu)和系統(tǒng)的研究等[2]。

1 研究區(qū)概況

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

圖1 洞峽子溝溝谷范圍圖

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

2 模擬方法

2.1 計(jì)算流體動(dòng)力學(xué)數(shù)值模擬方法

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

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

圖2 一維問題有限體積法計(jì)算網(wǎng)格

圖3 二維問題有限體積法計(jì)算網(wǎng)格

2.2 泥石流流變模型的選擇

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

式中：τ為剪切應(yīng)力；τb為屈服應(yīng)力；η為黏性系數(shù);du/dy為剪應(yīng)變。

剪切應(yīng)力主要是由泥石流中黏性粒子的凝集結(jié)構(gòu)和沉積物粒子之間的摩擦形成的，黏性粒子的凝集結(jié)構(gòu)是賓漢極限剪切應(yīng)力的主要組成部分，其影響因素與粒度和微粒物質(zhì)的含量有關(guān)，粒子越細(xì)，微粒含量越高，產(chǎn)生賓漢極限剪切應(yīng)力的臨界濃度就越低。實(shí)驗(yàn)證實(shí)，賓漢流體模型更適合一般泥石流流體的特性，因此本次模擬使用賓漢流體模型。

2.3 幾何建模

幾何建模選擇洞峽子溝主渠道大壩建模。結(jié)構(gòu)形式為普通重力塊壩。

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

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

2) 若使用雙向流固耦合需要設(shè)置移動(dòng)網(wǎng)格，難以確定移動(dòng)網(wǎng)格計(jì)算的準(zhǔn)確性，過于依賴手動(dòng)經(jīng)驗(yàn)，以及在研究體變形較小的前提下生成移動(dòng)網(wǎng)格的計(jì)算誤差可能會(huì)更大。

3) 雙向流固耦合計(jì)算對計(jì)算機(jī)性能要求很高，計(jì)算速度較慢。此模型最好使用單向流固耦合，因?yàn)樵诓捎脝蜗?固態(tài)耦合的100 s、步長設(shè)置為0.1 s的英特爾i7-7800k 32gb ddr4平面下，累積時(shí)間超過44 h。

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

圖4 流固耦合模型

圖5 流體網(wǎng)格FSI細(xì)節(jié)加密圖

需要特別指出的是，F(xiàn)SI面大小必須與壩體網(wǎng)格的大小相匹配。否則，將出現(xiàn)overflow(溢出)錯(cuò)誤。最終獲得壩網(wǎng)的節(jié)點(diǎn)為2 614 341，單位為623 280(圖6)。

圖6 流體網(wǎng)格細(xì)節(jié)展示圖

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

圖7 SYM流場邊界圖

2.4 流固耦合計(jì)算參數(shù)設(shè)置

見表1。

表1 流固耦合參數(shù)設(shè)置一覽表

3 模擬結(jié)果

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

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

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

將結(jié)果流場數(shù)據(jù)導(dǎo)入耦合計(jì)算后，獲得內(nèi)部實(shí)體的應(yīng)力應(yīng)變和位移等主要參數(shù)。當(dāng)泥漿流再次沉入整個(gè)水箱時(shí)，通道和水壩的液體最大壓力位置在水壩入口底部，壓力值為1.891 0×105Pa，水壩前面的最大壓力位置是負(fù)流沖擊區(qū)，沖擊壓力約為1.225 2×105Pa。

大壩的動(dòng)態(tài)響應(yīng)分析結(jié)果顯示，最大大壩應(yīng)力強(qiáng)度位于大壩進(jìn)水口表面的頂部，最大應(yīng)力值為4.304 2×105Pa，大壩入口表面的應(yīng)力強(qiáng)度也較大。塊壩的彈性變形強(qiáng)度分布基本上與壩的應(yīng)力強(qiáng)度分布相匹配，最大應(yīng)變值為2.224 1×105Pa。

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

4 結(jié) 論

泥石流災(zāi)害的防治大部分采用經(jīng)驗(yàn)和規(guī)范，很難對所有泥石流流體的管理采用相同的加固標(biāo)準(zhǔn)。因此,要對每條泥石流溝進(jìn)行單獨(dú)的數(shù)值模擬計(jì)算，以便對后防措施提供參考。在此次泥石流的數(shù)值模擬中，選擇合理的泥石流模擬參數(shù)，以便軟件模擬泥石流的起始、流動(dòng)和整個(gè)流域流場情況。本文研究成果可為泥石流流場和大壩流固耦合計(jì)算、泥石流流場的動(dòng)態(tài)特性和攔截工程的動(dòng)態(tài)響應(yīng)分析、泥石流流場動(dòng)力學(xué)研究以及泥石流防治工程的理論研究提供參考和借鑒。