首页/文章/ 详情

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

2年前浏览6266

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。成样如下图。

图片

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分8w11th Gen Intel(R) Core(TM) i7-11700F @ 2.50GHz
颗粒膨胀法2天8w11th Gen Intel(R) Core(TM) i7-11700F @ 2.50GHz

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

7.文件领取

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


代码&命令二次开发岩土水利离散元PFC
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2022-09-08
最近编辑:2年前
素墨
硕士 个人gzh:离散元数值模拟交流
获赞 53粉丝 56文章 9课程 0
点赞
收藏
作者推荐

周期性边界真三轴参数标定研究示例

真三轴试验是单元体尺度下土工试验中非常经典的试验,是研究土体力学响应的基本工具,它能够得到在特定应力路径下的土单元体应力应变关系。真三轴试验对岩土力学的相关理论发展(比如本构理论)具有重要的意义。周期性边界是离散元中的边界的一种,如果模型的上下边界为周期性边界,颗粒如果从上往下运动透过下边界,那么这个颗粒将会从上边界运动到模型域内,如下图。但是,目前基于离散元(DEM)的模拟大多采用刚性墙作为边界,从而控制土单元试样的应力路径。由于刚性墙的边界效应比较强,对于一些比较特殊颗粒级配、特殊颗粒形状的试样的模拟结果其实不太理想。目前,很多科研文章在用真三轴标定参数的时候,很多都采用周期性边界的真三轴试验。本文,将基于PFC6.0模拟低应力水平的三轴压缩,并复现了Ciantia[1]关于枫丹白露砂的参数研究,其中主要的难点在于编写周期性边界的应力伺服程序(参考了Help文件),试样内孔隙比、配位数、颗粒级配、应力的测量。·模型描述试样尺寸:3mm×3mm×3mm立方体边界:整个模型没有用wall,立面体边界都是周期性边界土样:模拟砂土,特定颗粒级配,采用赫兹接触模型颗粒级配:接触参数:(hertz接触模型)并禁止颗粒旋转!!!·建模流程首先是生成试样、然后在等向压力为10kPa下预压到制定孔隙比(通过调节颗粒的摩擦系数)、接着各向同性固结到围压为100kPa、最后在z方向施加偏压。·结果生成的试样并具有特定的级配:预压后得到了想要的孔隙率大概0.385:各向同性固结到100kPa,看看此时的力链,还是很均匀的,边界上并没有特别的应力集中:这里我们给出竖向应力(注意不是剪应力!)-竖向应变的关系图可以清晰的看到有一个小的应变软化的阶段,说明我们的试样处于一个稍微偏密的状态。然后是体应变与竖向应变的关系:也可以看出是先剪胀后剪缩的。最后是和试验对比的应力应变关系:(左本次模拟,右Ciantia文章中给出的结果)偏应力-竖向应变:左右体应变-竖向应变:·拓展对于复杂颗粒形状,比如PFC中由两个球粘接而成的clump,在标定其参数的时候,使用墙边界可能会导致很严重的边界效应,而周期性边界条件下的三轴压缩能够很好的解决这个问题!·说明本文复现了一篇文献中使用三轴压缩标定参数的研究,您可以借助这份代码进行您自己研究相关的参数标定。由于创作不易,原代码将有偿提供,并且承诺一次购买,全面答疑,如果使用过程中有任何疑问可以添加qq。代码链接跳转。代码中包含5个文件(前四个分别是制样、预压、施加围压、施加三轴压缩,最后一个是包含所有用到的fish函数文件)四个结果文件分别是以上每步生成的结果文件。参考文献[1]CiantiaMO,BoschiK,ShireT,etal.Numericaltechniquesforfastgenerationoflargediscrete-elementmodels[J].ProceedingsoftheInstitutionofCivilEngineers-EngineeringandComputationalMechanics,2018,171(4):147-161.

未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习 福利任务 兑换礼品
下载APP
联系我们
帮助与反馈