PFC2D大尺度离散元快速成样技术

前言

离散元的成样方法有很多，比如颗粒膨胀法、分层压实法、落雨法等。这些成样方法都能够很好地生成均匀的试样，但是面临大尺度的离散元模型试验模拟时，巨大的颗粒数目使得个人PC的计算能力日益窘迫。本文将介绍由 Matteo Oryem Ciantia 在文章 Numerical techniques for fast generation of large discrete-element models 提出的大尺度离散元模型快速成样方法。

1. 制作Brick单元

PFC提供了brick功能，能够将一小块brick试样copy到你想要的任意尺寸。因此成样第一步就是制作一块你想要的brick。

model new
model random 10001
model mechanical timestep scale

[hmin = 0]
[hmax = 0.1]
[r    = 0.6]

[rp   = 0.018]
[rv = rp * 5]
[rc   = rv +   0*rp]

[emodB = 1.0e8]
[kratio = 2.5]
[damp = 0.7]

[NumEnlarge = 2]
[grainPorosity = 0.2]
[density_ball = 2500.0]
[ballMin = 0.0006]
[ballMax = 0.0009]

model domain extent [0] [rc]  ...
                    [hmin] [hmax] ...
                    condition  periodic periodic


;generate balls
ball distribute porosity [grainPorosity]  ...
    radius [ballMin] [ballMax] box [0] [rc] [hmin] [hmax] ...
    resolution [NumEnlarge]    
ball attribute density [density_ball] damp [damp]     
contact cmat default model linear method deformability emod [emodB] kratio [kratio]


;expend method  
program call 'functions.fis'
[enlarge_ball] 

model solve
model save 'I0d20'

2. 进行Brick单元力学试验测试

这个步骤的目的是标定离散元参数，由于brick单元的边界都是周期性边界，因此单元试验也应当在周期性边界下进行。可进行的单元试验包括：三轴压缩（双轴压缩），三维单剪（二维单剪）。我写过一些关于周期性边界下的单元试验的帖子，可以参考：三轴压缩、双轴压缩（待开发，几天后推出）

3. 对Brick单元预压

预压的目的在于使得初始试样均匀地压实到制定的孔隙比，我们可以通过调节预压时的摩擦系数来得到我们想要的孔隙比试样。可以看到下图，力链还是比较均匀的。

model restore 'I0d20'
model domain tolerance 1e-3
program call 'servo_domain.fis'
model mechanical timestep automatic
model calm
ball fix spin
;zone gridpoint fix velocity
[do_xSevro = true]
[do_ySevro = true]

[txx = -5.0e3]
[tyy = -5.0e3]

[emodB = 1.0e8]
[kratio = 2.5]
[ballFrictionPre = 0]

contact cmat default model linear method deformability emod [emodB] kratio [kratio]

measure creat id 1 radius 0.01

;define ball  friction property
ball property 'fric' @ballFrictionPre


fish callback add @sevro_domain  -1.0

fish history name 41 @s_xx
fish history name 42 @s_yy
fish history name 43 @porosityHistory

[tol =  5e-3]
fish define stop_me
  if math.abs((s_yy - tyy)/tyy) > tol
    exit
  endif
  if math.abs((s_xx - txx)/txx) > tol
    exit
  endif
  if mech.solve("ratio-average") > 1e-5
    exit
  endif
  stop_me = 1
end

ball attribute displacement multiply 0.0
model solve fish-halt @stop_me

model save 'P0d20f0d3'

4. K0固结

预压之后，对Brick进行K0固结，需要定义竖向压力tzz和侧压力系数K0。成样如下图。

model restore 'P0d20f0d3'
model calm
ball fix spin

[do_xSevro = true]
[do_ySevro = true]

[k0 = 0.43]
;[K0 = 1.0]
[tyy = -100.0e3]
[txx = tyy*k0]


[ballFrictionC = 0.3]

ball property 'fric' @ballFrictionC




[tol =  5e-3]
[stop_me = 0]
fish define stop_me
  if math.abs((s_yy - tyy)/tyy) > tol
    exit
  endif
  if math.abs((s_xx - txx)/txx) > tol
    exit
  endif
  if mech.solve("ratio-average") > 1e-5
    exit
  endif
  stop_me = 1
end

ball attribute displacement multiply 0.0
model solve fish-halt @stop_me
;model save 'C0d20f0d3ISO'
model save 'C0d20f0d3K0'

5. 组装成大尺度试样并进行边界处理

组装成大尺度试样，由于原来是Domain边界，现在要变为墙边界，墙生成之后，墙与颗粒之间的重叠会很大，此时我们需要通过修改球墙边界的rgap属性，让球与墙的接触力到达一个合理的值，能减少之后模型的平衡时间。最后成样的效果如下图，具有8万颗粒。

model restore 'C0d20f0d3K0'
fish callback remove @sevro_domain

model domain condition periodic
brick make id 1
brick export id 1 skip-errors filename 'Brick0902'

ball delete
[emodBall = 1.0e8]
[kratioBall = 2.5]
[emodWall = 1.0e9]
[kratioWall = 2.5]

[bymax = domain.max.y()]
[bymin = domain.min.y()]
[bxmax = domain.max.x()]
[bxmin = domain.min.x()]
[nxbrick = 6]
[xmin = bxmin]
[xmax = bxmin + nxbrick*lx]
[nybrick = 10]
[ymin = bymin]
[ymax = bymin+nybrick*ly]

model domain extent [xmin-rp] [xmax + rp] ...
                    [-1*rp+ymin] [ymax + rp] ...
                    condition destroy destroy
brick assemble id 1 origin [xmin] [ymin] size [nxbrick] [nybrick]
model calm 
ball fix spin

wall creat id 11 name 'plane_left' group 'fricless' ...
                         vertices  ([xmin],[ymin-0.8*rp]) ([xmin],[ymax+0.8*rp])

wall creat id 12 name 'top' group 'fricless' ...
                         vertices  ([xmin-0.8*rp],[ymax]) ([xmax+0.8*rp],[ymax])
wall creat id 13 name 'bottom' group 'fricless' ...
                         vertices  ([xmin-0.8*rp],[ymin]) ([xmax+0.8*rp],[ymin])
wall creat id 14 name 'plane_right' group 'fricless' ...
                         vertices  ([xmax],[ymin-0.8*rp]) ([xmax],[ymax+0.8*rp])

program call 'servo_walls.fis'
fish callback add @sevro_walls  -1.0
[do_xSevro = true]
[do_ySevro = true]
history delete
history interval 100
fish history name 101 @s_xx
fish history name 102 @s_yy
fish history name 103 @wsxx
fish history name 104 @wsyy

program call 'Modify_gap.fis'
[targetForce = AverageContactForce]

contact cmat default type 'ball-facet' model linear method deformability emod [emodWall] kratio [kratioWall]

model cycle 1
[yPressure = 100e3]
[xPressure = k0 * yPressure]
[ytargetForce = yPressure*(xmax-xmin)]
[xtargetForce = xPressure*(ymax-ymin)]
[ Modify_BallWall(xtargetForce, 'plane_left')]
[ Modify_BallWall(ytargetForce, 'top')]
[ Modify_BallWall(ytargetForce, 'bottom')]
[ Modify_BallWall(xtargetForce, 'plane_right')]

[tol =  5e-3]
[stop_me = 0]
fish define stop_me
  if math.abs((wsyy - tyy)/tyy) > tol
    exit
  endif
  if math.abs((wsxx - txx)/txx) > tol
    exit
  endif
  if mech.solve("ratio-average") > 1e-5
    exit
  endif
  stop_me = 1
end
model solve fish-halt @stop_me

model save 'A0d20f0d3'


[ymin1 = wall.pos.y(wpBottom)]
[measure_ini_geostress(wlx/2, 0.08, wly,ymin1 , 'test6', 100000001)]

6. 总结与讨论

当颗粒数超过6w时，对于普通的PC电脑来说，一般的成样方法可能需要花费2-3天，而本文提出的成样方法对于二维大概一个小时、三维一天就能够算完。

这里对比一下这个案例里，颗粒膨胀法和本文的方法花费的时间：

成样方法	成样时间	颗粒数目	CPU型号
本文	<10分	8w	11th Gen Intel(R) Core(TM) i7-11700F @ 2.50GHz
颗粒膨胀法	2天	8w	11th Gen Intel(R) Core(TM) i7-11700F @ 2.50GHz

最后给出文章链接：Numerical techniques for fast generation of large discrete-element models

7.文件领取

以上代码需要使用几个包含fish函数的文件，关注公众号“离散元数值模拟交流”，并邀请3人，截图发后台领取全部代码。