基于CFD-DEM耦合的尾礦庫潰壩數值模擬

楊 柳 王金安 黎偉佳

(1.北京科技大學土木與資源工程學院,北京 100083;2.金屬礦山高效開采與安全教育部重點實驗室,北京 100083)

尾礦庫是用于堆放礦山尾礦的場所,是維持礦山開采所必須的工業設施。尾礦壩則指用于儲存尾礦的外圍壩體構筑物[1]。由于尾礦庫是一個具有高勢能的人造泥石流危險源,一旦發生危險,將對下游人民的生活和生產造成極大的災害。因此,尾礦庫的安全運營和尾礦壩的穩定性一直是礦山開采中的重要研究問題[2-4]。

當前,對于尾礦庫的穩定性研究通常采用有限元軟件[5-7]或者極限平衡法[8]、強度折減法[9]等基于連續介質假設的研究方法。然而,實際的尾礦庫中儲存的并不是連續、未擾動的土體,而是礦石選別后的松散堆積體,由于水力填充尾礦料沉積過程的隨機性,以及尾礦料自身復雜的固液耦合作用,使得尾礦壩的物理力學性質比一般土石壩更為復雜。因此,基于計算流體力學(CFD)和離散元(DEM)的流固耦合算法更適用于模擬尾礦庫潰壩過程,能更深入地解釋潰壩機理,更真實地反映潰壩過程中尾礦料的力學行為。

關于耦合算法在尾礦庫穩定性方面已經有諸多研究,李強等[10]將連續介質流固耦合法與強度折減法相結合,對尾礦壩滑移面進行了分析,獲得了尾礦壩內三維滑移面的空間分布。張鐸等[11]在連續介質力學有限差分的基礎上,選取尾礦壩在填充前后潛在滑移面上的典型局部區域進行離散—連續耦合分析,獲得了滑移帶形成過程的細觀力學特征。王君等[12]基于有效應力理論和動孔壓的應力應變模型,結合FLAC3D進行了地震作用下尾礦庫地基液化流固耦合分析,得到不同位置的孔壓、加速度、沉降位移等信息,提供了一種快速判別尾礦庫液化區域的方法。

上述關于尾礦庫的耦合算法研究主要集中于連續介質理論的流固耦合算法,而關于離散元理論的流固耦合算法鮮有報道。因此,本研究基于離散元(DEM)理論的SDEM軟件對尾礦庫潰壩過程進行模擬,并在此基礎上用 Fortran編寫計算流體力學(CFD)和CFD-DEM流固耦合程序,對某尾礦庫潰壩過程進行模擬,獲得潰壩的滑移面、位移場和速度場,根據速度場推斷受災區域,并將模擬結果與實際結果進行對比與驗證。

1 流固耦合計算基本理論

1.1 流體計算方程

流體介質需要同時滿足質量、動量以及能量的三大守恒方程,可采用經典的Navier-Stokes方程進行統一描述[13]。考慮孔隙率的三維不可壓縮粘性流動Navier-Stokes方程為

式中,n為孔隙率;V為流動速度矢量;φ為特征變量,在不同的方程中代表不一樣的變量;Sφ為源項;u、v、w為x、y、z3個方向上的速度。

將上述控制方程中變量φ、Γφ、源項Sφ按照表1取值,則分別獲得流體的連續方程、動量方程和能量方程。表1中Fx、Fy、Fz分別表示x、y、z方向上的單位質量流體所受的外部體力;Smx、Smy、Smz分別表示x、y、z方向上的單位質量流體的質量源項;μ為液體的動力粘性系數;k為熱導率;e為單位質量流體內能;p為流體壓強;?為耗散產生的熱能;Q是單位時間外部施加給單位體積流體的熱能。

表1 流體計算方程中各變量取值Table 1 Values of variables in fluid calculation equations

1.2 顆粒控制方程

1.3 流固耦合計算

2 流固耦合程序實現

采用SDEM軟件中的DEM模塊,并在此基礎上用Fortran語言編寫CFD模塊以及CFD-DEM交互模塊。其中CFD模塊通過N—S方程離散,再采用Issa[18]提出的壓力隱式算子分裂法(PISO算法)進行求解。該程序主要分為3個部分:流體的壓力場和速度場計算、顆粒運動過程計算、流體與顆粒相互作用計算。計算過程中主要涉及2個判斷:

(1)判斷流體計算程序與顆粒計算程序的物理時間是否相等。若兩者的物理時間相等,則通過流體與顆粒相互作用計算程序進行計算。其余時間內,流體與離散單元間的相互作用力保持不變。

(2)計算離散單元所在流體區域網格,根據流體區域的壓力大小來判斷顆粒是否處于流體當中。若處于流體中,則需計算顆粒與流體的相互作用力和方向。CFD-DEM總計算程序流程如圖1所示。

圖1 CFD-DEM計算程序流程Fig.1 Calculation program flow of CFD-DEM coupled

3 數值模擬參數確定

離散介質參數標定實質上是建立起細觀參數與宏觀參數的聯系,使數值分析結果更具有可靠性。離散單元的參數標定方法往往先確定會影響模擬結果的宏觀參數,而同一宏觀行為通常受到一個或多個的細觀參數的影響,一個細觀參數的選擇可對多種宏觀行為造成影響。因此,在建模之前,先建立一個尺寸為0.05 m×0.05 m×0.1 m的試樣,對不同細觀參數的試樣進行單軸壓縮,并與現場試樣進行對比,從而確定模擬參數。由于標定參數較多,此處僅以彈性模量、摩擦系數和粘結強度為例,介紹標定過程及結果。單軸壓縮模型如圖2所示。

圖2 單軸壓縮模型Fig.2 Uniaxial compression model

3.1 彈性模量標定

試樣的切向剛度和法向剛度決定了宏觀顆粒集合體的彈性模量,影響應力應變曲線的斜率,對破壞時的應變產生直接的影響。由圖3可知,細觀接觸剛度與顆粒介質集合體的宏觀彈性模量呈線性正相關關系,對抗壓強度的影響不大。細觀接觸剛度(x)與顆粒介質集合體的宏觀彈性模量(y)在數值上可以根據圖3(a)擬合表示為

圖3 接觸剛度對宏觀彈性模量、抗壓強度的影響Fig.3 Influence of contact stiffness on macroscopic elastic modulus and compressive strength

經室內土工試驗測得的宏觀彈性模量為18.46 MPa,則細觀接觸剛度取19.9 MPa。

3.2 摩擦系數和粘結強度標定

摩擦系數和粘結強度對試件的破壞形態會造成直接的影響,從而也影響力鏈的傳播方向,對壓縮強度造成一定的影響。不同摩擦系數條件下粘結強度和抗壓強度之間的關系如圖4所示。

圖4 不同摩擦系數條件下粘結強度和抗壓強度之間的關系Fig.4 Relationship between bond strength and compressive strength under different friction coefficients

對圖4進行擬合,得

式中,x為細觀粘結強度;μp為摩擦系數;y為宏觀壓縮強度。由于庫內水體的作用,摩擦系數一般取0.1,即當粘結強度為0.17 MPa時,壓縮強度為0.4 MPa。

綜上,在尾礦庫潰壩模擬試驗中,細觀參數的選取如表2所示。

表2 尾礦庫潰壩模擬試驗中的主要計算參數Table 2 Main calculation parameters in the simulation test of tailings pond dam break

4 尾礦庫潰壩過程的數值模擬

4.1 工程概況

以山西某尾礦庫潰壩作為工程背景,該尾礦庫剖面如圖5所示。尾礦庫總壩高為49.0 m,總庫容約39.6萬m3,儲存尾砂約31.7萬m3,屬于四級尾礦壩。下游子壩坡比為1∶3,壩頂標高為26.4 m;上游子壩坡比為1∶1.38,壩頂標高為49.0 m。初期壩高度為8.0m,為漿砌石壩。在潰壩事故發生前10 d,尾礦庫壩體表面出現滲水現象,這種情況極易出現滑坡型潰壩。現場采取的防治措施是在壩體表面覆蓋一層約4m厚的黃土以堵水,造成庫內水位升高,浸潤線貼近堵水隔層,黃土浸水而軟化。

圖5 尾礦庫剖面示意(單位:m)Fig.5 Profile schematic of tailings pond

4.2 尾礦庫模型建立

尾礦庫離散元模型及計算區域如圖6所示,1區為尾礦庫,2區表示水域,并在壩體上設置5個監測點,檢測潰壩過程中該點的運動軌跡。

圖6 計算模型Fig.6 Calculation model

尾礦庫的初期壩為漿砌石壩,在現場調研中發現,其并未遭受太大的破壞及位移,因此在數值模擬當中,固定漿砌石壩處的顆粒。模型的水壓力如圖7所示。由于坡面覆蓋了1層約4m厚的黃土,故在坡面上的流體網格的速度設置為0.0 m/s。由于黃土的隔水作用,壩體內水頭壓力基本保持水平,最大的水壓力位于壩體底部,壓力值為48.482 kPa。

圖7 水壓力場示意Fig.7 Schematic diagram of hydraulic force field

4.3 潰壩過程模擬

根據上述得到的細觀參數,對尾礦庫潰壩事故進行模擬,研究其潰壩形成與發展過程。

由圖8可知,當迭代2 000步(實際物理時間為3.6s)時,壩體形成貫通的滑動面,在滑移面后方存在1個受擠壓的區域。滑移面的劃分是按照在y方向位移量的正負值進行劃分的。滑移面進口位于接近壩頂的位置,滑面深且陡;出口位于坡比為1∶3的子壩中部位置,滑坡面平緩。位移最大的部位位于坡比為1∶1.38和坡比為1∶3的子壩的交接處,位移最大值約為2.6 m,符合該尾礦庫潰壩的事故調查報告中壩體中部位置首先隆起破壞的事實。

圖8 2 000迭代步時位移矢量圖Fig.8 Displacement vector graph at 2 000 steps

在運行至3 000迭代步時(實際物理時間為5.4 s),滑動面處的顆粒迅速崩落。在滑動面入口處產生貫通性的拉裂縫,拉裂縫后方物質由于失去支撐而處于崩落的臨界狀態,壩頂后方局部位置的拉裂縫也逐漸發展,導致壩頂后方尾砂也迅速崩落。如圖9所示,左上角為側視圖,拉裂縫貫穿了整個模型,而在拉裂縫后方的顆粒的位移值最大。

圖9 3 000迭代步時位移示意Fig.9 Displacement at 3 000 steps

在運行至4 000迭代步時(實際物理時間為7.2 s),原來的分段坡面已經演變成一個坡面,如圖10所示。從壩體的位移分布來看,壩體的位移呈現一種連續弧形的分布。這種位移的分布主要是由于z方向(重力方向)位移造成的,而不是y方向的位移。在原壩頂處的顆粒z方向的位移最大,在漿砌子壩處z方向的位移為正值,向下游運移,隨后尾礦砂越過漿砌子壩。y方向位移分布出現明顯的弧形分界面,深達模型的底部位置。

圖10 4 000迭代步時位移示意Fig.10 Displacement at 4 000 steps

潰壩過程中下泄的顆粒物質形態如圖11所示。當下泄物質越過初期壩向下游移動時,顆粒物質較為密實,下游顆粒速度迅速到達最大值17.267 m/s后緩慢減速。在潰壩發展過程中,滑動面上的尾砂具備更高的容積及動能,隨后的尾砂是由于失去下方物質支撐而逐漸補充至下游下泄物質當中。隨后由于相互碰撞,后方下泄物質給下游處物質提供持續的動能,而滯留于距離漿砌子壩不遠處,從而使兩端的物質更多,呈現出最終堆積狀態呈現兩頭密中間疏的現象。

圖11 潰壩過程中顆粒速度Fig.11 Particle velocity in the process of dam break

壩體處5個監測點的軌跡示意如圖12所示。由圖12可知,壩體底部顆粒的運移距離大于壩體上部顆粒。位于距離壩頂100 m處的1號監測點位移主要以垂直方向為主。而處于滑動面上的2號、3號、4號監測點,在潰壩過程中,其主要的位移方向依次為垂直—水平相互交替運移。而這種現象隨著高度增加而越加明顯。位于初期壩上的5號監測點位移主要為水平方向,并在潰壩初期在走向方向(即X方向)有一定波動。

圖12 監測點運移軌跡Fig.12 Trajectory of monitoring points

顆粒在Y方向的速度場如圖13所示。在迭代步數為10 000步時(實際物理時間為18 s),下泄物質達到Y方向最大速度17.267 m/s,位于距離子壩處約150 m處。隨著潰壩的演進,最大速度梯度所在位置的Y方向坐標逐漸增加。假設只有單個顆粒崩落的情況下,沒有顆粒間的相互碰撞,在重力作用下,理論上應當在壩底處達到最大速度,隨后在與邊界的摩擦作用下,速度逐漸降低。但在圖13中可以看出,速度最大值往往是處于下泄物質的前段,而不是在尾礦壩壩體底部。這是由于上游物質不斷補充,與下游物質相互碰撞,動能發生轉移造成的。

圖13 潰壩演進過程Y方向速度場Fig.13 Y-direction velocity field in the evolution process of dam break

在迭代步數大于60 000(實際物理時間為108 s)步時,由于上游補充的物質逐漸減少,龍頭處的下泄物質的速度衰減率比之前大,速度梯度減少。距離尾礦壩壩底越近的下泄物質越早停止運動。大部分物質停留在Y方向坐標為0～2 000m的位置,少部分物質由于沒有其他顆粒的阻擋作用,最終的運移距離可以到達2 000 m以上。由此可以看出,在一定范圍內,當庫容量越大時,下泄物質的最遠運移距離越大。在迭代步數為120 000步時(實際物理時間為216 s時),下泄物質基本停止。在不足4 min的時間內,最遠運移距離已經達到2 000余m,導致在實際過程中,附近的居民難以逃離災害區。

隨后峰值速度不斷減少,峰值速度隨時間變化曲線如圖14所示。變化曲線可由2條直線組成,加速階段時,平均加速度為0.959 3 m2/s;減速階段平均加速度為-0.063 95 m2/s,最遠運移距離為S=17.2672/(2×0.959 3)+17.2672/(2×0.063 95)=2 486.51 m,事故調查報告中潰壩影響范圍為2.5 km。實地考察中距離壩底約350 m處的磚混結構房屋已被整體沖毀,該處速度峰值為16.32 m/s。

圖14 潰壩過程中最大速度Fig.14 Maximum velocity in the process of dam break

4 結 語

(1)在離散元軟件SDEM的基礎上,用Fortran語言編寫了CFD-DEM流固耦合程序,流體與顆粒間的相互作用力僅考慮了拖曳力、浮力和壓力差,并闡述了流體計算模塊和離散元計算模塊的耦合過程。

(2)根據室內試驗的土樣數據和試樣的數值模擬結果,確定尾礦砂的摩擦系數,彈性模量和粘結強度等細觀參數,建立尾礦庫及相應水域的數值模型。數值模擬結果與潰壩事故報告結果相吻合,證明了該流固耦合算法的可靠性。

(3)潰壩形成初期,壩體中部位移最大,隨后形成貫通的滑動面,滑移面后方有一個受擠壓的區域;之后滑動面迅速崩塌,在滑動面入口處產生貫通性拉裂縫;尾礦壩潰壩后形成的堆積體呈兩頭密中間疏,在漿砌子壩前方形成明顯的分割帶。