基于线性接触模型的小麦三轴试验细观模拟

曾长女，于 航

（河南工业大学 土木建筑学院，河南 郑州 450001）

0 引言

在筒仓设计中，筒仓卸料引起的侧压力规律是关键要素之一.仓内储料的特性对筒仓侧压力具有较大的影响，其宏观特性取决于储料颗粒的细观特性.利用数值试验，可方便地获得需要的储料特性计算参数.离散元法用来处理散粒体流动特性具有很大的优势.PFC3D是ITASCA 公司基于Cundall 等[1]提出的离散元方法，是模拟不连续物体运动及相互作用的三维颗粒离散元模拟程序.该程序被广泛地应用于模拟大变形、散体流动特性等研究中.

土工三轴试验的离散元模拟开展得较多.罗勇等[2]利用PFC3D对无黏性土的室内常规三轴试验及其剪切带形成和发展进行了数值模拟.耿丽等[3]以PFC3D为工具，提出了考虑围压影响的粗粒土三轴试验细观模型，并研究了不同的颗粒细观参数对粗粒土强度的影响.邵磊等[4]通过生成不规则形状的数值颗粒，模拟三轴排水剪切试验.李识博等[5]利用PFC3D对陇西地区黄土三轴固结不排水剪切试验进行模拟，考虑了不同围压、不同含水率对应力、应变等的影响.Cheng 等[6]运用PFC3D对岩石的团块组合进行模拟，研究了三轴应力-应变关系.Hong-hua Zhao 等[7]采用PFC3D中的线性接触模型，研究了摩擦系数对三轴试验获得的密砂应力-应变的影响.Backstrom 等[8]利用颗粒离散元研究了花岗岩细观参数及其应力-应变曲线.现有的研究多集中于不同土体三轴试验，基于粮食散体细观参数的研究很少，而利用离散元法分析筒仓卸料时，细观参数多直接设定，没有更深入地说明.因此有必要获得小麦等粮食散体的细观参数，为今后系统地研究粮食散体的细观行为提供基础参数.

小麦的宏观特性必然取决于其颗粒的细观特性，离散元法从细观角度对小麦的力学特性进行数值模拟，分析小麦的细观参数变化对小麦的宏观力学行为的影响，再现小麦剪切试验曲线.在进行数值计算和理论研究时，小麦堆体的弹性模量等参数的确定对筒仓结构设计具有重要意义.作者基于小麦三轴试验，通过PFC3D模拟室内小麦三轴压缩试验，经反复试算，标定细观参数，并反算获得其宏观参数.

1 PFC3D 模拟过程

1.1 细观计算模型的选取

PFC3D中的接触模型包括3 个模型：接触刚度模型；滑移-分离模型；黏结模型.接触刚度模型描述了接触力和相对位移的弹性关系如图1 所示，采用线性接触刚度模型，该模型应用较多，细观参数较少，也可满足模拟结果的精度要求.

在线性接触刚度模型中，接触刚度为不变值，应力和相对位移是线性相关的.法向作用力满足式（1）：

式中：Kn为法向刚度，属割线模量；Un为法向位移.

切向作用力满足式（2）：

式中：Ks为切向刚度；ΔUs为切向位移增量.

图1 线性接触刚度模型Fig.1 Linear contact stistness model

在细观模拟中，有颗粒法向刚度Kn、颗粒切向刚度Ks、颗粒间摩擦系数μ、三轴试验加载速度v等细观参数需重点考虑.

1.2 三轴细观模型的建立

1.2.1 生成孔隙率为n 的试样

通过颗粒半径放大法使模型试样达到目标孔隙率.对所有颗粒采用相同的半径放大系数m，即：

式中：n0为初始孔隙率；n 为目标孔隙率；R0为初始颗粒半径.

模型试样生成时，先假定半径放大系数m，通过式（3）反算出初始孔隙率：

然后设定颗粒半径分布范围是（Rmin，Rmax），颗粒平均半径为：

试样尺寸为h×d，则生成的颗粒数目为：

根据前面步骤反复计算，直至达到设定的孔隙率.

1.2.2 墙体刚度的选取

生成的数值模型如图2 所示，包括圆柱形墙体和上、下两面墙体.圆柱形墙体模拟三轴试验中的橡皮膜，上、下两面墙体模拟三轴试验中的加载墙体.墙体刚度的选择对试验结果影响较大.已有研究表明，墙体刚度过大会导致不容易平衡，太小则导致侧压力过小，从而与实际不符.在数值模拟中，可采用圆柱形墙体的刚度为颗粒刚度0.1 倍时，可很好地模拟橡皮膜的柔性约束，上、下墙体的刚度为颗粒刚度10 倍以上时，模拟刚性加载面效果更好.因此，作者选择圆柱形墙体的刚度为颗粒刚度的0.1 倍，上、下加载墙体的刚度为颗粒刚度的10 倍[9].

图2 三轴试验试样数值模型Fig.2 Numerical model of triaxial test sample

1.2.3 三轴试验加载模拟

三轴试验模型加载时，圆柱形墙体用来模拟橡皮膜，可施加力模拟三轴试验中的围压.通过在上、下两墙施加力模拟三轴试验的轴向压力.具体计算公式如式（8）、（9）、（10）.

式中：r 为试样底面半径.

式中：h 为试样高度.

轴向应变和径向应变都可用下式计算：

式中：L 为当前式样的高度或宽度；L0为初始试样的高度或宽度.

试验加载步骤为：

（1）通过伺服机制给圆柱墙一个速度控制围压不变.

（2）给上、下墙大小相同、方向相反的速度压缩试样：

式（13）、（14），以使模型达到目标的围压和轴向荷载.

（4）计算参数G 的选择.

在一个计算时步中，墙体应力的增量为：

式中：Nc为 颗粒与墙体的接数目；为所有和墙体接触的平均刚度；A 为墙体的面积.

求解过程需要一定的稳定性，为确保加载过程的稳定，在实际操作中增加一个松弛因子α（一般为0.5）.

由式（11）、（12）、（15）和（16）可得：

为了获得稳定条件，把参数G 定义为：

在每一次循环开始时，墙的速度给定一个值满足式（11）、（12），而G 要满足式（18）.

2 细观模拟结果

PFC3D三轴模型参数的筛选与输入，是一个非常费时的过程，生成的颗粒数量对计算速度影响很大，经过多次试算，本次试验模拟采用的试样尺寸取值0.30 m×0.15 m，颗粒半径取值范围定为（0.004，0.008）m，孔隙率为0.3，该模型生成3 877个颗粒.

2.1 应力-应变曲线

根据已获得的室内试验的应力-应变曲线，不断调整细观参数并反复试算，模拟得出的应力-应变曲线与真实试验得出的应力-应变曲线对比如图3 所示.由图3 可知，试验曲线与数值模拟曲线拟合较好，由此可获得合适的细观参数，具体如表1 所示.

图3 试验与数值模拟的应力-应变曲线Fig.3 Stress-strain curve from test and numerital simulation

表1 材料参数Table 1 Material parameters used in simulations

2.2 三轴试验中小麦宏观参数推算

（1）弹性模量E

E 的确定如图4 所示，可计算获得小麦的宏观参数.对于三轴试验而言，采用式（19）可获得其弹性模量.初始阶段应力-应变曲线的斜率可取为弹性模量：

图4 E 的确定Fig.4 Stress-strain curve from simulation

（2）泊松比ν

由试验确定泊松比比较困难，通过如图5 所示体变-轴向应变曲线可确定泊松比.

（3）内摩擦角φ

该模拟试验时围压σ3保持不变，偏应力最大值σ1max，由图4 可知，利用公式（20）可计算得到颗粒的内摩擦角：

图5 泊松比确定Fig.5 Poisson ratio

（4）剪胀角ψ

体变-轴向应变曲线与x 轴有一交点，利用最小二乘法拟合出与该交点往后的曲线最接近的直线，由公式（21）可求出剪胀角ψ：

2.3 与文献结果对比

利用前面计算方法，小麦的宏观参数计算结果如表2 所示.为了验证模拟结果，对比了文献中关于小麦的弹性模量等宏观参数.文中试验模拟的是100 kPa 围压下的应力应变关系，由表2 可知，小麦的宏观参数如弹性模量等值差异性较大，一方面，采用不同测试方法测定会导致数值误差；另一方面，小麦本身的硬度、含水量等初始状态都会对弹性模型产生较大的影响.文中获得的弹性模量值在文献的数值范围之内.

表2 小麦的宏观参数对比Table 2 Comparison of wheat maoropameters

3 结论

利用PFC3D离散元方法模拟室内小麦三轴试验，获得小麦三轴试验的应力-应变曲线，通过对比宏观试验曲线，标定了小麦的细观参数，得出以下结论：

（1）PFC3D数值模拟能灵活改变三轴试验的边界、初始条件等，可有效地模拟室内试验，是进行三轴试验的有效数值手段.

（2）数值试验结果表明，小麦颗粒材料的宏观参数弹性模量、泊松比、内摩擦角和剪胀角与细观参数颗粒刚度、颗粒间摩擦系数等具有良好的相关性.

（3）通过数值试验，获得小麦的细观参数，对比分析得到的宏观参数符合要求.

［1］Cundall P A，Strack O D L.A discrete numerical model for granular assemblies［J］.Geote -chnique，1979，29（1）:47-65.

［2］罗勇，龚晓南，连峰.三维离散颗粒单元模拟无黏性土的工程力学性质［J］.岩土工程学报，2008，30（2）：292-297.

［3］耿丽，黄志强，苗雨.粗粒土三轴试验的细观模拟［J］.土木工程与管理学报，2011，28（4）：24-29.

［4］邵磊，迟世春，张勇，等.基于颗粒流的堆石料三轴剪切试验研究［J］.岩土力学，2013，34（3）:711-720.

［5］李识博，王常明，王念秦，等.黄土三轴试验的颗粒流数值模拟［J］.中国公路学报，2013，26（6）：22-29.

［6］Cheng Y P，Nakata Y，Bolton M D.Discrete element simulation of crushable soil ［J］.Geotechnique，2003，53（7）：633-641.

［7］Zhao Hong-hua，Shao Long-tan，Ji Shunying.Numerical simulation of triaxial test on the dense sand by DEM［J］.Geotechnical Special Publication，2011，220：233-241.

［8］Backstrom A，Antikainen J，Backers T.Numerical modelling of uniaxial compressive failure of granite with and without saline porewater［J］.International Journal of Rock Mechanics&Mining Sciences，2008，45（7）:1126-1142.

［9］罗勇.土工问题的颗粒流数值模拟及应用研究［D］.杭州：浙江大学，2007.

［10］程绪铎，安蓉蓉，曹阳，等.小麦粮堆弹性模量的实验测定与研究［J］.粮食储藏，2009，38（6）：22-25.

［11］程绪铎，陆琳琳，石翠霞.小麦摩擦特性的试验研究［J］.中国粮油学报，2012，27（4）：15-19.

［12］石翠霞.筒仓内小麦、玉米堆的压缩特性、仓壁压应力及储量总重量的研究［D］.南京：南京财经大学，2012.

［13］梁醒培，杨伯源，王西院.深筒仓动态压力数值模拟［J］.合肥工业大学学报：自然科学版，2008，31（3）：427-429.