使用FLAC自带FISH语言安装锚杆实例

合集下载

FLAC--FISH(最新)

卷首语要学好FLAC3D，必须学好FISH，FISH身为FLAC3D的内嵌语言，也就是是一种嵌入式编程语言，其编写短小精干实时，同时兼顾大家操作习惯，不但可以嵌入命令流文件里工作，而且还可以引用FLAC3D本身的任何命令，突破了一般标准程序代码的限制，实现了读者对FLAC3D软件的完美控制。

类似于ANSYS的APDL语言，FISH主要是用来处理现有命令程序难以实现（较难或较烦）的一些问题，对于没有编程经验的用户，可以编一些简单的FISH函数，最好是修改3DSHOP中的一些比较现成的函数以便使用；当然，既然是一门编程语言，也可以编复杂程序。

写这一卷的主要目的就是使大家对FISH有个较为充分的认识和理解，因为相关资料甚少，故卷中的有些资料是作者在Simwe,Yantubbs等论坛上收集后整理所得，在此表示感谢。

本卷分为4章，第一章为关于FISH语言法则、变量、函数等得一个综述，主要是想让大家从全局认识FISH编程语言；第二章介绍了FISH语句的类别和各自特点和用法；第三章则较为细致的向大家介绍了各变量，函数的具体含义，部分配合了用法实例；第四章补充了一些比较不错的FISH应用实例。

希望大家读后会有所收获，也欢迎大家批评指正，共同提高！Sunshinessw1216@QQ 61925252007-9-20于铁道科学研究院北京·中国第一章 Fish的语言法则一、四个概念函数与变量—Fish函数由define定义，end结束定义。

如下例为定义一个名叫abc 的函数。

上述函数中hh即为“变量”。

函数和变量是FISH语言中非常基本和重要的两个对象，区别如下：当我们调用一个FISH对象XXX时，如果XXX是函数，系统就会执行该函数；而如果XXX不是函数而是变量，系统则会简简单单的调用其当前置。

大家注意继续输以下命令行，注意输出结果，便可对这两个有个较为清晰的认识。

语句—在不少实际问题中有许多具有指示说明性的、控制性的、重复性操作，程序中需要执行他们，诸如Array, Define, end, Case, endcase, If ,else, endif 等即为语句。

FLAC_3D的锚杆拉拔数值模拟试验.kdh

第41卷第10期2009年10月哈尔滨工业大学学报JOURNAL OF HARBIN INSTITUTE OF TECHNOLOGYVol.41No.10Oct.2009FLAC－3D的锚杆拉拔数值模拟试验江文武1，2，徐国元1，马长年1（1.中南大学资源与安全工程学院，长沙410083，wenwujiang@；2.江西理工大学江西赣州341000）摘要：为研究锚杆锚固力主要影响因素、锚杆拉拔过程中的整体失稳的规律等，采用三维显式有限差分法，建立锚杆拉拔数值仿真模型，进行一系列计算机模拟.结果表明：数值模拟计算的结果和现场试验得到的结果基本吻合，表明数值模拟锚杆拉拔过程是可行的；锚杆拔出的过程是慢慢滑移到突然整体失稳过程；锚杆锚固效应随着锚固剂的摩擦角、粘结力和有效围压的增加而增强；锚固剂所受的剪切应力分布规律随锚杆拉拔过程而改变，在拉拔初始时，自由端至锚固端，自由端锚固剂处的剪切应力为最大，锚固端处的剪切应力为最小接近为零，随着锚杆拉拔的进行，自由端处锚固剂与锚杆的界面屈服点首先达到，造成自由端处锚固剂与锚杆出现滑移现象，而随着锚杆拉拔的进一步进行，锚固剂所受的剪切应力慢慢呈现均匀分布，均都达到了最大值，均达到了锚固剂与锚杆界面的屈服点，锚杆出现整体失稳；在锚杆的拉拔过程中，锚杆的变形规律是自由端处的变形量最大，自由端至锚固端锚杆的变形量逐渐变小，同时锚固剂变形的规律与锚杆的变形规律相同.关键词：FLAC－3D；有效围压；锚固力；摩擦角；数值模拟中图分类号：TD788.23文献标识码：A文章编号：0367－6234（2009）10－0129－05Numerical simulation on pull－tests of a cable by FLAC－3DJIANG Wen-wu1，2，XU Guo-yuan1，MA Chang-nian1（1.School of Resources＆Safety Engineering，Central South University，Changsha410083，China，wenwujiang@；2.Jiangxi University of Science and Technology，Ganzhou，341000，China）Abstract：To study major influential factors of cable anchorage force and the law of the whole instability in pull-tests for a cable，a numerical simulation model was established by3－D explicit finite difference method. It is proved that results of numerical simulation agree with field results of pull-tests.The pull-tests of the cable is a process in which the slow slippage turns abruptly into the whole distability.The anchorage effect of the ca-ble boosts up along with the augment of friction angle，adhesion stress and effective pressure.The shearing stress distribution regularities of anchorage agent change with the pull-tests process of the cable.At the initial stage of pull-tests，the shearing stress of anchorage agent in free end is maximum，while that in anchorage end is minimum.Cable grout yield points of free end are obtained firstly along with pull-tests，resulting in the slip-page of the cable and anchorage agent of free end.Shearing stress of anchorage agent takes on slowly homoge-neous distribution along with the further operation of pull-tests and reachs cable grout yield points from free end to anchorage end，so as that the whole cable becomes instable.The maximum deformation occurs in the free end during the pull-tests，the deformation amount of the cable is diminished gradually from free end to an-chorage end，and the deformation law of anchorage agent is uniform with that of a cable.Key words：FLAC-3D；effective pressure；anchorage force；friction angle；numerical simulation收稿日期：2006－07－14.基金项目：国家自然科学基金资助项目（59804007）；教育部“优秀青年教师资助计划”资助项目（EYTP－2134）.作者简介：江文武（1975—），男，博士研究生，讲师；徐国元（1965—），男，教授，博士生导师.岩体锚固技术在交通工程、矿业工程、隧道工程、水利工程等被广泛运用，适用于地质条件较差的岩体边坡稳定加固、隧洞的支护和混凝土结构的增强加固等.锚杆作为主要的加固手段，对岩体变形和强度起着重要的作用［1 4］.加锚岩体的数值模拟方法大都还是基于有限元法，但一般都过低估计锚固效果.然而FLAC －3D 即三维快速拉格朗日分析方法的出现，又为锚杆在岩体锚固机理提供了新的机遇.本文就锚固体的摩擦角、有效围压等对锚杆锚固性能的影响作了分析，对锚杆拉拔过程中锚杆锚固失效的特点进行了探索，并将现场试验与数值模拟计算进行了对比和分析.1数值模拟的平台安装了锚杆的锚固体包含岩体、锚杆和锚孔中的灌注砂浆几部分介质，对锚杆拔出过程进行数值仿真，应正确模拟锚固体的各个组成部分以及各部分之间的界面.数值模拟采用国际上著名的岩土工程分析软件FLAC －30［5］1.1轴向特性锚单元的轴向行为特性采用一维模型描述，轴向刚度K 为K =AE L.式中：A 为锚单元的横截面积，E 为弹性模量，L 为锚单元的长度.由轴向位移增量Δu t ，计算轴向力增量ΔF t：ΔF t =－K Δu t .其中，Δu t =（u ［b ］1－u ［a ］1）·t 1+（u ［b ］2－u ［a ］2）·t 2+（u ［b ］3－u ［a ］3）·t 3.式中：u ［b ］i ，u ［a ］i 为节点位移，i =1，2，3；［a ］，［b ］分别为节点a ，b ；t 1，t 2，t 3分别为锚单元轴线方向的方向余弦.锚单元可以指定其拉伸屈服强度F t 和压缩屈服强度F c ，单元轴力不能超过强度极限，如图1所示.拉伸屈服点F t拉力EA l轴向变形张拉压缩F c 压缩屈服点图1锚单元轴向力－变形特性曲线1.2界面模拟采用弹簧－滑块系统来描述锚杆（索）、砂浆及岩体之间的相互作用关系，如图2所示.该模型反映了锚杆（索）与岩体界面（浆体）之间的剪切特性.界面之间可能产生的最大剪切力取决于浆体的剪切刚度，单位厚度浆体的剪切刚度k g 为k g =2πGln （1+2t /D ）.式中：G 为浆体的剪切模量；t 为浆体环的厚度；D 为钢筋索的直径.浆体界面的剪应力τg 为τg =G （D /2+t ）Δuln （1+2t /D ）.式中：Δu 为浆体与岩体间的相对位移.浆体单位长度所能承受的最大剪应力为F max SL =c g +σm ˑtan （φg ）ˑp g .式中：c g 为浆体的粘结强度，φg 为浆体摩擦角，p g 为浆体与锚单元或岩体接触的实际周长，σm 为有效围压.界面材料特性假定为理想弹塑性，采用莫尔－库仑准则作为屈服准则.每个锚单元均允许沿轴向产生变形并发生屈服.若锚单元两节点与网格节点重合，锚单元与实体单元联结成整体，不产生相对位移；锚单元节点与网格节点不重合时，锚单元与实体元之间产生相对位移，其位移大小与界面模型参数有关.锚杆与岩体间的界面可以发生剪切屈服、产生滑动直至拉拔破坏.锚单元(钢筋)浆体环局部放大开挖面岩体钢筋的轴向刚度锚单元节点mmm滑片(浆体的剪切刚度)弹簧(浆体的剪切刚度)图2锚单元中的界面力学模型2计算模型的建立锚固体的计算模型取15m ˑ15m ˑ15m 的立方块，锚索布置在模型的中心，如图3、图4所示，模型由7200单元、7514个节点、20个结构体组成，锚杆的长度为5m ，数值模拟中，岩体采用8节点的六面体单元模拟，网格由锚索体向外呈·031·哈尔滨工业大学学报第41卷放射状逐渐由密变疏，如图3所示.在锚杆的自由端施加一个常速度v ，随着模型的计算，则锚杆的变形量等于计算的步速乘以常速度v ，同时在模拟锚杆拉拔的过程中，限制自由端所在的自由面Y 方向的位移，并且固定锚杆在X 、Z 方向速度为0.模型的本构关系采用莫尔库伦理想弹塑性理论.图3网格剖分图szzZXsxx 锚杆sxx X Y7.5mszz10m 5m 沿锚杆轴向施加固定的速度v限制Y 方向的位移图4锚杆拉拔数值模型示意图为了模拟锚杆拉拔过程中的影响因素，即影响锚杆锚固的效应的因素：1）模拟了在同样的外部条件下，唯有浆体的摩擦角（φg =00，100，200，300，400）不同的条件作用下沿着锚杆轴向、径向锚杆的应力与应变的分布规律以及锚杆的锚固力、浆体界面上的剪应力分布特征；2）模拟了在同样的外部条件下，唯有浆体有效围压（σm =0，2，4，6，8MPa ）不同的条件作用下沿着锚杆轴向、径向锚杆的应力与应变的分布规律以及锚杆的锚固力、浆体界面上的剪应力分布特征；同时还模拟了锚杆在拉拔过程中，锚索与岩体间的界面发生剪切屈服、产生滑动直至拉拔破坏具体过程.3数值模拟试验结果通过多种方案的数值模拟试验可知图5（a ）是现场试验得到一系列的力与位移之间的曲线，从图5（a ）中得知锚杆直径为15.2mm 的锚杆锚固力=17t /m.图5（b ）是根据现场的地质条件建模后计算得到的锚杆所受力与位移之间的曲线，图5（b ）中显示当锚杆自由端施加的力小于某一值时，力与位移基本成正比关系，当力达到一定值即锚固力时，力保持不变，而位移呈无限增大趋势，说明锚杆已经整体失稳，锚固作用失效，图6中显示锚固力=175kN /m =17.86t /m.图5（a ），图5（b ）对比来看可以得到在锚杆整体失稳、锚固作用失效之前力与位移之间的关系完全一致，并且得到的锚固力基本相等，而失稳后曲线的差异是由于FLAC 软件对锚杆单元力学特性的假设引起的.而工程中真正关心的是锚杆整体失效即力达到锚固力的过程，用FLAC 应用软件对锚杆拉拔试验进行模拟是切实可行的.30252015105015.2mm30.4mm 24mm 26mm51015202530位移/mm力/（t ·m -1)1.61.41.21.00.80.60.40.20.20.40.60.8 1.0×10-1(a)力与位移之间的曲线(b)锚杆所受力与位移之间的曲线图5现场试验结果与数值模拟结果的对比研究图6显示了当锚杆变形量各自为2.3、12.4、13.4mm 时锚杆的轴向力的分布情况以及锚杆和锚固剂的界面产生剪切屈服的情况，即锚杆和锚固体的界面产生滑动破坏情况.在锚杆位移为2.3mm 时，锚杆和锚固体的界面处刚刚产生剪切屈服点，即锚杆和锚固剂的界面刚刚产生了滑动或破坏；在锚杆位移为12.4mm 时，锚杆和锚固剂的界面进一步产生了滑动或破坏，但从图6中可以看出，锚杆和锚固剂界面产生滑动或破坏的过程比较缓慢；当位移为13.4mm 时，锚杆和锚固剂的界面全部剪切屈服，锚杆和锚固剂的界面迅速产生滑动，直到锚杆和锚固剂的界面整体产生失稳或滑动.所以锚杆在拉拔过程中，锚杆失效的过程是由量变到质变的过程，即缓慢滑动到突然整体失稳的过程.另一方面，在拉拔过程中，锚杆所受轴向力或应力的特点是锚杆的拉拔端受的力或应力最大，自由端或拉拔端至锚固端，锚杆所受的力或应力由大变到小直到接近为零.图7、图8显示了锚固体在拉拔过程中，锚固·131·第10期江文武，等：FLAC －3D 的锚杆拉拔数值模拟试验剂所受剪切应力的分布规律，在拉拔的初始阶段时，锚杆和锚固剂之间产生剪切应力，图7中显示(a)变形为2.3mm(b)变形为12.4mm(c)变形为13.4mm图6锚杆拉拔过程中锚杆轴向力分布与锚杆锚固剂屈服点图图7锚杆变形为2.3mm 时，锚固剂所受的剪切应力分布图沿着自由端至锚固端剪切应力呈现阶梯式分布，自由端处的剪切应力最大，沿着自由端直到固定端，剪切应力慢慢变小，接近为零；随着锚杆拉拔的进行，自由端处锚固剂与锚杆的界面屈服点首先达到，造成自由端处锚固剂与锚杆出现滑移现象；图8显示了锚杆变形量为12.4mm 时锚杆和锚固剂界面处的剪切应力分布规律，图中显示因为锚杆的变形量的增加，界面处剪切应力慢慢达到最大值，大概4m 的界面发生了剪切屈服；图9显示在锚杆拉拔结束阶段时，即锚杆变形量为13.4mm ，锚杆和锚固剂之间产生剪切应力分布规律，图9中显示剪切应力呈现均匀分布，剪切应力都达到了最大值，锚固剂与锚杆界面的屈服点均已达到，锚杆出现整体失稳.从锚固剂的剪切应力分析得知锚杆的失稳过程跟前面讲的失稳过程完全一致，锚杆失效过程也是慢慢滑动到突然整体失稳的过程.图8锚杆变形为12.4mm 时，锚固剂所受的剪切应力分布图图10是不同围岩应力下摩擦角与锚固力之间的曲线.锚杆的锚固力在摩擦角不同的情况下，而其它条件相同的情况时，当摩擦角为0ʎ时，尽管有效围压不同，但锚杆的锚固力都相等，它的值等于锚固剂和锚杆之间界面的粘结力，随着摩擦角的增大锚固力也随之增加，它与锚杆锚固作用机理理论是完全相符的.另一方面，锚杆的锚固力在摩擦角相同的情况下，而有效围压不同情况时，随着有效围压的增加，锚固力也按照一定的比例系数增大，即锚固力与围岩压力成正比关系；当摩擦角由0ʎ增加到40ʎ时，随着有效围压的增加，比例系数的值也跟着增加，也就是说，当有效围压改变量相等时，摩擦角值小时锚固力的增加量比摩擦角值大时锚固力的增加量要小，如当有效围压增加量为2MPa 时，摩擦角为10ʎ，锚固力的增加量约为30kN ，摩擦角为40ʎ，锚固力的增加量约为130kN ，增加量大了很多.从图11中可以得知，影响锚杆锚固作用效应除了锚杆周边岩体本身的地质条件外，主要有锚固剂和锚杆界面之间的摩擦角、粘结力以及有效围压等，从这些影响因素中，除了岩体的自身地质条件是天然的，不是人为能轻易改变的外，锚固剂与锚杆界面间的摩擦角、粘结力以及有效围压是人为控制的，根据实际的需要，选取适当匹配的锚固剂、锚杆的类型可以有效的控制界面间的摩擦角以及粘结力，采取一定的措施，改变锚杆锚固的·231·哈尔滨工业大学学报第41卷现场工艺，使锚固体在充填过程中密实，可以有效的提高锚固体所受的有效围压.图9锚杆变形为13.4mm 时，锚固剂所受的剪切应力分布图摩擦角/°锚固力/k N图10不同有效围压下锚杆锚固力与摩擦角的关系曲线摩擦角/°锚杆的变形量图11不同有效围压下锚杆变形量与摩擦角的关系曲线从摩擦角与锚杆变形量的曲线中，可以得到锚固力与摩擦角之间类似的结论，锚杆的变形量随着摩擦角的增加而增加，随着有效围压的增加而增加.锚杆的变形时的位移矢量如图12所示锚杆的自由端变形量最大，自由端至锚固端，锚杆的变形量也在慢慢变小直到接近为零.另一方面，在锚杆的拉拔过程中，锚固剂也发生变形，图13中在锚杆自由端处的锚固剂产生的变形量最大，从自由端至锚固端，锚固剂的变形量逐渐变小.图12锚固剂的位置锚固剂的变形量图13沿锚杆自由端至锚固端锚固剂的变形量曲线4结论1）数值模拟和现场试验的结果基本吻合，表明数值模拟锚杆拉拔过程是可行的；2）影响锚杆锚固效应的因素除了岩体自身的地质条件外，主要是锚固剂的摩擦角、粘结力和有效围压等，锚杆锚固效应随着锚固剂的摩擦角、粘结力和有效围压的增加而增强；3）锚固剂所受的剪切应力分布规律随锚杆拉拔过程而改变，4）在锚杆的拉拔过程中，锚杆的变形规律是自由端处的变形量最大，从自由端至锚固端锚杆的变形量逐渐变小，同时锚固剂也发生变形，变形的规律与锚杆的变形规律相同.参考文献：［1］杨强，任继承，张浩.岩石中锚杆拔出试验的数值模拟［J ］.水力学报，2002（12）：68－73.［2］杨松林，徐卫亚，黄启平.节理剪切过程中锚杆的变形分析［J ］.岩石力学与工程学报，2004，23（19）：3268－3273.［3］杜守继，职洪涛，翁慧俐，等.高速公路软岩隧道复合支护机理的FLAC 解析［J ］.中国公路学报，2003，16（2）：70－73，77.［4］丁秀丽，盛谦，韩军，等.预应力锚索锚固机理的数值模拟试验研究［J ］.岩石力学与工程学报，2002，21（7）：980－988.［5］Itasca Consulting Group Inc.FLAC3D （Version 2.0）us-ers manual ［R ］.USA ：Itasca Consulting Group Inc.，1997.（编辑张红）·331·第10期江文武，等：FLAC －3D 的锚杆拉拔数值模拟试验。

FLAC3D快速入门及简单实例

进阶书目
1. FLAC3D 在岩土工程中的应用【孙书伟，林杭，任连伟。中国水利水电出版社】
内容同书名一样，有一定的基础后，在处理实际问题时，可参考本书。 2. 各大论坛高手如云，卧虎藏龙，常去逛逛定会受益匪浅。
1
初识 FLAC3D
在这里说几句题外话，FLAC3D 的后处理功能虽然比较丰富，但却相对简单且不美观，若想做出漂亮的模型及结果，一般都要借助第三方软件。前处理可借用 ANSYS，网络上已经有高手写好了模型转换程序，可以将 ANSYS 模型导入 FLAC3D 进行计算，本文后面会有介绍；后处理软件选择比较多，常用的有 Excel、 Origin、Tecplot 等等，Origin 可以做出非常漂亮的曲线；借助网友编写的 fish 程序，还可以把 FLAC3D 计算结果导入 Tecplot 进行处理，关于这两个软件的使用，本文也会稍有介绍。
2
初识 FLAC3D
图 1.2 FLAC3D 基本菜单
基本菜单中，最常用的是 File 菜单。下面逐一介绍 File 菜单中各项的作用。 New：开始新工程，之前未保存的计算将被清除。命令：new 。 Call：导入命令流 txt 文件。命令：call filename.txt。 Model→Load：载入 dll 文件，用于载入用户开发的本构模型。 Restore：恢复保存的.sav 文件。命令：restore filename.sav。 Save：保存结算结果为.sav 文件。命令：save filename.sav。 Import Grid：导入.flac3d 模型文件，可以使 FLAC 自己生成的文件，也可是第三方软件(如 ANSYS)转成的符合 FLAC3D 模型数据格式的模型文件。命令： impgrid filename.flac3d。

flac命令及FISH小程序

1head1=nullp_z=zone_headloop while p_z#nullif ...new1=get_mem(2)mem(new1)=head1mem(new1+1)=p_zhead1=new1endifp_z=z_next(p_z)endloopend2这是fish中最常用的循环，即通过在内存中建立一链表，链表第一元素存储地址，第2个（或者更多）存储目标内容。

new = get mem(2)，即为在内存中建立一个两个元素的链表mem(new) = headmem(new+1) = p_gp为给这两个元素赋予内容3如果有whilestepping，此fish程序每时步执行一次，如果没有，fish程序只在读取到该程序时运行。

这里牵扯到链式结构和。

我根据我的理解解析一下。

def zs_topad = top_head；ad和top_head都是指针变量,此句将top_head的值；附给ad，以后ad代替top_head，指向get_mem创建的矩阵的第一个元素。

zftot = 0.0；顶点全部的力。

loop while ad # null；gp_pnt = mem(ad+1)；读取get_mem创建的矩阵第二个元素zf = gp_zfunbal(gp_pnt)；取出单元不平衡力zftot = zftot + zf；不平衡力叠加ad = mem(ad)；指针变量ad指向链表的下一个adendloopzs_top = zftot / 0.1414；力除以顶点处面积得顶点应力。

endballnew; p3 p4 是体对角线上的点，且一致。

def parmrad=10.0 ; radius of sphere.rad_size=5 ; radial zones.endparmgen zone pyramid p0 rad 0 0 p1 rad 0 rad p2 rad rad 0 p3 0 0 0 &p4 rad rad rad size rad_size rad_size rad_size group 1gen zone pyramid p0 0 rad 0 p1 rad rad 0 p2 0 rad rad p3 0 0 0 p4 rad rad rad & size rad_size rad_size rad_size group 2gen zone pyramid p0 0 0 rad p1 0 rad rad p2 rad 0 rad p3 0 0 0 p4 rad rad rad & size rad_size rad_size rad_size group 3def make_sphere; Loop over all GPs and remap their coordinates:p_gp=gp_headloop while p_gp#null; Get gp coordinate: P=(px,py,pz)px=gp_xpos(p_gp)py=gp_ypos(p_gp)pz=gp_zpos(p_gp)dist=sqrt(px*px+py*py+pz*pz)if dist>0 thenmaxp=max(px,max(py,pz))k=(maxp/rad)*(rad/dist)gp_xpos(p_gp)=k*pxgp_ypos(p_gp)=k*pygp_zpos(p_gp)=k*pzend_ifp_gp=gp_next(p_gp)end_loopendmake_spheregen zone ref ;z=0 planegen zone ref dip 90 ;y=0 planegen zone ref dip 90 dd 90 ;x=0 planeplot surfpl set back whpl bl grdataoflineforflac3d;; by tg0215@simwe 11/22/2006;; Function:Export gridpoint datas(disp,xdisp,ydisp,zdisp) and zone datas(sig1,sig2,sig3,sxx,syy,szz);; along a line defined by two point's coordinates.;; The exporting form is ten chart.The x-axes of the chart is the distance to the first point and;; the y-axes of the chart is gp data or zone data.;; Used for FLAC3D;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; Line Define by two points in spacedef Line_Definex1=in('Please input the x-coordinate of the first point:')y1=in('Please input the y-coordinate of the first point:')z1=in('please input the z-coordinate of the first point:')x2=in('Please input the x-coordinate of the second point:')y2=in('Please input the y-coordinate of the second point:')z2=in('Please input the z-coordinate of the second point:')numP=in('Please input the number of point on the line:')endLine_Define;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Initializationdef ini_Paraarray xx(numP)array yy(numP)array zz(numP)endini_Para;;;;;;;;;;;;;;; Calculate the coordinate of points on the linedef Line_CalculateLL=sqrt((x2-x1)^2+(y2-y1)^2+(z2-z1)^2)loop n(1,numP)if n=1 thenxx(n)=x1yy(n)=y1zz(n)=z1elseif n=nump thenxx(n)=x2yy(n)=y2zz(n)=z2elseif x2=x1 thenxx(n)=x1elsexx(n)=x1+LL*(n-1)/((numP-1)*sqrt(1+((y2-y1)/(x2-x1))^2+((z2-z1)/(x2-x1))^2)) endifif y2=y1 thenyy(n)=y1elseyy(n)=y1+LL*(n-1)/((numP-1)*sqrt(1+((x2-x1)/(y2-y1))^2+((z2-z1)/(y2-y1))^2)) endifif z2=z1 thenzz(n)=z1elsezz(n)=z1+LL*(n-1)/((numP-1)*sqrt(1+((y2-y1)/(z2-z1))^2+((x2-x1)/(z2-z1))^2)) endifendifendifendloopendLine_Calculate;; Calculate displacement magnitudedef get_gp_dispgp_disp=gp_xdisp(p_gp)*gp_xdisp(p_gp)gp_disp=gp_disp+gp_ydisp(p_gp)*gp_ydisp(p_gp)gp_disp=gp_disp+gp_zdisp(p_gp)*gp_zdisp(p_gp)gp_disp=sqrt(gp_disp)end;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; Write gp-related data to table,suck as disp,xdisp,ydisp,zdispdef write_gp_infoloop n(1,numP)p_gp=gp_near(xx(n),yy(n),zz(n))if p_gp # null thencaseof info_flagcase 1 ;dispxtable(1,n)=LL*(n-1)/float(nump-1)get_gp_dispytable(1,n)=gp_disp;aa=(xtable(1,n));bb=(ytable(1,n))case 2xtable(2,n)=LL*(n-1)/float(nump-1)ytable(2,n)=gp_xdisp(p_gp)case 3xtable(3,n)=LL*(n-1)/float(nump-1)ytable(3,n)=gp_ydisp(p_gp)case 4xtable(4,n)=LL*(n-1)/float(nump-1)ytable(4,n)=gp_zdisp(p_gp)endcaseendifcommand;print LL,n;print aa,bb;pauseendcommandendloopend;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Write zone-related data to table,such as sxx,syy,szz,sig1,sig2,sig3 def write_zone_infoloop n(1,numP)p_z=z_near(xx(n),yy(n),zz(n))if p_z # null thencaseof info_flagcase 5xtable(5,n)=LL*(n-1)/float(nump-1)ytable(5,n)=z_sig1(p_z)case 6xtable(6,n)=LL*(n-1)/float(nump-1)ytable(6,n)=z_sig2(p_z)case 7xtable(7,n)=LL*(n-1)/float(nump-1) ytable(7,n)=z_sig3(p_z)case 8xtable(8,n)=LL*(n-1)/float(nump-1) ytable(8,n)=z_sxx(p_z)case 9xtable(9,n)=LL*(n-1)/float(nump-1) ytable(9,n)=z_syy(p_z)case 10xtable(10,n)= LL*(n-1)/float(nump-1) ytable(10,n)=z_szz(p_z)endcaseendifendloopend;;;;;;;;;;;;;;;;;;;;;;;;;;;; Main function:draw curve chartdef draw_curveinfo_flag=1write_gp_infoinfo_flag=2write_gp_infoinfo_flag=3write_gp_infoinfo_flag=4write_gp_infoinfo_flag=5write_zone_infoinfo_flag=6write_zone_infoinfo_flag=7write_zone_infoinfo_flag=8write_zone_infoinfo_flag=9write_zone_infoinfo_flag=10write_zone_infocommandtable 1 name 'disp'table 2 name 'xdisp'table 3 name 'ydisp'table 4 name 'zdisp'table 5 name 'sig1'table 6 name 'sig2'table 7 name 'sig3'table 8 name 'sxx'table 9 name 'syy'table 10 name 'szz'set plot jpgpl create v-disppl set back whitepl set title text'Displacement along the line'pl add table 1 mark line xlabel 'distance along the line' ylabel 'displacement' set out v-disp.jpgpl hardcopy v-disppl create v-xdisppl set back whitepl set title text'X-Displacement along the line'pl add table 2 mark line xlabel 'distance along the line' ylabel 'x-displacement' set out v-xdisp.jpgpl hardcopy v-xdisppl create v-ydisppl set back whitepl set title text'Y-Displacement along the line'pl add table 3 mark line xlabel 'distance along the line' ylabel 'y-displacement' set out v-ydisp.jpgpl hardcopy v-ydisppl create v-zdisppl set back whitepl set title text'Z-Displacement along the line'pl add table 4 mark line xlabel 'distance along the line' ylabel 'z-displacement' set out v-zdisp.jpgpl hardcopy v-zdisppl create v-sig1pl set back whitepl set title text'Sig1 along the line'pl add table 5 mark line xlabel 'distance along the line' ylabel 'Sig1' set out v-sig1.jpgpl hardcopy v-sig1pl create v-sig2pl set back whitepl set title text'Sig2 along the line'pl add table 6 mark line xlabel 'distance along the line' ylabel 'Sig2' set out v-sig2.jpgpl hardcopy v-sig2pl create v-sig3pl set back whitepl set title text'Sig3 along the line'pl add table 7 mark line xlabel 'distance along the line' ylabel 'Sig3' set out v-sig3.jpgpl hardcopy v-sig3pl create v-sxxpl set back whitepl set title text'Sxx along the line'pl add table 8 mark line xlabel 'distance along the line' ylabel 'Sxx' set out v-sxx.jpgpl hardcopy v-sxxpl create v-syypl set back whitepl set title text'Syy along the line'pl add table 9 mark line xlabel 'distance along the line' ylabel 'Syy' set out v-syy.jpgpl hardcopy v-syypl create v-szzpl set back whitepl set title text'Szz along the line'pl add table 10 mark line xlabel 'distance along the line' ylabel 'Szz' set out v-szz.jpgpl hardcopy v-szzendcommandenddraw_curveenunew;set mech ratio 1e-6def _geoparmk=1.130 ;安全系数size1=40size2=10filecal0=string(size1)+'_'+string(size2)+'_cal'+'.sav'fosfile0=string(size1)+'_'+string(size2)+'_fos'+'.sav'bf=45 ; 改变此处坡角af=bf*pi/180 ;h1=20;x00=30.0x01=105.0x02=x00x03=x00+h1/tan(af)x04=x01x05=x03x06=105.0x07=x06;y00=0.0y01=y00y02=1.0y03=y00y04=y02y05=y02y06=y00y07=y02;z00=20.0z01=z00z02=z00z03=40.0z04=z00z05=z03z06=z03z07=z03end_geoparmgen zone brick p0 0 0 0 p1 30 0 0 p2 0 1 0 p3 0 0 20 size 12 1 8gen zone brick p0 30 0 0 p1 105 0 0 p2 30 1 0 p3 30 0 20 size size1 1 8gen zone brick p0 x00 y00 z00 p1 x01 y01 z01 p2 x02 y02 z02 p3 x03 y03 z03 &p4 x04 y04 z04 p5 x05 y05 z05 p6 x06 y06 z06 p7 x07 y07 z07 size size1 1 size2def modelMaterial; 弹性模量(Pa)flag=0E1=100e6; 泊松比poi1=0.35;; 粘结力(Pa)coh1=42e3/k;; 容重(N/m3)weight1=25e3; 膨胀角(度)dila1=0.0;; 内摩擦角(度)fri1=(atan((tan(20.0*pi/180))/k))*180/pibeta1=(sin(20.0*pi/180))/(1-2*poi1);afa1=E1*poi1poi1=0.5*(1-(sin(fri1*pi/180))/beta1)E1=afa1/poi1;抗拉强度ten1=0.01e6grav0=-9.80 ; 重力加速度(N/kg)dens1=-weight1/grav0K1=E1/(3*(1-2*poi1))G1=E1/(2*(1+poi1))endmodelMaterialmodel mohr ;采用莫尔库仑模型pro bulk K1 she G1 dens dens1 coh coh1 &friction fri1 dilation 0.fix x range x -0.1 0.1fix x range x 104.9 105.1fix yfix x y z range z -0.1 0.1set grav 0 0 grav0solvesave 1_130.savfish求出偶数和奇数def numberarray even(21),odd(21);even为偶数数组，odd为奇数数组m1=0m2=0loop n(0,41)a=n/2.0b=int(n/2)if a = b thenm2=m2+1even(m2)=nelsem1=m1+1odd(m1)=nendifendloopendnumberdef wriloop k(1,21)hh=out(even(k))endloopcommandpauseendcommandloop k(1,21)hh=out(odd(k))endloopendwrifish作图FISH作图有个方法仅供参考：因为应力都是与zone有关的，所以，可以编程求出该函数在每个zone中的值，然后用flac 重新运行出原网格，可以将你所关注的这个函数量作为flac中已经有的量（如szz)的初始条件来绘出，随便举了个例子set log onset logfile hua.loggen zone brick size 3 4 4;重新运行出的原网格mo elasprop bul=1e10 shear 6e9fix x range x -0.1 0.1fix x range x -2.9 3.1fix x y z range z -0.1 0.1fix y range y -0.1 0.1fix y range y 3.9 4.1;下面两个fish第一个读出原函数的zone值,第二个作为szz赋给zonedef fuzhiarray var(48)p_z=zone_headloop while p_z # nulln=z_id(p_z)var(n)=-10000*n;这var是已经读取的数组，这里随便给一个规律p_z=z_next(p_z)endloopendfuzhidef jialiloop m(1,48)azz=var(m)commandini szz azz range id m;赋值endcommandendloopendjialiini dens 2500set gra 0 0 -10plot contour szz outline onset log offget_plast(塑性区大小显示程序)def get_plastshearnow = 1tensionnow = 2shearpast = 4tensionpast = 8v_shear_now = 0v_tension_now = 0v_shear_past = 0v_tension_past = 0p_z = zone_headloop while p_z # nullif and(z_state(p_z,0),shearnow) = shearnow thenv_shear_now = v_shear_now + z_volume(p_z)endifif and(z_state(p_z,0),tensionnow) = tensionnow then v_tension_now = v_tension_now + z_volume(p_z)endifif and(z_state(p_z,0),shearpast) = shearpast thenv_shear_past = v_shear_past + z_volume(p_z)endifif and(z_state(p_z,0),tensionpast) = tensionpast then v_tension_past = v_tension_past + z_volume(p_z) endifp_z = z_next(p_z)endloopii = out('shear_now : ' + string(v_shear_now))ii = out('tension_now : ' + string(v_tension_now))ii = out('shear_past : ' + string(v_shear_past))ii = out('tension_past : ' + string(v_tension_past))endget_plastgroup遍列节点[1]gen zone brick size 6 8 8model mohrprop bulk = 1e8 shear = 0.3e8 fric = 35prop coh = 1e10 tens = 1e10set grav 0, 0, -9.81ini dens = 1000fix x range x -0.1 0.1fix x range x 5.9 6.1fix y range y -0.1 0.1fix y range y 7.9 8.1fix x y z range z -0.1 0.1group hua range x 2 4 y 3 5 z 3 5def loadgp_gp = gp_headload=0.0loop while p_gp#nullx=gp_xpos(p_gp)y=gp_ypos(p_gp)z=gp_zpos(p_gp)p_z=z_near(x,y,z)if z_group(p_z) = 'hua' thenload=load+gp_zfunbal(p_gp)endifp_gp = gp_next(p_gp)endlooploadg=loadendhist loadgplot his 1solveinrange遍列节点[1]gen zone brick size 6 8 8model mohrprop bulk = 1e8 shear = 0.3e8 fric = 35 prop coh = 1e10 tens = 1e10set grav 0, 0, -9.81ini dens = 1000fix x range x -0.1 0.1fix x range x 5.9 6.1fix y range y -0.1 0.1fix y range y 7.9 8.1fix x y z range z -0.1 0.1range name hua x 2 4 y 3 5 z 3 5def loadgp_gp = gp_headload=0.0loop while p_gp#nullif inrange('hua',p_gp)=1 thenload=load+gp_zfunbal(p_gp)endifp_gp = gp_next(p_gp)endlooploadg=loadendhist loadgplot his 1solvemov11; avi.dat; Create an AVI movie file (steptest.avi) using the step option. ; Snap a picture every 30 steps while solving.gen zone radcylinder size 25 1 25 25gen zone reflect normal -1 0 0gen zone reflect normal 0 0 -1gen merge 1e-5model mohrprop dens=1000 bu=1e8 sh=7e7 fric 25 coh 3.5e4 tens 1e10fix yfix x range x -24.9 -25.1fix x range x 24.9 25.1fix z range z -24.9 -25.1set grav 10ini szz -1e6 grad 0 0 1e6ini sxx -0.5e6 grad 0 0.5e6 0ini syy -0.5e6 grad 0 0.5e6 0plot create testviewplot current testviewplot add contour syy average outline onset movie avi step 1 file glj1.avimovie startsolve ratio 1e-2movie finishthn;================;导入网格数据以生成水面impgrid ww;为建立水面设置“虚拟界面”而分组为水上部分和水下部分group water_below range group 1 any group 2 anygroup water_above range group water_below notgen separate water_belowinterface 1 wrap water_below water_aboveset grav 0 -9.81 0water den 1000;==================================def CnumNum=0;p_i=i_headp_in=i_node_head(p_i)loop while p_in # nullNum=Num+1;p_in=in_next(p_in)endloopendCnum;++++++++++++++++++++++++++++++==def ini_vararray var(Num,3)i=0p_i=i_headp_in=i_node_head(p_i)loop while p_in # nulli=i+1var(i,1)=in_pos(p_in,1)var(i,2)=in_pos(p_in,2)var(i,3)=in_pos(p_in,3)p_in=in_next(p_in)endloopendini_var;删除单元组water_below和组water_above，为导入实际的材料分组做准备model null range group water_belowmodel null range group water_above;删除虚拟界面单元int 1 dele;===========================;建立实际材料分组部分;导入网格数据，重新分组impgrid wwgroup soil range group 2 any group 4 anygroup rock range group soil notattach face;生成水面def ini_ppminx=var(1,1)minz=var(1,3)maxx=var(1,1)maxz=var(1,3)loop t (2,Num)if var(t,1)<minxminx=var(t,1)endifif var(t,3)<minyminz=var(t,3)endifif var(t,1)>maxxmaxx=var(t,1)endifif var(t,3)>maxzmaxz=var(t,3)endifendlooppnt = gp_headloop while pnt # nullgpnum=gp_id(pnt)x=gp_xpos(pnt)y=gp_ypos(pnt)z=gp_zpos(pnt)dist_gp=x^2+z^2dist1=var(1,1)^2+var(1,3)^2if x<=maxx thenif x>=minx thenif z>=minz thenif z<maxz thenmak=0loop n (1,Num)disn=var(n,1)^2+var(n,3)^2if y<var(n,2) thenif disp_gp-disn<dist1 then mak=nendifendifendloopif mak>0 thenif y<var(mak,2) thengp_pp(pnt)=(var(mak,2)-y)*1.e4endifendifendifendifendifendifpnt = gp_next(pnt)end_loopendini_pp;建立真实的“界面”;gen separate rock;interface 2 sma 1.0 wrap rock soil;施加直边约束fix x y z rang y -0.1 0.1fix x z range x 0 600 z 0fix x z range x 0 z 0 700fix x z range x 0 370 z 700fix x z range x 600 z 0 220;施加斜边约束fix x z range plane normal 48 0 23 origin 600 0 220 fix x z range id 685fix x z range id 671fix x z range id 670fix x z range id 669fix x z range id 24fix x z range id 668fix x z range id 679fix x z range id 680fix x z range id 684fix x z range id 682fix x z range id 667fix x z range id 25fix x z range id 666fix x z range id 676fix x z range id 674fix x z range id 659fix x z range id 657fix x z range id 688fix x z range id 237 fix x z range id 235 fix x z range id 233 fix x z range id 3218 fix x z range id 3216 fix x z range id 3214 fix x z range id 3212 fix x z range id 3195 fix x z range id 3194 fix x z range id 3198 fix x z range id 3197 fix x z range id 3200 fix x z range id 3201 fix x z range id 3203 fix x z range id 3204 fix x z range id 3206 fix x z range id 3207 fix x z range id 3208 fix x z range id 2822 fix x z range id 2840 fix x z range id 2812 fix x z range id 2814 fix x z range id 2841 fix x z range id 80 fix x z range id 81 fix x z range id 82 fix x z range id 83 fix x z range id 84 fix x z range id 175 fix x z range id 178 fix x z range id 179 fix x z range id 2808 fix x z range id 2809 fix x z range id 2807 fix x z range id 2806 fix x z range id 2842 fix x z range id 2803 fix x z range id 2804 fix x z range id 2835 fix x z range id 2834 fix x z range id 2833 fix x z range id 2832 fix x z range id 2831fix x z range id 2827fix x z range id 2826fix x z range id 2844fix x z range id 2824fix x z range id 2823fix x z range id 2845fix x z range id 2846plot set cent 300 125.6 350 dist 2061 rotation 130 290 0 ang 22.5 mag 1.25 plo shosave th.savtxnewgen zone brick size 3 3 3 &p0 (0,0,0) p1 (3,0,0) p2 (0,3,0) p3 (0,0,1.5) &p4 (3,3,0) p5 (0,3,1.5) p6 (3,0,4.5) p7 (3,3,4.5)group Base; Create Top - 1 unit high for initial spacinggen zone brick size 3 3 3 &p0 (0,0,2.5) p1 (3,0,5.5) p2 (0,3,2.5) p3 (0,0,7) &p4 (3,3,5.5) p5 (0,3,7) p6 (3,0,7) p7 (3,3,7)group Top range group Base not;; Create interface elements on the top surface of the baseinterface 1 face range plane norm (-1,0,1) origin (1.5,1.5,3) dist 0.1ini z add -1.0 range group Topwater dens 1000set grav 0 0 -10def water_tablep_i=i_headp_ie=i_elem_head(p_i) ;界面1单元指针赋予p_ieloop while p_ie # null;返回界面单元的三个顶点的地址p_gp1=ie_vert(p_ie,1)p_gp2=ie_vert(p_ie,2)p_gp3=ie_vert(p_ie,3);以这三个界面单元顶点为水面的三个顶点，生成水面;========interface的节点和网格节点的地址是分别存放的x1=in_pos(p_gp1,1)y1=in_pos(p_gp1,2)z1=in_pos(p_gp1,3)x2=in_pos(p_gp2,1)y2=in_pos(p_gp2,2)z2=in_pos(p_gp2,3)x3=in_pos(p_gp3,1)y3=in_pos(p_gp3,2)z3=in_pos(p_gp3,3)commandwater table face x1,y1,z1 x2,y2,z2 x3,y3,z3 endcommandp_ie=ie_next(p_ie)endloopendmodel null range group Basemodel null range group Topgen zone brick size 3 3 3 &p0 (0,0,0) p1 (3,0,0) p2 (0,3,0) p3 (0,0,1.5) &p4 (3,3,0) p5 (0,3,1.5) p6 (3,0,4.5) p7 (3,3,4.5) group 1; Create Top - 1 unit high for initial spacinggen zone brick size 3 3 3 &p0 (0,0,2.5) p1 (3,0,5.5) p2 (0,3,2.5) p3 (0,0,7) &p4 (3,3,5.5) p5 (0,3,7) p6 (3,0,7) p7 (3,3,7)group 2 range group 1 not group Base not group Top not ini z add -1.0 range group 2water_tableint 1 delplo water安全度安全度问题1.计算之前 config zextra 12.根据计算结果，计算zone的点安全度f，赋给 zextra 1z_extra(p_z,1) = f3.绘图 pl con zextra 1p_z = zone_headloop while p_z # null计算安全度fz_extra(p_z,1) = fp_z = z_next(p_z)endloop测一下你的计算机的速度newdef raterate = (t1 - t0) / 6250.0enddef timet0 = clock / 100.0enddef time1t1 = clock / 100.0endgen zone brick size 5,5,5model elasprop bulk 1e8 shear 1e7fix y range y -0.1 0.1apply syy -1.0e6 range y 4.9,5.1time0step 50000time1print ratereturn查看各个时段不平衡力的具体数值查看各个时段不平衡力的具体数值采用his来记录你的计算，包括位移应力等命令his unbalhis gp(zone) zdis range (0 0 0) 或者id=？导出数据命令his write n vs m begin 时步 end 时步 file filename.his n表示你纪录的id m表示时步要导出不平衡力的具体数值his unbalstep 100000 or solvehis write 1 vs step begin 1 end 1000 file 123.his使用上述命令就可以查看各个时步下的不平衡力的具体数值得到主应力def getsigp_z=find_zone(3113)|s1=z_sig1(p_z)Ss2=z_sig2(p_z)s3=z_sig3(p_z)endgetsig地面曲面方程newgen zone brick p0 0,167.0072,0 p1 0,0,0 p2 242.6586,167.0072,0 p3 0,167.0072,32.0 size 56,81,11;;定义地面曲面方程DEF eqn537ARRAY c(45),cx(12),cy(12),v(70)c(1)=30.33465961749452c(2)=0.3768529280014758c(3)=-1.411874011262434c(4)=-0.4466983009380354c(5)=0.02498771591159658c(6)=-0.6225933592839830c(7)=-0.5572409253529970c(8)=0.2594590432443336c(9)=-0.2154254299775080c(10)=0.05611608056180144c(11)=-0.3573754488804990c(12)=0.02174566734976303c(13)=0.4051641051840550c(14)=0.1744265369646383c(15)=0.1147014441688443c(16)=-0.2088509864184064c(17)=0.1493365036211442c(18)=0.4311419368791974c(19)=0.008961362272148653c(20)=-0.07774920998867369c(21)=0.1811142002504400c(22)=-0.1519527001982752c(23)=0.07915472779033880c(24)=0.2559732910968667c(25)=0.1026065146630541c(26)=-0.2819625969135452c(27)=-0.09027129904676756c(28)=0.1130183438471391c(29)=-0.02169639756233838c(30)=0.1907291399047935c(31)=0.1297366637584324c(32)=0.08806739336824376c(33)=-0.2260336878126227c(34)=-0.2115733057775976c(35)=0.1148217058529207c(36)=-0.06508949301229160c(37)=-0.04511427839772801c(38)=0.1632228872577737c(39)=0.1931691979495782c(40)=0.08562312402296863c(41)=0.01872590692817278c(42)=-0.04198182565657024c(43)=0.04502760175694480c(44)=0.1160322072269608c(45)=-0.07952377661820480order%=44x=(x-(0.000000000000000))/(77.24063134751799) y=(y-(0.000000000000000))/(53.16004282387357) IF x<0.0 THENx=0.0END_IFIF x>3.14159265358979323846 THENx=3.14159265358979323846END_IFIF y<0.0 THENy=0.0END_IFIF y>3.14159265358979323846 THENy=3.14159265358979323846END_IFnc%=8cx(1)=1.0cy(1)=1.0LOOP j% (1,nc%)cx(j%+1)=COS(j%*x)cy(j%+1)=COS(j%*y)END_LOOPiv%=1LOOP j% (1,nc%+1)m%%=j%+1LOOP m% (1,j%)m%%=m%%-1v(iv%)=cx(m%%)*cy(j%-m%%+1)iv%=iv%+1END_LOOPEND_LOOPz=0.0LOOP j% (1,order%+1)z = z + c(j%)*v(j%)END_LOOPEND;;地面曲面方程定义完毕;;找出节点坐标后只要给出X,Y，然后调用eqn537就可以得到X,Y对应的Z;;定义第1层底面曲面方程DEF eqn566ARRAY c1(36),cx1(12),cy1(12),v1(70)c1(1)=28.62060766890477c1(2)=0.2866834946597056c1(3)=-0.8353575191107516c1(4)=-0.04915431878961194c1(5)=0.4039755560498526c1(6)=-0.3362136646742352c1(7)=-0.2665446036102640c1(8)=0.01228084747342192c1(9)=-0.2228703607987686c1(10)=-0.1056397207500379c1(11)=-0.01267143831082441c1(12)=-0.1996293677009942c1(13)=-0.06064874974162015c1(14)=-0.2080253561904090c1(15)=0.03108780356779321c1(16)=0.05965985858539676c1(17)=0.04285379002151423c1(18)=0.3260634131778066c1(19)=0.03617992023477324c1(20)=0.02800173263537663c1(21)=0.05300528967296362c1(22)=-0.07147038182113627c1(23)=0.09423171997760518c1(24)=-0.008003820958046695c1(25)=0.03654771159807278c1(26)=0.04195809554079616c1(27)=0.01776731215338432c1(28)=0.07995229437423655c1(29)=0.01184899721780239c1(30)=0.1914904626951018c1(31)=-0.03834379592630702c1(32)=-0.04935320463957614c1(33)=-0.04246583589991888c1(34)=-0.07508761240626857c1(35)=0.08511394861034457c1(36)=0.03158468483322942order%=35x=(x-(0.000000000000000))/(77.24063134751799)y=(y-(0.000000000000000))/(53.16004282387357)IF x<0.0 THENx=0.0END_IFIF x>3.14159265358979323846 THENx=3.14159265358979323846END_IFIF y<0.0 THENy=0.0END_IFIF y>3.14159265358979323846 THENy=3.14159265358979323846END_IFnc%=7cx1(1)=1.0cy1(1)=1.0LOOP j% (1,nc%)cx1(j%+1)=COS(j%*x)cy1(j%+1)=COS(j%*y)END_LOOPiv%=1LOOP j% (1,nc%+1)m%%=j%+1LOOP m% (1,j%)m%%=m%%-1v1(iv%)=cx1(m%%)*cy1(j%-m%%+1)iv%=iv%+1END_LOOPEND_LOOPz=0.0LOOP j% (1,order%+1)z = z + c1(j%)*v1(j%)END_LOOPEND;;地面曲面方程定义完毕;;找出节点坐标后只要给出X,Y，然后调用eqn536就可以得到X,Y对应的Z;;定义第2层底面曲面方程DEF eqn607ARRAY c2(45),px(12),py(12),v2(70) c2(1)=26.40374776302957c2(2)=-1.292234879166614c2(3)=-0.9987652673592940c2(4)=0.08847151489382983c2(5)=-6.524125629419486c2(6)=0.4002923561445450c2(7)=0.3031713620939287c2(8)=0.9293011989675050c2(9)=1.028374487465087c2(10)=0.7894526930446397c2(11)=-0.2710341955836033c2(12)=-1.833585849910791c2(13)=-0.6799558584481914c2(14)=1.202068992799922c2(15)=-0.6569131952237344c2(16)=1.079058946133437c2(17)=1.286040305186933c2(18)=0.3587835487607498c2(19)=-0.3055218324008147c2(20)=1.320108378629063c2(21)=1.360836777714370c2(22)=-1.086438562691549c2(23)=2.885522086015947c2(24)=-0.6574447894560952c2(25)=1.428461008364791c2(26)=0.8792263637632590c2(27)=0.3263370530986411c2(28)=-0.8544593089642435c2(29)=0.8177299018417186c2(30)=-0.8210160246224562c2(31)=0.2959394968085911c2(32)=-0.1081689772714739c2(33)=-1.546682781919654c2(34)=-0.8887755964444086c2(35)=-0.1826536331301314c2(36)=0.01080252224961026c2(37)=-0.2510335500991192c2(38)=1.172841958315430c2(39)=0.1775850816862299c2(40)=-2.437087642643127c2(42)=0.3321145625443073c2(43)=0.9669041422660922c2(44)=0.3462785144789184c2(45)=0.7677956032753429order%=44x=(x-(121.3293000000000))/(121.3293000000000)y=(y-(83.50360000000001))/(83.50360000000001)pcnt%=9px(1)=1.0py(1)=1.0px(2)=xpy(2)=yLOOP j% (3,pcnt%)ctr=-1.0+(j%-2)*(2.0/((pcnt%)-1.0))px(j%)=-1.0+2.0/(1.0+EXP(-(x-ctr)/0.12))py(j%)=-1.0+2.0/(1.0+EXP(-(y-ctr)/0.12))END_LOOPiv%=1LOOP j% (1,pcnt%)m%%=j%+1LOOP m% (1,j%)m%%=m%%-1v2(iv%)=px(m%%)*py(j%-m%%+1)iv%=iv%+1END_LOOPEND_LOOPz=0.0LOOP j% (1,order%+1)z = z + c2(j%)*v2(j%)END_LOOPEND;;地面曲面方程定义完毕;;找出节点坐标后只要给出X,Y，然后调用eqn607就可以得到X,Y对应的Z;;定义第3层底面曲面方程DEF eqn609ARRAY c3(66),px1(12),py1(12),v3(70)c3(1)=21.82428275547068c3(2)=2.221561754167763c3(3)=-1.581903359451976c3(4)=-2.316827507113655c3(5)=-2.703947398487063c3(6)=-1.246037280331843c3(8)=-2.339997159259053 c3(9)=-4.751446231799969 c3(10)=2.238112737137326 c3(11)=-2.307007727942070 c3(12)=2.704352656458496 c3(13)=6.218021412068429 c3(14)=6.082322362020687 c3(15)=-1.781025552909750 c3(16)=1.595374835381855 c3(17)=5.083409370236106 c3(18)=-12.05561353031727 c3(19)=-8.149539051933758 c3(20)=0.4629435785333617 c3(21)=3.854808615283558 c3(22)=-3.172213655564220 c3(23)=-7.493360842439101 c3(24)=18.03955724690294 c3(25)=15.92195138603906 c3(26)=5.276624705110507 c3(27)=4.834416374228492 c3(28)=-2.941631600267834 c3(29)=3.061011978625018 c3(30)=6.894707465545384 c3(31)=-19.55553399652331 c3(32)=-25.47011892985832 c3(33)=-8.508414420882716 c3(34)=-4.573107121770635 c3(35)=-2.159206741489122 c3(36)=3.630405781581091 c3(37)=-1.716300833677144 c3(38)=-1.710288020840512 c3(39)=13.39520120547237 c3(40)=26.62167851961098 c3(41)=9.133569050067894 c3(42)=2.465816025684572 c3(43)=5.760745683031747 c3(44)=3.387751985565470 c3(45)=-1.998506978790600 c3(46)=-0.1436671749659344 c3(47)=2.673753603392002 c3(48)=-3.414672834039968 c3(49)=-15.30717986499453 c3(50)=-5.830362201304671。

FLAC3D第一章

CHINA UNIVERSITY OF MINING AND TECHNOLOGY
FLAC/FLAC3D软件简介软件简介
基本特征
连续介质非线性，大应变模拟连续介质非线性，显式解题方案，为不稳定物理过程提供稳定解显式解题方案，界面或滑动面用来模拟可产生滑动或分离的离散面，从而模拟断层，节理或摩擦边界散面，从而模拟断层，内置材料模型：内置材料模型：
CHINA UNIVERSITY OF MINING AND TECHNOLOGY
岩土软件应用选择
总体而言，三个系列的程序在解决一些基本问题时具有一定的重迭区间，因此，这三个系列的程序或许都适用 FLAC/FLAC3D更适合于解决连续体的多介质多场耦合问题；UDEC/3DEC在解决结构面控制型问题时有独到的优势；PFC系列则特别适合脆性材料破裂发展和散粒体流动变形问题 Itasca软件是岩（土）工程问题的专用高级计算程序，不论使用哪种程序，研究者的专业知识是基本要求
三维网格图
坝体变形分布图
水平位移分布图
坝体应力分布图
CHINA UNIVERSITY OF MINING AND TECHNOLOGY
采用LS-DYNA模拟板材跌落水中激起浪花
LS-DYNA分析弹丸侵砌混凝土面板
基于强调折减法的金沙江某边坡稳定性分析（ABAQUS）
基于UDEC的金沙江某边坡稳定性分析
CHINA UNIVERSITY OF MINING AND TECHNOLOGY
Itasca软件的开发和应用历史软件的开发和应用历史
1984 – 应DNA and Corps. of Engineers开发UDEC以研究节理切割块体破坏的运动行为 1985 – 开发了岩土体连续介质分析的 FLAC 软件，帮助矿山咨询 1988 – 应加拿大安省Sudbury市 Falconbridge Ltd.要求开发3DEC，应用于深部矿山开采中遇到的岩爆问题 1993 – 由Flac扩展到FLAC3D 1994 – PFC2D/3D 应如下企业要求开发：Codelco Chile, Anglo American, Shell, Komatsu, Mitsubishi Heavy Industries, Taisei, Kajima, Pacific Consultants 应用于矿山开采的粒子流分析 1998 – 为FLAC 、UDEC增加了GIIC 2001 – FLAC/Slope 推出，计算边坡稳定 2001 – 应11家矿山公司要求开发REBOP，用于崩落法开采设计 2002 –开发BLO-UP 模拟岩石破裂和爆破飞石控制问题（基于PFC与 FLAC的耦合）

关于锚杆本构及其数值模拟分析总结与思考

纯拉：
纯拉
何礼理
纯剪拉剪
PILE
耦合
(1)塑性之前（OA 段），杆体拉力随着拉应变的增加而线性增长； (2)屈服后，杆体轴力保持不变（AB 段）；(3)当杆体的应变大于拉断破坏应变时，锚杆拉断,杆体拉力变为(BC 段)。
6
纯剪：
（1）在杆体受到的剪力达到抗剪极限能力之前，其剪力随着剪切位移的增加而线性增加（OA 段）；（2）当杆体达到剪切极限受力状态时，杆体产生剪断破坏。
拉剪耦合：
纯剪：拉剪耦合：
栾恒杰
拉剪
曹艳伟
PILE
6
耦合
蒋宇静
二、关于 FLAC3D 中材料及结构单元二次开发的思考
（1）深部巷道二次开发锚杆、锚索、锚网索支护稳定性控制模型上述文献[1-9]对 FLAC3D 锚杆单元抗拉、抗剪、拉剪耦合方面二次开发及破断理论进行了详细阐述，并在模拟实际工程取得良好效果。支护与围岩形成的相互协调的承载共同体是巷道围岩控制的核心，但是目前文章[1-18]仅仅是对锚杆、锚索支护结构单元自身进行研究，一是缺乏对恒阻吸能、恒阻大变形锚杆单元结构进行二次开发；二是缺乏二次开发锚杆单元与原结构单元在静态围岩应力场对比、动态应力场演变规律量化分析、函数拟合，静态围岩位移场对比、动态位移场演变规律量化分析、函数拟合；三是缺乏二次开发锚网索单元耦合与原结构单元在静态围岩应力场对比、动态应力场演变规律量化分析、函数拟合，静态围岩位移场对比、动态位移场演变规律量化分析、函数拟合，从而形成动、静应力作用下深部巷道二次开发结构单元锚杆、锚索、锚网索支护稳定性控制模型。围岩应力场稳定是巷道稳定的基础，国内外研究学者对围岩应力场演变进行了深入研究，其研究现状如下： Antonio, B.[19]对深部隧道围岩应力场进行深入研究，通过围岩应力和位移变化揭示应力场演变规律；Mohammad, R.Z.[20]研究了深部隧道围岩内弹塑性分区，对塑性损伤区给出了应力和位移的解析解；Srisharan, S.[21]将深部煤矿巷道围岩简化为等效连续体，并使用离散元法对围岩

基于flac3d隧道开挖的关键命令流

收稿日期：2019-12-03作者简介：包昊（1995-），男，硕士研究生，主要研究方向为隧道可靠度通信作者：方超（1990-），男，工程师，硕士，主要研究方向为岩土不确定性以及地下工程设计基于FLAC3D 隧道开挖的关键命令流包昊1，周旭辉1，葛彬1，方超2（1.河海大学土木与交通学院，南京210024； 2.安徽省综合交通研究股份有限公司，合肥230000）摘要：利用FLAC3D 软件开展数值模拟分析隧道工程中遇到的问题，以上海软黏土盾构为例，建立隧道开挖的有限差分模型，使用壳单元生成衬砌，锚索单元生成锚杆，用该方法得到的隧道模型与现实隧道更加贴切.运算结果表明：用该隧道模型诱发的地表沉降能用Peck 公式进行较好的拟合，基本符合高斯分布，可以有效地模拟隧道开挖引起的地表变形规律.关键词：FLAC3D ；隧道开挖；弹塑性求解；应力释放方法中图分类号：TK 730.2；Q 357.5文献标识码：ACommand Flow Analysis of Tunnel Excavation Based on FLAC3DBAO Hao 1，ZHOU Xuhui 1，GE Bin 1，FANG Chao 2（1.College of Civil and Transportation Engineering ，Hohai University ，Nanjing 210024，China ；2.Anhui Comprehensive Transportation Research Institute Co.Ltd.，Hefei 230000，China ）Abstract ：The numerical simulation of FLAC3D was used to analyze the tunnel engineering problems.Taking Shanghai soft clay shield as an example ，the finite difference model of tunnel excavation was established.The shell element was used to generate the lining ，and the cable element was used to generate the anchor rod.The tunnel model obtained by this method is more suitable to the real tunnel.The calculation results show that the ground settlement induced by the tunnel model can be well fitted by the Peck formula ，which basically conforms to the Gaussiandistribution and can effectively simulate the ground deformation law caused by tunnel excavation.Key words ：FLAC3D ；excavation of tunnel ；elastic and plastic method ；method of stress releasing随着城市化进程的不断加快，地表空间已经不能够满足人们的需求，地下空间的探索已然成为主流.盾构隧道施工技术是人们对于地下空间探索的最主要的方法之一.随着工程项目的增多，各地地形的复杂不一，因此隧道的开挖也遇到各种各样的问题，如何较为准确的预测隧道开挖过程中风险，成为国内外学者们尤为关注的问题.数值模拟分析成为如今分析隧道工程相关问题的重要方式之一.随着科技的发展，众多商业软件被开发出来用于模拟实际工程问题，与其他软件相较而言，FLAC3D 软件在用于隧道工程的模拟时具有多方面的优势［1］.首先，隧道开挖时，不同场地的物理力学参数不同，导致不同场地的本构模型之间具有差异，FLAC3D 软件内嵌多种本构模型，可以针对不同的场地进行合理的选择；第二，与实体（group ）单元不同的是，该软件本身拥有较多的结构单元用于模拟现实中的衬砌、锚索、梁等，在方便使用者的同时提高了模拟的准确性；第三，该软件为满足更多使用者的要求，其内置的FISH 语言使得参数的赋值以及数据的提取更加人性化.众所周知，隧道开挖的问题已经得到了广泛的分析［2-8］，而FISH 语言则是建模中必不可少的部分.第38卷第2期河南科学2020年2月本文简要介绍了隧道开挖模拟中关键部分的程序命令流，建立有限差分隧道模型，进行一次隧道开挖对地表位移影响的模拟运算，将结果与Peck ［9］经验公式拟合对比，证明其有效性.1关键命令1.1生成初始地应力场在采矿工程或者岩土工程领域中，必然存在着初始地应力场，它对于土体变形分析的影响不容小觑.传统初始地应力场的生成采用弹性求解法，而后再改换成塑性求解，忽略了土体的实际性质.而采用弹塑性求解法与前述方法相比可产生屈服的区域，相较而言，该方法初始地应力场的生成比前者更为合理.用简单例子更简易地表达弹塑性求解法的过程：newgen zone brick size 999model mohr ；采用摩尔库伦模型prop young …pois …fric …coh 3e10；将凝聚力设置为较大值fix xyz …；固定边界ini dens …set grav …；设置参数solve prop …coh 13e3；重新设置凝聚力solve需要要注意的是此简单例子只为说明生成初始地应力场的过程和方法，将其更为简明地展示出来，由于模型较为简单，故在自重作用之下并没有产生屈服区域，实际情形需要由使用者自己建立适用的模型进行观察分析.1.2应力释放应力释放法实际上就是应力的反向施加.张传庆等［10］分析了应力释放在隧道工程中的相关问题；程红战等［11-12］在将应力释放系数设为0.1的基础上建立隧道模型，分析了土体弹性模量的相关距离和变异系数对地表变形的影响；方超等［13］将围岩密度、弹性模量、内摩擦角视为三维正态随机场，研究围岩的相关距离对可靠度的影响，其中应力释放系数为0.30.应力释放方法的原理［14］是当土体开挖以后，在开挖边缘的单元节点上会失去原有的支持力，进行第一步计算（step 1）.此计算是为了获取其不平衡力P 0，将这些不平衡力以某一比例（应力释放系数a ）反向施加在原有的节点之上，紧接着添加shell 单元进行最后求解.应力释放后不能进行一次求解计算，必须添加衬砌后两者同时求解，否则隧道先变形后添加衬砌，其隧道掌子面变形量与衬砌变形量不相等，从而脱离实际.应力释系数的确定与当地的水文地质、开挖施工方法等都有一定关系，需综合分析确定.具体命令流为：def rel；定义FISH 语言coef=1.0-a c_x=41c_z=35--288引用格式：包昊，周旭辉，葛彬，等.基于FLAC3D隧道开挖的关键命令流［J］.河南科学，2020，38（2）：287-291.d=6.2；输入隧道中心点坐标以及隧道直径drelax1=gp_headloop while relax1#nullxa=gp_xpos（relax1）ya=gp_ypos（relax1）za=gp_zpos（relax1）dis=sqrt（（c_x-xa）^2+（c_z-za）^2）if dis<（d/2.+0.01）thenif dis>（d/2.-0.01）then；原理为隧道开挖掌子面距离隧道中心的距离等于半径，用该方法［15］来确定需要进行应力释放的节点xpow=-gp_xfunbal（relax1）*coefypow=-gp_yfunbal（relax1）*coefzpow=-gp_zfunbal（relax1）*coefkid=gp_id（relax1）commandapply xforce@xpow range id@kidapply yforce@ypow range id@kidapply zforce@zpow range id@kidendcommand；反向施加一定比例的不平衡力endifendifrelax1=gp_next（relax1）endloopend1.3设置锚索在隧道开挖工程的数值模拟中用锚杆和锚索的支护，常常用锚杆对岩石岩土工程进行加固，它的作用是利用水泥沿着长度方向提供的抗剪切能力，以生成局部阻力，借此抵御裂缝的位移变形，但是在目前已有的论文中极少有关于该命令流的介绍.由于隧道纵向长度远大于横向长度，将其视为平面应变情况，所以建模纵向距离取值1m，在该范围内的锚杆数量有限，用精确坐标的方法［16］即可完成锚杆布置的数值模拟.由于锚杆在圆形隧道四周呈放射状布置，故FLAC中的单元都以矩形方块为主，而放射状布置相对于单元形状是难以确定坐标的，可以使用CAD绘制准确图形，从CAD绘图软件中精确读取每根锚杆的坐标位置，并以其中一根锚杆的命令流为例，使用精确坐标法布置隧道四周的锚杆命令流：def cab_insloop iidx（1，6）y=iidx-0.5commandsel cable id=1begin39.32y43.16end36.88y42.62nseg4…endcommandendloopend-289-第38卷第2期河南科学2020年2月用该方法建立锚杆需输入每个点的坐标，故不适合较多坐标点的输入，否则既繁琐又容易出错，需慎重选择.锚杆一般与衬砌连用，共同作用于隧道掌子面，使得隧道变形降为最小值，就如上面添加衬砌一样，锚杆的添加也不是一步完成的，需要在此之前进行应力释放，由于前面已经定义了FISH 语言，这里不再复述.利用应力释放程序、衬砌（shell ）单元、锚杆（cable ）等建立较为完善的隧道开挖模拟.2隧道有限差分建模以在上海地区的软黏土中开挖隧道为原型，模拟在均质土体中开挖隧道对邻近建筑物的影响.有限差分模型的几何形状如图1所示，土体的宽度为120m ，深度为50m ，较大的边界有利于减少计算变形的误差.隧道的衬砌用shell 单元来模拟，由于在FLAC3D 中shell 单元为弹性连续环，这与实际衬砌的组装不相符，故需要对衬砌刚度进行折减，折减系数取0.7.隧道的轴线埋深为15m ，如表1所示，隧道的外径为6.2m ，内径为5.5m ，衬砌厚度为0.35m .由于隧道的纵向尺寸远远大于其截面的尺寸，故同前所述，纵向取值1m 假设为平面应变情形.假设衬砌为混凝土材料，其弹性模量、泊松比和重度分别为34.5GPa 、0.2和25kN/m 3.表1物理力学参数Tab.1Physico-mechanical parameters介质土体衬砌材料弹塑性材料线弹性材料衬砌厚度/m0.3重度/（kN·m -3）1825内摩擦角/（°）9.5黏聚力/kPa14泊松比0.320.22弹性模量/MPa121.5×104在各种数值模拟中，已有大量的土体模型用于模拟软土的非线性应力应变的特性.然而获得准确的土体参数进行合理的预测较为困难，故本研究采用Mohr-Coulomb 模型，这也是目前在模拟土体模型中应用最为广泛的数值模型之一.表1给出了土体参数，这是上海软黏土的常用值［17］.采用应力释放法开挖隧道，典型的上海软黏土的应力释放率为25%~30%，结合有限差分模型以及上海隧道的实际情形，本文的应力释放率取值0.25.均质土的地表位移如图2所示，该曲线能用Peck 的经验公式较好地拟合，基本服从高斯分布，也说明该有限差分模型可以有效地模拟隧道开挖引起的地表变形规律.图1有限差分模型Fig.1Finite difference model1201535土体深度/m土体宽度/m图2地表沉降与曲线拟合Fig.2Surface settlement and curve fitting地表沉降Peak 公式0-5-10-15-20地表沉降/m m-40-30-20-10010203040距隧道中心距离/m--290引用格式：包昊，周旭辉，葛彬，等.基于FLAC3D隧道开挖的关键命令流［J］.河南科学，2020，38（2）：287-291.3结语1）FLAC3D在分析隧道工程和采矿工程问题时具有较大的优势，内含的结构单元和本构模型可更为简便与准确地进行数值建模.但在模型较大、单元数量较多时，其计算过程较长，计算速度会显得较慢.2）传统的初始地应力场的生成虽然简便，但与实际有一定差距.通过改变强度参数的弹塑性法建立初始地应力场可优化这一过程，可反映土体的塑性区，但其计算速度亦会随着模型的增大减慢.3）使用应力释放法进行隧道开挖，可以在一定计算条件下较好地还原实地情形，但是由于应力释放率的确定与各种因素有关，故需综合确定.4）用FLAC3D建立有限差分模型，所得隧道开挖对地表位移的沉降曲线符合Peck经验公式，具有一定的有效性.参考文献：［1］陈育民.FLAC及FLAC3D基础与工程实例［M］.北京：中国水利水电出版社，2008.［2］ABID A，LYAMIN A V，HUANG J S，et al.Undrained stability of a single circular tunnel in spatially variable soil subjected to surcharge loading［J］.Computers&Geotechnics，2017，84：16-27.［3］MOLLON G，PHOON K K，DIAS D，et al.Validation of a new2D failure mechanism for the stability analysis of a pressurized tunnel face in a spatially varying sand［J］.Journal of Engineering Mechanics，2011，137（1）：8-21.［4］HUANG H W，GONG W P，KHOSHNEVISAN S，et al.Simplified procedure for finite element analysis of the longitudinal performance of shield tunnels considering spatial soil variability in longitudinal direction［J］.Computers&Geotechnics，2015，64：132-145.［5］BERNAT S，CAMBOU B.Soil-structure interaction in Shield Tunnelling in Soft Soil［J］.Computers&Geotechnics，1998，22（3）：221-242.［6］ORESTE P A.Probabilistic design approach for tunnel supports［J］.Computers&Geotechnics，2005，32（7）：520-534.［7］李志华，华渊，周太全，等.盾构隧道开挖面稳定的可靠度研究［J］.岩土力学，2008（S1）：315-319.［8］宋建康，刘春原.复杂地质条件下隧道支护体时效可靠性分析［J］.河北工业大学学报，2018，47（4）：113-118.［9］PECK R B.Deep excavation and tunnelling in soft ground［C］//Proceedings of the7th International Conference on Soil Mechanics and Foundation Engineering.Mexico：Sociedad Mexicana de Mecanica de Suelos，A.C.，1969：225-290.［10］张传庆，冯夏庭，周辉，等.应力释放法在隧洞开挖模拟中若干问题的研究［J］.岩土力学，2008（5）：1174-1180.［11］程红战，陈健，李健斌，等.基于随机场理论的盾构隧道地表变形分析［J］.岩石力学与工程学报，2016，35（S2）：4256-4264.［12］李健斌，陈健，罗红星，等.基于随机场理论的双线盾构隧道地层变形分析［J］.岩石力学与工程报，2018，37（7）：1748-1765.［13］方超，薛亚东.围岩空间变异性对隧道结构可靠度的影响分析［J］.现代隧术，2014，51（5）：41-47.［14］HUANG H W，XIAO L，ZHANG D M，et al.Influence of spatial variability of soil Young’s modulus on tunnel convergence in soft soils［J］.Engineering Geology，2017，228：327-370.［15］王涛，韩煊，赵先宇，等.FLAC3D数值模拟方法及工程应用——深入剖析FLAC3D5.0［M］.北京：中国建筑工业出版社，2015.［16］刘永乐.FLAC~（3D）5.0巷道锚杆安装命令流［J］.煤炭技术，2018，37（11）：124-126.［17］ZHANG Z G，HUANG M S.Geotechnical influence on existing subway tunnels induced by multiline tunneling in Shanghai soft soil［J］.Computers&Geotechnics，2014，56（3）：121-132.（编辑孟兰琳）-291-。

基于FLAC 3D的煤巷锚杆支护参数的模拟分析

全区断裂可分两组，一组Ｆ和岗德尔逆断层，为近南北走
（泥岩）和９的标志层。煤层顶底板岩性：直接顶为为９下＃
煤岩互层，厚度为２３ｍ，节理发育，抗压强度极低。老．０
顶为粗砂岩，厚度为１４ｍ，灰黑色，厚层状、较坚硬。．９
范围的破坏区Ｈ。锚杆支护的作用就在于保持破坏区范ｊ围内岩层的稳定性。锚杆支护设计关系到巷道锚杆支护工
Ｌ：Ｋ＋Ｈ１＋Ｌ２式中 — — 锚杆长度，ｍ；
（）１
日 —一软弱岩层厚度或冒落拱高度，ｍ；
— 一
安全系数，一般取２；
代人相关数据，由式（）：ａ＝１３ｍ。３得．２
３数值模拟分析３１模型建立．
针对神华集团平沟煤矿的地质采矿条件，根据采矿工
ｌ
程问题特点，建立相应的数值分析模型，考虑采空区的影
响，模型水平方向上取６ｍ，这样既可以模拟已开采煤层０对地表的影响和采空区周边围岩的响应状态，又可以模拟工作面巷道掘进的影响。根据岩层情况，竖直方向上取５ｍ，巷道设计高Ｘ宽＝２５４３．ｍＸｍ。采用应力边界条件，模型的上边界施加均匀的垂直压应力，模型的左右边界限制水平向移动，模型下边界固定，上边界自由］，模型如
结构简单，裂隙较发育，煤层中部含一层厚约０２ｍ夹石．８
１煤巷地质概况
平沟煤矿位于桌子山背斜西翼之中北部，因受岗德尔逆断层逆冲影响，使地层上挠而成不对称平缓向斜，全区的构造形态为向西倾斜的单斜地层为主，倾角为６。一１。Ｏ。

洞山隧道锚杆支护作用的FLAC模拟

ＣＡＮＮＤＧＥＯＴＥＣＨＮＩＬＥＮＧＩＥＥＲＩＧＷＯＲＬ
… … … … … … … … … … … … … ． … … … … … … … … …
…
ＶＯＬ１１
．
Ｎ
…
…
…
…
…
…
…
…
…
…
…
…
…
…
…
ｏ１．
３模型建立及ＦＡＣ数值计算Ｉ
洞山隧道断面支护模型见图１。
图５模型影响范围内岩土体ｌ方向应力等值线图，
由隧道模型影响范围内岩土体的向应力等
值线图（４和ｌ向应力等值线图（５可以看图），图）出，向较大应力发生在拱脚的一定范围内，即向应力主要集中在拱脚到拱腰的一定范围内，主要其
和灵活性。
１ＦＡＬＣ程序简介
ＦＡＦｓＬｇｎａｎｌｓｆｏｔｕｍ）Ｌｃ（ａｔａｒ百ｎＡａｉｏｎｉｕ分ａｙｓＣｎ析软件是美国ＩＡＣＴＳＡ咨询集团公司开发研制的一
种以牛顿第二定律为基础的有限差分法程序，以可
等值线图（３图）可以看出，本模型能够充分反应隧道因开挖而发生位移的空间效应。向位移由轴线向两边逐渐减小，较大位移发生在距隧道拱腰较且近的一定高度的范围内；Ｙ向最大位移发生在拱顶
附近，由拱顶向两边逐渐减小，明本分析比较符且说合实际。

基于FLAC模拟的基坑锚杆(土钉)支护分析

坑坡顶１ｍ。
因此，型中土钉和锚杆采用ＦＡ模ＬＣ中的锚杆单元（ａｌｃｂｅ
关键词：杆（锚土钉）护，坑，ＬＣ数值模拟支基ＦＡ
中图分类号：４５７Ｕ５．１
文献标识码：Ａ地表以下２ｍ，二排锚杆设置于地表以下６ｍ。基坑面层采用第
坑顶
０引言
问题之一。锚杆（土钉）护技术在深基坑工程中得到了较为广支泛的应用。一般情况下采用土钉支护即可，但是对于基坑的水平
２
形。本文以某城市基坑工程复合锚杆（土钉）护的设计与应用支为实例，使用ＦＡＬＣ数值模拟软件对复合锚杆（土钉）支护进行了
初步分析，为基坑复合锚杆（土钉）支护技术的设计和施工提供有益的参考。
图１预应力锚杆＋土钉复合支护设计剖面
粉质粘土中砂卵石
５５．５５
２０２．０３２．Ｏ６
２０００
２２３３４０
１８２５５０
０３．３０３０２．５
具体计算时根据变形模量Ｅ。利用公式（），１将其转化为Ｇ和Ｋ：
Ｋ，Ｇ（）１
１工程概况及方案设计
１１工程概况．
２数值模拟分析
对于锚杆（土钉）结构，数值模拟方法可以从宏观整体趋势上
ＦＡＣ适工程拟建场地平面形状为矩形，寸为１０ｒ×２基坑开分析岩土体开挖效应的变化规律。而众多数值方法中，Ｌ尺０ｎ８ｍ，用于大多数工程力学问题。因此，次数值模拟分析主要运用本挖深度为１０ｍ。基坑重要性等级为一级。基坑分５步开挖，每步ＦＡＬＣ程序进行计算。开挖２ｍ。基坑南侧为一条２宽城市道路，５ｍ道路南侧为多层住宅小

使用FLAC自带FISH语言安装锚杆实例

使用FLAC自带FISH语言安装锚杆实例define dingmaoganaa = 1xx = -3.7m1 = 2m2 = 14m3 = 15m4 = 24m5 = 1loop while aa < 7commandsel cable id aa beg xx 0.5 0.8 end xx 0.5 2.2 nseg 14sel cable id aa beg xx 0.5 2.2 end xx 0.5 3.2 nseg 10sel cable id aa prop emod 2e10 ytension 310e3 xcarea 0.0004906 gr_coh 1 gr_k 1 gr_per 0.0785 range cid m1 m2;sel cable id aa prop emod 2e10 ytension 310e3 xcarea 0.0004906 gr_coh 10e5 gr_k 2e7 range cid m3 m4;sel cable id aa prop emod 2e10 ytension 310e3 xcarea 0.0004906 gr_coh 10e8 gr_k 2e10 range cid m5 m5;sel cable id aa pretension 60e3 range cid m5 m2;end_commandaa = aa +1xx = xx + 0.7m1 = m1 + 24m2 = m2 + 24m3 = m3 + 24m4 = m4 + 24m5 = m5 + 24end_loopenddingmaogan;安装右帮锚杆define youmaoganaa = 7zz = -2.7m1 = 146m2 = 158m3 = 159m4 = 168m5 = 145loop while aa < 10commandsel cable id aa beg 0.8 0.5 zz end 2.2 0.5 zz nseg 14sel cable id aa beg 2.2 0.5 zz end 3.2 0.5 zz nseg 10sel cable id aa prop emod 2e10 ytension 310e3 xcarea 0.0004906 gr_coh 1 gr_k 1 gr_per 0.0785 range cid m1 m2;sel cable id aa prop emod 2e10 ytension 310e3 xcarea 0.0004906 gr_coh 10e5 gr_k 2e7 range cid m3 m4;sel cable id aa prop emod 2e10 ytension 310e3 xcarea 0.0004906 gr_coh 10e8 gr_k 2e10 range cid m5 m5;sel cable id aa pretension 60e3 range cid m5 m2;end_commandaa = aa +1zz = zz + 0.9m1 = m1 + 24m2 = m2 + 24m3 = m3 + 24m4 = m4 + 24m5 = m5 + 24end_loopendyoumaogan;安装左帮锚杆define zuomaoganaa = 10zz = -2.7m1 = 218m2 = 230m3 = 231m4 = 240m5 = 217loop while aa < 13commandsel cable id aa beg -4.6 0.5 zz end -6.0 0.5 zz nseg 14sel cable id aa beg -6.0 0.5 zz end -7.0 0.5 zz nseg 10sel cable id aa prop emod 2e10 ytension 310e3 xcarea 0.0004906 gr_coh 1 gr_k 1 gr_per 0.0785 range cid m1 m2;sel cable id aa prop emod 2e10 ytension 310e3 xcarea 0.0004906 gr_coh 10e5 gr_k 2e7 range cid m3 m4;sel cable id aa prop emod 2e10 ytension 310e3 xcarea 0.0004906 gr_coh 10e8 gr_k 2e10 range cid m5 m5;sel cable id aa pretension 60e3 range cid m5 m2;end_commandaa = aa + 1zz = zz + 0.9m1 = m1 + 24m2 = m2 + 24m3 = m3 + 24m4 = m4 + 24m5 = m5 + 24end_loopendzuomaogan;巷道右侧斜锚杆安装由上而下sel cable id 13 beg 0.6 0.5 0.6 end 1.6 0.5 1.6 nseg 14sel cable id 13 beg 1.6 0.5 1.6 end 2.3 0.5 2.3 nseg 10sel cable id 13 prop emod 2e10 ytension 310e3 xcarea 0.0004906 gr_coh 1 gr_k 1 gr_per 0.0785 range cid 290 302sel cable id 13 prop emod 2e10 ytension 310e3 xcarea 0.0004906 gr_coh 10e5 gr_k 2e7 range cid 303 312sel cable id 13 prop emod 2e10 ytension 310e3 xcarea 0.0004906 gr_coh 10e8 gr_k 2e10 range cid 289 289sel cable id 13 pretension 60e3 range cid 289 302sel cable id 14 beg 0.8 0.5 0.0 end 2.10.5 0.5 nseg 14sel cable id 14 beg 2.1 0.5 0.5 end 3.0 0.5 0.8 nseg 10sel cable id 14 prop emod 2e10 ytension 310e3 xcarea 0.0004906 gr_coh 1 gr_k 1 gr_per 0.0785 range cid 314 326sel cable id 14 prop emod 2e10 ytension 310e3 xcarea 0.0004906 gr_coh 10e5 gr_k 2e7 range cid 327 336sel cable id 14 prop emod 2e10 ytension 310e3 xcarea 0.0004906 gr_coh 10e8 gr_k 2e10 range cid 313 313sel cable id 14 pretension 60e3 range cid 313 326;巷道左侧斜锚杆安装由上而下sel cable id 15 beg -4.4 0.5 0.6 end -5.4 0.5 1.6 nseg 14sel cable id 15 beg -5.4 0.5 1.6 end -6.1 0.5 2.3 nseg 10sel cable id 15 prop emod 2e10 ytension 310e3 xcarea 0.0004906 gr_coh 1 gr_k 1 gr_per 0.0785 range cid 338 3500.0004906 gr_coh 10e5 gr_k 2e7 range cid 351 360sel cable id 15 prop emod 2e10 ytension 310e3 xcarea 0.0004906 gr_coh 10e8 gr_k 2e10 range cid 337 337sel cable id 15 pretension 60e3 range cid 337 350sel cable id 16 beg -4.6 0.5 0.0 end -5.9 0.5 0.5 nseg 14sel cable id 16 beg -5.9 0.5 0.5 end -6.8 0.5 0.8 nseg 10sel cable id 16 prop emod 2e10 ytension 310e3 xcarea 0.0004906 gr_coh 1 gr_k 1 gr_per 0.0785 range cid 362 374sel cable id 16 prop emod 2e10 ytension 310e3 xcarea 0.0004906 gr_coh 10e5 gr_k 2e7 range cid 375 384sel cable id 16 prop emod 2e10 ytension 310e3 xcarea 0.0004906 gr_coh 10e8 gr_k 2e10 range cid 361 361sel cable id 16 pretension 20e3 range cid 361 374;安装锚索define maosuoaa = 17xx = -3.4m2 = 414m3 = 415m4 = 426m5 = 385loop while aa < 20commandsel cable id aa beg xx 0.5 0.8 end xx 0.5 6.8 nseg 30sel cable id aa beg xx 0.5 6.8 end xx 0.5 8.0 nseg 12sel cable id aa prop emod 2e10 ytension 310e3 xcarea 0.0004906 gr_coh 10e5 gr_k 2e7 range cid m3 m40.0004906 gr_coh 10e8 gr_k 2e10 range cid m5 m5sel cable id aa pretension 60e3 range cid m5 m2 end_commandaa = aa + 1xx = xx + 1.5m3 = m3 + 42m4 = m4 + 42m5 = m5 + 42end_loopendmaosuo。

FLAC用于煤巷锚杆支护设计的二次开发_高富强

巷道工程与支护FL AC 用于煤巷锚杆支护设计的二次开发高富强1,张作樵2,康建东1(1.中国矿业大学能源与安全工程学院,江苏徐州221008;2.淮北矿业集团涡北煤矿,安徽淮北235000)[摘　要]　应用V is ua l Basic 进行了FLA C 用于煤巷锚杆支护设计的二次开发,实现了模型建立、运算、分析的自动化、便捷化,并以工程实例验证了该程序的处理功能和适用性。

[关键词]　FLAC ;二次开发;锚杆支护[中图分类号]TD311.52 [文献标识码]A [文章编号]1006-6225(2006)04-0054-03The Second D evelop m ent of FLAC for Bolti n g Supporti n g Design in Coal Road wayGAO Fu -qiang ,ZHANG Zuo -q iao ,KANG Jian -dong(1.Energy ＆S af ety C ollege ,Ch i na Un i versity ofM ining ＆Technol ogy ,Xu z hou 221008,C h i na ;2.Guobei Colli ery ,Hu ai beiM i n i ng G roup ,H uai bei 235000,C h i na )Ab strac t :T he second deve lopment o f FLAC fo r bo lting s upporti ng desi gn i n co al road w ay is made w it h V isual Basic .I t makes the m ode l buil ding ,runn i ng and anal y sis au t om atic ,fa st and conven i en.t The dea li ng func tions and its app licabilit y o f this prog ramm e are proved by eng inee ri ng cases .K ey word s :FLAC ;the second develop m en t ;bolting s upporti ng[收稿日期]2006-03-10[作者简介]高富强(1981-),男,河南周口人,硕士,主要从事采矿工程数值模拟方面的研究。

基于FLAC的煤巷锚杆支护设计系统开发与应用

ａＭｏｍｅｄ，ｔｅｒｏｕａｉｎｍｅｈｄａｄｎｍｅｃｌｉｌｔｎｎ￣￣ｏｌｏｙｃｍｐｔｔｔｏｎｕｒａｓｍｕａｉｍｅｈｄａｅｕｅｏｄｓｇｏｋｂｌｎｎｔｅｌｏｉｏｔｏｌｓｄｔｅｉｒｃｔｇｉｈｎｏｉ
构产生革命性的变化，使工业自动化系统朝着智能化、数字化、信息化、网络化、分散化的方向迈进。孔庄煤矿选
系统主要功能为：实时监视各用电设备的主要电量 ① 参数和故障信息，对故障状态作出报警或保护；② 实时监
视主要运输设备运行状态，对故障状态作出声光报警或保护；③实时记录主要设备用电参数和设备故障信息，并随
ｓｓｅＴｅｓｓｅｈｈｕｃｉｎｏｄｌｇｎｒｔｇａｔｍａｉａｌｈｃｎｂｅｎ —ｓｔｅｈｉｉｎｅｉｏｋｙｔｍ．ｈｙｔｍａｔｅｆｎｔｆｍｏｅｅｅａｉｕｏｔｌｙｗｉｈｅａｌｓｉｓｏｎｃｉｔｃｎｃａｓｔｄｓｇｒｃｕｏｎ
ｂｌｎｉＬｘｅｉｎｌＡｎａｏｒｈｎｉｅａｐｉａｉｎｉｏｇｈａｉｇｎｓｏｔｇｗｔＦＡＣｅｐｄｅｔｉｈｙ．ｄｈａｃｍｐｅｅｓｐｌｔｎＴｎｃｕｄｇｉｇ．ｓｖｃｏｎ
ｓｓｅｓｄｎｙｔｍｂａｅｏＦＬＡＣ
ＧＡＯＦｕ—ｑａｇｉｎ
（．ｅｉｅｅｒｈＩｓｔｔｏｏｌＭｉｉｇｅｔｌｏｅｅｒｈＩｓｉｔ，Ｂｉｎ００３Ｃｉａ１ＢｉｎＲｓａｃｎｔｕｅｆａｎｎ，ＣｎｒａＲｓａｃｔｕｅｅｉｇ１０１ｈｎ；ｊｇｉＣａＣｌｎｔｊ

基本FLAC 模拟与边际分析的锚杆方案选择

《装备维修技术》2020年第4期— 409 —基本FLAC 模拟与边际分析的锚杆方案选择张颖超1曹佳2孙伟3（1.青岛黄海学院山东青岛 266247；2.青岛海尔特种电冰柜有限公司山东省青岛市 2665553. 青岛黄海学院山东青岛 266555）摘要：本文通过建立巷道模型，通过使用计算机仿真软件FLAC,模拟描述锚杆直径变化对巷道受到应力的变化，形成变量关系；运用边际理论分析，结合支护投入，确定最优方案。

关键词：FLAC；锚杆直径；巷道应力；边际理论；最优方案引言锚网索支护是目前煤矿巷道支护采用的主要方法,它从成巷道作业成本,作业速度,断面的可重复利用率,人员劳动强度,安全作业等方面有巨大优势[1],对于锚网索支护的选择应用，除了岩层等外界条件外，其本身的性能也很重要，如影响巷道支护的具体参数有锚杆直径，长度，间距等。

本文从锚杆直径角度入手，基于模拟软件数值，运用边际成本分析法,从杆径，围岩变量，支护成本三方面，对影响巷道围岩稳定性的锚杆直径进行评价筛选，来确定合理最优的方案。

1 锚网索支护初步支护参数的确定在现场采矿过程中，受重的载荷主要来自于岩层上部，也就是上部岩层的自重；根据采集到的某矿一段采煤区域，简化为以下模型：1)只考虑上层地质与岩重，岩石材料与重量为厚砂岩26m 、，厚粉砂岩3.0m 、厚砂岩10.8m ，弱层1.2m;2)开采区域模型为长方形，采集数据为长65m,宽65米，高55m ;3)掘进巷道为梯形，位于其中：A 、B 、C 、D 是顶板上随机取的四个点，用来作为观察模拟持续受到相同同的岩重压力，不同的杆径产生的相应的变形量。

运用FLAC 软件模拟巷道的各个点的应力，随机取顶板上的A,B,C,D 四点位移变化，因此我们可以得到在FLAC 模拟软件下顶板四个点的不同变形量按照顺序为：150.4mm，173.3mm ，173.3mm ，167mm ，这四个数值可以作为选择影响围岩稳定性其中因子锚杆直径的参考对比值。

flac命令流

flac命令流newgen zone brick size 1 1 2model mohrprop bulk 3e7 shear 1e7 coh 10e3 fri 15 ten 0fix z ran z 0fix x ran x 0fix x ran x 1fix y ran y 0fix y ran y 1ini dens 2000 ran z 0 2ini szz -50e3 grad 0 0 20e3 ran z 0 1ini syy -30e3 grad 0 0 10e3 ran z 0 1ini sxx -30e3 grad 0 0 10e3 ran z 0 1ini pp 30e3 grad 0 0 -10e3 ran z 0 2app nstress -10e3 ran z 2set grav 0 0 -10solve；深埋工程初始应力场的S-B法生成newgen zone brick p0 0 0 0 p1 60 0 0 p2 0 60 0 p3 0 0 90 & p4 60 60 0 p5 0 60 90 p6 60 0 150 p7 60 60 150 & size 6 6 10model elasprop bulk 10e10 shear 10e10ini dens 2500apply sxx -1e9 grad 0 0 1.1111111e7 ran x -0.1 0.1 apply sxx -1e9 grad 0 0 7.6666666e6 ran x 59.9 60.1 apply syy -1e9 grad 0 0 8.3333333e6 ran y -0.1 0.1apply syy -1e9 grad 0 0 8.3333333e6 ran y 59.9 60.1apply szz -1e8 grad 0 0 8.3333333e5 ran z 0 120set grav 0 0 -10step 30000ini xdisp 0 ydisp 0 zdisp 0ini xvel 0 yvel 0 zvel 0plot cont szz常规方法生成newgen zone brick p0 0 0 0 p1 60 0 0 p2 0 60 0 p3 0 0 90 &p4 60 60 0 p5 0 60 90 p6 60 0 150 p7 60 60 150 &size 6 6 10model elasprop bulk 10e10 shear 10e10ini dens 2500apply sxx -1e9 grad 0 0 1.1111111e7 ran x -0.1 0.1apply sxx -1e9 grad 0 0 7.6666666e6 ran x 59.9 60.1apply syy -1e9 grad 0 0 8.3333333e6 ran y -0.1 0.1apply syy -1e9 grad 0 0 8.3333333e6 ran y 59.9 60.1apply szz -1e8 z -0.1 0.1fix x y z ran z -0.1 0.1set grav 0 0 -10solveini xdisp 0 ydisp 0 z disp 0ini xvel 0 yvel 0 zvel 0plot cont szz7.3 路基施工过程模拟gen zone brick p0 0 0 -50 p1 27.5 0 -50 p2 0 5 -50 p3 0 0 -10 size 8 1 10 group claygen zone brick p0 27.5 0 -50 p1 100 0 -50 p2 27.5 5 -50 p327.5 0 -10 &ratio 1.1 1 1 size 12 1 10 group claygen zone brick p0 0 0 -10 p1 27.5 0 -10 p2 0 5 -10 p3 0 0 0 ratio 1 1 0.8 &size 8 1 4 group soilgen zone brick p0 27.5 0 -10 p1 100 0 -10 p2 27.5 5 -10 p3 27.5 0 0 &ratio 1 1 0.8 size 12 1 4 group soilgen zone brick p0 0 0 0 p1 27.5 0 0 p2 0 5 0 p3 0 0 5 p4 27.5 5 0 p5 0 5 5 p6 20 0 5 &p7 20 5 5 size 8 1 5 group dam；边界条件fix x y z ran z -49.9 -50.1fix x ran x -.1 .1fix x ran x -99.9 -100.1fix ymodel mohr ran z -50 0model null ran z 0 5prop bulk 7.8e6 shear 3.0e6 coh 10e10 tension 1e10 ran group soilini dens 1500 ran group soilprop bulk 3.91e6 shear 1.5e6 coh 10e10 tension 1e10 ran group clayini dens 1800 ran group clayset grav 0 0 -9.8hist id=1 unbalsolveprop bulk 7.8e6 shear 3.0e6 coh 10e3 fric 15 ran group soil prop bulk 3.91e6 shear 1.5e6 coh 20e3 fric 20 ran group clay solveplot con szz outline onplot con sxx outline on;第一次填筑ini xdisp 0 ydisp 0 zdisp 0ini xvel 0 yvel 0 zvel 0hist id=2 gp zdis 0 0 0hist id=3 gp zdis 27.5 0 0hist id=4 gp xdis 27.5 0 0model elastic ran z 0 1prop bulk 7.8e6 shear 3.0e6 ran z 0 1 ini dens 1500 ran z 0 1solvesave fill-1.savplot con zdis outline onplot con xdis outlineonmodel elastic ran z 1 2prop bulk 7.8e6 shear 3.0e6 ran z 1 2 ini dens 1500 ran z 1 2solvesave fill-2.savmodel elastic ran z 2 3prop bulk 7.8e6 shear 3.0e6 ran z 2 3 ini dens 1500 ran z 2 3solvesave fill-3.savmodel elastic ran z 3 4prop bulk 7.8e6 shear 3.0e6 ran z 3 4 ini dens 1500 ran z 3 4save fill-4.savmodel elastic ran z 4 5prop bulk 7.8e6 shear 3.0e6 ran z 4 5ini dens 1500 ran z 4 5solvesave fill-5.savplot con zdis outline onplot con xdis outline onset log onset logfile 1.logrestore fill-1.savprint gp dis range id 517 any id 533 any restore fill-2.savprint gp dis range id 517 any id 533 any restore fill-3.savprint gp dis range id 517 any id 533 any restore fill-4.savprint gp dis range id 517 any id 533 any restore fill-5.savprint gp dis range id 517 any id 533 any set log off让土体的模量随小主应力变化的fish程序例子gen zon bri size 3 3 3model elasprop bulk 3e7 shear 1e7ini dens 2000fix x y z ran z -0.1 0.1fix x ran x -0.1 0.1fix x ran x 2.9 3.1fix y ran y -0.1 0.1fix y ran y 2.9 3.1set grav 10solveini xd 0 yd 0 zd 0 xv 0 yv 0 zv 0save 8-4.savrest 8-4.savdef E_modifyp_z=zone_headd_k=704d_n=0.38d_pa=101325.0loop while p_z# nullsigma_3 = -1.0*z_sig1(p_z)E_new = d_k*d_pa*(sigma_3/d_pa)^d_nz_prop(p_z,'young')=E_newp_z=z_next(p_z)endloopendE_modifyplot cont prop youngzone_head 与 gp_head:分别表示单元和节点的头指针z-next（）与gp_next():分别表示下一个单元（）节点单元遍历的程序框架为：p_z=zone headloop while p_z # null语句P_z=z_next(p_z)endloop节点遍历的程序框架p_gp=gp_headloop while p_gp # nullp_gp=gp_next(p_z)endloopgp_near(x,y,z):得到最靠近坐标（x,y,z）的节点地址xtable(n,s):对ID号为n的表的s行、X列进行赋值ytable(n,s):对ID号为n的表的s行、Y列进行赋值gp_zdisp(p_gp)：地址为p_gp的节点的z向变形string(n):将变量n转换成字符串的格式分级加载的施加与监测的FISH程序rest8-5.savtable 1 name load_settlementdef add_loadp_gp=gp_near(2,1,3)loop n (1 5)app_load = n*(-1000e3)file_name = '7-6_add_step'+string(n)+'.sav' commandapp nstress app_load ran z 2.9 3.1 ran x 12 y 12 solvesave file_nameendcommandxtable(1,n) = -1.0*app_loadytable(1,n) = gp_zdisp(p_gp)endloopendadd_loadsave 8-6.sav得最大位移的大小及发生位置的函数rest 8-6.savdei find_max_dispp_gp = gp_headmaxdisp_value = 0.0maxdisp_gpid = 0loop while p_gp # nulldisp_gp = sqrt(gp_xdisp(p_gp)^2 + gp_ydisp(p_gp)^2 + gp_zdisp(p_gp)^2)if disp_gp > maxdisp_valuemaxdisp_value = disp_gpmaxdisp_value = gp_id(p_gp)endifp_gp = gp_next(p_gp)endloopendfin_max_dispprint maxdisp_value maxdisp_gpid得到主应力差云图的fish程序rest 8-6.savconfig zextra 1def get_sigma_difp_z = zone_headloop while p_z_sig3(p_z)- z_sig1(p_z)z_extra(p_z,1) = sigma_difp_z = z_next(p_z)endloopendget_sigma_difplot con zextra 1例6-1简单网格的初始应力计算newgen zone brick size 3 3 3model mohrprop bulk 3e7 shear 1e7 con 10e3 fric 15 fix z ran z -0.1 0.1fix x ran x -0.1 0.1fix x ran x 2.9 3.1fix y ran y -0.1 0.1fix y ran y 2.9 3.1ini dens 2000hist unableset grav 10solve elasticsave 6-1.savrest 6-1.savini xd 0 yd 0 zd 0 xv 0 yv 0 zv 0app nstress -100e3 ran z 2.9 3.1 x 1 2 y 1 2 solvesave 6-2.savrest 6-1.savini xd 0 yd 0 zd 0 xv 0 yv 0 zv 0app nstress -100e3 ran z 2.9 3.1 x 1 2 y 1 2 hist id = 2 gp zdis 1 1 3hist id = 3 gp zdis 1 1 2hist id = 4 gp xdis 1 1 3hist id = 5 gp xdis 1 1 3hist id = 6 zone szz 1 1 3hist id 7 = zone szz 1.5 1.5 2.5hist id = 8 zone sxy 1.5 1.5 2.5solvesave 6-3.sav输出单元应力和接地啊位移rest 6-3.saavset log on &set logfile 6-2.log &print zone stress &print gp dis &set log off建立一个垂直于y轴的剖面plot set plane ori 1 1.5 0 norm 0 1 0ori(切面中一点的坐标）norm（切面的法线）生成影片文件rest 6-1.savini xd 0 yd 0 zd 0 xv 0 yv 0 zv 0app nstress -100e3 ren z 2.9 3.1 x 1 2 y 1 2 plot set rot 20 0 30plot con szz out on magf 10plot add histset movie aavi step 1 file 6-5.savmovie startsolve/////////////////////////////////prop bulk 7.8e6 shear 3.0e6 coh 10e3 fric 15 ran group soil prop bulk 3.91e6 shear 1.5e6 coh 20e3 fric 20 ran group clay solvesaveini xdisp 0 ydisp 0 zdisp 0ini xvel 0 yvel 0 zvel 0hist id=2 gp zdis 0 0 0hist id=3 gp zdis 27.5 0 0hist id=4 gp xdis 27.5 0 0model elastic ran z 0 1prop bulk 7.8e6 shear 3.0e6 ran z 0 1ini dens 1500 ran z 0 1step 5000save fill-1.savplot con zdis outline onplot con xdis outline onmodel elastic ran z 1 2prop bulk 7.8e6 shear 3.0e6 ran z 1 2 ini dens 1500 ran z 1 2step 5000save fill-2.savmodel elastic ran z 2 3prop bulk 7.8e6 shear 3.0e6 ran z 2 3 ini dens 1500 ran z 2 3step 5000save fill-3.savmodel elastic ran z 3 4prop bulk 7.8e6 shear 3.0e6 ran z 3 4 ini dens 1500 ran z 3 4step 5000save fill-4.savmodel elastic ran z 4 5prop bulk 7.8e6 shear 3.0e6 ran z 4 5 ini dens 1500 ran z 4 5step 5000save fill-5.sav。

PFC内置的FISH语言

PFC内置的FISH语言5.1 FISH语言的语法规则PFC帮助手册全文翻译。

QQ：708958918FISH是FLAC中自带的编程语言，被用户用于定义变量和新建函数，从而也拓展了FLAC的用途，增加了用户自定义功能。

例如，指定和输出新变量；编译特殊的网格生成器；将伺服控制用于数值实验；指定分布属性和自动进行参数分析等。

FISH是个“编译器”，而不仅是“翻译器”。

程序以FLAC数据文件的形式输入并转换成指令串(“伪代码”)，将指令串存储在FLAC的内存空间中，而源程序文件不存人FLAC中。

只要FISH函数被调用，伪代码就执行。

该编译代码(而不是源代码)的采用加速了程序的运行。

但FISH又不完全是编译器，变量名和变量数值可随时输出，而且用户可用FLAC中的SET命令修改数值。

FISH程序就简单的嵌在FLAC的普通数据文件中。

函数的命令行以DEFINE字符开始，END字符结束。

函数可以调用其他函数，相应的其他函数又可调用另外的函数。

只要函数在使用前都定义过，它们在函数文件中的顺序无关紧要。

因为编译后的FISH函数文件存储在FLAC的内存中，所以使用SAVE命令就可以保存函数和相关变量的当前值。

以下介绍FLAC的语法规则。

FISH程序能植入普通FLAC数据文件中执行或直接从键盘输入，FISH函数以DEFINE 字符开始，在END字符处结束，一个有效的FISH代码必须具有以下格式之一。

代码行用语句(statements)开始，比如，IP，LOOP等。

命令行含多个用户自定义函数时，函数名以空格间隔，如：fun—1 fun—2 fun—3以上名字和用户定义的函数名一致，函数即可依序执行。

FISH代码中函数定义无需先于函数声明，即，允许提前声明函数。

命令行包含赋值语句(将“：”右边的符号赋给左边的变量或函数名)。

代码行含有FLAC命令时，假定命令嵌在FISH代码中，则必须用COMMAND--ENDCOM-MAND命令与FISH隔开。

FLAC教程

(5)FLAC运动方程(含惯量项)的显式时间逼近解法允许进行岩体的渐进破坏与跨落。摩擦材料剪切带的形成与定位以及工程材料的大变形分析等。 (6)在求解过程中，FLAC又采用了离散元的动态松驰法，不需求解大型联立方程组(刚度矩阵)，便于微机上实现。 (7)它不但能处理一般的大变形问题，而且能模拟岩体沿某一软弱结构面产生的滑动变形。

2.1 FLAC简介

2.1.1 什么是FLAC？ FLAC其实是快速拉格朗日分析 (Fast Lagrangian Analysis of Continua)的简称，是一种用于工程力学计算的显式有限差分程序。该程序可模拟土、岩石等材料的力学行为。它利用拉格朗日有限差分原理编制，以牛顿第二定律为基础，可以对连续介质进行非连续大变形分析，特别适用于求解岩土工程中的大变形问题。 2.1 .2 FLAC的发展 (1).由美国ITASCA咨询集团公司开发，首先由Cundall 在20世纪80年代提出并将其程序化、实用化。 (2).在20世纪90年代，得以广泛应用。我国20世纪90 年代初才引进该软件。
FLAC中可以模拟的模型

(1)零空模型(Mull)：代表网格中的孔洞(开挖单元)； (2)应变硬化/软化模型(SS)：代表非线性，不可逆剪切破碎与压缩； (3)粘弹性蠕变模型； (4)界面模型(界面为平面，沿界面允许滑动和分开)：模拟断层、节理和摩擦边界； (5)水利模型：模拟可变形空隙体与粘性流动的全藕荷； (6)结构单元模型：模拟岩土体加固、衬砌、锚杆、混凝土喷层、可缩支柱及钢拱等。
生成网格
gen -3 3 -3 30 3 30 3 3
gen 3 -30 3 -3 30 -3 30 -30 gen 3 -3 3 3 30 3 30 -3 gen 3 3 3 30 30 30 30 3

FLAC若干问题的解答

FLAC若⼲问题的解答把结果中的信息输出到指定⽂件分为两种：1、把单元信息输出到指定⽂件set log onset log on finame.dat;可以输⼊⽂件路径，否则按当前⽬录处理print zoneset log off通过反复使⽤该命令，可以把不同信息输出到不同的⽂件。

2、把节点信息输出到指定⽂件set log onset log on finame.dat;可以输⼊⽂件路径，否则按当前⽬录处理print gpset log off请教⼀下print gp⾥的keyword该如何填。

⽐如想输出某点的zdis谢谢print gp dispprint gp positionlog ⽂件中有，然后处理，这些论坛⾥都有的关于flac3d—内置fish语⾔精讲FISH语⾔是FLAC3D程序的内置编程语⾔因为FLAC3D的最佳操作⽅式是命令流⽂件⽅式这⼀点与ANSYS很相似，⽽FISH就相当于ANSYS的APDL语⾔。

它包括循环、判断等结构。

如果你还⽤过其它⾼级语⾔，那么从形式上讲你也可以把它理解为⼦函数。

FISH语⾔的引⼊极⼤的⽅便了⽤户进⾏复杂的程序建模它不但可以嵌⼊命令流⽂件⾥⼯作⽽且还可以引⽤FLAC3D本⾝的任何命令所以说它实现了对FLAC3D的完全控制引⼊FISH语⾔的作法值得其它通⽤软件效仿。

以上仅是个⼈的使⽤体会请版主与管理员指教也供初学者⼊门之⽤。

谢谢！Fish函数增加了以下⼏种新的特性：1.增加了fish变量来获取结点、单元和界⾯变量2.fish提供了获取结构单元变量的途径3.休单元和⾯单元性质⽬前可以通过单元变量名z_prop(i_z string)和界⾯单元名i_prop(i_z string)分别加以识别4.fish函数可以获取单元应变和应⼒速率，还提供了全应变增量张量和应变速率张量5.提供了fish绘图⼦程序函数够⽣⽤户定义的图形内容6.fish函数已经增加了从⽂件读、写数据的能⼒FLAC3D是⼀个强⼤的软件，但是不得不承认的，它的界⾯没有ansys好⽤，商业运⾏不是很强的说也来凑个热闹，呵呵注意点：1.fish函数可以嵌套使⽤;2.以save命令保存模型时，fish函数和变量也同时保存；3.fish函数不⽀持缩写，这与flac3d命令不同，另外所有的fish函数或变量不区分⼤⼩写，程序同意转化为⼤写进⾏编译，当然也可以通过执⾏set case-sensitivity on来区分⼤⼩写；4.变量或函数名不能以⼀个数字开头也不能是下列字符：. , * / ^ = > < # ( ) [ ] @ ; “ '5.如果⽤命令set safe on 指定了编译安全模式，则⽤户调⽤fish函数时，函数名前必须加@；6.如果变量不曾赋值，则系统默认为零（整形），如果赋值，其类型由值的类型决定；7.fish函数的调⽤⽅法：.可以出项在其他fish函数的单独⾏中；.可以出现在其他fish函数的表达式中；.出现在flac3d的命令⾏中；.作为命令set,print,hist的参数。