FLAC3D数值模拟基础

FLAC3D数值模拟基础刘升贵中国矿业大学力学系liushg2002@163.com主要内容FLAC3D软件简介1、基坑开挖2、浅基础的稳定性分析3、网格的合并联结4、界面的生成5、隧道的生成6、模型材料问题7、Interface合并（联结）问题8、初始条件问题9、破坏问题分析10、综合实例-煤巷应力分析2主要内容FLAC3D软件简介1、基坑开挖2、浅基础的稳定性分析3、网格的合并联结4、界面的生成5、隧道的生成6、模型材料问题7、Interface合并（联结）问题8、初始条件问题9、破坏问题分析10、综合实例-煤巷应力分析3FLAC3D简介FastLagrangianAnalysisofContinua美国Itasca（依泰斯卡）咨询公司开发2D程序(1986)1990年代初引入中国有限差分法(FDM)DOS版→2.0→2.1→3.0Itasca其他软件4FLAC3D简介应用：岩土力学/岩石力学分析，例矿体滑坡、煤矿开采沉陷预测、水利枢纽岩体稳定性分析、采矿巷道稳定性研究等岩土工程、采矿工程、水利工程、地质工程特色：大应变模拟完全动态运动方程使得FLAC3D在模拟物理上的不稳定过程不存在数值上的障碍显示求解具有较快的非线性求解速度51承受荷载能力与变形分析：用于边坡稳定和基础设计2渐进破坏与坍塌反演：用于硬岩采矿和隧道设计3断层构造的影响研究：用于采矿设计4施加于地质体锚索支护所提供的支护力研究：岩锚和土钉的设计5排水和不排水加载条件下全饱和流体流动和孔隙压力扩散研究：挡土墙结构的地下水流动和土体固结研究6粘性材料的蠕变特性：用于碳酸钾盐矿设计7陡滑面地质结构的动态加载：用于地震工程和矿山岩爆研究8爆炸荷载和振动的动态响应：用于隧道开挖和采矿活动9结构的地震感应：用于土坝设计10由于温度诱发荷载所导致的变形和结构的不稳定：高辐射废料地下埋藏的性能评价12大变形材料分析：用于研究粮仓谷物流动及井巷和矿洞中材料的总体流动FLAC3D简介6基本原理开始生成网格并调整网格的形状；持续的运动和连续的物质属性；特定的边界条件和初始条件。到达平衡状态效果是否符合要求作如下改变：开挖模型的物质属性改变边界条件实施求解效果是否符合要求参数是否调整结束网格确定问题的几何尺寸；持续的运动和连续的物质属性决定了模型的扰动（如由于开挖引起的变形）形式；边界条件和初始条件确定了模型的初始状态（没有引起扰动或变形的状态）。建立FLAC计算模型，必须进行以下三个方面的工作：

1.有限差分网格

2.本构特性与材料性质

3.边界条件与初始条件

完成上述工作后，可以获得模型的初始平衡状态，也就是模拟开挖前的原岩应力状态。然后，进行工程开挖或改变边界条件来进行工程的响应分析。

7基本原理有限差分法Lagrangian网格空间混合离散技术Lagrangian格式动量平衡方程FLAC3D的求解过程FLAC3D的本构模型8有限差分法古老的方法(上世纪40年代)用差分格式转化控制方程中的微商格式流体力学；土工渗流问题；固结FDM&FEM的混合求解FDM的新进展9Lagrangian网格源自流体力学中的拉格朗日法跟踪流体质点的运动状态跟踪固体力学中结点，按时步用Lagrangian法研究网格节点的运动节点和单元随材料移动，边界和接触面与单元的边缘一致固体力学大变形理论法国数学家、物理学家拉格朗日

10FLAC3D的求解过程平衡方程(动量方程)应力—应变关系(本构模型)Gauss定律单元积分应变率速度节点力新的应力对所有的网格节点对所有单元11FLAC3D中的本构模型开挖模型null3个弹性模型各向同性弹性横观各向同性弹性正交各向同性弹性8个塑性模型（Drucker-Prager模型、Morh-Coulomb模型、应变硬化/软化模型、遍布节理模型、双线性应变硬化/软化遍布节理模型、修正剑桥模型和胡克布朗模型）12FLAC3D中的本构模型13FLAC3D的前后处理术语水平边界压力结构线模型边界格网点格网区域水平边界压力滚动底端边界内部开挖边界区域（Zone）……有限差分划分的带在几何上是最小的区域，在在这个区域里的每一个现象的变化，如应力应变都可以估计出，。各种形状的多面体（立方体、楔形、锥体、四面体等）可用来构造模型并可用plot显示出来。每一个多面体可能有一套或两套表层设置，这由5个四面体组成。默认的情况下，两个表层设置用在对计算精度要求高的情况下，区域的另外一种叫法是要素。栅格点（GridPoint）……栅格点是有限差分单元的角点。一个多面体可能有5个、6个、7个或8个网格点，主要取决于多面体的形状。给定每个节点的x，y和z值这样就具体确定了有限差分单元，。其他叫法有：节点，交点。有限差分栅格（FiniteDifferenceGrid）……有限差分网格是研究区域中一个或多个通过物理边界连接的有限差分单元的集合。另一个叫法是网格，有限差分网格也可以标识出模型中每个状态的存储位置，FLAC3D所生成的矢量都保存在节点上（如：受力、速度、位移）。标量和张量保存在单元的中心（如应力、材料属性）。14FLAC3D的前后处理术语水平边界压力结构线模型边界格网点格网区域水平边界压力滚动底端边界内部开挖边界模型边界（ModelBoundary）……即有限差分网格的外围，内部边界也同样是模型边界（如网格中的空洞）。边界条件（BoundaryCondition）……即模型边界的约束条件或控制条件的给定（如：限制位移、渗透条件、绝热条件）。初始条件（InitialConditions）……即在对模型加载或开挖等作用前的各种参数状态。基本模型（ConstitutiveModel)……基本模型（材料模型）即规定了FLAC3D模型中某一区域的变形或强度效应，可用大量基本模型去近视地质材料，可以单独定义FLAC3D模型中的基本模型和材料模型。空单元（NullZone）……空单元表示此区域为空（就象没有材料一样）。次级网格（SUB-GRID）……有限差分网格可由次级网格组成，它可用来在模型中创建不同形状的区域，次级网格是分别生成，可进行合并和连接。15FLAC3D的前后处理术语水平边界压力结构线模型边界格网点格网区域水平边界压力滚动底端边界内部开挖边界附属接触面ATTACHEDFACES……附属栅格面是由被划分的次级栅格组成的网格接触或合并的面栅格面，接触面必须是共面或接触，每个面的节点不一定一样，不同总密度的次级网格可以接触。接触面INTERFACE------即次级网格在计算过程中可以分开（滑动，开裂）的两部分之间的面，可表示不连续的物理特征，如，断层，节理面或材料性质突变的临界面。范围INTERFACE------范围是对FLAC3D模型空间值的一个描述，可给定一个命令的作用范围，即使模型发生运动---不影响模型中的区域和节点位置，一个范围或范围确定的单元也不发生改变。范围或范围内的单元所包含的区域，也可以用一个单元的ID号来表示，它与区域，节点，或结构单元密切联系。集合GROUP-----在FLAC3D模型中他们有共同的名称，由于限定具体命令的对象，如model命令对某一集合设置为某种材料，任何命令加于集合名称也就相当于作用于这一集合的所有区域。16FLAC3D的前后处理术语水平边界压力结构线模型边界格网点格网区域水平边界压力滚动底端边界内部开挖边界ID号码IDNUMBER-----FLAC3D模型中的单元以ID号加以区分，下面的单元有ID号；内部面、节点、区域、体积、历史、表格、显示项和结构单元的全部内容。这帮助用户确定模型中的单元，可用porint命令获得ID号，用户可给内部面、结构单元、历史等赋ID号。实体结构单元同样也有CID号，系统给每一个网格，单元都创建了一个CID号，这与梁，柱等不一样。结构单元STRUCTURALELEMENT------在FLAC3D中有两种结构单元。二结点，线性单元表示梁，柱作用。三结点，三角平面单元表示面状，结构单元用来模拟土体或岩体中结构支护的相互作用。非线性材料作用可用单元表示。每一结构单元实体（梁，柱，面体）包括三个内容：结点、单个单元（也叫sels）和网格连线，这些内容的不同可区别出梁、桩、面体的作用。步STEP------因为在FLAC3D是具体代码，问题的计算须分步进行，随步长的增加，现象的有关信息在研究区域传递。对于静态分析，需要给一个具体的步，让其达到平衡状态，典型的问题计算在2000-4000步之间，其他叫法有时间步，循环次。17FLAC3D的前后处理术语水平边界压力结构线模型边界格网点格网区域水平边界压力滚动底端边界内部开挖边界静态解答STATICSOLUTION-----如果模型中动量变化率小于了某一可忽略的值，就认为静态或类静态出现了，这通过限定运动方程实现，静态就是模型达到应力平衡，或流体材料受外力后从不稳定到稳定。这种分析方法在FLAC3D中是默认的分析方法，机械的静态分析也可与地下水渗透或热传递问题结合（通过特定设置后，动态问题可由带约束的静态分析代替）。非平衡力UNBALANCEDFORCE-------非平衡力标征静态分析达到机械平衡（或塑性变形前），严格的说平衡时每个节点上的应力矢量都为0。最大应力会自动被监测，当击活step或solve命令时，其值会显示在屏幕上。最大网格力也叫非平衡力或抗平衡力，非平衡力在数值上永远也不能达到0，当最大非平衡力相对加载的力很小时，我们就认为模型达到了平衡状态，如果非平衡力一直保持某一非0值，这就说明模型中可能发生了破坏或塑性变形。18FLAC3D的前后处理术语水平边界压力结构线模型边界格网点格网区域水平边界压力滚动底端边界内部开挖边界动态解答DYNAMICSOLUTION------在动态分析中，求解所有动态方程，动量的产生和消耗都将产生直接影响，在高频率发生或持续时间很短的过程中用到，如地震或爆炸。动态计算是FLAC3D的一个可选模块。（见附录K）大应变/小应变LARGESTRAIN/SMALLSTRAIN--------FLAC3D默认的都是小应变，也就是说，即使计算出来的位移相对通常区域尺寸很大，节点也不发生相应位移。大应变中，节点根据每一步计算出的位移量发生位移。几何非线性只有用大应变才能实现。19FLAC3D的前后处理命令驱动(推荐)程序控制图形界面接口计算模型输出指定本构模型及参数指定初始条件及边界条件，指定结构单元指定接触面指定自定义变量及函数(FISH)求解过程的变量跟踪进行求解模型输出20菜单驱动(计算模式)命令栏21菜单驱动(Plot)22一个最简单的例子genzonbrisize333;建立网格modelelas;材料参数propbulk3e8shear1e8inidens2000;初始条件fixzranz-.1.1;边界条件fixxranx-.1.1fixxranx2.93.1fixyrany-.1.1fixyrany2.93.1setgrav00-10solve;求解appnstr-10e4ranz3x12y12solveRUNFLAC3D23前后处理功能的优点多种zone类型后处理快捷、方便、丰富计算过程中的hist变量动态显示FISH可进行参数化模型设计单元状态的可编程计算暂停时的后处理与可保存24前后处理功能的缺点复杂模型的建模功能不强可以编程导入其他软件形成的网格(比如：Ansys、Adina、GeoCAD)无等值线的后处理功能(3D)可编程将.sav文件写入TecPlot等其他后处理软件全命令操作，学习困难鼠标功能单一(双击取击点坐标)25用tecplot绘制曲线

1.第一主应力

2.xdisp、ydisp、zdisp、disp

用excel做曲线

隧道

1做地表沉降槽（zdisp）

2地表横向位移（xdisp）

3隧道中线竖向沉降曲线（zdisp）

4提取位移矢量图，

5显示初期支护结构内力

6显示state（找塑性区）

基坑

1做地表沉降槽（zdisp）

2提取位移矢量图，

3显示初期支护结构内力

4显示state（找塑性区）边坡做安全系数和应变图后处理26需要掌握gen,ini,app,plo,solve等建模、初始条件、边界条件、后处理和求解的命令。

常用命令FLAC3D程序的编写步骤：1Config________2Grid________3Model________4求起始的应力平衡________

(1)建立x,y坐标与网格的关系，建议使用Gen指示：Genx1,y1x2,y2,x3,y3x4,y4i=i0,i1j=j0,j1

(2)设定材料性质：prop

(3)设定外力：SetGrav,ApplyPressure,inisxx,Syy

(4)设定边界条件：fix,free

(5)求起始的应力平衡：solve

(6)储存：Save5求工程的影响________求出区域内的应力分布情况后，再依工程的流程及步骤阶段执行各工程进行过程的影响，建议使用以下的步骤：

(1)调出起初的应力平衡：re_____.sav

(2)设定新的材料性质：model,prop

(3)设定新的支撑性质：struct

(4)设定新的外力

(5)设定边界条件

(6)求工程时的应力平衡

(7)储存27主要内容FLAC3D软件简介1、基坑开挖2、浅基础的稳定性分析3、网格的合并联结4、界面的生成5、隧道的生成6、模型材料问题7、Interface合并（联结）问题8、初始条件问题9、破坏问题分析10、综合实例-煤巷应力分析28Case-1：

土体中挖了一个长宽为2m×4m，深为4m的沟，挖土体的同时监测周围土体的变形情况（模型6m×8m×8m）初始化不同限定的格网，可以使用GENERATE（生成）命令：genzonebricksize688这个命令会建立以一个初始化的格网，这个格网在X方向上有6个分区，Y方向上有8个分区，Z方向有8个分区。所建模型的Z轴在垂直方向上。

第一步：初始模型的建立显示模型格网

命令plot

CreateTrenchAddsurfaceyellowAddaxesblackShow创建一个名为“Trench”的块视图并将块表面设为黄色，把视图块的轴设为黑色。show关键字是显示当前视图——自动地创建一个视图并设置为当前视图。模型的正面透视图会出现在窗口视图中。按住键盘的X、Y、Z键，图形能X、Y、Z方向旋转。使用M键可以放大视图，使用箭头键可以移动图形块（shift键然后按住相应的键可以反转和缩小图形

显示的图形块中的标题中给出了“Center:”，“Rotation:”，“Dist:”和“Mag:”的数据

29Case-1：

土体中挖了一个长宽为2m×4m，深为4m的沟，挖土体的同时监测周围土体的变形情况（模型6m×8m×8m）创建另一个视图——黑体部分显示了新增加的命令。（我们假定Plot>命令提示符仍然是被激活的，如果没有，在键入命令之前键入PLOT然后按回车键。）

CreateTrench2AddsurfaceyellowAddaxesblackSetrotation17.4959.84728.481Setcenter344Setdist26.12Setmagnification0.8Show第一步：初始模型的建立30Case-1：

土体中挖了一个长宽为2m×4m，深为4m的沟，挖土体的同时监测周围土体的变形情况（模型6m×8m×8m）使用莫尔-库仑(Mohr-Coulomb)准则弹塑性模型

（在MODEL命令中不指定区域的范围，FLAC3D假定所有的区域都是Mohr-Coulomb性质）

Propbulk=1e8shear=0.3e8fric=35coh=1e10tens=1e10Propcoh=1e10tens=1e10

注意：关键参数的值可能被空格键或“=”分开

体积、剪应力、内摩檫角、内聚力和抗压强度是一定的。我们所见的内聚力和抗压强度都给得很大，这仅仅是在重力作用阶段给材料的初始值。实际中，在初始应力阶段材料是弹性的。

setgrav0,0,-9.81inidens=1000//初始化（ini=initial）密度为1000Kg/m3在z的负方向加以9.81m/sec2的重力加速度（在坐标轴的正方向把重力视为正），为了给模型一个重力，材料密度需预先给定。命令INI用于把模型的每个区域的平均密度设为1000Kg/m3

第一步：初始模型的建立-模型属性31Case-1：

土体中挖了一个长宽为2m×4m，深为4m的沟，挖土体的同时监测周围土体的变形情况（模型6m×8m×8m）边界条件也可在FLAC3D的提示符下输入：

fixxrangex–0.10.1fixxrangex5.96.1fixyrangey–0.10.1fixyrangey7.98.1fixzrangez–0.10.11、前两个命令规定，沿着在x=0和x=6两个平面的栅格边界上的节点在x方向被“固定”，这两个边界面通过Fix命令在”range”限定的范围内下降；2、在y=0和y=8上的节点在y方向上被固定。第3、4条命令规定了他们的下降范围；3、沿着底面边界上的节点在z方向上被固定（z=0），这个平面被第5条命令固定了范围。

第一步：初始模型的建立-边界条件32Case-1：

土体中挖了一个长宽为2m×4m，深为4m的沟，挖土体的同时监测周围土体的变形情况（模型6m×8m×8m）在计算过程中我们想监测所选参数的变化，History命令可以帮助我们判断是处于平衡状态还是发生了不稳定的破坏。

histn=5//hist=history设置步长为5histunbal//最大非平衡力histgpzdisp448//测试栅格点448的位移我们设定监测参数的变化步长为5（默认步长是10），则每隔5步被选定的参数值就会自动保存在历史列表中，被保存的2个参数是：最大非平衡力和节点z方向上的位移（x=4,y=4,z=8点）；知道一个模型中的最大非平衡力很有用，如果最大非平衡力的值很小或位移很稳定，这说明模型已达到平衡状态。

为了考查模型的重力作用，先规定达到平衡的步数。Solve命令可自动判断平衡，*（计算时间，不是实际时间）当键入：

setmechforce=50solve计算开始，步数和最大非平衡力将在屏幕上显示；当最大非平衡力小于极限值（set命令给定的大小为50N的力）时，程序将停止运行。

第一步：初始模型的建立-监测网格点33Case-1：

土体中挖了一个长宽为2m×4m，深为4m的沟，挖土体的同时监测周围土体的变形情况（模型6m×8m×8m）可以根据最大主应力和z方向位移的历史记录判断它是否达到平衡。当键入：plothist1

可得到最大非平衡力的历史记录

最大非平衡力的历史记录图

z方向上位移历史记录图(在x=4,y=4,z=8这点上)

hist2

当最大非平衡力的历史值达到0，而且位移的历史记录稳定。这两个方面都可以反应出模型已达到了平衡状态。

得到z方向的位移历史记录;注意：通过hist命令，每一个历史都被连续地从一开始进行数值化，回到FLAC提示符下，键入；printhist得到一个历史列表和相应的数值。

第一步：初始模型的建立-监测网格点34Case-1：

土体中挖了一个长宽为2m×4m，深为4m的沟，挖土体的同时监测周围土体的变形情况（模型6m×8m×8m）同样可以在模型中考察重力效应。现在我们建立一个包含多个图形的复杂窗口，只需键入：plotcreateGravVplotsetplanedip=90dd=0origin=3,4,0plotsetrot15020;thiswouldbeachievedinteractivelyplotsetcenter;andisshowheresimplytoillustrateplotaddboundbehindplotaddbcontszzplaneplotaddaxesplotshow这样我们就建立了一个叫“GravV”的图形，将它设为当前图层，对此图形设置一个倾角为90°的面，（在x-y面上，约定z的负方向为下），倾向为0°。（在x-y面上，从y的正方向开始沿顺时针方向计量），过x-y面上一点（x=3,y=4,z=0）,在平面后加一个网格，垂直方向的应力Qzz，最后给模型加轴线，整体图形轮廓与添加图形相反，显示每个区域中心的应力计算值，每个区域的颜色直接对应基本区域的应力，当键入show,我们建立的图形被显示

在初始平衡时的重力应力等值线图

第一步：初始模型的建立-监测网格点35Case-1：

土体中挖了一个长宽为2m×4m，深为4m的沟，挖土体的同时监测周围土体的变形情况（模型6m×8m×8m）保存初始状态，你可以在进行参数分析的时候提取；保存此状态后，回到FLAC3D的提示符下键入；savetrench.Sav于是在默认路径下将建立一个叫“trench.Sav”的文件。

现在在模型中建立两个面，键入命令：plotprint就可以显示这些视图的列表，动态图将会用箭头标出，图形“base”为草稿视图。键入命令：plotprintview将显示当前图形的详细信息。在以上命令中加入图形名称或数字将产生一个图形的详细描述，包含所有设置和项目在图上。当我们把当前的“GravV”图形改变为“trench”图形时，只要通过输入命令：plotcurrentTrench//Trench视图必须存在（通过create命令生成）这样就让“Trench”接受项目或设置的改变，任何现存的视图都可直接显示出来，而不用先把它设为当前。plotshowTrench将在主菜单下显示所有可视的图形

第一步：初始模型的建立-监测网格点36Case-1：

土体中挖了一个长宽为2m×4m，深为4m的沟，挖土体的同时监测周围土体的变形情况（模型6m×8m×8m）第二步开挖巷道

在”plot”下，现在我们在土体内开挖巷道，先键入：propfric=35coh=1e3tens=1e3摩擦角35°粘聚力1×103抗拉强度1×103这就给整个区域设置了一个大小为1000Pa的内聚力和抗拉强度，这个强度足以防止在初始阶段发生破坏（即开挖前），但是你要选择一个步数，使它计算完后有破坏的可能性。进行开挖，键入：modelnullrangex=24y=26z=510开挖的巷道尺寸（即空材料），通过x-,y-,z-定义一个合适的范围。一个低内聚力和无垂直支护的巷道将发生跨塌，因为我们想仔细考察这一过程，应使用大应变模式。通过键入以下命令实现：setlarge为了得到更好的显示效果，我们只想看巷道开挖后的位移变化，而不想从重力加载就开始观察，所以我们去掉刚开始x,y,z方向上由于重力发生的位移；*inixdis=0ydis=0zdis=0我们有意设置足够小的内聚力让其破坏，我们现在不用“solve”命令计算，因为用他时要给一个失去平衡的应力极限（为了平衡而选择的），而我们的模拟永远都不会达到平衡；此外，我们可以给定一个时间步，当此时跨塌时显示其结果，这只是一个间接方法，模拟不必在每个计算过程中都达到平衡，因为我们用不着解决一系列数字式

37Case-1：

土体中挖了一个长宽为2m×4m，深为4m的沟，挖土体的同时监测周围土体的变形情况（模型6m×8m×8m）5、开挖巷道

用“step”命令step2000此时FLAC3D将计算2000步。这样考察计算主要通过显示一些结果。比如，为了产生出一个面的彩色等值线而重复“plot”。但要显示位移等值线图，键入plotcreateDispcontcopyGravVDispcontsettingaddcontdispplanebehindaddaxesshow

位移等值线图可显示出由于开挖引起的下沉量

38Case-1：

土体中挖了一个长宽为2m×4m，深为4m的沟，挖土体的同时监测周围土体的变形情况（模型6m×8m×8m）39主要内容FLAC3D软件简介1、基坑开挖2、浅基础的稳定性分析3、网格的合并联结4、界面的生成5、隧道的生成6、模型材料问题7、Interface合并（联结）问题8、初始条件问题9、破坏问题分析10、综合实例-煤巷应力分析40Case-2：浅基础的稳定性分析

条形基础是指基础长度远远大于宽度的一种基础形式。按上部结构分为墙下条形基础和柱下条形基础。基础的长度大于或等于10倍基础的宽度。条形基础的特点是，布置在一条轴线上且与两条以上轴线相交，有时也和独立基础相连在此，以条形基础作为研究对象进行说明41Case-2：浅基础的稳定性分析

42Case-2：浅基础的稳定性分析

可以用平面应力分析方法解决条形基础问题土体参数二维的条形基础的塑性流动(Tresca材料)

（模型20m×1m×10m）43Case-2：

二维的条形基础的塑性流动(Tresca材料)

（模型20m×1m×10m）44Case-2：

二维的条形基础的塑性流动(Tresca材料)

（模型20m×1m×10m）genzonebricksize20110;建立矩形的网格区域，其大小为20×1×10plotCreateTrenchAddsurfaceyellowAddaxesblackShow;显示模型格网pausemodelmohr;采用摩尔-库仑模型propbul2.e8shea1.e8cohesion1.e5;模型的材料性质：体积模量2.0×108;剪切模量1.0×108粘聚力1.0×105propfriction0.dilation0.tension1.e10;内摩擦角0°剪胀角0°抗拉强度1.0×1010fixxrangex-.1.1;在x方向上固定边界x=0（即只在x方向上受约束）fixxyzrangez-.1.1;在xyz方向上固定边界z=0（即为固定约束）fixxyzrangex19.920.1;在xyz方向上固定边界x=20fixy;约束所有y方向上运动fixxyzrangex-.13.1z9.910.1;在xyz方向上固定平面x=0，3z=10inizvel-0.5e-5rangex-.13.1z9.910.1;在平面x=0，3z=10上初始化z方向的速度为−0.5×10−545Case-2：

二维的条形基础的塑性流动(Tresca材料)

（模型20m×1m×10m）defp_cons;用fish语言定义函数p_conspdis1=gp_near(0.,0.,10.);靠近坐标（0，0，10）的节点的地址赋予pdis1pdis2=gp_near(0.,1.,10.);靠近坐标（0，0，10）的节点的地址赋予pdis2p_sol=(2.+pi);2+π赋予p_solend;用于fish语言的结尾p_cons;运行函数p_cons;---------------------------------------------------------------------;p_load:averagefootingpressure/c;场地的平均压力;c_disp:magnitudeofverticaldisplacementatfootingcenter/a;场地中心的垂直沉降量;---------------------------------------------------------------------defactLoad;定义函数actLoad（活荷载）（fish语言的起始句）46Case-2：

二维的条形基础的塑性流动(Tresca材料)

（模型20m×1m×10m）pnt=gp_head;把网格点指针的首个地址（gp_head）赋予变量pntpload=0.0;定义变量pload为0n=0;定义变量n为0loopwhilepnt#null;loop是循环命令，其条件是pnt≠0ifgp_zpos(pnt)>9.9then;当变量pnt所对应的z坐标大于9.9ifgp_xpos(pnt)<3.1then;且x坐标小于3.1时pload=pload+gp_zfunbal(pnt);变量pload为上次循环中的值与变量pnt所对应;的z方向上的不平衡力的和n=n+1;n自我加1endif;结束里面的ifendif;结束外面的ifpnt=gp_next(pnt);把网格点指针的下一个地址（gp_next）赋予变量pntendloop;结束循环actLoad=pload/(3.5*z_prop(zone_head,'cohesion'));这里函数actLoad的值，为后面的式子所赋予，;其中pload为循环结束后的最终;值z_prop(zone_head,‘cohesion’)为首个单元;的粘聚力的数值c_disp=-(gp_zdisp(pdis1)+gp_zdisp(pdis2))/7.0;同理，这里函数c_disp的值，;也为后面的式子所赋予，其中gp_zdisp(pdis1)为;变量pdis1所对应的网格点在z方向的位移，;gp_zdisp(pdis2)变量pdis2所对应;的网格点在z方向的位移end;fish语言的结束句47Case-2：

二维的条形基础的塑性流动(Tresca材料)

（模型20m×1m×10m）defp_err;定义函数p_errp_err=100.*(actLoad-p_sol)/p_sol;函数p_err的值为后面式子所定义end;fish语言的结束句histn50;对n做历史记录histactLoad;对actLoad做历史记录histp_sol;对p_sol做历史记录histc_disp;对c_disp做历史记录histunbal;对unbal（不平衡力）做历史记录plotsketch;绘出轮廓线plotaddvel;添加矢量plotaddhis12vs3;绘出历史记录n，actLoad，p_solcyc15000;运行15000时步printp_err;输出函数p_errsavepran.sav;形成sav文件，并保存为pran.sav48Case-2：

二维的条形基础的塑性流动(Tresca材料)

（模型20m×1m×10m）49Case-2：

二维的条形基础的塑性流动(Tresca材料)

（模型20m×1m×10m）50主要内容FLAC3D软件简介1、基坑开挖2、浅基础的稳定性分析3、网格的合并联结4、界面的生成5、隧道的生成6、模型材料问题7、Interface合并（联结）问题8、初始条件问题9、破坏问题分析10、综合实例-煤巷应力分析51Case-3：网格的合并联结plotcreategrid4plotaddblockgroupredyellow;把在group中的部分染成红色和黄色plotaddaxesblack;添加坐标轴线为黑色pausegenzonebricksize888p0-10,-10,-20p110,-10,-20&p2-10,10,-20p3-10,-10,0&grouplower_block;建立块体网格，大小为8×8×8，其尺寸是p0、p1、p2、p3来决定，其中p0为;起始点的坐标，p1为x方向，p2为y方向，p3位z方向的坐标，并把这些区域;建立成一个群，名为lower_block。genzonebricksize2,2,1p00,0,1p15,0,1&p20,5,1p30,0,3&groupupper_block;建立另外一个块体网格，其大小为2×2×1，其尺寸及位置是p0、p1、p2、p3;来决定，具体如上所述，并把这个区域归为一个群，名为upper_block。plotshow;在视图窗中显示pause;暂停inizadd-1rangegroupupper_block;群upper_block的所有单元，在z方向上向下移动1m;genzonebricksize2,2,1p00,0,0p15,0,0&;p20,5,0p30,0,2&;groupupper_block;重新建立另外一个块体网格，其大小为2×2×1，其尺寸及位置是p0、p1、;p2、p3来决定，具体如上所述，并把这个区域归为一个群，名为upper_block。genmerge1e-5rangez0;此命令是接触面单元合并成一个整体。52主要内容FLAC3D软件简介1、基坑开挖2、浅基础的稳定性分析3、网格的合并联结4、界面的生成5、隧道的生成6、模型材料问题7、Interface合并（联结）问题8、初始条件问题9、破坏问题分析10、综合实例-煤巷应力分析53Case-4：界面的生成;CreateBasegenzonebricksize333p0(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)groupBase;建立另外一个块体网格，其大小为3×3×3，其尺寸及位置是p0、p1、p2、p3、;p4、p5、p6和p7等楔型体的8个角点坐标来决定，并把这个区域归为一个群，;名为Base（即下面的底座）。;CreateTop-1unithighforinitialspacinggenzonebricksize333p0(0,0,1.5)p1(3,0,4.5)p2(0,3,1.5)p3(0,0,6)p4(3,3,4.5)p5(0,3,6)p6(3,0,6)p7(3,3,6)groupToprangegroupBasenot;建立另外一个块体网格，其大小为3×3×3，其尺寸及位置是p0、p1、p2、p3、;p4、p5、p6和p7等楔型体的8个角点坐标来决定，并把这个区域归为一个群，;名为Top（即上面的部分）。genseparateTop;使两部分的接触网格分离为两部分;Createinterfaceelementsonthetopsurfaceofthebase;在两部分之间添加交界面单元interface1wrapBaseTop;在这两部分之间添加接触单元;plotcreateview_int;显示，并创建标题view_intplotaddsurface;显示表面plotaddinterfacered;显示交互面为红色plotshow;打开图形saveint.sav;形成sav文件，并保存为int.sav54主要内容FLAC3D软件简介1、基坑开挖2、浅基础的稳定性分析3、网格的合并联结4、界面的生成5、隧道的生成6、模型材料问题7、Interface合并（联结）问题8、初始条件问题9、破坏问题分析10、综合实例-煤巷应力分析55Case-5：隧道的生成new;新建defsetup;定义函数setupnumy=8;定义常量numy为8depth=10.0;depth为10end;结束对函数的定义setup;运行函数setupgenzoneradcylsize5numy810&p0000p11002p20depth0p30010&p410depth2p50depth10p610010p710depth10;建立模型的单元为放射性圆柱，其内部圆柱大小尺寸为5×unmy×8，其总体的;具体的空间位置由p0,p1,p2,p3,p4,p5,p6,p7（它们的具;体意思参看generate命令）决定。plocreaqqqploaddsurfredploshopause;;暂停genzoneradcylsize5numy810&p0000p1100-10p20depth0p31002&p410depth-10p510depth2p6100-4p710depth-4;建立模型的单元为放射性圆柱，其内部圆柱大小尺寸为5×unmy×8，其中数值;10的意义是扩展网格的个数。其总体的具体的空间位置由;p0,p1,p2,p3,p4,p5,p6,p7（它们的具体意思参看generate命令）决定。plocreaqqqploaddsurfredploshopause;genzonebricksize5numy10&p000-10p1100-10p20depth-10p300-2p410depth-10&p50depth-2p620-2p72depth-2;建立块体单元网格，其大小为5×unmy×10，;空间位置由p0,p1,p2,p3,p4,p5,p6,p7;（它们的具体意思参看generate命令）决定。genzonereflectorig000norm100;对单元进行镜像，镜像面为以（0，0，0）为原点，;以（1，0，0）为法向向量。plocreaqqq;创建名为qqqploaddsurfred;表面颜色设置为红色ploshow;显示图形56主要内容FLAC3D软件简介1、基坑开挖2、浅基础的稳定性分析3、网格的合并联结4、界面的生成5、隧道的生成6、模型材料问题7、Interface合并（联结）问题8、初始条件问题9、破坏问题分析10、综合实例-煤巷应力分析57Case-6：模型材料问题genzonecylp0000p1100p2020p3001size454;建立圆柱的网格单元，大小为4×5×4genzonereflectnorm1,0,0;对单元进行镜像，法线为（1，0，0），原点为（0，0，0）genzonereflectnorm0,0,1;对以上所有单元再进行镜像，法线为（0，0，1），原点为（0，0，0）plotCreateTrenchAddsurfaceyellowAddaxesblackShow;显示模型格网pausemodelmohr;采用摩尔-库仑模型propbulk1.19e10shear1.1e10;为模型赋予参数，体积模量为1.19×1010，剪切模量为1.1×1010。propcoh2.72e5fric44ten2e5;内聚力为2.72×105，内摩擦角为44°，抗拉强的为2×105。fixxyzrangey-.1.1;边界y=0为固定支座fixxyzrangey1.92.1;边界y=2为固定支座iniyvel1e-7rangey-.1.1;在y=0的边界上，初始y方向上的速度为1×10−7iniyvel-1e-7rangey1.92.1;在y=2的边界上，初始y方向上的速度为-1×10−7;inipp1e5histgpydisp0,0,0;监测网格点（0，0，0）点的y方向上的位移histzonesyy0,1,0;监测单元（0，1，0）点的y方向上的应力histzonesyy1,1,0;监测单元（1，1，0）点的y方向上的应力step3000;运行3000步58Case-6：模型材料问题监测网格点（0，0，0）点的y方向上的位移59Case-6：模型材料问题监测单元（0，1，0）点的y方向上的应力60Case-6：模型材料问题监测单元（1，1，0）点的y方向上的应力61主要内容FLAC3D软件简介1、基坑开挖2、浅基础的稳定性分析3、网格的合并联结4、界面的生成5、隧道的生成6、模型材料问题7、Interface合并（联结）问题8、初始条件问题9、破坏问题分析10、综合实例-煤巷应力分析62Case-7：Interface合并（联结）问题new;新建;gluedinterfacegrid;粘合交互网格genzonebricksize444p00,0,0p14,0,0p20,4,0p30,0,2;建立块体网格单元，大小为4×4×4，其具体位置由p0,p1,p2,p3的坐标决定。genzonebricksize884p00,0,3p14,0,3p20,4,3p30,0,5;建立块体网格单元，大小为8×8×4，其具体位置由p0,p1,p2,p3的坐标决定。callplot.dat;调用plot.dat文件pause;暂停inter1facerangez2.9,3.1;在z=3的面上建立交互单元inter1propkn300e9ks300e9tens1e10SBRATIO=1;赋予交界面参数，法向刚度为3.0×1011，剪切刚度为3.0×1011，;抗拉强度为1.0×1010，网格的比率为1inizadd-1.0rangez2.9,5.1;z方向上在2.9到5.1的范围内，所有z的坐标都减1modelelas;材料模型为弹性模型propbulk8e9shear5e9;材料的参数为，体积模量为8×109，剪切模量为5×109fixzrangez-.1.1;固定边界z=0，为滚动支座fixxrangex-.1.1;固定边界x=0，为滚动支座fixxrangex3.94.1;固定边界x=4，为滚动支座fixyrangey-.1.1;固定边界y=0，为滚动支座fixyrangey3.94.1;固定边界y=4，为滚动支座applyszz-1e6rangez3.94.1x0,2y0,2;在z=4的平面中，x（0，2），y（0，2）的范围内施加z方向的应力为?1×106histunbal;监测不平衡力solve;计算，一般计算到ratio=1×10?5为止saveinter.sav;将文件保存为inter.savcallplot.dat;调用plot.dat文件pause;暂停63Case-7：Interface合并（联结）问题64Case-7：Interface合并（联结）问题new;新建文档;attachedgrid;连接网格genzonebricksize444p00,0,0p14,0,0p20,4,0p30,0,2;建立块体的网格单元，其大小为4×4×4，其空间位置由p0,p1,p2,p3决定genzonebricksize884p00,0,2p14,0,2p20,4,2p30,0,4;建立块体的网格单元，其大小为8×8×4，其空间位置由p0,p1,p2,p3决定attachfacerangez1.92.1;在z=2的平面上连接交界面modelelas;材料模型为弹性模型propbulk8e9shear5e9;体积模量为8×109，剪切模量为5×109fixzrangez-.1.1;固定边界z=0，为滚动支座fixxrangex-.1.1;固定边界x=0，为滚动支座fixxrangex3.94.1;固定边界x=4，为滚动支座fixyrangey-.1.1;固定边界y=0，为滚动支座fixyrangey3.94.1;固定边界y=4，为滚动支座applyszz-1e6rangez3.94.1x0,2y0,2;在z=4的平面中，x（0，2），y（0，2）的范围内施加z方向的应力为?1×106histunbal;监测不平衡力Solve;计算，一般计算到ratio=1×10?5为止saveatt.sav;将文件保存为att.savcallplot.dat;调用plot.dat文件pause;;暂停65Case-7：Interface合并（联结）问题new;新建文本;onegrid;一个网格genzonebricksize888p00,0,0p14,0,0p20,4,0p30,0,4;建立块体的网格单元，其大小为8×8×8，其空间位置由p0,p1,p2,p3决定modelelas;材料模型为弹性模型propbulk8e9shear5e9;材料的体积模量为8×109，剪切模量为5×109fixzrangez-.1.1;固定边界z=0，为滚动支座fixxrangex-.1.1;固定边界x=0，为滚动支座fixxrangex3.94.1;固定边界x=4，为滚动支座fixyrangey-.1.1;固定边界y=0，为滚动支座fixyrangey3.94.1;固定边界y=4，为滚动支座applyszz-1e6rangez3.94.1x0,2y0,2;在z=4的平面中，x（0，2），y（0，2）的范围内施加z方向的应力为?1×106histunbal;监测不平衡力solve;计算，一般计算到ratio=1×10?6为止savenoatt.sav;将文件保存为noatt.savcallplot.dat;调用plot.dat文件66主要内容FLAC3D软件简介1、基坑开挖2、浅基础的稳定性分析3、网格的合并联结4、界面的生成5、隧道的生成6、模型材料问题7、Interface合并（联结）问题8、初始条件问题9、破坏问题分析10、综合实例-煤巷应力分析67Case-8：初始条件问题APPLY命令可以在其后面选用关键词gradient，这时，我们就可以在一个指定的范围内线性地改变应力或应变。在关键词gradient后面的参数为gx，gy，gz，它们可以用来说明应力或者应变在x，y，z上的变化趋势。应力或者应变线性改变的距离，是以全局坐标系下的原点（0，0，0）为准的。例如：对于这个公式的理解，我们最好用下面的这得例子加以解释：applysxx-10e6gradient0,0,1e5rangez-100,0特别指出，我们可以应用此命令来模拟由重力引起的，随深度而有所增加的应力变化。我们在说明INITIAL命令时，确定应力场是随坐标而产生坡度变化的，并且重力加速度的值已有SETgravity命令给出。68Case-8：初始条件问题genzonebricksize8810ratio1.211;产生块体的网格单元，其到小为8×8×10，在x，y，z方向上的放大比率依;次为1.2，1，1modelmohr;材料模型为摩尔-库仑模型inidens2000;所有材料的初始密度为2000propbulk2e8shear1e8;材料的体积模量为2×108，剪切模量为1×108。propfric30;材料的内摩擦角为30°。fixxrangex-.1.1;固定边界x=0，为滚动支座fixxrangex7.98.1;固定边界x=8，为滚动支座fixyrangey-.1.1;固定边界y=0，为滚动支座fixyrangey7.98.1;固定边界y=8，为滚动支座fixzrangez-.1.1;固定边界z=0，为滚动支座setgrav10;设置重力加速度为10solve;运算saveini1.sav;形成sav文件，保存为ini1.savplotcreateszz_contour;创建文件名为szz_contour（z方向的应力云图）plotsetcent445;视图的中心坐标为（4，4，5）plotsetrot20030;视图的旋转角度为（20，0，30）plotsetmag1.0;视图的放大倍数为1.0plotaddcontszzoutonshadeon;打开szz的等值线图plotaddaxes;打开坐标轴线plotshow;显示窗口pause;暂停69Case-8：初始条件问题70Case-8：初始条件问题new;新建文本genzonebricksize8810ratio1.211;产生块体的网格单元，其到小为8×8×10，在x，y，;z方向上的放大比率依次为1.2，1，1modelmohr;材料模型为摩尔-库仑模型inidens2000;所有材料的初始密度为2000propbulk2e8shear1e8;材