PFC3D中文讲解例题

1、通过例子学习几种常见命令例二：PFC3D目录下的GuideStartfooting.dat，程序如下（注意：本程序与原程序不同，特加了plot set rotation (24.0,352.0,340.0)这句以使得模型转过一定角度，让用户看的更清楚，另外还将原程序一分为二，主要是为了更好地学习每个命令的作用），以下是程序代码。;fname: footing.DAT (tutorial example for PFC3D)new ; clear program state to begin new problemset random ; reset random-number generato

2、rtitle 'Tutorial Example'wall id=1 face ( 0, 0, 0) (10, 0, 0) (10, 0 -5) ( 0, 0,-5)wall id=2 face ( 0, 0, 0) ( 0, 0,-5) ( 0, 2,-5) ( 0, 2, 0)wall id=3 face (10, 0,-5) (10, 0, 0) (10, 2, 0) (10, 2,-5)wall id=4 face (10, 2,-5) (10, 2, 0) ( 0, 2, 0) ( 0, 2,-5)wall id=5 face ( 0, 0,-5) (10, 0,-5

3、) (10, 2,-5) ( 0, 2,-5)wall id=6 face ( 0, 0, 0) ( 0, 2, 0) (10, 2, 0) (10, 0, 0)gen id=1,750 rad 0.15,0.20 x=0,10 y=0,2 z=-5,0;plot create Footingplot set title text 'Basic collection of particles in a box'plot set rotation (24.0,352.0,340.0)plot add ball yellowplot add wall whiteplot add a

4、xes brownplot showpause;wall id=1 kn=1e8 ks=1e8wall id=2 kn=1e8 ks=1e8wall id=3 kn=1e8 ks=1e8wall id=4 kn=1e8 ks=1e8wall id=5 kn=1e8 ks=1e8wall id=6 kn=1e8 ks=1e8;prop density 2000 kn 1e8 ks 1e8ini rad mul 1.47;set hist_rep=5hist ball zvel 3 1 0hist diagnostic mufset dt dscalecycle 3000以下是程序中逐条命令的解释

5、：1、我们首先定义了模型的边界，使用WALL命令，如下：wall id=1 face ( 0, 0, 0) (10, 0, 0) (10, 0 -5) ( 0, 0,-5)wall id=2 face ( 0, 0, 0) ( 0, 0,-5) ( 0, 2,-5) ( 0, 2, 0)wall id=3 face (10, 0,-5) (10, 0, 0) (10, 2, 0) (10, 2,-5)wall id=4 face (10, 2,-5) (10, 2, 0) ( 0, 2, 0) ( 0, 2,-5)wall id=5 face ( 0, 0,-5) (10, 0,-5) (10

6、, 2,-5) ( 0, 2,-5)wall id=6 face ( 0, 0, 0) ( 0, 2, 0) (10, 2, 0) (10, 0, 0)这些命令创造了一个由六个墙面围成的箱子，每个墙面都赋予了一个独一无二的ID号。每个墙面四边形通过给出四边形的四个顶点坐标来定义，即四个顶点连成四条线，然后再围成一个四边形墙面。注意：顶点的顺序决定了墙面是否被激活，激活的墙面可以自动识别与小球的接触，不激活的话就会被小球穿透，2、然后在箱子里填充小球，用Generate命令，也可简写，如下：gen id=1,750 rad 0.15,0.20 x=0,10 y=0,2 z=-5,0这个命令规定了

7、小球填充的空间范围，小球的大小以及填充的总个数。Id=1,750的含义是：ID号从1到750其实就生成了750个小球，并给每个小球指定了各自的ID号。Rad 0.15,0.20的含义是：小球半径从0.15m到0.20m随机选择x=0,10 y=0,2 z=-5,0的含义是：小球填充的空间是0<x<10，0<y<2，-5<z<0的一个立方体空间。如果填充的空间太小或者小球的半径设置的太大，以至于这个空间容不下这么多的小球，那系统会自动减少生成的球的个数，并且系统会提示实际生成的球的个数。3、模型的显示用命令plot，程序如下：plot create Footi

8、ngplot set title text 'Basic collection of particles in a box'plot set rotation (24.0,352.0,340.0)plot add ball yellowplot add wall whiteplot add axes brownplot show从上面命令可以看出plot的用法规律，即与plot显示相关的命令的前面都用plot引出，然后后面再加相关的命令，如create Footing的意思就是创建一个名字为Footing的可视化窗口（这个在菜单栏点击“window”，下拉菜单可发现多了一项vi

9、ew 1Footing，用户可自己看下）；set title text是设定标题， set rotation 是设定模型的角度，以方便观察；add ball yellow是设定球的颜色为黄色；add wall white是设定墙面（也就是边界面）为白色；add axes brown是设定坐标轴为棕色； show是显示模型的命令。4、定义墙面的法向刚度（normal stiffness）和切向刚度（shear stiffness）wall id=1 kn=1e8 ks=1e8wall id=2 kn=1e8 ks=1e8wall id=3 kn=1e8 ks=1e8wall id=4 kn=1e

10、8 ks=1e8wall id=5 kn=1e8 ks=1e8wall id=6 kn=1e8 ks=1e8其中kn表示法向刚度，ks表示切向刚度，单位为N/m5、定义球的密度，法向刚度和切向刚度prop density 2000 kn 1e8 ks 1e8使用命令property，可简写为prop，定义密度density为2000kg/m3，法向刚度和切向刚度均为108N/m6、改变小球半径ini rad mul 1.47ini应该是initial的简写，用于初始化小球半径（有待考证）rad表示半径mul是multiply的简写，是倍数的意思。1.47是倍数因子，表示现在半径是原来半径的1.

11、47倍。7、设定监测变量和采样间隔时间set hist_rep=5hist ball zvel 3 1 0hist diagnostic mufhist是history历史记录是简写，set hist_rep=5表示每隔5个运算步记录一次变量的值（默认的是10个运算步记录一次）；本程序设定了两个监测对象，一个是hist ball zvel 3 1 0，表示的是中心坐标最接近x=3，y=1，z=0的球的z方向的速度被监测和记录（在这个例子中，282号球最接近这个坐标，因此作为监测球）；另一个是hist diagnostic muf，是用来监测系统的不平衡力，主要用于静力学分析，当这个值趋于很小时

12、，表示系统达到稳定状态。8、设定求解选项set dt dscale因为这是一个静态分析，因此我们开启density-scaling选项（简写为dscale）来优化求解效率，dt应该是这个选项中的某个模式（有待考证）。9、设定求解步数cycle 3000当写入这句程序后，计算就开始了，计算的步数，时间步，总时间，不平衡力平均值和最大不平衡力这些信息均显示在屏幕上，并且每隔5个计算步更新一次（这是set hist_rep=5的作用），图像动画每隔20个计算步更新一次（默认），这个更新时间可用命令set pinterval来改变。另外，阻尼已经在运动方程中被默认设定，用于吸收振动能量，更快地达到平衡

13、状态。10、查看结果求解完毕后，用户可输入命令plot his 1来查看历史记录中第一个变量的数值曲线。同理可输入命令plot his 2来查看历史记录中第二个变量的数值曲线。由曲线可知系统在cycle=1000时已趋于稳定，故可以将之前的3000步改为1000步，以减少计算时间（个人愚见，仅供参考）。先对plot补充点说明：1、plot create *命令可用于创建多个图像显示窗口，但每次只能看一种结果，可以通过菜单栏里的window选项卡来切换窗口，查看想看的结果2、如只想显示小球，则使用命令plot ball，如只想显示墙面，则使用命令plot wall，如二者都想显示则plot ba

14、ll wall，如运行完毕后想显示别的结果，如接触力则plot cforce（contact force的简写）。总而言之，plot后面加你想显示的东东即可，用户需灵活使用。此外，程序中的注释用英文的分号; ，后面加注释的内容即可。第二天的讲解中小球都是自由状态，无外力的作用，因此有必要讨论施加外力的情况。（1） 对小球施加重力作用;fname: footing.DAT (tutorial example for PFC3D)new ; clear program state to begin new problemset random ; reset random-number genera

15、tortitle 'Tutorial Example'wall id=1 face ( 0, 0, 0) (10, 0, 0) (10, 0 -5) ( 0, 0,-5)wall id=2 face ( 0, 0, 0) ( 0, 0,-5) ( 0, 2,-5) ( 0, 2, 0)wall id=3 face (10, 0,-5) (10, 0, 0) (10, 2, 0) (10, 2,-5)wall id=4 face (10, 2,-5) (10, 2, 0) ( 0, 2, 0) ( 0, 2,-5)wall id=5 face ( 0, 0,-5) (10, 0,

16、-5) (10, 2,-5) ( 0, 2,-5)wall id=6 face ( 0, 0, 0) ( 0, 2, 0) (10, 2, 0) (10, 0, 0)gen id=1,750 rad 0.15,0.20 x=0,10 y=0,2 z=-5,0;plot create Contact_Forceplot set title text 'Model at equilibrium'plot add ball yellowplot add wall whiteplot add cforce blackplot add axes brownplot show;wall i

17、d=1 kn=1e8 ks=1e8wall id=2 kn=1e8 ks=1e8wall id=3 kn=1e8 ks=1e8wall id=4 kn=1e8 ks=1e8wall id=5 kn=1e8 ks=1e8wall id=6 kn=1e8 ks=1e8;prop density 2000 kn 1e8 ks 1e8ini rad mul 1.47;set hist_rep=5hist ball zvel 3 1 0hist diagnostic mufset dt dscale;set grav 0 0 -9.81prop fric 1.0;cyc 1000此程序与第二天的程序几乎

18、相同，唯一的不同是多加了两句程序：set grav 0 0 -9.81prop fric 1.0第一句是设定重力加速度为9.81m/s2，方向为Z轴负方向，grav是重力gravity的简写；第二句是设定球与球之间的摩擦系数，fric是friction的简写，在该程序中设定摩擦系数为1.0。运行结果如下图所示：（2） 将上顶面换做一小块下降的板，压缩箱子里的小球程序如下：;fname: footing.DAT (tutorial example for PFC3D)new ; clear program state to begin new problemset random ; reset

19、random-number generatortitle 'Tutorial Example'wall id=1 face ( 0, 0, 0) (10, 0, 0) (10, 0 -5) ( 0, 0,-5)wall id=2 face ( 0, 0, 0) ( 0, 0,-5) ( 0, 2,-5) ( 0, 2, 0)wall id=3 face (10, 0,-5) (10, 0, 0) (10, 2, 0) (10, 2,-5)wall id=4 face (10, 2,-5) (10, 2, 0) ( 0, 2, 0) ( 0, 2,-5)wall id=5 fac

20、e ( 0, 0,-5) (10, 0,-5) (10, 2,-5) ( 0, 2,-5)wall id=6 face ( 0, 0, 0) ( 0, 2, 0) (10, 2, 0) (10, 0, 0)gen id=1,750 rad 0.15,0.20 x=0,10 y=0,2 z=-5,0;plot create Contact_Forceplot set title text 'Model at equilibrium'plot add ball yellowplot add wall blueplot show;wall id=1 kn=1e8 ks=1e8wall

21、 id=2 kn=1e8 ks=1e8wall id=3 kn=1e8 ks=1e8wall id=4 kn=1e8 ks=1e8wall id=5 kn=1e8 ks=1e8wall id=6 kn=1e8 ks=1e8;prop density 2000 kn 1e8 ks 1e8ini rad mul 1.47;set hist_rep=5hist ball zvel 3 1 0hist diagnostic mufset dt dscale;set grav 0 0 -9.81prop fric 1.0;delete wall 6wall id 10 face (0,0,0) (0,2

22、,0) (2,2,0) (2,0,0)wall id 10 ks 1e8 kn 1e8 fric 1.0wall id 10 zvel=-0.5e-3property n_bond=5e5 s_bond=5e5hist wall zforce id 10cyc 6000运行结果如图所示与只施加重力的情况相比，这个程序多加了如下几句代码：delete wall 6wall id 10 face (0,0,0) (0,2,0) (2,2,0) (2,0,0)wall id 10 ks 1e8 kn 1e8 fric 1.0wall id 10 zvel=-0.5e-3property n_bond

23、=5e5 s_bond=5e5hist wall zforce id 10其中，delete wall 6是删除编号为6的面，即箱子的顶面。wall id 10 face (0,0,0) (0,2,0) (2,2,0) (2,0,0)wall id 10 ks 1e8 kn 1e8 fric 1.0wall id 10 zvel=-0.5e-3是建立一个ID为10的墙面，设定面的法向刚度、切向刚度、摩擦系数和Z方向的速度，这个速度的单位是m/step，因为开启了density scaling（程序set dt dscale的作用），如果是完全的动力学分析的话，单位将会是m/s（这个地方可解释为

24、：静力学只考虑系统达到平衡稳态的情况，故单位与时间无关；动力学关注的是每一时刻的状态，故单位与时间有关）property n_bond=5e5 s_bond=5e5是为所有的接触面定义法向和切向强度，以此来模拟初始的接触力（应该是类似于预紧力的作用，值得商榷），单位为Nhist wall zforce id 10用于监测下降的板的Z方向的力的变化情况，计算完毕后如想查看，只需输入plot hist 3即可。结果如图所示这个力的曲线图与英文帮助手册上的图明显不同，但不知道是什么原因出的错，有待思考。下面附上英文帮助手册上的图。（1） 补充一句程序说明在每个程序的开头，经常有句set random

25、，这句话的意思是重置随机数，因为generate生成小球的时候，半径经常是个范围，需要随机生成，如果不加这句的话，随机得到的半径大小可能一直是某一组随机数，故需重置一下，这与vb编程中的随机数重置设定是一样的。（2）plot wall id=on 显示所有wall的ID（3）定义字符串变量命令：macro（手册上的例子均是将某字符串赋给某个变量，因此我取名字符串变量，如果是赋予具体数值，比如a=5之类的功能，应该有别的命令可实现，待后续）Macro xx x 0 10Macro yy y 0 10Macro zz z 0 10Macro size id 1 100 rad 0.5 0.5Mac

26、ro ball xx yy zz sizeMacro size id 1001 1100 rad 0.1 0.1 ;注意！Gen ball以上程序需注意两点：1、 macro命令的结构形式：macro+变量名+字符串，字符串要加单引号；2、 倒数第二句是重新定义了变量size，需要注意的是size重新定义的时候需加单引号，第一次定义不需要加。（4）一般分析流程这个流程应该是先建立模型，设定接触和材料属性，施加载荷和约束，然后再运行程序使达到稳态，最后根据需要改变条件继续求解，从而得出结果。（5）生成球体的两个命令命令generate和命令ballGenerate（简写为gen）在第二天的讲解中

27、已经说过，它可一次生成大量的球填充到某个空间内，ball命令是一次只生成一个球，这也是与generate的区别。Ball命令的格式为：ball rad * id * x * y * z *（*为各参数对应的值，需用户自己输入）（6）施加边界条件1、 给墙面施加速度和角速度Wall id 1 xvel 1e-3给ID号为1的墙面施加x方向速度，大小为1e-3Wall id 2 yspin 1e-3 x=2 y=2 z=0给ID号为2的墙面施加绕y轴转动的角速度，旋转中心坐标为（2,2,0）2、 给球施加恒定不变的速度或角速度为球施加速度或角速度，可用命令initial，change，property，但这种定义只能对球的速度进行初始化，不能保证以后球一直以恒定的速度运动下去，因此需要用到fix命令。程序如下：Ini xvel 5e-3 zspin 1e-3 range id 40Fix x zspin range id 40上面程序的意思就是先给ID为40的球赋予x方向的速度和z方向的角速度，之后用fix命令保证这两个速度恒定不变。3、 给球初始化运动参数Initial命令可初始化小球的力，