PFC-FLAC3D Coupling Examples

Punch Indentation of a Bonded Material

这个例子展示了一个粘合颗粒模型（BPM）被圆形冲头压缩时的压痕反应。冲床的应力/变形不能在PFC中评估，因为墙体只作为刚性边界。2

在创建一个耦合的PFC-FLAC3D模型之前，必须加载PFC模块，可以通过 " 工具 "►" 加载PFC "菜单项或者使用以下命令来加载模块。一旦一个项目被保存并加载了任何额外的模块，只要再次打开该项目，这些模块就会被加载。

PFC球是在FLAC3D的无效区域内创建的。由于PFC总是在大应变模式下运行，因此必须明确激活区域和结构元素的大应变模式。（the model largestrain command）

PFC球和团块的创建可以没有初始重叠，也可以有任意的初始重叠。在这种情况下，球的分布命令被用来在一个平行四边形体积内创建一个具有任意重叠的球云，并达到指定的体积分数。然后分配球的属性，并使用接触cmat命令来指定接触模型的行为。采取这些步骤后，可以允许球重新排列。模型循环命令中的calm关键字被用来周期性地消除系统中的所有动能。最后，所有溢出生成区域的球被删除。

;create a bed of zones for the balls to sit inside
zone create brick size 35 35 35 point 0 -4 -4 -4.5 point 1 4 -4 -4.5 ...
                                point 2 -4 4 -4.5 point 3 -4 -4 -1 group 'bed'
zone cmodel assign elastic
zone property young [materialMod] poisson 0.25
zone cmodel null ...
            range position-x -1.5 1.5 position-y -1.5 1.5 position-z -2 -1
;as the balls are packing into this region - 
;fix the gridpoint motions as the balls are settling
zone gridpoint fix velocity

;wrap the bed with walls
wall-zone create name 'bed' starting-zone 30013
;create a wall over the top 
wall generate name 'top' plane position 0 0 -1 

;create a cloud of balls that are arbitrarily overlapping  （任意重叠）
ball distribute porosity 0.4 box -1.5 1.5 -1.5 1.5 -2 -1 radius 0.04 0.06

;set the ball density and local damping coefficient
ball attribute density 2600 damp 0.7
;set the default contact behavior - 
;the deformability method sets properties of the linear portion 
;of the contact model
contact cmat default model linearpbond method deformability ...
                   emod [materialMod] kratio 1.0 
;allow the balls to rearrange, nulling the linear 
;and angular velocities every 100 cycles
model cycle 1000 calm 100
;delete any balls that flew out of the bed
ball delete range position -1.5 -1.5 -2  1.5 1.5 -1 not

把材料转化成BPM

用pb_deform方法设置线性平行键接触模型的pb_kn和pb_ks属性；此时，球-球和球-面接触的平行结合强度也被分配。

球-面的平行结合强度被设置得非常大，以确保边界处的球保持与邻近区域的结合。

;bond the balls together with the specified gap
contact method bond gap 2.0e-3
;set the parallel bond properties from the pb_deform method
contact method pb_deform emod [materialMod] kratio 1.0
;remove the linear portion of the contact properties
;and set the contact strength
;this will be a relatively weak material
contact property lin_mode 1 pb_ten 2e4 pb_coh 2e5 ...
        range contact type 'ball-ball'
;set the ball-facet contact strength to be very high
contact property lin_mode 1 pb_ten 1e100 pb_coh 1e100 ...
        range contact type 'ball-facet'
;remove the box
wall delete wall range name 'top'
;set the default behavior for new ball-ball and ball-facet contacts
contact cmat default model linear type ball-ball method deformability ...
                   emod [materialMod] kratio 1.0 property fric 0.3
contact cmat default model linear type ball-facet method deformability ...
                   emod [materialMod] kratio 1.0 property fric 0.3
;remove all velocities, fix the balls, 
;cycle 2 to remove any contact forces that previously existed from 
;the simple packing procedure
model calm
ball fix velocity spin
ball attribute force-contact multiply 0.0
contact property lin_force 0 0 0
model cycle 2
ball free velocity spin

大规模失效前位移图

大规模失效后位移图

加载命令

;push the punch into the balls with just minimal failure of the material
model cycle 1000
zone gridpoint fix velocity-z 0.0 range group "fixed"
;solve to an equilibrium state to show the current displacements
model solve ratio-average 1.0e-4 
model save 'beforeFailure'

;push the punch in and watch the zone displacement
;diminish after the onset of failure
zone gridpoint fix velocity-z -.0001 range group "fixed"
model cycle 1000
zone gridpoint fix velocity-z 0.0 range group "fixed"
;failure will continue and the granular specimen will rearrange,
;reducing the displacements
model solve ratio-average 1.0e-3 
model save 'afterFailure'

Sleeved Triaxial Test of a Bonded Material

这个例子展示了使用基于壳的结构元素对粘合颗粒模型（BPM）进行一套套筒式三轴试验。这种性质的模型可以单独使用PFC进行研究。然而，在这种情况下，为了使套筒壁模仿弹性膜的行为，将需要相当数量的FISH。

在创建一个耦合的PFC-FLAC3D模型之前，必须加载PFC模块。PFC球是在FLAC3D外壳组成的圆柱形套筒内创建的。为了在每次运行数据文件时创建相同的球的初始云，模型随机命令被用来初始化随机种子到一个指定的值。由于PFC总是在大应变模式下运行，因此必须明确激活区域和结构元素的大应变模式（见模型最大应变命令）。

命令流

弧线被创建为几何边缘，并被垂直挤压以创建圆柱形套筒。

;define relevant parameters for the cylinder
[rad = 1.0]
[height = 3.0]
[segments = 6]
[halfLen = height/2.0]
[freeRegion = height/2.0*0.8]

;create the cylinder with the geometry logic
;first create arcs making a circle of edges 首先创建弧线，形成一个圆的边缘
geometry edge create by-arc origin (0,0,[-halfLen]) ...
   start ([rad*(-1)],0,[-halfLen]) end (0,[rad*(-1)],[-halfLen]) ...
   segments [segments]
geometry edge create by-arc origin (0,0,[-halfLen]) ...
   start (0,[rad*(-1)],[-halfLen]) end ([rad],0,[-halfLen]) ...
   segments [segments]
geometry edge create by-arc origin (0,0,[-halfLen]) ...
   start ([rad],0,[-halfLen]) end (0,[rad],[-halfLen]) ...
   segments [segments]
geometry edge create by-arc origin (0,0,[-halfLen]) ...
   start (0,[rad],[-halfLen]) end ([rad*(-1)],0,[-halfLen]) ...
   segments [segments]
;extrude the edges to make a cylinder
geometry generate from-edges extrude (0,0,[height]) segments [segments*2]

结构外壳导入命令用于从创建的几何体中创建外壳。元素和节点组都被分配，以简化边界条件分配。代码片段还显示了一个FISH函数，用来设置每个节点的本地轴系，使其指向圆柱体的中心，这样就可以使用结构壳应用命令对模型中间的壳元素施加压力。z方向直接指向圆柱体的中心，y方向指向上方。节点的速度和旋转是固定的，这样球就可以在套筒内得到平衡。

;import the structural elements as a shell from the geometry logic
;and assign groups and properties
structure shell import from-geometry 'Default' element-type dkt-cst 
structure node group 'middle' range position-z [-freeRegion] [freeRegion]
structure node group 'top' range position-z [freeRegion] [freeRegion+1]
structure node group 'bot' range position-z [-freeRegion-1] [-freeRegion]
structure shell group 'middle' range position-z [-freeRegion] [freeRegion]
structure shell property isotropic (1e6, 0.0) thick 0.25 density 930.0
;initialize the structural element related data structures 
model cycle 0
;In order to use the STRUCTURE SHELL APPLY command 
;set the local system of each structural
;element to be pointing to the center of the triaxial cell
fish define setLocalSystem
    loop foreach local s struct.node.list()
        local p = struct.node.pos(s)
        local nid = struct.node.id.component(s)
        local mvec = vector(0,0,comp.z(p))
        zdir = math.unit(p-mvec)
        ydir = vector(0,0,1)
        command
            structure node system-local z [zdir] y [ydir] ...
                      range component-id [nid]
        endcommand
    endloop
	command
	  structure node fix system-local
	endcommand
end
@setLocalSystem
;make sure that local damping is active for the structural elements
structure damp local
;fix the structural element nodes for specimen creation
structure node fix velocity rotation

为套筒和压板创建墙体

;create a wall representation of the structural element
wall-structure create
;parameter to set the platen width relative to the cylinder radius
[rad2 = rad*1.3]
;create the platens
wall generate name 'platenTop' polygon ([-rad2],[-rad2],[halfLen]) ...
                                       ([rad2],[-rad2],[halfLen])  ...
                                       ([rad2],[rad2],[halfLen])   ...
                                       ([-rad2],[rad2],[halfLen])
wall generate name 'platenBottom' polygon ([-rad2],[-rad2],[-halfLen]) ...
                                          ([rad2],[-rad2],[-halfLen])  ...
                                          ([rad2],[rad2],[-halfLen])   ...
                                          ([-rad2],[rad2],[-halfLen])
;set the wall resolution strategy to full and set the cutoff angle for 
;proximity contacts                                          
wall resolution full
wall attribute cutoff-angle 20

设置接触

;set the ball modulus and generate a cloud of balls with arbitrary overlap
[ballemod = 1.0e8]
ball distribute box [-rad] [rad] [-rad] [rad] [-halfLen] [halfLen] ...
                porosity 0.3 radius 0.05 0.1 ...
                range cylinder end-1 (0,0,[-halfLen]) end-2 (0,0,[height]) ...
                      radius [rad*0.95]

;set the ball attributes
ball attribute density 2600 damp 0.8
;set the default contact behavior - 
;the deformability method sets properties of the 
;linear portion of the contact model
contact cmat default model linearpbond method deformability ...
                   emod [ballemod] kratio 1.0 
;allow the balls to rearrange, 
;nulling the linear and angular velocities every 100 cycles
model cycle 1000 calm 100

;provide some friction to kill additional energy 
;and apply this to all of the contacts
contact cmat default model linearpbond method deformability ...
                   emod [ballemod] kratio 1.0 property fric 0.3
contact cmat apply
;solve to a small average ratio
model solve ratio-average 1e-8

;bond the ball-ball contacts
contact method bond gap 1.0e-2 range contact type 'ball-ball'
;change the existing contact properties 
;so that the linear force is incremental and
;supply strength
contact property lin_mode 1 lin_force 0 0 0 pb_ten 2e5 pb_coh 2e6
;set the stiffness of the parallel-bond portion of the contact model 
contact method pb_deformability emod [ballemod] kratio 1.0 ...
        range contact type 'ball-ball'

;this set of operations removes all ball velocities 
;and all contact forces in the model
;so that the specimen is completely stress free and bonded
model calm
ball fix velocity spin
model cycle 2
ball free velocity spin

;free the nodes in the middle section of the sleeve 
;while keeping the top and bottom edges fixed
structure node free velocity rotation range group 'middle'

监测应力应变

;define the stress FISH function to measure the stress and strain
fish define stress
    local topForce = math.abs(comp.z(wall.force.contact(platenTop))) 
    local botForce = math.abs(comp.z(wall.force.contact(platenBottom)))
    currentStress = 0.5*(topForce+botForce)/area
    stress = currentStress
    strain = (height - (comp.z(wall.pos(platenTop)) - ...
              comp.z(wall.pos(platenBottom))))/height * 100
    if failureStress <= currentStress
        failureStress = currentStress
        failureStrain = strain
    endif
end

设置停止条件

一旦当前应力水平低于最大应力的85%，循环将停止。一个FISH函数可以用来终止循环，用模型解决fish-halt命令。一旦函数halt的返回值为1，循环将终止。

;define the halt FISH function to stop cycling 
;as the platens displace and the material fails
fish define halt
    halt = 0
    if currentStress < failureStress * 0.85
        halt = 1
    endif
end

利用FISH函数逐步增加压板和套筒的约束应力。此外，还使用了内置的壁面伺服机制，通过壁面伺服命令来提供压板上的应力。

;define a FISH function to increase the isotropic confining pressure 
;in increments so that the bonded material does not break during loading
fish define rampUp(beginIn,ending,increment)
    command
        ball attribute displacement (0,0,0)
        structure node initialize displacement (0,0,0)
    endcommand
    begin = beginIn
    loop while (math.abs(begin) < math.abs(ending))
        begin = begin + increment
        command
            ;apply the confining stress
            structure shell apply [begin] range group 'middle'
            ;apply the same confining stress on the platens
            wall servo force (0,0,[begin*area]) activate true ...
                 range name 'platenTop'
            wall servo force (0,0,[-begin*area]) activate true ...
                 range name 'platenBottom'
model cycle 200
            model calm
        endcommand
    endloop
    command
        ;once the stress state has been installed cycle and turn off the servo 
        ;mechanism on the walls
model cycle 1000
        wall servo activate false
        wall attribute velocity (0,0,0) range name 'platenTop'
        wall attribute velocity (0,0,0) range name 'platenBottom'
    endcommand
end