采場非膠尾砂充填料漿自重固結過程數值模擬

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

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

0 引 言

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

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

全尾砂細粒多，保水性好，滲透性差[4]，導致采場全尾砂料漿較難快速脫水完成固結。尤其是選用階段空場嗣后充填采礦法的大型礦山，其采場尺寸和綜合生產能力大，導致采空區充填料漿上升速度快、料漿滲流排水路徑長等問題突出，而且二步驟采場充填為節約成本往往不添加膠凝材料[5-6]，迫使礦山必須充分掌握采場全尾砂料漿的脫水固結過程。近期，國內外報道的大型采場尾砂料漿沖垮擋墻造成的安全生產事故時有發生[7]，因此，科學合理地評估采場非膠尾砂料漿的固結過程是礦山充填體力學領域的重要研究熱點之一。

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

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

1 采場充填料漿固結過程解析法

1.1 Gibson模型

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

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

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

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

(1)

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

(2)

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

1.2 充填料漿孔隙水壓力隨固結時間變化解析

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

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

續表1

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

γ=(ρd+nsρf)g

(3)

(4)

采場實際充填時，料漿高度隨時間的增長速率與礦山充填系統的充填能力、待充采場的水平截面面積有關，穩定充填時采場料漿高度隨時間近似線性變化，因此，可用式(2)Gibson模型中的m表征采場充填料漿的高度增長速率。考慮到該礦山充填能力80 m3/h和采場水平截面面積270 m2，計算得料漿增長速率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)，可解析計算得到不同時刻充填料漿孔隙水壓力在高度上的分布圖，如圖2所示，圖2中橫坐標、縱坐標分別為正則化的孔隙水壓力和標高。不同時刻用無量綱參數時間因數T表征(式(5))，依據慣例，選取T=1、2、4、8、16和32等值(也可根據需要選取64、128等更大的值)。時間因數T量化了以平均速率m將給定的料漿(固結特性由cv表征)充入采場空區且料漿高度等于h時采場的充填狀態，當料漿增長速率m和固結系數cv為常數時，T值越大對應采場料漿高度h越大。

(5)

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

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

為克服Gibson固結模型理論用于采場尾砂充填料漿固結過程計算的局限性，可采用數值計算方法模擬實際采場的固結過程，但需要首先校核所用數值模擬方法的可靠性。為此，基于FLAC3D軟件，首先數值模擬充填料漿的一維自重固結過程，并以Gibson模型解析解為標準，校核驗證FLAC3D模擬計算充填料漿固結過程的準確性和可靠性。

2 采場充填料漿固結過程數值模擬法

2.1 采場充填固結數值模型

采場空區實際充填過程中，料漿高度隨時間連續地均勻增大，但是在數值模擬中，受網格單元限制，模型無法連續變化達到最終充填高度，只能近似模擬料漿高度隨時間連續增長的實際充填過程。

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

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

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

2.2 模型參數及邊界條件

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

前文已選定時間因數T的典型值，即T=1、2、4、8、16和32，表征不同時刻的采場空區充填狀態。開展采場空區充填料漿固結過程的模擬計算時，模型高度通過網格加載由0逐漸增大，在分別達到h=1 m、2 m、4 m、8 m、16 m和32 m時，提取并對比相應時刻采場充填料漿的孔隙水壓力，直至達到采場充填總高度H。 本次為實現與Gibson模型理論解析解對比，選定采場充填料漿模型高度H最大為32 m。

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

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

2.3 網格尺寸的影響研究

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

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

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

圖6 T=1時，Lz取值對模擬結果的影響Fig.6 Influence of Lz on the simulation results when T=1

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

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

2.4 模型網格尺寸選取方法

在采場充填料漿固結過程模擬分析時，根據重點關注的充填料漿高度h與模型網格尺寸Lz的比值h/Lz，可平衡優化數值模型網格劃分方法與網格尺寸。基于時間因數T表征的采場充填狀態(不同時刻采場達到的料漿高度)，表2給出了前述模擬結果對應的h/Lz值。

對比表2和圖5，當h/Lz≥32時，分別能確保T=32(圖5(a))、T=16(圖5(b))、T=8(圖5(c))和T=4(圖5(d))時的FLAC3D軟件模擬結果與Gibson模型理論解析解良好吻合；而當h/Lz<32時，相應的模擬結果均與解析解有一定差距。對比表2和圖6，當關注的充填料漿高度較小時，即在T=1時，如果按照h/Lz=40的網格模型，則模擬結果與解析解仍有一定差異；但當h/Lz=100和200時，模擬結果與解析解吻合度較高。 表明關注的充填料漿模型高度越小時，數值計算結果對網格尺寸越敏感。基于上述分析可知，當數值模型網格劃分滿足h/Lz≥32時，可認為模型的網格劃分是足夠致密的，能夠模擬得出可靠穩定的數值解，但如果所關注的充填料漿模型高度較小或關注采場充填料漿早期的固結狀態時，為獲得相對準確的模擬結果，應使h/Lz不小于100。

表2 不同時間因數T對應的數值模型高度h與Lz的比值Table 2 Ratio of height h and Lz of the numericalmodel corresponding to different time factor

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

圖7 Lz=0.01 m時FLAC3D模擬結果與解析解對比Fig.7 Comparison between simulation results andGibson’s solution when Lz=0.01 m

3 采場充填料漿固結過程數值模擬方法工程應用

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

3.1 工程背景

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

圖8 KB金礦采場示意圖及對應的數值模型Fig.8 Stope diagram and corresponding numericalmodel of KB gold mine

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

圖9 KB金礦采場充填過程Fig.9 Filling sequence of KB gold mine

3.2 采場充填料漿固結過程數值模擬

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

3.3 數值模擬與現場監測結果對比驗證

圖10對比了KB金礦實際采場充填料漿固結過程孔隙水壓力與豎向總應力的數值模擬結果和現場監測值，分析該采場充填的3段充填固結過程可知。

圖10 孔隙水壓力與豎向總應力的模擬結果和現場監測結果對比Fig.10 Comparison of modeling results and monitoredvalues of PWP and vertical total stress

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

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

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

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

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

4 結 論

1) FLAC3D軟件可模擬非膠尾砂料漿固結過程，能夠實現Gibson固結理論模型問題的可靠模擬計算。

2) 網格尺寸對充填料漿固結過程的模擬結果影響顯著，且不同高度充填料漿固結過程的模擬結果對網格尺寸的敏感性不同。當充填料漿高度較小時，其模擬結果可靠性對網格尺寸更敏感，應采用更小網格尺寸才能獲得準確模擬結果。

3) 數值模擬非膠尾砂充填料漿的固結過程，需確保充填料漿模型高度與自重固結方向網格單元尺寸的比值不小于100，才能實現充填料漿不同時刻固結參數的可靠數值模擬計算。

4) 提出的數值方法可模擬采場非膠尾砂充填料漿孔隙水壓力和總應力隨固結時間變化，但是沒有考慮膠凝材料水化反應對固結過程的影響，導致同等條件下模擬采場膠結充填料漿固結過程時，在充填早期得出的數值解和現場監測值吻合性好，而采場充填后期的孔隙水壓力和豎向總應力模擬值均偏大。