1、,PFC及几种方法简介,201130410019,一、PFC简介,1、pfc背景,PFC系列软件是由ITASCA咨询集团（设有ITASCA中国公司）开发的颗粒流分析程序（Partical Flow Code)，分为PFC2D,PFC3D两种 特别用于模拟任意性状、大小的二维圆盘或三维球体集合体的运行及其相互作用的强大颗粒分析程序。除了模拟大体积流动和混合材料力学研究，程序更适合于描述在固体材料中细观/宏观裂纹扩展、破坏累积并断裂、破坏冲击和微震响应等高水平课题的深化研究。,PFC中的颗粒为刚性体，但在力学关系上允许重叠，以模拟颗粒之间的接触力。颗粒之间的力学关系非常简单，即牛顿第二定律。颗粒之
2、间的接触破坏可以为剪切和张开两种形式，当介质中颗粒间的接触关系（如断开）发生变化时，介质的宏观力学特性受到影响，随着发生破坏的接触数量增多，介质宏观力学特性可以经历从峰前线性到峰后非线性的转化，即介质内颗粒接触状态的变化决定了介质的本构关系。因此，在PFC计算中不需要给材料定义宏观本构关系和对应的参数，这些传统的力学特性和参数通过程序自动获得，而定义它们的是颗粒和水泥的几何和力学参数，如颗粒级配、刚度、摩擦力、粘结介质强度等微力学参数。,2、基本假设,1)颗粒单元为刚性体;2)接触发生在很小的范围内，即点接触;3)接触特性为柔性接触，接触处允许有一定的“重叠”量;4) “重叠”量的大小与接触力
3、有关，与颗粒大小相比，“ 重叠”量很小;5)接触处有特殊的连接强度;6)颗粒单元为圆盘形(或球形)。,3、优缺点,优点：第一、它有潜在的高效率。因为圆形物体间的接触探测比角状物体间的更简单。 第二、对可以模拟的位移大小实质上没有限制。 第三、由于它们是由粘结的粒子组成，块体可以破裂，不象UDEC和3DEC模拟的块体不能破裂。 缺点是：块体的边界不是平的，用户必须接受不平的边界以换取PFC2D提供的优点。,4、求解步骤,1)定义模拟对象 根据模拟意图定义模型的详细程序，假如只对某一力学机制的不同解释作出判断时，可以建立一个比较粗略的模型，只要在模型中能体现要解释的机制即可，对所模拟问题影响不大的
4、特性可以忽略。2)建立力学模型的基本概念首先对分析对象在一定初始特性形成初步概念。为此，应先提出一些问题，如系统是否将变为,不稳定系统、问题变形的大小、主要力学特性是否非线性、是否需要定义介质的不连续性、系统边界是实际边界还是无限边界、系统结构有无对称性等3)构造并运行简化模型 在建立实际工程模型之前，先构造并运行一系列简化的测试模型，可以提高解题效率。通过这种前期简化模型的运行，可对力学系统的概念有更深入的了解，有时在分析简化模型的结果后(例如所选的接触类型是否有代表性、边界条件对模型结果的影响程度等)，还需将第二步加以修改,4)补充模拟问题的数据资料 模拟实际工程问题需要大量简化模型运行的
5、结果，对于地质力学来说包括: a)几何特性，如地下开挖酮室的形状、地形地貌、坝体形状、岩土结构等； b)地质构造位置，如断层、节理、层面等； c)材料特性，如弹/塑性、后破坏特性等； d)初始条件，如原位应力状态、孔隙压力、饱和度等； e)外荷载，如冲击荷载、开挖应力等。,5)模拟运行的进一步准备 a)合理确定每一时步所需时间，若运行时间过长，很难得到有意义的结论，所以应该考虑在多台计算机上同时运行。 b)模型的运行状态应及时保存，以便在后续运行中调用其结果。例如如果分析中有多次加卸荷过程，要能方便地退回到每一过程，并改变参数后可以继续运行。 c)在程序中应设有足够的监控点(如参数变化 处、不
6、平衡等)，对中间模拟结果随时作出比较分析，并分析颗粒流动状态。,6)运行计算模型 在模型正式运行之前先运行一些检验模型，然后暂停，根据一些特性参数的试验或理论计算结果来检查模拟结果是否合理，当确定模型运行正确无误时，连接所有的数据文件进行计算。 7)解释结果 计算结果与实测结果进行分析比较。图形应集中反应要分析的区域如应力集中区，各种计算结果应能方便地输出分析。,二、 PFC2D计算模型的几种生成方法,1、有两个命令可用于生成颗粒流模型:BALL和GENER-ATE，其中，BALL命令是生成单个的颗粒，该命令生成的颗粒可与已存在的颗粒重叠，而GENERATE 可生成一系列指定数目的颗粒流，该命
7、令生成的颗粒是不允许重叠的。PFC2D里主要有两种类型的颗粒流：规则排列的和无规则排列的。尽管颗粒的排列是随机的，但在颗粒模型生成后，整个模型的结构特性还是可能会受影响的，比如弱的结构面或各向异性。对于无规律排列的颗粒流模型，一般不可能去描述它的初始接触力的量级大小，这必须在后期要经过一个压缩的过程才可能给予较好的评价。,1.1规律排列颗粒流,New def hex xc = x0 yc = y0 rc = radius idc = id_start r2 = 2.0 * radius yinc = radius * sqrt(3.0) loop row (1,n_row) loop col
8、(1,n_col) command ball id=idc x=xc y=yc rad=rc end_command idc = idc + 1 xc = xc + r2 end_loop yc = yc + yinc xc = x0 + radius * (row - (row/2) * 2),end_loop endset echo offset x0=0.2 y0=0.4 radius=0.1set id_start=100 n_col=7 n_row=8hexset echo onplot set cap size 20plot add axes blackplot add ball
9、yellowplot show,1.2不规则排列无规则排列，即：对一个给定空隙率的区域，采用颗粒来充填其中需要进行填充的空隙，并确保整个模型保持平衡。对于所能被填充的模型的初始空隙率，是有一个限制值，不能任意小。对于某些空隙率的模型，颗粒的填充可以无接触地排列，对于其它情况的空隙率，颗粒又可以重叠排列。 第一种方法，首先建立封闭区域的边界(简称墙体)，然后在封闭区域内任意生成一系列无接触的颗粒，最后移动区域的限制墙体，至所需要的空隙率。这种方法有三个缺点：1.区域的几何形状改变；2.收敛速度慢；3.最终的分布趋势是不均匀的,第二种方法：运用GENERATE命令生成颗粒体，同时配合关键词高斯分配
10、，即指定颗粒体半径的上下限，然后相应分配一个标准差，同时配合FISH函数来选择颗粒半径，最终生成我们所需要的模型。,半径扩展法,newdef expand ;- 输入数据 - n_stiff = 1e8 ; 法向连接刚度 s_stiff = 1e8 ; 剪切连接刚度 width = 10.0 ; 区域宽 height = 5.0 ; 区域高 tot_vol = width*height poros = 0.12 ; 最终目标空隙率 num = 300 ; 颗粒体数目 rat = 1.5 ; 最大最小半径比 ;- 导出所需数据 - mult = 1.6 ; 初始半径放大系数 n0 = 1.0 -
11、 (1.0 - poros) / mult2 r0 = sqrt(height*width*(1.0 - n0)/(pi*num) rlo = 2.0 * r0 / (1.0 + rat) rhi = rat * rlo _x1 = width*(1.0 + extend) _y1 = 0.0,command wall id=1 ks=s_stiff kn=n_stiff nodes (_x0,_y0) (_x1,_y1) end_command _x0 = width _y0 = -extend*height _x1 = width _y1 = height*(1.0 + extend) c
12、ommand wall id=2 ks=s_stiff kn=n_stiff nodes (_x0,_y0) (_x1,_y1) end_command _x0 = width*(1.0 + extend) _y0 = height _x1 = -extend*width _y1 = height command wall id=3 ks=s_stiff kn=n_stiff nodes (_x0,_y0) (_x1,_y1) end_command _x0 = 0.0 _y0 = height*(1.0 + extend) _x1 = 0.0 _y1 = -extend*height,com
13、mand wall id=4 ks=s_stiff kn=n_stiff nodes (_x0,_y0) (_x1,_y1) end_command ;- generate the balls and give them their properties command gen id=1,num rad=rlo,rhi x=0,width y=0,height prop dens=1000 ks=s_stiff kn=n_stiff end_command get_poros mult = sqrt(1.0 - poros) / (1.0 - pmeas) command ini rad mu
14、l mult cycle 1000 prop fric 0.2 cycle 2000 end_command enddef get_poros sum = 0.0 bp = ball_head loop while bp # null,sum = sum + pi * b_rad(bp)2 bp = b_next(bp) end_loop pmeas = 1.0 - sum / (width * height) end expand get_poros plot wall ball plot show print pmeas save expand.SAV,挤压排斥法 指定颗粒体的半径，不限制
15、颗粒的数目，使足够多的颗粒产生来达到所需要的空隙率。但这种方法所 带来的缺点是可能在局部区域造成大面积的颗粒重叠，这将会产生很大的挤压力，从而给予颗粒较大的初始速度，这就可能使得颗粒体脱离墙体的限制。为避免此情况的发生，可通过初始的有限步循环计算将颗粒的动能减至零，然后再计算至平衡态。,new set random ; 设定颗粒体数目无限制模式def explode n_stiff = 1e8 s_stiff = 1e8 width = 10.0 height = 5.0 poros = 0.12 n_max = 1000 rlo = 0.174 rhi = 0.261 command wal
16、l id 1 ks=s_stiff kn=n_stiff nodes (0,0) (width,0) wall id 2 ks=s_stiff kn=n_stiff nodes (width,0) (width,height) wall id 3 ks=s_stiff kn=n_stiff nodes (width, height) (0,height) wall id 4 ks=s_stiff kn=n_stiff nodes (0,height) (0,0) end_command pvol_sum = 0.0tot_vol = width * heightcount = 0,sectio
17、n loop n (1,n_max) r_ball = rlo + urand * (rhi - rlo) pvol_new = pvol_sum + pi * r_ball2 if (1.0 - pvol_new / tot_vol) poros then exit section ; (new porosity will be that specified) end_if rb2 = r_ball * 2.0 x_ball = r_ball + urand * (width - rb2) y_ball = r_ball + urand * (height - rb2) command ba
18、ll x=x_ball y=y_ball rad=r_ball end_command pvol_sum = pvol_new count = count + 1 end_loop end_section tot_poros = 1.0 - pvol_sum / tot_vol command prop dens=1000.0 ks=s_stiff kn=n_stiff end_command ii = out( +string(count)+ balls were created.) end,def quiet loop n (1,4) command cycle 5 ini xvel 0
19、yvel 0 end_command end_loop end explode quiet cyc 980 prop fric=0.2 cyc 250 plot wall ball print tot_poros save explode.SAV,三、边界条件,PFC2D中有三种边界条件，分别是：墙体边界、颗粒体边界和混合边界。其中颗粒体边界又分为速度边界和受力边界。 1).墙体边界 建模过程中，墙体可作为颗粒体的生成范围约束，但同时也可以将墙体作为边界来施加约束。对于墙体，我们只能施加速度约束，而不能直接对其施加外力，因为运动定律对墙体是不适用的。其速度由以下三个参数控制:线速度、角速度和旋
20、转中心。墙的运动是通过不断更新定义墙的基点的位置来描述。采用WALL命令设置，如： wall id=1 x=1.0 y=1.0 spin=10.0,(2)颗粒体边界 PFC2D中的模型可以将一连串的颗粒体作为边界条件。基本方法：在模型紧密压缩至平衡后，我们通过FISH函数将与墙体相接触的颗粒体逐个提取，将这一系列的颗粒体采用共同的边界条件限制，最后删除初始的限制墙体，即实现了以颗粒体代替墙体来作为边界条件。,颗粒体外力边界,其前期操作与速度边界类似，差异处是先将边界颗粒体的速度初始化为零，再删除限制墙体，采用FISH函数反向施加颗粒体所受的不平衡力，运行一定的计算步至平衡态。在此过程中，区域内
21、部的颗粒体会因轻微的扰动而导致其自身产生运动而脱离边界颗粒体，此时前面将边界颗粒体的速度固定为零。,混合边界条件,混合边界即速度边界与外力边界同时存在，对于双轴压缩试验，其两侧端面不允许存在明显的几何变形，否则其边界力将会无效。因此，此时需要用到混合边界，两侧端面采用外力边界来约束，确保其不发生几何变形，侧面Y向采用速度边界，使其速度线性增加，以此来模拟试块轴向变形。,四、本构模型,在PFC2D中，材料的本构特性是通过接触本构模型来模拟的。每一颗粒的接触本构模型有:1)接触刚度模型;2)滑动模型;3)连接模型。接触刚度模型提供了接触力和相对位移的弹性关系，滑动模型则强调切向和法向接触力使得接触
22、颗粒可以发生相对移动，而连接模型是限制总的切向和法向力使得在连接强度范围内发生接触。,五、赋予材料属性,PFC2D程序里，除了颗粒体联结相关的属性(联结刚度、强度以及平行联结半径)外，其它颗粒体的属性赋予均采用PROPERTY 命令来执行。 同时，在执行命令时，对与不同区域不同位置的颗粒体，其属性的变化可通过关键词 gradient和group来配合使用，达到我们所需设置的要求和变化。,六、节理面的生成及属性设置,PFC2D中，节理面可以用来模拟颗粒组间的滑移和分离，还可以模拟节理、断层和层理等，同时还能来模拟两种不同材料间的接触面，如基础和土体。 我们采用JSET命令来设置，可生成单独的一条
23、或一组节理面。但只有存在法向力的颗粒间的接触，才能被设置为节理面接触，因此，节理面命令必须在颗粒体模型生成后并达到平衡后才能执行。,如若生成一条45度角的单节理面，我们可以采用下面的命令： JSET id=1 dip=-45.0 origin=(0.0,2.0),七、加载,PFC2D荷载中的类型一般分为两类：主动荷载和被动荷载。被动荷载如重力，采用SET GRAVITY 设置。主动荷载则包括施加颗粒体速度和外力。颗粒体可以直接施加速度和外力，但外力却不能直接施加在墙体上，因为运动方程对于墙体是无效的，但可以通过其他等效的方法来模拟墙体受集中力或应力。,加载的操作主要从这四方面来执行：1).墙体
24、的控制；2).颗粒体外力的控制；3).颗粒体速度控制；4).混合受力的控制。 上述几方面，其具体编程方法与前面所介绍的边界条件与初始条件类似，在这里就不详细介绍，具体的现实问题或工程所需的特殊的加载方式，需采用FISH来进行编写和操作，同时需从简单的模型来验证其实际控制效果。,八、求解,PFC2D中所谓的求解即为执行许多个计算循环步(使用CYCLE和STEP)，同时观察计算后的模型响应。对于静力问题，模型系统要么计算达到平衡，要么破坏坍塌。对于动力问题，就显得较复杂些，特别需要注意如何消除个别已脱离模型颗粒体。虽然程序会对这些已脱离模型的颗粒体持续计算，但他们的变化对整个模型的影响会越来越小，因为只有当颗粒体处于计算的约束区域内才会产生其应有的效应。一种简易的消除逃逸颗粒体的方法是每隔几千步计算就在一定范围内执行DELETE BALL 命令。,1).静力求解,对于求解命令CYCLE或STEP，用户也可以用SOLVE命令来进行静力求解。该命令是给定某个不平衡力准则，当达到这个标准时，SOLVE命令即停止。同时可配合一些关键词共同执行。如Average，允许指定平均不平衡力与平均接触力的比值；maximum则允许指定最大不平衡力与最大接触力的比值， 默认的值为0.01(1%)。 solve average=0.005 maximum=1e-10,
