不同细料含量土石混合料塑性行为离散元模拟1)

王涛 朱俊高,2) 刘斯宏

* (河海大学岩土力学与堤坝工程教育部重点实验室,南京 210098)

† (河海大学水利水电学院,南京 210098)

引言

土石混合料是指由大粒径块石和作为填充成分的细粒土组成的混合料[1-2].土石混合料在自然界中分布广泛,可作为基础填筑材料应用于水利、铁路、公路等工程建设之中.作为典型的二元混合颗粒材料,土石混合料的力学特性与细料含量(fc)密切相关.细料含量较低时,土石混合料由石料承担骨架,其力学特性主要由石颗粒间的摩擦特性决定;细料含量超过阈值细料含量时,土石混合料由土料承担骨架,石颗粒悬浮在土料骨架基质中,其力学特性主要由土颗粒间的摩擦特性决定[3].

目前,国内外学者针对土石混合料已开展了大量的室内试验[4-6]、现场试验[7-8]与数值模拟[9-11]研究,发现土石混合料的抗剪强度、剪胀特性、临界状态线与破坏模式等均与细料含量密切相关.然而,有关细料含量对土石混合料塑性行为的影响却鲜有研究,从细观层面阐明细料含量对土石混合料塑性行为的影响机制,对建立合理的弹塑性本构模型具有重要意义.

为了模拟土石混合料的应力应变关系,近年来众多唯象类本构模型相继被提出[12-13],这些本构模型是在总结试验规律基础上提出的,能够较好地拟合试验结果,但是缺乏物理机制的支撑.土石混合料是离散颗粒组成的集合体,其宏观应力应变由集合体内部接触的几何分布与接触力大小决定[14-15],但是在唯象类本构模型中其离散的本质并未被显性表达,在此层面上,离散单元法由于易获得代表体积单元中颗粒位置和接触信息,是研究细观参数如何影响宏观特征的有效方法.颗粒集合体由无数个细观结构单元组成,这些细观结构的几何与力学特性直接决定了集合体宏观的强度变形特性[16].近年来,针对细观结构的研究成果有效解释了颗粒材料的破坏模式[17]、失稳特性[18]与剪切带的形成过程[19].因此,探究加载过程中不同细料含量土石混合料内细观结构的演化规律可揭示细料含量对土石混合料塑性行为的影响机制.

本文目的在于探究细料含量对土石混合料塑性行为的影响规律及其细观机制.需要说明的是,细料含量对石料骨架结构与土料骨架结构塑性行为的影响机制不同,因此本文重点讨论石料骨架结构土石混合料.本文基于二维离散元数值模拟,研究细料含量对石料骨架结构土石混合料(fc=0,0.05,0.1)屈服面、塑性流动方向与失稳破坏等特性的影响规律,并从细观结构角度对数值模拟结果进行解释.

1 DEM 数值试验与分析方法

1.1 双轴压缩数值模拟

离散单元法最早是由Cundall 和Strack[20]提出,之后成功应用于岩土不连续介质的非线性力学分析中.离散单元法的基本思想是将不连续介质离散化为独立的单元,通过对各个独立结构之间的相互作用进行分析,进而研究整体的力学性质.离散单元法一般采用显式时步算法,假设单个时步内颗粒速度和加速度不变,求解颗粒的运动方程和力-位移方程,并在下一步中更新颗粒、墙体和接触的信息,如此不断更新交替,直至整个计算结束.

本研究采用开源软件Yade[21]进行计算,Yade主要应用于岩土工程领域,对土力学细观机理、山体滑坡、地震等工程问题有很好的支持,具有计算效率高的优点.

法向与切向接触力与接触位移可由下式确定

式中,kn为法向刚度,kt为切向刚度,δn为法向重叠量,δt为切向重叠量,ϕ为颗粒间摩擦角.

法向刚度kn可由接触颗粒的尺寸和材料模量E决定,并假定切向刚度与法向刚度的比值取定值r,kn与kt可由下式确定

式中,Rp和Rq分别为接触颗粒p和q的半径.

本次研究所采用的模型参数见表1,表1 的参数参照文献[17-19]进行选取,其中设置颗粒-墙体摩擦角为0°可确保墙体的应力为主应力(上下墙体为大主应力,左右墙体为小主应力).

表1 DEM 模型参数取值Table 1 Parameters for DEM tests

为确保数值试样为石料骨架结构,制备低细料含量试样(fc=0,0.05,0.1),采用间断级配,土颗粒在0.4～0.8 mm 粒径范围内均匀分布,石颗粒在4～8 mm 粒径范围内均匀分布.试样制备分三个步骤.(1)生成试样:首先在长方形墙体内生成相互无接触的颗粒群(见图1),将颗粒球心的z方向坐标固定以形成二维情况.(2)固结:首先采用内部加载法,通过增大颗粒粒径使墙体处达到90 kPa 的围压,随后停止增大颗粒粒径,利用墙体的移动来施加设计的围压(100 kPa 或200 kPa).(3)剪切:保持侧面墙体的围压不变,以足够小恒定的速率移动顶部和底部墙体以剪切试样,以保证剪切过程为准稳态[22].

图1 DEM 双轴压缩模型Fig.1 DEM biaxial model

本研究采用相对密实度指标,保持不同细料含量的试样初始相对密实度相同(Dr=10%),相对密实度的计算公式为

式中emax和emin分别为试样的最大与最小孔隙比.

在DEM 模拟中,通常通过设置固结过程中不同颗粒-颗粒摩擦角控制试样的密度.本研究将颗粒间摩擦角设为0°获得试样的最小孔隙比,这与文献[23-24]中所采用的方法相同.对于最大孔隙比,文献[25-26]研究发现,设置颗粒间摩擦角为35°可获得试样最大孔隙比,因此本研究通过设置颗粒间摩擦角为35°制备最大孔隙比试样.获取不同细粒含量试样的最大、最小孔隙比后,可在固结过程中不断试算颗粒间摩擦角直至其相对密实度达到10%.不同细料含量土石混合料试样的最大孔隙比、最小孔隙比与制样孔隙比列于表2.

表2 不同细料含量试样最大最小孔隙比、制样孔隙比与相对密实度Table 2 Minimum and maximum void ratios,void ratio for sample preparing and relative density for samples with different fines contents

图2 给出了fc=0 试样剪切过程中应力比η(η=q/p)与体积应变εv的变化曲线.可见,η与εv均随轴向应变增大而增大,随后达到某一恒定值,试样达到临界状态,体现出典型的松散试样的特性.加载初期,应力比急速增长,随后波动上升,在波动上升阶段,准稳态假设在某些加载步中可能不再成立.为了研究试样在不同应力状态下的稳定性,在双轴压缩过程中保存处于不同应力比的试样(如图2(a)中蓝色实心点),这些应力状态对应的η∈{0.15,0.22,0.33,0.43,0.46,0.49,0.52}.

图2 双轴压缩过程中应力比η 与体积应变εv 随轴向应变变化曲线(fc=0)Fig.2 Evoluitons of stress ratio η and volumetric strain εv with axial strain εyy (fc=0)

1.2 分岔失稳点分析

作为典型的非关联流动材料,岩土材料在达到塑性破坏极限之前便可能发生失稳破坏,通常采用二阶功失稳准则来描述失稳破坏.二阶功失稳准则最早由Hill[27]提出,并被广泛应用于岩土工程计算中.Darve等[28-29]详细描述了颗粒材料分岔失稳特性,并建立了颗粒材料失稳分析的框架.

二阶功失稳准则是指,若存在一个加载路径,在此应力路径下二阶功W2=dσ:dε 为负,在该应力路径下土体失稳.dσ 和 dε 分别为应力增量矢量和应变增量矢量.

利用二阶功失稳准则,采用经典的应力控制的方向加载法研究不同应力状态下试样的稳定性.方向加载分为两步:预稳定和方向加压.

(1)预稳定:1.1 节中试样是在加载过程中保存的,需要让其在恒定的应力下达到稳定,稳定的判别标准是试样的不平衡力Funb＜ 10-5.不平衡力是判断试样是否平衡的指标,其计算公式为

式中,Nc为颗粒Np的接触总数;Fcp为颗粒Np所受到的合外力,Fc为颗粒Np所受所有接触力的平均值.

1.1 节中保存的试样在预稳定过程中产生的体积应变如图2(b)所示,蓝色实心点为预稳定前的体积应变,蓝色空心点为预稳定后的体积应变.可见,当η ≤0.46 时,预稳定过程中的体积应变很小,可忽略;当 η=0.49,0.52,由于颗粒重排布预稳定过程中产生了不可忽略的体积应变.

(2)方向加压:预稳定结束后,施加不同方向的应力增量 dσ,如图3(a) 所示,应力增量的数值为,加载方向角为α.可根据式(5)计算施加在侧向墙和竖直墙上的应力 dσxx与 dσyy.当围压为100 kPa 时,dσ 值设置为5 kPa,当围压为200 kPa 时,由于在固结过程中试样储存了更多的弹性能,需要采用更大的 dσ 来释放这些弹性能,因此当围压为200 kPa 时,dσ 值设置为10 kPa

图3 应力增量dσ 及其应变响应增量dεFig.3 Stress probe dσ and incremental strain response dε

对于施加的应力增量dσ,可以得到相应的应变增量响应 dε,dε 的方向角为β,如图3(b)所示.可由dσ和 dε 计算归一化的二阶功[30]

图4为fc=0 的试样在不同初始应力比下的归一化二阶功,图4 中的红色圆圈代表为0,曲线在红圈内时,为负,反之,为正.可以看出对于所有的应力比,主要在第一和第三象限降低,且最小值出现在第三象限.图4(a)是η≤0.46 的情况,即预稳定过程中的体积应变可以忽略的情况,可以看出当应力比较小时,对于所有的加载方向均不存在负的,即试样在低应力比时是稳定的.随着应力比增大,最小值逐渐降低,当η=0.43,0.46 时开始出现负的,产生失稳锥(所有产生负的的加载方向的集合).图4(b)是η＞ 0.46 的情况,可见在η=0.49 时无负值,η=0.52 时存在负值,即随着η=0.46 增大至0.49 时,试样由不稳定状态转变为稳定状态,出现这一现象的原因是η=0.49 的试样在预稳定过程中产生了较大的体积变形,试样的细观结构因此产生了较大改变,即进行方向加载的试样与前期加载过程中保存的试样并不相同,因此研究失稳特性时,应重点关注预稳定过程中不产生较大体积变形的试样.

图4 在不同初始应力比下的归一化二阶功 (fc=0)Fig.4 Circular diagrams of the normalized second-order work at different stress ratio η (fc=0)

2 细料含量对土石混合料塑性行为的影响

2.1 失稳特性

图5 给出了不同细料含量的土石混合料在相同应力比下(η=0.43)的归一化二阶功.可以看出,随着细料含量增大,试样的最小值逐渐增大.当fc为0 时,存在一个较宽的失稳锥;当fc为0.05 时,失稳锥变窄;当fc增大至0.1 时,失稳锥消失,试样在所有加载方向上的均大于0.这说明对于石料骨架结构,增大细料含量有利于试样的稳定.

图5 不同细料含量的试样(fc=0,0.05,0.1)的归一化二阶功(η=0.43)Fig.5 Circular diagrams of the normalized second-order work computed for the specimen with different fines content (fc=0,0.05,0.1)under the same stress ratio

如前所述,岩土材料在达到莫尔库伦破坏极限之前就有可能发生破坏,即应力空间中出现分岔失稳区域[28-30].分岔失稳区域仅存在于非关联流动材料,对于关联流动材料,莫尔-库伦极限线与分岔失稳区域的下边界线重合,无分岔失稳区域.分岔区域越大,表示材料的非关联流动性越强.

为了进一步探究细料含量对土石混合料分岔失稳区域的影响,图6 给出了不同细料含量试样莫尔-库伦极限线与分岔失稳区域的下边界线的夹角θ,两条直线之间的区域即为分岔失稳区,可见随着细料含量增加,θ逐渐降低,说明细料含量的增加会降低材料的非关联流动特性.

关于细料含量增大会降低材料的非关联流动特性这一发现,将在2.2 节中重点讨论.

2.2 非关联流动性

在弹塑性力学框架中,应变增量可写成弹性应变增量(dεe)和塑性应变增量(dεp)之和

为了获得弹性应变,可在加载过程中设置颗粒间摩擦角为90°,此时颗粒间不会发生滑移,产生的应变全部为弹性应变.获得弹性应变增量之后,可根据式(7)计算塑性应变增量.

图7 给出了fc=0.1 的试样在应力比为0.43 时不同加载方向下的应变增量.可见弹性应变增量数据点形成了一个中心位于坐标原点的椭圆,表现出纯弹性行为[31].塑性应变增量数据点形成一条直线,表明塑性应变的方向与加载方向无关,证明了流动法则的存在.

图8 给出了不同应力加载方向角α下的塑性应变增量方向角βp及其数值大小,可见βp保持不变,根据流动法则,塑性应变方向与加载方向无关且与该应力状态下的塑性势面法向一致,因此βp即为塑性势面法方向.在弹塑性力学框架中,当应力在屈服面以内变化,只产生弹性应变,应力若超过屈服面,将产生塑性应变.当α ∈[60,240]时,试样产生塑性应变,可知α=60和α=240为屈服面法方向的两个切方向,因此屈服面的法方向为=150.

图9 比较了不同细料含量试样塑性应变大小,可见在弹性区(0≤α≤60与240≤α≤360),不产生塑性变形;在塑性区(60≤α≤240),产生塑性变形.由于W2=dσ:dε=dσ:dεe+dσ:dεp,且dσ:dεe≥0,因此塑性应变很大程度上决定了二阶功的正负.由图9 可见增大细料含量可以减小塑性应变大小,这也解释了为什么增大细料含量有利于试样的稳定.

图9 不同细料含量试样(fc=0,0.05,0.1)在相同应力比下(η=0.43)的塑性应变大小Fig.9 Norms of plastic strain for specimens with fc=0,0.05,0.1 at the same η

图10 给出了屈服面f与塑性势面g的示意图,及其法方向n和m随细料含量的变化情况.可见,随着细料含量增加,屈服面法方向n保持不变,而塑性势面法方向m逐渐靠近n,m与n之间夹角减小,试样的非关联流动特性降低,这与2.1 节中的结果一致.图11 为有细粒和无细粒情况下塑性变形流动方向的示意图.集合体变形主要是来自力链在载荷作用下崩塌引起的颗粒重新排布,图11(a)所示,初始状态下,力链的方向与加载方向(y方向)一致,一旦力链发生崩塌竖向会产生较大的竖向变形,同时水平向会产生一定的膨胀,需要注意的是竖直方向的变形远大于水平向的变形.当试样中含有细颗粒时,如图11(b)所示,由于细颗粒会侧向支撑力链,因此试样竖向的变形相较于无细颗粒时更小,所以加入细颗粒后会使塑性流动方向发生逆时针的旋转,这与图10 所展示的m方向的变化一致.

图10 屈服面f 与塑性势面g 的示意图及其法方向随细料含量的变化情况(η=0.43)Fig.10 Schematic diagram of yield surface f and plastic potential surface g (on the left) and their normal directions m and n (on the right)for specimens with different fc at the same η=0.43

图11 有细粒和无细粒情况下流动方向示意图Fig.11 Illustration of the plastic flow direction with and without fine grain

3 细观机理

3.1 基本概念回顾

(1)各向异性

各向异性是颗粒材料的一个重要特性[32],可分为几何各向异性和力学各向异性[33].几何各向异性可用标量ac定量描述

(2)力链

力链是颗粒集合体最基本的细观结构,力链在颗粒结合体中负责传递较大的接触力并具有准直线的几何特征.Peters等[36]对产生力链的三个条件归纳如下:(1) 力链上的颗粒的最大主应力大于所有颗粒的平均最大主应力;(2) 力链颗粒与其相邻颗粒的圆心连线和最大主应力方向夹角小于45°;(3) 一条力链至少含有3 个颗粒.

力链的屈曲会使集合体承载力降低,甚至导致试样发生崩塌[37].如图12 所示,考虑3 个颗粒组成的力链,力链的偏角θ定义为中心颗粒与两端颗粒圆心连线的夹角.假设在t时刻力链偏角为 θt,经过时间 dt后,力链偏角为θt+dt,则力链屈曲角θb为θb=θt-θt+dt.需要定义一个阈值θc,当θb＜θc时,认为力链未发生屈曲,反之,力链发生屈曲[37].本研究中,作者经过反复试算,发现阈值取为1 度时试验规律较好,因此选取θc=1°.

图12 力链偏角示意图Fig.12 Illustration of force chain bending angle

3.2 细料含量对内部几何结构与力学状态的影响

为探究细料含量对试样几何与力学各向异性的影响,图13 绘制了不同细料含量试样在双轴压缩过程中ac(几何各向异性)和an(力学各向异性)随应力比的变化过程.可见,不同细料含量试样在相同应力比时an相同,但ac不同,这意味着在相同应力比时不同细粒含量试样具有相似的力学状态,但内部几何结构不同.

图13 不同细料含量试样 ac (几何各向异性)和 an (力学各向异性)随应力比的变化过程Fig.13 Evolution of ac (filled symbols) and an (empty symbols) with respect to η in specimens with different fc during biaxial tests

在外载荷作用下,试样的塑性变形主要取决于颗粒间的滑移、滚动、原有接触的打开与新接触的生成,这些都与试样内部的几何结构相关,而不同细料含量试样具有不同的内部几何结构,因此细料含量会影响试样的塑性流动方向(塑性势面法向).这一结果与一些学者发现的剪胀与试样内部几何结构密切相关这一结论一致[38-39].然而,屈服面的方向与集合体的应力状态相关,随着细料含量增加,集合体力学各向异性不发生改变,因此细料含量不改变屈服面的方向.

材料的屈服对应着承载力的降低,细观上表现为集合体内部的力链难以抵抗外载荷,发生屈曲变形,并伴随着大量的接触滑动.以下从力链屈曲、接触滑动角度进一步解释为何细料含量不改变屈服面方向.图14 为不同细料含量试样在双轴压缩过程中力链屈曲角的概率密度函数.可见,不同细料含量的试样屈曲角的概率密度函数几乎重合,这说明试样的内部力学状态基本相同,不受细料含量的影响.

图14 不同细料含量试样(fc=0,0.05,0.1)在双轴压缩过程中力链屈曲角的概率密度函数(应力比为0.43)Fig.14 Buckling angle probability density functions of specimens with fc=0,0.05,0.1(at stress ratio η=0.43)

为了进一步探究细料含量对试样力学状态的影响规律,定义接触滑动因子指标Ip

式中,ϕ 为颗粒-颗粒间摩擦角;Fn和Ft分别为法向和切向接触力.Ip介于0 至1 之间,反映了某一接触是否接近滑动,Ip=1 说明接触滑动.

图15 为不同细料含量试样的所有接触数Ip概率密度函数,可见,不同细料含量的试样Ip具有相似的概率密度函数,这一结果进一步佐证了细颗粒的加入并不会很大程度影响颗粒集合体的内部力学状态,因此,细颗粒的加入不影响屈服面的方向.

图15 不同细料含量试样的所有接触数 Ip 概率密度函数(应力比为0.43)Fig.15 Sliding index probability density functions of specimens with fc=0,0.05,0.1 at the same stress ratio η=0.43

4 结论与展望

通过二维离散元数值模拟,研究了细料含量对石料骨架结构整体稳定性和塑性行为的影响规律,主要结论如下.

(1)细颗粒有利于颗粒集合体的整体稳定,细颗粒能够限制集合体的塑性变形,提高试样的稳定性.

(2)细颗粒控制颗粒集合体塑性变形的方向,随着细料含量增加,塑性势面法方向和屈服面法方向之间的夹角减小,材料的非关联流动性降低,材料分岔失稳区域变窄.

(3)尽管加入到石颗粒中的部分细颗粒与石颗粒共同承担骨架作用,但是细颗粒的加入不影响颗粒集合体的力学状态,不改变试样的屈服面方向.

细观分析结果表明不同细料含量的试样在经历相同加载路径至同一应力比时,具有相似的力学状态和不同的内部几何结构.从实用性角度来说,这一发现对构建考虑细料含量影响的土石混合料弹塑性本构模型具有指导意义,例如对于石料骨架结构土石混合料,屈服函数应与细粒含量无关,塑性势函数应为细粒含量相关,且保证材料非关联流动性随细粒含量增大而减弱.

本文主要研究目的在于揭示细粒含量影响土石混合料塑性行为的细观机制,模拟时仅采用了常规二维球形颗粒,拟在后续研究中拓展至三维情况并考虑颗粒形状的影响