一套可行的PFC力控制实现循环加载方法（附同步视频教程+离散元学习资料）

详细信息

作者优秀: 优秀教师/意见领袖/博士学历/特邀专家/独家讲师
内容充实: 工程案例
平台推荐: 主编推荐/内容稀缺/全网独家

7月前浏览11461

导读：

我认识Lobby老师（仿真秀优秀讲师）是通过他连续在仿真秀平台发布原创文章，并且深受用户喜爱（持续有用户订阅），这让我很惊奇！于是我翻看了他发布的所有文章，毫不犹豫发起对他约稿。

后来进一步获悉：他毕业于上海同济大学，一直在使用离散元技术搞科研和解决复杂工程问题，尽管工作十分饱和，工作之余在仿真秀平台分享他这些年积累的离散元经验和学习心得，希望能对理工科学子和研发工程师成长有些许帮助，如有不当，欢迎批评指正，共同进步！

以下是正文：

一、基本原理

PFC中的计算元素包括ball、clump、wall，其中ball和clump是参与力学计算的，也就是说可以由接触力来计算其速度位移。但wall是不可以进行力学计算的，只能使用位移进行控制，我们当然可以通过伺服原理来实现位移和力的相互转化，但是使用wall直接进行力加载始终有些不顺畅，原因在于力的传递需要时间，而伺服原理根本而言给的是试样内部的力，而不是边界的力，这个和现实的加载又有些出入。并且利用伺服实现力的加载，对伺服系数敏感度较高，有可能会出现力的波动情况。

clump是PFC内置的可以模拟不规则形状的单元，其力学原理和ball一样，只是形状不是圆形而已。这里便采用矩形的clump来模拟加载板，这样和实际也可以对的上。当然模拟也可以高出现实，在于我们可以设置加载板间是可以没有力的，这样便可以实现加载板的互相穿透，摆脱了现实真三轴试验的仪器限制。最终的加载板如图所示，加载板中的pebble数量一共800个，相对于模拟的上万颗粒数而言不算很多了，所以对计算效率的影响不会很大。

图1 clump加载板

二、clump实现预压和围压

Clump直接可以施加外力，所以根据尺寸和应力计算出力就可以直接施加在加载板上面了。但是有时候如果应力变化比较大的话，对于散体材料肯定会出现很大的扰动，所以可以运用分级加载来实现应力缓慢的变化。这里采用二次函数实现，如果是指数函数的话可能还有点偏大。

图2 分级加载曲线

这里可以根据当前应力wsxx00和目标应力mubiaoxx来计算系数a:

a= (wsxx00-mubiaoxx)/ 10^2

则当前加载阶段的应力计算为：

txx=mubiaoxx ((math.abs(n-10))^10)*a

成样后模型的状态为：

图3 模型图

则在预压阶段应力的变化为：

图4 预压应力变化图

围压阶段应力的变化为：

图5 围压应力变化图

三、等偏应力幅的循环加载

首先我们一下加载路径，如图6所示，对试样的纵向施加余弦应力振动。则纵向应力的值为：

tyy=txx fengzhi*0.5-fengzhi*0.5*math.cos(2*math.pi*pinlv*mech.age)

其中：tyy为竖向应力，txx为横向应力，fengzhi为偏应力峰值，pinlv为加载频率，即1s内有多少个余弦波，mech.age为当前时间。

图6 余弦应力路径

我们监测加载过程中tyy与wsyy的变化，如图7所示，蓝色线为tyy，红色线为wsyy。其中tyy是我们预设的应力，也就是施加给加载板的力，wsyy为加载板和试样的接触力。这里给的围压为100kPa，偏应力峰值为50kPa，则按照假设，纵向应力应该在100-150kPa间作余弦波动。从图中可以看出，给加载板施加的力与加载板监测到的力之间，有一定的滞后，在数值上也不是完全一致，材料的动力特性以及加载频率影响比较大，图8和图9为频率为2以及50的曲线，可以看出加载频率越小，滞后性越小。

图7 频率为10的应力曲线

图8 频率为2的应力曲线

图9 频率为50的应力曲线

图10为频率为10 的加载频率的横向应力与竖向应力的变化，可以看到由于加载的滞后性以及加载速度的问题，横向应力wsxx并不是一直维持在100kPa，而是有一定的波动，这个波动随着频率的减小也会不怎么明显。这也是一样的道理，wsxx为试样内部力，而txx为仪器给的力，所以txx可以维持不变，而wsxx由于惯性肯定会有所变化。

图10 横向与纵向的应力曲线

试样循环加载应力应变曲线如图11所示，注意这里的应力用的是wsyy，因为wsyy有一定的滞后性，所以采用tyy的话效果会好很多，tyy对应的也是现实中仪器给的应力值，而wsyy对应的是试样内部的应力值。由图可以看出，频率越小，这种应力的滞后性越不明显，频率比较大的时候，不光应力会出现滞后性，应力的值也会有所区别。

图11 tyy与wsyy曲线

在做不排水的循环试验时，我们更关注其累计变形。图12为频率为10的轴向应变随时间的关系，即粉色的曲线，首先曲线增大，随后在某一个区间内往复运动。循环加载中的变形可以分为两个部分，一部分为不可恢复的塑性变形，即蓝色的区域，一部分为可恢复的弹性变形，即粉色曲线的区域。从图中可以看出累计塑性变形先增大，最后稳定在一个数值。总的变形先是大部分转化为塑性变形，之后都是转化为弹性变形。

图12 轴向应变曲线

四、参数分析

1、加载频率的影响

对试样进行2Hz、10Hz、50Hz的加载，整理后可得下图的曲线对比，图13为tyy曲线的对比，可以看到频率越小，第一次加载产生的应变越大，其最终的应变也是越大的。10Hz表现出的滞回圈最大，这是因为频率比较小的时候，试样的变形比较彻底，而频率比较大时，试样来不及大的变形。

图13 tyy曲线

图14为wsyy曲线，这里也可以看出为什么大频率下产生的应变会小，这是因为试样内部的应力没有达到仪器所给的应力，而低频率下应力基本差不多，所以最后的变形值也不会有太大的区别。相比于wsyy曲线，tyy曲线更加圆润一点，但是wsyy曲线更加能够体现试样内部的应力应变情况。

图14 wsyy曲线

图15为不同加载频率的轴向变形曲线，可以看出2Hz下的曲线能够更快的达到一个平稳阶段，这个和应力曲线也是相对应的。虽然2Hz下的累计变形最大，但是累计塑性变形却最小，而50Hz的加载频率可以发生更大的累计塑性变形。

图15 轴向变形曲线

2、阻尼系数的影响

在做单元试验的时候，我们通常会采用damp去耗散系统中的动能，使其能够更快的平衡。这个对于准静态条件来说是完全可以成立的，但是对于动力试验我们需要仔细去考虑damp的取值。一般来讲damp模拟的是介质对颗粒的粘滞效应，也就是说damp的值决定了你的材料是放在空气中，还是放在水里，或者是放在黏液中。

这里采用三种不同的damp值对其动力特性进行研究，即取damp为0.2、0.5、0.7。

图16和17分别为试样的tyy应力曲线以及wsyy应力曲线。可以看出，阻尼系数越小，其产生的累计变形越大。这个也是可以理解的，颗粒变形的速度会受到阻尼的影响，阻尼越大，颗粒速度减小的越快，于是其产生的变形越小。

图16 tyy应力曲线

图17为不同阻尼系数下的累计变形曲线，可以看出不同的阻尼系数其累计塑性变形不会有很大的区别，但是可恢复的弹性变形随着阻尼系数的减小而增大。

图17 累计变形曲线

3、CSR循环应力比的影响

CSR为偏应力与正应力的比值，这里做了三组试验，分别为CSR=0.1、0.2、0.5。图18为三种CSR的应力应变曲线，可以看出CSR越大，其产生的应变越大。图19也能看出一样的效果。

图18 tyy图

图19 累计曲线图

五、PFC模拟仿真实践与技巧

本文使用clump模拟加载板，直接对加载板施加力来实现对试样的加载。提出了一套可行的模拟方法，并进行等偏应力的循环加载方法。对模拟进行参数分析，分析了加载频率、阻尼系数和CSR对模拟的敏感性，也证明了模型的正确性。为此笔者在仿真秀平台同步发布了视频教程，感兴趣的朋友可以关注我的专栏。点击下图即可查看视频教程