摘要
非均质地层参数的获取与表征方法是浅层地质体改造的重要依据。基于非均质地层的地质建模,表征具有随机不确定性和模糊不确定性的土体表观参数和物理力学参数。利用 Weibull 分布统计描述局部地层信息,并应用Diamond‑Square分形插值方法,结合差分盒子维计算非均质地层分形维数,进一步合理演算整体地层参数。采用PFC3D数值模拟软件,结合马尔科夫链蒙特卡洛方法采样数据,对非均质地层块石颗粒及土颗粒进行建模。通过非均质地层模拟结果,发现其破坏过程与应力应变曲线要比均质土体复杂得多,由于块石大颗粒的含量、分布特性对土样力学性质影响显著。本文建模方法有助于非均质地层参数表征研究。
关键词
非均质地层是一种由均质土体与岩块共同构成的特殊地层,主要由具有一定尺寸的强度较高岩块、强度相对较低的土体及孔隙等组
现阶段对非均质地层表征的研究主要集中在现场及室内试验和数值模拟两方
非均质地层具有复杂性、多因素性、分布的高度非均质性。有限的地层数据对土体物性参数进行描述和预测,已有研究表明运用单一的方法成效甚微,需要结合随机概率方法、随机分形技术、数字图形技
正态分布或者泊松分布在表征非均质地层时过于理想化,而韦伯分布从概率密度函数来看一般具有长尾分布,即右偏分布的特点相对来说更贴近真实地层表征。韦伯分布的概率密度函数及分布均质与方差为:
(1) |
(2) |
(3) |
式中:——随机变量;——比例参数;——形状参数;——伽马函数。
从

图1 形状参数对Weibull分布影响

图2 比例参数对Weibull分布影响
获取局部非均质地层土样,对局部土样进行表观、力学参数分析,发掘并归纳非均质土样的规律性和变异性。但对局部土样的分析归纳不能很好地表征整个非均质地层参数。为满足从局部到整体的预测,最终表征实际非均质地层整体参数值,结合插值算法,对地层获取的局部土样特征进行加密处理,在上述研究的基础上进行分形插值计算。
以地层孔隙率参数为例,首先对待表征地层分点取样,测定地层孔隙率大小并记录其所处空间位置,生成孔隙率参数的控制点,根据控制点的完整程度灵活采用四边形和三角形连线。再根据分形插值算法生成各控制点连线中点孔隙率参数,获得孔隙率参数的分形分布场。相关其他非均质地层参数的分形分布场也可用相同方式进行处理。
分形插值算法采用Diamond‑Square算法,该方法选取4个控制点构成正方形作为基本单元,Diamond在正方形对角线连线上生成中点o,中点参数取值方法如
(4) |
(5) |
式中:——网格间距;——迭代次数;——赫斯特指数;——控制点数据的标准差;——服从的高斯分布随机数。
Square取相邻4个点构成棱形,于棱形4个点连线中点再生成一个中点,计算方式同Diamond。经过1次Diamond‑Square算法后,方形数目变为原先4倍,经过次分形插值算法后,方形数目将达到原先数量的倍,Diamond-Square算法原理如

图3 Diamond‑Square算法原理
为得到分形插值算法中的重要参数赫斯特指数,还需要识别获得非均质地层的分形维数。分形维数不仅对分析、识别和解释图像纹理粗糙度有重要的意义,同时利用数据定义及试验的手段来近似计算复杂地层的不规则性的程度。引入差分盒子维DBC(differential box‑counting)计算非均质地层图像的分形维数。
差分盒子维方法首先得到非均质地层局部土样的现场照片,对已有照片进行灰度处理。


图4 图像灰度曲面示意
构成的灰度曲面后,将曲面嵌入空间划分边长的盒子集中,计算出曲面所占据的盒子数。盒子数度量灰度曲面分形维数原理见
(6) |
式中:——分割立方体边长;——含有灰度曲面格子数。

图5 曲面盒子维计量原理
在完成一种盒子边长后,再换用边长更小的盒子对灰度曲面进行度量,重复上述过程,得到与分布图,两参数所得曲线斜率即为分形维数,其精确值采用最小二乘法可得。实例地层分形维数计算结果如

图6 盒维数分形维计算结果
由地层局部信息演算出非均质地层整体表观参数后,非均质地层表观参数可表示为一系列服从随机变量的概率分布函数:
(7) |
结合前文所得表观参数的分布密度函数曲线,对于任意区间,随机变量落在该区间的概率为:
(8) |
因此非均质地层的各个表观参数在进行数值模拟计算前,可以对已有的非均质地层表观参数进行采样,将采样参数对PFC模拟单元进行赋值,从而构建出所需的PFC土样模型参数,用以后续对非均质土样物性参数求解。
采样方式选用马尔科夫链蒙特卡洛方法(MCMC)。该方法的适用性在于:当得到稳定概率分布对应的马尔科夫状态转移矩阵,对于任意的初始概率分布,经过序列转换最终可以得到稳定概率分布样本。同时该性质对于离散状态及连续状态均成立。采样方法如下:
(1)输入马尔科夫状态转移矩阵,设定达到稳定概率分布的转移次数,所需样本个数;
(2)设任意初始概率分布为,初始状态值,对每个分布:
(9) |
(3)
①条件概率分布中采样得到;
②如果,则接受转移,赋值;
③否则不接受转移,最终得到平稳分布对应的样本集。
非均质地层主要由强度较弱的黏性基粒与大颗粒不规则块石胶结构成,这种混合体的宏观力学性质不仅与地层黏粒的级配、孔隙率、块石的大小和含量有关,同时与黏粒与块石之间的胶结性质、接触本构模型以及摩擦状态密切相关。PFC程序(Particle Flow Code)是基于离散元代码的软件,允许颗粒产生平移和转动,在颗粒计算过程中自动识别新的接触,每个时间步中更新颗粒的内力、力矩、位移状态。同时颗粒球可以组合为一个刚性簇或柔性簇,用于模拟非均质地层中的大颗粒成分,最终模拟黏粒与不规则块石间的接触及其相互作用过程。
土体中黏粒与大颗粒块石间的界限是相对的。非均质地层的土体特征随研究尺度的变化固然会发生改变,地层中大颗粒块石相对于黏粒的尺寸也会随研究尺度而发生相对变化,故在对土体特征进行分析前,首先应确定研究土样范围内块石与黏粒间的阈值大小。在地层块石颗粒建模前,通过现场取样及室内筛分,初步筛选土体中大颗粒块石进行数字图像处理,导入CAD软件后,标识各颗粒平均粒径大小,通过
(10) |
式中:——非均质土中土/石阈值;——工程特征尺寸;——研究区断面面积。
将外部软件生成块石颗粒模型导入PFC3D,PFC对块石的模拟可采用刚性簇和柔性簇 2种方式,2种方法生成的块石模型均由数个小球体通过一定强度粘结在一起,区别在于刚性簇中的小球体单元没有相对位移,刚性簇只会产生平移和旋转,在外力作用下不会发生大的变形和破坏。柔性簇内小球体单元的接触情况随计算时间步不断变化更新,在外力加载足够大时,柔性簇内部颗粒的粘结将会发生破坏,这种方式在研究簇即块石破坏问题时更为适用。模拟中为探究非均质地层样整体力学特征,不考虑块石颗粒自身强度对整体的影响,故将外部模型导入PFC生成刚性簇模板集模拟土样块石。刚性簇投放过程包含以下几个步骤:
(1)确定各形态参数刚性簇数量比例及刚性簇体积占土样总体积大小;
(2)确定刚性簇模板投放时绕x轴、y轴、z轴的旋转角度,即确定各块石颗粒的方向和倾角;
(3)于整个计算区域中投放刚性簇,为简化控制块石所在土样内空间位置,块石质心位置随机;
(4)定义fish函数,搜索计算区域内所有刚性簇,对每个刚性簇赋予指针,判别其所处位置,删除模拟实验加载边墙外的刚性簇以及与边墙相交的刚性簇;
(5)赋予刚性簇质量、密度等力学参数,设定表面摩擦角及接触性质。
将块石刚性簇成功投放后,于实验加载范围内均匀投放一定粒径范围小球体,小球体尺寸遵循高斯分布。同样赋予小球密度、摩擦角、接触性质等参数,对试样整体循环少量时间步,弹开刚性簇与小球体生成后的重叠部分。因重叠量的不同,各颗粒应力状态差异很大,同时无法保证颗粒之间具备良好的接触关系。为方便后续加载实验正常进行,采用伺服机制能有效应对该情况。
对边界墙伺服一方面用于释放刚生成非均质土样中颗粒间存在的大量不可控应变能,逐渐对整体土样进行加压达到设定预加应力大小,模拟地层真实地应力环境。另一方面逐渐减弱各颗粒生成后速度失量,不断调整从而最终实现非均质土样均匀空隙率、均匀应力。伺服过程中边墙对土样施加力的过程通过动态改变边墙的法向速度实现,每个时间步得到边墙应力大小与目标应力差值,将差值乘以相应系数后得到速度赋予下一个时间步。随应力差值越小,边墙施加力也会愈发精细,最终完成伺服。三维状态下的计算原理如下:
(11) |
(12) |
式中:——边墙法向速度;——边墙当前应力值;——边墙目标应力值;——伺服参数;——边墙与颗粒接触面积;,,——墙与颗粒在x,y,z方向上接触力;,,——墙在x,y,z方向的单位向量;,,——第个边墙在x,y,z方向上的伺服速度。
PFC3D非均质模型建立流程如

图7 PFC3D非均质地层模型建立流程
PFC模拟材料的整体力学性质依靠模型中每个颗粒与颗粒、颗粒与簇之间的接触本构模型实现,根据颗粒接触特征可以将被测材料定义为线性模型、接触黏结模型、平行黏结模型和平节理模型。模拟实验中土颗粒与碎块石间的接触特性是影响土石混合体性质的重要因素,为更好地反映土石混合体的力学行为,选用接触黏结模型。以边墙作为模型边界,模型尺寸100 mm×100 mm。土体颗粒实际尺寸一般<0.075 mm,以实际土颗粒尺寸构件模型虽然对实际地层反映更为精确但计算量十分巨大,因此模型构造是土颗粒球大小进行了等比缩放,将其颗粒半径限定于1~3 mm范围,共生成约34000个颗粒模拟非均质地层土样。PFC3D模拟中生成非均质地层粗颗粒及土颗粒如

图8 生成地层粗颗粒

图9 填充地层土颗粒
PFC中颗粒间的相互作用状态可以通过接触力链来表示,

图10 地层原始接触力

图11 直剪后接触力
从

(a) 2%直剪试验

(b) 5%直剪试验

(c) 8%直剪试验

(d) 2%剪切力-位移曲线

(e) 5%剪切力-位移曲线

(f) 8%剪切力-位移曲线

(g) 10%直剪试验

(h) 15%直剪试验

(i) 20%直剪试验

(j) 10%剪切力-位移曲线

(k) 15%剪切力-位移曲线

(l) 20%剪切力-位移曲线
图12 块石颗粒不同体积比的非均质地层直剪实验对比
由
(1)对非均质地层表观参数及物性参数统计描述时,相对于正态分布曲线,Weibull分布更适用于描述地层的非均质性。采用拟合方法得到相应参数的分布函数将有助于非均质地层模型的精确构建。
(2)基于分形插值算法,采用Diamond‑Square算法对有限的地层资料进行加密处理,考虑不确定因素的影响,将局部土体信息合理扩展到整体地层。采用差分盒子维DBC(differential box‑counting)对地层土样颗粒分布的分形维数进行分析和计算,结合分形插值算法对局部数据高效利用。
(3)采用马尔科夫链蒙特卡洛方法(MCMC)对PFC3D的非均质地层土样进行赋值,构建在真实地层中随机取样模型,直剪实验结果表明,非均质地层土样与均质土样相比应力-应变曲线更为复杂,在大颗粒含量较少且变异性较大时,地层土样规律性不强,随大颗粒含量上升,剪切力曲线更为平滑。模拟结果能有效反映非均质地层的特定参数,基于本文建模思路,提供了一种非均质地层研究及表征方法。
参考文献
徐文杰,胡瑞林.土石混合体概念、分类及意义[J].水文地质工程地质,2009,36(4):50-56,70. [百度学术]
徐文杰,胡瑞林,岳中琦.基于数字图像处理的土石混合体细观结构[J].辽宁工程技术大学学报(自然科学版),2008,27(1):51-53. [百度学术]
罗荣,曾亚武.岩石矿物细胞元随机性参数赋值方法研究[J].岩土力学,2012,33(7):2221-2228. [百度学术]
吴振君.土体参数空间变异性模拟和土坡可靠度分析方法应用研究[D].北京:中国科学院大学,2009. [百度学术]
Lumb P. The variability of natural soils[J]. Canadian Geotechnical Journal,1966,3(2):74-97. [百度学术]
雷进生,武增琳,姚奇,等.基于Weibull分布的非均质地层物性参数随机模型[J].地下空间与工程学报,2017,13(3):684-691. [百度学术]
丁秀丽,张宏明,黄书岭,等.基于细观数值试验的非饱和土石混合体力学特性研究[J]. 岩石力学与工程学报,2012,31(8):1553-1566. [百度学术]
刘新荣,涂义亮,王林枫,等.土石混合体的剪切面分形特征及强度产生机制[J].岩石力学与工程学报,2017,36(9):2260-2274. [百度学术]
XU W J, XU Q, HU R L. Study on the shear strength of soil-rock mixture by large scale direct shear test[J]. International Journal of Rock Mechanics and Mining Sciences, 2011,48(8):1235-1247. [百度学术]
钟祖良,涂义亮,何晓勇,等.土石混合体物理指标及强度特性研究进展[J].地下空间与工程学报,2016,12(4):952-961. [百度学术]
车灿辉.巨厚非均质潜水含水层抽水试验及参数计算[J].探矿工程(岩土钻掘工程),2018,45(11):60-64. [百度学术]
Vanmarcke E H. Probabilistic modeling of soil profiles [J]. Journal of the Geotechnical Engineering, 1997,103(11):1227-1246. [百度学术]
Elishakoff I.Probabilistic methods in the theory of structure[M]. New York: John Wiley & Sons, Incc, 1983. [百度学术]
郭红,李波,张博.基于完全样本的两参数Weibull分布的参数估计[J].华中师范大学学报(自然科学版),2007,41(3):348-351. [百度学术]
薛亚东,岳磊,李硕标.含水率对土石混合体力学特性影响的试验研究[J].工程地质学报,2015,23(1):21-29. [百度学术]
雷进生,夏磊,王乾峰,等.非均质地层地质剖面的随机预测模型[J].地下空间与工程学报,2016,12(1):84-89. [百度学术]
唐俊伟,沈平平,黄平中.分形插值模拟渗透率及孔隙度平面分布(I)———理论分析[J].石油勘探与开发,1997,24(3):66-69. [百度学术]
张颖,陈晨,王彧佼,等.基于FLAC3D对大连某工程桩基承载力数值模拟研究[J].探矿工程(岩土钻掘工程),2019,46(5):65-71,85. [百度学术]
LI X, LIAO Q L, HE J M. In-situ test and a stochastic structural model of rock and soil aggregate in the Three Gorges reservoir area[J]. International Journal of Rock Mechanics and Mining Sciences, 2004,41(S1):702-707. [百度学术]
廖秋林,李晓,郝钊,等.土石混合体的研究现状及研究展望[J].工程地质学报,2006,14(6):800-807. [百度学术]
油新华,李晓,贺长俊.土石混合体实测结构模型的自动生成技术[J].岩土工程界,2003,6(8):60-62. [百度学术]
杨传俊.数字图像处理技术在土石混合体颗分中的应用[J].云南水力发电,2009,25(3):24-25,32. [百度学术]
舒志乐,刘新荣,刘保县.基于分形理论的土石混合体强度特征研究[J].岩石力学与工程学报,2009,28(S1):2652-2656. [百度学术]
石崇,王盛年,刘琳,等.基于数字图像分析的冰水堆积体结构建模与力学参数研究[J].岩土力学,2012,33(11):3393-3399. [百度学术]
丁秀丽,李耀旭,王新.基于数字图像的土石混合体力学性质的颗粒流模拟[J].岩石力学与工程学报,2010,29(3):477-484. [百度学术]
董启朋,卢正,詹永祥,等.土石混合体原位试验的颗粒流数值模拟分析[J].上海交通大学学报,2013,47(9):1382-1389. [百度学术]