PFC的一些整理资料.doc

PFC2D (Particle Follow Code 2 Dimension)即二维颗粒流程序，是通过离散单元方法来模拟圆形颗粒介质的运动及其相互作用。最初，这种方法是研究颗粒介质特性的一种工具，它采用数值方法将物体分为有代表性的数百个颗粒单元，期望利用这种局部的模拟结果来研究边值间题连续计算的本构模型。以下两种因素促使PFC2D方法产生变革与发展:(1)通过现场实验来得到颗粒介质本构模型相当困难:(2)随着微机功能的逐步增强，用颗粒模型模拟整个问题成为可能，一些本构特性可以在模型中自动形成。因此，PFC2D便成为用来模拟固体力学和颗粒流问题的一种有效手段。 2、 颗粒流方法的基本假设 ： 颗粒流方法在模拟过程中作了如下假设:1)颗粒单元为刚性体;2)接触发生在很小的范围内，即点接触;3)接触特性为柔性接触，接触处允许有一定的“重叠”量;4) “重叠”量的大小与接触力有关，与颗粒大小相比，“重叠”量很小;5)接触处有特殊的连接强度;6)颗粒单元为圆盘形(或球形)。3、颗粒流方法的特点： PFC2D可以直接模拟圆形颗粒的运动和相互作用问题。颗料可以代表材料中的个别颗粒，例如砂粒，也可以代表粘结在一起的固体材料，例如混凝土或岩石。当粘结以渐进的方式破坏时，它能够破裂。粘结在一起的集合体可以是各向同性，也可以被分成一些离散的区域或块体。这类物理系统可以用处理角状块体的离散单元程序UDEC和3DEC来模拟。PFC2D有三个优点： 第一、它有潜在的高效率。因为圆形物体间的接触探测比角状物体间的更简单。 第二、对可以模拟的位移大小实质上没有限制。 第三、由于它们是由粘结的粒子组成，块体可以破裂，不象UDEC和3DEC模拟的块体不能破裂。 用PFC2D模拟块体化系统的缺点是，块体的边界不是平的，用户必须接受不平的边界以换取PFC2D提供的优点。PFC2D能模拟任意大小圆形粒子集合体的动态力学行为。 粒子生成器根据粒子的指定分布规律自动概率地生成。粒子半径按均匀分布或按高斯分布规律分布。 初始孔隙度一般比较高，但通过控制粒子半径的扩大可以获得密度压实。在任何阶段任何因素都可以改变半径。所以不需反复试验就可以获得指定孔隙度的压实状态。 属性与各个粒子或接触有关，而不是与“类型号”有关。 因此，可以指定属性和半径的连续变化梯度。“节理生成器”用来修改沿指定轨迹线的接触特性。假定这些线叠加在颗粒集合体上。用这种方法，模型可以被成组的弱面，如岩石节理切割。 粒子颜色也是一种属性，用户可以指定各种标记方案。 PFC2D模型中为了保证数据长期不漂移，用双精度数据存储坐标和半径。接触的相对位移直接根据坐标而不是位移增量计算。接触性质由下列单元组成： 1）线性弹簧或简化的Hertz-Mindlin准则； 2）库仑滑块； 3）粘结类型：粘结接触可承受拉力，粘结存在有限的抗拉和抗剪强度。 可设定两种类型的粘结，接触粘结和平行粘结。这两种类型粘结对应两种可能的物理接触：接触粘结再现了作用在接触点一个很小区域上的附着作用；平行粘结再现了粒子接触后浇注其它材料的作用（如水泥灌浆）。平行粘结中附加材料的有效刚度具有接触点的刚度。 块体逻辑支持附属粒子组或块体的创建，促进了程序的推广普及。块体内粒子可以任意程度的重叠，作为刚性体具有可变形边界的每一个块体，可作为一般形状的超级粒子。通过指定墙的速度、混合的粒子速度、施加外力和重力来给系统加载。“扩展的FISH库”提供了在集合体内设置指定应力场或施加应力边界条件的函数。时步计算是自动的，包括因为Hertz接触模型刚度变化的影响。模拟过程中，根据每个粒子周围接触数目和瞬间刚度值，时步也在变化。基于估计的粒子数，单元映射策略采用最佳的单元数目，自动调整单元的外部尺寸来适应粒子缺失和指定的新对象。单元映射方案支持接触探测算法以保证求解时间随粒子数目线性增加，而不是二次方增加。 类似于FLAC，PFC提供了局部无粘性阻尼。这种阻尼形式有以下优点： 1）对于匀速运动，体力接近于零，只有加速运动时才有阻尼； 2）阻尼系数是无因次的； 3）因阻尼系数不随频率变化，集合体中具有不同自然周期的区域被同等阻尼，采用同样的阻尼系数。 PFC2D可以在半静态模式下运行以保证迅速收敛到静态解，或者在完全动态模式下运行。 PFC2D包含功能强大的内嵌式程序语言FISH，允许用户定义新的变量和函数使数值模型适合用户的特殊需求。例如，用户可以定义特殊材料的模型和性质、加载方式、实验条件的伺服控制、模拟的顺序以及绘图和打印用户定义的变量等。2. FISH语言简介 FISH是一种内置于Itasca软件内的编程语言，使用FISH 用户可以定义新的变量和函数，从而使得这些函数被用来扩展Itasca软件的用法或者增加用户自定义的特性。例如，用户可以绘制(PLOT)或打印(PRINT)新的变量，也能够改进特殊的颗粒体模型(网格)生成器，可以对数值试验进行伺服控制，可以设定一些性质的特殊分布以及自动化参数研究。 Itasca已经写了一些简单但非常有用的FISH函数作为库文件包含在各个具体的软件中，这些FISH函数一方面方便了一些没有编程经验的用户写一些简单的FISH函数，另一方面也能使用户在这些提供的简单函数的基础上作进一步的改进。不过，FISH语言象其它编程语言一样，也可以编制出非常复杂的程序 利用FISH语言进行编程，应该首先编一些简单的函数，然后仔细检查函数的功能，测试是否有错误。如果没有发现错误，再逐渐增加其功能，增加一项功能检查一下，直至发展到最后比较复杂的程序。这是因为虽然FISH是一种编译型语言，但它没有自己独立的编译器，不象VC+或VB能够实时全面地检查错误，FISH检查错误的能力很差，因此在使用他们到真实的应用之前，一定要用一些简单的数据（假如可能的话）来检查所有定义的函数。 FISH函数内置于标准的Itasca软件的数据文件中，函数的格式必须以DEFINE开始，以END结束。函 数可以嵌套调用，但定义函数的次序没有关系，只要在使用之前全部定义就行。由于FISH函数的编译格式储存在Itasca软件的内存中，因此可以用SAVE命令保存函数以及相关变量的当前值。 FISH也可以用来改进用户写的本构模型，如例1： DEF abc abc=22*3+5 END Print abc 对上例子稍作改进(例2)： new def abc hh=22 abc=hh*3+5 end1).稍有编程常识的人可以看出，执行上面的例子（PRINT abc），其结果与例1相同：abc=71.在这个函数中，我们首先把22赋值给变量hh，然后把这个变量带入abc的表达式中，因此二者的结果相同。2).FISH的执行过程如下：当在程序命令中使用一个FISH符号名时（例如执行PRINT 符号名），如果符号名也是一个函数名，那么执行这个函数（例如abc）；如果符号名不是函数名，那么使用符号目前的值(例如hh)。3).在输入完例2的各行后，如果我们执行命令：PRINT hh，此时hh=0,因为在这个时候没有执行FISH函数，因此hh的初始值为0；我们接着执行PRINT abc，结果显示abc=71;再次执行PRINT hh,此时结果为hh=22，这是因为我们首先运行了abc函数，在这个过程中hh已被赋值。4). 下面的试验将进一步解释函数与变量之间的差别。注意：Itasca软件的SET命令可以用来设置任何用户定义的FISH符号的值，与在FISH中使用的符号无关。下面的例3建立在例2的基础之上，我们不使用NEW命令来清除内存中的值，因为我们想继续使用那些值: set abc=0 hh=0 print hh print abc print hh5).在这个例子中，我们首先把abc和hh都赋值为0，由于hh是一个变量，第一个Print命令显示当前hh的值，hh=0;第二个Print命令由于abc是一个函数名，因此执行abc函数，先前定义的abc=0不起作用，重新计算了hh和abc的值，因此第三个Print命令显示的值是它在abc函数内指定的值，即hh=22,例4是这个试验完整的命令。 象其它高级编程语言一样，FISH有执行循环命令的功能，标准的格式如下： LOOP var(expr1,expr2) END_LOOP其中LOOP和END_LOOP是FISH语句，符号var代表循环变量，expr1和pxpr2代表表达式或者单个变量，下面的例子用循环命令计算从1到10的和以及乘积,见下例： new def abc sum = 0 prod = 1 loop n (1,10) sum = sum + n prod = prod * n end_loop end abc print sum, prod 在这个例子中，首先给两个变量赋于初始值，sum用来保存和的结果，prod用来保存积的结果，然后执行循环，最后分别打印出这两个变量的最后结果。循环变量n(1,10)表示从1开始，连续计算到10结束。关于LOOP的注意事项：1).FISH接受END_LOOP和ENDLOOP的写法，但不接受END LOOP这样中间有空格的写法，其它类似的命令有着同样的规则，如END_IF，END_COMMAND等命令；2).在上面的例子中，如果执行Print n或Print fish命令，你会看到n=11而不是10，注意：这不是FISH的错误，这是一个基本的计算机指令存储规则，当循环结束后，计数器的值保存的是n+1而不是n,所有的高级编程语言有着相同的规则。1.DEFINE function END2.CASEOF expr Case n endcase3. IF expr1 test expr2 THEN ELSE ENDIF4. LOOP var (expr1, expr2) ENDLOOP5. LOOP WHILE expr1 test expr2ENDLOOP6. COMMAND ENDCOMMAND7. HISTORY var PRINT var SET var value PLOT add .sh fname另外，在FISH中还有许多其它的预定义对象,其中一类是尺度变量(scalar variables),它们是单个的数字,下面是总的尺度变量: clock-时钟时间,单位是秒的100倍. unbal-最大不平衡力 pi-圆周率 step-目前的时步数目 urand-0.0-1.0之间均匀分布的随机变量 这仅是其中的一小部分,完全的列表以后再述 另一类非常有用的内置对象是固有函数(intrinsic functions),这些函数能在FISH内进行一些比较高级的数学运算,完整的列表见FISH 手册,下面给出其中的一部分: abs(a)-a的绝对值 cos(a)-a的余玄(a为弧度) log(a)-a的底数为10的对数 max(a,b)-返回a,b中的最大值 sqrt(a)-a的平方根3. PFC2D计算模型的生成方法 有两个命令可用于生成颗粒流模型:BALL和GENER-ATE，其中，BALL命令是生成单个的颗粒，该命令生成的颗粒可与已存在的颗粒重叠，而GENERATE 可生成一系列指定数目的颗粒流，该命令生成的颗粒是不允许重叠的。PFC2D里主要有两种类型的颗粒流：规则排列的和无规则排列的。一系列规则排列的颗粒流可以用来模拟模拟结构部分，如梁，而不规则排列的颗粒流可用来模拟实体或内部结构无规则的颗粒材料，如岩石内部所包含的胶结颗粒。 尽管颗粒的排列是随机的，但在颗粒模型生成后，整个模型的结构特性还是可能会受影响的，比如弱的结构面或各向异性。对于无规律排列的颗粒流模型，一般不可能去描述它的初始接触力的量级大小，这必须在后期要经过一个压缩的过程才可能给予较好的评价。3.1规律排列颗粒流： 生成规则排列的颗粒流，主要采用FISH语言配合BALL命令，循环生成一系列的颗粒，如下例: 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 (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 end set 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 yellowplot show3.2不规则排列 ； 无规则排列，即：对一个给定空隙率的区域，采用颗粒来充填其中需要进行填充的空隙，并确保整个模型保持平衡。对于所能被填充的模型的初始空隙率，是有一个限制值，不能任意小。对于某些空隙率的模型，颗粒的填充可以无接触地排列，对于其它情况的空隙率，颗粒又可以重叠排列目前没有一个普遍的方法来将紧密充填的圆形颗粒体限制在一个任意封闭的面积内。这里来介绍两种方法如何充填特定半径的颗粒体达到我们所需要的空隙率。 第一种方法，首先建立封闭区域的边界(简称墙体)，然后在封闭区域内任意生成一系列无接触的颗粒，最后移动区域的限制墙体，至所需要的空隙率。这种方法有三个缺点：1.区域的几何形状改变；2.收敛速度慢；3.最终的分布趋势是不均匀的第二种方法：运用GENERATE命令生成颗粒体，同时配合关键词高斯分配，即指定颗粒体半径的上下限，然后相应分配一个标准差，同时配合FISH函数来选择颗粒半径，最终生成我们所需要的模型。 下面来介绍采用这种方法是如何进行操作的，重点介绍2种方法。3.2.1 半径扩展法 首先需要介绍一个概念：半径放大系数m，它会引起模型空隙率的改变。空隙率的定义： 颗粒体的总面积； 封闭区域的总面积因此： 我们定义初始空隙率为 ，新生成的空隙率为 ,那么，则有实例：区域宽：10单位 区域高：5单位 目标空隙率：0.12 颗粒体数目：300 最大最小半径比：1.5 初始假设一个m之值，可以求出初始的 为从而得到平均半径大小为：由初始给定的颗粒体最大直径和最小直径的比例r可得到颗粒体半径的上下限：通过上面的计算，同时编写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 - (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) command 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 command 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_commandget_poros mult = sqrt(1.0 - poros) / (1.0 - pmeas) command ini rad mul 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.SAV3.2.2 挤压排斥法： 第2种方法：指定颗粒体的半径，不限制颗粒的数目，使足够多的颗粒产生来达到所需要的空隙率。但这种方法所 带来的缺点是可能在局部区域造成大面积的颗粒重叠，这将会产生很大的挤压力，从而给予颗粒较大的初始速度，这就可能使得颗粒体脱离墙体的限制。为避免此情况的发生，可通过初始的有限步循环计算将颗粒的动能减至零，然后再计算至平衡态。 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 wall 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_commandpvol_sum = 0.0tot_vol = width * heightcount = 0section 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 = id1 then if b_id(bp) = id2 then sum = sum + pi * b_rad(bp)2 end_if end_if bp = b_next(bp) end_loop pmeas = 1.0 - sum / tot_volenddef final_poros tot_vol = width * height id1 = 1 id2 = 1200 get_poros final_poros = pmeasend set x1=0.0 x2=5.0 y1=0.0 y2=5.0 id1=1 id2=50make_block set mult_a=mult set x1=5.0 x2=10.0 y1=0.0 y2=5.0 id1=1001 id2=1200make_block set mult_b=mult ini rad mul=mult_a c_index 0 range id 1,50 ini rad mul=mult_b c_index 1 range id 1001,1200 plo create the_assembly plot add ball lgreen lorange plot add wall black plot showcycle 1000 prop fric 0.2 cycle 500 print final_porossave expand2.SAV3.4 运用模型生成“过滤器” 有些情况需要我们建立复杂区域形状的颗粒流模型，如右图，此时我们可运用模型生 成过滤器来获得所需要的模型， 即filter命令，其后由用户定义 FISH函数来控制,其中,颗粒的半 径通过fc_arg(0)进行检验,x和y 的坐标位置分别通过fc_arg(1) 和fc_arg(2)进行检验。如果颗 粒满足要求，则FISH函数值设为 0，否则为1。详见下例。newdef ff_rect ; - 用户定义生成过滤器生成方形环状颗粒流模型 ; 中心 (ff_x, ff_y), 内径 ff_r1 and外径 ff_r2. _brad = fc_arg(0) _bx = fc_arg(1) _by = fc_arg(2) _skip = 0 _rx = abs( _bx - ff_x ) - _brad _ry = abs( _by - ff_y ) - _brad if _rx ff_r1 then if _ry ff_r1 then _skip = 1 end_if end_if ff_rect = _skipenddef gen_balls _xlo = ff_x - ff_r2 _xhi = ff_x + ff_r2 _ylo = ff_y - ff_r2 _yhi = ff_y + ff_r2 command generate x=(_xlo, _xhi) y=(_ylo, _yhi) & rad=(0.09, 0.11) & filter=ff_rect & id=(1,250) end_commandend set ff_x=1.0 ff_y=1.0 ff_r1=2.0 ff_r2=3.0gen_balls property dens=1000 kn=1e8 ks=1e8 wall id=1 nodes (-1.0,-1.0) (-1.0, 3.0) (3.0,3.0) (3.0,-1.0) close wall id=2 nodes (-2.0,-2.0) ( 4.0,-2.0) wall id=3 nodes ( 4.0,-2.0) ( 4.0, 4.0) wall id=4 nodes ( 4.0, 4.0) (-2.0, 4.0) wall id=5 nodes (-2.0, 4.0) (-2.0,-2.0) wall id=1 kn=1e8 ks=1e8 wall id=2 kn=1e8 ks=1e8 wall id=3 kn=1e8 ks=1e8 wall id=4 kn=1e8 ks=1e8 wall id=5 kn=1e8 ks=1e8 plot create the_view plot add ball yellow plot add axes black plot add wall blue id=on plot show pause property rad mul 1.5 plot add cf greencycle 5004. 边界条件： PFC2D中有三种边界条件，分别是：墙体边界、颗粒体边界和混合边界。其中颗粒体边界又分为速度边界和受力边界。 1).墙体边界 建模过程中，墙体可作为颗粒体的生成范围约束，但同时也可以将墙体作为边界来施加约束。对于墙体，我们只能施加速度约束，而不能直接对其施加外力，因为运动定律对墙体是不适用的。其速度由以下三个参数控制:线速度、角速度和旋转中心。墙的运动是通过不断更新定义墙的基点的位置来描述。采用WALL命令设置，如： wall id=1 x=1.0 y=1.0 spin=10.0(2)颗粒体边界PFC2D中的模型可以将一连串的颗粒体作为边界条件。基本方法：在模型紧密压缩至平衡后，我们通过FISH函数将与墙体相接触的颗粒体逐个提取，将这一系列的颗粒体采用共同的边界条件限制，最后删除初始的限制墙体，即实现了以颗粒体代替墙体来作为边界条件。颗粒体边界条件分速度边界和外力边界，下例为速度边界程序实现restore expand.sav def bound bp = ball_head loop while bp # null section cp = b_clist(bp) loop while cp # null if c_nforce(cp) # 0.0 then b2 = c_ball2(cp) if pointer_type( b2 ) = 101 then ; b2 is a wall b_xfix(bp) = 1 ; fix original ball in x,y b_yfix(bp) = 1 b_color(bp) = 1 exit section ; all done for this ball end_if end_if if c_ball1(cp) = bp cp = c_b1clist(cp) else cp = c_b2clist(cp) end_if end_loop end_section bp = b_next(bp) end_loopendboundini xvel=0.0 grad=(-0.1, 0)ini yvel=0.0del wall 1del wall 2del wall 3del wall 4plot create assemblyplot add ball red greenplot add vel blackplot showplot add cf whitecycle 200颗粒体外力边界： 其前期操作与速度边界类似，差异处是先将边界颗粒体的速度初始化为零，再删除限制墙体，采用FISH函数反向施加颗粒体所受的不平衡力，运行一定的计算步至平衡态。在此过程中，区域内部的颗粒体会因轻微的扰动而导致其自身产生运动而脱离边界颗粒体，此时前面将边界颗粒体的速度固定为零。ini xvel=0 yvel=0 spin=0 del wall 1 del wall 2 del wall 3 del wall 4 cycle 1def replace ; Replace unbalanced forces with applied forces bp = ball_head loop while bp # null if b_xfix(bp) = 1 b_xfap(bp) = -b_xfob(bp) b_yfap(bp) = -b_yfob(bp) end_if bp = b_next(bp) end_loopendreplacefree x y spincycle 200(3) 混合边界条件： 混合边界即速度边界与外力边界同时存在，对于双轴压缩试验，其两侧端面不允许存在明显的几何变形，否则其边界力将会无效。因此，此时需要用到混合边界，两侧端面采用外力边界来约束，确保其不发生几何变形，侧面Y向采用速度边界，使其速度线性增加，以此来模拟试块轴向变形。5. 初始条件 : 首先介绍一下应力的概念。在连续介质中，应力是一个连续量，应力一般定义为作用于单位面积上的力，对于二维问题，应力则定义为作用于单位线段上的合力。然而，对于松散介质这种定义是不合适的，因为颗粒介质是离散的，在颗粒模型中没有存在于每一个点上连续的应力，而且变化幅度很大。对于应变也存在这样的问题。在PFC模型中，计算颗粒间接触力和颗粒的位移，这些量对于在微观上研究材料的特性很有意义，但不能将这些量直接联系到连续模型中，需要经过一个平均的过程，才能从微观传到连续模型中。所以,在PFC2D中采用平均应力的概念来表示5.1 施加等向应力 def expand_radii_sum = 0.0 bp = ball_head loop while bp # null ;找到与bp相连的下一个颗粒的物理地址cp = b_clist(bp)loop while cp # nullif c_nforce(cp) # 0.0 then _rcp = sqrt(c_x(cp) - b_