4/5/2025, 9:57:58 PM 星期六
网刊加载中。。。

使用Chrome浏览器效果最佳,继续浏览,你可能不会看到最佳的展示效果,

确定继续浏览么?

复制成功,请在其他浏览器进行阅读

基于非均质地层参数表征的颗粒流模型构建  PDF

  • 夏余宏烨 1,2
  • 孙平贺 1,2
  • 伍新民 3
  • 盛汉洋 3
  • 王心龙 1,2
  • 张绍和 1,2
1. 有色金属成矿预测与地质环境监测教育部重点实验室,湖南 长沙 410083; 2. 中南大学地球科学与信息物理学院,湖南 长沙 410083; 3. 中南建设集团有限公司,湖南 长沙 4l0015

中图分类号: P634

最近更新:2021-12-06

DOI:10.12143/j.ztgc.2021.S1.009

  • 全文
  • 图表
  • 参考文献
  • 作者
  • 出版信息
EN
目录contents

摘要

非均质地层参数的获取与表征方法是浅层地质体改造的重要依据。基于非均质地层的地质建模,表征具有随机不确定性和模糊不确定性的土体表观参数和物理力学参数。利用 Weibull 分布统计描述局部地层信息,并应用Diamond‑Square分形插值方法,结合差分盒子维计算非均质地层分形维数,进一步合理演算整体地层参数。采用PFC3D数值模拟软件,结合马尔科夫链蒙特卡洛方法采样数据,对非均质地层块石颗粒及土颗粒进行建模。通过非均质地层模拟结果,发现其破坏过程与应力应变曲线要比均质土体复杂得多,由于块石大颗粒的含量、分布特性对土样力学性质影响显著。本文建模方法有助于非均质地层参数表征研究。

0 引言

非均质地层是一种由均质土体与岩块共同构成的特殊地层,主要由具有一定尺寸的强度较高岩块、强度相对较低的土体及孔隙等组

1-2。非均质地层的物质组成、组织结构、含水率状态和孔隙变化具有随机3,其性质也随空间和时间表现出很大的变异4-5,土颗粒及块石大颗粒物性参数相差悬殊,共同影响非均质土性质,因此该类介质在应力-应变特征、变形破坏机制等方面不同于一般的土体。这导致了非均质地层的复杂6

现阶段对非均质地层表征的研究主要集中在现场及室内试验和数值模拟两方

7-8,现场多运用大型水平推剪和直剪试验对得到土石混合体的变形特点和抗剪强度特性,而数值模拟方法则偏向土石混合体随机结构模型的自动生成技9-10,FLAC、PFC和ABAQUS等许多有限元或有限差分程序被广泛用于非均质地层分11。但对于高非均质、非均匀的介质,其结构模型的建立仍是一个难点,需要通过其他方法辅助进行。随机概率方法及随机分形技12-13被应用于土体参数的分布模型研究,其中Weibull分布被广泛应用于如孔隙率、土体强度等表征参数的分14。从而进一步研究非均质地层土性参数的分布规律及它们之间的相关15。同时为克服地层参数获取存在的“局部预测整体,整体代表局部”的局限性,分形理论被运用于建立地质变量的空间函数关系,进而建立符合地质实际的地层概念模型。已成为研究非均质性地层参数空间分布的重要工16-17。用以描述许多包括地层参数、测井数据、地震数据和岩石断裂特性等地质现象。另外对非均质地层中块石的表征Monte Carlo方法被用于模拟块石在土石混合体中的分布特征,包括块石的空间位置、含量、大小、形状和方位18-20。数码摄像、自动图像识别和灰度处理等处理技术也作为建立非均质体的精细结构力学模型的方21-22

非均质地层具有复杂性、多因素性、分布的高度非均质性。有限的地层数据对土体物性参数进行描述和预测,已有研究表明运用单一的方法成效甚微,需要结合随机概率方法、随机分形技术、数字图形技

23-24和数值模25-26等多种方法。本文采用数字图像处理技术及Weibull分布,采集和统计非均质土体中黏粒及块石颗粒的具体参数。运用随机分形理论,对有限的地层数据进行加密处理。采用马尔科夫链蒙特卡洛方法(MCMC)将采样数据投放至PFC3D颗粒流程序中,生成与真实地层相近的地层模型进行数值模拟试验,研究并表征非均质地层表观及物理力学性质。

1 非均质地层参数分布的统计描述

正态分布或者泊松分布在表征非均质地层时过于理想化,而韦伯分布从概率密度函数来看一般具有长尾分布,即右偏分布的特点相对来说更贴近真实地层表征。韦伯分布的概率密度函数及分布均质与方差为:

f(x,λ,k)=kλxλk-1e-(x/λ)k0x0x<0 (1)
F(x)=1-e-(λx)k (2)
E(x)=λΓ(1+1k)D(x)=λ2Γ(1+2k)-Γ(1+1k)2Γ(z)=0tz-1etdt (3)

式中:x——随机变量;λ——比例参数;k——形状参数;Γ——伽马函数。

图1图2可以看出,当k=1时,韦伯分布服从指数分布。当λ=1时,则称为最小化的韦伯分布。不论增加比例参数,还是增加形状参数,分布都更加类似正态分布,峰值也会趋向于随λk值的增加而增加。其中形状参数主要决定密度曲线的形状。如图1所示,当比例参数保持不变时,随形状参数的增大Weibull分布由指数分布形式逐渐过渡到正态分布形式。而比例参数则控制Weibull分布曲线的尺度大小,如图2所示,当保持形状参数不变,随形状参数的增加,Weibull分布密度函数曲线峰值右偏的同时由尖峰分布渐变为平峰分布。在形状参数和比例参数共同控制下,Weibull分布相比与正态分布形式更为多样,适用范围更广,更能有效地表征地层孔隙特征。

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

图2 比例参数对Weibull分布影响

2 非均匀地层参数的分形插值方法

获取局部非均质地层土样,对局部土样进行表观、力学参数分析,发掘并归纳非均质土样的规律性和变异性。但对局部土样的分析归纳不能很好地表征整个非均质地层参数。为满足从局部到整体的预测,最终表征实际非均质地层整体参数值,结合插值算法,对地层获取的局部土样特征进行加密处理,在上述研究的基础上进行分形插值计算。

以地层孔隙率参数为例,首先对待表征地层分点取样,测定地层孔隙率大小并记录其所处空间位置,生成孔隙率参数的控制点,根据控制点的完整程度灵活采用四边形和三角形连线。再根据分形插值算法生成各控制点连线中点孔隙率参数,获得孔隙率参数的分形分布场。相关其他非均质地层参数的分形分布场也可用相同方式进行处理。

分形插值算法采用Diamond‑Square算法,该方法选取4个控制点构成正方形作为基本单元,Diamond在正方形对角线连线上生成中点o,中点参数取值方法如公式(4),为4个角点的平均值加上服从高斯分布的随机值。

ho=(ha+hb+hc+hd)/4+Δo (4)
Δi=(d/2n)Hσ1-22H-2gauss() (5)

式中:d——网格间距;n——迭代次数;H——赫斯特指数;σ——控制点数据的标准差;gauss()——服从N(0,1)的高斯分布随机数。

Square取相邻4个点构成棱形,于棱形4个点连线中点再生成一个中点,计算方式同Diamond。经过1次Diamond‑Square算法后,方形数目变为原先4倍,经过i次分形插值算法后,方形数目将达到原先数量的4i倍,Diamond-Square算法原理如图3所示。

图3 Diamond‑Square算法原理

为得到分形插值算法中的重要参数赫斯特指数H,还需要识别获得非均质地层的分形维数。分形维数不仅对分析、识别和解释图像纹理粗糙度有重要的意义,同时利用数据定义及试验的手段来近似计算复杂地层的不规则性的程度。引入差分盒子维DBC(differential box‑counting)计算非均质地层图像的分形维数。

差分盒子维方法首先得到非均质地层局部土样的现场照片,对已有照片进行灰度处理。图4为图像灰度曲面,图像灰度曲面的复杂程度越高时相应的分形维数也越高,度量图像的粗糙程度也越大。

图4 图像灰度曲面示意

构成的灰度曲面后,将曲面嵌入空间划分边长ε的盒子集中,计算出曲面所占据的盒子数N(ε)。盒子数度量灰度曲面分形维数原理见图5。通过公式(6)即可计算出当前盒子边长ε下所得分形维数。分形维数的原理计算公式如下:

d=limε0[lgN(ε)/lg(1/ε)] (6)

式中:ε——分割立方体边长;N(ε)——含有灰度曲面格子数。

图5 曲面盒子维计量原理

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

图6 盒维数分形维计算结果

3 MCMC参数采样方法

由地层局部信息演算出非均质地层整体表观参数后,非均质地层表观参数可表示为一系列服从随机变量Xi的概率分布函数:

Fi(x)=PiXx (7)

结合前文所得表观参数的分布密度函数曲线,对于任意区间[x1,x2],随机变量Xi落在该区间的概率为:

Pix1<Xix2=PiXix2-PiXix1 (8)

因此非均质地层的各个表观参数在进行数值模拟计算前,可以对已有的非均质地层表观参数进行采样,将采样参数对PFC模拟单元进行赋值,从而构建出所需的PFC土样模型参数,用以后续对非均质土样物性参数求解。

采样方式选用马尔科夫链蒙特卡洛方法(MCMC)。该方法的适用性在于:当得到稳定概率分布对应的马尔科夫状态转移矩阵,对于任意的初始概率分布,经过序列转换最终可以得到稳定概率分布样本。同时该性质对于离散状态及连续状态均成立。采样方法如下:

(1)输入马尔科夫状态转移矩阵P,设定达到稳定概率分布的转移次数n1,所需样本个数n2

(2)设任意初始概率分布为π0(x),初始状态值x0,对每个分布πi(x)

πi(x)=πi-1(x)P=πi-2(x)P2...=π0(x)Pi (9)

(3)t=0n1+n2-1

①条件概率分布Q(xxt)中采样得到x*

②如果0,1<π(x*)Q(X*,xt),则接受转移,赋值xt+1=x*

③否则不接受转移,最终得到平稳分布对应的样本集(xn1,xn2,...,xn1+n2-1)

4 非均质地层土样的PFC模型构建

非均质地层主要由强度较弱的黏性基粒与大颗粒不规则块石胶结构成,这种混合体的宏观力学性质不仅与地层黏粒的级配、孔隙率、块石的大小和含量有关,同时与黏粒与块石之间的胶结性质、接触本构模型以及摩擦状态密切相关。PFC程序(Particle Flow Code)是基于离散元代码的软件,允许颗粒产生平移和转动,在颗粒计算过程中自动识别新的接触,每个时间步中更新颗粒的内力、力矩、位移状态。同时颗粒球可以组合为一个刚性簇或柔性簇,用于模拟非均质地层中的大颗粒成分,最终模拟黏粒与不规则块石间的接触及其相互作用过程。

土体中黏粒与大颗粒块石间的界限是相对的。非均质地层的土体特征随研究尺度的变化固然会发生改变,地层中大颗粒块石相对于黏粒的尺寸也会随研究尺度而发生相对变化,故在对土体特征进行分析前,首先应确定研究土样范围内块石与黏粒间的阈值大小。在地层块石颗粒建模前,通过现场取样及室内筛分,初步筛选土体中大颗粒块石进行数字图像处理,导入CAD软件后,标识各颗粒平均粒径大小,通过公式(10)的计算,剔除小于土/石阈值的颗粒,将粒径高于土/石阈值的大颗粒进行3D模型绘制。

dthr=0.05LcLc=A (10)

式中:dthr——非均质土中土/石阈值;Lc——工程特征尺寸;A——研究区断面面积。

将外部软件生成块石颗粒模型导入PFC3D,PFC对块石的模拟可采用刚性簇和柔性簇 2种方式,2种方法生成的块石模型均由数个小球体通过一定强度粘结在一起,区别在于刚性簇中的小球体单元没有相对位移,刚性簇只会产生平移和旋转,在外力作用下不会发生大的变形和破坏。柔性簇内小球体单元的接触情况随计算时间步不断变化更新,在外力加载足够大时,柔性簇内部颗粒的粘结将会发生破坏,这种方式在研究簇即块石破坏问题时更为适用。模拟中为探究非均质地层样整体力学特征,不考虑块石颗粒自身强度对整体的影响,故将外部模型导入PFC生成刚性簇模板集模拟土样块石。刚性簇投放过程包含以下几个步骤:

(1)确定各形态参数刚性簇数量比例及刚性簇体积占土样总体积大小;

(2)确定刚性簇模板投放时绕x轴、y轴、z轴的旋转角度,即确定各块石颗粒的方向和倾角;

(3)于整个计算区域中投放刚性簇,为简化控制块石所在土样内空间位置,块石质心位置随机;

(4)定义fish函数,搜索计算区域内所有刚性簇,对每个刚性簇赋予指针,判别其所处位置,删除模拟实验加载边墙外的刚性簇以及与边墙相交的刚性簇;

(5)赋予刚性簇质量、密度等力学参数,设定表面摩擦角及接触性质。

将块石刚性簇成功投放后,于实验加载范围内均匀投放一定粒径范围小球体,小球体尺寸遵循高斯分布。同样赋予小球密度、摩擦角、接触性质等参数,对试样整体循环少量时间步,弹开刚性簇与小球体生成后的重叠部分。因重叠量的不同,各颗粒应力状态差异很大,同时无法保证颗粒之间具备良好的接触关系。为方便后续加载实验正常进行,采用伺服机制能有效应对该情况。

对边界墙伺服一方面用于释放刚生成非均质土样中颗粒间存在的大量不可控应变能,逐渐对整体土样进行加压达到设定预加应力大小,模拟地层真实地应力环境。另一方面逐渐减弱各颗粒生成后速度失量,不断调整从而最终实现非均质土样均匀空隙率、均匀应力。伺服过程中边墙对土样施加力的过程通过动态改变边墙的法向速度实现,每个时间步得到边墙应力大小与目标应力差值,将差值乘以相应系数后得到速度赋予下一个时间步。随应力差值越小,边墙施加力也会愈发精细,最终完成伺服。三维状态下的计算原理如下:

u˙n=G(σc-σr)=GΔσσc=f+wx2f+wy2fwz2/A (11)
vix=u˙nnxviy=u˙nnyviz=u˙nnz (12)

式中:u˙n——边墙法向速度;σc——边墙当前应力值;σr——边墙目标应力值;G——伺服参数;A——边墙与颗粒接触面积;fwxfwyfwz——墙与颗粒在xyz方向上接触力;nxnynz——墙在xyz方向的单位向量;vixviyviz——第i个边墙在xyz方向上的伺服速度。

PFC3D非均质模型建立流程如图7所示。

图7 PFC3D非均质地层模型建立流程

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

图8 生成地层粗颗粒

图9 填充地层土颗粒

5 PFC模型试验结果及分析

PFC中颗粒间的相互作用状态可以通过接触力链来表示,图10图11所示为颗粒间的受力接触,力链的方向表示颗粒接触力的方向,力链越粗表示接触力越大。以不同颜色表示接触力的大小。

图10 地层原始接触力

图11 直剪后接触力

图10图11中也可以看出,非均质地层大颗粒及土颗粒生成并完成边墙伺服后,土颗粒间的接触力分布较为均匀,但在大颗粒与土颗粒、大颗粒与大颗粒、大颗粒与边墙间的接触应力状态较为集中且接触力较大。进行直剪实验后,颗粒间接触应力分布发生了改变,接触力链更多的集中于与剪切方向成约45°角的斜长带中。垂直于剪切方向上的接触力明显小于平行剪切方向上的力。图12中(a)~(l)为大颗粒块石所占体积比分别为2%、5%、8%、10%、15%、20%时数值模拟直剪实验图及相应剪切力-位移曲线。

(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 块石颗粒不同体积比的非均质地层直剪实验对比

图12各块石大颗粒体积比剪切实验对比图所示,非均质地层土样在进行直剪实验时剪切力波动值较大,特别是当大颗粒体积比较小时,土样的整体抗剪强度受大颗粒空间位置及空间方位角度影响更大,抗剪强度指标不具有规律性。当大颗粒体积比增大,土样剪切力值波动范围相对减少,最终数据也更为平稳,试验结果更加接近均质土体抗剪强度曲线。

6 结论

(1)对非均质地层表观参数及物性参数统计描述时,相对于正态分布曲线,Weibull分布更适用于描述地层的非均质性。采用拟合方法得到相应参数的分布函数将有助于非均质地层模型的精确构建。

(2)基于分形插值算法,采用Diamond‑Square算法对有限的地层资料进行加密处理,考虑不确定因素的影响,将局部土体信息合理扩展到整体地层。采用差分盒子维DBC(differential box‑counting)对地层土样颗粒分布的分形维数进行分析和计算,结合分形插值算法对局部数据高效利用。

(3)采用马尔科夫链蒙特卡洛方法(MCMC)对PFC3D的非均质地层土样进行赋值,构建在真实地层中随机取样模型,直剪实验结果表明,非均质地层土样与均质土样相比应力-应变曲线更为复杂,在大颗粒含量较少且变异性较大时,地层土样规律性不强,随大颗粒含量上升,剪切力曲线更为平滑。模拟结果能有效反映非均质地层的特定参数,基于本文建模思路,提供了一种非均质地层研究及表征方法。

参考文献

1

徐文杰胡瑞林.土石混合体概念、分类及意义[J].水文地质工程地质2009364):50-56,70. [百度学术

2

徐文杰胡瑞林岳中琦.基于数字图像处理的土石混合体细观结构[J].辽宁工程技术大学学报(自然科学版)2008271):51-53. [百度学术

3

罗荣曾亚武.岩石矿物细胞元随机性参数赋值方法研究[J].岩土力学2012337):2221-2228. [百度学术

4

吴振君.土体参数空间变异性模拟和土坡可靠度分析方法应用研究[D].北京中国科学院大学2009 [百度学术

5

Lumb P. The variability of natural soils[J]. Canadian Geotechnical Journal196632):74-97 [百度学术

6

雷进生武增琳姚奇.基于Weibull分布的非均质地层物性参数随机模型[J].地下空间与工程学报2017133):684-691. [百度学术

7

丁秀丽张宏明黄书岭.基于细观数值试验的非饱和土石混合体力学特性研究[J]. 岩石力学与工程学报2012318):1553-1566. [百度学术

8

刘新荣涂义亮王林枫.土石混合体的剪切面分形特征及强度产生机制[J].岩石力学与工程学报2017369):2260-2274. [百度学术

9

XU W JXU QHU 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 Sciences2011488):1235-1247. [百度学术

10

钟祖良涂义亮何晓勇.土石混合体物理指标及强度特性研究进展[J].地下空间与工程学报2016124):952-961 [百度学术

11

车灿辉.巨厚非均质潜水含水层抽水试验及参数计算[J].探矿工程(岩土钻掘工程)20184511):60-64 [百度学术

12

Vanmarcke E H. Probabilistic modeling of soil profiles [J]. Journal of the Geotechnical Engineering199710311):1227-1246. [百度学术

13

Elishakoff IProbabilistic methods in the theory of structure[M]. New YorkJohn Wiley & Sons, Incc1983 [百度学术

14

郭红李波张博.基于完全样本的两参数Weibull分布的参数估计[J].华中师范大学学报(自然科学版)2007413):348-351. [百度学术

15

薛亚东岳磊李硕标.含水率对土石混合体力学特性影响的试验研究[J].工程地质学报2015231):21-29. [百度学术

16

雷进生夏磊王乾峰.非均质地层地质剖面的随机预测模型[J].地下空间与工程学报2016121):84-89. [百度学术

17

唐俊伟沈平平黄平中.分形插值模拟渗透率及孔隙度平面分布(I)———理论分析[J].石油勘探与开发1997243):66-69. [百度学术

18

张颖陈晨王彧佼.基于FLAC3D对大连某工程桩基承载力数值模拟研究[J].探矿工程(岩土钻掘工程)2019465):65-71,85. [百度学术

19

LI XLIAO Q LHE 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 Sciences200441S1):702-707. [百度学术

20

廖秋林李晓郝钊.土石混合体的研究现状及研究展望[J].工程地质学报,200614(6):800-807. [百度学术

21

油新华李晓贺长俊.土石混合体实测结构模型的自动生成技术[J].岩土工程界200368):60-62. [百度学术

22

杨传俊.数字图像处理技术在土石混合体颗分中的应用[J].云南水力发电2009253):24-25,32. [百度学术

23

舒志乐刘新荣刘保县.基于分形理论的土石混合体强度特征研究[J].岩石力学与工程学报200928S1):2652-2656. [百度学术

24

石崇王盛年刘琳.基于数字图像分析的冰水堆积体结构建模与力学参数研究[J].岩土力学20123311):3393-3399. [百度学术

25

丁秀丽李耀旭王新.基于数字图像的土石混合体力学性质的颗粒流模拟[J].岩石力学与工程学报2010293):477-484. [百度学术

26

董启朋卢正詹永祥.土石混合体原位试验的颗粒流数值模拟分析[J].上海交通大学学报2013479):1382-1389. [百度学术