0 引言
PFC的结果往往都是离散化的结果,我们对于结果大部分都是采用ball的一些显示,比如速度场、位移场等,或者是接触的一些东西。当然我是支持利用离散化的结果去进行描述的,这点也符合离散元的理论特性。
但是当某些情况下,颗粒的离散化比较大的时候,我们需要将一部分颗粒集 合体作为单元来显示,作为某一部分的结果。这时候就是我们测量圆的概念,这点在很多地方也会遇到,比如式样内部的应力分布、空隙分布等。这个是单个或者少部分颗粒无法描述的性状。
这时候或许用有限元那种单元的可视化效果会更加好一点。
1 案例
以一个土体自重沉降的案例为例:
成样:
model newdef parwidth = 0.2hight = width*2rdmin = 0.009rdmax = 0.006poro = 0.4end@pardomain extent [-width*2.0] [width*2.0] ...[-hight*2.0] [hight*2.0]model random 10001wall generate box [-width*0.5] [width*0.5] ...[-hight*0.5] [hight*0.5] expand 1.5cmat default model linear method deform emod 100e6 kratio 1.5ball distribute porosity @poro ...radius [rdmin] [rdmax] ...box [-width*0.5] [width*0.5] ...[-hight*0.5] [hight*0.5] ...[-hight*0.5] [hight*0.5]ball attribute density 2200 damp 0.7model cycle 2000 calm 10model cycle 1model solveball delete range pos-x [-width] [-width*0.5]ball delete range pos-x [width*0.5] [width]ball delete range pos-y [-hight] [-hight*0.5]ball delete range pos-y [hight*0.5] [hight]ball delete range pos-z [-hight] [-hight*0.5]ball delete range pos-z [hight*0.5] [hight]model save "ball_sample"
沉降:
model restore "ball_sample"model gravity 9.8wall delete walls range id 2model cycle 1model solvemodel save "zizhong"
沉降后的结果如图;

可以从力链上分析这个应力是上小下大的,是一个比较典型的力的分布。
3 绘制
这里使用的概念是,先把计算范围网格化,然后以网格中心为测中心,最小的长度作为半径生成测量圆。后面把测量圆的相关数据存到rblock的extra中,ball、clump和rblock这些基础元素都有这种extra属性用来存储用户自定义的数据。
model restore "zizhong"def create_measure(x_min,x_max,y_min,y_max,z_min,z_max,n_x,n_y,n_z)x_length=(x_max-x_min)/float(n_x)y_length=(y_max-y_min)/float(n_y)z_length=(z_max-z_min)/float(n_z)rad=math.min(x_length,y_length)rad=math.min(rad,z_length)loop local n(1,n_x)x_pos=x_min x_length*(n-1) x_length*0.5loop local m(1,n_y)y_pos=y_min y_length*(m-1) y_length*0.5loop local k(1,n_z)z_pos=z_min z_length*(k-1) z_length*0.5commandmeasure create position [x_pos] [y_pos] [z_pos] radius [rad*0.5]rblock create box [x_pos-x_length*0.5] [x_pos x_length*0.5] ...[y_pos-y_length*0.5] [y_pos y_length*0.5] ...[z_pos-z_length*0.5] [z_pos z_length*0.5]endcommandendloopendloopendloopend@create_measure([-width*0.5],[width*0.5],[-hight*0.5],[hight*0.5],[-hight*0.5],[hight*0.5],2,4,4)def get_dataloop foreach mp measure.listid=measure.id(mp)rb=rblock.find(id)rblock.extra(rb,1)=measure.stress.zz(mp)endloopend@get_data
这里存储了z向的应力,测量圆分布如图:

把rblock的相关显示设置为:

与其相对应的rblock的效果为:

隐藏掉measure:

换种显示方式可以得到如下的结果:

配位数的分布如图:

本案例以方形的简单算例进行了计算,针对于复杂的模型,也可以先进行三角形或者多边形的rblock划分,然后再在rblock位置处生成测量圆的方式去做。理论应该是一致的,这里留给大家自己去拓展。