PFC2D大尺度離散元快速成樣技術(shù)

前言

離散元的成樣方法有很多，比如顆粒膨脹法、分層壓實(shí)法、落雨法等。這些成樣方法都能夠很好地生成均勻的試樣，但是面臨大尺度的離散元模型試驗(yàn)?zāi)M時(shí)，巨大的顆粒數(shù)目使得個(gè)人PC的計(jì)算能力日益窘迫。本文將介紹由 Matteo Oryem Ciantia 在文章 Numerical techniques for fast generation of large discrete-element models 提出的大尺度離散元模型快速成樣方法。

1. 制作Brick單元

PFC提供了brick功能，能夠?qū)⒁恍Kbrick試樣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. 進(jìn)行Brick單元力學(xué)試驗(yàn)測(cè)試

這個(gè)步驟的目的是標(biāo)定離散元參數(shù)，由于brick單元的邊界都是周期性邊界，因此單元試驗(yàn)也應(yīng)當(dāng)在周期性邊界下進(jìn)行。可進(jìn)行的單元試驗(yàn)包括：三軸壓縮（雙軸壓縮），三維單剪（二維單剪）。我寫過一些關(guān)于周期性邊界下的單元試驗(yàn)的帖子，可以參考：三軸壓縮、雙軸壓縮

3. 對(duì)Brick單元預(yù)壓

預(yù)壓的目的在于使得初始試樣均勻地壓實(shí)到制定的孔隙比，我們可以通過調(diào)節(jié)預(yù)壓時(shí)的摩擦系數(shù)來(lái)得到我們想要的孔隙比試樣。可以看到下圖，力鏈還是比較均勻的。

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固結(jié)

預(yù)壓之后，對(duì)Brick進(jìn)行K0固結(jié)，需要定義豎向壓力tzz和側(cè)壓力系數(shù)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. 組裝成大尺度試樣并進(jìn)行邊界處理

組裝成大尺度試樣，由于原來(lái)是Domain邊界，現(xiàn)在要變?yōu)閴吔纾瑝ι芍螅瑝εc顆粒之間的重疊會(huì)很大，此時(shí)我們需要通過修改球墻邊界的rgap屬性，讓球與墻的接觸力到達(dá)一個(gè)合理的值，能減少之后模型的平衡時(shí)間。最后成樣的效果如下圖，具有8萬(wàn)顆粒。

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. 總結(jié)與討論

當(dāng)顆粒數(shù)超過6w時(shí)，對(duì)于普通的PC電腦來(lái)說，一般的成樣方法可能需要花費(fèi)2-3天，而本文提出的成樣方法對(duì)于二維大概一個(gè)小時(shí)、三維一天就能夠算完。

這里對(duì)比一下這個(gè)案例里，顆粒膨脹法和本文的方法花費(fèi)的時(shí)間：

成樣方法	成樣時(shí)間	顆粒數(shù)目	CPU型號(hào)
本文	<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

文章來(lái)源：離散元數(shù)值模擬交流