FLAC,FLAC3D基础与应用-陈育民

1、1河海大学研究生课程岩土数值分析 2002年本科毕业于中国矿业大学 2007年博士毕业于河海大学 2013.7-2014.7美国普渡大学访问学者 2015.22015.5日本东京大学访问学者 研究方向：土动力学与岩土地震工程土木工程防灾减灾岩土工程数值分析 联系方式：32005-11-29 河海土木院研究生会组织2006-10-13 同济大学土木工程学院2006-10-26 河海大学金水节2007-04-15 东南大学交通学院2007-07-18 同济大学土木工程学院2007-11-03 河海大学岩土所组织FLAC学术沙龙2007-11-29 河南工业大学2008-11-15 河海大学河海金水

2、节培训2010-11-10 河海大学校庆报告2011-06-18 河海大学举办ITASCA技术与应用专题（南京）研讨会2011-10-16 河南理工大学2011-11-03 南京工业大学交通学院2011-11-24 河海大学土木与交通学院研究生会2011-06-18_ITASCA技术与应用专题（南京）研讨会2012-08-31_解放军理工大学FLAC讲座4 什么是FLAC？ 为什么要用FLAC？ FLAC能做什么？ FLAC为何这么流行？ 怎么学FLAC？ 研究生课程岩土数值分析上课学生调查67 Fast Lagrangian Analysis of Continua89检索期刊：岩土工程学报

3、岩土力学岩石力学与工程学报关键词： 软件软件论文数量论文数量FLAC650ABAQUS278ANSYS206PLAXIS61ADINA41PFC171UDEC73Updated on June 2, 201510检索期刊：Journal of Geotechnical and Geoenvironmental EngineeringGeotechniqueCanadian geotechnical JournalSoils and foundations 关键词：Updated on June 2, 2015 软件软件论文数量论文数量FLAC31ABAQUS23ANSYS5PLAXIS14AD

4、INA0PFC2UDEC911检索期刊：GEO*SOILS*关键词：Updated on June 2, 2015 软件软件论文数量论文数量FLAC88ABAQUS50ANSYS27PLAXIS74ADINA3PFC12UDEC16 岩土工程中的绝大多数问题土力学、岩石力学、防灾减灾、隧道、地下空间等 采矿工程中的大部分问题 水工结构中的部分问题 结构工程国际通用的岩土工程国际通用的岩土工程专业分析程序专业分析程序12 Charles Fairhurst美国工程院、瑞典皇家工程院院士，国际岩石力学学科和岩石力学学会创始人之一，历任国际岩石力学学会主席和副主席，国际岩石力学学会MULLER奖、美

5、国岩石力学学会终生成就奖获得者。 Peter Cundall美国工程院、英国皇家工程院院士，国际资深计算岩石力学学家。13 第一讲：基本介绍、静力分析、前后处理 第二讲：接触面、FISH语言、流固耦合分析 第三讲：动力分析、自定义本构、结构单元 第四讲：FLAC（2D）基本介绍与应用实例 讨论14 了解FLAC的基本概念 了解FLA的适用范围和缺陷 了解深入学习FLAC的方法 会用FLAC分析简单的路堤填筑问题（考题）复杂的网格往往会把核心问题覆盖掉，经复杂的网格往往会把核心问题覆盖掉，经常会造成计算结果无法解释。常会造成计算结果无法解释。现实中的桩现实中的桩虚拟的桩虚拟的桩冰碛土体结构模拟结

6、果冰碛土体结构模拟结果台阶坡面上的砾石产出状态台阶坡面上的砾石产出状态 数值分析要对实际工程进行大量的、细心的简化，从效率和结果两方面保证数模分析的成果 结果的判断，需要深入扎实的理论功底及“丰富”的工程经验 数值分析的作用从“锦上添花”到“雪中送炭”，在于各位的努力，指日可待FLAC3D基本介绍、静力分析、前后处理24 Fast Lagrangian Analysis of Continua 美国Itasca咨询公司开发2D程序(1986) 1990年代初引入中国 有限差分法(FDM) 3D版本：DOS版2.0 2.1 3.0 3.14.05.0 2D版本：DOS版4.05.06.07.02

7、5 大应变、小应变计算模式. 丰富的本构模型、提供自定义的本构模型功能 接触面可以模拟不同材料的接触 流固耦合实现土体的固结与渗流 拥有各种功能的结构单元类型，模拟土与结构的相互作用 强大的动力分析功能. 流变分析，拥有粘弹性模型和粘塑性模型 热力学分析.26Shear strainrate contours27- 与FLAC类似，是FLAC的三维版本- 与 FLAC拥有相同的优点upstreamdownstream 内置材料模型连续介质非线性，大应变模拟显式解题方案，为不稳定物理过程提供稳定解界面或滑动面用来模拟可产生滑动或分离的离散面，从而模拟断层，节理或摩擦边界内置材料模型丰富：零模型,

8、 三个弹性模型 (各向同性，横观各向同性和正交各向异性), 八个朔性模型 (德鲁克-布拉格, 摩尔-库伦, 应变硬化/软化，单一节理，双线性应变硬化/软化单一节理, 双屈服，修正剑桥粘土，霍克-布朗)28隧道工程 可选模块包括: 热力学,热-力学耦合,热-流体-力学耦合包括热传导和对流; 粘弹,粘朔性(蠕变)材料模型;动力学分析,并可以模拟静边界和自由域 使用C+定义自己的模型29核废料储存中的热力学研究问题核废料储存中的热力学研究问题 FLAC/FLAC3D利用有限差分，显示方案，动态松弛方法模拟连续体的非线性力学行为：即使对准静态问题，程序仍然求解完整的动力学方程。这种方法的好处在于可以为

9、物理非稳定过程例如塌方提供稳定解；在 “松弛”方法中，使用阻尼来吸收动能以模拟系统的“静态”反应。 这种方法可以用比其它方案如解矩阵法更为真实有效地模拟塌方问题。30 源自流体力学中的拉格朗日法跟踪流体质点的运动状态跟踪固体力学中结点，按时步用Lagrangian法研究网格节点的运动 节点和单元随材料移动，边界和接触面与单元的边缘一致 固体力学大变形理论31法国数学家、物理学家拉格朗日 32+/2=每个每个 为常应力为常应力/应变应变:体积应变由整个四边形算出体积应变由整个四边形算出 . 应变偏量则有两个三角形应变偏量则有两个三角形 和和 分别算出分别算出(混合离散混合离散 过程过程)解题过程

10、中网格坐标按照解题过程中网格坐标按照“拉格朗日方式更新拉格朗日方式更新” (网格随材料移动网格随材料移动), 且为且为显式显式 (一个时步内局部变化不会影响邻域一个时步内局部变化不会影响邻域)FLAC3D混和离散混和离散33+/2=34n结构域离散为可由四面体单元组合形成的五面体或六面体等结构域离散为可由四面体单元组合形成的五面体或六面体等单元单元; ;n以以 为基本单元为基本单元( (常应力、常应变常应力、常应变););n体应变的计算：体应变的计算： ; ;n偏应变的计算：偏应变的计算： . .动态松弛法在动态松弛法中，网格点根据牛顿运动定律运动. 网格点的速度与该点的不平衡力呈正比. 这种

11、求解方法所决定的一系列位移将把系统带入平衡状态，或表明破坏模式.在动态松弛法中有两个因素很重要:时步的选择阻尼效应 35显式解与隐式解的比较显式解与隐式解的比较36显式显式,逐时推进逐时推进隐式隐式, 静态静态1.无需进行反复迭代来实现非线无需进行反复迭代来实现非线性本构关系性本构关系 .2. 类似问题求解时间呈类似问题求解时间呈 N3/2 规规律增长律增长3. 物理非稳定性不会引起数值物理非稳定性不会引起数值不稳定性不稳定性.4. 因为无需储存矩阵，用较小因为无需储存矩阵，用较小内存即可模拟大尺度问题内存即可模拟大尺度问题.5.对大位移、大应变问题同样适对大位移、大应变问题同样适合，无需额外

12、的计算合，无需额外的计算 .1.需进行反复迭代来实现非线性需进行反复迭代来实现非线性本构关系本构关系 2.类似问题求解时间呈类似问题求解时间呈 N2 甚至甚至 N3规律增长规律增长.3.难以模拟物理非稳定性问题难以模拟物理非稳定性问题.4.需存储刚度矩阵，需克服相关需存储刚度矩阵，需克服相关的带宽问题，需要的内存较大的带宽问题，需要的内存较大 .5.对大位移、大应变问题需进行对大位移、大应变问题需进行大量的计算大量的计算 . 使用Intel Fortran compiler拥有更快的计算速度 自动网格重画功能，解决 bad-geometry 问题. 新的模拟颗粒土材料的硬化模型 更新的通用网格

13、生成工具37多处理器的并行计算功能 新结构单元类型 “Embedded Liner” 提供两个方向的接触作用，可以很好地模拟挡土墙对四面体单元采用新的混合离散方法 “Nodal Mixed Discretization” 提供塑性问题更精确的解答64位程序包含命令手册、FISH手册和应用实例的帮助38 模拟颗粒状材料的硬化模型 自动网格重画功能，解决 bad-geometry 问题. 改进的interface 更快的渗流计算 更新的动力计算功能3940F(t)duFm amdtijiijdugdtx, ,u u u m牛顿运动定律对于连续体在静力平衡条件下，加速度项为0，方程变为平衡方程41G

14、 = mgS = 1/2gt2 = 20m命令流：命令流：config dyngen zon bri size 1 1 1ini x mul 0.1 y m 0.1 z m 0.1model elasprop bulk 3e8 shear 1e8ini dens 1000set grav 0 0 -10solve age 24243节点gridpoint：节点zone：单元boundary：边界44平衡方程(动量方程)应力应变关系(本构模型)Gauss定律单元积分应变率速度节点力新的应力对所有的网格节点对所有单元 开挖模型null 3个弹性模型各向同性弹性横观各向同性弹性正交各向同性弹性 8个

15、弹塑性模型Drucker-Prager模型、Morh-Coulomb模型、应变硬化/软化模型、遍布节理模型、双线性应变硬化/软化遍布节理模型、修正剑桥模型和胡克布朗模型4546gen zon bri size 3 3 3 ;建立网格(前处理)model elas ;材料参数prop bulk 3e6 shear 1e6ini dens 2000 ;初始条件fix z ran z -.1 .1 ;边界条件fix x ran x -.1 .1fix x ran x 2.9 3.1fix y ran y -.1 .1fix y ran y 2.9 3.1set grav 0 0 -10solve ;

16、求解app nstr -10e4 ran z 3 x 1 2 y 1 2solveplo con zd ;后处理切片功能47RUN FLAC3D48建立网格建立网格初始条件初始条件边界条件边界条件初始应力平衡初始应力平衡外荷载外荷载求解求解前处理前处理后处理后处理n保存文件保存文件 (*.sav) 含有所有状态变量和用户定义条件的二进制文件 n数据文件数据文件 (*.dat) 数据文件由用户创建的一种ASC格式的文件，它包括一系列的用于描述所分析问题的FLAC3D命令 nFISH文件文件(*.fis) FISH程序文件nFLAC3D文件文件(*.flac3d) FLAC3D的网格信息文件n历史

17、记录文件历史记录文件 (*.his) 记录输入输出历史值的文件 n图形文件图形文件 图形文件(各种标准格式)n电影文件电影文件 (*.dcx) AVI或PCX图像文件，这些图像文件可以当作电影放映49 为什么要单独列出？分析过程中出现的很多问题都与初始应力是否合理有关手册中的例子五花八门是所有后续分析的基础！ 生成方法弹性求解更改强度参数的弹塑性求解设置初始应力的弹塑性求解存在水压力的初始应力生成水下建筑的初始应力生成50模型尺寸单元数量密度KGu112 (m3)112200030MPa10MPa0.3551gen zon bri size 1 1 2m elasprop bulk 3e7 s

18、hear 1e7fix z ran z 0fix x ran x 0fix x ran x 1fix y ran y 0fix y ran y 1ini dens 2000set grav 0 0 -10solveStep = 162 z = -40e3 x = -21.54e3模型尺寸单元数量密度KGcfu112 (m3)112200030MPa10MPa10kPa150.3552gen zon bri size 1 1 2model mohrprop bulk 3e7 shear 1e7 c 1e10 f 15 ten 1e10fix z ran z 0fix x ran x 0fix x

19、 ran x 1fix y ran y 0fix y ran y 1ini dens 2000set grav 0 0 -10solveprop bulk 3e7 shear 1e7 c 10e3 f 15 ten 0solveStep = 163 z = -40e3 x = -21.54e3Or: solve elastic模型尺寸单元数量密度KGcfu112 (m3)112200030MPa10MPa10kPa150.3553gen zon bri size 1 1 2model mohrprop bulk 3e7 shear 1e7 c 10e3 f 15 ten 0fix z ran

20、z 0fix x ran x 0fix x ran x 1fix y ran y 0fix y ran y 1ini dens 2000ini szz -40e3 grad 0 0 20e3 ran z 0 2ini syy -20e3 grad 0 0 10e3 ran z 0 2ini sxx -20e3 grad 0 0 10e3 ran z 0 2set grav 0 0 -10solveStep = 0 z = -40e3 x = -20e3模型尺寸单元数量饱和密度KGcfu水位线孔隙率112(m3)112200030MPa10MPa10kPa150.351m0.554gen zon

21、e 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 1ini dens 1500 ran z 1 2ini szz -35e3 grad 0 0 20e3 ran z 0 1ini syy -22.5e3 grad 0 0 15e3 ran z 0 1ini sxx -22.5e3 grad 0 0 15e3 ran z 0 1i

22、ni szz -30e3 grad 0 0 15e3 ran z 1 2ini syy -15e3 grad 0 0 7.5e3 ran z 1 2ini sxx -15e3 grad 0 0 7.5e3 ran z 1 2ini pp 10e3 grad 0 0 -10e3 ran z 0 1set grav 0 0 -10solveStep = 0 z = -35e3 x = -22.5e3 d = s ns f模型尺寸单元数量饱和密度KGcfu水位线孔隙率112(m3)112200030MPa10MPa10kPa150.351m0.555config fluidgen zon bri s

23、ize 1 1 2model elasprop bu 3e7 sh 1e7ini dens 1500model fl_isoini fdens=1000 fmod 0prop por 0.5set grav 0 0 -10water dens 1500water table face 0 0 1, 0 1 1, 1 1 1, 1 0 1 ini pp 10e3 grad 0 0 -10e3 ran z 0 1ini szz -30e3 grad 0 0 15e3 ran z 1 2.set fluid offsolveStep = 142 z = -40e3 x = -24e3单元数较少产生的

24、误差单元数较少产生的误差模型尺寸单元数量饱和密度KGcfu水位线112(m3)112200030MPa10MPa10kPa150.353m56gen zon bri size 1 1 2model mprop bulk 3e7 shear 1e7 c 10e10 f 15 ten 1e10fix 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 2ini syy -40e3 grad 0 0 15e3 ra

25、n z 0 2ini sxx -40e3 grad 0 0 15e3 ran z 0 2ini pp 30e3 grad 0 0 -10e3 ran z 0 2app nstress -10e3 ran z 2set grav 0 0 -10solveStep = 0 z = -50e3 x = -40e3FLAC3D 3.00Itasca Consulting Group, Inc.Minneapolis, MN USAStep 1 Model Perspective17:30:28 Wed Nov 02 2011Center: X: 5.000e-001 Y: 5.000e-001 Z:

26、1.000e+000Rotation: X: 0.000 Y: 0.000 Z: 0.000Dist: 6.030e+000Mag.: 1Ang.: 22.500Contour of SXX Magfac = 0.000e+000 Gradient Calculation-4.0000e+004 to -3.7500e+004-3.7500e+004 to -3.5000e+004-3.5000e+004 to -3.2500e+004-3.2500e+004 to -3.0000e+004-3.0000e+004 to -2.7500e+004-2.7500e+004 to -2.5000e

27、+004-2.5000e+004 to -2.2500e+004-2.2500e+004 to -2.0000e+004-2.0000e+004 to -1.7500e+004-1.7500e+004 to -1.5000e+004-1.5000e+004 to -1.2500e+004-1.2500e+004 to -1.0000e+004-1.0000e+004 to -1.0000e+004 Interval = 2.5e+00357 基本前后处理命令操作菜单操作 dd & dip attach & merge 外界模型的导入 复杂模型的网格检查58 命令驱动(推荐)程序

28、控制图形界面接口计算模型输出指定本构模型及参数指定初始条件及边界条件，指定结构单元指定接触面指定自定义变量及函数(FISH)求解过程的变量跟踪进行求解模型输出5960命令栏6162FLAC3D网格生成的关键特征：网格生成的关键特征：FLAC3D 是命令驱动.使用FLAC3D内置基元进行形状组合可形成复杂网格.用户自定义FISH 函数可以用来修改基元网格以创建更为复杂的网格.第三方软件导入.63块体块体退化块体退化块体楔体楔体金字塔金字塔四面体四面体柱体柱体径向块体径向块体径向隧道径向隧道径向柱体径向柱体柱状壳体柱状壳体柱状交叉柱状交叉隧道交叉隧道交叉6465块体块体径向柱体径向柱体柱状交叉柱状

29、交叉66676869707172TecplotTecplot中云图的效果与中云图的效果与flac3dflac3d的的shade onshade on的效果差的效果差不多不多, ,但是其出三维等值线的功能是卓越的但是其出三维等值线的功能是卓越的. .7374TecplotTecplot的的sliceslice切片切片功能功能. .和和CADCAD一样一样, ,可以任意可以任意切剖面出切剖面出图图. .最大优最大优点是可以点是可以几个剖面几个剖面同时出图同时出图. .75整体和截面的整体和截面的数据以单元形数据以单元形式存在列表中式存在列表中, ,可自由选择和可自由选择和组合多个单元组合多个单元出

30、图出图. . Simwe.Com上有近百的讨论贴 地质上的概念，倾向和倾角 建议用ori和norm代替 后处理的切片功能plo set plane ori (*,*,*) norm (*,*,*)plo con zdis plane概念清晰，使用快捷 其他应用reflect网格建立水面76xyzorinorm attach用于连接次节点 (sub-grid)允许网格存在一定随意性要求成比例(1:2, 1:3等)可用于模型的检查慎用 merge用于节点的合并外来模型导入的精度差异tolerance的设置77sub-grid 4:25:2tol121merg 复杂网格的生成难度大 接口编写不同软件

31、之间的精度差异Group的定义采用.flac3d的文件格式 .flac3d文件的格式G 1 1.0e+00 1.0e+00 1.0e+00Z B8 1 2 3 4 5 6 7 8ZGROUP Soil1 2 3 impgrid & expgrid仅限于网格783w zones 3sec FLAC3D生成的复杂网格attach face无接触面时检查整体模型是否存在无接触面时检查整体模型是否存在sub-grid有接触面时给定范围进行检查有接触面时给定范围进行检查 其他软件导入的复杂模型网格划分的检查gen merge弹性模型model elastic求解“独立独立”节点节点“畸形畸形”单

32、元单元79 FLAC本身的Check功能十分有限错误提示很少十分开放的工作平台 (自由落体) 检查的基本步骤网格检查(如前所述)边界条件检查速度约束条件速度约束条件 plo gpfix red sk荷载条件荷载条件 plo fap red sk模型检查模型赋值模型赋值plo block model参数赋值参数赋值plo block prop *初始应力检查(如后所述)80 初始应力的计算时间不会“非常长” 经常检查模型的响应plo con szz (syy, sxx)应力场plo con zdis (ydis, xdis)位移场plo blo sta屈服状态plo gpfix red sk速度

33、约束条件plo fap red sk体力plo hist (unbal)不平衡力plo interface nstress (sstress)接触面单元81FLAC3D接触面、FISH语言、流固耦合分析 应用范围 原理 建模方法推荐方法复杂内部接触面的设置方法 参数选择 单桩承载力分析 挡土墙的接触面设置思路：未知问题的分析方法 网格不连续 岩体介质中的解理、断层、岩层面 地基与土体的接触 箱、槽及其内充填物的接触 空间中无变形的固定“障碍” 三角形单元 (无厚度!) 参数较多 三种工作模式粘结界面粘接滑移库伦滑动 关键要形成同一位置的两个节点(面) “移来移去”(推荐)建两个分开的模型建立接

34、触单元通过INI * add使模型接触注意dist的含义NO merge, NO attach!接触面dist1234 “导来导去”利用expgrid, impgrid命令进行网格导出与导入配合DELETE命令适于内部接触面的建立，或其他前处理工具建立的网格 GEN separate INTERFACE wrap指定正确的group 虚构的为了合并节点而设置的接触面Kn=ks=10*max(K+4/3G)/Dzmin 真实的刚性接触面如料仓下料c,D,Tension重要，kn,ks不重要 真实的柔性接触面断层；水力劈裂材料试验得到参数对于kn,ks：岩石断层10100MPa/m(粘土); 10

35、0GPa(岩石)反分析方法：通过断层中岩石的变形与原岩的变形 软土地基bulk 1.6878E6 shear 3.6167E5 coh 15E3 fric 12dens 1.73E3 桩体bulk 5e9shear 3.75e9dens 2.5e30.5m 8m 10m 20m 施加桩顶荷载施加桩顶荷载计算结果计算结果7.50E+03fricgk3coh/0.70.7fric10g10k23.00E+0420100e100e1cohfrickskn水平因素1233(9)3123(8)2313(7)2132(6)1322(5)3212(4)3331(3)2221(2)1111(1)cohfric

36、kskn 水平因素工况Ks取取1 Ks取取2 Ks取取3 96.9%1233(9)22.4%3123(8)54.1%2313(7)96.9%2132(6)7.1%1322(5)37.8%3212(4)96.9%3331(3)42.9%2221(2)37.8%1111(1)Ra差异度差异度cohfrickskn水平因素水平因素2222111122101122最优方最优方案案17.3%6.8%72.8%11.9%极差极差52.4%52.7%96.9%57.8%k364.6%59.2%24.1%47.3%k247.3%52.4%43.2%59.2%k1cohfrickskn水平因水平因素素单桩分析单

37、桩分析简单网格简单网格接触面参数接触面参数多次试算多次试算理论、实测理论、实测加密网格加密网格接触参数接触参数理想结果理想结果群桩分析群桩分析Pile结构单元结构单元单元参数单元参数理想结果理想结果Pile结构单元结构单元单元参数单元参数理想结果理想结果对于未知问题的分析思路3个独立的接触面3个同ID的接触面2个独立的接触面2个独立的接触面并进行底部merge挡墙挡墙 土体土体 123WallSoilZ = 5.16cmZ = 1.14cmZ = 不能初始平衡不能初始平衡 Z = 0.35cm？X = 1.95cmX = 1.91cmX = 1.47cm？ 前两种方法的差别实质3个ID的独立接

38、触面在相同位置产生互不影响的两个节点共同ID的接触面在相同位置自动设置为1个节点 最终的结论需要您自己去判断！3 interfaces, 2 IDs 1 interface, 1 ID 软件自带的编程语言是否一定要学？视情况而定，需要时查询FISH变量即可语法简单xxxend_xxx注意事项与FLAC本身的关键字冲突保留字不可缩写变量可不定义，因此注意检查程序print fishtable, extra等命令使用def abc abc = 1 + 2 * 3 abcd = 1.0 / 2.0endabcprint fish数据格式 都可以在FISH函数中进行赋值，赋值操作与常规的编程语言类似，

39、按照运算符的优先级先后顺序来执行。 函数和变量的赋值遵守数据类型的规则，即整型的计算结果为整型，浮点型的计算结果为浮点型，因此读者在进行除法运算、开方运算时都需要将数据类型设置为浮点型，数字尽量使用小数点以保证运算正确。 变量和函数名的命名规则是不能以数字开头，不能含有中文，并且不能包含如下的字符。. , * / + - = # ( ) ; 变量和函数名不能与FLAC3D、FISH的保留字相冲突，不要采用过于简单的单词，比如a，hist等，这些都与保留字相冲突。 即使程序中存在与保留字相冲突的变量，FLAC3D也不会提供任何提示，所以提醒读者在编制FISH程序时尽量使用较长的、复杂的变量和函数

40、名。对变量进行赋值时，不能使用当前函数的函数名放在“=”的右边，比如采用下面的定义abcd = abc + 1.0在FISH程序执行时会提出错误，因为这样会形成递归调用，这种调用方式在FISH程序中是不允许的。变量和函数的作用是全局的，在命令中的任何地方修改变量的值都会立即生效，因此在实际应用中尽量避免不同的函数中含有相同的变量，因为这样可能会造成赋值错误，并难以检查。在FLAC3D中可以用如下的命令来引用FISH函数和变量PRINT用于查看函数和变量的数值；HISTORY命令可以对函数和变量的数值进行记录；SET命令用于变量的赋值。 选择语句CASEOF 表达式默认语句CASE n1表达式的

41、值为n1时的语句CASE n2表达式的值为n2时的语句ENDCASE 条件语句IF 条件表达式 THENELSEENDIFFISH中条件运算符没有“并”、“或”、“否”这样的符号 if aa 1.0 if aa 2.0 执行语句 endifendif表达“1aa maxdisp_value maxdisp_value = disp_gp maxdisp_gpid = gp_id(p_gp) endif p_gp = gp_next(p_gp) endloop endfind_max_dispprint maxdisp_value maxdisp_gpid第一步第一步第二步第二步def abce

42、ndabcdef abc p_gp = gp_head loop while p_gp # null p_gp = gp_next(p_gp) endloopendabc第三步第三步第四步第四步def abc p_gp = gp_head loop while p_gp # null command endcommand p_gp = gp_next(p_gp) endloopendabcdef abc p_gp = gp_head loop while p_gp # null command app nstress endcommand p_gp = gp_next(p_gp) endloo

43、pendabcTunnel.txt两个圆形隧道的连接部分两个圆形隧道的连接部分 变直径的隧道部分变直径的隧道部分 主要采用PRINT fish命令 查看变量的赋值是否合理，主要检查值为0的函数和变量，因为FISH程序中一般定义的变量都有实际的意义，输出为0的变量很可能是与保留字相冲突的变量（如a就是apply的保留字） 或者由于编写笔误产生的变量（如数字0与大写字母O）。 很难！ 基本功能 理论框架 计算模式 渗流边界条件，初始条件 单渗流计算及渗流耦合计算 渗流各向同性、各向异性不同的渗流模型和属性流体压力，涌入量，渗漏量和不渗水边界抽水井、点源、体积源饱和渗流可采用显式差分法、隐式差分法

44、非饱和渗流采用显式差分法渗流-固体-热的耦合流体和固体的耦合程度依赖于土体颗粒(骨架)的压缩程度，用Biot系数表示颗粒的可压缩程度。 循环荷载引起的动水压力变化和土体液化 准静态Biot理论 多孔介质中遵循Darcy定律的单相渗流 描述多孔介质中流体渗流的变量孔隙水压力，饱和度，特定排水向量的三个分量 质量守恒定律 达西定律 本构定律考虑流体响应孔隙水压力改变，饱和度改变，体积应变改变和温度改变 无渗流模式不设置CONFIG Fluid孔压不改变设置孔压分布INITIAL ppWATER table WATER density WATER table face SET gravity手动设置

45、干湿密度F 渗流模式 设置CONFIG fluid 设置土体干密度 渗流模型 MODEL fl_isotropic MODEL fl_anisotropic MODEL fl_null SET fluid off set WATER bulk = 0 默认的边界条件是不透水边界 孔隙压力自由(不透水边界 ) 固定孔隙水压力(透水边界) 如：井 孔隙压力，孔隙率，饱和度和流体属性的初始分布可以用INITIAL命令或者PROPERTY命令定义。 时间比例 完全耦合分析方法 孔压固定分析(有效应力分析) 单渗流得到孔压分布 无渗流计算孔压的力学响应 流-固耦合计算 用途：排水沟；抽水井；耦合计算 计

46、算步骤CONFIG fluid SET mech offSET fluid implicit on/offMODEL fl_; PROP STEP; SOLVE age; SET fluid ratioSET fluid off mech onPROP biot_c 0 (or INI fmod 0) 不排水短期响应 两种分析方法：干法和湿法干法：Ku=K+a2M两种破坏形式两种破坏形式 WATER或INI获得常孔压，不排水的c, (孔压改变较小) =0,c=cu (MK+4/3G)湿法：耦合体系的短期行为使用排水的使用排水的K, c,若若SET fluid off, Biot_mod(fmo

47、d)真实真实力学过程的特征时间力学过程的特征时间 流体扩散过程的特征时间流体扩散过程的特征时间 时间比例短期行为 (不排水) ts(分析时间分析时间)tc施加扰动的属性流体扰动：渗流可不与力学过程耦合力学扰动：耦合等级取决于流固刚度比流固刚度比4/3mccutLKG24/3kMRKG2fccLtc CONFIG fluid; M(Kf); K(渗透系数) 真实，则FLAC3D默认耦合计算pevevp 预估流/力特征时间 耦合计算前先达到一个平衡状态SET fluid on mech off; SET fluid off mech on; STEPSET mech force; SET mech

48、 substep n auto; SET fluid substep m (=1)STEP:渗流步足够小 手动调整的STEP求解SET fluid on mech offSTEPSET fluid off mech onSTEP 主从进程的SOLVE求解SET mech forceSET mech substep n auto (从进程)SET mech substep m(主进程)SOLVE age 自动STEP求解STEP 时间比例(ts, tc)稳态不排水状态相当 扰动类型力学扰动孔压扰动 流固刚度比Rk是否1 完全耦合模式时间比例相当；力学扰动newconfig fluid set f

49、luid offgen zon brick p0 0 0 -10 size 20 1 10gen zon brick p0 5 0 0 p1 15 0 0 p2 5 1 0 p3 9 0 5 p4 15 1 0 p5 9 1 5 p6 11 0 5 p7 11 1 5 size 10 1 5group soilgroup dam ran x 5 7 z -5 0 group dam ran id 201 a id 211 a id 221 a id 231 a id 241 a group dam ran id 202 a id 212 a id 222 a id 232 a id 242 a

50、m eprop bu 3e7 sh 1e7ini pp 0 grad 0 0 -10e3 ran z 0 -10ini dens 2000model fl_isoprop por 0.5 perm 1e-10 ini fden 1000 ften -1e10ini sat 0.0 ran z 0 5model fl_null ran gro dam;ini pp 0 ran gro damfix z ran z -10fix x ran x 0fix x ran x 20fix yset grav 10solvesave elastic.sav网格模型初始孔压rest elastic.savi

51、ni xd 0 yd 0 zd 0 xv 0 yv 0 zv 0app nstress -40e3 grad 0 0 10e3 ran z 0 4 x 0 9solvesave pressure.sav竖向应力沉降rest pressure.savset fluid on mech offini fmod 2e3 ften 0.0 ran gro soilini xd 0 yd 0 zd 0 xv 0 yv 0 zv 0app pp 40e3 grad 0 0 -10e3 ran z 0 4 x 0 9app pp 0 ran z 0 x 15 20hist id=10 zone pp id

52、215solveRatio = 1Ratio = 1E-510 m3 mK = 500 MPa, G = 300 MPa, c = 10 kPa f = 1540 kPa20 m; - apply load slowly -def ramp ramp = min(1.0,float(step)/200.0)endapply nstress = -40e3 hist ramp range x -.1 3.1 z 9.9 10.1; - fluid flow model -model fl_isoini fmod 2e9; - pore pressure fixed at zero at the

53、surface -fix pp 0 range z 9.9 10.1; - settings -set fl off; - test -step 75011因为本例中没有设置初始应力，这里只进行了750步的求解。FLAC3D 3.00Itasca Consulting Group, Inc.Minneapolis, MN USAStep 750 Model Perspective16:45:13 Wed Apr 09 2008Center: X: 1.000e+001 Y: 5.000e-001 Z: 5.000e+000Rotation: X: 0.000 Y: 0.000 Z: 0.000

54、Dist: 5.580e+001Mag.: 1Ang.: 22.500Plane Origin: X: 0.000e+000 Y: 5.000e-001 Z: 0.000e+000Plane Normal: X: 0.000e+000 Y: 1.000e+000 Z: 0.000e+000Job Title: 荷载引起的超孔隙水压力View Title: Contour of Pore Pressure Plane: on Magfac = 0.000e+000-8.7557e+002 to 0.0000e+000 0.0000e+000 to 5.0000e+003 5.0000e+003

55、to 1.0000e+004 1.0000e+004 to 1.5000e+004 1.5000e+004 to 2.0000e+004 2.0000e+004 to 2.5000e+004 2.5000e+004 to 3.0000e+004 3.0000e+004 to 3.5000e+004 3.5000e+004 to 4.0000e+004 4.0000e+004 to 4.5000e+004 4.5000e+004 to 4.6350e+004 Interval = 5.0e+003SET mech force 设置一个不平衡力的大小，达到这个不平衡力系统认为暂时达到平衡状态；SE

56、T mech substep n auto设置力学进程为从进程，在主进程每执行一步中必须执行n步，当系统达到平衡时也可以少于n步；SET fluid substep m设置流体进程为主进程。讨论讨论1：对收敛准则进行对比分析：对收敛准则进行对比分析 不平衡力（force）：1E3、5E3、1E4、5E4不平衡力比（ratio）：1E-4、1E-3、1E-2、1E-1采用1E-3的收敛准则既可以满足计算流固耦合过程中的计算精度要求，同时又具有较高的计算效率 讨论讨论2：子步数的影响：子步数的影响 set mech sub 100 fluid sub 10set mech sub 10 fluid

57、 sub 10set mech sub 1 fluid sub 1设置合理的子步数也很重要。设置过大，则会导致计算时间大大增加，过小又会造成计算结果的误差。 孔压边界条件tstc长期分析(排水)Rk1骨架很软孔压扰动进行biot_mod调整砂层软土层粘土层PVD2m8m10m 第一讲：基本介绍、静力分析、前后处理 第二讲：接触面、FISH语言、流固耦合分析 第三讲：动力分析、自定义本构、结构单元 第四讲：FLAC（2D）基本介绍与应用实例 讨论139140非常复杂！Said by Prof. Peter Cundall FLAC 可以模拟体系（土，岩石，结构，流体）受到的外部动力荷载（比如地震

58、）或内部动力荷载（比如基础振动、爆炸）。 可以计算塑性引起的永久变形以及孔隙水压力的消散。 土动力学中常用的等效线性方法无法直接处理上述问题。 141动力荷载与边界条件材料响应与阻尼土体液化142 动力输入的类型加速度时程速度时程应力(压力)时程力时程 APPLY INTERIOR (内部)TABLEFISH143 静态(quiet,粘性)边界Lysmer and Kuhlemeyer(1969) 模型边界法向和切向设置独立的阻尼器 性能对于法向p波和s波能很好的吸收对于倾斜入射的波和Rayleigh波也有所吸收，但存在反射人工边界仍应当足够远144 内部振动(如隧道中的列车振动问题)动力荷载

59、直接施加在节点上使用Quiet边界减小人工边界上的反射不需要FF边界 外部荷载的底部边界软土地基上的地震荷载不适合用加速度或速度边界条件使用应力条件t = -2Csvs 地震底部输入的侧向边界扭曲了入射波145quietquietquiet146Cundall et al. (1980)自由场网格与主体网格的耦合粘性阻尼器，自由场网格的不平衡力施加到主体网格边界上设置条件底部水平，重力方向为z向侧面垂直，法向分别为x, y向其他边界条件在APPLY ff之前相当于一个阻尼器 APPLY ff将边界上单元的属性、条件和变量全部转移ff单元上； 设置以后主体网格上的改动将不会被FF边界所响应 可存

60、在任意的本构模型以及流体耦合(仅竖向) FF边界进行小变形计算，主体网格可大变形，FF边界上的变形要相对较小 存在attach的边界将不能设置FF边界 边界上的Interface将不能连续 动力边界设置需在FF边界设置之前1472SsCvt 148模型底部边界fix施加速度或加速度荷载刚性边界Free施加应力时程荷载柔性边界对于软弱的地基不适合施加速度（加速度荷载），而应当施加应力荷载Note that there is a factor of 2 because the input energy divides into a downward- & upward-propagating wave.149连续的非线性，表