堆石料颗粒破碎的连续时间Markov链模型

石修松，程展林

（长江科学院水利部岩土力学与工程重点实验室，武汉 430010）

堆石料颗粒破碎的连续时间Markov链模型

石修松，程展林

（长江科学院水利部岩土力学与工程重点实验室，武汉 430010）

对各级颗粒含量进行等效分级（按粒径由小到大分为四级）并确定其破坏模式，依据颗粒破碎的Markov性建立模拟堆石料颗粒破碎的连续时间Markov链模型。根据室内三轴试验数据的统计分析确定稳态的判别模式，给出计算流程和仿真实例，得到了各级颗粒含量的动态变化曲线，并对模型的灵敏度进行了分析。结果表明：随着时间的增加，粒径较大的第三级和第四级颗粒的含量单调下降，下降的速率随时间的增加而减缓，其中第四级颗粒下降的幅度较大。第二级颗粒的含量先随时间增长，到第15 d左右后随时间下降。通过对模型的灵敏度分析，该模型对试验参数的变化是稳定的。

堆石料；颗粒破碎；破坏模式；Markov链；灵敏度分析

堆石料是由大小颗粒彼此充填而呈粒状结构的散粒体，堆石体的颗粒破碎对堆石体的强度和变形特性有着重大的影响，对高堆石坝而言，其变形特性与堆石颗粒破碎明显相关。根据一些已建的高堆石坝变形观测成果，大坝的变形较大，主要原因在于堆石颗粒在高压下重新分布以及颗粒破碎引起的孔隙压缩［1］。堆石料的应力-应变关系、强度、体变和渗透性等都与应力变化引起的颗粒的完整性或颗粒的破碎程度有关。颗粒破碎对堆石体物理力学特性的影响已经为国内外学者所重视，并开展了一系列的研究工作［1，2］，但大都以试验为手段针对颗粒破碎现象、颗粒破碎的影响因素以及颗粒破碎对堆石体强度等进行研究［3－5］，对堆石料粒径再分布的动态变化研究较少。本文通过对堆石料的级配的等效分级建立了其破坏模式，在此基础上得到连续时间的Markov链模型，根据室内三轴试验数据的统计分析确定稳态的判别模式，并通过计算机进行了仿真。

1 堆石料级配的等效分级与破坏模式

1.1 堆石料级配的等效分级

堆石料级配的等效分级是对其粒径级配的虚拟化。一级颗粒的等效粒径定义为小于某粒径的颗粒粒径的加权平均值。同理可定义更高级颗粒的等效粒径。第i级颗粒在形式上由i个一级颗粒组成。二级颗粒与其它一级颗粒簇的不同点在于，组成它的一级颗粒仅是形式上的，它本身作为一个整体参与受力和变形计算，颗粒内部是均匀的，组成它的一级颗粒在接触和运动计算中不再具有个体特性。只有在二级颗粒破碎时，组成它的一级颗粒的个体特性才恢复：即二级颗粒破碎前作为一个连续体处理不考虑空隙及一级颗粒接触面的影响。二级颗粒被模型化为单个颗粒，由此可定义更高级的颗粒。

1.2 颗粒的破坏模式

在上述模型中定义的一级颗粒根据其存在状态分为独立颗粒和非独立颗粒。所谓独立颗粒是指那些没有参与更高级别组成的一级颗粒，非独立颗粒是指那些构成了更高级颗粒一级颗粒，它们不具有作为独立颗粒时的物理力学性质。一级等效颗粒的定义具有重要的意义，它为可破碎的更高级颗粒建立了模型，高级颗粒作为一个整体参与破碎过程，它只能破碎成低级的颗粒，不具有可逆性。假设在极短的时间内高级颗粒破碎只能破碎成两部分，即一个高级颗粒破碎时只有两个低级颗粒生成，低级颗粒再参与下一次的破碎。这里以四级等效颗粒为例探讨其破碎模式，各级颗粒的破碎模式见表1，类似可得到更高级颗粒的破坏模式。可见第四级颗粒对应2种破坏模式。其他低级的颗粒破坏模式唯一。

表1 各级颗粒的破坏模式Table 1 Damage patterns of every class particles

2 连续时间Markov链模型的建立

2.1 模型的建立

Markov链是利用某变量的现在状态和动向预测该变量未来状态及其动向的一种分析手段，若X＝｛X｛t｝，t≥0｝为Markov过程［6］，它应满足下述性质：随机过程在时刻ti所处的状态x（ti）已知时，则随机过程在时刻ti＋1＞ti所处的状态x（ti＋1）仅与ti时刻所处的状态有关，与ti时刻之前发生的状态无关。显然，如果堆石料颗粒的破碎过程为脆性的话，其破碎过程符合Markov性。设各级颗粒中一级颗粒所处的状态为i，则Markov链对应的状态空间为S＝｛1，2，…，m｝，其中m为颗粒等效分级数。

设给定连续时间Markov链对应的状态转移矩阵为~P，对应的转移率矩阵［7］为

对应的分量形式为

假定在t0时刻一级颗粒到达状态i（i∈S），以τi记颗粒状态在离开i之前在i停留的时间，则τi服从指数分布［8］P｛τi＞t／X（t0）＝i｝＝exp（－λit），对应颗粒级别A－D，各级颗粒中一级颗粒所处的状态为1－4，各级颗粒含量占总颗粒的百分比ui（t），时间调整参数λi（其中状态4为吸收态，无时间调整参数），见表2。

表2 各级颗粒占总颗粒的百分比和时间调整参数Table 2 The content percent of every class particles to the total and the tim e ad justing parameters

第四级颗粒对应2种破坏模式。设其对应破坏模式1的概率分别为w。则当时间增量Δt很小时，如果一级颗粒在时刻t处于状态1，而在时刻t＋Δt变为状态2，4，此时只要求在时间（t，t＋Δt）内状态恢复。

初始条件：

由式（1）至（3），可解得pij（t），

各级颗粒中一级颗粒在占总的一级颗粒的百分比v（t0）＝［v4（t0），v3（t0），v2（t0），v1（t0）］，则

由（4）即可得v（t），从而由v（t）和u（t）的对应关系得各级颗粒含量占总颗粒的百分比u（t）。

2.2 平衡状态的判别

由式（4）可知，v（t）＝lin v（t0）P~＝0，由于破碎

t→∞的时间不可能为无限长，这显然是不合实际的。那么究竟怎样判别到达稳定状态的时刻呢，实际的情况是，由于初始的粒径级配与受力状态不适应，在外力的作用下颗粒产生破碎，最终会趋于一个稳定值。Marsal［9］在研究堆石体粒间接触应力和破坏机制的基础上，建议用破碎率BM来表征相应应力下颗粒破碎的程度，BM定义为试验前后颗粒在曲线上各粒组含量之差的正值之和。

式中，ΔWk＝Wki－Wkf，Wki为原级配曲线上某级粒组的含量，Wkf为试验后级配曲线上相同粒组的含量。上式反应了粒径破坏程度的大小，本文采用BM的定义［10］，即

式中n表示粒组数，BM的取值范围为（0～1）。

式中t0为试验的开始时刻，tf为试验终止的时刻。

考虑到颗粒破碎率的增量一般是前期较大后期逐渐减小并趋于稳定，文献［3］采用如下的指数型衰减函数来定义相对颗粒破碎率。

2.3 模型参数的获取

模型中一共有6个参数，各状态的时间调整参数λi（i＝1～3）；第四级颗粒对应破坏模式1的概率w；平衡判据的拟合参数a，b。其中w，a，b均可通过堆石料三轴试验后的粒径级配曲线的统计和拟合得到。设时刻tf对应的u（tf）＝（u4（tf），u3（tf），u2（tf），u1（tf）），通过换算得对应的v（tf）＝（v4（tf），v3（tf），v2（tf），v1（tf））。由式（4）可得：

其中：

由（6）式仅含λ1，λ2，λ3三个未知数，而且P~*为满秩矩阵，故解是唯一的。解得λi（i＝1～3）后，将其代入p4ivi（0）p4i＝v1（t）自动满足。这样转移概率矩阵中的参数均已知。

3 模拟计算

3.1 模拟步骤

首先参考相关文献［1，9］给出平衡状态的计算参数，a＝0.126 4，b＝0.106 2，σ2＝σ3＝200 kPa，σ1＝1 000 kPa，pa＝101.3 kPa，令t0＝0，t＝20 d，Δt＝1 d。u（0）＝（0.450，0.250，0.200，0.100），u（20）＝（0.145，0.129，0.242，0.484）。

计算流程如图1所示：

图1 模型的计算流程图Fig.1 The calcu lation flow chart of themodel

3.2 模拟结果

各级颗粒的初始含量百分比、试验结束时（t＝20 d）的百分比和平衡状态的百分比见表3。取w＝0.80，0.60，0.40，0.20，算得λi见表4。

表3 不同时刻对应的各级颗粒含量变化Table 3 Particle content percentage varying w ith time

表4 w＝0.80，0.60，0.40，0.20对应λi和Δt的值Table 4 Corresponding values ofλiand N to w＝0.80，0.60，0.40，0.20

由上表可知，随着第四级颗粒破坏模式为1的概率w的减少，λ1没有变化，λ2下降，λ3增加，这是显然的；因为无论第四级颗粒为那种破坏模式，颗粒均受到破坏。破坏模式为1时，有利于第三级颗粒的增加，故其时间调整参数λ2减小，λ3增加。

取w＝0.8，对应λ1＝0.033 8，λ2＝0.048 1，λ3＝0.021 4。用Matlab绘制各级颗粒含量百分比随时间的动态变化如图2（a）至图2（d）所示。

图2 颗粒含量百分比随时间的变化Fig.2 Particle content percentage varying w ith time

由上图可知，随着时间的增加，第三级和第四级颗粒的含量单调下降，下降的速率随时间的增加逐渐减缓，其中第四级颗粒含量下降的幅度较大，第三级颗粒下降幅度相对较少。第二级的含量先随时间增长，到第15 d左右后随时间下降，这是可以解释的，由于刚开始第二级颗粒的含量相对较少，它转化为一级颗粒的速度较缓，而更高级颗粒含量较高，向其转化速度较快，故其含量增加。但是随着时间的延续，更高级别颗粒的含量下降，转化速率降低，二级颗粒含量增加后向一级颗粒转化的速率增加，故其含量下降。由于一级颗粒处于吸收态，故其含量一直增加。

3.3 模型的灵敏度分析

由于a，b，w是试验参数，故其肯定存在误差，这些误差会反映在计算结果上。如果试验参数发生微小的变动，反映到模型得出的计算结果发生成倍的扩大，则该模型是不稳定的。下面就对本文模型的灵敏度进行分析。由于试验参数a，b仅对平衡判据有影响，而平衡判据是决定系统从开始状态到达稳态时间的因素。所以它们仅对持续时间有影响，这里不予分析。下面仅对w进行灵敏度分析。w取不同值时各级颗粒含量百分比随时间的变化如图3所示。

由图3（a）、图3（c）可知w对第一级和第三级颗粒含量的动态变化影响不明显，各种情况下颗粒含量的变化趋势基本相同。第二级颗粒受w的影响较大，但误差范围是可以接受的。这里有一个有趣的现象，不论w为何值，各级颗粒的动态变化曲线均交于同一点，在该点的左侧第一级和第三级颗粒含量随w的减小而减小，在右侧随w的减小而增加；第二级颗粒含量随w的变化则截然相反。

图3 w取不同值时颗粒含量百分比随时间的变化Fig.3 Particle content percentage varying w ith time by different values of w

4 结 论

本文建立了模拟堆石料颗粒破碎的连续时间Markov链模型，并给出仿真实例，得到了各级颗粒含量的动态变化曲线。研究结果表明：①随着时间的增加，第三级和第四级的含量单调下降，下降的速率随时间的增加减缓，其中第四级颗粒下降的幅度较大，第三级颗粒下降幅度相对较少。第二级的含量先随时间增长，到第15 d左右后随时间下降；②随着w的减少，系统到达稳定状态所需的时间Δt时间是上升的；③通过对模型的灵敏度分析，该模型对试验参数的变化是稳定的。不论w为何值，各级颗粒的动态变化曲线均交于同一点，在该点的两侧w的影响效应是不同的。如果围压和材料的试验参数已知，则可由该模型估计各级颗粒含量的动态变化；各级颗粒相对含量的变化将对堆石等堆石料的强度和变形特性产生影响。

该模型可在以下方面予以改进：①由于颗粒所处的状态是复杂的，可将模糊数学与连续时间Markov链模型建立相应的数学模型，这样与实际也会更加吻合；②本文仅是在理论上对颗粒破碎动态模拟的一次尝试，尚待缺乏有力的试验验证，故可通过试验对该模型进行验证；③作为该模型的应用，可在本构模型中考虑颗粒动态破碎的影响。

［1］ 刘汉龙，秦红玉，高玉峰，等.堆石堆石料颗粒破碎试验研究［J］.岩土力学，2005，26（4）：562－566.（LIU Han-long，QIN Hong-yu，GAO Yu-feng，et al.Experimental study on particle breakage of rockfill and coarse aggregates［J］.Rock and Soil Mechanics，2005，26（4）：562－566.（in Chinese））

［2］ 郭熙灵，胡 辉，包承纲.堆石料颗粒破碎对剪胀性及抗剪强度的影响［J］.岩土工程学报，1997，19（3）：83－88.（GUO Xi-ling，HU Hui，BAO Cheng-gang.Experimental studies of the effects of grain breakage on the dilatancy and shear strength of rock fill［J］.Chinese Journal of Geotechnical Engineering，1997，19（3）：83－88.（in Chinese））

［3］ 迟世春，贾宇峰.土颗粒破碎对罗维剪胀模型模型的修正［J］.岩土工程学报，2005，27（11）：1266－1269.（CHI Shi-chun，JIA Yu-feng.Rowe’s stress-dilatancy modelmodified for energy dissipation of particle breakage［J］.Chinese Journal of Geotechnical Engineering，2005，27（11）：1266－1269.（in Chinese））

［4］ 楚锡华，李锡夔.离散颗粒多尺度分级模型与破碎模拟［J］.大连理工大学学报，2006，46（3）：319－326.（CHU Xi-hua，LIXi-xie.Hierarccalmulti-scale discrete particle model and crushing simulation［J］.Journal of Dalian University of Technology，2006，46（3）：319－326.（in Chinese））

［5］ 张 柱，杨云川，晋艳娟.单颗粒破碎机理分析［J］.太原科技大学学报，2005，26（4）：306－308.（ZHANG Zhu，YANG Yun-chuan，JIN Yan-juan.Persistent currents inmesoscopic ring driven by time-dependentelectric fields［J］.Journal of Taiyuan University of Science and Technology，2005，26（4）：306－308.（in Chinese））

［6］ 张 波，张景肖.应用随机过程［M］.北京：清华大学出版社，2004：115－125.（ZHANG Bo，ZHANG Jing-xiao.Applied Stochastic Processes［M］.Beijing：Tsinghua University Press，2004：115－125.（in Chinese））

［7］ CARLE SF，FOGGG E.Modeling spatial variability with one andmulti-dimensional continuousmarkov chains［J］.Mathematical Geology，1997，29（7）：891－917.

［8］ 林元烈.应用随机过程［M］.北京：清华大学出版社，2002：223－235.（LIN Yuan-lie.Applied Stochastic Processes［M］.Beijing：Tsinghua University Press，2002：223－235.（in Chinese））

［9］ MARSAL R J，Mechanical Properties of Rockfill［M］.New York：John Wiley＆Sons，1973：109－200.

［10］米占宽，李国英，陈铁林.考虑颗粒破碎的堆石体本构模型［J］.岩土工程学报，2007，29（12）：1865－1869.（MI Zhan-kuan，LIGuo-ying，CHEN Tie-1in.Constitutive model for rockfill material considering grain crushing［J］.Chinese Journal of Geotechnical Engineering，2007，29（12）：1865－1869.（in Chinese））

（编辑：刘运飞）

Continual M arkov Chain M odel for Particle Breakage of Rock fill

SHIXiu-song，CHENG Zhan-lin
（Key Laboratory of Geotechnical Mechanics and Engineering of The Ministry ofWater Resources，Yangtze River Scientific Research Institute，Wuhan 430010，China）

Carrying on the equivalent graduation to all levels of particle content and determining its destruction mode，the author constructs a continual Markov chainmodel based on the Markov nature of the particle breakage to simulate the particle breakage of coarse aggregates.According to data of the indoor triaxial tests，the criterion of the steady-state is obtained and then a simulation example is given.The results show that，as the time increases，the particle contents of the fourth and third classes aremonotonically decreased，and their rate slowed down with time，of which the fourth class is dropped largely，and that the particle content of the second class is first grown with the time and then declined after 15 days.Through the sensitivity analysis of themodel，it is stablewith the test parameters.

particle breakage；damage pattern；Markov chain；sensitivity analysis

TU41

A

1001－5485（2010）06－0038－05

2009-09-03；

2010-01-13

石修松（1985-），男，安徽淮北人，硕士研究生，主要从事堆石料的力学特性与本构关系的研究，（电话）15971436907（电子信箱）qingsongsaint＠126.com。