网上经典-38pfc汇总

1、PFC2D (Particle Follow Code 2 Dimen)即二维颗粒流程序，是通过离散单元方法来模拟圆形颗粒介质的运动及其相互作用。最初，这种方法是颗粒介质特性的一种工具，它采用数值方法将物体分为有代表性的数百个颗粒单元，期望利用这种局部的模拟结果来边值间题连续计算的本构模型。以下两种促使 PFC2D 方法产生与发展:(1)通过现场实验来得到颗粒介质本构模型相当:(2)随着微机功能的逐步增强，用颗粒模型模拟整个问题成为可能，一些本构特性可以在模型中自动形成。因此，PFC2D 便成为用来模拟固体力学和颗粒流问题的一种有效。2、颗粒流方法的基本假设 ：颗粒流方法在模拟过程中作了如下假

2、设: 1)颗粒单元为刚性体;2)接触发生在很小的范围内，即点接触;3)接触特性为柔性接触，接触处允许有一定的“”量;4) “”量的大小与接触力有关，与颗粒大小相比，“”量很小;接触处有特殊的连接强度;颗粒单元为圆盘形(或球形)。3、颗粒流方法的特点：PFC2D 可以直接模拟圆形颗粒的运动和相互作用问题。颗料可以代表材料中的个别颗粒，例如砂粒，也可以代表粘结在一起的固体材料，例如混凝土或岩石。当粘结以渐进的方式破坏时，它能够破裂。粘结在一起的集合体可以是各向同性，也可以被分成一些离散的区域或块体。这类物理系统可以用处理角状块体的离散单元程序 UDEC 和 3DEC 来模拟。 PFC2D 有三个优

3、点：第一、它有潜在的高效率。因为圆形物体间的接触探测比角状物体间的更简单。第二、对可以模拟的位移大小实质上没有限制。第三、由于它们是由粘结的粒子组成，块体可以破裂，不象 UDEC 和 3DEC 模拟的块体不能破裂。用 PFC2D 模拟块体化系统的缺点是，块体的边界不是平的，用户必须接受不平的边界以换取 PFC2D 提供的优点。PFC2D 能模拟任意大小圆形粒子集合体的动态力学行为。粒子分布规律分布。根据粒子的指定分布规律自动概率地生成。粒子半径按均匀分布或按初始孔隙度一般比较高，但通过控制粒子半径的扩大可以获得密度压实。在任何阶段任何都可以改变半径。所以不需反复试验就可以获得指定孔隙度的压实状

4、态。属性与各个粒子或接触有关，而不是与“类型号”有关。因此，可以指定属性和半径的连续变化梯度。“节理”用来修改沿指定轨迹线的接触特性。假定这些线叠加在颗粒集合体上。用这种方法，模型可以被成组的弱面，如岩石节理切割。粒子颜色也是一种属性，用户可以指定各种标记方案。PFC2D 模型中为了保证数据长期不漂移，精度数据坐标和半径。接触的相对位移直接根据坐标而不是位移增量计算。接触性质由下列单元组成：1）线性弹簧或简化的 Hertz-Mindlin 准则；2）滑块；3）粘结类型：粘结接触可承受拉力，粘结存在有限的抗拉和抗剪强度。可设定两种类型的粘结，接触粘结和平行粘结。这两种类型粘结对应两种可能的物理接

5、触：接触粘结再现了作用在接触点一个很小区域上的附着作用；平行粘结再现了粒子接触后浇注其它材料的作用（刚度。泥灌浆）。平行粘结中附加材料的有效刚度具有接触点的块体逻辑支持附属粒子组或块体的创建，促进了程序的推广普及。块体内粒子可以任意程度的，作为刚性体具有可变形边界的每一个块体，可作为一般形状的超级粒子。通过指定墙的速度、混合的粒子速度、施加外力和重力来给系统加载。“扩展的 FISH 库”提供了在集合体内设置指定应力场或施加应力边界条件的函数。时步计算是自动的，包括因为 Hertz接触模型刚度变化的影响。模拟过程中，根据每个粒子周围接触数目和瞬间刚度值，时步也在变化。基于估计的粒子数，单元来适应

6、粒子缺失和指定的新对象。单元目线性增加，而不是二次方增加。策略采用最佳的单元数目，自动调整单元的外部尺寸方案支持接触探测算法以保证求解时间随粒子数类似于 FLAC，PFC 提供了局部无粘性阻尼。这种阻尼形式有以下优点：对于匀速运动，体力接近于零，只有阻尼系数是无因次的；运动时才有阻尼；3）因阻尼系数不随频率变化，集合体中具有不同自然周期的区域被同等阻尼，采用同样的阻尼系数。PFC2D 可以在半静态模式下运行以保证迅速收敛到静态解，或者在完全动态模式下运行。PFC2D 包含功能强大的内嵌式程序语言 FISH，允许用户定义新的变量和函数使数值模型适合用户的特殊需求。例如，用户可以定义特殊材料的模型

7、和性质、加载方式、实验条件的伺服控制、模拟的顺序以及绘图和打印用户定义的变量等。2.FISH 语言简介FISH 是一种内置于 Itasca内的编程语言，使用 FISH用户可以定义新的变量和函数，从而使得这些函数被用来扩展 Itasca的用法或者增加用户自定义的特性。例如，用户可以绘制(PLOT)或打印(PR)新的变量，也能够改进特殊的颗粒体模型(网格)，可以对数值试验进行伺服控制，可以设定一些性质的特殊分布以及自动化参数。Itasca 已经写了一些简单但非常有用的 FISH 函数作为库文件包含在各个具体的中，这些 FISH 函数一方面方便了一些没有编程经验的用户写一些简单的 FISH 函数，另

8、一方面也能使用户在这些提供的简单函数的基础上作进一步的改进。不过，FISH 语言象其它编程语言一样，也可以编制出非常复杂的程序利用 FISH 语言进行编程，应该首先编一些简单的函数，然后仔细检查函数的功能，测试是否有错误。如果没有发现错误，再逐渐增加其功能，增加一项功能检查一下，直至发展到最后比较复杂的程序。这是因为虽然 FISH 是一种编译型语言，但它没有自己独立的编译器，不象 VC+或 VB 能够实时全面地检查错误，FISH 检查错误的能力很差，因此在使用他们到真实的应用之前，一定要用一些简单的数据（假如可能的话）来检查所有定义的函数。FISH 函数内置于标准的 Itasca的数据文件中，

9、函数的格式必须以 DEFINE 开始，以END 结束。函 数可以嵌套调用，但定义函数的次序没有关系，只要在使用之前全部定义就行。由于 FISH 函数的编译格式数以及相关变量的当前值。在 Itasca的内存中，因此可以用 SAVE 命令保存函FISH 也可以用来改进用户写的本构模型，如例 1：DEF abcabc=22*3+5ENDPr对上例子稍作改进(例 2)：abcnewdef abchh=22abc=hh*3+5end的人可以看出，执行上面的例子（PR1).稍有编程abc），其结果与例 1 相同：abc=71.在这个函数中，首先把 22 赋值给变量 hh，然后把这个变量带入 abc 的表达

10、式中，因此二者的结果相同。2).FISH 的执行过程如下：当在程序命令中使用一个 FISH 符号名时（例如执行 PR符号名），如果符号名也是一个函数名，那么执行这个函数（例如 abc）；如果符号名不是函数名，那么使用符号目前的值(例如 hh)。3).在输入完例 2 的各行后，如果执行命令：PRhh，此时 hh=0,因为在这个时候没有执行 FISH 函数，因此 hh 的初始值为 0；接着执行 PRabc，结果显示 abc=71;再次执行 PR被赋值。hh,此时结果为 hh=22，这是因为首先运行了 abc 函数，在这个过程中 hh 已4). 下面的试验将进一步解释函数与变量之间的差别。注意：It

11、asca的 SET 命令可以用来设置任何用户定义的 FISH 符号的值，与在 FISH 中使用的符号无关。下面的例 3 建立在例 2 的基础之上，不使用 NEW 命令来清除内存中的值，因为想继续使用那些值:set abc=0 hh=0pr prprhh abchh5).在这个例子中，首先把 abc 和 hh 都赋值为 0，由于 hh 是一个变量，第一个 Pr命令显示当前 hh 的值，hh=0;第二个 Pr命令由于 abc 是一个函数名，因此执行 abc 函数，先前定义的 abc=0 不起作用，重新计算了 hh 和 abc 的值，因此第三个 Pr命令显示的值是它在 abc 函数内指定的值，即 h

12、h=22,例 4 是这个试验完整令。象其它高级编程语言一样，FISH 有执行循环命令的功能，标准的格式如下：LOOPvar(expr1,expr2)END_LOOP其中 LOOP 和 END_LOOP 是 FISH 语句，符号 var 代表循环变量，expr1 和 pxpr2 代表表达式或者单个变量，下面的例子用循环命令计算从 1 到 10 的和以及乘积,见下例：new def abcsum= 0prod = 1loop n (1,10)sum= sum + n prod = prod * nend_loopend abcprsum, prod在这个例子中，首先给两个变量赋于初始值，sum 用

13、来保存和的结果，prod 用来保存积的结果，然后执行循环，最后分别打印出这两个变量的最后结果。循环变量 n(1,10)表示从 1 开始，连续计算到 10 结束。关于 LOOP 的注意事项：1).FISH 接受 END_LOOP 和 ENDLOOP 的写法，但不接受 END LOOP 这样中间有空格的写法，其它类似令有着同样的规则，如 END_IF，MAND 等命令；2).在上面的例子中，如果执行 Prn 或 Prfish 命令，你会看到 n=11 而不是 10，注意：这不是 FISH 的错误，这是一个基本的计算机指令规则，当循环结束后，计数器的值保存的是 n+1 而不是 n,所有的高级编程语言

14、有着相同的规则。 1.DEFINEfunctionENDCASEOFexpr Case nendcaseIFexpr1 test expr2 THEN ELSEENDIF4. LOOPvar (expr1, expr2)ENDLOOP5. LOOP WHILE expr1 test expr2 ENDLOOP6.DMAND7. HISTORY varPRvarSET var value PLOT add .sh fname另外，在 FISH 中还有许多其它的预定义对象,其中一类是尺度变量(scalar variables),它们是单个的数字,下面是总的尺度变量:clock时钟时间,unbal最

15、大不平衡力pi圆周率是秒的 100 倍.step目前的时步数目urand0.0-1.0 之间均匀分布的随量这仅是其中的一小部分,完全的列表以后再述另一类非常有用的内置对象是固有函数(rinsic 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 命令

16、是生成单个的颗粒，该命令生成的颗粒可与已存在的颗粒，而 GENERATE 可生成一系列指定数目的颗粒流，该命令生成的颗粒是不允许的。PFC2D 里主要有两种类型的颗粒流：规则排列的和无规则排列的。一系列规则排列的颗粒流可以用来模拟模拟结构部分，如梁，而不规则排列的颗粒流可用来模拟实体或内部结构无规则的颗粒材料，如岩石内部所包含的胶结颗粒。尽管颗粒的排列是随机的，但在颗粒模型生成后，整个模型的结构特性还是可能会受影响的，比如弱的结构面或各向异性。对于无规律排列的颗粒流模型，一般不可能去描述它的初始接触力的量级大小，这必须在后期要经过一个压缩的过程才可能给予较好的评价。3.1 规律排列颗粒流：生成

17、规则排列的颗粒流，主要采用 FISH 语言配合 BALL 命令，循环生成一系列的颗粒，如下例:Newdef hexxc = x0 yc = y0rc = radius idc = id_startr2 = 2.0 * radiusyinc = radius * sqrt(3.0) loop row (1,n_row)loop col (1,n_col)dball id=idc x=xc y=yc rad=rc mandidc = idc + 1 xc = xc + r2 end_loopyc = yc + yincxc = x0 + radius * (row - (row/2) * 2) e

18、nd_loopendset echo off？set x0=0.2 y0=0.4 radius=0.1set id_start=100 n_col=7 n_row=8 hexset echo onplot set cap size 20 plot add axes black plot add ball yellow plot show3.2 不规则排列；无规则排列，即：对一个给定空隙率的区域，采用颗粒来充填其中需要进行填充的空隙，并确保整个模型保持平衡。对于所能被填充的模型的初始空隙率，是有一个限制值，不能任意小。对于某些空隙率的模型，颗粒的填充可以无接触地排列，对于其它情况的空隙率，颗粒又

19、可以排列目前没有一个普遍的方法来将紧密充填的圆形颗粒体限制在一个任意封闭的面积内。这里来介绍两种方法如何充填特定半径的颗粒体达到所需要的空隙率。第法，首先建立封闭区域的边界(简称墙体)，然后在封闭区域内任意生成一系列无接触的颗粒，最后移动区域的限制墙体，至所需要的空隙率。这种方法有三个缺点：1.区域的几何形状改变；2.收敛速度慢；3.最终的分布趋势是不均匀的第二种方法：运用 GENERATE 命令生成颗粒体，同时配合分配，即指定颗粒体半径的上下限，然后相应分配一个标准差，同时配合 FISH 函数来选择颗粒半径，最终生成所需要的模型。下面来介绍采用这种方法是如何进行操作的，重点介绍 2 种方法。

20、3.2.1 半径扩展法首先需要介绍一个概念：半径放大系数 m，它会引起模型空隙率的改变。空隙率的定义：颗粒体的总面积； 封闭区域的总面积因此：定义初始空隙率为，新生成的空隙率为,那么，则有实例：区域宽：10区域高：5目标空隙率：0.12颗粒体数目：300最大最小半径比：1.5初始假设一个 m 之值，可以求出初始的为从而得到平均半径大小为：由初始给定的颗粒体最大直径和最小直径的比例 r到颗粒体半径的上下限：通过上面的计算，同时编写 FISH 函数来生成满足需要的空隙率的模型，并能自动计算空隙率的大小，以便用户进行检验。newdef expand;- 输入数据 - n_stiff = 1e8 s_

21、stiff = 1e8 width = 10.0height = 5.0; 法向连接刚度; 剪切连接刚度; 区域宽; 区域高tot_vol = width*heightporos = 0.12num = 300rat = 1.5导出所需数据mult = 1.6; 最终目标空隙率; 颗粒体数目; 最大最小半径比-; 初始半径放大系数;-n0 = 1.0 - (1.0 - poros) / mult2r0 = sqrt(height*width*(1.0 - n0)/(pi*num) rlo = 2.0 * r0 / (1.0 + rat)rhi = rat * rlo_x1 = width*(1

22、.0 + extend)_y1 = 0.0dwall id=1 ks=s_stiff kn=n_stiff nodes (_x0,_y0) (_x1,_y1) mand_x0 = width_y0 = -extend*height_x1 = width_y1 = height*(1.0 + extend) dwall id=2 ks=s_stiff kn=n_stiff nodes (_x0,_y0) (_x1,_y1) mand_x0 = width*(1.0 + extend)_y0 = height_x1 = -extend*width_y1 = height dwall id=3 ks

23、=s_stiff kn=n_stiff nodes (_x0,_y0) (_x1,_y1)mand_x0 = 0.0_y0 = height*(1.0 + extend)_x1 = 0.0_y1 = -extend*height dwall id=4 ks=s_stiff kn=n_stiff nodes (_x0,_y0) (_x1,_y1) mand;- generate the balls and give them their properties dgen id=1,num rad=rlo,rhi x=0,width y=0,height prop dens=1000 ks=s_st

24、iff kn=n_stiffmandget_porosmult = sqrt(1.0 - poros) / (1.0 - pmeas) dini rad mul mult cycle 1000prop fric 0.2cycle 2000mandenddef get_porossum = 0.0bp = ball_headloop whip # nullsum = sum + pi * b_rad(bp)2bp = b_next(bp) end_looppmeas = 1.0 - sum / (width * height)end expand get_poros plot wall ball

25、plot showprpmeassave expand.SAV3.2.2 挤压排斥法：第 2 种方法：指定颗粒体的半径，不限制颗粒的数目，使足够多的颗粒产生来达到所需要的空隙率。但这种方法所 带来的缺点是可能在局部区域造成大面积的颗粒，这将会产生很大的挤压力，从而给予颗粒较大的初始速度，这就可能使得颗粒体脱离墙体的限制。为避免此情况的发生，可通过初始的有限步循环计算将颗粒的动能减至零，然后再计算至平衡态。newset random ; 设定颗粒体数目def exploden_stiff = 1e8 s_stiff = 1e8 width = 10.0height = 5.0poros = 0.

26、12n_max = 1000rlo = 0.174rhi = 0.261 d模式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)mand pvol_sum = 0.0tot_vo

27、l = width * height count = 0sectionloop n (1,n_max)r_ball = rlo + urand * (rhi - rlo) pvol_new = pvol_sum + pi * r_ball2if (1.0 - pvol_new / tot_vol) poros thenexit section ; (new porosity will be = id1 thenif b_id(bp) = id2 thensum = sum + pi * b_rad(bp)2 end_ifend_ifbp = b_next(bp) end_looppmeas =

28、 1.0 - sum / tot_volenddef final_porostot_vol = width * height id1 = 1id2 = 1200get_poros final_poros = pmeasendset x1=0.0 x2=5.0 y1=0.0 y2=5.0 id1=1 id2=50make_blockset mult_a=multset x1=5.0 x2=10.0 y1=0.0 y2=5.0 id1=1001 id2=1200make_blockset mult_b=multini rad mul=mult_a c_index 0 range id 1,50in

29、i rad mul=mult_b c_index 1 range id 1001,1200 plo create the_assemblyplot add ball lgreen lorange plot add wall blackplot show cycle 1000prop fric 0.2cycle 500prfinal_porossave expand2.SAV3.4 运用模型生成“过滤器”有些情况需要建立复杂区域形状的颗粒流模型，如右图，此时可运用模型生成过滤器来获得所需要的模型，即 filter 命令，其后由用户定义 FISH 函数来控制,其中,颗粒的半径通过 fc_arg(0

30、)进行检验,x 和y的坐标位置分别通过 fc_arg(1) 和 fc_arg(2)进行检验。如果颗粒满足要求，则 FISH 函数值设为0，否则为 1。详见下例。newde_rect;用户定义生成过滤器生成方形环状颗粒流模型; 中心 (ff_x, ff_y), 内径 ff_r1 and 外径_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 thenif _ry ff_r1 then_s

31、kip = 1 end_ifend_ifff_rect = _skipff_r2.enddef gen_balls_xlo = ff_x - ff_r2_xhi = ff_x + ff_r2_ylo = ff_y - ff_r2_yhi = ff_y + ff_r2dgenerate x=(_xlo, _xhi) y=(_ylo, _yhi) & rad=(0.09, 0.11) &filter=ff_rect & id=(1,250)mandendset ff_x=1.0 ff_y=1.0 ff_r1=2.0 ff_r2=3.0gen_ballsproperty dens=1000 kn=1

32、e8 ks=1e8wall id=1 nodes (-1.0,-1.0) (-1.0, 3.0) (3.0,3.0) (3.0,-1.0) closewall 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=1e8wall id=2 kn=1e8 ks=1e8 wall id=3 kn=1e8 ks=1e

33、8wall id=4 kn=1e8 ks=1e8wall id=5 kn=1e8 ks=1e8 plot create the_viewplot add ball yellow plot add axes blackplot add wall blue id=on plot showpauseproperty rad mul 1.5 plot add cf green cycle 5004.边界条件：PFC2D 中有三种边界条件，分别是：墙体边界、颗粒体边界和混合边界。其中颗粒体边界又分为速度边界和受力边界。1).墙体边界建模过程中，墙体可作为颗粒体的生成范围约束，但同时也可以将墙体作为边界来

34、施加约束。对于墙体，只能施度约束，而不能直接对其施加外力，因为运动定律对墙体是不适用的。其速度由以下三个参数控制:线速度、角速度和旋转中心。墙的运动是通过不断更新定义墙的基点的位置来描述。采用 WALL 命令设置，如：wall id=1 x=1.0 y=1.0 spin=10.0(2)颗粒体边界PFC2D 中的模型可以将一连串的颗粒体作为边界条件。基本方法：在模型紧密压缩至平衡后，通过 FISH 函数将与墙体相接触的颗粒体逐个提取，将这一系列的颗粒体采用共同的边界条件限制，最后删除初始的限制墙体，即实现了以颗粒体代替墙体来作为边界条件。颗粒体边界条件分速度边界和外力边界，下例为速度边界程序实现

35、restore expand.sav def boundbp = ball_headloop whisectionp # nullcp = b_clist(bp) loop while cp # nullif c_nforce(cp) # 0.0 then b2 = c_ball2(cp)if poer_type( b2 ) = 101 then ; b2 is a wall b_xfix(bp) = 1 ; fix original ball in x,y b_yfix(bp) = 1b_color(bp) = 1exit section ; all done for this ball e

36、nd_ifend_ifif c_ball1(cp) = bpcp = c_b1clist(cp) elsecp = c_b2clist(cp) end_ifend_loop end_sectionbp = b_next(bp) end_loopend boundini xvel=0.0 grad=(-0.1, 0) ini yvel=0.0del wall 1del wall 2del wall 3del wall 4plot create assembly plot add ball red green plot add vel blackplot showplot add cf white

37、 cycle 200颗粒体外力边界：其前期操作与速度边界类似，差异处是先将边界颗粒体的速度初始化为零，再删除限制墙体，采用 FISH 函数反向施加颗粒体所受的不平衡力，运行一定的计算步至平衡态。在此过程中，区域内部的颗粒体会因轻微的扰动而导致其自身产生运动而脱离边界颗粒体，此时前面将边界颗粒体的速度固定为零。ini xvel=0 yvel=0 spin=0 del wall 1del wall 2del wall 3del wall 4cycle 1def replace ; Replace unbalanced forbp = ball_headwippd forloop whip # nu

38、llif b_xfix(bp) = 1b_xfap(bp) = -b_xfob(bp) b_yfap(bp) = -b_yfob(bp)end_ifbp = b_next(bp)end_loopend replacefree x y spincycle 200(3)混合边界条件：混合边界即速度边界与外力边界同时存在，对于双轴压缩试验，其两侧端面不允许存在明显的几何变形，否则其边界力将会无效。因此，此时需要用到混合边界，两侧端面采用外力边界来约束，确保其不发生几何变形，侧面 Y 向采用速度边界，使其速度线性增加，以此来模拟试块轴向变形。5.初始条件 : 首先介绍一下应力的概念。在连续介质中，应力

39、是续量，应力一般定义线段上的合力。然而，为作用于面积上的力，对于二维问题，应力则定义为作用于对于松散介质这种定义是不合适的，因为颗粒介质是离散的，在颗粒模型中没有存在于每一个点上连续的应力，而且变化幅度很大。对于应变也存在这样的问题。在 PFC 模型中，计算颗粒间接触力和颗粒的位移，这些量对于在微观上材料的特性很有意义，但不能将这些量直接联系到连续模型中，需要经过一个平均的过程，才能从微观传到连续模型中。所以,在 PFC2D 中采用平均应力的概念来表示5.1 施加等向应力def expand_radii_sum = 0.0bp= ball_headloop whip # null;找到与 bp

40、 相连的下一个颗粒的物理地址cp = b_clist(bp) loop while cp # nullif c_nforce(cp) # 0.0 then_rcp = sqrt(c_x(cp) - b_x(bp)2 + (c_y(cp) - b_y(bp)2) if _rcp = b_rad(bp) then; 刚好接触或if c_ball1(cp) = bp then bp_other = c_ball2(cp)elsebp_other = c_ball1(cp) end_ifif poer_type(bp_other) = 101 then;21, 100, 101, 102, 103 o

41、r 104_phi = b_rad(bp)else; memory, ball,wall,contact,clump,circle_phi = b_rad(bp) + b_rad(bp_other) end_if_sum = _sum + _rcp * c_kn(cp) * _phi end_ifend_ifif c_ball1(cp) = bp thencp = c_b1clist(cp) elsecp = c_b2clist(cp) end_ifend_loopbp = b_next(bp) end_loop_alpha = -1.0 * _lambda * tot_vol * _diso

42、 / _sumbp = ball_headloop whip # nullb_rad(bp) = (1.0 + _alpha)*b_rad(bp) bp = b_next(bp)end_loop enddef achieve_isostrloop while 1 # 0;始终执行，除非遇到 exit 命令_diso = req_isostr - isostrif abs(_diso/req_isostr) = req_isostr_tol then exitend_ifoo = out(Current isotropic stress = + string(isostr) expand_rad

43、iidsolvemandend_loopenddef isostrii= measure(mp,1); 计算应力 1 表示应力，2 表示应变isostr = (m_s11(mp) + m_s22(mp) / 2.0enddef get_mp；返回 measure id=1 mp = find_meas(1)endmeas id=1 x 5 y 2.5 rad 2.0 get_mphistory id=1 isostrSET req_isostr = -4e5req_isostr_tol = 0.01 cyc 20prprmeas 1isostrachieve_isostr6.接触本构模型在 P

44、FC2D 中，材料的本构特性是通过接触本构模型来模拟的。每一颗粒的接触本构模型有:1)接触刚度模型;2)滑动模型;3)连接模型。接触刚度模型提供了接触力和相对位移的弹性关系，滑动模型则强调切向和法向接触力使得接触颗粒可以发生相对移动，而连接模型是限制总的切向和法向力使得在连接强度范围内发生接触。1).接触刚度模型接触刚度通过以下两式把接触力与相对位移联系起来，即:式中，是法向刚度，它为割线刚度，联系总的法向力与位移。式中，是切向刚度，它为切线刚度，通过增量形式联系切向力与位移。在 PFC2D 中有两种接触刚度模型，即线性接触刚度模型与简化的 Hertz-Minlin 接触刚度模型,不同的接触刚

45、度模型有不同的接触刚度值。.线性接触模型:分析颗粒间接触特性时，将这种接触的颗粒想象为一点在颗粒中心弹性梁，受力或力矩相当于作用于颗粒的中心。这种以下特征参数描述:1)几何参数:长度(L)、断面积(A)、惯性矩(I);2)变形参数:杨氏模量 E、泊松比 v ; 3)强度参数:法向强度、切向强度 ;颗粒间接触杨氏模量用 表示，平行连接杨氏模量表示。同时假定 PFC2D 中所有颗粒均为厚度为 t 的圆盘。若颗粒 A 和 B 接触，则梁的半径为:梁的长度为:式中:和分别为接触颗粒的半径。相互接触颗粒之间的受力特性，相当于弹性致，此时梁断面各 A 与惯性矩 I 为:受纯轴向或纯切向荷载时的情况一，式中

46、:t 为假设的颗粒圆盘厚度。关于接触变形，不仅指接触连接模型的接触连接变形，同样适用于颗粒之间以及颗粒与墙之间的接触变形，所以，下面介绍的计算方法即可用于接触连接材料也可用于非接触连接材料(砂土)。对于纯轴向荷载或纯切向荷载，其法向接触刚度与切向接触刚度分别为：；式中， 是接触杨氏模量，不同于材料整体杨氏模量(通常大于整体杨氏模量)。对于线性接触模型，其接触刚度 k是假设相互接触两球为串联：式中，分别代表切向和法向刚度。若两球有相同的刚度，即：此时，通过公式转换，在接触连接时，接触模量与球的刚度间的关系式:所以在描述变形微观参数时，先给定颗粒一颗粒接触的接触模量 EC 以及颗粒法向刚度与切向刚

47、度比值 kn/ks。根据式 1.计算kn。再据给定比值求 ks。线性接触模型是通过两个接触体(颗粒颗粒，颗粒一墙体)的切向刚度 ks 和法向刚度kn 定义的。计算线性接触模型的接触刚度时假定接触体为串联，根据式2.计算。.Hertz-Mindlin 接触模型Hertz-Mindlin 接触模型是 Mindlin&Deresiewich (1953)和 Cundall(1988)理论近似得出的一种非线性接触模型。它只能严格用于颗粒体接触，而且不能再现剪切过程中的连续非线性(特别使用初始剪切模量时)。Hertz-Mindlin 模型由两个接触球体的参数剪切模量 G 和泊松比 v 来定义。该模型不适

48、用于采用接触连接的两球体的接触，因为这种模型没有定义球体受张力的情况。对于颗粒颗粒相互接触，其弹性常数为平均值，对于颗粒墙体接触，假定墙为刚性，所以直接取球的弹性常数。接触法向割线刚度和切向切线刚度为:式中 Un 为颗粒体接触量，Fin 为法向接触力其它的系数为两接触体的几何与材料特性的函数。对于颗粒一颗粒相互接触，这些系数可表示为:对于颗粒墙相互接触，这些系数可表示为:式中，G 为弹性剪切模量，v 为泊松比，R 为颗粒的半径，A, B为相互接触的两球体。2).滑动模型滑动模型是相互接触球体的一种固有特性，它没有法向抗拉强度，允许颗粒在抗剪强度范围内发生滑动，这种模型在接触连接模型发生作用之前

49、一直有效。同时，滑动模型也可与半行连接模型同时起作用。滑动模型是通过两接触体间最小摩擦系数 u 定义的，若颗粒间则令法向和切向接触力等于零。发生滑动的判别条件为:量小于或等于零，若,则可以发生滑动，其中：3).连接模型：PFC2D 模型允许相互接触颗粒连接在一起，有两种连接模型:即接触连接与平行连接模型。接触连接认为连接只发生在接触点很小范围内，而平行连接发生在接触颗粒间圆形或方形有限范围内。接触连接只能传递力，而平行连接同时能传递力和力矩。两种类型的接触可以同时存在，直到超过接触强度。连接模型只能是颗粒之间的连接，颗粒与墙之间的连接不能采用连接模型。.接触连接模型接触连接可以想象为一对有恒定

50、法向刚度与切向刚度的弹簧作用于颗粒接触点处，同时，这些弹簧设定一定的抗拉与抗剪强度。只要接触连接存在就没有颗粒间的滑动发生。接触连接在颗粒间量小于零时，允许出现张力，但法向接触张力过接触连接强度。在 PFC2D 中，接触连接是有法向连接强度 Fcn 定义。当法向抗拉接触力大于或等于法向接触连接强度时，连接破坏并把法向、切向接触力赋值为零。当切向接触力大于或等于切向连接强度时，连接也发生破坏，但是接触力不发生变化，并假设切向力没有超过摩擦极限。颗粒接触点处接触力与相对位移的关系的本构特性见下图。接触力的法向分量接触力的切向分量.平行连接模型平行连接模型可以描述颗粒之间有限范围内有夹层材料的木构特

51、性。相互连接的两个颗粒可以看作是球体或柱体。这种连接建立颗粒间一种弹性相互作用关系，可与前面所述的滑动模型或接触连接模型同时存在。平行连接可以想象为一组有恒定法向刚度与切向刚度的弹簧均匀分布于接触平面内，这些弹簧作用的本构关系类似与点接触弹簧模拟颗粒刚度的本构特性。接触的相对运动在平行连接处产生力和力矩，作用于相互连接的颗粒上，且与连接材料的最向、切向应力有关。平行连接模型是由法向刚度、切向刚度、法向强度、切向强度和连接半径五个参数定义的。与平行连接相对应的总接触力和力矩用和表示，根据，总的接触力和力矩表示平行连接在颗粒 B 上的作用，如下图所示。可以将总的接触力沿接触面分解为切向分量和法向分

52、量：式中，表示法向分量， 表示切向分量。7.赋予材料属性 ： PFC2D 程序里，除了颗粒体联结相关的属性(联结刚度、强度以及平行联结半径)外，其它颗粒体的属性赋予均采用 PROPERTY 命令来执行。同时，在执行命令时，对与不同区域不同位置的颗粒体，其属性的变化可通过gradient 和group 来配合使用，达到所需设置的要求和变化。如颗粒体间接触剪切强度 随着 X 方向线形变化，其表达式为：可以这样来编写命令流文件：prop s_bond 1e6 grad -1e4, 0 range x 0 50具体求解为：x=15m 处，s_bond=1e6+(-1e4)*15同时还可将某个区域内的颗

53、粒体均定义为一个组，用 group 来执行，然后再给这个组赋予属性，具体操作为：group strong_rock range x 0 10 y 5 10group weak_rock range x 0 10 y 0 5 prop fric 1.0 range group strong_rock prop fric 0.1 range group weak_rock如果颗粒体的材料属性是一个非线形变化的量，对于颗粒间摩擦系数表达式为：就需要采用 FISH 语言来赋予属性,如其中 d 为颗粒体距模型顶面距离，h 为模型区域的宽度和高度,此时赋予属性时, FISH语言可这样来表述：def var

54、y_friction bp = ball_headloop whip # nully_= b_y(bp)depth = height - abs(y_)b_fric(bp) = fric0 * (1.0 + (depth/height)2) bp = b_next(bp)end_loopendSET height=5 fric0=0.2vary_friction8.节理面的生成及属性设置 ： PFC2D理面可以用来模拟颗粒组间的滑移和分离，还可以模拟节理、断层和层理等，同时还能来模拟两种不同材料间的接触面，如基础和土体。采用 JSET 命令来设置，可生成单独的一条或一组节理面。但只有存在法向力的颗粒间的接触，才能被设置为节理面接触，因此，节理面命令必须在颗粒体模型生成后并达到平衡后才能执行。dip 和 origin可配合使用来生成所需位置的节理面，dip 是指节理面的倾角，origin是指节理面所经过的某一点的坐标。 如若生成一条 45 度角的单节理面， JSET id=1 dip=-45.0 origin=(0.0,2.0)可以采用下面令：number 和对于多组节理面，可采用spacing 来控制。JSET id=2 dip=0 origin=(3.0,2.0) & number=2 spacing=5.0JSET i