戈 凱

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

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

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

1 充填體的力學(xué)性能多場(chǎng)耦合過(guò)程

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

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

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

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

在實(shí)際情況中，3 個(gè)物理過(guò)程會(huì)互相影響，因此需要對(duì)其進(jìn)行耦合，三者的耦合過(guò)程如圖1所示。

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

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

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

充填料漿填充到采場(chǎng)中后，由于黏合劑的水合反應(yīng)，水合產(chǎn)物將逐漸積聚，充填體將逐漸顯示出多孔介質(zhì)的特性。COMSOL Multiphysics 提供了一個(gè)描述多孔介質(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為速度場(chǎng)，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）。

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

式中，ki為CPB 的固有滲透率，m/d；kr為CPB 對(duì)水的相對(duì)滲透率，m/d；μw為CPB 中水流的速度場(chǎng)，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。

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

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

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

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

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

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

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

2.3 模擬結(jié)果分析

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

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

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

2.3.2 充填體縱斷面位移變化

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

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

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

圖8 為10，25 和40 m 三維截線所對(duì)應(yīng)的應(yīng)力分布，不同標(biāo)記的線代表模擬時(shí)間。對(duì)比3條線相同高度相同模擬時(shí)間所對(duì)應(yīng)的應(yīng)力大小，如高度為10 m、模擬時(shí)間為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)護(hù)時(shí)間越長(zhǎng)，其應(yīng)力越大，這是由于養(yǎng)護(hù)時(shí)間越久，充填體中膠結(jié)劑的水化反應(yīng)越充分。養(yǎng)護(hù)時(shí)間固定時(shí)，距離底板越遠(yuǎn)，應(yīng)力越低，這是因?yàn)槌涮铙w下部封閉，熱量的集中加速了膠結(jié)劑的水化反應(yīng)以及料漿自重應(yīng)力的作用。

2.3.4 充填體橫斷面位移變化

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

3 結(jié) 論

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

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

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