PFC培训PPT.ppt

PFC2D颗粒流软件培训中仿科技有限公司 主要内容 第一部分PFC2D颗粒流程序简介第二部分有限差分法基础介绍第三部分离散元法基础介绍第四部分PFC2D的使用 第一部分PFC2D颗粒流程序简介 1 理论背景2 颗粒流方法的基本假设3 颗粒流方法的特点4 可选特性5 应用领域6 求解步骤 作为离散元的一种 二维颗粒流程序 ParticleFollowCodePFC2D 数值模拟新技术 其理论基础是Cundall 1979 提出的离散单元法 用于颗粒材料力学性态分析 如颗粒团粒体的稳定 变形及本构关系 专门用于模拟固体力学大变形问题 它通过圆形 或异型 离散单元来模拟颗粒介质的运动及其相互作用 由平面内的平动和转动运动方程来确定每一时刻颗粒的位置和速度 作为研究颗粒介质特性的一种工具 它采用有代表性的数百个至上万个颗粒单元 通过数值模拟实验可以得到颗粒介质本构模型 1 理论背景 PFC2D ParticleFollowCode2Dimension 即二维颗粒流程序 是通过离散单元方法来模拟圆形颗粒介质的运动及其相互作用 最初 这种方法是研究颗粒介质特性的一种工具 它采用数值方法将物体分为有代表性的数百个颗粒单元 期望利用这种局部的模拟结果来研究边值间题连续计算的本构模型 以下两种因素促使PFC2D方法产生变革与发展 1 通过现场实验来得到颗粒介质本构模型相当困难 2 随着微机功能的逐步增强 用颗粒模型模拟整个问题成为可能 一些本构特性可以在模型中自动形成 因此 PFC2D便成为用来模拟固体力学和颗粒流问题的一种有效手段 2 颗粒流方法的基本假设 颗粒流方法在模拟过程中作了如下假设 1 颗粒单元为刚性体 2 接触发生在很小的范围内 即点接触 3 接触特性为柔性接触 接触处允许有一定的 重叠 量 4 重叠 量的大小与接触力有关 与颗粒大小相比 重叠 量很小 5 接触处有特殊的连接强度 6 颗粒单元为圆盘形 或球形 其中 颗粒为刚性体的假设 对于模拟介质运动为只沿相互接触面的表面发生的问题非常重要 比如象砂土或粮食这种颗粒组合体材料 利用这种假设在总体上来讲是比较恰当的 因为这种材料的变形是来自于颗粒刚性体间的滑动和转动以及接触面处的张开和闭锁 而不是来自于每个刚性颗粒本身的变形 对于这种特殊材料 没有必要采用非常精确的数值模型 来得到对材料特性的近似 3 颗粒流方法的特点 PFC2D可以直接模拟圆形颗粒的运动和相互作用问题 颗料可以代表材料中的个别颗粒 例如砂粒 也可以代表粘结在一起的固体材料 例如混凝土或岩石 当粘结以渐进的方式破坏时 它能够破裂 粘结在一起的集合体可以是各向同性 也可以被分成一些离散的区域或块体 这类物理系统可以用处理角状块体的离散单元程序UDEC和3DEC来模拟 PFC2D有三个优点 第一 它有潜在的高效率 因为圆形物体间的接触探测比角状物体间的更简单 第二 对可以模拟的位移大小实质上没有限制 第三 由于它们是由粘结的粒子组成 块体可以破裂 不象UDEC和3DEC模拟的块体不能破裂 用PFC2D模拟块体化系统的缺点是 块体的边界不是平的 用户必须接受不平的边界以换取PFC2D提供的优点 PFC2D中几何特征 物理特性和解题条件的说明不如FLAC和UDEC程序那样直截了当 例如用连续介质程序 创建网格 设置初始压力 设置固定或自由边界 在象PFC2D这样的颗粒程序中 由于没有唯一的方法在一个指定的空间内组合大量的粒子 粒子紧密结合的状态一般不能预先指定 必须跟踪类似于物体压实的过程 直到获得要求的孔隙率 由于颗粒相对位置变化产生接触力 初始应力状态的确定与初始压密有关 由于边界不是由平面组成 边界条件的设定比连续介质程序更复杂 当要求满足有实验室实际测试的模拟物体的力学特性时 出现了更大的困难 在某种程度上 这是一个反复试验的过程 因为目前还没有完善的理论可以根据微观特性来预见宏观特性 然而 给出一些准则应该有助于模型与原型的匹配 如哪些因素对力学行为的某些方面产生影响 哪些将不产生影响 应该意识到 由于受现有知识的限制 这样的模拟很难 然而 用PFC2D进行试验 对固体力学 特别是对断裂力学和损伤力学 可以获得一些基本认识 PFC2D能模拟任意大小圆形粒子集合体的动态力学行为 粒子生成器根据粒子的指定分布规律自动概率地生成 粒子半径按均匀分布或按高斯分布规律分布 初始孔隙度一般比较高 但通过控制粒子半径的扩大可以获得密度压实 在任何阶段任何因素都可以改变半径 所以不需反复试验就可以获得指定孔隙度的压实状态 属性与各个粒子或接触有关 而不是与 类型号 有关 因此 可以指定属性和半径的连续变化梯度 节理生成器 用来修改沿指定轨迹线的接触特性 假定这些线叠加在颗粒集合体上 用这种方法 模型可以被成组的弱面 如岩石节理切割 粒子颜色也是一种属性 用户可以指定各种标记方案 PFC2D模型中为了保证数据长期不漂移 用双精度数据存储坐标和半径 接触的相对位移直接根据坐标而不是位移增量计算 接触性质由下列单元组成 1 线性弹簧或简化的Hertz Mindlin准则 2 库仑滑块 3 粘结类型 粘结接触可承受拉力 粘结存在有限的抗拉和抗剪强度 可设定两种类型的粘结 接触粘结和平行粘结 这两种类型粘结对应两种可能的物理接触 接触粘结再现了作用在接触点一个很小区域上的附着作用 平行粘结再现了粒子接触后浇注其它材料的作用 如水泥灌浆 平行粘结中附加材料的有效刚度具有接触点的刚度 块体逻辑支持附属粒子组或块体的创建 促进了程序的推广普及 块体内粒子可以任意程度的重叠 作为刚性体具有可变形边界的每一个块体 可作为一般形状的超级粒子 通过指定墙的速度 混合的粒子速度 施加外力和重力来给系统加载 扩展的FISH库 提供了在集合体内设置指定应力场或施加应力边界条件的函数 时步计算是自动的 包括因为Hertz接触模型刚度变化的影响 模拟过程中 根据每个粒子周围接触数目和瞬间刚度值 时步也在变化 基于估计的粒子数 单元映射策略采用最佳的单元数目 自动调整单元的外部尺寸来适应粒子缺失和指定的新对象 单元映射方案支持接触探测算法以保证求解时间随粒子数目线性增加 而不是二次方增加 类似于FLAC PFC提供了局部无粘性阻尼 这种阻尼形式有以下优点 1 对于匀速运动 体力接近于零 只有加速运动时才有阻尼 2 阻尼系数是无因次的 3 因阻尼系数不随频率变化 集合体中具有不同自然周期的区域被同等阻尼 采用同样的阻尼系数 PFC2D可以在半静态模式下运行以保证迅速收敛到静态解 或者在完全动态模式下运行 PFC2D包含功能强大的内嵌式程序语言FISH 允许用户定义新的变量和函数使数值模型适合用户的特殊需求 例如 用户可以定义特殊材料的模型和性质 加载方式 实验条件的伺服控制 模拟的顺序以及绘图和打印用户定义的变量等 4 可选特性 1 热学分析2 并行处理技术3 能写用户定义接触模型4 用户写C 程序的C 编程 热学选项用来模拟材料内热量的瞬间流动和热诱导位移和力的顺序发展 热学模型可以独立运行或耦合到力学模型 通过修改粒子半径和平行粘结承受的力 产生热应变来解释粒子和粘结材料的受热 用户定义的接触本构模型可以用C 语言来编写 并编译成动态链接库文件 一旦需要就可以加载 用户写的C 程序选项允许用户用C 语言写自己的程序 创建可执行的PFC2D个人版本 这个选项可以用来代替FISH函数 大大提高运行的速度 并行处理技术允许将一个PFC2D模型分成几个部分 每个部分可以在单独的处理器上平行运行 与一个PFC2D模型在一个处理器上运行相比 平行处理在内存容量和计算速度方面得到大大提高 5 应用领域 PFC2D既可解决静态问题也可解决动态问题 既可用于参数预测 也可用于在原始资料详细情况下的实际模拟 PFC2D模拟试验可以代替室内试验 在岩石与土体中开挖问题的研究与设计方面 实测资料相对较少 关于初始应力 不连续性等问题也只能部分了解 而在松散介质流动问题中 影响流动介质不规律分布的影响因素很难定量描述 因此 应用PFC2D初步研究影响整个系统的一些参数的特性 对整个系统的特性有所了解后 就可以方便地设计模型模拟整个过程 PFC2D可以模拟颗粒间的相互作用问题 大变形问题 断裂问题等 适用于以下领域 1 在槽 管 料斗 筒仓中松散物体的流动问题 2 矿山冒落法开采中的岩体断裂 坍塌 破碎和岩块的流动问题 3 铸模中粉料的压实问题 4 由粘结粒子组成物体的碰撞及其动态破坏 5 梁结构的地震响应及垮塌 6 颗粒材料的基本特性研究 如屈服 流动 体积变化等 7 固体的基本特性研究 如累积破坏 断裂 6 求解步骤 1 定义模拟对象根据模拟意图定义模型的详细程序 假如只对某一力学机制的不同解释作出判断时 可以建立一个比较粗略的模型 只要在模型中能体现要解释的机制即可 对所模拟问题影响不大的特性可以忽略 2 建立力学模型的基本概念首先对分析对象在一定初始特性形成初步概念 为此 应先提出一些问题 如系统是否将变为不稳定系统 问题变形的大小 主要力学特性是否非线性 是否需要定义介质的不连续性 系统边界是实际边界还是无限边界 系统结构有无对称性等 综合以上内容来描述模型的大致特征 包括颗粒单元的设计 接触类型的选择 边界条件的确定以及初始平衡状态的分析 3 构造并运行简化模型在建立实际工程模型之前 先构造并运行一系列简化的测试模型 可以提高解题效率 通过这种前期简化模型的运行 可对力学系统的概念有更深入的了解 有时在分析简化模型的结果后 例如所选的接触类型是否有代表性 边界条件对模型结果的影响程度等 还需将第二步加以修改 4 补充模拟问题的数据资料模拟实际工程问题需要大量简化模型运行的结果 对于地质力学来说包括 a 几何特性 如地下开挖酮室的形状 地形地貌 坝体形状 岩土结构等 b 地质构造位置 如断层 节理 层面等 c 材料特性 如弹 塑性 后破坏特性等 d 初始条件 如原位应力状态 孔隙压力 饱和度等 e 外荷载 如冲击荷载 开挖应力等 因为一些实际工程性质的不确定性 特别是应力状态 变形和强度特性 所以必须选择合理的参数研究范围 第三步简化模型的运行有助于这项选择 从而为更进一步的试验提供资料 5 模拟运行的进一步准备a 合理确定每一时步所需时间 若运行时间过长 很难得到有意义的结论 所以应该考虑在多台计算机上同时运行 b 模型的运行状态应及时保存 以便在后续运行中调用其结果 例如如果分析中有多次加卸荷过程 要能方便地退回到每一过程 并改变参数后可以继续运行 c 在程序中应设有足够的监控点 如参数变化处 不平衡等 对中间模拟结果随时作出比较分析 并分析颗粒流动状态 6 运行计算模型在模型正式运行之前先运行一些检验模型 然后暂停 根据一些特性参数的试验或理论计算结果来检查模拟结果是否合理 当确定模型运行正确无误时 连接所有的数据文件进行计算 7 解释结果计算结果与实测结果进行分析比较 图形应集中反应要分析的区域如应力集中区 各种计算结果应能方便地输出分析 第二部分有限差分法基础介绍 连续介质三维快速拉格朗日有限差分计算方法 FLAC3D 是近20年来逐步成熟完善起来的一种新型数值计算方法 它基于显式差分法来求解运动方程和动力方程 可模拟岩土或其他材料的三维力学行为 其求解时首先将计算区域离散化 分成若干三维单元 单元之间由节点联结 节点受荷载作用后 其平衡方程 运动方程 可以写成时间步长为 t的有限差分形式 由于采用动态应力松弛显式差分求解技术 在某一微小的时段内 作用于该节点的荷载只对周围若干节点有影响 根据单元节点的速度变化和时段 t 可求出单元之间的相对位移 进而求出单元应变 利用单元材料的本构关系即可求出单元应力 在此基础上 求出单元之间的不平衡力 将此不平衡力重新作用到节点上 再进行下一步的迭代过程 直到整个系统不平衡力足够小或节点位移趋于平衡为止 FLAC3D可以解决诸多的有限元程序难以模拟的复杂的工程问题 例如分布开挖 大变形 非线性及非稳定系统 甚至大面积屈服 失稳或完全塌方 第三部分离散元法基础介绍 离散单元法是一种模拟非连续介质的计算方法 自Cundall在70年代提出以来 在岩石力学 土力学 结构分析等领域的数值模拟中得到广泛应用 是一种新兴的非连续体分析方法 离散单元法允许单元间的相对运动 不一定满足位移连续和变形协调条件 计算速度快 所需存储空间小 特别适用于节理岩体的大位移 大变形分析 离散单元法自问世以来有了长足的发展 已经成为解决岩石力学问题的一种重要的数值方法 因为工程中所见到的岩体其形态呈非连续结构 所形成的岩石块体运动和受力情况多是几乎或材料非线性问题 所以很难用解决连续介质力学问题的有限单元法或边界单元法等 数值方法来进行求解 而离散单元法正是充分考虑到岩体结构的不连续性 适用于解决节理岩石力学问题 近年来 离散元法的应用领域又扩展到求解连续介质向非连续介质转化的力学问题 混凝土等脆性材料在冲击 侵彻等动荷载作用下产生的损伤和破坏 其实质是力学模型从连续体到非连续体的转变过程 建立在传统的连续介质力学基础上的有限元法等数值计算方法难以直接用于计算和模拟材料具体的破坏形式和破坏的整个过程 而离散元法在这一方面显示出巨大的生命力 第四部分PFC2D的使用 1 对PFC软件的使用界面 菜单功能及作用进行介绍 2 FISH语言简介3 PFC2D分析模型的生成方法4 边界条件的设置方法5 初始条件的设置6 接触本构模型 接触刚度模型 滑动模型 连接模型7 赋予材料属性 相关命令的使用方法介绍 第四部分PFC2D的使用 8 节理面的生成及属性设置9 加载方法 主动荷载和被动荷载 10 求解过程 静力求解 动力求解 11 流体与热分析简介 12 介绍PFC2D软件的用户自定义本构模块的相关功能 操作等 13 常用命令使用方法及相关的重要概念 14 讲述PFC2D工程应用的实例 1 使用界面 菜单功能介绍 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 DEFabcabc 22 3 5ENDPrintabc 对上例子稍作改进 例2 newdefabchh 22abc hh 3 5end1 稍有编程常识的人可以看出 执行上面的例子 PRINTabc 其结果与例1相同 abc 71 在这个函数中 我们首先把22赋值给变量hh 然后把这个变量带入abc的表达式中 因此二者的结果相同 2 FISH的执行过程如下 当在程序命令中使用一个FISH符号名时 例如执行PRINT符号名 如果符号名也是一个函数名 那么执行这个函数 例如abc 如果符号名不是函数名 那么使用符号目前的值 例如hh 3 在输入完例2的各行后 如果我们执行命令 PRINThh 此时hh 0 因为在这个时候没有执行FISH函数 因此hh的初始值为0 我们接着执行PRINTabc 结果显示abc 71 再次执行PRINThh 此时结果为hh 22 这是因为我们首先运行了abc函数 在这个过程中hh已被赋值 4 下面的试验将进一步解释函数与变量之间的差别 注意 Itasca软件的SET命令可以用来设置任何用户定义的FISH符号的值 与在FISH中使用的符号无关 下面的例3建立在例2的基础之上 我们不使用NEW命令来清除内存中的值 因为我们想继续使用那些值 setabc 0hh 0printhhprintabcprinthh 5 在这个例子中 我们首先把abc和hh都赋值为0 由于hh是一个变量 第一个Print命令显示当前hh的值 hh 0 第二个Print命令由于abc是一个函数名 因此执行abc函数 先前定义的abc 0不起作用 重新计算了hh和abc的值 因此第三个Print命令显示的值是它在abc函数内指定的值 即hh 22 例4是这个试验完整的命令 现在我们总结一下 Itasca软件的三个重要的命令PRINT SET HISTORY可以直接操作简单的FISH变量或函数 如下图所示 其中var代表变量名或函数名 HISTORY命令的用法在这不作重复介绍 即定义了FISH函数后 对其中的变量可以进行跟踪 即HISTTORYVar 象其它高级编程语言一样 FISH有执行循环命令的功能 标准的格式如下 LOOPvar expr1 expr2 END LOOP其中LOOP和END LOOP是FISH语句 符号var代表循环变量 expr1和pxpr2代表表达式或者单个变量 下面的例子用循环命令计算从1到10的和以及乘积 见下例 newdefabcsum 0prod 1loopn 1 10 sum sum nprod prod nend loopendabcprintsum prod在这个例子中 首先给两个变量赋于初始值 sum用来保存和的结果 prod用来保存积的结果 然后执行循环 最后分别打印出这两个变量的最后结果 循环变量n 1 10 表示从1开始 连续计算到10结束 关于LOOP的注意事项 1 FISH接受END LOOP和ENDLOOP的写法 但不接受ENDLOOP这样中间有空格的写法 其它类似的命令有着同样的规则 如END IF END COMMAND等命令 2 在上面的例子中 如果执行Printn或Printfish命令 你会看到n 11而不是10 注意 这不是FISH的错误 这是一个基本的计算机指令存储规则 当循环结束后 计数器的值保存的是n 1而不是n 所有的高级编程语言有着相同的规则 1 DEFINEfunctionEND2 CASEOFexprCasenendcase3 IFexpr1testexpr2THENELSEENDIF4 LOOPvar expr1 expr2 ENDLOOP 5 LOOPWHILEexpr1testexpr2ENDLOOP6 COMMANDENDCOMMAND7 HISTORYvarPRINTvarSETvarvaluePLOTadd shfname 另外 在FISH中还有许多其它的预定义对象 其中一类是尺度变量 scalarvariables 它们是单个的数字 下面是总的尺度变量 clock 时钟时间 单位是秒的100倍 unbal 最大不平衡力pi 圆周率step 目前的时步数目urand 0 0 1 0之间均匀分布的随机变量这仅是其中的一小部分 完全的列表以后再述 另一类非常有用的内置对象是固有函数 intrinsicfunctions 这些函数能在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命令 循环生成一系列的颗粒 如下例 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 r2 end loopyc yc yincxc x0 radius row row 2 2 end loopendsetechooffsetx0 0 2y0 0 4radius 0 1setid start 100n col 7n row 8hexsetechoonplotsetcapsize20plotaddaxesblackplotaddballyellowplotshow 3 2不规则排列 无规则排列 即 对一个给定空隙率的区域 采用颗粒来充填其中需要进行填充的空隙 并确保整个模型保持平衡 对于所能被填充的模型的初始空隙率 是有一个限制值 不能任意小 对于某些空隙率的模型 颗粒的填充可以无接触地排列 对于其它情况的空隙率 颗粒又可以重叠排列 目前没有一个普遍的方法来将紧密充填的圆形颗粒体限制在一个任意封闭的面积内 这里来介绍两种方法如何充填特定半径的颗粒体达到我们所需要的空隙率 第一种方法 首先建立封闭区域的边界 简称墙体 然后在封闭区域内任意生成一系列无接触的颗粒 最后移动区域的限制墙体 至所需要的空隙率 这种方法有三个缺点 1 区域的几何形状改变 2 收敛速度慢 3 最终的分布趋势是不均匀的 第二种方法 运用GENERATE命令生成颗粒体 同时配合关键词高斯分配 即指定颗粒体半径的上下限 然后相应分配一个标准差 同时配合FISH函数来选择颗粒半径 最终生成我们所需要的模型 下面来介绍采用这种方法是如何进行操作的 重点介绍2种方法 3 2 1半径扩展法 首先需要介绍一个概念 半径放大系数m 它会引起模型空隙率的改变 空隙率的定义 颗粒体的总面积 封闭区域的总面积因此 我们定义初始空隙率为 新生成的空隙率为 那么 则有 如果整个模型使用相同的半径放大系数 则 或 实例 区域宽 10单位区域高 5单位目标空隙率 0 12颗粒体数目 300最大最小半径比 1 5初始假设一个m之值 可以求出初始的为 从而得到平均半径大小为 由初始给定的颗粒体最大直径和最小直径的比例r可得到颗粒体半径的上下限 下限 上限 通过上面的计算 同时编写FISH函数来生成满足需要的空隙率的模型 并能自动计算空隙率的大小 以便用户进行检验 newdefexpand 输入数据 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 commandwallid 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 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 cycle2000end 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 SAV 模型 接触力 第2种方法 指定颗粒体的半径 不限制颗粒的数目 使足够多的颗粒产生来达到所需要的空隙率 但这种方法所带来的缺点是可能在局部区域造成大面积的颗粒重叠 这将会产生很大的挤压力 从而给予颗粒较大的初始速度 这就可能使得颗粒体脱离墙体的限制 为避免此情况的发生 可通过初始的有限步循环计算将颗粒的动能减至零 然后再计算至平衡态 3 2 2挤压排斥法 newsetrandom 设定颗粒体数目无限制模式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 sectionloopn 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 0 x 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 stiff end commandii out string count ballswerecreated enddefquietloopn 1 4 commandcycle5inixvel0yvel0end commandend loopendexplodequietcyc980propfric 0 2cyc250plotwallballprinttot porossaveexplode SAV 平衡前模型 平衡后模型 速度向量图 颗粒间接触 颗粒体间的连接关系图 3 3不同尺度颗粒的分布 如果一个模型的不同区域需要不同尺寸的颗粒体分布时 我们就可以运用方法二 挤压排斥法 来生成目标模型 具体操作 将区域分成两部分 左侧随机充填大直径颗粒体 右册随机充填小直径颗粒体 左侧和右侧的区分通过x和y坐标来控制 最终生成后 通过平衡计算使两侧颗粒体自由扩散直到目标模型 详见下例 newsetrandomsetdiskon 将球体看作单位厚度的圆盘状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 wallid1ks 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 end defget 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 pmeasend setx1 0 0 x2 5 0y1 0 0y2 5 0id1 1id2 50make blocksetmult a multsetx1 5 0 x2 10 0y1 0 0y2 5 0id1 1001id2 1200make blocksetmult b multiniradmul mult ac index0rangeid1 50iniradmul mult bc index1rangeid1001 1200plocreatethe assemblyplotaddballlgreenlorangeplotaddwallblackplotshowcycle1000propfric0 2cycle500printfinal porossaveexpand2 SAV 平衡前模型 平衡后模型 法向接触力 切向接触力 颗粒体间的连接关系图 3 4运用模型生成 过滤器 有些情况需要我们建立复杂区域形状的颗粒流模型 如右图 此时我们可运用模型生成过滤器来获得所需要的模型 即filter命令 其后由用户定义FISH函数来控制 其中 颗粒的半径通过fc arg 0 进行检验 x和y的坐标位置分别通过fc arg 1 和fc arg 2 进行检验 如果颗粒满足要求 则FISH函数值设为0 否则为1 详见下例 newdefff 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 defgen 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 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 5plotaddcfgreencycle500 4 边界条件 PFC2D中有三种边界条件 分别是 墙体边界 颗粒体边界和混合边界 其中颗粒体边界又分为速度边界和受力边界 1 墙体边界建模过程中 墙体可作为颗粒体的生成范围约束 但同时也可以将墙体作为边界来施加约束 对于墙体 我们只能施加速度约束 而不能直接对其施加外力 因为运动定律对墙体是不适用的 其速度由以下三个参数控制 线速度 角速度和旋转中心 墙的运动是通过不断更新定义墙的基点的位置来描述 采用WALL命令设置 如 wallid 1x 1 0y 1 0spin 10 0 2 颗粒体边界PFC2D中的模型可以将一连串的颗粒体作为边界条件 基本方法 在模型紧密压缩至平衡后 我们通过FISH函数将与墙体相接触的颗粒体逐个提取 将这一系列的颗粒体采用共同的边界条件限制 最后删除初始的限制墙体 即实现了以颗粒体代替墙体来作为边界条件 颗粒体边界条件分速度边界和外力边界 下例为速度边界程序实现 restoreexpand 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 end ifend loopend sectionbp b next bp end loopendboundinixvel 0 0grad 0 1 0 iniyvel 0 0delwall1delwall2delwall3delwall4plotcreateassemblyplotaddballredgreenplotaddvelblackplotshowplotaddcfwhitecycle200 放大 放大 颗粒体速度边界示意图 其前期操作与速度边界类似 差异处是先将边界颗粒体的速度初始化为零 再删除限制墙体 采用FISH函数反向施加颗粒体所受的不平衡力 运行一定的计算步至平衡态 在此过程中 区域内部的颗粒体会因轻微的扰动而导致其自身产生运动而脱离边界颗粒体 此时前面将边界颗粒体的速度固定为零 颗粒体外力边界 inixvel 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 3 混合边界条件 混合边界即速度边界与外力边界同时存在 对于双轴压缩试验 其两侧端面不允许存在明显的几何变形 否则其边界力将会无效 因此 此时需要用到混合边界 两侧端面采用外力边界来约束 确保其不发生几何变形 侧面Y向采用速度边界 使其速度线性增加 以此来模拟试块轴向变形 5 初始条件 首先介绍一下应力的概念 在连续介质中 应力是一个连续量 应力一般定义为作用于单位面积上的力 对于二维问题 应力则定义为作用于单位线段上的合力 然而 对于松散介质这种定义是不合适的 因为颗粒介质是离散的 在颗粒模型中没有存在于每一个点上连续的应力 而且变化幅度很大 对于应变也存在这样的问题 在PFC模型中 计算颗粒间接触力和颗粒的位移 这些量对于在微观上研究材料的特性很有意义 但不能将这些量直接联系到连续模型中 需要经过一个平均的过程 才能从微观传到连续模型中 所以 在PFC2D中采用平均应力的概念来表示 5 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 程序8 ifpointer type bp other 101then 21 100 101 102 103or104 phi b rad bp memory ball wall contact clump circleelse 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 loop 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 end commandend loopenddefisostrii measure mp 1 计算应力1表示应力 2表示应变isostr m s11 mp m s22 mp 2 0enddefget mp 返回measureid 1mp find meas 1 endmeasid 1x5y2 5rad2 0get mphistoryid 1isostrSETreq isostr 4e5req isostr tol 0 01cyc20printmeas1printisostrachieve isostr 6 接触本构模型 在PFC2D中 材料的本构特性是通过接触本构模型来模拟的 每一颗粒的接触本构模型有 1 接触刚度模型 2 滑动模型 3 连接模型 接触刚度模型提供了接触力和相对位移的弹性关系 滑动模型则强调切向和法向接触力使得接触颗粒可以发生相对移动 而连接模型是限制总的切向和法向力使得在连接强度范围内发生接触 1 接触刚度模型接触刚度通过以下两式把接触力与相对位移联系起来 即 式中 是法向刚度 它为割线刚度 联系总的法向力与位移 式中 是切向刚度 它为切线刚度 通过增量形式联系切向力与位移 在PFC2D中有两种接触刚度模型 即线性接触刚度模型与简化的Hertz Minlin接触刚度模型 不同的接触刚度模型有不同的接触刚度值 线性接触模型 分析颗粒间接触特性时 将这种接触的颗粒想象为一梁端点在颗粒中心弹性梁 梁端部受力或力矩相当于作用于颗粒的中心 这种梁由以下特征参数描述 1 几何参数 长度