颗粒物质局部变形的离散元模拟

中国科学 G辑: 物理学 力学 天文学 2008年 第38卷 第6期: 692 ~ 703

www.scichina.com phys.scichina.com

《中国科学》杂志社

SCIENCE IN CHINA PRESS

颗粒物质局部变形的离散元模拟

王等明, 周又和*

兰州大学西部灾害与环境力学教育部重点实验室, 兰州 730000 * 联系人, E-mail: zhouyh@lzu.edu.cn 收稿日期: 2007-06-06; 接受日期: 2008-04-08 国家自然科学基金重点项目资助(批准号: 10532040)

摘要 在外加载荷作用下的变形分析是目前颗粒物质研究的热点问题之一. 利用改进的离散单元法, 对随机分布的颗粒物质在平面应变状态下的变形进行了研究, 并分析了初始空隙比、围压、摩擦系数等因素对颗粒物质局部变形特征和整体稳定性的影响. 数值模拟结果表明: 在平面应变状态下, 颗粒物质内部会出现两个相互交叉的剪切面, 并且局部变形特征会随着竖向应变的增大而发生变化, 这些特征和目前存在的实验结果具有较好的一致性. 同时, 沿着两个剪切滑移面的切向速度廓线表现为两种不同的模式, 一种呈类流态模式, 另一种呈类固态模式. 最后对各种材料参数和因素对颗粒物质局部变形特征的影响进行了详细讨论.

关键词

颗粒物质 局部变形 改进的离散单元法

颗粒物质是由大量离散粒子组成的集合体, 在集合体内部, 各个粒子可以的移动且粒子之间只有在接触点产生作用力, 这种离散的特性导致颗粒物质表现出一些不同于一般固体和液体的复杂行为. 为了了解和掌握颗粒物质的这些独特性质, 近十年来很多学者从不同领域对这个问题进行了相关研究. 在这些研究中, 在外部荷载下的变形由于具有非常重要的工业应用背景而成为目前颗粒物质研究中的热点问题之一, 这些应用主要包括土木工程中的地基处理、边坡和堤坝的稳定行分析、颗粒物质的输运和处理、颗粒路面的设计和维护等[1,2]. 目前这方面的研究主要集中在实验、理论和数值模拟三方面.

在实验方面, 迄今为止, 对颗粒材料在不同边界和加载条件下的变形特征已开展了大量研究[3~5]. 这些实验主要是观测颗粒物质的抗压强度、破坏模式以及宏观的应力应变关系等, 但是导致颗粒物质发生局部变形的微观机理目前还不清楚. 一个主要的原因就是要通过实验手段观测颗粒物质内部的微观过程, 目前还存在技术上的挑战[3].

理论分析是目前研究密集颗粒物质变性特征的一种比较常用的方法. 这种方法主要是从连续体的本构理论出发, 借助有限元和其它数值方法来解决颗粒物质在不同边界和加载条件下的变形问题[6~9]. 在此方法中, 选取的本构模型对最后的结果起着非常重要的作用. 但是,

692

中国科学 G辑: 物理学 力学 天文学 2008年 第38卷 第6期

目前还没有一个合适的本构模型可以完全描述颗粒物质的性质[7]. 现有的关于颗粒物质的本构模型, 大多数都是以弹塑性本构理论为基础, 很少考虑颗粒物质内部粒子之间的微观作用, 而粒子之间的微观作用对颗粒物质的整体变形起着决定性的作用. 故理论方面的研究并不能完全真实的反映颗粒物质的变形特征.

为了能满足当前工程的需要, 许多研究者借助于数值模拟来解决目前存在的涉及到颗粒物质的问题[5], 其中离散单元方法就是一种最常用的数值方法. 最早利用数值方法来模拟颗粒物质中局部变形和剪切带的是Iwashita等人[9,10], 他们发现, 在准静态加载下, 颗粒物质基本的微观变形机理是在硬化过程中柱形结构的产生和软化过程中它的破坏. 虽然传统的滑动模型在这个研究领域已经被广泛接受, 但一些研究也显示, 滚动效应在颗粒物质的微观变形中起也起着非常重要的作用[9~11]. Hu等人[12]研究了密集规则排列的颗粒物质剪切带的形成机理和物理参数对剪切带的影响.

虽然目前开展了很多关于颗粒物质局部变形特征的相关研究, 但不管是实验、理论还是数值模拟, 都缺乏一种理解, 即各种物理因素和材料参数如何联合起来影响颗粒物质的局部变形特征、剪切面的出现及其演化规律等. 本文基于改进的离散单元法, 研究了随机排列的颗粒物质中剪切面的形成机理及其变形模式, 并分析了空隙比、围压、摩擦系数等对局部变形特征和整体稳定性的影响规律.

1 改进的离散单元法

离散单元法最早是由Cundall于1971年提出的一种不连续的数值方法. 它的基本思想是把研究对象离散为刚性单元的集合体, 使各个刚性单元满足各自的运动方程, 用时步迭代法求解运动方程就可得到不连续体的整体运动形态[13]. 离散单元法允许单元间的相对运动, 不一定要满足位移连续和变形协调条件, 故它可以比较有效的模拟诸如大变形、非线性以及多物理场作用下颗粒材料的复杂运动学和力学特性.

1.1 颗粒的运动方程

离散单元方法是基于牛顿第二运动定理和颗粒在接触点上的力与位移之间的本构关系所建立的数值模拟方法. 对任意一个粒子i, 其运动的控制方程可以表示为[14,15]

Ii

mi

dvic

=∑Fji+Fi, (1) dt

dωic

=∑(Ri×Fji+Mθji)+Mi, (2) dt

其中mi和Ii分别为第i个粒子的质量和转动惯量; vi和ω i为第i个粒子的平动速度矢量和转动

c

角速度; Fi和Mi为外部施加在第i个粒子上的力和力矩; Fji为第j个粒子对第i个粒子的接触

力; Ri为第i个粒子的中心到接触点的矢量, Mθji为颗粒i和颗粒j之间由于接触面上法向力的分布不均匀而产生的滚动抵抗力矩.

693

王等明等: 颗粒物质局部变形的离散元模拟

1.2 颗粒之间的接触力分析

颗粒物质内部两个粒子之间的相互作用, 可由两个粒子之间的相对距离来描述, 如图1所示.

图1 颗粒之间相互作用示意图

如果两个粒子之间发生接触, 即满足以下条件:

c

其中Dij为粒子i,j质心之间的距离, 则两个粒子之间会产生接触力Fji, 此接触力可以分解为cncs

和平行于接触面的切向接触力Fji, 即 垂直于接触面的法向接触力Fji

δn=(Ri+Rj)−Dij>0, (3)

ccncscncsFji=Fji+Fji=Fjinji+Fjisji, (4)

其中nji和sji分别为粒子i, j之间的接触面的法向矢量和切向矢量. 颗粒之间的接触模型是离散单元法的核心, 目前应用比较广泛和成熟的是在颗粒之间施加变形元件, 例如弹簧、阻尼器

cn

等[16,17]. 第j个粒子施加在第i个粒子上的法向接触力Fji可以表示为

cn

Fji=knδn−cnme(vr⋅nij), (5)

其中kn为法向弹性常数, cn为法向阻尼系数, vr为粒子之间的相对速度, me为有效质量, me=mimj

()(mi+mj). 这里粒子之间的接触模型采用了弹簧-阻尼器模型. 弹簧的弹性力模

拟粒子相互接触时的排斥力, 而相应的阻尼力模拟粒子碰撞过程中能量的损失. 第j个粒子施

cs

加在第i个粒子上的切向接触力Fji可以表示为

cscn

Fji=−sign(δs)min(ksδs,µFji)−csmevr⋅sij, (6)

()其中ks为切向弹性常数, cs为切向阻尼系数, µ为摩擦系数, δ s为接触过程中的切向位移. 当切

cn

时, 滑动摩擦力将会取代切向弹性力. 向弹性力ksδs大于滑动摩擦力µFji

694

中国科学 G辑: 物理学 力学 天文学 2008年 第38卷 第6期

1.3 颗粒之间的滚动抵抗力矩

传统的离散单元法中都假定粒子的旋转是由切向接触力单独决定的, 然而最近一些实验和数值模拟结果表明, 颗粒之间的滚动效应尤其是滚动抵抗力矩, 会对颗粒物质的微观变形产生重要的影响[9~11]. 当颗粒i,j相互接触时, 在接触点附近会产生接触力F, 两个颗粒在接触点附近会发生变形, 这样会使得两粒子之间的接触点逐渐演化为接触面. 由于切向接触力不为零, 法向接触力关于接触面的中心是不对称分布的, 从而会导致一个力矩产生, 如图2所示, 这个力矩称为滚动抵抗力矩[11].

图2 滚动抵抗力矩产生示意图

目前对于粒子之间的滚动抵抗的模拟主要采用弹簧-阻尼器模型, 即在粒子之间的接触点上添加扭转弹簧和阻尼器. 滚动抵抗力矩的表达式, 采用和切向方向类似的Mohr-Coulomb准则给出, 即

⎧⎪Mθ=Kr∆θ,⎨nn⎪⎩Mθ=ηF=αBF,

Kr∆θ<ηFn,

else,

(7)

其中Kr为扭转刚度, ∆θ 为粒子之间的相对转角, η 为滚动摩擦系数. 可以看出, η 具有长度单位, 由于滚动抵抗力矩和粒子之间的接触面积有关系, 故可以方便的将η 分解为两部分, 一部分为接触宽度B, 另一部分为无量纲的参数α, 此参数用来控制滚动抵抗力矩的大小. α 主要取决于接触力在接触表面的分布, 如果粒子的横截面是圆形, 则从物理意义上来说α 应当小于1. 通

过理论分析[11], α 的值可近似取为1/3.

2 数值模拟及结果分析

2.1 数值模型

为了模拟颗粒物质的局部变形特征, 必须选择一个合适的数值模型. 本文建立图3所示的平面应变状态下的加载变形模型. 矩形颗粒集合体的变形通过4个边界来控制: 下边界为固定的刚性平板, 左右边界施加大小恒定的静水压力(围压)来确保颗粒物质初始的稳定性. 顶端施加一个加载平板, 平板以给定的速度竖直向下运动. 当竖向的应变或应力达到其极限值后, 颗粒物质内部会产生局部的变形, 在某些条件下会产生明显的剪切面或剪切带.

695

王等明等: 颗粒物质局部变形的离散元模拟

整个颗粒集合体包含4000颗左右的圆形颗粒

(颗粒数目会随着孔隙比的不同而有所变化), 粒径在1~3 mm之间满足均匀型概率分布, 中间粒子的粒径用D表示. 初始颗粒样本的生成采用粒子沉积法, 即初始时刻让所有粒子在一定的尺度范围内随机分布, 并给所有粒子都赋予大小和方向都随机分布的初速度, 让其在重力作用下自由沉降. 当颗粒集合体达到稳定后, 再在其上部施加一定的均布压力, 压力的大小可以根据所需颗粒物质的初始空隙比来调节. 当上部施加的压力使颗粒集合体达到所需的高度和空隙比后, 需要通过对颗粒样本的逐步放松来消除内部的残余应力[18].

在本文中, 模型初始的几何尺寸为: 宽度l = 88 mm、高度h = 178 mm. 数值模拟中各物理参数取为[16,17]: 颗粒的密度ρ = 2850 kg/m3、粒子之间的

图3 平面应变状态下颗粒物质受压模型

摩擦系数µ = 0.3、粒子与平板之间的摩擦系数µ1 = 0.35、法向的弹簧刚度kn = 1.0×107 N/m、切向的弹簧刚度ks = 2/3·kn、扭转弹簧刚度kr = 1.3×102 Nm/rad、粒子之间的法向阻尼系数cn = 8.33×10−2 N·s/m、切向阻尼系数cs = 8.33×10−2 N·s/m、扭转阻尼系数cr = 1.0×10−2 N·ms/rad、计算时间步∆t = 10−6 s.

2.2 颗粒物质局部变形的典型特征

与一般的固体材料的变形不同, 颗粒物质在外部载荷作用下的变形具有明显的局部特性, 即在某些局部的位置会产生较大的变形, 而在其余位置其变形非常微小, 这种特殊的性质会使得颗粒物质内部产生明显的剪切面(一侧变形速率较大而另一侧变形速率非常小的面)或剪切带(变形速率的变化集中在一个相对较窄的区域内). 不同于一般的环形剪切或平面剪切[19], 颗粒物质在顶部受压时一般会沿着顶端产生两个相互交叉的剪切面. 定义竖向应变为

ε=

∆h

, (8) h

其中∆h为上部加载板的竖向位移, h为颗粒集合体的原始高度. 图4为中间粒径D = 2 mm, 围压σc=15 kPa、顶部加载速率v = 0.3 m/s时, 颗粒物质在不同竖向应变阶段的速度矢量图. 图4中每个矢量的方向代表了单个粒子的运动方向, 而矢量的长度则代表了粒子运动速度的相对

大小. 从图4中可以看出, 在这种加载模式和边界条件下, 沿着颗粒集合体顶端会产生两个相互交叉的剪切面, 并且剪切变形的形式和Khalid等人[3]通过实验得到的局部变形特征具有较好的一致性.

本文采用的加载方式是给颗粒样本顶端施加竖向的应变. 当此应变达到其极限值后, 颗粒物质内部会产生明显的局部变形, 这将直接导致颗粒集合体对上部加载平板反作用力的变化. 将此作用力经过统计平均后, 就可以得到加载平板上的应力分布. 由于这种应力分布是与

696

中国科学 G辑: 物理学 力学 天文学 2008年 第38卷 第6期

图4 颗粒物质在不同竖向应变时的速度矢量图

颗粒物质内部的局部变形直接相关的, 故可用其来反映颗粒物质内部局部变形的出现及演化

规律等. 这里采用的预测参数η, 取为顶端平板上的平均应力与围压的比值, 即

η=

σv

, (9) σc

其中σv为顶部平板上的平均正应力, σc为水平围压. 图5为一个典型的应力比随着竖向应变的变化关系图. 从图5中可以看出, η 刚开始逐渐增大, 当达到其峰值后开始下降. 数值模拟结果表明: 应力比峰值对应的竖向应变, 刚好是颗粒物质内部出现局部变形时的竖向应变. 在此之前颗粒物质有一个整体的压缩过程, 在这个过程中并未出现明显的局部变形, 这个阶段称为颗粒材料的硬化过程, 而应力比下降阶段为颗粒材料的软化阶段, 在这个阶段颗粒物质内部会出现明显的局部变形, 并伴随着剪切面的形成和发展. 此现象也得到了实验的验证[3]. 图6为ε = 0.017时颗粒物质内部的切向力链图, 其中链的粗细代表了粒子之间切向作用力的相对大小. 从图中可以看出, 在剪切面附近尤其是两个剪切面的交界部分, 切向力明显增强. 数值模拟结果中并未发现法向接触力在剪切面附近有明显增强的现象, 这说明平面应变状态下, 颗粒物质发生局部变形并形成明显的剪切面时, 颗粒之间的剪切力起着非常

图5 应力比随竖向应变的变化

697

重要的作用.

王等明等: 颗粒物质局部变形的离散元模拟

为了分析剪切面的几何特征和演化规律, 图7给出了颗粒物质变形时剪切面的几何示意

图. 由于宏观剪切面是由于样本中部分微观粒子的集体运动导致的, 故可以根据颗粒物质中局部粒子的运动特征来反映剪切面的演化特征, 并可以定义剪切面的倾斜角等. 本文采用沿某个方向满足一定条件粒子(沿x正向的速度大于所有粒子沿x方向的速度的平均值的粒子)的

运动方向与x轴所夹角度的平均值, 来定义剪切面的倾斜角为

1θ=

Nc

∑wciθci,

i=1

Nc

这里θci为满足条件粒子的运动方向和x轴所夹的角度, wci为θci对应的权值, 可用来区分不同运动粒子的速度在剪切面倾斜角定义中所占的比重. 在定义了两个方向不同的剪切面倾斜角以后, 可以分别分析两个剪切面对应的变形特征、剪切面的倾斜角随各种因素的变化规律等.

图6 切向力链示意图 图7剪切面的几何示意图

按照上面剪切面倾斜角的定义, 可将整个速度矢量图沿两个剪切面的切线方向进行分解, 图8(a)为D = 2 mm, ε = 0.017时颗粒物质整体剪切变形时的速度矢量图, 图8(b)和(c)为相同条件下沿两个剪切面切线方向分解后的速度矢量图, 从图中可以明显的看到, 在此应变下, 颗粒物质内部形成了两个明显的剪切面.

为了分析剪切面在颗粒物质变形中的演化特征, 图9(a)和(b)分别给出了沿两个剪切面的法

线方向上, 剪切速度廓线随竖向应变的变化规率. 图9中y1和y2为两个剪切面法线方向的局部

坐标, 这里对剪切速度和局部坐标都进行了归一化. 可以看出, 两组速度廓线反映出了两种不同的剪切变形模式. 首先从图9(a)可以看出, 沿着第一个剪切面法线方向, 从加载位置开始,

剪切速度衰减很快, 类似负指数分布. 而从图9(b)可以看出, 除去加载附近一个很小的区域外, 剪切速度类似于线性分布, 呈扩散剪切. 在某些时候会出现明显的剪切带, 即速度廓线中间会

698

中国科学 G辑: 物理学 力学 天文学 2008年 第38卷 第6期

图8 剪切变形沿两个剪切面方向的分解

(a) 整体变形时的速度矢量; (b) 沿第一剪切面方向分解后的速度矢量; (c) 沿第二剪切面方向分解后的速度矢量

图9 不同竖向应变下的切向速度廓线图

(a) 类流态模式; (b) 类固态模式

出现一段速度变化很剧烈的区域. Aharanov等人[19]模拟平面剪切流时也发现存在两种典型的

剪切变形模式, 分别类似于流体变形的剪切模式(低围压, 高剪切速度)和类似于固体变形的剪切模式(高围压, 低剪切速度). 为了区分起见, 将速度呈负指数形式衰减的剪切模式称作类流态模式, 简称模式F. 将速度呈线性形式分布的剪切模式称作类固态模式, 简称模式S. 本文通过数值模拟发现: (ⅰ) 在平面应变状态下, 竖向加载时颗粒物质的剪切面不是对称分布的, 而

699

王等明等: 颗粒物质局部变形的离散元模拟

会像平面剪切流一样, 出现两种不同的局部变形模式, 即类流态模式F和类固态模式S, 但与平面剪切流也有不同之处, 即这两种变形模式是同时出现在颗粒物质的整体变形中; (ⅱ) 在整个竖向加载过程中, 两种变形模式保持一定的持续性, 其中类流态的剪切变形在整个变形过程中起主导作用, 一方面, 从剪切面的倾斜角看, 类流态的剪切面的倾斜角始终要大于类固态的剪切面的倾斜角, 而倾斜角越大, 颗粒物质越容易发生流动; 另一方面, 只要颗粒物质发生局部流动, 类流态的剪切面就会出现. 同时变形在单剪切面和双剪切面模式之间演化时, 类流态模式伴随着整个变形的过程.

以上结论说明在平面应变状态下, 颗粒物质在局部变形时呈现出了流体的特性. 但是由于围压的存在, 在此过程中, 颗粒物质又会表现出一些固体物质变形的特征. 2.3 各种材料参数和因素对颗粒物质变形的影响 2.3.1 初始空隙比的影响

颗粒物质的初始空隙比对其局部变形有很大的影响, 图10给出了不同初始空隙比下应力

比随竖向应变的变化关系, 从图10中可以看出, 初始空隙比越小, 即初始的颗粒物质越密集, 其内部开始出现局部变形所需的竖向应变越小. 同时, 初始空隙比越小, 应力比的峰值就越大. 这说明越密集的颗粒物质, 在相同的围压下, 只需要较小的竖向应变就能发生明显的局部变形, 但是可以承受更大的加载力. 这个现象与Khalid等人[3]的实验结果在定性上是一致的. 出现这种现象的原因, 可能是颗粒物质比较疏松(空隙比越大)时, 在发生局部变形前存在整体的竖向压缩. 在相同的围压下, 只有当空隙比达到其极限值后, 颗粒物质内部才会出现明显的局部变形. 数值模拟结果表明, 在其他条件完全相同的情况下, 颗粒物质发生局部变形时的空隙比是基本不变的. 若初始空隙比越大, 在加载过程中要达到其极限值所需要的竖向应变也越大. 正是由于发生局部变形前存在整体的竖向压缩, 所以会出现颗粒物质初始空隙比越大, 其发生局部变形所需的竖向应变越大的现象.

图11给出了初始空隙比对剪切面倾斜角的影响关系. 可以看出, 颗粒物质越密集(空隙比越小), 剪切面的倾斜角越小. 同时也可以看出, 类流态模式对应的剪切面倾斜角要大于类固

图10 不同初始孔隙比下应力比的变化 图11 初始孔隙比对剪切带倾斜角的影响 700

中国科学 G辑: 物理学 力学 天文学 2008年 第38卷 第6期

态对应的剪切面倾斜角, 这种现象与实验结果在定性上也是一致的[3]. 2.3.2 围压的影响

当不考虑颗粒之间的黏性时, 需要施加一定的围压来保持颗粒物质的整体稳定性. 图12

给出了不同围压下, 顶部平板的应力随着竖向应变的变化关系. 从图中可以看出, 围压越大, 峰值应力越大, 也就是颗粒物质要发生局部变形所需要施加在顶部平板上的应力越大. 这充分反映出围压对颗粒物质的抗压强度有着直接的影响, 围压越大, 颗粒物质抗压强度越好.

从图12同样可以看出, 当围压σc

≥50 kPa后, 随着围压的增大, 不仅应力峰值增大, 而且应力峰值对应的竖向应变值也在增大. 这说明当围压较大时, 颗粒物质有一个整体压缩的过程, 在此过程中颗粒物质不会发生局部变形. 当竖向应变达到其极限值后, 颗粒物质开始发生局部变形. 同时, 围压越大, 颗粒物质的极限应变值越大, 颗粒物质也越稳定. 当围压较小时, 如σc ≤15 kPa时, 应力峰值的出现不是很明显, 通过数值模拟结果发现, 这是由于围压较小

图12 不同围压下应力的变化关系图

时, 颗粒物质发生变形时主要以整体的水平扩散为主. 这也说明若不考虑颗粒

之间的黏性, 围压是颗粒物质内部是否会出现局部变形和剪切面的重要条件. 当围压较小时, 颗粒物质的变形是整体的, 其内部观测不到明显的局部变形和剪切面, 这个结论与Iwashita等人[10]的模拟结果是一致的. 同时, 通过数值结果发现, 随着围压的增大, 类固态的剪切模式逐渐减弱, 双剪切面模式逐渐向单个剪切面模式演化.

2.3.3 摩擦系数的影响

通过上面的力链图可以看出, 颗粒之间的摩擦力对颗粒物质剪切面的产生有直接的影响. 下面分析摩擦系数的大小对于颗粒物质变形特征的影响. 通过数值模拟结果发现, 当摩擦系数µ = 0时, 颗粒之间无摩擦力, 其内部产生了一些类似环形的粒子运动而没有出现明显的剪切面. 随着摩擦系数的增大, 颗粒物质内部逐渐出现了明显的剪切面, 如图13所示. 当µ 较小时, 利用切向速度的分布廓线分析可以发现, 颗粒物质内部主要是类流态的局部变形模式. 随着µ 的增大, 类固态的剪切模式逐渐出现并增强, 剪切面由一个逐渐演化为两个. Hu等人[12]在对密集规则排列的等粒径颗粒集合体的模拟中, 也得到了相似的模拟结果: 即随着摩擦系数的增大, 单个剪切面会逐渐演化成双剪切面形式. 由于双剪切面模式的抗压强度要高于单剪切面模式, 故摩擦系数越大, 颗粒物质的抗压强度越高, 整体稳定性也越好.

图14为不同摩擦系数下, 应力比随竖向应变的变化关系. 从图14中可以看出, 当摩擦系

701

王等明等: 颗粒物质局部变形的离散元模拟

图13 不同摩擦系数下颗粒物质在某一瞬时的速度矢量图

数较小时, 由于粒子之间的摩擦力较小, 其在空间的重新排列相对比较容易, 应力比随竖向应变有较大的波动. 在此过程中伴随着颗粒物质交替的应变强化和弱化, 颗粒物质抗压强度较差. 当摩擦系数较大时, 经过峰值应力以后, 应力比随着竖向应变波动较小, 颗粒物质内部一般不会出现较大的竖向位移, 颗粒物质抗压强度较好. 由于摩擦系数较小时, 颗粒物质内部只有单个的剪切面, 故图15仅仅给出了单个剪切面的倾斜角随摩擦系数的变化关系, 可以看出,

随着摩擦系数的增大, 剪切面的倾斜角也会逐渐增大. 这是由于当摩擦系数较大时, 在相同的围压下, 颗粒物质侧向运动的摩擦力变大导致的.

图14 摩擦系数对应力比变化的影响 图15 摩擦系数对剪切面倾斜角的影响

3 结论

本文在平面应变状态下, 对颗粒物质的局部变形特性进行了数值模拟, 同时对初始空隙比、围压、摩擦系数等因素对颗粒物质稳定性及剪切面的影响进行了分析, 得到以下结论.

(ⅰ) 平面应变状态下颗粒物质在竖向压缩下顶端会出现两个相互交叉的剪切面. (ⅱ) 相互交叉的两个剪切面沿着各自的法线方向具有不同形式的速度廓线分布, 分别表

702

中国科学 G辑: 物理学 力学 天文学 2008年 第38卷 第6期

现为类流体和类固体的剪切变形模式. 两种剪切变形模式具有一定的持续性, 并且在整体变形中类流体模式占主导地位.

(ⅲ) 初始空隙比对颗粒物质的变形有较大的影响. 初始空隙比越小, 在相同的围压下颗粒物质抗压强度越好, 并且随着空隙比增大, 剪切面的倾斜角逐渐增大.

(ⅳ) 围压会影响颗粒物质的抗压强度和变形模式. 围压越大, 应力峰值越大, 颗粒物质抗压强度越好. 随着围压的增大, 类固态的剪切模式逐渐减弱, 双剪切面模式逐渐向单个剪切面模式演化.

(ⅴ) 当摩擦系数µ 较小时, 颗粒物质内部主要是类流态的局部变形模式. 随着µ 的增大, 类固态的剪切模式逐渐出现并增强, 单剪切面模式逐渐演化为双剪切面模式, 颗粒物质的整体稳定性也逐渐增强.

参考文献

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

Tordesillas A, Walsh D C D. Incorporating rolling resistance and contact anisotropy in micromechanical models of granular media. Powder Technol, 2002, 124: 106—111

Kuhn M R. Structured deformation in granular materials. Mech Mater, 1999, 31: 407—429

Khalid A A, Stein S. Shear band formation in plane strain experiments of sand. J Geotech Geoenviron, 2000, 126(6): 495—503

Labuz J F, Dai S T. Residual strength and fracture energy from plain-strain testing. J Geotech Geoenviron, 2000, 126(10): 882—8

Wang Q, Lade P V. Shear banding in true triaxial tests and its effect on failure in sand. J Eng Mech, 2001, 127(8): 7—761

Nübel K, Huang W X. A study of localized deformation pattern in granular media. Comput Meth Appl Mech Eng, 2004, 193: 2719—2743

Tejchman J. Influence of a characteristic length on shear zone formation in hypoplasticity with different en-hancements. Comput Geotechnol, 2004, 31: 595—611

Gardiner B S, Tordesillas A. Micromechanics of shear bands. Int J Solids Struct, 2004, 41: 5885—5901

Iwashita K, Oda M. Rolling resistance at contacts in simulation of shear band development by DEM. J Eng Mech, 1998, 124(3): 285—292

Iwashita K, Oda M. Micro-deformation mechanics of shear banding process based on modified distinct element method. Powder Technol, 2000, 109: 192—205

Jiang M J, Yu H S, Harris D. A novel discrete model for granular material incorporating rolling resistance. Com-put Geotechnol, 2005, 32: 340—357

Hu N, Molinari J F. Shear bands in dense metallic granular materials. J Mech Phys Solids, 2004, 52: 499—531 Cundall P A, Strack O D L. A discrete numerical model for granular assemblies. Geotechnique, 1979, 29(1): 47—65

Zhu H P, Yu A B. A theoretical analysis of the force models in discrete element method. Powder Technol, 2006, 161: 122—129

Nezami E G, Hashash Y M. A fast contact detection algorithm for 3-D discrete element method. Comput Geo-technol, 2004, 31: 575—587

Hsiau S S, Yang S C. Numerical simulation of self-diffusion and mixing in a vibrated granular bed with the cohe-sive effect of liquid bridges. Chem Eng Sci, 2003, 58: 339—351

Bertrand F, Leclaire L A, Levecque G. DEM-based models for the mixing of granular materials. Chem Eng Sci, 2005, 60: 2517—2531

Aharanov E, Sparks D W. Rigidity phase transition in granular packing. Phys Rev E, 1999, 60: 60—66

Aharanov E, Sparks D W. Shear profiles and localization in simulations of granular materials. Phys Rev E, 2002, 65(5): 051302-1—12

703