椭球颗粒体系剪切过程中自由体积的分布与演化1)

2021-11-10邹宇雄李易奥邱焕峰

力学学报 2021年9期

关键词：元胞椭球剪切

邹宇雄马刚,2) 李易奥王頔邱焕峰周伟

* (武汉大学水资源与水电工程科学国家重点实验室,武汉 430072)

† (水工岩石力学教育部重点实验室,武汉 430072)

** (中国电建集团贵阳勘测设计研究院有限公司,贵阳 550081)

引言

颗粒物质是由大量离散固体相互作用而形成的复杂体系,在自然界和工业界广泛存在.以水利和岩土工程中的堆石体为例,由于取材方便,且对地形、地质条件有较强的适应性,在堆石坝、路堤、机场高填方地基等工程建设中得到广泛应用.这些工程一旦出现问题和隐患,将会严重威胁人民生命财产安全,因此迫切需要深入研究颗粒材料的物理力学特性.

由于颗粒材料具有结构异质性、各向异性和多尺度结构等特点,其表现出复杂的物理力学特性,如剪胀[1]、应变局部化[2],固−液转变[3-4]等.对岩土类颗粒材料,自1963 年Roscoe 等[5]提出剑桥模型以来,许多学者根据颗粒材料的宏观应力变形特性发展了各种各样的本构理论[6-8].但目前还没有统一的理论能够解释颗粒材料所有的应力变形特性.蒋明镜[9]指出这些基于唯象的土力学本构理论虽然在岩土工程实践中发挥了重要作用,但是难以从机理上刻画岩土颗粒系统的非连续性、大变形、破坏等问题.因此,越来越多学者尝试从细观尺度研究颗粒材料复杂力学行为的起源,并试图将颗粒形状、表面摩擦、粒径分布等颗粒尺度的特性与宏观力学响应建立联系.

岩土颗粒材料的形状各异,大量研究表明颗粒形状对颗粒材料的剪切强度[10-13]、塑性行为[14]以及堆积特性[15-16]等有显著影响.深入研究颗粒形状对颗粒材料物理力学特性的影响,剖析这些影响的细观机理,对颗粒材料的理论研究和工程应用有重要意义.受试验技术限制,目前只能通过光弹试验技术[17]、X-ray 断层扫描技术[18]等研究颗粒材料宏观力学响应的细观机理.与此同时,细观数值模拟方法在颗粒材料的宏细观力学特性研究方面得到广泛应用,如离散单元法(discrete element method,DEM)[19]和连续离散耦合分析方法(combined finite and discrete element method,FDEM)[20-21].其中FDEM 中将颗粒离散为有限元网格,理论上可以考虑任意颗粒形状,在模拟复杂颗粒形状方面具有显著优势.

虽然颗粒材料的结构特性与其物理力学性质的相关性已经得到大量研究证明[22-26],但是具有复杂颗粒形状和宽粒径分布的颗粒材料的结构表征,以及建立宏观力学响应和细观结构的联系仍是一个重大挑战.基于该现状,颗粒统计力学得到了快速发展,例如Edwards 等[27-28]建立了颗粒系统的体积系综理论(Edwards' volume ensemble),用统计物理的框架来研究非平衡态的颗粒堆积问题.该理论强调了局部自由体积(孔隙)在描述静态颗粒材料阻塞行为中的重要作用.在工程中,常采用孔隙率、孔隙比、体积分数等表征颗粒材料的密实程度,但是这些表征量仅能提供宏观尺度的自由体积信息,忽视了颗粒尺度自由体积的形态、分布和演化.

已有研究表明,颗粒尺度的自由体积与颗粒材料的力学性能、变形特征等密切相关[29-30],所以有必要探索颗粒材料的自由体积的统计分布特性及其在剪切过程中的演化规律.将颗粒与其周围的自由体积组成一个局部元胞(以Voronoi 元胞较为典型)用以定量分析颗粒材料自由体积特征是一种常见思路.不少学者研究了颗粒材料在特定状态下的局部自由体积[31-32].但是颗粒材料在剪切过程中会产生颗粒重排列现象,颗粒材料的细观结构和自由体积会不断演化.因此,探索剪切过程非球颗粒材料的局部自由体积的演变规律,可以为研究颗粒材料的宏细观多尺度力学特性提供新的思路.

本文将选择具有不同非球度的椭球颗粒体系为研究对象,采用连续离散耦合分析方法进行三轴剪切数值模拟,分析不同非球度的椭球颗粒试样在剪切过程中的宏细观力学响应,并基于Voronoi 元胞对椭球颗粒体系的局部自由体积进行定量分析,探讨局部自由体积的分布特征和演化规律以及受颗粒形状的影响.

1 椭球颗粒体系三轴剪切数值试验

1.1 FDEM 基本理论

FDEM 最早由Munjiza 等[20]提出,充分融合了有限元和离散元方法的优势.颗粒的应力变形计算采用有限单元法,颗粒运动和接触模型采用离散单元法.FDEM 中颗粒均由有限元网格构成,相比其他不连续数值模拟方法而言,最吸引人的优点是能够考虑各种复杂形状[33-36].

有限元网格不仅定义了颗粒的形状和边界,而且可以很方便进行接触检索以及引入接触模型.在FDEM 中,罚函数法被用于颗粒间的接触分析,假定接触力偶相互侵入,根据重叠量在节点间产生分布式接触力.切向力的最大值根据库仑摩擦定律控制.本文采用线性接触模型计算颗粒间的接触力,基于两个接触颗粒的相对速度,通过阻尼力定律确定阻尼力.接触力的计算公式如下

式中,和分别为节点法向与切向接触刚度,定义knkt为kn=pnAnode,kt=ptAnode,pn和pt分别为法向和切向罚刚度,Anode为接触面上的节点控制区的面积,为滑移率,Δt为时间增量,δn为节点侵入量,βn和βt为法向与切向临界阻尼系数,ζn和 ζt为法向和切向恢复系数,m为节点质量.

1.2 颗粒形状与数值试验

Barrett[37]将颗粒形状量化分为形态、圆度和表面纹理3 个层面.本文重点关注形态,故选取椭球颗粒体系作为研究对象.共生成7 种主轴长度不同的椭球颗粒,颗粒的各种形态参数如图1 所示.其中L,I,S分别代表椭球长轴、中轴和短轴的长度.α 为Domokos系数[38],定义为α=(1/S+1/I+1/L)表征颗粒主轴各向异性,在椭球形颗粒中可以反映颗粒的非球度.ψ为球度,定义为表征颗粒形状偏离圆球的程度,式中A和V分别代表颗粒的表面积和体积.后文分析中选取 α 作为椭球颗粒的形状量化指标.

图1 颗粒形状描述Fig.1 Particle shape descriptors

在立方体容器中生成等效半径(定义为等体积球体的半径)在6.8～ 10.2 mm 范围内均匀分布的6394 个颗粒,颗粒位置和倾向均随机以避免产生初始各向异性.采用各向同性压缩制样,制样过程中颗粒间摩擦系数设为0.1,六面无摩擦的刚性墙体缓慢压缩试样直至达到目标围压0.5 MPa.图2 为α=3.394的数值试样,试样尺寸为0.42 m × 0.42 m ×0.42 m.采用常规三轴应力路径加载,围压保持σ2=σ3=0.5 MPa,控制顶部和底部两块刚性墙体以恒定速度相向移动.

图2 数值试样Fig.2 Numerical sample

1.3 细观参数率定

数值试验的可靠度很大程度取决于参数取值.本试验中颗粒密度、弹性模量和泊松比参考天然砂砾石的自然属性.根据岩石材料的回弹性质[39],临界阻尼系数取0.03,对应的恢复系数为0.9.Tatone 和Grasselli[40]验证了当接触罚刚度设为10 倍弹性模量时所得到的弹性响应与室内物理试验匹配较好,故本试验接触罚刚度设为弹性模量的10 倍.通过颗粒柱坍塌物理试验和数值试验对比分析来率定颗粒摩擦系数.具体步骤如图3 所示:选取粒径在6～ 10 mm均匀分布的砂砾石装入无底板圆筒,随后缓慢抽离圆筒,重复10 次试验,测得休止角均值为29.19°.对物理试验中的砂砾石进行扫描,通过主成分分析得到主轴长度,在此基础上重构出形状相近的椭球颗粒.调整颗粒摩擦系数进行颗粒柱坍塌数值试验,使得休止角逼近物理试验的结果,最终摩擦系数设为0.5 时休止角为29.5°.FDEM 数值模拟所需细观参数汇总于表1.

表1 FDEM 数值模拟细观参数Table 1 Parameters used in the FDEM simulation

图3 颗粒柱坍塌试验Fig.3 Granular column collapse tests

1.4 宏观力学响应

图4 和图5 分别为7 组不同非球度的椭球颗粒试样的偏应力−轴向应变和体积应变−轴向应变关系曲线.不同试样的宏观响应呈现相似的规律,偏应力在加载初期迅速上升达到峰值,随后开始下降,呈现出明显的应变软化现象.试样在加载初期有少量体积压缩,随后发生明显的体胀.在较大轴向应变时偏应力和体积应变都趋于稳定或呈缓慢变化趋势.这种应力和体积响应在密实的无黏性颗粒材料中非常典型[6].随形状参数 α 增大,即颗粒非球度增大,试样的剪切强度明显增大,并伴随着更加显著的剪胀.目前普遍认为这种影响源自非球颗粒的互锁效应,其增强了颗粒抵抗转动的能力[12,41].

图4 不同非球度的椭球颗粒偏应力−轴向应变Fig.4 Curves of deviatoric stress versus axial strain for ellipsoidal particles with different asphericity

图5 不同非球度的椭球颗粒体积应变−轴向应变Fig.5 Curves of volumetric strain versus axial strain for ellipsoidal particles with different asphericity

2 非球形颗粒的Voronoi 剖分方法

Voronoi 剖分作为一种基本几何结构划分方法,将离散的颗粒在空间上划分,利用所得到的空间结构来表征材料的细观结构,目前已成为结构表征的一种常用手段.其中圆球颗粒的Voronoi 剖分被广泛应用,计算软件也较为成熟(如Voro++[42]).相比而言,非球颗粒的Voronoi 剖分由于实现困难,在颗粒材料的计算分析中鲜有报道.为实现非球颗粒的Voronoi 剖分,Schaller 等[43]提出了Set Voronoi 剖分方法,现已被成功运用于非球颗粒体系的细观结构分析[44-45].具体步骤如图6 所示:(1)在颗粒表面均匀分布足够多的点,得到每个颗粒表面的离散点集;(2)计算所有颗粒表面离散点的Voronoi 元胞;(3)将属于同一颗粒的Voronoi 元胞合并,形成该颗粒的Voronoi 元胞.

图6 Set Voronoi 剖分二维示意图Fig.6 Two-dimensional illustration of Set Voronoi tessellation

上述步骤适用于完全分离的非球颗粒体系.然而在FDEM 中,由于颗粒间会相互侵入,会造成局部计算失准.为了避免这种情况,在进行离散点的Voronoi 剖分之前,将每个颗粒表面离散点沿其法向方向收缩一定距离,可在不影响Voronoi 剖分结果的同时避免局部误差[43].图7 为圆球颗粒和椭球颗粒的Voronoi 元胞示意图,可以看出圆球的Voronoi元胞是凸多面体,而椭球的Voronoi 元胞表面会存在凹面,形态更加复杂.

图7 圆球和椭球颗粒的Voronoi 元胞Fig.7 Voronoi cells of spherical and ellipsoidal particles

3 Voronoi 元胞的各向异性

对于颗粒材料而言,即使是圆球颗粒体系,在剪切过程中仍然会产生各向异性[46].由于各向异性与颗粒材料的力学性质密切相关,相关研究受到广泛关注.例如Ouadfel 和Rothenburg[24]基于组构和接触力的各向异性张量,推导了(stress-force-fabric,SFF)关系,建立了组构各向异性、接触力各向异性与宏观应力比之间的关系.组构和接触力是对接触状态的描述,相比而言,Voronoi 元胞是对颗粒局部空间的描述,其形状的不规则程度可以用来量化局部空间各向异性.由于Voronoi 元胞本质上是一个多面体,在量化其形状特性上可以借鉴于颗粒形状的描述指标.引入最常见的球度,描述Voronoi 元胞形状偏离圆球(各向同性)的程度,定义其为与元胞等体积球体的表面积与元胞表面积的比值

式中,Vc和Sc分别代表Voronoi 元胞的体积与表面积.

为了避免边界效应,本文分析排除了在边界2 倍平均粒径范围内的颗粒.图8 为不同非球度的椭球颗粒试样在不同加载阶段(εa为轴向应变)的Voronoi 元胞球度 ψc的概率密度分布函数(probability density function,PDF).可以发现 ψc呈现高斯分布,图中实线为高斯分布函数拟合曲线.在加载过程中 ψc分布逐渐向左移动(均值降低),且在加载到较大应变时趋于稳定.这说明颗粒材料在剪切过程中,局部空间的各向异性会逐渐增强直至系统达到临界状态,且随颗粒非球度增大各向异性增强越明显,这与组构各向异性和接触力各向异性的演化类似[41,47].在剪切过程中,ψc分布的标准差也逐渐增大.

图8 不同非球度的椭球颗粒在不同加载阶段的Voronoi 元胞球度概率密度分布Fig.8 Probability density distributions of sphericity of Voronoi cells at the different loading states for ellipsoidal particles with different asphericity

当颗粒材料充分剪切后,会达到一个与初始状态无关的力学稳定态,即临界状态.颗粒材料在该状态时的力学性质受到广泛关注[7,11].本文试样虽未完全达到理想的临界状态,但是在加载后期应力状态和体积已经趋近于稳定,故选取轴向应变40%时为试样的“临界状态”,代表充分剪切后的状态.图9 为不同非球度的椭球颗粒试样在临界状态的Voronoi元胞球度概率密度分布.随形状参数 α 增大 ψc分布向左平移,分布更广.图10 展示了临界状态下Voronoi元胞球度与颗粒形状参数 α 的线性关系.这说明局部各向异性的分布强烈依赖颗粒形状.相比圆球,椭球颗粒可以产生更为复杂的局部结构.

图9 不同非球度的椭球颗粒在临界状态的Voronoi 元胞球度概率密度分布Fig.9 Probability density distributions of sphericity of Voronoi cells at critical state for ellipsoidal particles with different asphericity

4 局部自由体积分布及演化规律

为了量化颗粒尺度自由体积的大小,引入局部孔隙比eloc=(Vc−Vp)/Vp,其中Vp和Vc分别为颗粒体积和其对应的Voronoi 元胞体积.图8 展示了不同非球度的椭球颗粒试样在剪切过程中局部孔隙比eloc的概率密度分布.局部孔隙比eloc呈现类高斯分布.针对局部自由体积的统计分布特性已有一定报道,Schaller 等[45]发现扁椭球颗粒局部体积分数1/(1+eloc)近似服从高斯分布,Dong 等[32]发现椭球颗粒和柱状颗粒局部孔隙比服从对数正态分布,Zhao 等[30]认为局部孔隙率eloc/(1+eloc)服从一个修改的对数正态分布.此外,还有一些研究发现圆球颗粒堆积体的局部体积分数的倒数(1+eloc) 服从k−Γ 分布[48-49].上述研究中以k−Γ 分布的拟合参数最少,其函数为

式中,k为表示分布函数形状的参数,Γ 为伽马函数,emin为当前堆积状态下最小局部孔隙比,loc为局部孔隙比的均值.

图11 中实线为k−Γ 函数拟合曲线.虽然k−Γ 分布是针对圆球颗粒体系提出的,但是对于椭球颗粒体系仍然适用.在加载过程中eloc分布逐渐向右移动,峰值减小并且分布范围更宽,这意味着颗粒在剪切过程中局部自由体积变大且分布更广.在加载后期这种演变逐渐趋于稳定,与试样先剪胀后达到临界状态的现象相吻合.随非球度增大这种演变越明显,这与图5 中剪胀随 α 增强相对应.

图11 不同非球度的椭球颗粒在不同加载阶段的局部孔隙比概率密度分布Fig.11 Probability density distributions of local void ratio at the different loading states for ellipsoidal particles with different asphericity

图12 为全局孔隙比eglo=0.638时不同非球度的椭球颗粒试样的局部孔隙比的概率密度分布,此时7 组试样剪切状态不同,轴向应变分别为9.0%,12.7%,15.2%,21.0%,27.1%,31.8%和37.7%.可以发现不同非球度试样的局部孔隙比分布基本一致,这说明颗粒局部孔隙比的分布仅与全局孔隙比相关,不受颗粒形态和剪切状态的影响.图13 为不同非球度的椭球颗粒试样的局部孔隙比分布标准差σ(eloc) 与全局孔隙比eglo的关系.σ (eloc)与eglo呈现明显的线性关系,且这个线性关系不受颗粒形态影响,进一步证明了颗粒局部孔隙比的分布仅受全局孔隙比控制.

图12 不同非球度的椭球颗粒在全局孔隙比相同时局部孔隙比概率密度分布Fig.12 Probability density distributions of local void ratio at the states as same global void ratio for ellipsoidal particles with different asphericity

图13 局部孔隙比的标准差与全局孔隙比的关系Fig.13 Relationship between standard deviation of local void ratio and global void ratio

从局部孔隙比分布的演化可知剪切过程中局部自由体积的变化是非均匀的,探索局部孔隙比波动可以帮助了解颗粒体系剪胀的细观机理.选取了5 个加载阶段(εa=4%,10%,20%,30%和40%),分析不同非球度的椭球颗粒试样在应变窗口 Δ εa≈0.8%的局部孔隙比波动.图14 为不同非球度的椭球颗粒试样在不同加载阶段的局部孔隙比波动的概率密度分布.局部孔隙比的波动呈现非对称拉普拉斯分布(asymmetric laplace distribution,ALD),其函数为

图14 不同非球度的椭球颗粒试样在不同加载阶段的局部孔隙比波动概率密度分布Fig.14 Probability density distributions of local void ratio fluctuations at the different loading states for particles with different shapes

式中,m为分界点,λ 为尺度参数,κ 为非对称参数.

Δeloc的分布可以分为左右两个非对称的指数函数(在单对数坐标系下表现为斜率不一致的两段线性分布),分界点在 Δeloc=0附近.左右两侧分布的非对称性刻画了局部自由体积收缩和膨胀的博弈.在加载初期(εa=4%)体积膨胀明显处于优势地位(图中图14(a)斜率绝对值小于图14(b)),试样表现出明显的剪胀行为.随形状参数 α 增大,两侧斜率差异增大,这意味着有更显著的体积膨胀,与宏观体积响应相符合(图5).随着剪切进行,代表体积膨胀的图14(a)体积波动分布变化较小,而图14(b)分布的斜率有明显下降,两侧趋于对称,即剪胀会随剪切过程逐渐停止.

ALD 中非对称参数 κ 为分布两侧斜率比值的平方,描述了ALD 的非对称程度,也刻画了体积膨胀的强度.κ 越小体积膨胀越明显,κ=1时体积无变化.如图15 所示,κ 与全局孔隙比eglo呈现明显的线性关系,并在 κ=1附近截断.这意味着体积膨胀的强度很大程度上取决于试样当前全局孔隙与临界全局孔隙比的差值,故随 α 增大体胀更显著可以归因于椭球颗粒在相同固结条件下可以形成更密实的堆积.