戈 凱

（西藏華泰龍礦業(yè)開發(fā)有限公司）

礦山充填將地表堆積的固體廢棄物回填至井下采空區(qū)，兼具平衡地壓和處置固體廢棄物的雙重作用，是實現(xiàn)無廢、綠色礦山建設(shè)的一項關(guān)鍵技術(shù)［1-3］，廣泛應(yīng)用于國內(nèi)外礦產(chǎn)資源的開采。充填體強度的穩(wěn)定性是影響采場安全及礦山企業(yè)收益的關(guān)鍵性因素，充填質(zhì)量的好壞本決定了采場能否穩(wěn)定生產(chǎn)，而充填配比是否合理又直接決定了礦山能否取得較好的經(jīng)濟效益。充填體的穩(wěn)定性主要由其力學(xué)性能反映出來，充填體的力學(xué)性能受多個過程影響，這些過程相互影響，共同控制著充填體的力學(xué)行為，包括充填體內(nèi)膠凝材料水化引起的熱過程、充填體內(nèi)外水流動引起的滲流過程、水化反應(yīng)料漿凝固形成充填體引起的力學(xué)過程，三者會互相影響。

目前，多場耦合研究巖土工程問題被學(xué)者們廣泛應(yīng)用，比如丁金華等［4］提出一種適用于膨脹土邊坡的濕度場—膨脹變形場—應(yīng)力場的多場耦合數(shù)值分析方法，闡述膨脹土邊坡淺層漸進性破壞機制的本質(zhì)是水力邊界條件變化引起的膨脹土膨脹變形作用。黃康橋等［5］采用改進的熱—水—力耦合模型與物理試驗進行驗證分析，將該耦合模型應(yīng)用于研究凍融作用下混凝土三維細觀尺度的性能演化規(guī)律。張培森等［6］基于Rock Top 多場耦合試驗儀，展開靜水壓力條件和三軸壓縮條件下的滲流試驗。礦山充填相關(guān)的多場耦合研究較晚，甘永生等［7］研究了流變場和溫度場耦合條件下充填料漿輸送數(shù)值模擬，結(jié)果表明在考慮溫度場條件下，模擬結(jié)果與其實際情況一致，可為料漿大倍線輸送設(shè)計提供理論依據(jù)。加拿大渥太華大學(xué)Fall 等［8-10］通過綜合化學(xué)、熱力學(xué)及力學(xué)過程得出了一個耦合模型，并利用試驗數(shù)據(jù)對模型進行了驗證。關(guān)于采場充填體的研究多集中在其力學(xué)行為分析，在理論基礎(chǔ)上研究充填體成拱應(yīng)力發(fā)展過程以及其應(yīng)力分布和發(fā)展過程對采場圍巖、充填擋墻的力學(xué)影響。為研究高海拔高寒環(huán)境下充填體穩(wěn)定性，本研究通過建立合理有效的多場耦合模型，采用熱—流—力耦合方式定量研究充填體穩(wěn)定性。

1 充填體的力學(xué)性能多場耦合過程

充填體的力學(xué)性能受熱過程、滲流過程和力學(xué)過程的影響。

（1）熱過程。由于充填體內(nèi)部膠結(jié)劑遇水發(fā)生化學(xué)反應(yīng)，放出相應(yīng)的熱量，這會使得充填體溫度升高，加快水化反應(yīng)，進而提高充填體熱量的產(chǎn)生，并且溫度升高導(dǎo)致的熱膨脹效應(yīng)也是充填體力學(xué)行為變化的組成部分。

（2）滲流過程。CPB（膠結(jié)充填料漿）中的水一部分會與黏結(jié)劑發(fā)生水化反應(yīng)進行消耗，另一部分會由于水的自重應(yīng)力作用發(fā)生滲流現(xiàn)象，這也會對充填體強度產(chǎn)生影響，且水的流動會帶動熱量的流動，進而影響充填體中膠結(jié)劑水化反應(yīng)的不均勻分布。

（3）力學(xué)過程。水化反應(yīng)產(chǎn)生的水化產(chǎn)物積聚，充填體固相增多，骨架形成，強度產(chǎn)生，且礦壓、水壓等也會對充填體的穩(wěn)定性產(chǎn)生影響，充填體強度形成后，會導(dǎo)致充填體空隙率的改變，進而影響料漿的滲流過程。

在實際情況中，3 個物理過程會互相影響，因此需要對其進行耦合，三者的耦合過程如圖1所示。

2 充填體熱—流—力耦合模擬

2.1 充填體熱—流—力耦合數(shù)學(xué)模型

本研究采用COMSOL Multiphysics 數(shù)值模擬軟件，其具有耦合多物理場的能力。COMSOL Multiphysics 允許用戶在內(nèi)置的物理界面中編寫偏微分方程，以實現(xiàn)耦合多物理場的目標(biāo)。在COMSOL Multiphysics 中，解決特定問題的整個過程是構(gòu)建幾何模型、設(shè)置邊界條件、輸入材料參數(shù)、網(wǎng)格化并計算可視化結(jié)果。

充填料漿填充到采場中后，由于黏合劑的水合反應(yīng)，水合產(chǎn)物將逐漸積聚，充填體將逐漸顯示出多孔介質(zhì)的特性。COMSOL Multiphysics 提供了一個描述多孔介質(zhì)傳熱的通用數(shù)學(xué)模型：

keq=φ( )ksolid-kfluid+kfluid， （3）

式中，ρ為充填體的密度，kg/m3；Cp為充填體的比熱容，J/（kg·K）；keq為熱導(dǎo)率，W/（m·K）；( )ρCp eq為恒定壓力下的等效體積熱容量，J/（m3·K）；u為速度場，m/s；QH為熱源，W/m3；φ為充填體的孔隙率；ρsolid為CPB 中固體基質(zhì)的密度，kg/m3；Csolid為CPB 中固體基質(zhì)的熱容量，J·K；ρfluid為CPB中流體的密度，kg/m3；Cfluid為CPB中流體的熱容量，J·K；ksolid為CPB中固體基質(zhì)的熱導(dǎo)率，W/（m·K）；kfluid為CPB 中流體的熱導(dǎo)率，W/（m·K）。

充填體的滲流過程可由式（4）表示。

式中，ki為CPB 的固有滲透率，m/d；kr為CPB 對水的相對滲透率，m/d；μw為CPB 中水流的速度場，m/s；pw為孔隙水壓力，Pa；Qw為排水量，m3。

充填體中的總應(yīng)力由有效應(yīng)力以及水和空氣的壓力組成，總應(yīng)力矢量σt為

σt=σe+βδijpa， （5）

pa=Sfluidpfluid+Sairpair， （6）

式中，σe為有效應(yīng)力，Pa；β為比奧系數(shù)；δij為克羅尼克參數(shù)；pa為水和空氣的平均壓力，Pa；Sfluid為CPB中流體的飽和度；pfluid為CPB 中流體的壓力，Pa；Sair為CPB 中空氣的飽和度；pair為CPB 中空氣的孔隙水壓力，Pa。

有效應(yīng)力σe是熱膨脹、孔隙水壓力和充填體重力引起的總應(yīng)力，如式（7）所示。

式中，σt為總矢量壓力，Pa；ξij為克朗尼參數(shù)，ξii= 1。

綜上，充填體在井下會受到多種物理場的影響，由于受到料漿內(nèi)部如膠結(jié)劑的水化反應(yīng)和外部如高寒高海拔的影響，充填體會表現(xiàn)出復(fù)雜的多物理場耦合的現(xiàn)象。本研究建立的多場耦合模型主要包括以下幾個方面。

（1）熱學(xué)過程。料漿充入采場之后，膠結(jié)劑與水發(fā)生水化反應(yīng)產(chǎn)生熱量，這部分熱量會在充填體之間以及外界發(fā)生熱傳遞現(xiàn)象，主要運用傅里葉定律和COMSOL 內(nèi)置的導(dǎo)熱公式闡明熱傳遞過程。在熱傳遞過程中，充填體中的一些力學(xué)參數(shù)會隨著溫度的變化而發(fā)生改變，闡述了充填體力學(xué)參數(shù)（如流體密度）與溫度的關(guān)系變化。

（2）滲流過程。料漿中的水由于充填體水化反應(yīng)后的密實過程擠出充填體，水在流動過程中，受到膠結(jié)劑水化反應(yīng)熱的影響，水的黏度等參數(shù)會發(fā)生改變。

（3）力學(xué)過程。充填體中的總應(yīng)力由有效應(yīng)力以及水和空氣的壓力組成，基于Drucker—Prager 屈服準(zhǔn)則的演化屈服函數(shù)定義了初始屈服和荷載函數(shù)。充填體的力學(xué)過程受到熱學(xué)過程和流體過程的雙重影響，膠結(jié)劑水化反應(yīng)放出的熱量造成了充填體的熱膨脹效應(yīng)，而流體的流動帶動了熱量的分布，進而影響充填體的力學(xué)行為，最終充填體力學(xué)行為的改變又會反過來影響熱學(xué)過程和流體過程。

2.2 充填體數(shù)值模型建立

建立幾何模型的過程及模擬過程包括添加物理場接口、建立幾何模型、添加初始條件、設(shè)置邊界條件、材料添加、網(wǎng)格剖分和后處理。本研究所選用的物理場接口為“固體力學(xué)模塊”“達西定律”和“傳熱模塊”，由于模擬結(jié)果是隨著時間變化的，因此選擇“瞬態(tài)”模擬。在數(shù)值模擬軟件中，選擇2 個工作平面，分別為XY平面和ZX平面，將采場平面圖和剖面圖分布導(dǎo)入2個工作平面中，并根據(jù)采空區(qū)實際長度進行拉長，拉伸的長度應(yīng)大于采空區(qū)實際長度。之后通過改變坐標(biāo)使2個拉伸實體進行重合，最終通過布爾操作中的交集得到符合實際的充填體模型。建立的幾何模型網(wǎng)格剖分結(jié)果如圖2所示。

為了方便后續(xù)更直觀地觀測分析充填體應(yīng)力及位移的分布，在充填體幾何模型中布置了縱斷面及橫斷面2 個斷面，縱斷面距離Y軸零點坐標(biāo)6 m，橫斷面距離充填體底面邊界10 m，斷面布置如圖3所示。

2.3 模擬結(jié)果分析

2.3.1 充填體縱斷面應(yīng)力變化

為了研究充填體在X軸方向上的水平應(yīng)力和位移變化，在后處理過程中設(shè)置了3條三維截線，如圖4所示。3條三維截線分別距底板3，10和20 m，所處的Y軸坐標(biāo)與充填體縱斷面一致。

圖5 為3，10 和20 m 截線各個位置的應(yīng)力分布，每一條不同標(biāo)記的線代表模擬的養(yǎng)護時間。圖5 表明，隨著模擬時間的增加，充填體應(yīng)力也在增加，充填體的應(yīng)力分布自上而下是逐漸增加的。單獨分析圖5（a），由于圖4 a 線0 點坐標(biāo)處受到邊界條件的支撐作用，其應(yīng)力分布較低，隨著向線的中心位置推進，應(yīng)力不斷增加，這是因為中心受到兩側(cè)充填體的擠壓作用，使其應(yīng)力增加明顯。而線的終端位置出現(xiàn)了明顯的應(yīng)力集中現(xiàn)象，這可能是因為所建立的幾何模型為不規(guī)則體，邊界的不規(guī)則性導(dǎo)致了應(yīng)力集中的現(xiàn)象，圖4 b線的基本情況與其類似。圖4 c線由于邊界較為規(guī)則，線的兩端并沒有出現(xiàn)應(yīng)力集中的現(xiàn)象，出現(xiàn)應(yīng)力集中現(xiàn)象的位置為線的中心位置，這可能是由于充填體的頂部為拱結(jié)構(gòu)，拱的兩側(cè)位置的料漿對圖4c 線中心處產(chǎn)生擠壓作用，導(dǎo)致了其中心位置的應(yīng)力集中現(xiàn)象。

2.3.2 充填體縱斷面位移變化

圖6 為3，10 和20 m 截線各個位置的位移分布。圖6表明隨著模擬時間的增加，充填體發(fā)生的位移逐漸加大，充填體的位移自上而下逐漸降低，且由于充填體下方受到底板的支撐作用，很難發(fā)生位移現(xiàn)象，因此位移為零，充填體上方與外界接觸，為自由面，位移為最大。分析圖6（b）和圖6（c）可知，三維截線的中心位置的位移較兩端大，這是由于在幾何模型邊界條件的設(shè)定中，四周邊界為輥支撐，使兩端發(fā)生的位移并不明顯，而線的中心位置由于受到兩側(cè)的擠壓和上覆料漿的重力作用，產(chǎn)生的位移更大。而圖4 a線產(chǎn)生此種形態(tài)可能是由于幾何模型的不規(guī)則性導(dǎo)致的，0 點坐標(biāo)處幾何模型邊界的不規(guī)則性使其發(fā)生了更加明顯的位移現(xiàn)象。

2.3.3 充填體橫斷面應(yīng)力變化

為了研究充填體在Z軸方向的應(yīng)力及位移變化，在幾何模型中布置3 條豎直方向的三維截線，如圖7所示，分別距X軸0點10，25和40 m。

圖8 為10，25 和40 m 三維截線所對應(yīng)的應(yīng)力分布，不同標(biāo)記的線代表模擬時間。對比3條線相同高度相同模擬時間所對應(yīng)的應(yīng)力大小，如高度為10 m、模擬時間為8 d，圖7 a 線的應(yīng)力約為0.12 MPa，圖7 b線的應(yīng)力約為0.14 MPa，圖7 c 線的應(yīng)力約為0.115 MPa，表明同一高度上，應(yīng)力從中心向兩端逐漸降低。分析單條三維截線上的應(yīng)力分布，如圖8（a）可以看出，在同一高度上，養(yǎng)護時間越長，其應(yīng)力越大，這是由于養(yǎng)護時間越久，充填體中膠結(jié)劑的水化反應(yīng)越充分。養(yǎng)護時間固定時，距離底板越遠，應(yīng)力越低，這是因為充填體下部封閉，熱量的集中加速了膠結(jié)劑的水化反應(yīng)以及料漿自重應(yīng)力的作用。

2.3.4 充填體橫斷面位移變化

圖9 為10，25 和40 m 三維截線不同高度、不同模擬時間的位移變化。對比3條三維截線相同高度、相同模擬時間的位移情況，如高度為10 m、模擬時間為12 d，圖7 a線的位移約為6.1 cm，圖7 b線的位移約為7 cm，圖7 c 線的位移約為6 cm，表明相同高度、相同養(yǎng)護時間的條件下，充填體的位移自中心向兩端遞減，這是由于邊界條件的限制作用影響以及采場幾何拱兩端擠壓。分析單條三維截線，如圖9（a）可以看出，同一高度條件下，充填體養(yǎng)護時間越久，其位移越大，這是由于料漿中膠結(jié)劑水化反應(yīng)不斷消耗水分以及料漿自重應(yīng)力作用。當(dāng)養(yǎng)護時間一致時，距底板越遠，其位移越大，這是由于料漿中的水自下而上排出，壓縮了料漿自重應(yīng)力作用。

3 結(jié) 論

（1）運用驗證過的數(shù)學(xué)模型，對采場充填后的力學(xué)行為變化進行熱—流—力多場耦合模擬預(yù)測，其結(jié)果與實際相符合，驗證了數(shù)值模擬技術(shù)手段研究充填體力學(xué)行為的適用性。

（2）采場中的應(yīng)力自充填體上表面至底板逐漸增加。一方面是膠結(jié)劑與水發(fā)生水化反應(yīng)，水化反應(yīng)產(chǎn)生的熱量在采場下部封閉區(qū)域積聚，反過來加速了膠結(jié)劑的水化反應(yīng)，使充填體下部充填體強度更大；另一方面是水自下而上滲出，使充填體上部積聚了更多水分，降低了上部料漿的濃度，使上部充填體強度難以進一步提高，以及料漿的自重應(yīng)力作用加大了下部充填體的應(yīng)力分布。

（3）采場中的位移自充填體上表面至底板逐漸減小。由于在模擬軟件中，采場模型的上表面為自由面，其余面由于圍巖的支撐作用，表現(xiàn)在模型中為輥支撐和固定約束，使得充填體邊界難以發(fā)生位移變化，且由于自重應(yīng)力作用，充填體上部不斷發(fā)生沉降現(xiàn)象，而下部沉降現(xiàn)象并不明顯，因此位移自上而下逐漸減小。