PFC初期简单练习程序

PFC简单程序1
Pfc2d程序1：用ball生成球
;fname: raft.dat
ball rad 0.25 id=1501 x=0.25 y=0.25
ball rad 0.25 id=1502 x=0.75 y=0.25
ball rad 0.25 id=1503 x=1.25 y=0.25
ball rad 0.25 id=1504 x=1.75 y=0.25
macro Raft 'id=1501,1504';组群
property density 2000 kn=1e9 ks=1e9 range Raft property n_bond 1e10 s_bond 1e10 range Raft
property pb_rad 1.0 pb_s 1e20 pb_n 1e20 range Raft property pb_kn 1e10 pb_ks 1e10 c_index 1 range Raft
ini xv 0 range Raft
fix x range Raft
ini yforce -5.0e5 range Raft
pl add ball red
pl show
Pfc2d程序2：用generate生成球
;fname: gen1.dat
new
set random;生成随机数
macro big 'id 1,20';用macro定义颗粒组
macro small 'id 51,100'
gen big rad=.6,.6 x=0,10 y=0,10
gen small rad=.3,.3 x=0,10 y=0,10
change c_index 0 range big
change c_index 1 range small
plot create fred
plot add ball lblue lmag
plot show
扩展：生成两种不同的球
new
set random;生成随机数
macro big 'id 1,2000';用macro定义颗粒组
macro small 'id 5001,10000'
gen big rad=.15,.18 x=0,40 y=0,20
gen small rad=.02,.03 x=0,40 y=20,50
change c_index 0 range big
change c_index 1 range small
plot create fred
plot add ball lblue lmag
plot show
wall id 1 ks=1e8 kn=1e8 node 0, 0 40, 0 wall id 2 ks=1e8 kn=1e8 node 40,0 40,20
wall id 3 ks=1e8 kn=1e8 node 40,20 0,20 wall id 4 ks=1e8 kn=1e8 node 0,20 0,0 wall id 5 ks=1e8 kn=1e8 node 40, 20 40,50 wall id 6 ks=1e8 kn=1e8 node 40, 50 0,50 wall id 7 ks=1e8 kn=1e8 node 0,50 0,20
wall id 8 ks=1e8 kn=1e8 node 40,50 40,60 wall id 9 ks=1e8 kn=1e8 node 40,60 0,60 wall id 10 ks=1e8 kn=1e8 node 0,60 0,50 wall id 11 ks=1e8 kn=1e8 node -2,-2 -2,-32 wall id 12 ks=1e8 kn=1e8 node -2,-32 42,-32 wall id 13 ks=1e8 kn=1e8 node 42,-32 42,-2
Pfc2d程序3：用wall生成一组墙
new
def make_walls
command
wall id=1 nodes=(0, 0) (5, 0)
wall id=2 nodes=(5, 0) (5, 3)
wall id=3 nodes=(5, 3) (0, 3)
wall id=4 nodes=(0, 3) (0, 0) endcommand
end
make_walls
plot wall id=on red
Pfc2d程序4：用wall生成一组墙
new
wall id 1 ks=1e8 kn=1e8 node 0, 0 10, 0 wall id 2 ks=1e8 kn=1e8 node 10, 0 10,10 wall id 3 ks=1e8 kn=1e8 node 10,10 0,10 wall id 4 ks=1e8 kn=1e8 node 0,10 0,0
gen id=1,500 rad=0.12,0.17 x=0,10 y=0,10;生成球
prop dens=1000 kn=1e8 ks=1e8;给球赋属性
ini rad mul 1.60;初始半径乘以1.6倍
cyc 8000;运行8000步
plot create the_view ;创建显示框
plot add wall black;显示墙为黑色
plot add ball yellow;显示球为黄色
plot add meas blue;显示测量圈为蓝色
plot show;显示
measure id 1 x 5 y 5 rad 2.5;测量圈1以（5,5）为中心半径为2.5 print measure 1 ;打印测量圈1
这里解释一下，这里测量的是孔隙率，滑动摩擦系数，应力和应变
Pfc2d程序5：扩展墙与修改墙的刚度
;fname: biax_1.DAT biaxial sample
new
SET random ;生成随机数
SET disk on ;二维中生成的圆盘
; ----------------------------------------------------
def make_walls ; 创建扩展的墙壁
extend = 1.0;扩展系数
_x0 = -extend*width
_y0 = 0.0
_x1 = width*(1.0 + extend)
_y1 = 0.0
command
wall id=1 kn=w_stiff nodes (_x0,_y0) (_x1,_y1)
end_command
_x0 = width
_y0 = -extend*height
_x1 = width
_y1 = height*(1.0 + extend)
command
wall id=2 kn=w_stiff nodes (_x0,_y0) (_x1,_y1)
end_command
_x0 = width*(1.0 + extend)
_y0 = height
_x1 = -extend*width
_y1 = height
command
wall id=3 kn=w_stiff nodes (_x0,_y0) (_x1,_y1)
end_command
_x0 = 0.0
_y0 = height*(1.0 + extend)
_x1 = 0.0
_y1 = -extend*height
command
wall id=4 kn=w_stiff nodes (_x0,_y0) (_x1,_y1)
end_command
end
def assemble ; assemble sample
s_stiff = 0.0 ; 初始刚度
n_stiff = 5e8 ;法向刚度
w_stiff = 5e8;墙的法向刚度
tot_vol = height * width * 1.0;总体积
rbar = 0.5 * (rlo + rhi);
num = int((1.0 - poros) * tot_vol / (pi * rbar^2)) ;产生的颗粒数目mult = 1.6 ; initial radius multiplication factor
rlo_0 = rlo / mult
rhi_0 = rhi / mult
make_walls ;运行前面的墙命令
command
gen id=1,num rad=rlo_0,rhi_0 x=0,width y=0,height
prop dens=1000 ks=s_stiff kn=n_stiff
end_command
ii = out(string(num)+' particles were created') ;输出字符串
sum = 0.0 ; get actual porosity
bp = ball_head
loop while bp # null
sum = sum + pi * b_rad(bp)^2
bp = b_next(bp)
end_loop
pmeas = 1.0 - sum / tot_vol ;最终孔隙率
mult = sqrt((1.0 - poros) / (1.0 - pmeas))
command
ini rad mul mult
cycle 1000
prop ks=5e8 fric 0.25
cycle 250
end_command
end
; ----------------------------------------------------
def cws ;改变侧壁刚度
command
wall id 2 kn=w_stiff
wall id 4 kn=w_stiff
end_command
end
; ----------------------------------------------------
macro zero 'ini xvel 0 yvel 0 spin 0'
SET height=12.0 width=6.0 rlo=0.075 rhi=0.100 poros=0.14
assemble
SET w_stiff= 5e7 ; make lateral wall stiffness=1/10 of ball stiffness
cws
cyc 500
zero
plot create assembly
plot add ball lorange wall black
save bt_ass.SAV
return
;可以生成节理
Jset id=2 dip=0 origin=(3.0,2.0) number=2 spacing=5.0
plot add contact red range jset 2
Plot show
Pfc2d程序6：塌落隧道
;fname: tun1.DAT
new ; clear program state to begin new problem
set random ; reset random-number generator
set disk on ; treat balls as disks of unit thickness set hist_rep 5;每5步记录一次
wall id=1 kn=1e8 ks=1e8 node 0 -5 10 -5
wall id=2 kn=1e8 ks=1e8 node 10 -5 10 0
wall id=3 kn=1e8 ks=1e8 node 10 0 0 0
wall id=4 kn=1e8 ks=1e8 node 0 0 0 -5
gen id 1 500 rad 0.08 0.13 x 0 10 y -5 -0.5
property density 2000 kn 1e8 ks 1e8
ini rad mul 1.50;颗粒初始半径乘以1.5
hist diag muf;记录不平衡力
hist diag mcf
hist ball yvel id 442;记录442球y方向的速度
hist ball ypos 3 0;记录y方向位置在（3,0）附近
set dt dscale;设计时步
cycle 3500
set grav 0 -9.81
prop fric 1.0
cycle 2500
print info;
pl create Model
pl set cap size 25
pl add ball yellow
pl add wall black
pl add cfor black
plot show
plot current 0
plot hist 1
pause
plot hist 1 2 begin 3000
pause
plot hist 3
pause
plot hist 4
pause
save tun1.sav
;------------
;fname: tun2.dat
rest tun1.sav;恢复tun1.sav
prop s_bond=2e5 n_bond=2e5 ; add contact bonds
; delete top wall and add footing raft made of balls del wall 3
call raft.dat
set dt auto ; turn off density scaling
prop xdisp 0.0 ydisp 0.0 ; set displacements to zero hist reset
hist diag muf
hist diag mcf
hist ball ypos id 1501
hist ball yvel id 1501
solve
plot current 0
plot hist 3 ymin -0.02
pause
plot create disp_view
plot add wall black
plot add disp black
plot show
pause
plot show model
save tun2.sav
;-------
;fname: tun3_unl.DAT
rest tun2.sav
;excavate tunnel beneath and to right of footing macro Hole ’circle center 3 -2.5 rad 1.3’
del ball range Hole
prop xdisp 0.0 ydisp 0.0 ; set displacements to zero cycle 8000
plot show disp_view
pause
plot show model
pause
save tun3_unl.sav
生成边界的命令流
☐Wall type cylinder end1 0 0 0, end2 0 0 1, rad 1,1, id=1, kn=1e6, ks=1e6, fric=0.2 ;生成半径为1的圆柱面；
☐
☐Wall type cylinder end1 0 0 0, end2 0 0 1, rad 0.0,1, id=1, kn=1e6, ks=1e6, fric=0.2 ;生成底面半径为1的圆锥面；
☐
☐Wall type cylinder end1 0 0 0, end2 0 0 1, rad 1,2, id=1, kn=1e6, ks=1e6, fric=0.2 ;生成上底面半径1，下底面半径为2的圆台面；
☐
☐Wall type spiral end1 0 0 0, end2 10 0 0, radin 1 radout 2, pitch=2, id=1, kn=1e6, ks=1e6, fric=0.2；生成5个螺纹的螺旋面；
☐
☐wall id=5 face (2,1,2) (5,1,2) (5,0,2) (2,0,2)
wall id=5 face (5,0,2) (5,1,2) (5,1,5) (5,0,5)
wall id=5 face (2,0,5) (5,0,5) (5,1,5) (2,1,5)
wall id=5 face (2,1,2) (2,0,2) (2,0,5) (2,1,5)
注意：当Generate命令使用no_shadow关键词时，颗粒只会在有效侧与阴影区生成，非阴影区则不能生成颗粒。

当在较复杂的非方体空间内生成颗粒时，No_shadow关键
词（结合边界函数）极为有用。

Pfc2d程序7：颗粒规则排列
;fname:GZKLSC.DAT Generates a sheet of close-packed balls new
title 'Tutorial example'
def hex
xc = x0
yc = y0
rc = radius
idc = id_start
r2 = 2.0 * radius
yinc = radius * sqrt(3.0)
loop row (1,n_row) ;每排1到n_row
loop col (1,n_col)
command
ball id=idc x=xc y=yc rad=rc
end_command
idc = idc + 1
xc = xc + r2
end_loop
yc = yc + yinc
xc = x0 + radius * (row - (row/2) * 2)
end_loop
end
set echo off
set x0=0.15 y0=0.15 radius=0.1
set id_start=100 n_col=7 n_row=8
hex
set echo on
非规则排列之半径扩大法
非规则排列之颗粒排斥法
颗粒半径非均匀分区分布
当每个区域都生成颗粒组并分别达到平衡后，使用命令DEL WALL ID 将隔墙删除。

也可以不设隔墙，采用坐标控制颗粒的生成区域，但这种方法生成颗粒组后，若经过CYCLE 平衡后，颗粒将不能完全停留在各自的原始区域。

;fname
new
set random
gen x 0 5 y 0 5 rad .07 .07 id=1,120 annulus 2.5,0 4,5
gen x 0 5 y 0 5 rad .1 .1 id=121,240 annulus 2.5,0 2.5,4
plot create the_annulus
plot add ball white shade on;添加球白色阴影
plot add axes mag;显示坐标轴为紫色
plot show
Pfc2d程序8：圆柱体中生成球
;FNAME filter_cylinder.dat
new;
set random;
wall id=1 type cylinder end1=(0.0 0.0 0.0) end2=(0.0 0.0 3.0) &
rad=(1.0 1.0) kn=1.0e8 ks=1.0e8 fric=0.2;
wall id=2 normal (0 0 1) origin (0.0 0.0 0.0);
wall id=3 normal (0 0 -1) origin (0.0 0.0 3.0);
;
def ff_cylinder
_brad = fc_arg(0);
_bx = fc_arg(1);
_by = fc_arg(2);
_bz = fc_arg(3);
_rlim = 1.0;
_rxy = sqrt(_bx*_bx + _by*_by);
if (_rxy + _brad )<= _rlim
ff_cylinder = 0;
else
ff_cylinder = 1;
end_if
end
gen id=1,500 rad=0.1 0.1 x=-1.0 1.0 y=-1.0 1.0 z=0.0,3.0 filter ff_cylinder tries 400000000;
prop kn=1.0e6 ks=1.0e6 dens=1.0e3 fric=0.4
plot creat filter_test;
plot add wall white id on;
plot add ball yellow;
plot add axes;
plot set back white;
plot show;
☐本例中定义了一个名为ff_cylinder的过滤器(filter)，用于在一个半径为1.0m，高3.0m的圆柱形空间内生成500个半径等于0.1m的颗粒。

☐需注意的是，ff_cylinder(即fname)必须是一个有返回值的函数，且其返回值只有0或1，返回值必须赋予与函数名同名的变量。

☐第4—7行生成一个由圆柱壁面和2个无限平面组成的封闭圆柱形空间。

☐第9—21行即为过滤器ff_cylinder函数的定义过程，其中第11—13行的作用是将“试产球”的半径和位置信息通过fc_arg(0), fc_arg(1), fc_arg(2), fc_arg(3) （注：fc_arg(N)是PFC3D的内嵌效用函数）分别传递给变量_brad, _bx, _by, _bz。

Generate生成颗粒的过程为不断尝试在特定位置生成特定大小的球，每次尝试时都会生成“试产球”的半径和位置的x, y, z坐标值，这些值分别保存在fc_arg(0), fc_arg(1), fc_arg(2), fc_arg(3)之中（必须注意对应关系），调用这些函数将“试产球”的状态信息传递给变量。

☐第14—20行的作用为设置过滤条件，是过滤器的关键语句。

_rxy为“试产球”在XOY平面上离圆柱壁面轴线的距离，当(_rxy + _brad) 小于_rlim=1.0m(圆柱半径)时，ff_cylinder = 0(接受)，即当“试产球”处于圆柱壁面内部时，该球被接受并生成，否则不予生成。

Fish:需注意的是FISH 中的变量名、函数名以及语句，都要写全，不建议缩写。

此外，代码行不能续行，不过可以通过引入新变量解决此问题。

还有一点，FISH 不区分大小写。

Pfc2d程序9：记录曲线
new
wall id 5 kn 1e8 nodes -5,0 5,0
ball id 1 x 0 y 1.1 rad 1
prop dens 1000 kn 1e8
set grav 0 -10
set dt max=2e-3;设置最大时步
def down_force
down_force = b_yfob(ball_head)
end
hist down_force ;记录时间
cyc 500
plot his 1 ;显示时间曲线
new
def box_geometry
n_stiff = 2e8
s_stiff = 1e8
xx = 5.0
yy = 3.0
end
box_geometry
macro wall1nodes 'nodes ( 0, 0) (xx, 0)' macro wall2nodes 'nodes (xx, 0) (xx, yy)' macro wall3nodes 'nodes xx yy 0 yy' macro wall4nodes 'nodes 0 yy 0 0' wall id=1 ks=s_stiff kn=n_stiff wall1nodes wall id=2 ks=s_stiff kn=n_stiff wall2nodes wall id=3 ks=s_stiff kn=n_stiff wall3nodes wall id=4 ks=s_stiff kn=n_stiff wall4nodes plot wall red id on
print wall prop
print macro
定义一个一维矩阵，并进行赋值
代码1
def fuzhi
array num(20)
loop i(1,20)
num(i)=i^2+1
"wanjiabaodian"=out(num(i))
Endloop
End
Fuzhi
代码2
def number
array even(21),odd(21);even为偶数数组，odd为奇数数组m1=0
m2=0
loop n(0,41)
a=n/2.0
b=int(n/2)
if a = b then
m2=m2+1
even(m2)=n
else
m1=m1+1
odd(m1)=n
endif
endloop
end
number
def wri
loop k(1,21)
hh=out(even(k))
endloop
command
pause
endcommand
loop k(1,21)
hh=out(odd(k))
endloop
end
wri。