采場(chǎng)非膠尾砂充填料漿自重固結(jié)過程數(shù)值模擬

馬青海，劉光生，楊小聰，郭利杰

(1.北京礦冶研究總院，北京 100160； 2.金屬礦綠色開采國(guó)家國(guó)際聯(lián)合研究中心，北京 102628； 3.礦冶科技集團(tuán)有限公司，北京 102628)

0 引 言

充填采礦法既能改善地表尾砂堆存問題，又能處理井下采空區(qū)，已成為地下金屬非金屬礦山主流采礦方法[1-3]。選廠尾砂是采場(chǎng)空區(qū)充填的主要骨料，尾砂充填工藝以水為載體，通過管道將尾砂料漿輸送至井下空區(qū)。尾砂料漿充入采場(chǎng)后，須盡快排出自由水并達(dá)到一定固結(jié)程度才能發(fā)揮其地壓控制作用，否則在采空區(qū)長(zhǎng)期堆積的難脫水固結(jié)的流態(tài)狀尾砂料漿反而會(huì)構(gòu)成礦山井下重大危險(xiǎn)源。

以往礦山充填前，大多先將選廠全尾砂中的細(xì)粒分出(多排放至尾礦庫(kù))，使用粒級(jí)組成較粗的分級(jí)尾砂作為骨料，其滲透性較好，分級(jí)尾砂料漿的自由水可快速排出采場(chǎng)完成固結(jié)。但近年來隨著選礦工藝發(fā)展，選廠磨礦產(chǎn)出的全尾砂顆粒越來越細(xì)，而且全尾砂細(xì)粒占比也越來越高，加上礦區(qū)環(huán)保要求趨緊，細(xì)粒尾砂地表排放難，促使礦山盡量選用不分離細(xì)顆粒的全尾砂作為充填骨料。

全尾砂細(xì)粒多，保水性好，滲透性差[4]，導(dǎo)致采場(chǎng)全尾砂料漿較難快速脫水完成固結(jié)。尤其是選用階段空?qǐng)鏊煤蟪涮畈傻V法的大型礦山，其采場(chǎng)尺寸和綜合生產(chǎn)能力大，導(dǎo)致采空區(qū)充填料漿上升速度快、料漿滲流排水路徑長(zhǎng)等問題突出，而且二步驟采場(chǎng)充填為節(jié)約成本往往不添加膠凝材料[5-6]，迫使礦山必須充分掌握采場(chǎng)全尾砂料漿的脫水固結(jié)過程。近期，國(guó)內(nèi)外報(bào)道的大型采場(chǎng)尾砂料漿沖垮擋墻造成的安全生產(chǎn)事故時(shí)有發(fā)生[7]，因此，科學(xué)合理地評(píng)估采場(chǎng)非膠尾砂料漿的固結(jié)過程是礦山充填體力學(xué)領(lǐng)域的重要研究熱點(diǎn)之一。

國(guó)內(nèi)學(xué)者目前對(duì)采場(chǎng)尾砂充填料漿脫水固結(jié)方面的研究主要集中在采場(chǎng)排水措施應(yīng)用上，包括脫水管脫水[8]、負(fù)壓脫水[9]和電滲脫水[10]等，較少關(guān)注采場(chǎng)尾砂漿固結(jié)過程的力學(xué)機(jī)理。國(guó)外學(xué)者針對(duì)膏體充填料漿固結(jié)過程開展了較多原位監(jiān)測(cè)和模擬分析工作，BELEM等[11]在Doyon金礦的待充填采場(chǎng)內(nèi)預(yù)先布置傳感器，監(jiān)測(cè)分析了充填過程中及充填完成后膏體料漿孔隙水壓力和總應(yīng)力隨時(shí)間變化過程；THOMPSON等[12]以及DOHERTY等[13]也相繼開展了類似的現(xiàn)場(chǎng)監(jiān)測(cè)試驗(yàn)，對(duì)特定礦山的充填料漿固結(jié)過程進(jìn)行了研究。然而，現(xiàn)場(chǎng)監(jiān)測(cè)存在成本高、數(shù)據(jù)離散性大、不確定影響因素干擾多等問題，且只能得到給定礦山充填料漿的固結(jié)特性，料漿固結(jié)分析的可重復(fù)性不強(qiáng)，但也為揭示實(shí)際采場(chǎng)充填料漿固結(jié)過程的力學(xué)機(jī)理提供了重要數(shù)據(jù)支撐。而數(shù)值模擬方法具備成本低且效率高、模擬采場(chǎng)靈活性好、結(jié)果可重復(fù)性強(qiáng)等優(yōu)勢(shì)，一些學(xué)者主要采用數(shù)值模擬手段研究充填料漿固結(jié)特性，ELMKADMI等[14]基于SIGMA/W軟件數(shù)值模擬了采場(chǎng)充填料漿的排水固結(jié)過程，研究了充填速度、采場(chǎng)尺寸等因素對(duì)應(yīng)力分布狀態(tài)的影響；LI等[15]基于GeoStudio軟件模擬了采場(chǎng)脫水管對(duì)連續(xù)充填過程中料漿排水的影響；JAOUHAR等[16]基于SIGMA/W軟件模擬了采場(chǎng)連續(xù)充填過程，分析了充填體固結(jié)狀態(tài)對(duì)采場(chǎng)擋墻荷載的影響。然而，在國(guó)內(nèi)充填采礦領(lǐng)域普遍用于充填體力學(xué)問題計(jì)算的FLAC3D軟件，僅有VEENSTRA[17]基于此軟件研究了早期充填料漿的應(yīng)力狀態(tài)發(fā)展過程。

考慮國(guó)內(nèi)礦山推行使用全尾砂充填骨料的發(fā)展趨勢(shì)，以及對(duì)采場(chǎng)全尾砂充填料漿固結(jié)機(jī)理認(rèn)識(shí)欠缺的問題，有必要探究非膠尾砂充填料漿的固結(jié)過程。為此，本文提出一種基于FLAC3D軟件的采場(chǎng)非膠尾砂充填料漿固結(jié)過程數(shù)值模擬方法，以孔隙水壓力為定量表征參數(shù)，對(duì)比了尾砂料漿固結(jié)過程的數(shù)值模擬結(jié)果和Gibson理論解析解，校核了FLAC3D軟件模擬采場(chǎng)充填料漿固結(jié)過程的可靠性，并給出了FLAC3D數(shù)值模型中網(wǎng)格參數(shù)的選取方法，進(jìn)而示例模擬了某礦采場(chǎng)尾砂料漿應(yīng)力狀態(tài)隨固結(jié)時(shí)間演變過程，與該礦采場(chǎng)充填料漿孔隙水壓力和豎向總應(yīng)力的實(shí)際監(jiān)測(cè)數(shù)據(jù)進(jìn)行對(duì)比驗(yàn)證，進(jìn)一步校驗(yàn)了該數(shù)值模擬方法計(jì)算實(shí)際采場(chǎng)尾砂料漿固結(jié)參數(shù)的有效性。 本文研究可為采場(chǎng)非膠尾砂充填料漿自重固結(jié)研究提供一種可靠的數(shù)值計(jì)算方法。

1 采場(chǎng)充填料漿固結(jié)過程解析法

1.1 Gibson模型

采場(chǎng)空區(qū)實(shí)際充填過程中，料漿高度隨時(shí)間連續(xù)地均勻增大，與此同時(shí)充填料漿在自重作用下發(fā)生脫水固結(jié)，通常使用孔隙水壓力定量表征充填料漿固結(jié)過程。在礦山采場(chǎng)尾砂充填領(lǐng)域，目前沒有專用于描述尾砂充填料漿孔隙水壓力隨固結(jié)時(shí)間變化的理論模型，對(duì)此可借鑒土力學(xué)領(lǐng)域Gibson提出的飽和土體固結(jié)模型及其孔隙水壓力解析計(jì)算方法。

圖1展示了Gibson一維固結(jié)模型示意圖，假設(shè)存在厚度隨時(shí)間連續(xù)均勻增大的飽和土體，其在水平方向無限延伸，并在自重作用下脫水固結(jié)，圖示時(shí)刻土體厚度為h(h=m×t，m為土體厚度均勻增大的速率，t為時(shí)間)，最終可增加至H，在此過程中土體的頂面和底面保持為透水邊界。

圖1 Gibson一維固結(jié)模型示意圖Fig.1 Gibson’s one-dimensional consolidation model

為得到土體內(nèi)任一位置z在任一時(shí)刻t的孔隙水壓力，Gibson模型提出滲流固結(jié)控制方程，見式(1)。

(1)

解式(1)可得孔隙水壓力pw，見式(2)。

(2)

式中：pw為孔隙水壓力，kPa；z為土體內(nèi)任一點(diǎn)標(biāo)高(0≤z≤h)，m；γ為土體飽和容重，kN/m3；cv為固結(jié)系數(shù)(m2/s)；h為當(dāng)前土體厚度，m；t為時(shí)間(t=0時(shí)，h=0)，s；m為土體厚度增加速率(常數(shù))；ξ為積分變量。

1.2 充填料漿孔隙水壓力隨固結(jié)時(shí)間變化解析

位于西澳大利亞的Kanowna Belle(KB)金礦采用膏體充填，采場(chǎng)長(zhǎng)18 m、寬15 m、高40 m[18]。為評(píng)估該礦山采場(chǎng)充填料漿孔隙水壓力隨固結(jié)時(shí)間的變化過程，嘗試選用Gibson一維固結(jié)模型計(jì)算孔隙水壓力分布，其中須假設(shè)料漿頂面和底面均為透水邊界，僅在自重作用下脫水固結(jié)且不考慮水平方向的位移和滲流。該礦山采場(chǎng)充填使用其選廠產(chǎn)出的全尾砂，表1為充填料漿的力學(xué)參數(shù)和滲流參數(shù)[18]，用于Gibson模型解析計(jì)算及后續(xù)數(shù)值模擬。此外，其他用于模型計(jì)算的常量參數(shù)也列于表中，包括水密度、水壓縮模量和重力加速度。

表1 KB金礦膏體充填料漿的力學(xué)參數(shù)和滲流參數(shù)Table 1 Mechanical and fluid-flow parameters ofpaste backfill from KB gold mine

續(xù)表1

根據(jù)式(3)和式(4)和表1數(shù)據(jù)計(jì)算得出充填體飽和容重γ和固結(jié)系數(shù)cv：γ=17.2 kN/m3，cv=9.67×10-5m2/s。

γ=(ρd+nsρf)g

(3)

(4)

采場(chǎng)實(shí)際充填時(shí)，料漿高度隨時(shí)間的增長(zhǎng)速率與礦山充填系統(tǒng)的充填能力、待充采場(chǎng)的水平截面面積有關(guān)，穩(wěn)定充填時(shí)采場(chǎng)料漿高度隨時(shí)間近似線性變化，因此，可用式(2)Gibson模型中的m表征采場(chǎng)充填料漿的高度增長(zhǎng)速率。考慮到該礦山充填能力80 m3/h和采場(chǎng)水平截面面積270 m2，計(jì)算得料漿增長(zhǎng)速率m=0.3 m/h=8.3×10-5m/s。將γ=17.2 kN/m3，cv=9.67×10-5m2/s和m=8.3×10-5m/s代入式(2)，可解析計(jì)算得到不同時(shí)刻充填料漿孔隙水壓力在高度上的分布圖，如圖2所示，圖2中橫坐標(biāo)、縱坐標(biāo)分別為正則化的孔隙水壓力和標(biāo)高。不同時(shí)刻用無量綱參數(shù)時(shí)間因數(shù)T表征(式(5))，依據(jù)慣例，選取T=1、2、4、8、16和32等值(也可根據(jù)需要選取64、128等更大的值)。時(shí)間因數(shù)T量化了以平均速率m將給定的料漿(固結(jié)特性由cv表征)充入采場(chǎng)空區(qū)且料漿高度等于h時(shí)采場(chǎng)的充填狀態(tài)，當(dāng)料漿增長(zhǎng)速率m和固結(jié)系數(shù)cv為常數(shù)時(shí)，T值越大對(duì)應(yīng)采場(chǎng)料漿高度h越大。

(5)

由圖2可知，由于料漿頂面和底面為透水邊界，兩處的孔隙水壓力始終為零。另外，采場(chǎng)連續(xù)充填過程中，隨時(shí)間增加，孔隙水壓力逐漸增大，在采場(chǎng)下部區(qū)域(z/h=0.3)孔隙水壓力增幅尤其明顯，這與現(xiàn)場(chǎng)監(jiān)測(cè)的孔隙水壓力變化趨勢(shì)吻合，表明Gibson一維固結(jié)模型理論的解析計(jì)算解可近似反映采場(chǎng)內(nèi)孔隙水壓力隨固結(jié)過程變化。Gibson一維自重固結(jié)理論適用于描述厚度逐漸增大的飽和土體的自重固結(jié)過程，多用于港口、航道和海洋工程中淤泥等類型土體堆場(chǎng)固結(jié)研究。由圖2可知，該理論模型也可在一定假設(shè)條件下解析計(jì)算礦山采場(chǎng)非膠尾砂充填料漿的孔隙水壓力隨自重固結(jié)過程的演變規(guī)律，但仍存在一定局限性。因?yàn)槲采俺涮盍蠞{在三維尺寸有限的采場(chǎng)空區(qū)內(nèi)發(fā)生固結(jié)時(shí)，采場(chǎng)側(cè)壁圍巖、排水管、充填擋墻等均可作為排水邊界，料漿內(nèi)自由水能夠在水平方向和豎直方向發(fā)生滲流，這與Gibson固結(jié)模型中僅設(shè)定飽和土頂面和底面為透水邊界、滲流僅發(fā)生在豎直方向有一定差別；另外，采場(chǎng)充填料漿自重固結(jié)過程中會(huì)發(fā)生水平位移和豎直位移。因此，不能將采場(chǎng)充填料漿自重固結(jié)過程計(jì)算簡(jiǎn)化為荷載方向、滲流方向和位移方向保持一致的一維固結(jié)問題。

圖2 不同時(shí)刻尾砂料漿孔隙水壓力在高度上的分布圖(Gibson的解析解)Fig.2 Values of pore water pressure along height atdifferent time (Gibson’s solution)

為克服Gibson固結(jié)模型理論用于采場(chǎng)尾砂充填料漿固結(jié)過程計(jì)算的局限性，可采用數(shù)值計(jì)算方法模擬實(shí)際采場(chǎng)的固結(jié)過程，但需要首先校核所用數(shù)值模擬方法的可靠性。為此，基于FLAC3D軟件，首先數(shù)值模擬充填料漿的一維自重固結(jié)過程，并以Gibson模型解析解為標(biāo)準(zhǔn)，校核驗(yàn)證FLAC3D模擬計(jì)算充填料漿固結(jié)過程的準(zhǔn)確性和可靠性。

2 采場(chǎng)充填料漿固結(jié)過程數(shù)值模擬法

2.1 采場(chǎng)充填固結(jié)數(shù)值模型

采場(chǎng)空區(qū)實(shí)際充填過程中，料漿高度隨時(shí)間連續(xù)地均勻增大，但是在數(shù)值模擬中，受網(wǎng)格單元限制，模型無法連續(xù)變化達(dá)到最終充填高度，只能近似模擬料漿高度隨時(shí)間連續(xù)增長(zhǎng)的實(shí)際充填過程。

如圖3(a)所示，將實(shí)際充填過程簡(jiǎn)化為間斷式地由下至上“添加”一系列分層。在FLAC3D軟件內(nèi)建立預(yù)定高度的充填體模型后，首先將整個(gè)模型的本構(gòu)關(guān)系選擇為“空模型”，此時(shí)相當(dāng)于未充填的采場(chǎng)采空區(qū)。之后，在t=t1時(shí)刻，將模型最底部厚度為h0的分層重新選擇為線彈性模型，相當(dāng)于加載了一個(gè)分層，緊接著令現(xiàn)有的h0充填料漿模型計(jì)算t0時(shí)間，到達(dá)t=t2時(shí)刻，在上一時(shí)步的模型頂部再加載一厚度為h0的分層，同樣計(jì)算t0時(shí)間，到達(dá)t=t3時(shí)刻。重復(fù)進(jìn)行這一步驟，直至模型達(dá)到預(yù)定的采場(chǎng)充填高度，這就近似地模擬了以平均充填速率m=h0/t0的采場(chǎng)采空區(qū)充填過程。 通過調(diào)整t0，可以模擬不同充填速率及各層充填料漿脫水固結(jié)時(shí)間。

按照上述模擬方法，在FLAC3D軟件內(nèi)建立圖3(b)所示的一維充填體模型，其中X軸、Y軸沿水平方向，Z軸沿豎直方向，模型高度為H，將其劃分為一系列立方體網(wǎng)格單元，邊長(zhǎng)為L(zhǎng)z。

圖3 充填料漿自重固結(jié)過程模型化處理Fig.3 Modeling of self-weight consolidation of backfill

2.2 模型參數(shù)及邊界條件

基于FLAC3D軟件模擬采場(chǎng)充填料漿的自重固結(jié)過程，需分別為力學(xué)過程和滲流過程計(jì)算指定本構(gòu)關(guān)系和材料參數(shù)。本次數(shù)值模擬的力學(xué)過程本構(gòu)關(guān)系選為線彈性模型，滲流過程選為各向同性滲流模型。數(shù)值模擬的力學(xué)及滲流過程所需的材料參數(shù)從表1選取。為方便建模，假設(shè)充填料漿增長(zhǎng)速率m在數(shù)值上等于cv，即m=9.67×10-5m/s。 在m=cv條件下，根據(jù)式(5)時(shí)間因數(shù)的定義：T=hm/cv=h，此時(shí)，時(shí)間因數(shù)T和對(duì)應(yīng)的充填料漿模型高度h在數(shù)值上相等。

前文已選定時(shí)間因數(shù)T的典型值，即T=1、2、4、8、16和32，表征不同時(shí)刻的采場(chǎng)空區(qū)充填狀態(tài)。開展采場(chǎng)空區(qū)充填料漿固結(jié)過程的模擬計(jì)算時(shí)，模型高度通過網(wǎng)格加載由0逐漸增大，在分別達(dá)到h=1 m、2 m、4 m、8 m、16 m和32 m時(shí)，提取并對(duì)比相應(yīng)時(shí)刻采場(chǎng)充填料漿的孔隙水壓力，直至達(dá)到采場(chǎng)充填總高度H。 本次為實(shí)現(xiàn)與Gibson模型理論解析解對(duì)比，選定采場(chǎng)充填料漿模型高度H最大為32 m。

圖4為充填料漿一維固結(jié)模型的邊界條件。 對(duì)于力學(xué)邊界條件，模型底面限制水平和豎直方向位移，側(cè)面限制水平方向位移。對(duì)于滲流邊界條件，每次加載新充填分層后，模型頂面、底面設(shè)置為透水邊界(孔隙水壓力為0)，側(cè)面為不透水邊界，與Gibson模型的排水邊界保持一致。需要注意的是，每次加載新充填分層后，充填料漿模型的頂面位置上移，在設(shè)置頂面透水邊界時(shí)，需刪除上一時(shí)刻的頂面透水邊界。

圖4 數(shù)值模型的邊界條件Fig.4 Boundary conditions for numerical model

2.3 網(wǎng)格尺寸的影響研究

為了使數(shù)值模擬的充填過程更接近于實(shí)際情況，圖3中每個(gè)加載時(shí)步添加的充填料漿分層厚度h0應(yīng)盡量取較小值，但h0受限于網(wǎng)格單元尺寸不能無限減小，而且如果h0取為一個(gè)單元大小(即h0=Lz)，F(xiàn)LAC3D軟件在模擬計(jì)算過程中會(huì)報(bào)錯(cuò)。因此，每個(gè)加載時(shí)步添加的充填料漿分層厚度h0的最小值為兩個(gè)網(wǎng)格單元高度，即h0=2Lz。在此基礎(chǔ)上， 只需開展網(wǎng)格單元高度Lz取值的敏感性分析，就能相應(yīng)得出數(shù)值模擬充填料漿加載分層厚度h0的可靠值，實(shí)現(xiàn)充填料漿固結(jié)過程模擬結(jié)果數(shù)值穩(wěn)定，為此，在選定采場(chǎng)充填總高度H=32 m的前提下，取網(wǎng)格單元尺寸Lz=1 m、0.5 m、0.25 m、0.1 m，分別在FLAC3D軟件內(nèi)建模分析充填料漿固結(jié)過程，并將數(shù)值模擬結(jié)果與Gibson模型理論解析解對(duì)比，探究模型網(wǎng)格尺寸對(duì)料漿固結(jié)過程模擬結(jié)果的影響。

圖5對(duì)比了不同Lz取值時(shí)FLAC3D軟件模擬的和Gibson模型理論解析計(jì)算的不同時(shí)刻充填料漿孔隙水壓力沿高度分布狀態(tài)。從圖5可以看出，隨網(wǎng)格尺寸減小，孔隙水壓力的數(shù)值解逐漸趨近于解析解。當(dāng)T=32(或充填高度h=32 m)時(shí)，即使網(wǎng)格尺寸Lz取較大值(如1 m)時(shí)數(shù)值解也能與解析解很好地吻合，并且減小Lz，數(shù)值解不發(fā)生明顯變化，這表明模擬結(jié)果已達(dá)到了數(shù)值穩(wěn)定。時(shí)間因數(shù)T=16、8、4時(shí)，數(shù)值解與解析解均有一定差異，但隨Lz減小，兩者的差距逐漸縮小。 值得注意的是，從圖5(d)可以看出，在T=1或2時(shí)，選取了相對(duì)較小的網(wǎng)格尺寸Lz=0.1 m，數(shù)值解與解析解仍有較大差距，表明此網(wǎng)格單元尺寸構(gòu)建的數(shù)值模型無法模擬得出準(zhǔn)確的固結(jié)狀態(tài)，還需進(jìn)一步網(wǎng)格細(xì)化。

為獲得T=1或2時(shí)準(zhǔn)確的數(shù)值解，繼續(xù)加密網(wǎng)格，使Lz=0.025 m、0.010 m和0.005 m，模擬完成后對(duì)比T=1時(shí)的FLAC3D模擬結(jié)果和Gibson模型理論解析解，如圖6所示，減小Lz至0.01 m后，數(shù)值解與解析解吻合度較高，但Lz從0.01 m再降低至0.005 m時(shí)，模擬結(jié)果不發(fā)生明顯變化，表明充填料漿固結(jié)過程的模擬結(jié)果已達(dá)到了數(shù)值穩(wěn)定。

圖6 T=1時(shí)，Lz取值對(duì)模擬結(jié)果的影響Fig.6 Influence of Lz on the simulation results when T=1

前述分析表明，采用提出的FLAC3D數(shù)值模擬方法，選取合適網(wǎng)格尺寸，可以得到充填料漿固結(jié)參數(shù)的準(zhǔn)確模擬結(jié)果，而網(wǎng)格尺寸對(duì)模擬結(jié)果的影響與所關(guān)注的充填料漿高度有關(guān)。 在圖5中，當(dāng)所關(guān)注充填料漿高度較大時(shí)(如T=32時(shí)，h=32 m)，即使網(wǎng)格尺寸相對(duì)較大(如Lz=1 m)，模擬結(jié)果也能與Ginson模型理論解析解相吻合。 然而，當(dāng)所關(guān)注的充填料漿高度較小時(shí)(如T=1時(shí)，h=1 m)，網(wǎng)格尺寸須根據(jù)當(dāng)前的模型高度設(shè)定得足夠小(如Lz=0.01 m)，才能確保模擬結(jié)果達(dá)到準(zhǔn)確的Gibson模型理論解析解。因此，如何根據(jù)重點(diǎn)關(guān)注的采場(chǎng)充填料漿高度合理選取網(wǎng)格尺寸十分重要，為確保充填料漿固結(jié)參數(shù)模擬結(jié)果的準(zhǔn)確性并優(yōu)化計(jì)算效率，分析推薦數(shù)值模型網(wǎng)格尺寸的優(yōu)化選取方法如下所述。

圖5 不同Lz取值的模擬結(jié)果與Gibson解析計(jì)算解對(duì)比Fig.5 Comparison between the simulation results of different Lz and Gibson’s solution

2.4 模型網(wǎng)格尺寸選取方法

在采場(chǎng)充填料漿固結(jié)過程模擬分析時(shí)，根據(jù)重點(diǎn)關(guān)注的充填料漿高度h與模型網(wǎng)格尺寸Lz的比值h/Lz，可平衡優(yōu)化數(shù)值模型網(wǎng)格劃分方法與網(wǎng)格尺寸?；跁r(shí)間因數(shù)T表征的采場(chǎng)充填狀態(tài)(不同時(shí)刻采場(chǎng)達(dá)到的料漿高度)，表2給出了前述模擬結(jié)果對(duì)應(yīng)的h/Lz值。

對(duì)比表2和圖5，當(dāng)h/Lz≥32時(shí)，分別能確保T=32(圖5(a))、T=16(圖5(b))、T=8(圖5(c))和T=4(圖5(d))時(shí)的FLAC3D軟件模擬結(jié)果與Gibson模型理論解析解良好吻合；而當(dāng)h/Lz<32時(shí)，相應(yīng)的模擬結(jié)果均與解析解有一定差距。對(duì)比表2和圖6，當(dāng)關(guān)注的充填料漿高度較小時(shí)，即在T=1時(shí)，如果按照h/Lz=40的網(wǎng)格模型，則模擬結(jié)果與解析解仍有一定差異；但當(dāng)h/Lz=100和200時(shí)，模擬結(jié)果與解析解吻合度較高。 表明關(guān)注的充填料漿模型高度越小時(shí)，數(shù)值計(jì)算結(jié)果對(duì)網(wǎng)格尺寸越敏感?；谏鲜龇治隹芍?，當(dāng)數(shù)值模型網(wǎng)格劃分滿足h/Lz≥32時(shí)，可認(rèn)為模型的網(wǎng)格劃分是足夠致密的，能夠模擬得出可靠穩(wěn)定的數(shù)值解，但如果所關(guān)注的充填料漿模型高度較小或關(guān)注采場(chǎng)充填料漿早期的固結(jié)狀態(tài)時(shí)，為獲得相對(duì)準(zhǔn)確的模擬結(jié)果，應(yīng)使h/Lz不小于100。

表2 不同時(shí)間因數(shù)T對(duì)應(yīng)的數(shù)值模型高度h與Lz的比值Table 2 Ratio of height h and Lz of the numericalmodel corresponding to different time factor

圖7對(duì)比了Lz=0.01 m時(shí)(h/Lz≥100) FLAC3D模擬的和Gibson模型理論解析計(jì)算的不同時(shí)刻充填料漿孔隙水壓力沿高度分布狀態(tài)，可以看出數(shù)值解和解析解良好吻合。Lz=0.01 m時(shí)，在所有關(guān)注的充填狀態(tài)中，T=1時(shí)刻的充填料漿模型高度最小(即h=1 m)，此時(shí)h/Lz=100，而其他時(shí)刻的h/Lz分別為200、400、800、1 600、3 200，均大于100。最終表明，采用本文提出的方法模擬充填料漿固結(jié)過程，充填體模型網(wǎng)格劃分時(shí)，將自重固結(jié)方向的網(wǎng)格尺寸設(shè)定為此時(shí)刻充填料漿高度的1/100以內(nèi)，就能獲得準(zhǔn)確的模擬計(jì)算結(jié)果。

圖7 Lz=0.01 m時(shí)FLAC3D模擬結(jié)果與解析解對(duì)比Fig.7 Comparison between simulation results andGibson’s solution when Lz=0.01 m

3 采場(chǎng)充填料漿固結(jié)過程數(shù)值模擬方法工程應(yīng)用

前文基于FLAC3D軟件和Gibson模型理論，從一維角度提出并驗(yàn)證了采場(chǎng)非膠尾砂充填料漿固結(jié)過程的可靠數(shù)值模擬方法，然而充填料漿在三維尺寸的采場(chǎng)內(nèi)發(fā)生固結(jié)時(shí)不能簡(jiǎn)單地簡(jiǎn)化為一維固結(jié)問題，因此，有必要進(jìn)一步評(píng)估所提出的數(shù)值方法模擬實(shí)際三維采場(chǎng)非膠充填料漿固結(jié)過程的有效性。為此，采用本文的數(shù)值模擬方法，模擬計(jì)算了KB金礦采場(chǎng)充填料漿的孔隙水壓力和豎向總應(yīng)力隨固結(jié)時(shí)間的演變過程，再與礦山采場(chǎng)實(shí)際監(jiān)測(cè)的相應(yīng)應(yīng)力值進(jìn)行綜合對(duì)比，進(jìn)一步校驗(yàn)評(píng)估本文提出的數(shù)值模擬方法的有效性。

3.1 工程背景

為研究礦山實(shí)際采場(chǎng)充填料漿的孔隙水壓力和總應(yīng)力隨固結(jié)時(shí)間的變化，HELINSKI等[18]在KB金礦某典型待充采場(chǎng)底板中心位置布設(shè)了孔隙水壓計(jì)和土壓力傳感器，開展充填料漿固結(jié)過程中應(yīng)力狀態(tài)演變現(xiàn)場(chǎng)監(jiān)測(cè)，該采場(chǎng)尺寸及現(xiàn)場(chǎng)監(jiān)測(cè)點(diǎn)位如圖8(a)所示。 采場(chǎng)采用全尾砂充填，灰砂比為1∶32，全尾砂料漿的力學(xué)參數(shù)和滲流參數(shù)見表1。

圖8 KB金礦采場(chǎng)示意圖及對(duì)應(yīng)的數(shù)值模型Fig.8 Stope diagram and corresponding numericalmodel of KB gold mine

為規(guī)避連續(xù)充填時(shí)采場(chǎng)充填擋墻承受過大側(cè)向推力而失穩(wěn)，KB金礦實(shí)際采用間斷式充填，即先充填10 m，隔1 d后，再連續(xù)充填剩余空區(qū)至40 m。圖9展示了該采場(chǎng)空區(qū)實(shí)際充填過程的3個(gè)時(shí)段，即：0～34.5 h時(shí)段，連續(xù)充填采場(chǎng)空區(qū)底部10 m高度， 邊充填邊排水；34.5～58.5 h時(shí)段，停止充填，等待24 h使第一時(shí)段已充的10 m料漿固結(jié)排水；58.5～184 h時(shí)段，連續(xù)充填剩余采空區(qū)高度直至充滿采場(chǎng)總高度40 m，邊充填邊排水。

圖9 KB金礦采場(chǎng)充填過程Fig.9 Filling sequence of KB gold mine

3.2 采場(chǎng)充填料漿固結(jié)過程數(shù)值模擬

假定采場(chǎng)滿足平面應(yīng)變條件，在FLAC3D軟件內(nèi)建立采場(chǎng)充填料漿固結(jié)數(shù)值模型，如圖8(b)所示，其中X軸沿采場(chǎng)寬度方向、Z軸沿采場(chǎng)高度方向。依照前文所述的網(wǎng)格尺寸選取方法，以實(shí)際采場(chǎng)第一時(shí)段完成的充填高度h=10 m為基準(zhǔn)劃分網(wǎng)格，則數(shù)值模型中豎直方向的網(wǎng)格單元尺寸取為0.1 m，確保Z方向單元尺寸與此時(shí)刻充填料漿模型高度之比等于100。數(shù)值模擬計(jì)算的力學(xué)本構(gòu)關(guān)系選為莫爾庫(kù)倫模型，滲流本構(gòu)關(guān)系選為各向同性的滲流模型，力學(xué)參數(shù)和滲流參數(shù)依據(jù)表1選取。對(duì)于力學(xué)邊界條件，除頂面外其余邊界面均限制水平和豎直方向位移。對(duì)于滲流邊界條件，模型頂面和充填擋墻處設(shè)置為透水邊界，其他邊界為不透水邊界。采用提出的方法模擬采場(chǎng)充填料漿的固結(jié)過程，并提取充填過程中不同固結(jié)時(shí)刻采場(chǎng)底板中心處的孔隙水壓力和豎向總應(yīng)力，模擬完成后，綜合對(duì)比數(shù)值模擬結(jié)果和現(xiàn)場(chǎng)監(jiān)測(cè)結(jié)果。

3.3 數(shù)值模擬與現(xiàn)場(chǎng)監(jiān)測(cè)結(jié)果對(duì)比驗(yàn)證

圖10對(duì)比了KB金礦實(shí)際采場(chǎng)充填料漿固結(jié)過程孔隙水壓力與豎向總應(yīng)力的數(shù)值模擬結(jié)果和現(xiàn)場(chǎng)監(jiān)測(cè)值，分析該采場(chǎng)充填的3段充填固結(jié)過程可知。

圖10 孔隙水壓力與豎向總應(yīng)力的模擬結(jié)果和現(xiàn)場(chǎng)監(jiān)測(cè)結(jié)果對(duì)比Fig.10 Comparison of modeling results and monitoredvalues of PWP and vertical total stress

1) 0～34.5 h時(shí)段，數(shù)值模擬結(jié)果與現(xiàn)場(chǎng)監(jiān)測(cè)值吻合度較高，兩應(yīng)力參數(shù)(孔隙水壓力和豎向總應(yīng)力)均隨時(shí)間逐漸增大，且豎向總應(yīng)力始終大于孔隙水壓力，表明充填料漿內(nèi)已開始發(fā)展有效應(yīng)力。數(shù)值模擬結(jié)果和現(xiàn)場(chǎng)監(jiān)測(cè)值吻合度較高，主要是因?yàn)镵B金礦膏體充填料漿的膠凝材料含量較低(灰砂比=1∶32)，水化反應(yīng)短期內(nèi)(34.5 h)對(duì)料漿固結(jié)過程影響有限，該膏體充填料漿表現(xiàn)出近似非膠結(jié)充填料漿的固結(jié)特性。這表明本文提出的數(shù)值模擬方法可以有效模擬計(jì)算實(shí)際采場(chǎng)非膠充填料漿固結(jié)過程中的應(yīng)力演變。

2) 34.5～58.5 h時(shí)段，采場(chǎng)充填料漿高度保持10 m不變，可以看出孔隙水壓力和豎向總應(yīng)力隨時(shí)間逐漸減小，其中，孔隙水壓力的數(shù)值模擬結(jié)果與現(xiàn)場(chǎng)監(jiān)測(cè)值能夠較好吻合；豎向總應(yīng)力均表現(xiàn)出下降趨勢(shì)，但數(shù)值模擬結(jié)果下降幅度小于現(xiàn)場(chǎng)監(jiān)測(cè)值下降幅度，該差異來源于實(shí)際充填料漿中膠凝材料水化產(chǎn)物的影響。

膠凝材料的水化產(chǎn)物可填充尾砂顆粒間隙，細(xì)化孔隙結(jié)構(gòu)，引起充填體剛度、強(qiáng)度隨水化反應(yīng)逐漸增大，這使得阻止充填料漿固結(jié)過程中下沉趨勢(shì)的能力提高，進(jìn)而降低了充填體內(nèi)的豎向總應(yīng)力。然而在數(shù)值模擬中，不考慮膠凝材料水化反應(yīng)對(duì)充填體強(qiáng)度、剛度的影響，充填體的力學(xué)參數(shù)保持為初始充填狀態(tài)的固定值，導(dǎo)致模擬的充填料漿下沉趨勢(shì)更明顯，故豎向總應(yīng)力較大。數(shù)值模擬結(jié)果與現(xiàn)場(chǎng)監(jiān)測(cè)值的差異，也反映了本文提出的采場(chǎng)充填料漿固結(jié)過程模擬方法主要適用于非膠充填料漿，后續(xù)還須進(jìn)一步研究考慮膠凝材料水化反應(yīng)的充填料漿固結(jié)過程數(shù)值模擬方法。

3) 58.5～184 h時(shí)段，從采場(chǎng)剩余30 m高度范圍的連續(xù)充填過程可看出，數(shù)值模擬結(jié)果和現(xiàn)場(chǎng)監(jiān)測(cè)值的差異隨時(shí)間逐漸增大?，F(xiàn)場(chǎng)監(jiān)測(cè)的豎向總應(yīng)力在58.5～85.0 h時(shí)段逐漸增大，達(dá)到峰值后隨時(shí)間逐漸減小，在第120 h后趨于穩(wěn)定，孔隙水壓力監(jiān)測(cè)結(jié)果也表現(xiàn)出與豎向總應(yīng)力相似的變化趨勢(shì)，且在120 h后采場(chǎng)底部監(jiān)測(cè)的孔隙水壓力接近于零。然而，數(shù)值模擬的兩應(yīng)力參數(shù)均隨時(shí)間持續(xù)增大。該差異的主要原因是：在實(shí)際采場(chǎng)中，經(jīng)過長(zhǎng)時(shí)間的自重固結(jié)過程，并受膠凝材料水化反應(yīng)影響，采場(chǎng)底部充填體的強(qiáng)度、剛度逐漸發(fā)展，滲透系數(shù)逐漸減小，形成具有一定力學(xué)承載能力的致密固結(jié)體，在其頂部新充入的料漿對(duì)采場(chǎng)底部監(jiān)測(cè)點(diǎn)應(yīng)力狀態(tài)的影響逐漸減弱，因此現(xiàn)場(chǎng)監(jiān)測(cè)的兩應(yīng)力值表現(xiàn)出先增大、后減小、最終趨于穩(wěn)定的變化趨勢(shì)。然而在數(shù)值模擬中，不考慮膠凝材料對(duì)固結(jié)過程的影響，且充填料漿數(shù)值模型的剛度和滲透系數(shù)等參數(shù)設(shè)置為定值，新充入料漿產(chǎn)生的影響，經(jīng)過一定時(shí)間固結(jié)計(jì)算便會(huì)從充填料漿頂部傳導(dǎo)至采場(chǎng)底部，故底部應(yīng)力狀態(tài)持續(xù)受充填過程的影響，因此數(shù)值模擬的采場(chǎng)底部應(yīng)力值持續(xù)增大。

通過綜合對(duì)比KB金礦典型采場(chǎng)充填固結(jié)過程的數(shù)值模擬和現(xiàn)場(chǎng)監(jiān)測(cè)結(jié)果，可看出在0～58.5 h時(shí)段兩者吻合度高，進(jìn)一步驗(yàn)證了本文所述方法模擬采場(chǎng)非膠充填料漿固結(jié)過程的有效性；58.5 h之后兩者差異逐漸增大，體現(xiàn)了當(dāng)前的數(shù)值模型無法考慮膠凝材料水化反應(yīng)的局限性，本文作者及合作者將對(duì)此問題繼續(xù)開展研究。

4 結(jié) 論

1) FLAC3D軟件可模擬非膠尾砂料漿固結(jié)過程，能夠?qū)崿F(xiàn)Gibson固結(jié)理論模型問題的可靠模擬計(jì)算。

2) 網(wǎng)格尺寸對(duì)充填料漿固結(jié)過程的模擬結(jié)果影響顯著，且不同高度充填料漿固結(jié)過程的模擬結(jié)果對(duì)網(wǎng)格尺寸的敏感性不同。當(dāng)充填料漿高度較小時(shí)，其模擬結(jié)果可靠性對(duì)網(wǎng)格尺寸更敏感，應(yīng)采用更小網(wǎng)格尺寸才能獲得準(zhǔn)確模擬結(jié)果。

3) 數(shù)值模擬非膠尾砂充填料漿的固結(jié)過程，需確保充填料漿模型高度與自重固結(jié)方向網(wǎng)格單元尺寸的比值不小于100，才能實(shí)現(xiàn)充填料漿不同時(shí)刻固結(jié)參數(shù)的可靠數(shù)值模擬計(jì)算。

4) 提出的數(shù)值方法可模擬采場(chǎng)非膠尾砂充填料漿孔隙水壓力和總應(yīng)力隨固結(jié)時(shí)間變化，但是沒有考慮膠凝材料水化反應(yīng)對(duì)固結(jié)過程的影響，導(dǎo)致同等條件下模擬采場(chǎng)膠結(jié)充填料漿固結(jié)過程時(shí)，在充填早期得出的數(shù)值解和現(xiàn)場(chǎng)監(jiān)測(cè)值吻合性好，而采場(chǎng)充填后期的孔隙水壓力和豎向總應(yīng)力模擬值均偏大。