PFC培训ppt课件.ppt

PFCPFC2D2D颗粒流软件培训颗粒流软件培训1PFC培训主要内容主要内容第一部分第一部分 PFC2D颗粒流程序粒流程序简介介 第二部分第二部分 有限差分法基有限差分法基础介介绍第三部分第三部分 离散元法基离散元法基础介介绍第四部分第四部分￿ ￿PFC2D的使用的使用2PFC培训第一部分第一部分 PFCPFC2D2D颗粒流程序简介颗粒流程序简介ü1 1 、理论背景ü2 2、颗粒流方法的基本假设ü3 3、颗粒流方法的特点ü4 4、可选特性ü5 5、应用领域ü6 6、求解步骤3PFC培训 作为离散元的一种，二维颗粒流程序(Particle (Particle Follow Code PFC2D)Follow Code PFC2D)数值模拟新技术，其理论基础是Cundall [1979]Cundall [1979]提出的离散单元法，用于颗粒材料力学性态分析，如颗粒团粒体的稳定、变形及本构关系，专门用于模拟固体力学大变形问题它通过圆形( (或异型) )离散单元来模拟颗粒介质的运动及其相互作用由平面内的平动和转动运动方程来确定每一时刻颗粒的位置和速度作为研究颗粒介质特性的一种工具，它采用有代表性的数百个至上万个颗粒单元，通过数值模拟实验可以得到颗粒介质本构模型。

1.1.理论背景4PFC培训 PFC2D (Particle Follow Code 2 Dimension)PFC2D (Particle Follow Code 2 Dimension)即二维颗粒流程序，是通过离散单元方法来模拟圆形颗粒介质的运动及其相互作用最初，这种方法是研究颗粒介质特性的一种工具，它采用数值方法将物体分为有代表性的数百个颗粒单元，期望利用这种局部的模拟结果来研究边值间题连续计算的本构模型以下两种因素促使PFC2DPFC2D方法产生变革与发展:(1):(1)通过现场实验来得到颗粒介质本构模型相当困难:(2):(2)随着微机功能的逐步增强，用颗粒模型模拟整个问题成为可能，一些本构特性可以在模型中自动形成因此，PFC2DPFC2D便成为用来模拟固体力学和颗粒流问题的一种有效手段 5PFC培训2 2、颗粒流方法的基本假设、颗粒流方法的基本假设 颗粒流方法在模拟过程中作了如下假设: :1)1)颗粒单元为刚性体; ;2)2)接触发生在很小的范围内，即点接触; ;3)3)接触特性为柔性接触，接触处允许有一定的““重叠””量; ;4) “4) “重叠””量的大小与接触力有关，与颗粒大小相比，““重叠””量很小; ;5)5)接触处有特殊的连接强度; ;6)6)颗粒单元为圆盘形( (或球形) )。

6PFC培训 其中，颗粒为刚性体的假设，对于模拟介质运动为只沿相互接触面的表面发生的问题非常重要，比如象砂土或粮食这种颗粒组合体材料，利用这种假设在总体上来讲是比较恰当的，因为这种材料的变形是来自于颗粒刚性体间的滑动和转动以及接触面处的张开和闭锁，而不是来自于每个刚性颗粒本身的变形，对于这种特殊材料，没有必要采用非常精确的数值模型，来得到对材料特性的近似7PFC培训3 3、颗粒流方法的特点、颗粒流方法的特点 PFC2DPFC2D可以直接模拟圆形颗粒的运动和相互作用问题颗料可以代表材料中的个别颗粒，例如砂粒，也可以代表粘结在一起的固体材料，例如混凝土或岩石当粘结以渐进的方式破坏时，它能够破裂粘结在一起的集合体可以是各向同性，也可以被分成一些离散的区域或块体这类物理系统可以用处理角状块体的离散单元程序UDECUDEC和3DEC3DEC来模拟8PFC培训 PFC2D PFC2D有三个优点： 第一、它有潜在的高效率因为圆形物体间的接触探测比角状物体间的更简单 第二、对可以模拟的位移大小实质上没有限制。

第三、由于它们是由粘结的粒子组成，块体可以破裂，不象UDECUDEC和3DEC3DEC模拟的块体不能破裂 用PFC2DPFC2D模拟块体化系统的缺点是，块体的边界不是平的，用户必须接受不平的边界以换取PFC2DPFC2D提供的优点9PFC培训 PFC2D PFC2D中几何特征、物理特性和解题条件的说明不如FLACFLAC和UDECUDEC程序那样直截了当 例如用连续介质程序，创建网格、设置初始压力、设置固定或自由边界在象PFC2DPFC2D这样的颗粒程序中，由于没有唯一的方法在一个指定的空间内组合大量的粒子，粒子紧密结合的状态一般不能预先指定必须跟踪类似于物体压实的过程，直到获得要求的孔隙率 由于颗粒相对位置变化产生接触力，初始应力状态的确定与初始压密有关由于边界不是由平面组成，边界条件的设定比连续介质程序更复杂10PFC培训 当要求满足有实验室实际测试的模拟物体的力学特性时，出现了更大的困难在某种程度上，这是一个反复试验的过程，因为目前还没有完善的理论可以根据微观特性来预见宏观特性。

然而，给出一些准则应该有助于模型与原型的匹配，如哪些因素对力学行为的某些方面产生影响，哪些将不产生影响应该意识到，由于受现有知识的限制，这样的模拟很难然而，用PFC2DPFC2D进行试验，对固体力学，特别是对断裂力学和损伤力学，可以获得一些基本认识 11PFC培训 PFC2D PFC2D能模拟任意大小圆形粒子集合体的动态力学行为 粒子生成器根据粒子的指定分布规律自动概率地生成粒子半径按均匀分布或按高斯分布规律分布 初始孔隙度一般比较高，但通过控制粒子半径的扩大可以获得密度压实在任何阶段任何因素都可以改变半径所以不需反复试验就可以获得指定孔隙度的压实状态￿12PFC培训 属性与各个粒子或接触有关，而不是与属性与各个粒子或接触有关，而不是与“类型号类型号”有关 因此，可以指定属性和半径的连续变化梯度因此，可以指定属性和半径的连续变化梯度节理生成器节理生成器”用来修改沿指定轨迹线的接触特性假用来修改沿指定轨迹线的接触特性假定这些线叠加在颗粒集合体上用这种方法，模型可定这些线叠加在颗粒集合体上。

用这种方法，模型可以被成组的弱面，如岩石节理切割以被成组的弱面，如岩石节理切割 粒子颜色也是一种属性，用户可以指定各种标记方粒子颜色也是一种属性，用户可以指定各种标记方案13PFC培训 PFC2DPFC2D模型中为了保证数据长期不漂移，用双精度数据存储坐标和半径接触的相对位移直接根据坐标而不是位移增量计算接触性质由下列单元组成： 1 1）线性弹簧或简化的Hertz-MindlinHertz-Mindlin准则； 2 2）库仑滑块； 3 3）粘结类型：粘结接触可承受拉力，粘结存在有限的抗拉和抗剪强度 可设定两种类型的粘结，接触粘结和平行粘结这两种类型粘结对应两种可能的物理接触：①①接触粘结再现了作用在接触点一个很小区域上的附着作用；②②平行粘结再现了粒子接触后浇注其它材料的作用（如水泥灌浆）平行粘结中附加材料的有效刚度具有接触点的刚度 14PFC培训 块体逻辑支持附属粒子组或块体的创建，促进了程序的推广普及块体内粒子可以任意程度的重叠，作为刚性体具有可变形边界的每一个块体，可作为一般形状的超级粒子。

通过指定墙的速度、混合的粒子速度、施加外力和重力来给系统加载扩展的FISHFISH库””提供了在集合体内设置指定应力场或施加应力边界条件的函数时步计算是自动的，包括因为HertzHertz接触模型刚度变化的影响模拟过程中，根据每个粒子周围接触数目和瞬间刚度值，时步也在变化基于估计的粒子数，单元映射策略采用最佳的单元数目，自动调整单元的外部尺寸来适应粒子缺失和指定的新对象单元映射方案支持接触探测算法以保证求解时间随粒子数目线性增加，而不是二次方增加15PFC培训 类似于FLACFLAC，PFCPFC提供了局部无粘性阻尼这种阻尼形式有以下优点： 1 1）对于匀速运动，体力接近于零，只有加速运动时才有阻尼； 2 2）阻尼系数是无因次的； 3 3）因阻尼系数不随频率变化，集合体中具有不同自然周期的区域被同等阻尼，采用同样的阻尼系数 PFC2DPFC2D可以在半静态模式下运行以保证迅速收敛到静态解，或者在完全动态模式下运行 PFC2DPFC2D包含功能强大的内嵌式程序语言FISHFISH，允许用户定义新的变量和函数使数值模型适合用户的特殊需求。

例如，用户可以定义特殊材料的模型和性质、加载方式、实验条件的伺服控制、模拟的顺序以及绘图和打印用户定义的变量等￿16PFC培训4 4、可选特性、可选特性 1 1）热学分析2 2）并行处理技术3 3）能写用户定义接触模型4 4）用户写C++C++程序的C++C++编程 17PFC培训 热学选项用来模拟材料内热量的瞬间流动和热诱导位移和力的顺序发展热学模型可以独立运行或耦合到力学模型通过修改粒子半径和平行粘结承受的力，产生热应变来解释粒子和粘结材料的受热 用户定义的接触本构模型可以用C++C++语言来编写，并编译成动态链接库文件，一旦需要就可以加载 用户写的C++C++程序选项允许用户用C++C++语言写自己的程序，创建可执行的PFC2DPFC2D个人版本这个选项可以用来代替FISHFISH函数，大大提高运行的速度 并行处理技术允许将一个PFC2DPFC2D模型分成几个部分，每个部分可以在单独的处理器上平行运行与一个PFC2DPFC2D模型在一个处理器上运行相比，平行处理在内存容量和计算速度方面得到大大提高18PFC培训5 5、应用领域、应用领域 PFC2D PFC2D既可解决静态问题也可解决动态问题，既可用于参数预测，也可用于在原始资料详细情况下的实际模拟。

PFC2D PFC2D 模拟试验可以代替室内试验在岩石与土体中开挖问题的研究与设计方面，实测资料相对较少，关于初始应力、不连续性等问题也只能部分了解而在松散介质流动问题中，影响流动介质不规律分布的影响因素很难定量描述因此，应用PFC2D PFC2D 初步研究影响整个系统的一些参数的特性，对整个系统的特性有所了解后，就可以方便地设计模型模拟整个过程 19PFC培训 PFC2DPFC2D可以模拟颗粒间的相互作用问题、大变形问题、断裂问题等，适用于以下领域：  （1 1）在槽、管、料斗、筒仓中松散物体的流动问题；   （2 2）矿山冒落法开采中的岩体断裂、坍塌、破碎和岩块的流动问题；   （3 3）铸模中粉料的压实问题；   （4 4）由粘结粒子组成物体的碰撞及其动态破坏；   （5 5）梁结构的地震响应及垮塌；   （6 6）颗粒材料的基本特性研究，如屈服、流动、体积变化等；   （7 7）固体的基本特性研究，如累积破坏、断裂￿20PFC培训6 6、求解步骤、求解步骤 1) 1)定义模拟对象 根据模拟意图定义模型的详细程序，假如只对某一力学机制的不同解释作出判断时，可以建立一个比较粗略的模型，只要在模型中能体现要解释的机制即可，对所模拟问题影响不大的特性可以忽略。

21PFC培训2)2)建立力学模型的基本概念建立力学模型的基本概念 首先对分析对象在一定初始特性形成初步概念首先对分析对象在一定初始特性形成初步概念为此，应先提出一些问题，如系统是否将变为不稳为此，应先提出一些问题，如系统是否将变为不稳定系统、问题变形的大小、主要力学特性是否非线定系统、问题变形的大小、主要力学特性是否非线性、是否需要定义介质的不连续性、系统边界是实性、是否需要定义介质的不连续性、系统边界是实际边界还是无限边界、系统结构有无对称性等际边界还是无限边界、系统结构有无对称性等 综合以上内容来描述模型的大致特征，包括颗综合以上内容来描述模型的大致特征，包括颗粒单元的设计、接触类型的选择、边界条件的确定粒单元的设计、接触类型的选择、边界条件的确定以及初始平衡状态的分析以及初始平衡状态的分析22PFC培训 3) 3)构造并运行简化模型 在建立实际工程模型之前，先构造并运行一系列简化的测试模型，可以提高解题效率通过这种前期简化模型的运行，可对力学系统的概念有更深入的了解，有时在分析简化模型的结果后( (例如所选的接触类型是否有代表性、边界条件对模型结果的影响程度等) )，还需将第二步加以修改。

23PFC培训 4) 4)补充模拟问题的数据资料 模拟实际工程问题需要大量简化模型运行的结果，对于地质力学来说包括: : a) a)几何特性，如地下开挖酮室的形状、地形地貌、坝体形状、岩土结构等； b)b)地质构造位置，如断层、节理、层面等； c)c)材料特性，如弹/ /塑性、后破坏特性等； d)d)初始条件，如原位应力状态、孔隙压力、饱和度等； e)e)外荷载，如冲击荷载、开挖应力等 因为一些实际工程性质的不确定性( (特别是应力状态、变形和强度特性) )，所以必须选择合理的参数研究范围第三步简化模型的运行有助于这项选择，从而为更进一步的试验提供资料24PFC培训 5) 5)模拟运行的进一步准备 a)a)合理确定每一时步所需时间，若运行时间过长，很难得到有意义的结论，所以应该考虑在多台计算机上同时运行 b)b)模型的运行状态应及时保存，以便在后续运行中调用其结果例如如果分析中有多次加卸荷过程，要能方便地退回到每一过程，并改变参数后可以继续运行。

c)c)在程序中应设有足够的监控点( (如参数变化 处、不平衡等) )，对中间模拟结果随时作出比较分析，并分析颗粒流动状态25PFC培训 6) 6)运行计算模型 在模型正式运行之前先运行一些检验模型，然后暂停，根据一些特性参数的试验或理论计算结果来检查模拟结果是否合理，当确定模型运行正确无误时，连接所有的数据文件进行计算 7)7)解释结果 计算结果与实测结果进行分析比较图形应集中反应要分析的区域如应力集中区，各种计算结果应能方便地输出分析26PFC培训第二部分第二部分 有限差分法基础介绍有限差分法基础介绍 连续介质三维快速拉格朗日有限差分计算方法( ( FLACFLAC3D3D) ) 是近2020年来逐步成熟完善起来的一种新型数值计算方法, ,它基于显式差分法来求解运动方程和动力方程, ,可模拟岩土或其他材料的三维力学行为其求解时首先将计算区域离散化, ,分成若干三维单元, ,单元之间由节点联结, ,节点受荷载作用后, ,其平衡方程( (运动方程) ) 可以写成时间步长为Δt Δt 的有限差分形式, ,由于采用动态应力松弛显式差分求解技术, , 在某一微小的时段内, , 作用于该节点的荷载只对周围若干节点有影响。

27PFC培训 根据单元节点的速度变化和时段Δt ,Δt ,可求出单元之间的相对位移, ,进而求出单元应变, ,利用单元材料的本构关系即可求出单元应力在此基础上, ,求出单元之间的不平衡力, ,将此不平衡力重新作用到节点上, ,再进行下一步的迭代过程, ,直到整个系统不平衡力足够小或节点位移趋于平衡为止 FLACFLAC3D3D可以解决诸多的有限元程序难以模拟的复杂的工程问题，例如分布开挖、大变形、非线性及非稳定系统（甚至大面积屈服/ /失稳或完全塌方）28PFC培训第三部分第三部分 离散元法基础介绍离散元法基础介绍 离散单元法是一种模拟非连续介质的计算方法，自CundallCundall在7070年代提出以来，在岩石力学、土力学、结构分析等领域的数值模拟中得到广泛应用，是一种新兴的非连续体分析方法离散单元法允许单元间的相对运动，不一定满足位移连续和变形协调条件，计算速度快，所需存储空间小，特别适用于节理岩体的大位移，大变形分析 离散单元法自问世以来有了长足的发展，已经成为解决岩石力学问题的一种重要的数值方法，因为工程中所见到的岩体其形态呈非连续结构，所形成的岩石块体运动和受力情况多是几乎或材料非线性问题，所以很难用解决连续介质力学问题的有限单元法或边界单元法等。

29PFC培训 数值方法来进行求解，而离散单元法正是充分考虑到岩体结构的不连续性，适用于解决节理岩石力学问题 近年来，离散元法的应用领域又扩展到求解连续介质向非连续介质转化的力学问题混凝土等脆性材料在冲击、侵彻等动荷载作用下产生的损伤和破坏，其实质是力学模型从连续体到非连续体的转变过程建立在传统的连续介质力学基础上的有限元法等数值计算方法难以直接用于计算和模拟材料具体的破坏形式和破坏的整个过程，而离散元法在这一方面显示出巨大的生命力30PFC培训第四部分第四部分 PFCPFC2D2D的使用的使用Ø1.1.对PFCPFC软件的使用界面、菜单功能及作用进行介绍；Ø2.FISH2.FISH语言简介Ø3.PFC2D3.PFC2D分析模型的生成方法Ø4.4.边界条件的设置方法Ø5.5.初始条件的设置Ø6.6.接触本构模型：①①接触刚度模型②②滑动模型③③连接模型Ø7.7.赋予材料属性：相关命令的使用方法介绍31PFC培训第四部分第四部分 PFCPFC2D2D的使用的使用Ø8.8.节理面的生成及属性设置Ø9.9.加载方法：主动荷载和被动荷载; ;Ø10.10.求解过程：静力求解、动力求解; ;Ø11.11.流体与热分析简介; ;Ø12.12.介绍PFCPFC2D2D软件的用户自定义本构模块的相关功能、操作等; ;Ø13.13.常用命令使用方法及相关的重要概念; ;Ø14.14.讲述PFCPFC2D2D工程应用的实例32PFC培训1.使用界面、菜单功能介绍33PFC培训34PFC培训2.FISH语言简介 FISHFISH是一种内置于是一种内置于ItascaItasca软软件内的件内的编编程程语语言，使用言，使用FISH FISH 用用户户可以定可以定义义新的新的变变量和函数，从而使得量和函数，从而使得这这些函数被用来些函数被用来扩扩展展ItascaItasca软软件的用法或者增加用件的用法或者增加用户户自定自定义义的特性。

例如，用的特性例如，用户户可可以以绘绘制制(PLOT)(PLOT)或打印或打印(PRINT)(PRINT)新的新的变变量，也能量，也能够够改改进进特殊的特殊的颗颗粒体模型粒体模型( (网格网格) )生成器，可以生成器，可以对对数数值试验进值试验进行伺服控制，可以行伺服控制，可以设设定一些性定一些性质质的特殊分布以及自的特殊分布以及自动动化参数研究化参数研究 ItascaItasca已已经经写了一些写了一些简单简单但非常有用的但非常有用的FISHFISH函数作函数作为库为库文文件包含在各个具体的件包含在各个具体的软软件中，件中，这这些些FISHFISH函数一方面方便了一些函数一方面方便了一些没有没有编编程程经验经验的用的用户户写一些写一些简单简单的的FISHFISH函数，另一方面也能使函数，另一方面也能使用用户户在在这这些提供的些提供的简单简单函数的基函数的基础础上作上作进进一步的改一步的改进进不过过，，FISHFISH语语言象其它言象其它编编程程语语言一言一样样，也可以，也可以编编制出非常复制出非常复杂杂的程序 35PFC培训 利用FISH语言进行编程，应该首先编一些简单的函数，然后仔细检查函数的功能，测试是否有错误。

如果没有发现错误，再逐渐增加其功能，增加一项功能检查一下，直至发展到最后比较复杂的程序这是因为虽然FISH是一种编译型语言，但它没有自己独立的编译器，不象VC++或VB能够实时全面地检查错误，FISH检查错误的能力很差，因此在使用他们到真实的应用之前，一定要用一些简单的数据（假如可能的话）来检查所有定义的函数 36PFC培训 FISH函数内置于标准的Itasca软件的数据文件中，函数的格式必须以DEFINE开始，以END结束函￿数可以嵌套调用，但定义函数的次序没有关系，只要在使用之前全部定义就行由于FISH函数的编译格式储存在Itasca软件的内存中，因此可以用SAVE命令保存函数以及相关变量的当前值￿￿￿￿FISH也可以用来改进用户写的本构模型，如例1：￿ DEF abc abc=22*3+5 END Print abc 37PFC培训对上例子稍作改进(例2)： new def abc hh=22 abc=hh*3+5 end1).稍有编程常识的人可以看出，执行上面的例子（PRINT abc），其结果与例1相同：abc=71.在这个函数中，我们首先把22赋值给变量hh，然后把这个变量带入abc的表达式中，因此二者的结果相同。

38PFC培训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已被赋值39PFC培训4). 下面的试验将进一步解释函数与变量之间的差别注意：Itasca软件的SET命令可以用来设置任何用户定义的FISH符号的值，与在FISH中使用的符号无关下面的例3建立在例2的基础之上，我们不使用NEW命令来清除内存中的值，因为我们想继续使用那些值: set abc=0 hh=0print hh print abcprint hh40PFC培训5).在这个例子中，我们首先把abc和hh都赋值为0，由于hh是一个变量，第一个Print命令显示当前hh的值，hh=0;第二个Print命令由于abc是一个函数名，因此执行abc函数，先前定义的abc=0不起作用，重新计算了hh和abc的值，因此第三个Print命令显示的值是它在abc函数内指定的值，即hh=22,例4是这个试验完整的命令。

41PFC培训 现在我们总结一下：Itasca软件的三个重要的命令PRINT,SET,HISTORY可以直接操作简单的FISH变量或函数如下图所示，其中var代表变量名或函数名￿￿￿￿￿HISTORY命令的用法在这不作重复介绍，即定义了FISH函数后，对其中的变量可以进行跟踪行跟踪，即HISTTORY Var 42PFC培训 象其它高级编程语言一样，FISH有执行循环命令的功能，标准的格式如下：￿￿￿￿￿￿￿￿￿￿￿LOOP var(expr1,expr2) END_LOOP 其中LOOP和END_LOOP是FISH语句，符号var代表循环变量，expr1和pxpr2代表表达式或者单个变量，下面的例子用循环命令计算从1到10的和以及乘积,见下例：43PFC培训 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结束44PFC培训关于LOOP的注意事项：1).FISH接受END_LOOP和ENDLOOP的写法，但不接受END LOOP这样中间有空格的写法，其它类似的命令有着同样的规则，如END_IF，END_COMMAND等命令；2).在上面的例子中，如果执行Print n或Print fish命令，你会看到n=11而不是10，注意：这不是FISH的错误，这是一个基本的计算机指令存储规则，当循环结束后，计数器的值保存的是n+1而不是n,所有的高级编程语言有着相同的规则45PFC培训1.DEFINE function END2.CASEOF expr Case n endcase3. IF expr1 test expr2 THEN ELSE ENDIF4. LOOP var (expr1, expr2) ENDLOOP46PFC培训5. LOOP WHILE expr1 test expr2ENDLOOP6. COMMAND ENDCOMMAND7. HISTORY var PRINT var SET var value PLOT add .sh fname47PFC培训 另外，在FISH中还有许多其它的预定义对象,其中一类是尺度变量(scalar variables),它们是单个的数字,下面是总的尺度变量: clock----时钟时间,单位是秒的100倍. unbal----最大不平衡力￿￿￿￿￿pi----圆周率￿￿￿￿￿step----目前的时步数目￿￿￿￿￿urand----0.0-1.0之间均匀分布的随机变量￿￿这仅是其中的一小部分,完全的列表以后再述。

48PFC培训 另一类非常有用的内置对象是固有函数(intrinsic functions),这些函数能在FISH内进行一些比较高级的数学运算,完整的列表见FISH 手册,下面给出其中的一部分: abs(a)----a的绝对值￿￿￿￿￿￿cos(a)----a的余玄(a为弧度) log(a)----a的底数为10的对数￿￿￿￿￿￿max(a,b)----返回a,b中的最大值￿￿￿￿￿￿sqrt(a)----a的平方根49PFC培训3.PFC2D计算模型的生成方法 有两个命令可用于生成颗粒流模型:BALL和GENER-ATE，其中，BALL命令是生成单个的颗粒，该命令生成的颗粒可与已存在的颗粒重叠，而GENERATE 可生成一系列指定数目的颗粒流，该命令生成的颗粒是不允许重叠的PFC2D里主要有两种类型的颗粒流：规则排列的和无规则排列的一系列规则排列的颗粒流可以用来模拟模拟结构部分，如梁，而不规则排列的颗粒流可用来模拟实体或内部结构无规则的颗粒材料，如岩石内部所包含的胶结颗粒50PFC培训 尽管颗粒的排列是随机的，但在颗粒模型生成后，整个模型的结构特性还是可能会受影响的，比如弱的结构面或各向异性。

对于无规律排列的颗粒流模型，一般不可能去描述它的初始接触力的量级大小，这必须在后期要经过一个压缩的过程才可能给予较好的评价 51PFC培训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 + r252PFC培训end_loopyc = yc + yincxc = x0 + radius * (row - (row/2) * 2)end_loopendset 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 show53PFC培训54PFC培训3.2不规则排列 无规则排列，即：对一个给定空隙率的区域，采用颗粒来充填其中需要进行填充的空隙，并确保整个模型保持平衡。

对于所能被填充的模型的初始空隙率，是有一个限制值，不能任意小对于某些空隙率的模型，颗粒的填充可以无接触地排列，对于其它情况的空隙率，颗粒又可以重叠排列 55PFC培训 目前没有一个普遍的方法来将紧密充填的圆形颗粒体限制在一个任意封闭的面积内这里来介绍两种方法如何充填特定半径的颗粒体达到我们所需要的空隙率￿￿￿￿￿第一种方法，首先建立封闭区域的边界(简称墙体)，然后在封闭区域内任意生成一系列无接触的颗粒，最后移动区域的限制墙体，至所需要的空隙率这种方法有三个缺点：1.区域的几何形状改变；2.收敛速度慢；3.最终的分布趋势是不均匀的56PFC培训 第二种方法：运用GENERATE命令生成颗粒体，同时配合关键词高斯分配，即指定颗粒体半径的上下限，然后相应分配一个标准差，同时配合FISH函数来选择颗粒半径，最终生成我们所需要的模型￿￿￿￿下面来介绍采用这种方法是如何进行操作的，重点介绍2种方法57PFC培训3.2.1 半径扩展法 首先需要介绍一个概念：半径放大系数m，它会引起模型空隙率的改变空隙率的定义： ——颗粒体的总面积；￿——封闭区域的总面积因此：58PFC培训 我们定义初始空隙率为￿￿￿￿，新生成的空隙率为￿￿￿￿,那么，则有： 如果整个模型使用相同的半径放大系数，则： 或59PFC培训实例：区域宽：10单位￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿区域高：5单位￿￿￿￿￿￿目标空隙率：0.12 颗粒体数目：300 最大最小半径比：1.5 初始假设一个m之值，可以求出初始的￿￿￿为：从而得到平均半径大小为： 60PFC培训 由初始给定的颗粒体最大直径和最小直径的比例r可得到颗粒体半径的上下限：￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿下限：￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿上限：￿￿￿￿￿￿￿￿￿￿￿￿￿￿通过上面的计算，同时编写FISH函数来生成满足需要的空隙率的模型，并能自动计算空隙率的大小，以便用户进行检验。

61PFC培训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) / mult^2 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程序程序2 2：：62PFC培训 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 63PFC培训 _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 64PFC培训 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.SAV65PFC培训 模型模型接触力接触力66PFC培训 第2种方法：指定颗粒体的半径，不限制颗粒的数目，使足够多的颗粒产生来达到所需要的空隙率。

但这种方法所￿带来的缺点是可能在局部区域造成大面积的颗粒重叠，这将会产生很大的挤压力，从而给予颗粒较大的初始速度，这就可能使得颗粒体脱离墙体的限制为避免此情况的发生，可通过初始的有限步循环计算将颗粒的动能减至零，然后再计算至平衡态￿￿￿￿3.2.2 挤压排斥法67PFC培训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 = 0程序程序3 3：：68PFC培训section loop n (1,n_max) r_ball = rlo + urand * (rhi - rlo) pvol_new = pvol_sum + pi * r_ball^2 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 ball x=x_ball y=y_ball rad=r_ball end_command pvol_sum = pvol_new count = count + 1 end_loopend_section tot_poros = 1.0 - pvol_sum / tot_vol command prop dens=1000.0 ks=s_stiff kn=n_stiff69PFC培训 end_command ii = out(' '+string(count)+' balls were created.') enddef quiet loop n (1,4) command cycle 5 ini xvel 0 yvel 0 end_command end_loopend explode quiet cyc 980 prop fric=0.2 cyc 250 plot wall ball print tot_poros save explode.SAV70PFC培训平衡前模型平衡前模型平衡后模型平衡后模型71PFC培训速度向量速度向量图颗粒粒间接触接触72PFC培训73PFC培训颗粒体粒体间的的连接关系接关系图74PFC培训3.3 不同尺度颗粒的分布 如果一个模型的不同区域需要不同尺寸的颗粒体分布时，我们就可以运用方法二(挤压排斥法)来￿生成目标模型。

具体操作：将区域分成两部分，左侧随机充填大直径颗粒体，右册随机充填小直径颗粒体，左侧和右侧的区分通过x和y坐标来控制，最终生成后，通过平衡计算使两侧颗粒体自由扩散直到目标模型，详见下例75PFC培训new set random set disk on ; 将球体看作单位厚度的圆盘状def setup n_stiff = 1e8 s_stiff = 1e8 width = 10.0 height = 5.0 poros = 0.12 rat = 1.5 mult_0 = 1.6 ; initial radius multiplication factor mult_a = 0 ; 左侧半径放大系数 mult_b = 0 ;右侧半径放大系数 id1 = 0 id2 = 0endsetup程序程序4 4：：76PFC培训 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); ---------------------------------------------def make_block tot_vol = (x2 - x1) * (y2 - y1) num = id2 - id1 + 1 n0 = 1.0 - (1.0 - poros) / mult_0^2 r0 = sqrt(tot_vol*(1.0 - n0)/(pi*num)) rlo = 2.0 * r0 / (1.0 + rat) rhi = rat * rlo command ; Generate reduced-radius particles gen id=id1,id2 rad=rlo,rhi x=x1,x2 y=y1,y2 prop dens=1000 ks=s_stiff kn=n_stiff range id=id1,id2 end_command get_poros mult = sqrt((1.0 - poros) / (1.0 - pmeas))end77PFC培训def get_poros sum = 0.0 bp = ball_head loop while bp # null if b_id(bp) >= 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 = pmeasend78PFC培训 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.SAV79PFC培训平衡前模型平衡前模型平衡后模型平衡后模型80PFC培训法向接触力法向接触力切向接触力切向接触力81PFC培训颗粒体粒体间的的连接关系接关系图82PFC培训3.4 运用模型生成“过滤器” 有些情况需要我们建立复杂区域形状的颗粒流模型，如右图，此时我们可运用模型生￿￿成过滤器来获得所需要的模型，￿￿即filter命令，其后由用户定义￿￿FISH函数来控制,其中,颗粒的半￿￿径通过fc_arg(0)进行检验,x和y 的坐标位置分别通过fc_arg(1) 和fc_arg(2)进行检验。

如果颗￿￿粒满足要求，则FISH函数值设为￿￿0，否则为1详见下例83PFC培训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 = _skipend程序程序5 5：：84PFC培训def 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) 85PFC培训 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 50086PFC培训87PFC培训88PFC培训4.边界条件 PFC2D中有三种边界条件，分别是：墙体边界、颗粒体边界和混合边界。

其中颗粒体边界又分为速度边界和受力边界￿￿1).墙体边界￿￿￿￿￿￿建模过程中，墙体可作为颗粒体的生成范围约束，但同时也可以将墙体作为边界来施加约束对于墙体，我们只能施加速度约束，而不能直接对其施加外力，因为运动定律对墙体是不适用的其速度由以下三个参数控制:线速度、角速度和旋转中心墙的运动是通过不断更新定义墙的基点的位置来描述采用WALL命令设置，如： wall id=1 x=1.0 y=1.0 spin=10.089PFC培训(2)颗粒体边界￿￿￿￿PFC2D中的模型可以将一连串的颗粒体作为边界条件基本方法：在模型紧密压缩至平衡后，我们通过FISH函数将与墙体相接触的颗粒体逐个提取，将这一系列的颗粒体采用共同的边界条件限制，最后删除初始的限制墙体，即实现了以颗粒体代替墙体来作为边界条件颗粒体边界条件分速度边界和外力边界，下例为速度边界程序实现90PFC培训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)程序程序6 6：：91PFC培训 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 20092PFC培训放大放大放大放大 颗粒体速度边界示意图93PFC培训 其前期操作与速度边界类似，差异处是先将边界颗粒体的速度初始化为零，再删除限制墙体，采用FISH函数反向施加颗粒体所受的不平衡力，运行一定的计算步至平衡态。

在此过程中，区域内部的颗粒体会因轻微的扰动而导致其自身产生运动而脱离边界颗粒体，此时前面将边界颗粒体的速度固定为零颗粒体外力边界94PFC培训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程序程序7 7：：95PFC培训(3)混合边界条件 混合边界即速度边界与外力边界同时存在，对于双轴压缩试验，其两侧端面不允许存在明显的几何变形，否则其边界力将会无效。

因此，此时需要用到混合边界，两侧端面采用外力边界来约束，确保其不发生几何变形，侧面Y向采用速度边界，使其速度线性增加，以此来模拟试块轴向变形96PFC培训5.初始条件 首先介绍一下应力的概念在连续介质中，应力是一个连续量，应力一般定义为作用于单位面积上的力，对于二维问题，应力则定义为作用于单位线段上的合力然而，对于松散介质这种定义是不合适的，因为颗粒介质是离散的，在颗粒模型中没有存在于每一个点上连续的应力，而且变化幅度很大对于应变也存在这样的问题在PFC模型中，计算颗粒间接触力和颗粒的位移，这些量对于在微观上研究材料的特性很有意义，但不能将这些量直接联系到连续模型中，需要经过一个平均的过程，才能从微观传到连续模型中所以,在PFC2D中采用平均应力的概念来表示97PFC培训5.1 施加等向应力def expand_radii_sum = 0.0 bp = ball_head loop while bp # null ;找到与bp相连的下一个颗粒的物理地址 cp = b_clist(bp) loop while cp # null if 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) else bp_other = c_ball1(cp) end_if 程序程序8 898PFC培训if pointer_type(bp_other) = 101 then ; 21, 100, 101, 102, 103 or 104 _phi = b_rad(bp) ; memory, ball,wall,contact,clump,circle else _phi = b_rad(bp) + b_rad(bp_other) end_if _sum = _sum + _rcp * c_kn(cp) * _phi end_if end_if if c_ball1(cp) = bp then cp = c_b1clist(cp) else cp = c_b2clist(cp) end_if end_loop bp = b_next(bp) end_loop99PFC培训 _alpha = -1.0 * _lambda * tot_vol * _diso / _sum bp = ball_head loop while bp # null b_rad(bp) = (1.0 + _alpha)*b_rad(bp) bp = b_next(bp) end_loopenddef achieve_isostrloop while 1 # 0 ;始终执行，除非遇到exit命令 _diso = req_isostr - isostr if abs(_diso/req_isostr) <= req_isostr_tol then exit end_if oo = out('Current isotropic stress = ' + string(isostr)) expand_radii command solve 100PFC培训 end_command end_loop enddef isostr ii = 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.0get_mphistory id=1 isostrSET req_isostr = -4e5 req_isostr_tol = 0.01cyc 20 print meas 1print isostrachieve_isostr101PFC培训6.接触本构模型 在PFC2D中，材料的本构特性是通过接触本构模型来模拟的。

每一颗粒的接触本构模型有:1)接触刚度模型;2)滑动模型;3)连接模型接触刚度模型提供了接触力和相对位移的弹性关系，滑动模型则强调切向和法向接触力使得接触颗粒可以发生相对移动，而连接模型是限制总的切向和法向力使得在连接强度范围内发生接触102PFC培训1).接触刚度模型￿￿￿￿￿接触刚度通过以下两式把接触力与相对位移联系起来，即: 式中，￿￿是法向刚度，它为割线刚度，联系总的法向力与位移￿￿￿￿￿￿￿￿￿￿式中，￿￿是切向刚度，它为切线刚度，通过增量形式联系切向力与位移￿￿￿￿￿在PFC2D中有两种接触刚度模型，即线性接触刚度模型与简化的Hertz-Minlin接触刚度模型,不同的接触刚度模型有不同的接触刚度值103PFC培训<1>.线性接触模型 分析颗粒间接触特性时，将这种接触的颗粒想象为一梁端点在颗粒中心弹性梁，梁端部受力或力矩相当于作用于颗粒的中心这种梁由以下特征参数描述:1)几何参数:长度(L)、断面积(A)、惯性矩(I);2)变形参数:杨氏模量E、泊松比v ; 3)强度参数:法向强度￿￿、切向强度￿￿;颗粒间接触杨氏模量用￿￿表示，平行连接杨氏模量￿￿￿表示。

同时假定PFC2D中所有颗粒均为厚度为t的圆盘若颗粒A和B接触，则梁的半径为: 104PFC培训梁的长度为: 式中: 和￿￿￿￿分别为接触颗粒的半径￿￿￿￿￿￿相互接触颗粒之间的受力特性，相当于弹性梁端部受纯轴向或纯切向荷载时的情况一致，此时梁断面各A与惯性矩I为: ，式中:t为假设的颗粒圆盘厚度 105PFC培训 关于接触变形，不仅指接触连接模型的接触连接变形，同样适用于颗粒之间以及颗粒与墙之间的接触变形，所以，下面介绍的计算方法即可用于接触连接材料也可用于非接触连接材料(砂土)对于纯轴向荷载或纯切向荷载，其法向接触刚度与切向接触刚度分别为：￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿；￿￿￿￿式中，￿￿是接触杨氏模量，不同于材料整体杨氏模量(通常大于整体杨氏模量)对于线性接触模型，其接触刚度k‘是假设相互接触两球为串联：￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿￿式中， 分别代表切向和 法向刚度。

式式1.106PFC培训 若两球有相同的刚度，即： ￿￿￿￿￿￿此时，通过公式转换，可得在接触连接时，接触模量与球的刚度间的关系式: 所以在描述变形微观参数时，先给定颗粒一颗粒接触的接触模量EC以及颗粒法向刚度与切向刚度比值kn/ks根据式1.计算kn再据给定比值求ks￿￿￿￿￿线性接触模型是通过两个接触体(颗粒—颗粒，颗粒一墙体)的切向刚度ks和法向刚度kn定义的计算线性接触模型的接触刚度时假定接触体为串联，根据式￿￿2.计算式式￿ ￿2.107PFC培训<2>.Hertz-Mindlin接触模型￿￿￿￿Hertz-Mindlin接触模型是Mindlin&Deresiewich (1953)和Cundall(1988)理论近似得出的一种非线性接触模型它只能严格用于颗粒体接触，而且不能再现剪切过程中的连续非线性(特别使用初始剪切模量时)￿￿￿￿￿￿￿Hertz-Mindlin模型由两个接触球体的参数剪切模量G和泊松比v来定义该模型不适用于采用接触连接的两球体的接触，因为这种模型没有定义球体受张力的情况对于颗粒—颗粒相互接触，其弹性常数为平均值，对于颗粒—墙体接触，假定墙为刚性，所以直接取球的弹性常数。

108PFC培训接触法向割线刚度和切向切线刚度为:式中Un为颗粒体接触重叠量，Fin为法向接触力109PFC培训 其它的系数为两接触体的几何与材料特性的函数对于颗粒一颗粒相互接触，这些系数可表示为: 对于颗粒—墙相互接触，这些系数可表示为: 式中，G为弹性剪切模量，v为泊松比，R为颗粒的半径，[A], [B]为相互接触的两球体110PFC培训2).滑动模型 滑动模型是相互接触球体的一种固有特性，它没有法向抗拉强度，允许颗粒在抗剪强度范围内发生滑动，这种模型在接触连接模型发生作用之前一直有效同时，滑动模型也可与半行连接模型同时起作用￿￿￿￿￿滑动模型是通过两接触体间最小摩擦系数u定义的，若颗粒间重叠量小于或等于零，则令法向和切向接触力等于零发生滑动的判别条件为: 若 ,则可以发生滑动，其中：￿￿￿￿￿￿￿￿￿￿￿￿￿￿111PFC培训3).连接模型 PFC2D模型允许相互接触颗粒连接在一起，有两种连接模型:即接触连接与平行连接模型接触连接认为连接只发生在接触点很小范围内，而平行连接发生在接触颗粒间圆形或方形有限范围内。

接触连接只能传递力，而平行连接同时能传递力和力矩两种类型的接触可以同时存在，直到超过接触强度连接模型只能是颗粒之间的连接，颗粒与墙之间的连接不能采用连接模型112PFC培训<1>.接触连接模型 接触连接可以想象为一对有恒定法向刚度与切向刚度的弹簧作用于颗粒接触点处，同时，这些弹簧设定一定的抗拉与抗剪强度只要接触连接存在就没有颗粒间的滑动发生接触连接在颗粒间重叠量小于零时，允许出现张力，但法向接触张力不能超过接触连接强度在PFC2D中，接触连接是有法向连接强度Fcn定义当法向抗拉接触力大于或等于法向接触连接强度时，连接破坏并把法向、切向接触力赋值为零当切向接触力大于或等于切向连接强度时，连接也发生破坏，但是接触力不发生变化，并假设切向力没有超过摩擦极限颗粒接触点处接触力与相对位移的关系的本构特性见下图113PFC培训 接触力的法向分量 接触力的切向分量114PFC培训<2>.平行连接模型 平行连接模型可以描述颗粒之间有限范围内有夹层材料的木构特性相互连接的两个颗粒可以看作是球体或柱体这种连接建立颗粒间一种弹性相互作用关系，可与前面所述的滑动模型或接触连接模型同时存在。

￿￿￿￿￿￿平行连接可以想象为一组有恒定法向刚度与切向刚度的弹簧均匀分布于接触平面内，这些弹簧作用的本构关系类似与点接触弹簧模拟颗粒刚度的本构特性接触的相对运动在平行连接处产生力和力矩，作用于相互连接的颗粒上，且与连接材料的最大法向、切向应力有关115PFC培训 平行连接模型是由法向刚度￿￿￿、切向刚度￿￿￿、法向强度￿￿￿、切向强度￿￿￿和连接半径￿￿￿五个参数定义的与平行连接相对应的总接触力和力矩用￿￿￿￿￿和￿￿￿表示，根据习惯，总的接触力和力矩表示平行连接在颗粒B上的作用，如下图所示可以将总的接触力沿接触面分解为切向分量和法向分量：￿￿￿￿￿￿￿￿￿￿￿式中，￿￿￿表示法向分量，￿￿表示切向分量116PFC培训 当连接形成时，￿￿￿和￿￿￿均初始化为零，以后在接触处由位移增量和旋转增量引起的弹性力和力矩的增量将迭加在当前值中117PFC培训7.赋予材料属性 PFC2D程序里，除了颗粒体联结相关的属性(联结刚度、强度以及平行联结半径)外，其它颗粒体的属性赋予均采用PROPERTY 命令来执行￿￿￿￿￿同时，在执行命令时，对与不同区域不同位置的颗粒体，其属性的变化可通过关键词￿￿gradient和group来配合使用，达到我们所需设置的要求和变化。

￿￿￿118PFC培训 如颗粒体间接触剪切强度￿￿随着X方向线形变化，其表达式为：￿￿￿￿￿￿￿￿￿￿￿我们可以这样来编写命令流文件：￿￿￿prop s_bond 1e6 grad -1e4, 0 range x 0 50具体求解为：￿￿x=15m处，s_bond=1e6+(-1e4)*15119PFC培训 同时还可将某个区域内的颗粒体均定义为一个组，用group来执行，然后再给这个组赋予属性，具体操作为： group strong_rock range x 0 10 y 5 10 group weak_rock range x 0 10 y 0 5 prop fric 1.0 range group strong_rock prop fric 0.1 range group weak_rock120PFC培训 如果颗粒体的材料属性是一个非线形变化的量，我们就需要采用FISH语言来赋予属性,如对于颗粒间摩擦系数表达式为：￿￿￿￿￿￿￿￿其中d为颗粒体距模型顶面距离，h为模型区域的宽度和高度,此时赋予属性时, FISH语言可这样来表述：￿￿￿￿121PFC培训def vary_friction bp = ball_head loop while bp # null y_pos = b_y(bp) depth = height - abs(y_pos) b_fric(bp) = fric0 * (1.0 + (depth/height)ˆ2) bp = b_next(bp) end_loopend SET height=5 fric0=0.2vary_friction122PFC培训8.节理面的生成及属性设置 PFC2D中，节理面可以用来模拟颗粒组间的滑移和分离，还可以模拟节理、断层和层理等，同时还能来模拟两种不同材料间的接触面，如基础和土体。

￿￿￿￿￿我们采用JSET命令来设置，可生成单独的一条或一组节理面但只有存在法向力的颗粒间的接触，才能被设置为节理面接触，因此，节理面命令必须在颗粒体模型生成后并达到平衡后才能执行123PFC培训 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.0124PFC培训125PFC培训JSET id 1 dip 50,2 number 5 spacing 50,3 上述命令是生成5条间距为单位50的节理面，第2至第5条节理面生成的过程中，倾角允许偏差为2度，间距值允许偏差为单位3￿￿赋予节理面属性时，同样采用PROPERTY命令，如：property n_bond= 0 s_bond= 0 range jset 1或：￿￿￿￿property fric 0.25 range jset 1,2126PFC培训9.加载 PFC2D荷载中的类型一般分为两类：主动荷载和被动荷载。

被动荷载如重力，采用SET GRAVITY 设置主动荷载则包括施加颗粒体速度和外力颗粒体可以直接施加速度和外力，但外力却不能直接施加在墙体上，因为运动方程对于墙体是无效的，但可以通过其他等效的方法来模拟墙体受集中力或应力127PFC培训 加载的操作主要从这四方面来执行：1).墙体的控制；2).颗粒体外力的控制；3).颗粒体速度控制；4).混合受力的控制￿￿￿￿￿上述几方面，其具体编程方法与前面所介绍的边界条件与初始条件类似，在这里就不详细介绍，具体的现实问题或工程所需的特殊的加载方式，需采用FISH来进行编写和操作，同时需从简单的模型来验证其实际控制效果128PFC培训10.求解 PFC2D中所谓的求解即为执行许多个计算循环步(使用CYCLE和STEP)，同时观察计算后的模型响应对于静力问题，模型系统要么计算达到平衡，要么破坏坍塌对于动力问题，就显得较复杂些，特别需要注意如何消除个别已脱离模型颗粒体虽然程序会对这些已脱离模型的颗粒体持续计算，但他们的变化对整个模型的影响会越来越小，因为只有当颗粒体处于计算的约束区域内才会产生其应有的效应一种简易的消除逃逸颗粒体的方法是每隔几千步计算就在一定范围内执行DELETE BALL 命令。

129PFC培训1).静力求解 对于求解命令CYCLE或STEP，用户也可以用SOLVE命令来进行静力求解该命令是给定某个不平衡力准则，当达到这个标准时，SOLVE命令即停止同时可配合一些关键词共同执行如Average，允许指定平均不平衡力与平均接触力的比值；maximum则允许指定最大不平衡力与最大接触力的比值，￿￿默认的值为0.01(1%) solve average=0.005 maximum=1e-10130PFC培训 SOLVE ratio命令可以设置模型计算不平衡力的比率大小，默认值为1×10-5，ratio的￿值可通过SET solve_ratio 命令来设置，一般该命令在热力学求解模式下使用较多￿￿￿￿￿Solve Clock t(以分钟计) Solve time t(以秒计) …… 对于静力计算的时间步，是与颗粒体的大小和接触刚度有关的，如果模型中颗粒体的大小和接触刚度变化幅度较大，时间步将不统一，为了能较快的收敛平衡，可以设置set dt dscale，该命令使颗粒体所有的自由度均采用相等的时间相应，这样有效的加快了模型的收敛。

131PFC培训2).动力求解 动力求解和静力求解的差异主要是动力求解存在一个现实的阻尼值但一般动力求解之前总是要进行静力求解，其目的是为了获得初始应力平衡分配的计算模型在此基础上来研究其进一步的动力上的力学响应 对于阻尼值的设定，我们需采用实际值，它的获取通过试验室测定或现场试验结果，程序默认的值为0.7采用PROPERTY命令执行132PFC培训def startup old_time = time omega = freq * 2.0 * pienddef shake while_stepping real_time = time - old_time if real_time < tlength vel = ampl * sin(omega*real_time) else vel = 0.0 end_if bp = ball_head loop while bp # null if b_xfix(bp) = 1 b_xvel(bp) = vel end_if133PFC培训 bp = b_next(bp) end_loopendproperty damp=0.157 set freq=2.0 ampl=15.0 tlength=1.0startup hist id=1 ball xvel 0 0 hist id=2 ball xvel 0 6 hist id=3 ball xvel 0 12 hist id=4 real_timeplot cur 0plot hist 1 2 3 v 4solve time 0.64134PFC培训135PFC培训136PFC培训137PFC培训 PFC2D中的热选项能够模拟材料热传导和存储，并且能够模拟由热效应导致的位移和力的发展。

热材料被热容器网（与每个粒子有关）和热管（与接触有关）所表示￿￿￿￿￿￿￿通过传导，在连接热容器的活跃的热管中出现热流动可辐射的和对流的热传递目前没有公式两个颗粒和连接材料的热效应由平行连接来说明，通过修改每个平行连接中的粒子的半径和力载荷可以产生热导致的应变 11.流体与热分析简介138PFC培训 热管与每个连接有关，单是这些管子中只有一些是活跃的如果接触的两个粒子是重叠的或者存在连接，那么这个管子就默认是活跃的（可以通过选项￿THERMAL pipe 来修改管子是否活跃）一旦热的宏观特性被设定，接下来的加载和破坏将改变活跃管的数目，并且改变材料的热传导能力例如，一个连接材料受压，新热管生成，热传导率将增加，但是当连接破坏，相连的管子变得不活跃，宏观热传导率又将减小￿139PFC培训常用命令：CONFIG thermalTHERMAL apply (deltemp, power, psource)THERMAL pipe active check off/onPROP densityTHERMAL prop sheat 指定热体积量[J/kg℃]THERMAL prop exp 指定线性膨胀率[1/℃]THERMAL make conductiv指定热传导率[W/m℃]THERMAL prop resista 指定单位热阻抗[℃/W m]140PFC培训 PFC2D中固定网格流选项被执行来模拟流体耦合。

这个设置用来解决连续性和笛卡儿坐标系中的N-S等式，考虑每个cell中存在的ball和他们的特性，然后得出每个固定网格的力和速度向量每个cell包含10或20个ball，流体的推动力被作为体力施加给ball相同的力被加到流体等式上，通过流向的压力剃度改变反映出来￿￿￿￿￿还有一些附加的特性可以用fish来实现例如压力，空隙率，每个流体细胞的速度向量用户还可以用C++语言或fish来自定义作用条件代替默认的模式 141PFC培训常用命令：Config fluidFluid model xl,xh,yl,yh size xs,ysFluid property (density, viscosity)Fluid boundary (slip,nonslip,pressure,velocity)set perm=1.0 crit=1e-5(1e-6) p_relax=pr 142PFC培训143PFC培训12.用户自定义本构 PFC2D程序允许用户用C++语言自己写代码在PFC2D中执行，自己写的函数可以在PFC命令行中运行，也可以在其他的地方被调用￿￿￿￿￿￿￿用户写的C++程序选项允许用户用C++语言写自己的程序，创建可执行的PFC2D个人版本。

这个选项可以用来代替FISH函数，大大提高运行的速度144PFC培训 FISH语言可用于编写新的本构模型如果编译成功的话，对用户而言，新模型就和PFC内置的模型一样，即可用命令来使用和停用模型自定义模型中也可拥有类似PFC内置参数一样的用户自己命名的特征参数￿￿￿￿￿用户自定义本构模型仅仅是包含特定语句，涉及与某个区域网格相关联的特定变量的FISH函数￿￿￿￿当新的自定义模型建好后，必须所有结点都固定的单一区域进行测试，给定的应变路径用于不同单元以分析其应力响应采用该方法，新模型的调试非常有效145PFC培训 用户自定义模型的运行速度比内置模型要慢模型优化后，正常的运行速度为内置模型的1/4到1/3但是，因为模型往往只发生在局部网格，所以用户自定义模型常常只需要在PFC网格的小范围数以用，其他大范围的材料都能用标准模型描述自定义模型虽然增加了运行时间，但是却节省了人力，克服了内置模型的不适用问题146PFC培训 同时，用户还可以运用C++语言来编写本构模型在C++语言中,重点是以面向对象的方法规划结构，利用类去表示对象同对象联系的书记由对象封装，且在对象外部是不可见的。

同对象的联系通过在封装数据中运行的成员函数实现另外，对于对象的分级制度有强大的支持—新对象类型可以由基类衍生，并且基类的成员函数可以由衍生类重载￿￿￿￿利用这些优点，我们编译好程度代码，并将其转化为DLL文件(动态数据库),它可以在任何需要的时候载入 147PFC培训 登陆Visual Studio .NET(Visual C++ 7.0)，选￿￿￿择 Solution，选择到udm.sln文件￿1. 编译工程￿￿￿￿￿￿￿￿选择Build/Rebuild solution，将在Release文件夹中产生一个名为userssoft.dll文件￿2. 创建一个新的工程￿￿￿￿￿￿￿对用户定义的模型要重新建立一个工程，将现有的工程(*.vcproj)拷贝到一个新的文件夹，具有有三个步骤：￿（1）拷贝udm.vcproj到mymodel.vcproj（mymodel为自己取得模型名字）￿（2）在udm solution图表上右键，选择Add/Existing project 148PFC培训（3）选择mymodel.vcproj 可以看到在solution explorer中有两个相同名字的工程（都成为udm）,右键选择第2个工程，将其名字改成mymodel。

3.选择Release/Debug编译选项（1）选择Build/Configuration Manager （2）选择Release或者Debug 注1：vcmodels.lib是一个Release编译文件（不包含任何调试信息），但是通过Debug选项也可以建立DLL文件并读入pfc2D，虽然它不能全部使用且可能会降低一些优化选项的运行（具体参照Debug和Release的区别）￿￿￿￿￿注2：下面的摄制过程仅对当前编译环境有效149PFC培训4. 改变定义模型的输出文件名（1）在mymodel的工程上右键，选择Properties / Configuration Properties / Linker / General （2）在Output File输入"Release/mymodel.dll" 5. 添加自定义模型的资源和头文件到工程中￿（1）右键usersoft.cpp，选择Remove，同样对￿￿￿￿￿￿￿￿usersoft.h进行移出操作（2）在上述的位置填入新的.cpp和.h文件（选择￿￿￿￿￿￿￿Add / Add Existing Item）￿￿￿（3）对输出dll（上述的第4步）进行更名，重新￿￿对第1步进行编译即可以得到所需的自定义模型文件.DLL。

150PFC培训 注意：DLL文件必须用Microsoft Visual C++ Version6.0(sp4)或更新版本编写才能在PFC2D3.0版本上运行和载入用户所需要的DLL模型 运行和装载： CONFIG cppudm Model load usercm.dll151PFC培训13.实用和常用的命令 Group Hist Range Movie如：Group name range x *,* y #,#Group name2 range group name1Group name2 range group name1 not Group name3 range group name1 any group & name2 any152PFC培训Hist ball xvel(yvel xpos ypos) *,* (id *)Hist ball s11(s12 s22) *,* (id *)Hist ball energy body(bond bound fric kinetic strain)Hist write 1 skip 20 Hist write 1 v 3 skip 20 Hist write 1 v 3 begin * end * 153PFC培训range name master_region x=(20,90) & y=(20,90) circle center=(20,20) rad=80.0Group master range master_region等效于：Group master range x=(20,90) y=(20,90) & circle center=(20,20) rad=80.0 若对之前已定义的区域再次定义name，以最新定义的name为准。

154PFC培训Avi格式的影音文件： set plot avi size 640 480 movie avi open movie step 10 1(所指窗口id) cycle 1000 movie avi close step2d.aviDCX格式的影音文件： set plot dcx size 640 480 movie step 10 1 step2d.dcx演示155PFC培训实例一. 矿区垮落过程模拟 各种不同组成的矿块在其破碎垮落过程中相互影响，其垮落过程以及相互作用情况均可使用PFC2D合理模拟这些不同的组成包括：矿块的整体跨度、岩体强度、变形能力、节理结构特性(如方向、间距、强度)以及初始应力等这些方面性质的改变均可造成结果的改变￿￿￿￿￿本例呈现了一理想化的垮落计算模型，模型为一包含两组相交节理的矿体14.工程应用实例工程应用实例156PFC培训节理组 倾角 倾向 间距 摩擦系数 1 30° 90° 40m 0.0 2 0° 90° 25m 0.0节理面性理面性质157PFC培训模拟过程 岩块是由统一半径的2000个颗粒组成的，整个矿区的模型大小为：水平×垂直＝300m × 200m。

通过球体的扩展计算至压密状态达到所需要的应力状态(约70MPa)平衡状态时应力组成为SXX、SYY和SXY为-67.2MPa、－66.8MPa和￿￿0.61MPa在平衡后，颗粒体聚集收敛后执行JSET命令设置节理面由矿区底部某部分矿块的松动、塌落直至最后影响周边相临块体的运动，本模拟将整个塌落过程完全模拟，同时监控各块体的变化情况158PFC培训初始生成颗粒体平衡计算模型159PFC培训160PFC培训161PFC培训162PFC培训163PFC培训164PFC培训165PFC培训166PFC培训 实例二. 梁柱结构的动力分析问题描述：￿￿￿￿￿￿此次计算是采用PFC2D分析一两跨五层的梁柱结构在地表存在一周期变化的正弦波的情况下，其结构变形、内力及某些关键点的速度变化情况主要以此例来演示PFC2D的动力分析方法结构初始在重力下的自稳以及后期的动力影响甚至破坏，均可以较好的模拟分析在重力状态下的初始状态，结构的反应是弹性的，随着动力荷载及初始静力荷载，颗粒体间的连接逐渐破坏，直至最后的完全垮落167PFC培训初始初始计算模型算模型168PFC培训初始颗粒接触力平行连接关系169PFC培训平行连接力初始位移向量170PFC培训1.4秒后的变形及平行连接力171PFC培训1.4秒后的速度及位移172PFC培训173PFC培训174PFC培训 0.33秒 0.53秒175PFC培训 1.1秒 1.4秒176PFC培训 2.2秒 2.9秒177PFC培训 4.3秒 38秒178PFC培训179PFC培训动 画 演 示播 放180PFC培训实例三.颗粒体由漏斗内倾卸的模拟 如右图，漏斗内紧密充填颗粒体，上部伺服板处于活动状态，可均匀向下运动。

漏斗两侧的墙体固定，当底板拆除后，可限制颗粒体的运动同时，上部伺服板向下运动，自动调整运动速度以确保颗粒体与两侧墙体的接触力保持不变181PFC培训模拟过程： 第一阶段：建立所指定空隙率和颗粒半径的模型， 计算至平衡压密阶段 第二阶段：伺服板下降，压密颗粒直至达到所需的 颗粒间的接触压力 第三阶段：底板拆除，颗粒体流动，同时伺服板向 下运动以保证颗粒体与漏斗两侧墙体的 压力保持衡值 在这样的过程中，伺服板的压力和速度的变化情况，以及颗粒接触间的摩擦能量都会被记录，以此来更好的了解整个系统182PFC培训第一第一阶段段183PFC培训第二第二阶段段伺服板速度变化曲线184PFC培训伺服板应力变化曲线第二第二阶段段185PFC培训摩擦耗能变化曲线图186PFC培训187PFC培训伺服板应力变化曲线图188PFC培训摩擦耗能变化曲线图189PFC培训190PFC培训 将球体-墙体和球体-球体间的摩擦系数设为2.6，计算将得到另一种结果。

191PFC培训192PFC培训193PFC培训194PFC培训伺服板压力变化曲线图195PFC培训摩擦耗能变化曲线图196PFC培训THE END 谢谢大家！谢谢大家！197PFC培训。