PFC中文详解教育课件

PFC中文详解PPT讲座PFC中文详解PPT讲座1主要内容第一部分PFC2D颗粒流程序简介第二部分有限差分法基础介绍第三部分离散元法基础介绍第四部分PFC2D的使用主要内容第一部分PFC2D颗粒流程序简介2第一部分PFC2D颗粒流程序简介1、理论背景2、颗粒流方法的基本假设3、颗粒流方法的特点4、可选特性5、应用领域6、求解步骤第一部分PFC2D颗粒流程序简介1、理论背景3

作为离散元的一种，二维颗粒流程序(ParticleFollowCodePFC2D)数值模拟新技术，其理论基础是Cundall[1979]提出的离散单元法，用于颗粒材料力学性态分析，如颗粒团粒体的稳定、变形及本构关系，专门用于模拟固体力学大变形问题。它通过圆形(或异型)离散单元来模拟颗粒介质的运动及其相互作用。由平面内的平动和转动运动方程来确定每一时刻颗粒的位置和速度。作为研究颗粒介质特性的一种工具，它采用有代表性的数百个至上万个颗粒单元，通过数值模拟实验可以得到颗粒介质本构模型。1.理论背景作为离散元的一种，二维颗粒流程序(Particle4

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

PFC2D(ParticleFollowC52、颗粒流方法的基本假设

颗粒流方法在模拟过程中作了如下假设:1)颗粒单元为刚性体;2)接触发生在很小的范围内，即点接触;3)接触特性为柔性接触，接触处允许有一定的“重叠”量;4)“重叠”量的大小与接触力有关，与颗粒大小相比，“重叠”量很小;5)接触处有特殊的连接强度;6)颗粒单元为圆盘形(或球形)。2、颗粒流方法的基本假设颗粒流方法在模拟过程中作了如下假设6

其中，颗粒为刚性体的假设，对于模拟介质运动为只沿相互接触面的表面发生的问题非常重要，比如象砂土或粮食这种颗粒组合体材料，利用这种假设在总体上来讲是比较恰当的，因为这种材料的变形是来自于颗粒刚性体间的滑动和转动以及接触面处的张开和闭锁，而不是来自于每个刚性颗粒本身的变形，对于这种特殊材料，没有必要采用非常精确的数值模型，来得到对材料特性的近似。其中，颗粒为刚性体的假设，对于模拟介质运动为只73、颗粒流方法的特点

PFC2D可以直接模拟圆形颗粒的运动和相互作用问题。颗料可以代表材料中的个别颗粒，例如砂粒，也可以代表粘结在一起的固体材料，例如混凝土或岩石。当粘结以渐进的方式破坏时，它能够破裂。粘结在一起的集合体可以是各向同性，也可以被分成一些离散的区域或块体。这类物理系统可以用处理角状块体的离散单元程序UDEC和3DEC来模拟。3、颗粒流方法的特点PFC2D可以直接模拟8

PFC2D有三个优点：

第一、它有潜在的高效率。因为圆形物体间的接触探测比角状物体间的更简单。

第二、对可以模拟的位移大小实质上没有限制。

第三、由于它们是由粘结的粒子组成，块体可以破裂，不象UDEC和3DEC模拟的块体不能破裂。用PFC2D模拟块体化系统的缺点是，块体的边界不是平的，用户必须接受不平的边界以换取PFC2D提供的优点。PFC2D有三个优点：9PFC2D中几何特征、物理特性和解题条件的说明不如FLAC和UDEC程序那样直截了当。例如用连续介质程序，创建网格、设置初始压力、设置固定或自由边界。在象PFC2D这样的颗粒程序中，由于没有唯一的方法在一个指定的空间内组合大量的粒子，粒子紧密结合的状态一般不能预先指定。必须跟踪类似于物体压实的过程，直到获得要求的孔隙率。由于颗粒相对位置变化产生接触力，初始应力状态的确定与初始压密有关。由于边界不是由平面组成，边界条件的设定比连续介质程序更复杂。PFC2D中几何特征、物理特性和解题条件的说明不如F10

当要求满足有实验室实际测试的模拟物体的力学特性时，出现了更大的困难。在某种程度上，这是一个反复试验的过程，因为目前还没有完善的理论可以根据微观特性来预见宏观特性。然而，给出一些准则应该有助于模型与原型的匹配，如哪些因素对力学行为的某些方面产生影响，哪些将不产生影响。应该意识到，由于受现有知识的限制，这样的模拟很难。然而，用PFC2D进行试验，对固体力学，特别是对断裂力学和损伤力学，可以获得一些基本认识。

当要求满足有实验室实际测试的模拟物体的力学特性时11PFC2D能模拟任意大小圆形粒子集合体的动态力学行为。粒子生成器根据粒子的指定分布规律自动概率地生成。粒子半径按均匀分布或按高斯分布规律分布。初始孔隙度一般比较高，但通过控制粒子半径的扩大可以获得密度压实。在任何阶段任何因素都可以改变半径。所以不需反复试验就可以获得指定孔隙度的压实状态。

PFC2D能模拟任意大小圆形粒子集合体的动态力学行为。12

属性与各个粒子或接触有关，而不是与“类型号”有关。因此，可以指定属性和半径的连续变化梯度。“节理生成器”用来修改沿指定轨迹线的接触特性。假定这些线叠加在颗粒集合体上。用这种方法，模型可以被成组的弱面，如岩石节理切割。粒子颜色也是一种属性，用户可以指定各种标记方案。属性与各个粒子或接触有关，而不是与13

PFC2D模型中为了保证数据长期不漂移，用双精度数据存储坐标和半径。接触的相对位移直接根据坐标而不是位移增量计算。接触性质由下列单元组成：

1）线性弹簧或简化的Hertz-Mindlin准则；

2）库仑滑块；

3）粘结类型：粘结接触可承受拉力，粘结存在有限的抗拉和抗剪强度。可设定两种类型的粘结，接触粘结和平行粘结。这两种类型粘结对应两种可能的物理接触：①接触粘结再现了作用在接触点一个很小区域上的附着作用；②平行粘结再现了粒子接触后浇注其它材料的作用（如水泥灌浆）。平行粘结中附加材料的有效刚度具有接触点的刚度。PFC2D模型中为了保证数据长期不漂移，用双14

块体逻辑支持附属粒子组或块体的创建，促进了程序的推广普及。块体内粒子可以任意程度的重叠，作为刚性体具有可变形边界的每一个块体，可作为一般形状的超级粒子。通过指定墙的速度、混合的粒子速度、施加外力和重力来给系统加载。“扩展的FISH库”提供了在集合体内设置指定应力场或施加应力边界条件的函数。时步计算是自动的，包括因为Hertz接触模型刚度变化的影响。模拟过程中，根据每个粒子周围接触数目和瞬间刚度值，时步也在变化。基于估计的粒子数，单元映射策略采用最佳的单元数目，自动调整单元的外部尺寸来适应粒子缺失和指定的新对象。单元映射方案支持接触探测算法以保证求解时间随粒子数目线性增加，而不是二次方增加。

块体逻辑支持附属粒子组或块体的创建，促进了程15

类似于FLAC，PFC提供了局部无粘性阻尼。这种阻尼形式有以下优点：

1）对于匀速运动，体力接近于零，只有加速运动时才有阻尼；

2）阻尼系数是无因次的；

3）因阻尼系数不随频率变化，集合体中具有不同自然周期的区域被同等阻尼，采用同样的阻尼系数。

PFC2D可以在半静态模式下运行以保证迅速收敛到静态解，或者在完全动态模式下运行。

PFC2D包含功能强大的内嵌式程序语言FISH，允许用户定义新的变量和函数使数值模型适合用户的特殊需求。例如，用户可以定义特殊材料的模型和性质、加载方式、实验条件的伺服控制、模拟的顺序以及绘图和打印用户定义的变量等。

类似于FLAC，PFC提供了局部无粘性阻尼。这种阻164、可选特性

1）热学分析2）并行处理技术3）能写用户定义接触模型4）用户写C++程序的C++编程。

4、可选特性1）热学分析17

热学选项用来模拟材料内热量的瞬间流动和热诱导位移和力的顺序发展。热学模型可以独立运行或耦合到力学模型。通过修改粒子半径和平行粘结承受的力，产生热应变来解释粒子和粘结材料的受热。

用户定义的接触本构模型可以用C++语言来编写，并编译成动态链接库文件，一旦需要就可以加载。

用户写的C++程序选项允许用户用C++语言写自己的程序，创建可执行的PFC2D个人版本。这个选项可以用来代替FISH函数，大大提高运行的速度。

并行处理技术允许将一个PFC2D模型分成几个部分，每个部分可以在单独的处理器上平行运行。与一个PFC2D模型在一个处理器上运行相比，平行处理在内存容量和计算速度方面得到大大提高。热学选项用来模拟材料内热量的瞬间流动和热诱导位移和185、应用领域

PFC2D既可解决静态问题也可解决动态问题，既可用于参数预测，也可用于在原始资料详细情况下的实际模拟。PFC2D模拟试验可以代替室内试验。在岩石与土体中开挖问题的研究与设计方面，实测资料相对较少，关于初始应力、不连续性等问题也只能部分了解。而在松散介质流动问题中，影响流动介质不规律分布的影响因素很难定量描述。因此，应用PFC2D初步研究影响整个系统的一些参数的特性，对整个系统的特性有所了解后，就可以方便地设计模型模拟整个过程。

5、应用领域PFC2D既可解决静态问题也可解决19

PFC2D可以模拟颗粒间的相互作用问题、大变形问题、断裂问题等，适用于以下领域：

（1）在槽、管、料斗、筒仓中松散物体的流动问题；

（2）矿山冒落法开采中的岩体断裂、坍塌、破碎和岩块的流动问题；

（3）铸模中粉料的压实问题；

（4）由粘结粒子组成物体的碰撞及其动态破坏；

（5）梁结构的地震响应及垮塌；

（6）颗粒材料的基本特性研究，如屈服、流动、体积变化等；

（7）固体的基本特性研究，如累积破坏、断裂。

PFC2D可以模拟颗粒间的相互作用问题、大变206、求解步骤1)定义模拟对象根据模拟意图定义模型的详细程序，假如只对某一力学机制的不同解释作出判断时，可以建立一个比较粗略的模型，只要在模型中能体现要解释的机制即可，对所模拟问题影响不大的特性可以忽略。6、求解步骤1)定义模拟对象212)建立力学模型的基本概念首先对分析对象在一定初始特性形成初步概念。为此，应先提出一些问题，如系统是否将变为不稳定系统、问题变形的大小、主要力学特性是否非线性、是否需要定义介质的不连续性、系统边界是实际边界还是无限边界、系统结构有无对称性等。综合以上内容来描述模型的大致特征，包括颗粒单元的设计、接触类型的选择、边界条件的确定以及初始平衡状态的分析。2)建立力学模型的基本概念223)构造并运行简化模型在建立实际工程模型之前，先构造并运行一系列简化的测试模型，可以提高解题效率。通过这种前期简化模型的运行，可对力学系统的概念有更深入的了解，有时在分析简化模型的结果后(例如所选的接触类型是否有代表性、边界条件对模型结果的影响程度等)，还需将第二步加以修改。3)构造并运行简化模型234)补充模拟问题的数据资料模拟实际工程问题需要大量简化模型运行的结果，对于地质力学来说包括:a)几何特性，如地下开挖酮室的形状、地形地貌、坝体形状、岩土结构等；

b)地质构造位置，如断层、节理、层面等；

c)材料特性，如弹/塑性、后破坏特性等；

d)初始条件，如原位应力状态、孔隙压力、饱和度等；

e)外荷载，如冲击荷载、开挖应力等。因为一些实际工程性质的不确定性(特别是应力状态、变形和强度特性)，所以必须选择合理的参数研究范围。第三步简化模型的运行有助于这项选择，从而为更进一步的试验提供资料。4)补充模拟问题的数据资料245)模拟运行的进一步准备

a)合理确定每一时步所需时间，若运行时间过长，很难得到有意义的结论，所以应该考虑在多台计算机上同时运行。

b)模型的运行状态应及时保存，以便在后续运行中调用其结果。例如如果分析中有多次加卸荷过程，要能方便地退回到每一过程，并改变参数后可以继续运行。

c)在程序中应设有足够的监控点(如参数变化处、不平衡等)，对中间模拟结果随时作出比较分析，并分析颗粒流动状态。5)模拟运行的进一步准备256)运行计算模型在模型正式运行之前先运行一些检验模型，然后暂停，根据一些特性参数的试验或理论计算结果来检查模拟结果是否合理，当确定模型运行正确无误时，连接所有的数据文件进行计算。

7)解释结果计算结果与实测结果进行分析比较。图形应集中反应要分析的区域如应力集中区，各种计算结果应能方便地输出分析。6)运行计算模型26第二部分有限差分法基础介绍

连续介质三维快速拉格朗日有限差分计算方法(FLAC3D)是近20年来逐步成熟完善起来的一种新型数值计算方法,它基于显式差分法来求解运动方程和动力方程,可模拟岩土或其他材料的三维力学行为。其求解时首先将计算区域离散化,分成若干三维单元,单元之间由节点联结,节点受荷载作用后,其平衡方程(运动方程)可以写成时间步长为Δt的有限差分形式,由于采用动态应力松弛显式差分求解技术,在某一微小的时段内,作用于该节点的荷载只对周围若干节点有影响。第二部分有限差分法基础介绍连续介质三维快速27

根据单元节点的速度变化和时段Δt,可求出单元之间的相对位移,进而求出单元应变,利用单元材料的本构关系即可求出单元应力。在此基础上,求出单元之间的不平衡力,将此不平衡力重新作用到节点上,再进行下一步的迭代过程,直到整个系统不平衡力足够小或节点位移趋于平衡为止。

FLAC3D可以解决诸多的有限元程序难以模拟的复杂的工程问题，例如分布开挖、大变形、非线性及非稳定系统（甚至大面积屈服/失稳或完全塌方）。根据单元节点的速度变化和时段Δt,可求出单元之间28第三部分离散元法基础介绍

离散单元法是一种模拟非连续介质的计算方法，自Cundall在70年代提出以来，在岩石力学、土力学、结构分析等领域的数值模拟中得到广泛应用，是一种新兴的非连续体分析方法。离散单元法允许单元间的相对运动，不一定满足位移连续和变形协调条件，计算速度快，所需存储空间小，特别适用于节理岩体的大位移，大变形分析。离散单元法自问世以来有了长足的发展，已经成为解决岩石力学问题的一种重要的数值方法，因为工程中所见到的岩体其形态呈非连续结构，所形成的岩石块体运动和受力情况多是几乎或材料非线性问题，所以很难用解决连续介质力学问题的有限单元法或边界单元法等。第三部分离散元法基础介绍离散单元法是一种模拟非29

数值方法来进行求解，而离散单元法正是充分考虑到岩体结构的不连续性，适用于解决节理岩石力学问题。近年来，离散元法的应用领域又扩展到求解连续介质向非连续介质转化的力学问题。混凝土等脆性材料在冲击、侵彻等动荷载作用下产生的损伤和破坏，其实质是力学模型从连续体到非连续体的转变过程。建立在传统的连续介质力学基础上的有限元法等数值计算方法难以直接用于计算和模拟材料具体的破坏形式和破坏的整个过程，而离散元法在这一方面显示出巨大的生命力。数值方法来进行求解，而离散单元法正是充分考虑到岩30第四部分PFC2D的使用1.对PFC软件的使用界面、菜单功能及作用进行介绍；2.FISH语言简介3.PFC2D分析模型的生成方法4.边界条件的设置方法5.初始条件的设置6.接触本构模型：①接触刚度模型②滑动模型③连接模型7.赋予材料属性：相关命令的使用方法介绍第四部分PFC2D的使用1.对PFC软件的使用界面、菜单31第四部分PFC2D的使用8.节理面的生成及属性设置9.加载方法：主动荷载和被动荷载;10.求解过程：静力求解、动力求解;11.流体与热分析简介;12.介绍PFC2D软件的用户自定义本构模块的相关功能、操作等;13.常用命令使用方法及相关的重要概念;14.讲述PFC2D工程应用的实例第四部分PFC2D的使用8.节理面的生成及属性设置321.使用界面、菜单功能介绍1.使用界面、菜单功能介绍33PFC中文详解教育课件342.FISH语言简介

FISH是一种内置于Itasca软件内的编程语言，使用FISH用户可以定义新的变量和函数，从而使得这些函数被用来扩展Itasca软件的用法或者增加用户自定义的特性。例如，用户可以绘制(PLOT)或打印(PRINT)新的变量，也能够改进特殊的颗粒体模型(网格)生成器，可以对数值试验进行伺服控制，可以设定一些性质的特殊分布以及自动化参数研究。

Itasca已经写了一些简单但非常有用的FISH函数作为库文件包含在各个具体的软件中，这些FISH函数一方面方便了一些没有编程经验的用户写一些简单的FISH函数，另一方面也能使用户在这些提供的简单函数的基础上作进一步的改进。不过，FISH语言象其它编程语言一样，也可以编制出非常复杂的程序。

2.FISH语言简介FISH是一种内置于Ita35

利用FISH语言进行编程，应该首先编一些简单的函数，然后仔细检查函数的功能，测试是否有错误。如果没有发现错误，再逐渐增加其功能，增加一项功能检查一下，直至发展到最后比较复杂的程序。这是因为虽然FISH是一种编译型语言，但它没有自己独立的编译器，不象VC++或VB能够实时全面地检查错误，FISH检查错误的能力很差，因此在使用他们到真实的应用之前，一定要用一些简单的数据（假如可能的话）来检查所有定义的函数。

利用FISH语言进行编程，应该首先编一些简单36

FISH函数内置于标准的Itasca软件的数据文件中，函数的格式必须以DEFINE开始，以END结束。函数可以嵌套调用，但定义函数的次序没有关系，只要在使用之前全部定义就行。由于FISH函数的编译格式储存在Itasca软件的内存中，因此可以用SAVE命令保存函数以及相关变量的当前值。

FISH也可以用来改进用户写的本构模型，如例1：

DEFabcabc=22*3+5ENDPrintabcFISH函数内置于标准的Itasca软件的数据文件中37对上例子稍作改进(例2)：

newdefabchh=22abc=hh*3+5end1).稍有编程常识的人可以看出，执行上面的例子（PRINTabc），其结果与例1相同：abc=71.在这个函数中，我们首先把22赋值给变量hh，然后把这个变量带入abc的表达式中，因此二者的结果相同。对上例子稍作改进(例2)：382).FISH的执行过程如下：当在程序命令中使用一个FISH符号名时（例如执行PRINT符号名），如果符号名也是一个函数名，那么执行这个函数（例如abc）；如果符号名不是函数名，那么使用符号目前的值(例如hh)。3).在输入完例2的各行后，如果我们执行命令：PRINThh，此时hh=0,因为在这个时候没有执行FISH函数，因此hh的初始值为0；我们接着执行PRINTabc，结果显示abc=71;再次执行PRINThh,此时结果为hh=22，这是因为我们首先运行了abc函数，在这个过程中hh已被赋值。2).FISH的执行过程如下：当在程序命令中使用一个FISH394).下面的试验将进一步解释函数与变量之间的差别。注意：Itasca软件的SET命令可以用来设置任何用户定义的FISH符号的值，与在FISH中使用的符号无关。下面的例3建立在例2的基础之上，我们不使用NEW命令来清除内存中的值，因为我们想继续使用那些值:setabc=0hh=0printhhprintabcprinthh4).下面的试验将进一步解释函数与变量之间的差别。注意：I405).在这个例子中，我们首先把abc和hh都赋值为0，由于hh是一个变量，第一个Print命令显示当前hh的值，hh=0;第二个Print命令由于abc是一个函数名，因此执行abc函数，先前定义的abc=0不起作用，重新计算了hh和abc的值，因此第三个Print命令显示的值是它在abc函数内指定的值，即hh=22,例4是这个试验完整的命令。5).在这个例子中，我们首先把abc和hh都赋值为0，由于h41

现在我们总结一下：Itasca软件的三个重要的命令PRINT,SET,HISTORY可以直接操作简单的FISH变量或函数。如下图所示，其中var代表变量名或函数名。

HISTORY命令的用法在这不作重复介绍，即定义了FISH函数后，对其中的变量可以进行跟踪，即HISTTORYVar。

现在我们总结一下：Itasca软件的三个重要的命令42

象其它高级编程语言一样，FISH有执行循环命令的功能，标准的格式如下：

LOOPvar(expr1,expr2)

END_LOOP

其中LOOP和END_LOOP是FISH语句，符号var代表循环变量，expr1和pxpr2代表表达式或者单个变量，下面的例子用循环命令计算从1到10的和以及乘积,见下例：象其它高级编程语言一样，FISH有执行循环命令的功能43newdefabcsum=0prod=1loopn(1,10)sum=sum+nprod=prod*nend_loopendabcprintsum,prod

在这个例子中，首先给两个变量赋于初始值，sum用来保存和的结果，prod用来保存积的结果，然后执行循环，最后分别打印出这两个变量的最后结果。循环变量n(1,10)表示从1开始，连续计算到10结束。

44关于LOOP的注意事项：1).FISH接受END_LOOP和ENDLOOP的写法，但不接受ENDLOOP这样中间有空格的写法，其它类似的命令有着同样的规则，如END_IF，END_COMMAND等命令；2).在上面的例子中，如果执行Printn或Printfish命令，你会看到n=11而不是10，注意：这不是FISH的错误，这是一个基本的计算机指令存储规则，当循环结束后，计数器的值保存的是n+1而不是n,所有的高级编程语言有着相同的规则。关于LOOP的注意事项：451.DEFINEfunctionEND2.CASEOFexprCasenendcase3.IFexpr1testexpr2THENELSEENDIF4.LOOPvar(expr1,expr2)ENDLOOP1.DEFINEfunction465.LOOPWHILEexpr1testexpr2ENDLOOP6.COMMANDENDCOMMAND7.HISTORYvarPRINTvarSETvarvaluePLOTadd.shfname5.LOOPWHILEexpr1testexpr247

另外，在FISH中还有许多其它的预定义对象,其中一类是尺度变量(scalarvariables),它们是单个的数字,下面是总的尺度变量:clock----时钟时间,单位是秒的100倍.unbal----最大不平衡力

pi----圆周率

step----目前的时步数目

urand----0.0-1.0之间均匀分布的随机变量这仅是其中的一小部分,完全的列表以后再述。另外，在FISH中还有许多其它的预定义对象,其中48

另一类非常有用的内置对象是固有函数(intrinsicfunctions),这些函数能在FISH内进行一些比较高级的数学运算,完整的列表见FISH手册,下面给出其中的一部分:abs(a)----a的绝对值

cos(a)----a的余玄(a为弧度)log(a)----a的底数为10的对数

max(a,b)----返回a,b中的最大值

sqrt(a)----a的平方根另一类非常有用的内置对象是固有函数(intrin493.PFC2D计算模型的生成方法

有两个命令可用于生成颗粒流模型:BALL和GENER-ATE，其中，BALL命令是生成单个的颗粒，该命令生成的颗粒可与已存在的颗粒重叠，而GENERATE可生成一系列指定数目的颗粒流，该命令生成的颗粒是不允许重叠的。PFC2D里主要有两种类型的颗粒流：规则排列的和无规则排列的。一系列规则排列的颗粒流可以用来模拟模拟结构部分，如梁，而不规则排列的颗粒流可用来模拟实体或内部结构无规则的颗粒材料，如岩石内部所包含的胶结颗粒。3.PFC2D计算模型的生成方法有两个命令可用于50

尽管颗粒的排列是随机的，但在颗粒模型生成后，整个模型的结构特性还是可能会受影响的，比如弱的结构面或各向异性。对于无规律排列的颗粒流模型，一般不可能去描述它的初始接触力的量级大小，这必须在后期要经过一个压缩的过程才可能给予较好的评价。

尽管颗粒的排列是随机的，但在颗粒模型生成后，整个513.1规律排列颗粒流

生成规则排列的颗粒流，主要采用FISH语言配合BALL命令，循环生成一系列的颗粒，如下例:Newdefhexxc=x0yc=y0rc=radiusidc=id_startr2=2.0*radiusyinc=radius*sqrt(3.0)looprow(1,n_row)loopcol(1,n_col)commandballid=idcx=xcy=ycrad=rcend_commandidc=idc+1xc=xc+r23.1规律排列颗粒流生成规则排列的颗粒流，主要52end_loopyc=yc+yincxc=x0+radius*(row-(row/2)*2)end_loopendsetechooffsetx0=0.2y0=0.4radius=0.1setid_start=100n_col=7n_row=8hexsetechoonplotsetcapsize20plotaddaxesblackplotaddballyellowplotshowPFC中文详解教育课件53PFC中文详解教育课件543.2不规则排列

无规则排列，即：对一个给定空隙率的区域，采用颗粒来充填其中需要进行填充的空隙，并确保整个模型保持平衡。对于所能被填充的模型的初始空隙率，是有一个限制值，不能任意小。对于某些空隙率的模型，颗粒的填充可以无接触地排列，对于其它情况的空隙率，颗粒又可以重叠排列。

3.2不规则排列无规则排列，即：对一个给定55

目前没有一个普遍的方法来将紧密充填的圆形颗粒体限制在一个任意封闭的面积内。这里来介绍两种方法如何充填特定半径的颗粒体达到我们所需要的空隙率。第一种方法，首先建立封闭区域的边界(简称墙体)，然后在封闭区域内任意生成一系列无接触的颗粒，最后移动区域的限制墙体，至所需要的空隙率。这种方法有三个缺点：1.区域的几何形状改变；2.收敛速度慢；3.最终的分布趋势是不均匀的目前没有一个普遍的方法来将紧密充填的圆形颗粒体限制56

第二种方法：运用GENERATE命令生成颗粒体，同时配合关键词高斯分配，即指定颗粒体半径的上下限，然后相应分配一个标准差，同时配合FISH函数来选择颗粒半径，最终生成我们所需要的模型。下面来介绍采用这种方法是如何进行操作的，重点介绍2种方法。第二种方法：运用GENERATE命令生成颗粒体，同573.2.1半径扩展法

首先需要介绍一个概念：半径放大系数m，它会引起模型空隙率的改变。空隙率的定义：

——颗粒体的总面积；——封闭区域的总面积因此：3.2.1半径扩展法首先需要介绍一个概念：半径放58

我们定义初始空隙率为，新生成的空隙率为,那么，则有：如果整个模型使用相同的半径放大系数，则：

或我们定义初始空隙率为，新生成的空隙率为59实例：区域宽：10单位区域高：5单位目标空隙率：0.12

颗粒体数目：300

最大最小半径比：1.5

初始假设一个m之值，可以求出初始的为：从而得到平均半径大小为：

实例：区域宽：10单位60

由初始给定的颗粒体最大直径和最小直径的比例r可得到颗粒体半径的上下限：

下限：

上限：

通过上面的计算，同时编写FISH函数来生成满足需要的空隙率的模型，并能自动计算空隙率的大小，以便用户进行检验。由初始给定的颗粒体最大直径和最小直径的比例r可得61newdefexpand;---输入数据---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^2r0=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：new程序2：62commandwallid=1ks=s_stiffkn=n_stiffnodes(_x0,_y0)(_x1,_y1)end_command_x0=width_y0=-extend*height_x1=width_y1=height*(1.0+extend)commandwallid=2ks=s_stiffkn=n_stiffnodes(_x0,_y0)(_x1,_y1)end_command_x0=width*(1.0+extend)_y0=height_x1=-extend*width_y1=heightcommandwallid=3ks=s_stiffkn=n_stiffnodes(_x0,_y0)(_x1,_y1)end_command

command63_x0=0.0_y0=height*(1.0+extend)_x1=0.0_y1=-extend*heightcommandwallid=4ks=s_stiffkn=n_stiffnodes(_x0,_y0)(_x1,_y1)end_command;---generatetheballsandgivethemtheirpropertiescommandgenid=1,numrad=rlo,rhix=0,widthy=0,heightpropdens=1000ks=s_stiffkn=n_stiffend_commandget_porosmult=sqrt((1.0-poros)/(1.0-pmeas))commandiniradmulmultcycle1000propfric0.2

_x0=0.064cycle2000end_commandenddefget_porossum=0.0bp=ball_headloopwhilebp#nullsum=sum+pi*b_rad(bp)^2bp=b_next(bp)end_looppmeas=1.0-sum/(width*height)endexpandget_porosplotwallballplotshowprintpmeassaveexpand.SAVcycle200065

模型接触力模型接触力66

第2种方法：指定颗粒体的半径，不限制颗粒的数目，使足够多的颗粒产生来达到所需要的空隙率。但这种方法所带来的缺点是可能在局部区域造成大面积的颗粒重叠，这将会产生很大的挤压力，从而给予颗粒较大的初始速度，这就可能使得颗粒体脱离墙体的限制。为避免此情况的发生，可通过初始的有限步循环计算将颗粒的动能减至零，然后再计算至平衡态。

3.2.2

挤压排斥法第2种方法：指定颗粒体的半径，不限制颗粒的数目，使67newsetrandom;设定颗粒体数目无限制模式defexploden_stiff=1e8;法向刚度s_stiff=1e8;切向刚度width=10.0height=5.0poros=0.12n_max=1000rlo=0.174rhi=0.261commandwallid1ks=s_stiffkn=n_stiffnodes(0,0)(width,0)wallid2ks=s_stiffkn=n_stiffnodes(width,0)(width,height)wallid3ks=s_stiffkn=n_stiffnodes(width,height)(0,height)wallid4ks=s_stiffkn=n_stiffnodes(0,height)(0,0)end_commandpvol_sum=0.0tot_vol=width*heightcount=0程序3：new程序3：68sectionloopn(1,n_max)r_ball=rlo+urand*(rhi-rlo)pvol_new=pvol_sum+pi*r_ball^2if(1.0-pvol_new/tot_vol)<porosthenexitsection;(newporositywillbe<thatspecified)end_ifrb2=r_ball*2.0x_ball=r_ball+urand*(width-rb2)y_ball=r_ball+urand*(height-rb2)commandballx=x_bally=y_ballrad=r_ballend_commandpvol_sum=pvol_newcount=count+1end_loopend_sectiontot_poros=1.0-pvol_sum/tot_volcommandpropdens=1000.0ks=s_stiffkn=n_stiffsection69end_commandii=out(''+string(count)+'ballswerecreated.')enddefquietloopn(1,4)commandcycle5inixvel0yvel0end_commandend_loopendexplodequietcyc980propfric=0.2cyc250plotwallballprinttot_porossaveexplode.SAVPFC中文详解教育课件70平衡前模型平衡后模型平衡前模型平衡后模型71速度向量图颗粒间接触速度向量图颗粒间接触72PFC中文详解教育课件73颗粒体间的连接关系图颗粒体间的连接关系图743.3不同尺度颗粒的分布

如果一个模型的不同区域需要不同尺寸的颗粒体分布时，我们就可以运用方法二(挤压排斥法)来生成目标模型。具体操作：将区域分成两部分，左侧随机充填大直径颗粒体，右册随机充填小直径颗粒体，左侧和右侧的区分通过x和y坐标来控制，最终生成后，通过平衡计算使两侧颗粒体自由扩散直到目标模型，详见下例。3.3不同尺度颗粒的分布如果一个模型的不同区域75newsetrandomsetdiskon;将球体看作单位厚度的圆盘状defsetupn_stiff=1e8s_stiff=1e8width=10.0height=5.0poros=0.12rat=1.5mult_0=1.6;initialradiusmultiplicationfactormult_a=0;左侧半径放大系数

mult_b=0;右侧半径放大系数

id1=0id2=0endsetup程序4：new程序4：76wallid1ks=s_stiffkn=n_stiffnodes(0,0)(width,0)wallid2ks=s_stiffkn=n_stiffnodes(width,0)(width,height)wallid3ks=s_stiffkn=n_stiffnodes(width,height)(0,height)wallid4ks=s_stiffkn=n_stiffnodes(0,height)(0,0);---------------------------------------------defmake_blocktot_vol=(x2-x1)*(y2-y1)num=id2-id1+1n0=1.0-(1.0-poros)/mult_0^2r0=sqrt(tot_vol*(1.0-n0)/(pi*num))rlo=2.0*r0/(1.0+rat)rhi=rat*rlocommand;Generatereduced-radiusparticlesgenid=id1,id2rad=rlo,rhix=x1,x2y=y1,y2propdens=1000ks=s_stiffkn=n_stiffrangeid=id1,id2end_commandget_porosmult=sqrt((1.0-poros)/(1.0-pmeas))endwallid1ks=s_stiffkn=n77defget_porossum=0.0bp=ball_headloopwhilebp#nullifb_id(bp)>=id1thenifb_id(bp)<=id2thensum=sum+pi*b_rad(bp)^2end_ifend_ifbp=b_next(bp)end_looppmeas=1.0-sum/tot_volenddeffinal_porostot_vol=width*heightid1=1id2=1200get_porosfinal_poros=pmeasenddefget_poros78setx1=0.0x2=5.0y1=0.0y2=5.0id1=1id2=50make_blocksetmult_a=multsetx1=5.0x2=10.0y1=0.0y2=5.0id1=1001id2=1200make_blocksetmult_b=multiniradmul=mult_ac_index0rangeid1,50iniradmul=mult_bc_index1rangeid1001,1200plocreatethe_assemblyplotaddballlgreenlorangeplotaddwallblackplotshowcycle1000propfric0.2cycle500printfinal_porossaveexpand2.SAVsetx1=0.0x2=5.0y1=0.079平衡前模型平衡后模型平衡前模型平衡后模型80法向接触力切向接触力法向接触力切向接触力81颗粒体间的连接关系图颗粒体间的连接关系图823.4

运用模型生成“过滤器”

有些情况需要我们建立复杂区域形状的颗粒流模型，如右图，此时我们可运用模型生成过滤器来获得所需要的模型，即filter命令，其后由用户定义

FISH函数来控制,其中,颗粒的半径通过fc_arg(0)进行检验,x和y

的坐标位置分别通过fc_arg(1)

和fc_arg(2)进行检验。如果颗粒满足要求，则FISH函数值设为

0，否则为1。详见下例。3.4运用模型生成“过滤器”有些情况需要我们建83newdefff_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)-_bradif_rx<ff_r1thenif_ry<ff_r1then_skip=1end_ifend_ifff_rect=_skipend程序5：new程序5：84defgen_balls_xlo=ff_x-ff_r2_xhi=ff_x+ff_r2_ylo=ff_y-ff_r2_yhi=ff_y+ff_r2commandgeneratex=(_xlo,_xhi)y=(_ylo,_yhi)&rad=(0.09,0.11)&filter=ff_rect&id=(1,250)end_commandendsetff_x=1.0ff_y=1.0ff_r1=2.0ff_r2=3.0gen_ballspropertydens=1000kn=1e8ks=1e8wallid=1nodes(-1.0,-1.0)(-1.0,3.0)(3.0,3.0)(3.0,-1.0)closewallid=2nodes(-2.0,-2.0)(4.0,-2.0)wallid=3nodes(4.0,-2.0)(4.0,4.0)

defgen_balls85

wallid=4nodes(4.0,4.0)(-2.0,4.0)wallid=5nodes(-2.0,4.0)(-2.0,-2.0)wallid=1kn=1e8ks=1e8wallid=2kn=1e8ks=1e8wallid=3kn=1e8ks=1e8wallid=4kn=1e8ks=1e8wallid=5kn=1e8ks=1e8plotcreatethe_viewplotaddballyellowplotaddaxesblackplotaddwallblueid=onplotshowpausepropertyradmul1.5plotaddcfgreencycle500wallid=4nodes(4.0,4.86PFC中文详解教育课件87PFC中文详解教育课件884.边界条件

PFC2D中有三种边界条件，分别是：墙体边界、颗粒体边界和混合边界。其中颗粒体边界又分为速度边界和受力边界。

1).墙体边界建模过程中，墙体可作为颗粒体的生成范围约束，但同时也可以将墙体作为边界来施加约束。对于墙体，我们只能施加速度约束，而不能直接对其施加外力，因为运动定律对墙体是不适用的。其速度由以下三个参数控制:线速度、角速度和旋转中心。墙的运动是通过不断更新定义墙的基点的位置来描述。采用WALL命令设置，如：

wallid=1x=1.0y=1.0spin=10.04.边界条件PFC2D中有三种边界条件，分别是：墙89(2)颗粒体边界

PFC2D中的模型可以将一连串的颗粒体作为边界条件。基本方法：在模型紧密压缩至平衡后，我们通过FISH函数将与墙体相接触的颗粒体逐个提取，将这一系列的颗粒体采用共同的边界条件限制，最后删除初始的限制墙体，即实现了以颗粒体代替墙体来作为边界条件。颗粒体边界条件分速度边界和外力边界，下例为速度边界程序实现。(2)颗粒体边界90restoreexpand.savdefboundbp=ball_headloopwhilebp#nullsectioncp=b_clist(bp)loopwhilecp#nullifc_nforce(cp)#0.0thenb2=c_ball2(cp)ifpointer_type(b2)=101then;b2isawallb_xfix(bp)=1;fixoriginalballinx,yb_yfix(bp)=1b_color(bp)=1exitsection;alldoneforthisballend_ifend_ififc_ball1(cp)=bpcp=c_b1clist(cp)elsecp=c_b2clist(cp)程序6：restoreexpand.sav程序6：91end_ifend_loopend_sectionbp=b_next(bp)end_loopendboundinixvel=0.0grad=(-0.1,0)iniyvel=0.0delwall1delwall2delwall3delwall4plotcreateassemblyplotaddballredgreenplotaddvelblackplotshowplotaddcfwhitecycle200end_if92放大放大

颗粒体速度边界示意图放大放大颗粒体速度边界示意图93

其前期操作与速度边界类似，差异处是先将边界颗粒体的速度初始化为零，再删除限制墙体，采用FISH函数反向施加颗粒体所受的不平衡力，运行一定的计算步至平衡态。在此过程中，区域内部的颗粒体会因轻微的扰动而导致其自身产生运动而脱离边界颗粒体，此时前面将边界颗粒体的速度固定为零。颗粒体外力边界其前期操作与速度边界类似，差异处是先将边界颗粒体的94inixvel=0yvel=0spin=0delwall1delwall2delwall3delwall4cycle1defreplace;Replaceunbalancedforceswithappliedforcesbp=ball_headloopwhilebp#nullifb_xfix(bp)=1b_xfap(bp)=-b_xfob(bp)b_yfap(bp)=-b_yfob(bp)end_ifbp=b_next(bp)end_loopendreplacefreexyspincycle200程序7：inixvel=0yvel=0spin=0程序7：95(3)混合边界条件

混合边界即速度边界与外力边界同时存在，对于双轴压缩试验，其两侧端面不允许存在明显的几何变形，否则其边界力将会无效。因此，此时需要用到混合边界，两侧端面采用外力边界来约束，确保其不发生几何变形，侧面Y向采用速度边界，使其速度线性增加，以此来模拟试块轴向变形。(3)混合边界条件混合边界即速度边界与外力边界965.初始条件

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

cp=b_clist(bp)loopwhilecp#nullifc_nforce(cp)#0.0then_rcp=sqrt((c_x(cp)-b_x(bp))^2+(c_y(cp)-b_y(bp))^2)if_rcp<=b_rad(bp)then;刚好接触或重叠

ifc_ball1(cp)=bpthenbp_other=c_ball2(cp)elsebp_other=c_ball1(cp)end_if

程序85.1施加等向应力defexpand_radii程序898ifpointer_type(bp_other)=101then;21,100,101,102,103or104

_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)*_phiend_ifend_ififc_ball1(cp)=bpthencp=c_b1clist(cp)elsecp=c_b2clist(cp)end_ifend_loopbp=b_next(bp)end_loopifpointer_type(bp_other)=1099

_alpha=-1.0*_lambda*tot_vol*_diso/_sumbp=ball_headloopwhilebp#nullb_rad(bp)=(1.0+_alpha)*b_rad(bp)bp=b_next(bp)end_loopenddefachieve_isostrloopwhile1#0;始终执行，除非遇到exit命令

_diso=req_isostr-isostrifabs(_diso/req_isostr)<=req_isostr_tolthenexitend_ifoo=out('Currentisotropicstress='+string(isostr))expand_radiicommandsolve

_alpha=-1.0*_lambda*tot100end_commandend_loopenddefisostri