采场非胶尾砂充填料浆自重固结过程数值模拟

2021-07-27马青海刘光生杨小聪郭利杰

中国矿业 2021年7期

马青海，刘光生，杨小聪，郭利杰

(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之后两者差异逐渐增大，体现了当前的数值模型无法考虑胶凝材料水化反应的局限性，本文作者及合作者将对此问题继续开展研究。