离散元软件

2024-07-31

离散元软件（通用7篇）

离散元软件篇1

1 工程背景

随着工程规模的日益扩大, 作为建筑材料以及地基介质的岩土体一直以来都是人们研究的热点。目前, 解决岩土类工程问题的主要研究手段有三大类, 分别是理论分析、室内实验和数值模拟。其中, 数值模拟又因其便捷、高效以及可控性成为当今最为常用的研究方法。常用的数值模拟方法主要基于有限元分析理论或离散元分析理论, 运用计算机软件对实际的工程问题进行模拟, 从而发现问题、得到内在规律并作出研究结论, 对于指导工程实践具有重大意义。

2 研究方法

针对不连续介质问题以及大变形问题[1], 近年来新兴的离散元理论逐步代替传统的有限元理论被人们采用。其主要理论依据是牛顿第二定律和力与位移的关系, 求解思路为:将求解空间离散为分离的单元阵, 相邻两个单元之间用合理的元件相连接;基本变量是单元间的相对位移, 通过位移和力的关系可以分别得到单元之间的法向、切向作用力;求单元上受到的合力及合力矩, 由牛顿运动定律计算出单元体的加速度;最后, 积分得到单元的速度和位移。通过以上计算过程可以得到任意时刻所有单元体的位移、速度和加速度。

在基于离散元分析理论开发出来的计算机软件中, 最为广泛使用的是美国ITASCA公司开发的UDEC (universal distinct element code) 、3DEC (3-dimensional distinct element code) 块体离散元程序, 和PFC2D (particle flow code in 2 dimensions) 、PFC3D (particle flow code in 3 dimensions) 软件。前者用于模拟岩石的力学过程, 后者则可以模拟颗粒流问题, 如粗粒土、砂土的力学现象等。现阶段, 国内使用PFC软件进行离散元数值模拟研究的情况并不多见, 目前国内已出版的相关中文书籍也仅有《离散元颗粒流软件 (PFC) 在道路工程中的应用》和软件自带的全英文版帮助手册。

PFC软件[2]从微观结构角度着手, 将土体看成是由土颗粒组成, 其宏观力学特性取决于颗粒与颗粒之间粘结 (接触) 方式的几何、物理特性。在PFC中用等厚度的刚性圆盘代表颗粒, 颗粒之间允许重叠, 同样遵循牛顿第二定律。颗粒的破坏主要有剪切破坏和张拉破坏两种方式, 当颗粒发生破碎时土体的宏观力学特性也会发生相应变化, 即介质内颗粒接触状态的变化决定了介质的本构关系。因此, 在PFC计算中无需给定材料的宏观本构关系和对应的参数[3], 这些传统的参数和力学特性在程序中可以自动得到, 需要使用者自行定义的是颗粒的几何力学参数, 包括颗粒级配、刚度、摩擦力、粘结介质强度等细观力学参数。

3 参数标定方法及常见问题

需要注意的一点是, 运用PFC软件进离散元数值模拟时, 可以自行定义颗粒的物理力学参数, 也可以根据实际情况标定出较为真实的颗粒参数[4]。而通常情况下, 为了使模拟得到的结果更具有说服力, 研究者通常会采用后一种方法获取颗粒基本参数。因此, 数值模拟的首要任务是进行参数标定, 而在标定过程中也会遇到一些共性的问题, 下面以笔者的实际模拟过程为例进行简要说明。笔者的主要研究内容是运用PFC3.0软件进行砂土的颗粒破碎模拟, 选用真实的室内三轴剪切试验, 在PFC中对其进行三轴试验的数值模拟, 将模拟得到的宏观应力应变曲线与真实的室内实验得到的应力应变曲线进行比对, 如果二者在线性、斜率和峰值强度上都能够很好的定量吻合, 则说明数值模拟中采用的颗粒参数是相对真实可靠的。在PFC中如何根据宏观参数确定细观参数是一个困难的问题。在标定阶段, 影响宏观应力应变曲线的因素主要有颗粒之间的摩擦系数、颗粒与墙体间的摩擦系数、颗粒的切向及法向刚度、试样的初始孔隙率和剪切速率。如何选取适当的参数成为数值模拟的一个关键性问题, 参数标定的具体方案是:先拟定好以上参数, 再通过控制变量法逐一调整每个因素, 直至最后的曲线与实际曲线定量地对应上。值得注意的一点是, 先要保证参数的数量级是正确的, 在此基础上再微调参数的具体数值。

4 研究结论

通过笔者近一个月的标定过程, 最后总结出一些关于三轴试验中参数标定的规律:

(1) 在三轴试验条件下, 容易发生颗粒破碎的砂土的摩擦系数在0.7~0.9, 二维模拟情况下选取的摩擦系数尽量控制在0.9以下, 否则模拟将失真。

(2) 切向刚度、法向刚度的数量级在e9~e10, 也可以调整切向、法向的刚度比;

(3) 二维情况下的初始孔隙率在0.2~0.35之间, 曲线的线性则主要取决于初始孔隙率;

(4) 应力应变曲线的峰值强度很大程度上由加载速率和摩擦系数决定, 刚度影响并不是很大。

对于密砂和中密砂, 应力应变曲线在加载初期陡升, 后期处于平缓状态, 会出现应力峰值;而松砂的应力应变曲线则一直处于递增状态, 是否会出现峰值由围压和砂的种类决定。

摘要：简要介绍了运用离散元理论及颗粒流软件PFC进行岩土类工程问题数值模拟的研究思路和方法, 简述了运用PFC软件进行数值模拟时通常会遇到的关键性问题, 通过结合笔者在研究过程中的实践经验, 针对数值模拟的首要步骤即参数标定这一步, 给出了相应的解决方案及研究结论。

关键词：颗粒流软件PFC,离散元,数值模拟,参数标定

参考文献

[1]王卫华, 李夕兵.离散元法及其在岩土工程中的应用综述[J].岩土工程技术, 2005 (19) :177-181.

[2]陈俊, 张东, 黄晓明.离散元颗粒流软件 (PFC) 在道路工程中的应用[J].北京:人民交通出版社, 2015, 1ISBN 978-7-114-11667-4.

[3]赵吉坤, 李骅, 张慧清.基于离散元法的岩土细观破坏及参数影响研究[J].防灾减灾工程学报, 2013 (33) :218-224.

[4]杨洋, 唐寿高.颗粒流的离散元模拟及其进展[J].中国粉末技术, 2006 (05) :38-43.

[5]PFC3.0 help.美国, ITASCA.

离散元软件篇2

单轴贯人试验就是一种满足这种要求的试验方法。一般来说, 单轴贯入试验有两种特点, 在试件内部, 剪应力的分布与在实际路面的车轮荷载下的剪应力分布基本一致。由于压头直径大大小于试件的直径, 在加载过程中周围的材料将形成对压头下圆柱体的侧向约束, 侧向约束的大小与沥青的性能和混合料的性能密切相关。在这个过程中, 竖向压力是主动的, 侧向压力是被动的, 竖向压力不同侧向约束也不同, 这个过程与实际路面相同。

试件直径×高度=150×100 mm, 底端竖直方向固定, 试件四周受侧向限制作用, 经计算, 得其中心点处各应力沿试件高度方向和径向的应力路径, 之所以采用压头直径为28.5 mm, 是因为在该尺寸下, 试件内部的应力值和道路上在实测轮载作用下的应力值最为接近, 且较常用面层混合料的集料尺寸大1.5倍以上。

目前, 关于沥青混合料的单轴贯入试验还停留在传统试验阶段, 这就使得一个项目或者课题中, 完成沥青混合料单轴贯入试验需要大量的重复劳动, 而且在试验中由于种种原因有时会导致试验数据的分布离散, 不能准确反应试件的性能和材料真实的结果。在这种情况下, 科研人员研究并运用PFC3D软件建立试验模型, 模拟传统单轴贯入实验, 并得到相近的试验数据, 并将计算机模拟的结果与实际实验结果对比, 确保实验结果真实可信。

二、PFC3D软件简介

PFC (Particle Flow Code) 是利用显式差分算法和离散元理论开发的微、细观力学程序, 它是从介质的基本粒子结构的角度考虑介质的基本力学特性, 并认为给定介质在不同应力条件下的基本特性主要取决于粒子之间接触状态的变化, 适用研究粒状集合体的破裂和破裂发展问题, 以及颗粒的流动 (大位移) 问题。

三、模拟条件

本文将应用PFC 3D离散元计算软件来模拟沥青混合料的单轴贯入试验, 试验条件如下。

(1) 沥青混合料类型:AC16;

(2) 试件尺寸:按照标准马歇尔试件尺寸, 直径101.6 mm, 高度63.5 mm;

(3) 压头形式:圆形;

(4) 压头尺寸:30 mm;

(5) 侧限:有;

(6) 材料参数:颗粒间刚度kn=ks=1e8;强度nstren=sstren=500;摩擦系数fric=0.8。

考虑到沥青路面的上面层、中面层和下面层是分开铺装, 虽然在铺装过程中层间设有联接层, 但不是一次成型, 路面使用过程中会存在一定的缺陷。64 mm与上面层和中面层的厚度相似, 符合路面实际受力情况, 故选用该尺寸的标准马歇尔试件。除此之外, 沥青混合料中的沥青、矿粉以及纤维材料等作为沥青胶浆为骨架颗粒提供粘结作用, 本次数值模拟采用接触粘结模型以及滑动模型。接触粘结模型认为两个球体颗粒之间的接触为点接触, 可把接触力分为带有法向和切向两个方向的刚度。这种模型只能传递力而不能传递力矩。接触粘结和滑动粘结不能同时存在。

四、建模过程

在采用PFC3D建模过程中, 我们假设颗粒全部为球形刚性体, 接触只发生在很小的范围内, 即点接触。接触为柔性接触, 即在接触处允许有一定的重叠, 重叠量的大小由接触力的大小决定, 与颗粒尺寸相比可以忽略。调用PFC3D中的WALL命令生成一个半径为R0, 高度为H的圆柱形墙体, 其上下底面圆心坐标分别为 (X0, Y0, Z0+H) 和 (X0, Y0, Z0) , 以此达到对颗粒集合体的轴向或者周向压缩。墙体之间不会产生接触力。因此在PFC3D中存在两种接触类型:颗粒—颗粒接触和颗粒—墙体接触。在试验中要求圆柱体墙面单元的刚度足够大, 不易产生变形。设每档集料的半径范围Rmin

(1) 首先生成一个球形颗粒, 然后检查该球体是否完全位于该圆柱体范围之内。如果球颗粒满足上述条件, 则继续生成下一个颗粒;若不满足, 则将该颗粒删除。

(2) 建立一个数据链并随时更新。若发生嵌入现象, 则将该颗粒及其相关数据一并删除;重复之前的步骤生成新的颗粒。若生成球体颗粒的尝试次数达到一定的数量, 则认为该空间内已无法再加入新的颗粒, 随即保存相关数据。

此时生成的球形颗粒排列较为松散, 因此, 我们采用压缩成型法来提高颗粒排列的密度。现在模型顶部墙体上施加一定的速度, 使得墙体自上而下运动, 在墙体移动的过程中会推动一些颗粒向下移动。在压缩过程中, 墙体的速度不宜过大, 否则会导致结果与实际情况不符。同时还需要实时监测球体颗粒对上底墙和下底墙的作用力Ft和Fd的比值, 理论上当Ft/Fd=1时可认为颗粒被压缩密实, 但此时会模型内部会产生较大的内应力。之后撤除该速度, 使球体颗粒能充分调整相互之间的位置, 以消除颗粒间的不平衡力, 当不平衡力趋于平衡时模型压缩成型完成。最终得到底面直径为101.6 mm, 高为64 mm的标准马歇尔试件模型。

五、数据分析

模型建好之后, 采用直径为60 mm的圆柱形压头对马歇尔试件加压。利用PFC3D-crack调用实验模拟程序, 得到马歇尔试件的单轴贯入试验的轴向应力—应变曲线。

由此可知, AC-16沥青混凝土在直径60 mm的压头作用下进行单轴贯入试验, AC-16沥青混凝土在60mm直径压头作用下的单轴贯入试验, 试件所能承受的最大轴向应力为:

最大压力为:

试件轴向应变为:

试件的轴向位移为:

由于是在有侧限的条件下进行的试验, 因此试件的径向位移为0。

六、结语

本次实验在标准马歇尔试件上对有无侧线, 不同压头尺寸的情况下, 试件破坏的情况进行了详细的建模和分析。利用在PFC 3D软件下模拟的结果, 可以对实验结果进行对比分析, 得到沥青混合料的断裂过程和剪切强度的变化, 从而分析结构的破坏过程, 得到提高沥青混合料强度的有效方法, 并在实际工程中使用, 检验试验结果的准确性, 具有很好的科研价值和经济效益。

参考文献

[1]刘玉.沥青混合料习惯结构模型的离散元数值分析[D].郑州:河南工业大学, 2005.

[2]鲍飞剑.沥青混合料抗剪强度影响因素的试验研究[J].科技信息, 2011.

离散元软件篇3

长期以来常规球磨机是细磨的主要设备,由于它主要采用冲击粉磨方式进行磨矿,能量虚耗在转动笨重的筒体上,利用率低。搅拌球磨机是一种有发展前途而且能量利用率高的超细粉碎设备,其特点是采用搅拌方法直接输入能量。球磨机介质颗粒运动学分析首先是由E.W.Davis于1919年提出创立的,他建立了著名的球磨机介质运动学分析理论,即Davis理论[1]。我国的周思浦[2]也做了大量工作,在实验球磨机中观察了介质运动状态。在搅拌磨机中,介质的运动学分析同样重要,国内外学者进行了一些研究工作,但是介质球的运动学分析及能量转换关系较为复杂,研究者较少[3]。本文引入离散元方法[4],以真实的物理参数作为初始条件对搅拌磨系统进行模拟分析。

搅拌磨机主要应用于超细粉磨作业,在实际运行过程中,内部的介质由钢球颗粒和物料组成,但物料微粒的粒径很小,大部分分布在几十到几百微米之间,钢球直径为5mm,此外,钢球介质的填充率在80%以上,在如此高的填充率下,钢球之间的间隙非常小。根据以上分析,物料对仿真计算影响较小,为了简化模型,我们以钢球介质为模拟对象。

1 搅拌磨机系统数学模型

1.1 物理模型

图1为搅拌磨机的三维图[5],在粉磨过程中,电机经过减速器带动螺旋搅拌器,筒体静止,搅拌器在填充有钢球和物料的磨筒内做旋转运动,钢球和物料在搅拌器的驱动下做复杂的自转和循环运动,物料在冲击﹑摩擦和剪切综合作用下破碎[6]。

根据计算机仿真步骤,首先要选定所要模拟的搅拌磨机类型。为了更好地模拟工作情况,本文以某公司试验用搅拌磨机为研究对象,并进行实体比例尺寸建模,磨机的几何参数见表1。

在建立搅拌磨离散元模型前,除了要规定搅拌磨的规格,还要选好合适的模型输入参数,即需要设置介质参数来进行仿真计算,钢球的关键参数信息见表2[7]。

搅拌磨系统的建模过程如下:

(1)生成筒体和螺旋搅拌器结构。实际磨机筒体为一圆柱体,下端有端盖封闭,螺旋搅拌器为实体结构。在模拟中,筒体用圆柱墙体表示,端盖用平面墙体表示,搅拌轴用圆柱墙体表示,叶片用螺旋状墙体表示,生成的模型见图2。

(2)建立钢球模型。在PFC3D中,钢球模型表示为颗粒集合,颗粒集合充满在筒壁和搅拌器之间的特定区域里,编程中要用到表2的相关信息。

(3)自然堆积过程。在整个筒体内,各钢球之间没有接触,还处于非平衡状态。对钢球施加重力,钢球在重力作用下以自由落体运动下落,呈自然堆积状态,最终达到初始平衡状态,见图3。

1.2 接触力学模型

模型的建立是基于单个介质是刚体的假设基础之上的,刚体在接触点允许重叠,重叠量的大小与接触力的大小相关。两颗粒单元的接触作用通常分解为法向和切向来讨论,颗粒受力如图4所示。

根据Hertz Mindlin接触模型[8],法向力计算公式为

式中,E*为等效弹性模量,Pa;δn为法向重叠量,m;R*为等效半径,m。

R*和E*计算公式如下:

式中,R1和R2分别为两接触颗粒的半径,在本模型中,两颗粒半径相同,都是R;E1和μ1﹑E2和μ2分别为两颗粒的弹性模量和泊松比,在本模型中,两颗粒的弹性模量和泊松比相同,都是E和μ。

1.3 计算模型

介质系统内颗粒间的相互作用由接触力学原理得出,在每一时步扫描到的颗粒瞬时运动由该颗粒所受的接触力确定,边界条件由空间条件(螺旋搅拌器或圆筒壁)限定。接触模型是拟静态的,但整个系统是动态的,在每一时步对所有的单元进行扫描迭代计算,确定各个单元的接触状态,从而得到整个颗粒系统的运动形态。

离散元程序在计算时,对接触力的计算需要寻找颗粒间的接触点。Cundall提出了时间差分格式的计算方法来循环计算接触力增量和随后颗粒的运动[9]。对每一次循环,每一个颗粒的位移和加速度由牛顿第二运动定律得出,通过一个很小的时步下的数值积分得出更新后的速度和位移。每个颗粒的速度用于寻找接触颗粒之间的叠合量,以便根据力与位移叠合量之间的关系来计算法向接触力增量。解出接触力后便可以得到各个颗粒上的不平衡力,由此根据下个时步来计算各个颗粒的加速度。整个颗粒系统的动态循环过程示意图见图5。

对于颗粒在一个时步Δt内,外力和加速度都假定为常数,因此由牛顿第二定律,可得颗粒的运动方程:

式中,i=1,2,3,表示沿x,y,z三个坐标轴方向的分量;Fi为球颗粒的不平衡力分量,N;β为黏滞阻尼系数,N/(m/s);vi为平动速度分量,m/s;m为球颗粒的质量,kg;Mi为由接触力引起的不平衡力矩分量,N·m;βg为旋转(黏性)阻尼系数N·m·s/rad;ωi为转动速度分量,rad/s;I为球的转动惯量,kg·m2。

用差分法表示t时刻的平衡方程:

整理式(6)、式(7)可得速度迭代的表达式:

得到速度的表达式后,可以根据下式计算相应的每个颗粒新的位置:

式中,ui和φi分别为颗粒的线位移和角位移。

此循环不断继续下去,直到力或位移达到平衡状态为止,或是循环达到一定次数为止。

2 仿真计算及结果分析

搅拌磨系统的离散元模型生成后,赋予螺旋搅拌器240r/min的转速。程序运行结束后,我们得到了一系列的结果并对结果进行了分析。

2.1 介质运动形态

速度矢量图能够直观地观察介质的运动形态,从图6中可以看出,螺旋搅拌器对速度矢量影响较大,决定了速度矢量的分布,紧贴螺旋壁面的钢球速度矢量为最大值,随着钢球位置远离壁面,速度值逐渐减小,上层介质球和紧贴圆柱筒体壁面的速度值较小。可以得出,在搅拌器周围,钢球和物料运动充分,粉磨效果好,远离搅拌器壁面处粉磨效果较差。

2.2 系统平衡状态判定

平均接触应力代表了介质与介质之间﹑介质与筒体和搅拌器之间接触应力的平均值,平均不平衡应力是评估模型状态的一个重要参数,它的值需要与平均接触应力做对比。如图7所示,平均接触应力达到峰值后基本稳定在一定值,并上下波动,平均不平衡力近似恒定,且值较小,可以判定系统呈近似平衡状态。由此可见,搅拌磨机在粉磨过程中近似处于平衡状态,这便于在特定的运行时间内,得到理想中的物料细度。

2.3 系统能量耗散

搅拌磨机系统的能量包括几个部分,在这里只分析摩擦功和动能。摩擦分为钢球之间、钢球与螺旋搅拌器或筒体之间的摩擦,摩擦均为滑动摩擦,各个接触处的滑动摩擦力均做功。系统的动能是总体钢球的动能,包括平移时的动能和转动时的动能。功能曲线图见图8和图9。

从图8摩擦功曲线中可以看出,摩擦功随着时步呈近似线性增加;图9动能曲线表明,介质的动能在搅拌器刚开始转动时急剧增加,增加至某一平衡值后上下波动,且随着时步稳定在一定的数值范围内。随着时步的增加,总动能围绕某一定值上下波动,摩擦功近似线性增加。钢球的总动能越大,钢球之间的运动越剧烈,冲击力越大,搅拌磨机的粉磨效率越高。

3 结论

(1)在全部钢球中,螺旋面附近的钢球速度值较大,在这个区域里物料的粉磨效果好,同时也加剧了钢球的磨损,这对于研究合适的钢球速度有贡献意义。

(2)系统的平均不平衡应力较小,系统在刚开始时就能达到近似平衡状态,且以后运行一直保持近似平衡状态。

(3)搅拌磨系统的摩擦功随着时步线性增加,但钢球的总动能在较短时步内达到峰值后在平均值上下波动,表明介质球的复杂运动影响动能,但对摩擦功的线性增加基本没有影响。

摘要：针对搅拌磨机应用广泛但对其内部介质运动信息了解甚少的现状,以某公司立式搅拌磨机为研究对象,运用三维离散元方法跟踪磨机内的介质颗粒,考虑其重力、摩擦力和碰撞力的综合作用,建立了搅拌磨机系统的数学模型,对磨机内介质颗粒的复杂运动全过程进行了数值模拟。研究结果表明,在螺旋搅拌器附近的介质球速度值较大,粉磨效率较高,搅拌磨机系统在运行中达到近似平衡状态,系统摩擦功随时间线性增加,但总动能在固定值上下波动。

关键词：搅拌磨机,离散元法,介质运动,数值模拟

参考文献

[1]陈炳辰.磨矿原理[M].北京:冶金工业出版社,1989.

[2]周恩浦.粉碎机械的理论与应用[M].长沙:中南工业大学出版社,2004.

[3]张国旺.超细搅拌磨机的流场模拟和应用研究[D].长沙:中南大学,2005.

[4]王泳嘉,刑纪波.离散单元法及其在岩土中的应用[M].沈阳:东北工学院出版社,1991.

[5]Jankovic A,Sinclair S.The Shape of Product SizeDistributions in Stirred Mills[J].Minerals Engi-neering,2006,19:1528-1536.

[6]何亮,董为民.超临速球磨机介质运动三维离散元数值模拟[D].昆明:昆明理工大学,2007.

[7]史国军.基于三维离散单元法的球磨机介质工作参数研究[D].昆明:昆明理工大学,2008.

[8]朴香兰,王国强,张英爽.转送站的离散元模拟研究[J].矿山机械,2009,37(23):57-61.

离散元胞蚂蚁算法及其收敛性篇4

组合优化问题在现实领域中有着许多广泛的应用, 因而寻找其实际而有效的算法就显得颇为重要。遗传算法、模拟退火算法、蚂蚁算法等一系列来自自然界的演化型算法, 在一系列困难的组合优化问题 (TSP, QAP, JSP等) 求解中已取得了成效。

蚂蚁算法是由意大利学者Colorni等人[1]于1991年通过模拟自然界中蚂蚁的寻径行为提出的一种基于种群的启发式仿生进化算法, 许多学者从不同的角度对其进行了研究[2,3,4,5,6,7,8,9], 并将其应用于许多实际领域[10,11,12,13] 。

元胞自动机是冯诺伊曼最早提出, 沃尔夫勒姆等人使其成为一种可模拟复杂结构和过程的模型[14], 并已应用于许多实际领域[15,16,17,18,19]。离散元胞蚂蚁算法是一种新型的仿生算法, 它利用元胞在离散元胞空间的演化规律和蚂蚁寻优的特点, 为解决实际问题提供了一种优化方法。

以下首先给出离散元胞蚂蚁算法, 然后对算法的收敛性进行了论证, 最后给出一些数值仿真。

1 元胞自动机原理

冯诺伊曼最初提出元胞自动机是为研究复杂系统的行为提供一种构架。元胞自动机是一个时间、空间离散的动力系统。元胞自动机由规则网格中的元胞组成, 每个元胞的状态是k个有限状态之一。元胞自动机中的元胞状态依据一个已定义的局部规则, 即按照相邻元胞和元胞本身的前一时刻的值同步更新。

从集合论角度元胞自动机有着严格地描述和定义[20]。用数学符号来表示, 标准的元胞自动机是一个四元组:

A= (Ld, S, N, f ) 。

这里, A代表一个元胞自动机系统;L表示元胞空间, d是一正整数, 表示元胞自动机内元胞空间的维数; S是元胞的有限、离散状态集合; N表示一个邻域内所有元胞的组合 (包括中心元胞) ; f表示将SN映射到S上的一个局部转换函数.可以记为:

f:S $_{t}^{Ν}$ → St+1。

对元胞空间内的所有元胞, 独立施加上述局部函数, 则可得到全局的演化。

2 离散元胞蚂蚁算法原理

元胞是元胞自动机的最基本的组成部分。元胞空间是元胞的集合。以TSP问题为例, 阐述离散元胞蚂蚁算法的数学描述。

G= (V, E) 为赋权图, V 顶点集, E为边集, arc (i, j) 连结i和j节点, 并且有赋值dij = d (i, j) 。

经典的TSP问题可写为如下的数学规划模型

min Z= $\sum_{i = 1}^{n} \sum_{j = 1}^{n} d_{i j} x_{i j}$

$s . t {\begin{cases} \sum_{j = 1}^{n} x_{i j} = 1, i \in V \\ \sum_{i = 1}^{n} x_{i j} = 1, j \in V \\ \sum_{i \in S} \sum_{j \in S} x_{i j} \leq | S 1 | - 1, \forall S 1 \subset V \\ x_{i j} \in {0, 1} \end{cases}$

定义1 给定城市元素的集合C = { c1, c2, …, cn }, 则C中的任意排序组合的集合为元胞空间, 可表为:

L={CellX= (c1, …, ci, …, cj, …, cn) |, ci∈C, ci≠cj, i, j=1, 2, …, n}。

其中每一个组合CellX为元胞。

定义2 元胞邻居采用扩展Moore邻居类型

NMoore={CellY|diff (CellY-CellX) ≤r;CellX, CellY∈L}其中diff (CellY-CellX) ≤r为两个组合排序的差异, 若无差异为0, 有差异时, 最小为2。r为差异的程度, 现r为2。

定义3 蚂蚁的相邻结点的转移概率定义为

$Ρ_{i j} = \frac{[τ_{i j}]^{α} [η_{i j}]^{β}}{\sum_{k} [τ_{i k}]^{α} [η_{i k}]^{β}}$ 。

其中, ηij为边弧 (i, j) 的能见度 (visibility) , 即1/dij, τij为边弧 (i, j) 的轨迹强度 (intensity) , α为轨迹的相对重要性 (α≥0) , β为能见度的相对重要性 (β≥0) 。

定义4 信息素强度更新方程为:τ $_{i j}^{n e w}$ =ρτ $_{i j}^{o l d}$ +∑k △τ $_{i j}^{k}$ , 其中, △τ $_{i j}^{k}$ 为蚂蚁k于边弧 (i, j) 上留下的单位长度轨迹信息素数量,

$Δ τ_{i j}^{k} = {\begin{cases} Q / d_{i j}, 若 (i, j) 在最优路径上； \\ 0 ‚ 其他。 \end{cases}$

Q为体现蚂蚁所留轨迹数量的一个常数, ρ体现强度的持久性 (0≤ρ<1) 。

定义5 元胞演化规则, 依据元胞邻居的定义计算其邻居的目标解, 比较元胞和其邻居的差异, 选择最好的目标解。

于是, 元胞蚂蚁算法主要步骤可用以下流程图描述叙述如下:

Step 1 nc←0; (nc为迭代步数或搜索次数) 各τij和△τij的初始化, 将m个蚂蚁置于n个顶点上;

Step 2 将各蚂蚁的初始出发点置于当前解集中;对每个蚂蚁k (k=1, …, m) , 按概率P $_{i j}^{k}$ 移至相邻结点cj, 将顶点j置于当前解集;

Step 3 计算各蚂蚁的目标函数值Zk (k=1, …, m) , 记录当前的最好解;

Step 4 按元胞邻居的定义, 在邻居范围内演化, 记录最好解;

Step 5 按更新方程修改轨迹强度;

Step 6 对各边弧 (i, j) , 置△τij←0; nc←nc+1;

Step 7 若nc<预定的迭代次数且无退化行为 (即找到的都是相同解) , 则转步骤2;

Step 8 输出目前的最好解。

对于元胞蚂蚁算法实现的描述, 可知, 当选择差异的程度为2时, 邻居内演化其运算时间的复杂度为O (n2nc) , 因此元胞蚂蚁算法整个过程的时间复杂度为O (ncn2m) , 其中n为城市节点数目, m为蚂蚁数, nc为运算迭代次数。如果选取m≈n, 则算法的时间复杂度为O (ncn3) 。应用上, 这个复杂度在计算时间上是可接受的。

3 算法收敛性分析

设 (Ω, A, p) 表示一个完全的概率测度空间, 这里p (Ω) =1, a = {a1, a2, …ai, …}, 为算法每次迭代输出的最优解序列, 其中ai= (c1, …, ci, …cn) , ci∈C, ai∈L。L为元胞空间, ai为一个元胞, C为元胞的状态空间。

定义6 转移算子Tp是指按概率蚂蚁从结点i到结点j的转移过程, 它是一种元胞的状态空间到状态空间的映射Tp :Ω × C → C, 定义转移算子Tp为

p (ω: $Τ_{p} (ω, (c_{i}, c_{j})) = p_{i j}^{} = \frac{[τ_{i j}]^{α} [η_{i j}]^{β}}{\sum_{k} [τ_{i k}]^{α} [η_{i k}]^{β}}$ 。

定义7 演化算子Tf是指元胞区域的演化规律, 它是元胞空间到元胞空间的映射Tf:Ω×L→L, 可定义为

p (ω:Tf (ω, Cell

$X)) = {\begin{cases} = C e l l X ‚ Ζ_{C e l l X} \geq Ζ_{C e l l Y} ； \\ = C e l l Y ‚ Ζ_{C e l l X} < Ζ_{C e l l Y} 。 \end{cases}$

由离散元胞蚂蚁算法的求解过程可知, 它是一个迭代过程, 在每1次迭代有若干个转移和演化算子, 因此可以进一步抽象一为个映射T, 即

$Τ = \prod^{n} Τ_{p} \prod^{n} Τ_{f}$ ;T:Ω×a→a。

离散元胞蚂蚁算法在求解TSP问题过程中, 每次都保留最好解, 整个求解过程中目标函数是一个非增序列。由于优化过程只关心满意解的变化过程, 为了分析方便, 每次迭代用最优解来代替该次迭代的目标值。于是, 可定义影射为:

an+1=T (ω) an;an, an+1∈a。

其中an为第n次迭代输出的最优解序列, an+1为第n+1次迭代输出的最优解序列。则对应的Zn和Zn+1为输出的最优目标解。

定义8 度量d:a×a→R, 其中d的表达式定义如下

$d (a_{i}, a_{j}) = {\begin{cases} | Ζ_{i} - Μ | + | Ζ_{j} - Μ |, a_{i} \neq a_{j} ； \\ 0, a_{i} = a_{j} 。 \end{cases}$

其中M是Z的下界, 对于一个最小值问题, M必存在。

定理1 (a, d) 是完备可分的度量空间。

证明a为非空集合, d为a×a上的实值函数, 对∀ai, aj∈a, d (ai, aj) 满足:

1) d (ai, aj) >0, 当且仅当ai=aj, d (ai, aj) =0, 满足非负性;

2) d (ai, aj) =d (aj, ai) , 满足对称性;

3) 满足三角不等式

d (ai, aj) =|Zi-M|+|Zj-M|≤|Zi-M|+|Zk-M|+|Zk-M|+|Zj-M|=d (ai, ak) +d (ak, aj) 。

所以 (a, d) 是度量空间。

a是一有限状态空间, 即a中解集的数目是有限的。对任意柯西列ai, 及任意ε>0, 存在自然数N, 当自然数n, k>N时, d (an, ak) <ε, 当n→∞时, an→a*, 且a*∈a。因此, (a, d) 是完备的度量空间。

设G是a的子集, 由于a为有限集合, 因此G为可数子集。又G的闭包包含a中所有元素, 所以G在a中稠密, 这就证明了 (a, d) 是可分的。

因此 (a, d) 是完备可分的度量空间。

定理2 元胞蚂蚁算法形成的影射T为的随机压缩算子。

证明如果存在非负实值随机变量

k (ω) <1 a.s. ∀x, y∈a。

使得 p ({ω:d (T (ω) x, T (ω) y) ≤k (ω) d (x, y) }) =1。

随机算子称为随机压缩算子[21,22]。

根据算法原理, 每次都保留最好解, 整个求解过程中目标函数是一个非增序列。则有:

Zi-1≥Zi≥Zi+1;

d (T (ω, ai-1) , T (ω, ai) ) =d (ai, ai+1) =|Zi-M|+|Zi+1-M|≤|Zi-M|+|Zi-1-M|=d (ai-1, ai) 。

则存在一个非负的随机变量 0 ≤ k (ω) <1 使得

d (T (ω, ai-1) , T (ω, ai) ) ≤k (ω) d (ai-1, ai) 。

且 (a, d) 是完备的度量空间。所以元胞蚂蚁算法形成的影射T为随机压缩算子。

定理3[21,22] 设随机算子T:a→a满足对几乎所有的ω∈Ω, T (ω) 均为压缩算子, 即存在Ω0⊆Ω, p (Ω0) =1, 使任一ω∈Ω0有

d (T (ω, ai-1) , T (ω, ai) ) ≤k (ω) d (ai-1, ai)

∀ai-1, ai∈a其中0≤k (ω) <1。

对任一ω∈Ω0, 则T (ω) 有唯一随机不动点g (ω) , 即T (g (ω) , ω) =g (ω) 。

结论元胞蚂蚁算法的求解迭代过程是一个随机压缩影射, 根据定理3可知, 该迭代过程存在唯一的不动点, 则元胞蚂蚁算法的求解迭代过程是收敛的。

4 数值实验

本文选用“http://www.iwr.uni-heidelberg.de/groups/comopt/software/TSPLIB95/”公布的标准问题库TSPLIB中的实例。在Windows 2000下用Delphi 7.0编制程序进行了计算, 各种参数的设定如下, α=1, β=1, Q =100, ρ=0.7, m=n (蚂蚁数等于城市数) 。并使用2-OPT算法进行邻域优化。下面是部分问题得到的结果, 其结果如表1所示 (在计算边长时四舍五入取整) , L1为本算法求得的最优解, L2为TSP库中公布的最好解, b1= (L1-L2) / L2, 表中问题名称后的数字为结点总个数。

从表1的数据中可以看到元胞蚂蚁算法在TSP问题求解过程总体效果良好, 算法最优解与TSP库中的最好解的误差基本控制在1%以下, 说明该算法可以作为一种解决组合优化问题的方法。

5 结束语

结合TSP问题给出了离散元胞蚂蚁算法的描述, 并指出此算法的搜索求解实际是一种随机压缩映射, 依据随机不动点理论证明了离散元胞蚂蚁算法具有唯一随机不动点, 即存在唯一随机解, 同时通过数值仿真证明算法的有效性。

基于离散元对放煤方式的研究篇5

我国引进综采放顶煤技术已有二十几年, 特别是最近几年我国大型煤矿的综采放顶煤技术得到了迅速发展, 并且在一些矿区取得了较好的技术经济效果。但由于我国厚煤层储量丰富, 分布较广, 而且各地条件各异, 因此在不同条件下采用的方法和工艺也将有所不同。

不同的煤层及开采条件选择合理的放顶煤工艺模式是实现工作面高产高效的重要保证[1]。本文在总结放顶煤研究的基础上[2,3], 运用离散元PFC3D软件、结合前人的理论基础, 对放煤方式进行模拟。由于篇幅有限, 集中对单轮间隔放煤、单轮顺序放煤和多轮间隔放煤3种放煤方式加以研究。

1 模拟方案

图1为工作面模型, 模拟顶煤厚度为5.5 m, 共有12个支架, 模型两侧各留有1架或2架不进行放煤, 支架参数设置:放煤口大小为1.5 m, 放煤步距为0.6 m, 模拟中低位放煤。分别模拟以下3个方案, 对不同放煤方式的放煤效果进行分析比较。

(1) 单轮间隔放煤:

先打开1#、3#、5#、7#、……单号支架上的放煤口放煤, 且见矸石关闭放煤口, 留下一定的脊背煤, 然后放双号支架, 将留下的脊背煤放出。

(2) 单轮顺序放煤:

放煤顺序按 1#、2#、3#、……放煤口顺序放煤, 每个支架放煤后见矸后关闭放煤口。

(3) 多轮间隔放煤:

先打开 1#、3#、5#、7#、……单号支架上的放煤口, 将1/2顶煤放出, 然后放双号支架, 同样也是放出1/2, 之后再打开单号支架放煤口放出剩余的1/2, 见矸石后关闭放煤口, 然后再次打开一次双号支架放煤, 都按见矸石后关闭进行。

2 模拟结果分析

离散元PFC3D软件是高级非连续介质程序, 适用于任何需要考虑大应变、破裂、破裂发展以及颗粒流动问题的场合, 可以用来研究结构开裂、堆石材料特性和稳定性、矿山崩落开采、边坡解体、散体崩落等一系列传统数值方法难以解决的问题。本文运用该软件模拟煤层下落过程, 形象地反映煤层在不同放煤方式下出现的不同跨落效果。

(1) 单轮间隔放煤

图2为单轮间隔放煤过程。从图2 (a) 可看出, 单轮间隔放煤过程是先打开单号支架放煤口, 放完单号支架上方顶煤后, 双号支架放煤口上方的余煤形态类似于三角形。单号支架上方放煤后形成规则的类似三角形的空间由上部红色球体填充, 从红色球体下落效果来看, 符合椭球体放矿理论。相隔2个单号支架上方存在一定间距, 形成双号支架上留置的脊背煤。在双号支架放煤口打开后, 双号支架放煤口上方三角状顶煤受到两侧矸石的水平力作用, 大部分从放煤口放出, 少部分由于受到上方压力产生垂直位移, 混入临近的单号支架放煤口上方而不能放出, 造成了顶煤的损失。从图2可看出, 单号支架放煤后形成了的空间没有被直接顶矸石填满, 当双号支架打开后, 上方三角状煤层的两翼就容易被挤压而填补单号支架上的空间, 这是间隔放煤造成煤损的原因之一。建议在放煤完毕, 对放煤口上方积煤多的支架可单独进行放煤, 这样可减少煤的损耗。

(2) 单轮顺序放煤

图3为单轮顺序放煤过程。图3 (a) 为打开第一个放煤口进行放煤结束时的煤层状态, 放煤口上方的矸石在顶煤落下后填补了留下的空间, 从图3 (b) 中看到在关闭第一个放煤口、打开第二个放煤口后, 在放煤口上方的煤层因为受到侧向矸石的挤压导致支撑能力下降, 表现出更明显的流动性, 由放煤口放出, 经过一段时间放煤, 相邻的第一个放煤口上方部分矸石也下落至放煤口。最后, 所有放煤口放煤结束后, 放出煤中含矸较少, 这主要是由于见矸石关闭放煤口造成的。由模拟得知, 采用这种方式放煤, 建议在顶煤硬度低的煤层中进行高位放煤, 对于低位放煤, 顶煤硬度可适当大一些。

(3) 多轮间隔放煤

多轮间隔放煤方式现在采用得较少。这种方式理论上煤岩界面下降平稳, 可以减少混矸, 顶煤成拱后可以破坏, 放出率较高。但操作过程较繁琐, 混入量和混矸层均有增加, 减慢了放煤速度。多轮间隔放煤过程如图4所示, 图4 (a) 中单号支架上方形成了类似三角状的空隙, 这时上覆煤岩压力主要集中在双号支架上方, 而双号支架上方的煤主要受重力作用, 所以打开双号支架放煤口后放煤速度较快。图4 (b) 为双号支架放完一半煤的效果, 经过这次放煤, 单号支架上方的煤层压紧。观察得知这是由于上方煤的下沉和相邻支架上方煤层的运移形成的。继续放煤直至最终放煤完毕 (见图4 (c) ) , 由于每次放煤时, 煤岩接触面都有移动, 临近顶煤放完时已经有矸石混入。这种方式在实际生产中不可取。

3 结语

(1) 本文对3种不同的放煤方式依据模拟结果加以讨论, 分析了3种不同放煤方式的优缺点, 提出了一些加以改进的意见, 达到了运用模拟手段指导生产的目的。实际生产中, 可针对矿山实际情况, 采取较为合适的放煤方式。

(2) 在本文的研究中, 模拟过程中煤层设置为均匀的散体, 这与实际放煤过程有些差异, 实际生产过程中煤层不可能破碎到均匀体下落, 需要进一步改进。

参考文献

[1]郭金刚.提高综放采出率的理论与技术[M].北京:煤炭工业出版社, 2002.

[2]孙二伟, 王洪胜, 王淑玲.综采放顶煤采放参数优化研究[J].现代矿业, 2009 (4) :26-28.

[3]吴健.我国放顶煤开采的理论研究与实践[J].煤炭学报, 1991 (3) :34-36.

[4]许延春, 张玉卓.应用离散元法分析采矿引起厚松散层变形的特征[J].煤炭学报, 2002 (3) :68-71.

离散元软件篇6

为简化地基处理, 加快施工进度, 工程中经常采用柔性堰体结构, 如格宾挡土墙、加筋土挡土墙等结构型式。石笼结构, 即格宾 (gabion) 结构, 是将符合粒径要求的石料填入具有柔性的铁丝笼中达到一定的孔隙率, 逐层砌筑的一种新型的柔性挡土构筑物[1]。这种结构可以承受苛刻环境和气候条件的破坏, 它具有透排水良好、坡面美观、柔性结构、贴地性佳、耐久性好等优点, 在国外得到了广泛的应用, 国内也有所应用, 例如作为护坡成功用于长江干堤[2]。

这种柔性堰体其基本工作单元为填石的石笼, 由单个石笼堆积而成, 每一个石笼独立工作, 具有很强的离散特性, 他们之间通过一定的连接而形成一个整体。工程上常用的数值计算方法主要是有限元方法, 然而它是以连续介质为基础的数值模拟方法, 在处理不连续面等方面有很大的困难, 因此如何建立合理描述具有不连续性的岩体的力学行为, 引起了很多学者的关注, 并提出了不同的数值计算方法, 如离散单元法 (DEM) 、块体元、边界元和DDA等, 这些方法大部分都能得到较满意的结果, 在应用方面以离散元方法最为广泛。

UDEC (Universal Distinct Element Code) 是一种二维通用离散元显示有限差分程序[3], 使用显式时间递步法, 适应大变形和旋转的工程计算, 主要模拟静载或动载条件下非连续介质 (如节理块体) 的力学行为特征。UDEC在计算时不用形成刚度矩阵, 具有计算稳定、效率高等优点, 因此得到了广泛的应用, 包括对地下洞室围岩稳定、重力坝抗滑稳定以及三峡船闸高边坡的静动力稳定性进行了分析[4,5,6,7,8,9,10]。

1离散元计算原理

在计算时采用动态过程描述系统的响应, 使用时间递步算法, 把演化过程分成很多个时间步, 假定每一个时间步足够小, 并且速度和加速度为常数, 在一个单独的时间步内, 扰动来不及在离散元块体之间传播。在分析连续体时采用显示有限差分算法, 时间步同时受到接触和块体的限制, 对于刚性块体, 通过块体质量和界面刚度定义时间步大小;对于变形体, 通过网格的大小和系统的刚度决定时间步的大小, 其中系统的刚度同时受到接触岩块的模量和接触刚度的影响。

在计算循环中交替应用牛顿第2定律和力-位移定律。牛顿第2定律用来计算块体的运动, 更新块体的位置坐标;而力-位移定律用来计算接触力, 使块体保持平衡, 如果块体是变形体, 计算的将是块体内节点的运动, 块体材料本构将给出新的单元应力。

2工程实例

向家坝水电站2期纵向围堰大坝上游段直接与2期横向土石围堰相接, 堰体主要起导水和挡土墙作用。采用格宾挡土墙结构型式, 堰体结构由单个格宾堆砌而成, 格宾之间纵横绑扎以形成整体柔韧结构。目前工程上常用的格宾是由低碳钢丝经机器编制而成的六边形双绞合金属网面构成, 钢丝表面采用热厚镀锌防蚀处理方式。单个格宾尺寸为1 m×1 m×0.5 m或1 m×1 m×1 m (长×宽×高) , 格宾内填充石料。堰体结构形式如下:底部宽度38.00 m, 堰顶宽度3.00 m, 两侧坡比为1∶0.50, 最大堰高35.00 m, 堰体底部设置1.50 m厚的加筋碎石垫层、1.50 m厚的碎石垫层[11]。

利用柔性堰体结构对地基的适应能力强, 承载力要求低的特点, 堰基底部回填采用碎石分层碾压, 提高承载能力, 减少不均匀沉降。

采用离散元对2期纵向围堰大坝坝体进行模拟, 每个格宾都被离散为一个单独的块体, 块体采用弹性材料, 块体之间采用库伦摩擦模型。每个块体采用三角形差分单元进行细分, 块体的接触检索在块体表面三角形差分单元中进行, 变形体的计算采用大变形有限差分方法。计算模型如图1、图2所示。

3计算结果及分析

此次计算的目的是了解坝体自身由于自重的变形情况, 以及受土压力荷载后坝体的应力和变形情况, 因此计算时将坝体从整个模型中分出来单独计算, 并考虑施工顺序分层计算, 共分35层, 依次激活坝体单元, 计算得到坝体的应力及变形结果;然后再在边界上等效施加土压力荷载, 计算考虑土压力对坝体应力变形的影响。

取块体的密度ρ=1 750 kg/m3, 弹性模量E=60 MPa, 泊松比μ=0.3, 每个块体之间将接触面作为边界条件相互作用, 摩擦符合库伦强度准则, 其摩擦系数取0.7, 节理法向和切向刚度均取0.1 GPa/m。

3.1自重条件下坝体计算结果

由于采用离散元方法, 块体之间可以脱开, 从而导致等值线不连续, 因此选用云图出示坝体计算结果。图3～图6分别显示了坝体水平位移、铅直位移、水平应力及铅直应力计算结果。

由图3～图6可知, 坝体应力位移与EB模型计算结果规律一致, 水平位移为底部朝外鼓, 其中关键点1和2为水平位移值最大处;铅直位移发生在大约一半坝高处, 即关键点4所在位置;铅直正应力最大值发生在坝体底部中间, 即关键点5所在位置, 其值约为0.4 MPa。

图7为坝体剪切位移图, 从图中可以看到坝体的剪切位移主要发生在坝坡的格宾结构, 由于不受外荷载作用, 由自重产生的位移呈左右对称, 最大值在1/3坝高处, 最大值为1.073 cm。由此可知坝体的最危险部分在两侧坝坡处, 这是一般连续力学方法难以考虑的, 体现了离散元的适用性及准确性。

3.2施加土压力后坝体计算结果

向家坝水电站2期纵向围堰堰体主要起导水和挡土墙作用。为了能够更准确地了解土压力对柔性挡土坝应力变形的影响, 计算时在坝右侧边界上施加等效土压力荷载, 其值由EB模型计算得到。顶部水平应力为1 814.2 Pa, 铅直应力为9 094.2 Pa;底部水平应力为134 200 Pa, 铅直应力为646 850 Pa, 中间以线性差值计算得到。

计算得到施加土压力后坝体的变形图, 如图8和图9所示。由图可知土压力对坝体铅直位移影响不大, 对水平位移有较大影响。在坝体右下侧, 由于受到土压力的影响, 其铅直位移明显增大, 但坝体铅直位移最大值位置不变, 依然在坝体中部。坝体水平位移由对称变为非对称形式, 其值基本为负值, 说明坝体受到土压力作用向左边滑移。

图10为坝体受土压力后剪切位移图, 从图10中可以看到坝体在受土压力后的剪切位移不再对称, 左侧坝坡的位移明显大于右侧坝坡, 最大值仍在1/3坝高处, 最大值为1.262 cm, 跟考虑土压力之前相比, 增大了1.89 mm, 发生在左侧坝坡处, 右侧坝坡由于受压, 格宾之间的摩擦力增大导致剪切位移减小, 因此对于左侧坝坡的格宾应加大之间的黏结力以减小其相互错动。

4结论

(1) 施加土压力前后坝体的沉降变化不大, 坝体水平位移发生较大变化, 从1.5 cm增大到5 cm, 发生在坝顶位置, 说明土压力主要起水平推力作用。

(2) 坝体水平位移为块体变形及块体之间剪切滑移共同作用的结果, 剪切滑移最大值为1.262 cm, 发生在左侧坝坡, 为保证坝体的稳定应对坝坡部分格宾加强黏结力。

(3) 使用离散元程序UDEC对石笼挡墙这种具有离散性的填石挡土结构进行分析, 可以比较准确地抓住这种挡土形式的柔性特点, 结果更加合理。

摘要：针对柔性挡土坝的结构特性, 采用离散元程序UDEC对向家坝水电站2期纵向围堰挡土坝进行了数值模拟, 得到了坝体在挡土前后的应力位移规律, 通过对关键点部位的应力及位移变化情况的记录, 更清楚地了解坝体的变形情况和受力特点。计算结果表明使用离散元分析这种柔性挡土结构是可行的, 可以比较准确地抓住这种挡土形式的柔性特点, 并能较好地模拟挡土坝发生的较大的法向变形和切向滑移, 能够对构筑物的稳定性做出合理判断, 对于施工应采取的措施提出了有益的参考。

关键词：UDEC,离散元,柔性挡土坝,石笼结构,向加坝,土石围堰

参考文献

[1]孟云伟, 肖文, 柴贺军, 等.柔性石笼挡墙变形受力的PFC (2D) 数值模拟[J].地下空间与工程学报, 2007, (2) .

[2]张焕洲, 谢平, 戴秋红.格宾网材在黄石长江干堤合兴堤段的应用[J].人民长江, 2002, 33 (9) :38-40.

[3]侯艳丽.砼坝-地基破坏的离散元方法与断裂力学的耦合模型研究[D].北京:清华大学, 2005.

[4]Zhang Chuhan, Pekau O A, Jin Feng, et al.Application of distinctelement method in dynamic analysis of high rock slopes and blockystructures[J].Soil Dynamics and Earthquake Engineering, 1997, 16 (6) :385-394.

[5]Iofis IM, Maksimov AV, Mironov VV.Some Practical aspects ofnumerical simulation of jointed rock mass by distinct elementmethod[R].Proceedings of the International Conference on Me-chanics of Jointed and Faulted Rock, 1990:479-486.

[6]雷远见, 王水林.基于离散元的强度折减法分析岩质边坡稳定性[J].岩土力学, 2006, (10) .

[7]王泳嘉.离散单元法一种使用于节理分析的数值方法[C]∥第一届全国岩石力学数值计算及模型试验讨论会论文集.成都:西南交通大学出版社, 1986:32-37.

[8]张楚汉.论岩石、混凝土离散-接触-断裂分析[J].岩石力学与工程学报, 2008, (2) .

[9]朱焕春, Andrieux P, 钟辉亚.节理岩体数值计算方法及其应用 (二) :工程应用[J].岩石力学与工程学报, 2005, (1) .

[10]谢和平, 陈忠辉, 周宏伟, 等.基于工程体与地质体相互作用的两体力学模型初探[J].岩石力学与工程学报, 2005, (9) .

[11]常晓林, 周伟.向家坝水电站二期纵向围堰结构计算报告[R].武汉:武汉大学水利水电学院, 2007.

离散元软件篇7

砌块是由各种散体如粉煤灰、炉渣、砂石、浮石、火山渣、天然煤矸石等作为骨料和胶凝材料经机械成型、养护而成的块材,近年来得到了广泛的应用。砌块的一个主要指标就是抗压强度,而抗压强度在很大程度上取决于砌块的密实度。如何提高砌块的密实度,砌块在振动加压下是如何密实的,颗粒在振动状态下的受力和运动规律如何,到目前为止很少有学者在理论上做过深入的研究,主要是因为传统的方法很难对大量的散体进行受力和运动规律计算和模拟。近来,利用离散元法对散体进行数值模拟的研究工作取得了很大进展,离散元法的基本思想最初于20世纪70年代由Cundal1[1]等人提出的,用来研究散体或非连续体的力学问题。Zhou[2,3]等用改造的DEM程序对玻璃珠的堆积进行了数值模拟并与实验结果对比。周健[4]等对砂土的室内双轴试验和砂土剪切带形成与发展进行了数值模拟,研究了模型试样颗粒粒径、颗粒刚度和摩擦系数等细观参数的变化,对剪切带的形成与发展过程的影响。离散元法在求解散体的力学问题,特别是对于研究散体颗粒的运动,显示出巨大的生命力。本文基于离散元理论,对砌块成型过程中颗粒物料的运动规律和受力情况进行了模拟,分析了不同参数下颗粒物料在振动作用下的孔隙度的变化。

2 散体颗粒振动压实机理分析

目前振动压实技术的基础理论研究可归纳为以下几种学说[5]:

a.内部摩擦减小学说:由于振动作用使被振压材料的内部摩擦阻力急剧减小,剪切强度降低,抗压阻力变得很小,材料在重力作用下易于压实。

b.共振学说:当激振频率与压实材料的固有频率一致时,振动压实最为有效。实践证明利用共振现象进行压实效果显著。但由于材料的固有频率是变化的,要求激振器的频率做相应的变化比较困难。

c.反复载荷学说:振动所产生的周期性压缩运动的作用,使被压材料受反复载荷作用而达到振动压实的目的。实验表明:该学说在低频率范围内有一定的现实性,而在高频率范围内并无充足的理论依据,在高频范围内,振动压实的效果远远超过反复载荷的作用效果。

d.交变剪应变学说:利用土力学交变剪应变原理,振动使颗粒产生剪应变,使被压材料的颗粒重新排列而达到密实效果。

以上各种学说从不同的角度对振动压实机理进行解释,然而砌块骨料种类繁多,实际工况各异,且被压实材料多为非匀质异性材料,影响压实的因素很多,因此各学说都有待进一步的试验研究和探讨。

对于砌块而言,散体颗粒的摩擦系数、级配、振幅对于振动压实效果有明显的影响,从波动方程来看振动压实对散体中剪切应力的影响与被压实散体的种类无关。振动作用使被振压材料的内部摩擦阻力急剧减小,剪切强度降低,抗压阻力变得很小,材料在压力作用下易于压实。散体的抗剪强度τf可由库仑定律公式表示为:

式中c—颗粒凝聚力;σ—法向应力;准—颗粒内摩擦角。

对于砌块骨料这样的理想散粒而言,其c值很小,可以忽略不记,它在剪切过程中形成的剪阻力主要是颗粒间的摩擦力,其抗剪强度可以简化为:

从细观的角度而言,其颗粒间的摩擦还可以细分为颗粒间的滑动摩擦阻力和颗粒间的咬合摩擦阻力。滑动摩擦力是颗粒接触面粗糙不平形成的微细咬合作用。颗粒间距离的微弱增长,会使微细咬合作用产生很大的衰减。如果振动能使颗粒质点间的距离产生微弱的增长,就会使滑动摩擦力减小。振动过程中,散体振动加速度a为:

式中A—振幅;ω—激振频率;β—相位角;t—时间。

散体颗粒的惯性力I为

式中mk—散体颗粒质量。令e=Aω2

称e为振动强度。则有

由式(4)可以看出,散体颗粒惯性力I与颗粒质量和振动强度成正比。当振动强度较小时,或散体颗粒质量较小时,散体颗粒的惯性力也较小,它将在自己原来的位置振动。当振动强度和散体颗粒质量足够大时,散体颗粒的惯性力足以克服周围其他散体颗粒凝聚力的作用,使散体颗粒偏离自己原来的位置,具有良好级配的散体颗粒,相邻散体颗粒间的粒径大小不同,即它们的颗粒质量不同,因此相邻散体颗粒在具有相同振动强度时,惯性力大小不同。这种差别必然会使颗粒质点间的距离发生微小的变化,对颗粒间的微细咬合作用产生很大的衰减,导致内摩擦角的减小,即内摩擦力tan准减小,当减小到某个值的时候,散体颗粒间的摩阻力不足以抗拒振动带来的惯性力,使散体结构发生重排列。

散体颗粒质量大,凝聚力小,振动可以显著减小散体颗粒的滑动摩擦,对散体颗粒的抗剪强度影响很大,对于颗粒均一的散体,在同样的振动强度时,颗粒受到的转动惯性力大体相同,使得颗粒会出现定向排列,定向排列的结果是在某个方向上颗粒均为面-面接触,颗粒在这个方向上接触应力很小,抵抗外力的能力大为加强,这也是振冲产生预振效应的原因。对于颗粒级配良好的砌块散体颗粒而言,由于其颗粒大小形状比较复杂,颗粒排列在宏观统计上不易出现定向性,但是由于粗颗粒骨架的影响,细小颗粒多与粗颗粒骨架形成面-面接触,一样表现出很明显的预振效应,而且级配良好的颗粒相对于颗粒均一的颗粒而言更易达到较高密实度。

3 振动压实过程模拟

3.1 基本假设

密实度是指材料体积内,被固体物质充实的程度以D表示,孔隙度是指材料体积内,孔隙体积所占有的比例.以P表示,孔隙度与密实度从两个不同侧面来反映材料的致密度,D+P=1,为了仿真过程中便于测量,采用空隙度来表示砌块的致密程度。

为了模拟颗粒散体在振动压实下的孔隙度变化规律,在模拟过程中做了如下假设:

a.将砌块骨料颗粒等效为两种粒径颗粒,粒径比为1:6。

b.颗粒单元为刚性体不存在压溃。

c.颗粒之间的接触为点接触,颗粒之间的接触存在重叠量,但对于自身尺寸而言很小,重叠量的大小与接触力有关。

d.颗粒单元为球形颗粒。

3.2 振幅与砌块不同高度孔隙度关系模拟

以标准砌块尺寸190 mm×190 mm×390 mm,颗粒由直径大小分别为6 mm和0.5 mm组成,颗粒与颗粒之间以及颗粒与墙壁之间的摩擦系数为0.6,颗粒比重为2.65,切向和法向刚度都为1.0 E+09 N/m,图1给出了在颗粒摩擦系数、级配和振动频率一定的情况下,通过调整振幅的大小来观察砌块在不同高度孔隙度的变化规律。振动的施加引起抗剪强度急剧下降,随着振幅的增加,抗剪强度逐渐降低,使颗粒重排,孔隙度减小,由于振波在传递过程受到空间散射、颗粒间内摩擦、颗粒的黏滞性等因素形象衰减很严重。振幅越大,振动能量越大,传递的高度也就越高,克服周围凝聚力颗粒越多,颗粒重排力度就越大,孔隙度变得越小;反之振幅越小,只使最底层颗粒产生颗粒重排作用,对25 mm以上的颗粒作用微乎其微。

3.3 摩擦系数与孔隙度关系模拟

在颗粒振动频率、级配和振幅一定的情况下,通过调整颗粒间的摩擦系数来模拟砌块在不同高度孔隙度的变化规律。从图2可以看出,摩擦系数为0.5时,使底层颗粒出现了振动液化[6],颗粒对振动的阻尼减小,增加了振动传递的距离,从图2可以明显看出孔隙度低于其他曲线,其他曲线在60 mm以下孔隙度的变化趋势几乎一致只是大小差异较大。在60 mm以上,孔隙度随摩擦系数的变化不太明显。

3.4 颗粒级配与孔隙度关系模拟

在颗粒摩擦系数、振幅和振动频率一定的情况下,通过调整级配来观察砌块在不同高度孔隙度的变化规律。图3给出了A、B、C、D等四种级配在振动条件下不同高度的孔隙度变化,B、C两种级配的孔隙度明显低于其他两种情况,良好的级配可以增加砌块的密实度。

图4给出了颗粒自动堆积状态,通过测量图中1,2位置的孔隙度在振动条件下随时间的变化可以得出不同级配与孔隙度的关系,对比图5、图6可以看出图4中位置2孔隙度明显小于位置1,因为散体颗粒的惯性力与其质量成正比,当散体颗粒质量较小时,散体颗粒的惯性力也较小,它将在自己原来的位置振动,而质量大颗粒足以克服周围其他颗粒内聚力的作用,将偏离自己原来的位置,使颗粒重新排列,小颗粒填充孔隙,提高了砌块的密实度。同时也可以看出振动时间越长未必孔隙度越小。

3.5 双向振压与孔隙度关系模拟

由前面的分析可知,单向振动随高度的增加,孔隙度变得越来越大,也就是说离开振源较远的地方,振动密实的效果就越差。基于此,我们提出采用双向振动双向加压的方式,试图提高砌块的密实度,从而生产高强度的砌块。表2为模拟时所取参数,图7为双向振压的模拟效果图,孔隙度采用同一高度多点测量取均值。从图8中可以看出双向振压的方式孔隙度明显低于其他成型方式。由于振波在散体中传播的距离有限,采用下部振动不能够使上部颗粒产生足够大的惯性力来突破周围散体的凝聚力而产生滑动。从图8中单纯的压实虽然是孔隙度有所降低,但是相对于振动后压实来说孔隙度还是偏高的,采用两个方向振动方式使颗粒产生错位重排,填充孔隙,最后加压得到更小的孔隙度。

4 结论

通过对砌块振动成型方式几个重要参数离散元法模拟,至少得到以下4点结论:

a.振幅的大小对砌块孔隙度影响较大,大振幅可以使离振源比较远的颗粒获得足够的惯性力来突破周围散体颗粒的凝聚力,使颗粒产生滑动或转动,促使颗粒重排,在压力的作用下得到更加小的孔隙度,从而提高砌块的强度。

b.振动压实条件下摩擦系数越大,颗粒滑动需要的能量就越大,同时吸收的振动能也就越多,只有通过加大振幅和加大压实力才能得到较小的孔隙度。

c.级配在振压条件下对砌块的孔隙度影响比较明显,合理的级配更加容易得到较小的孔隙度。级配不合理即使振幅再大也很难得到理想砌块密实度。

d.采用双向振压方式,由于振源从上下两个方向对颗粒进行能量传递,使更多的颗粒在惯性力下发生了颗粒重排,得到了更加理想的孔隙度,进一步提高了砌块的强度,用离散元模拟砌块成型过程中振动参数及组分参数对孔隙度的影响,是一种非常有效的方法,取得了良好的模拟效果,分析结果对设备参数的选择具有一定的指导意义。

参考文献

[1]Cundall P A,Strack O L.A discrete numerical model for granular assembles.Geotechnique,1979,29(1).

[2]Zhou Y C,Xu B H,Yu A B,et a1.An experimenta1and numerica1study of the ang1e of repose of coarse spheres[J].Powder Technology,2002,125(1).

[3]Zhou Y C,Wright B D,Yang R Yet al.Rolling friction in the dynamic simulation of sand pile for-mation[J].Physica A,1999,269(2-4).

[4]周健,池永.砂土力学性质的细观模拟[J].岩土力学,2003,24(6).

[5]王戈,王贵慎,张世英.压实机械〔M」.北京:中国建筑工业出版社,1992.