FLAC,FLAC3D基础与应用ppt课件.ppt

FLAC/FLAC3D基础与应用1关于教师●2002年本科毕业于中国矿业大学●2007年博士毕业于河海大学●2013.7-2014.7美国普渡大学访问学者●2015.2~2015.5日本东京大学访问学者●研究方向：○土动力学与岩土地震工程○土木工程防灾减灾○岩土工程数值分析●联系方式：ymchenhhu@2关于教材3关于课程●2005-11-29 河海土木院研究生会组织●2006-10-13 同济大学土木工程学院●2006-10-26 河海大学金水节●2007-04-15 东南大学交通学院●2007-07-18 同济大学土木工程学院●2007-11-03 河海大学岩土所组织FLAC学术沙龙●2007-11-29 河南工业大学●2008-11-15 河海大学河海金水节培训●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？5课堂“作业”●研究生课程《岩土数值分析》上课学生调查6课堂“作业”7什么是FLAC？●Fast Lagrangian Analysis of Continua连续介介质拉格朗日算法拉格朗日算法有限差分法有限差分法快速快速8为什么要用FLAC？检索期刊：•岩土工程学报•岩土力学•岩石力学与工程学报•关键词：××× 软件件论文数量文数量FLAC650ABAQUS278ANSYS206PLAXIS61ADINA41PFC171UDEC73Updated on June 2, 20159为什么要用FLAC？检索期刊：•Journal of Geotechnical and Geoenvironmental Engineering•Geotechnique•Canadian geotechnical Journal•Soils and foundations 关键词：×××Updated on June 2, 2015 软件件论文数量文数量FLAC31ABAQUS23ANSYS5PLAXIS14ADINA0PFC2UDEC910为什么要用FLAC？检索期刊：•GEO*•SOILS*关键词：×××Updated on June 2, 2015 软件件论文数量文数量FLAC88ABAQUS50ANSYS27PLAXIS74ADINA3PFC12UDEC1611FLAC能做什么？●岩土工程中的绝大多数问题○土力学、岩石力学、防灾减灾、隧道、地下空间等●采矿工程中的大部分问题●水工结构中的部分问题●结构工程国际通用的岩土工程国际通用的岩土工程专业分析程序专业分析程序12FLAC为何这么流行？●Charles Fairhurst美国工程院、瑞典皇家工程院院士，国际岩石力学学科和岩石力学学会创始人之一，历任国际岩石力学学会主席和副主席，国际岩石力学学会MULLER奖、美国岩石力学学会终生成就奖获得者。

●Peter Cundall美国工程院、英国皇家工程院院士，国际资深计算岩石力学学家13课程安排●第一讲：基本介绍、静力分析、前后处理●第二讲：接触面、FISH语言、流固耦合分析●第三讲：动力分析、自定义本构、结构单元●第四讲：FLAC（2D）基本介绍与应用实例●讨论14课程要求●了解FLAC的基本概念●了解FLA的适用范围和缺陷●了解深入学习FLAC的方法●会用FLAC分析简单的路堤填筑问题（考题）15题外话——“虚拟”与“现实”16现实的复杂与虚拟的简化17是否一定要“像”18我们也做过很像的东西19认清土体数值分析的影响因素参数本构网格软件复杂的网格往往会把核心问题覆盖掉，经复杂的网格往往会把核心问题覆盖掉，经常会造成计算结果无法解释常会造成计算结果无法解释20虚拟是简单的世界现实中的桩现实中的桩虚拟的桩虚拟的桩21虚拟是简单的世界冰碛土体结构模拟结果冰碛土体结构模拟结果台阶坡面上的砾石产出状态台阶坡面上的砾石产出状态22小结●数值分析要对实际工程进行大量的、细心的简化，从效率和结果两方面保证数模分析的成果●结果的判断，需要深入扎实的理论功底及“丰富”的工程经验●数值分析的作用从“锦上添花”到“雪中送炭”，在于各位的努力，指日可待23第一讲FLAC3D基本介绍、静力分析、前后处理24软件介绍●Fast Lagrangian Analysis of Continua●美国Itasca咨询公司开发2D程序(1986)●1990年代初引入中国●有限差分法(FDM)●3D版本：DOS版→2.0 →2.1 →3.0 →3.14.05.0●2D版本：DOS版4.05.06.07.025FLAC●大应变、小应变计算模式.●丰富的本构模型、提供自定义的本构模型功能●接触面可以模拟不同材料的接触●流固耦合实现土体的固结与渗流●拥有各种功能的结构单元类型，模拟土与结构的相互作用●强大的动力分析功能.●流变分析，拥有粘弹性模型和粘塑性模型●热力学分析.Shear strainrate contours26FLAC3D- 与FLAC类似，是FLAC的三维版本- 与 FLAC拥有相同的优点upstreamdownstream 27基本特点●内置材料模型●连续介质非线性，大应变模拟●显式解题方案，为不稳定物理过程提供稳定解●界面或滑动面用来模拟可产生滑动或分离的离散面，从而模拟断层，节理或摩擦边界●内置材料模型丰富：○零模型, ○三个弹性模型 (各向同性，横观各向同性和正交各向异性), ○八个朔性模型 (德鲁克-布拉格, 摩尔-库伦, 应变硬化/软化，单一节理，双线性应变硬化/软化单一节理, 双屈服，修正剑桥粘土，霍克-布朗)隧道工程28可选模块●可选模块包括: ○ 热力学,热-力学耦合,热-流体-力学耦合包括热传导和对流; ○ 粘弹,粘朔性(蠕变)材料模型;○动力学分析,并可以模拟静边界和自由域 ○使用C++定义自己的模型核废料储存中的热力学研究问题核废料储存中的热力学研究问题29FLAC/FLAC3D基本原理●FLAC/FLAC3D利用有限差分，显示方案，动态松弛方法模拟连续体的非线性力学行为：○即使对准静态问题，程序仍然求解完整的动力学方程。

这种方法的好处在于可以为物理非稳定过程例如塌方提供稳定解；○在 “松弛”方法中，使用阻尼来吸收动能以模拟系统的“静态”反应 这种方法可以用比其它方案如解矩阵法更为真实有效地模拟塌方问题30Lagrangian法●源自流体力学中的拉格朗日法○跟踪流体质点的运动状态○跟踪固体力学中结点，按时步用Lagrangian法研究网格节点的运动●节点和单元随材料移动，边界和接触面与单元的边缘一致●固体力学大变形理论法国数学家、物理学家拉格朗日 31混合离散技术+/2=每个每个 为常应力为常应力/应变应变:体积应变由整个四边形算出体积应变由整个四边形算出 . 应变偏量则有两个三角形应变偏量则有两个三角形 和和 分别算出分别算出(混合离散混合离散 过程过程)解题过程中网格坐标按照解题过程中网格坐标按照“拉格朗日方式更新拉格朗日方式更新” (网格随材料移动网格随材料移动), 且为且为显式显式 (一个时步内局部变化不会影响邻域一个时步内局部变化不会影响邻域)32混合离散技术FLAC3D混和离散混和离散+/2=33FLAC3D混和离散n结构域离散为可由四面体单元组合形成的五面体或六面体等单结构域离散为可由四面体单元组合形成的五面体或六面体等单元元; ;n以以 为基本单元为基本单元( (常应力、常应变常应力、常应变););n体应变的计算：体应变的计算： ; ;n偏应变的计算：偏应变的计算： . .34动态松弛●动态松弛法○在动态松弛法中，网格点根据牛顿运动定律运动. 网格点的速度与该点的不平衡力呈正比. 这种求解方法所决定的一系列位移将把系统带入平衡状态，或表明破坏模式.●在动态松弛法中有两个因素很重要:○时步的选择○阻尼效应 35显式算法显式解与隐式解的比较显式解与隐式解的比较显式显式,逐时推进逐时推进隐式隐式, 静态静态1.无需进行反复迭代来实现非线无需进行反复迭代来实现非线性本构关系性本构关系 .2. 类似问题求解时间呈类似问题求解时间呈 N3/2 规规律增长律增长3. 物理非稳定性不会引起数值物理非稳定性不会引起数值不稳定性不稳定性.4. 因为无需储存矩阵，用较小因为无需储存矩阵，用较小内存即可模拟大尺度问题内存即可模拟大尺度问题.5.对大位移、大应变问题同样适对大位移、大应变问题同样适合，无需额外的计算合，无需额外的计算 .1.需进行反复迭代来实现非线性需进行反复迭代来实现非线性本构关系本构关系 2.类似问题求解时间呈类似问题求解时间呈 N2 甚至甚至 N3规律增长规律增长.3.难以模拟物理非稳定性问题难以模拟物理非稳定性问题.4.需存储刚度矩阵，需克服相关需存储刚度矩阵，需克服相关的带宽问题，需要的内存较大的带宽问题，需要的内存较大 .5.对大位移、大应变问题需进行对大位移、大应变问题需进行大量的计算大量的计算 .36New Features in FLAC Version 6.0●使用Intel Fortran compiler拥有更快的计算速度●自动网格重画功能，解决 bad-geometry 问题. ●新的模拟颗粒土材料的硬化模型●更新的通用网格生成工具37New Features in FLAC3D Version 3.11.多处理器的并行计算功能 2.新结构单元类型 “Embedded Liner” 提供两个方向的接触作用，可以很好地模拟挡土墙3.对四面体单元采用新的混合离散方法 “Nodal Mixed Discretization” 提供塑性问题更精确的解答4.64位程序5.包含命令手册、FISH手册和应用实例的帮助38New Features in FLAC3D Version 4.0●模拟颗粒状材料的硬化模型●自动网格重画功能，解决 bad-geometry 问题. ●改进的interface●更快的渗流计算●更新的动力计算功能39Lagrangian格式动量平衡方程F(t)m牛顿运动定律对于连续体在静力平衡条件下，加速度项为0，方程变为平衡方程40自由落体的模拟G = 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 241自由落体的模拟(movie)42FLAC3DFLAC3D中模型术语中模型术语节点gridpoint：节点zone：单元boundary：边界43FLAC3D的求解过程平衡方程(动量方程)应力—应变关系(本构模型)Gauss定律单元积分应变率速度节点力新的应力对所有的网格节点对所有单元44FLAC3D中的本构模型●开挖模型null ●3个弹性模型○各向同性弹性○横观各向同性弹性○正交各向同性弹性 ●8个弹塑性模型○Drucker-Prager模型、Morh-Coulomb模型、应变硬化/软化模型、遍布节理模型、双线性应变硬化/软化遍布节理模型、修正剑桥模型和胡克布朗模型45FLAC3D中的本构模型46一个最简单的例子gen 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 ;求解app nstr -10e4 ran z 3 x 1 2 y 1 2solveplo con zd ;后处理切片功能RUN FLAC3D47分析问题的过程建立网格建立网格初始条件初始条件边界条件边界条件初始应力平衡初始应力平衡外荷载外荷载求解求解前处理前处理后处理后处理48FLAC3DFLAC3D的文件格式的文件格式n保存文件保存文件 (*.sav) – 含有所有状态变量和用户定义条件的二进制文件 n数据文件数据文件 (*.dat) – 数据文件由用户创建的一种ASCⅡ格式的文件，它包括一系列的用于描述所分析问题的FLAC3D命令 nFISH文件文件(*.fis) –FISH程序文件nFLAC3D文件文件(*.flac3d) –FLAC3D的网格信息文件n历史记录文件历史记录文件 (*.his) – 记录输入输出历史值的文件 n图形文件图形文件 – 图形文件(各种标准格式)n电影文件电影文件 (*.dcx) – AVI或PCX图像文件，这些图像文件可以当作电影放映49初始应力的生成●为什么要单独列出？○分析过程中出现的很多问题都与初始应力是否合理有关○手册中的例子五花八门○是所有后续分析的基础！！●生成方法○弹性求解○更改强度参数的弹塑性求解○设置初始应力的弹塑性求解○存在水压力的初始应力生成○水下建筑的初始应力生成50弹性求解模型尺寸单元数量密度KGu1×1×2 (m3)1×1×2200030MPa10MPa0.35gen zon bri size 1 1 2m elasprop bulk 3e7 shear 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 = 162s sz = -40e3s sx = -21.54e351更改强度参数的弹塑性求解模型尺寸单元数量密度KGcfu1×1×2 (m3)1×1×2200030MPa10MPa10kPa150.35gen 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 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 = 163s sz = -40e3s sx = -21.54e3Or: solve elastic52设置初始应力的弹塑性求解模型尺寸单元数量密度KGcfu1×1×2 (m3)1×1×2200030MPa10MPa10kPa150.35gen zon bri size 1 1 2model mohrprop bulk 3e7 shear 1e7 c 10e3 f 15 ten 0fix z ran 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 = 0s sz = -40e3s sx = -20e353存在水压力的初始应力生成(1)模型尺寸单元数量饱和密度KGcfu水位线孔隙率1×1×2(m3)1×1×2200030MPa10MPa10kPa150.351m0.5gen 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 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 1ini 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 = 0s sz = -35e3s sx = -22.5e3r rd = r rs – n×s×r rf54存在水压力的初始应力生成(2)模型尺寸单元数量饱和密度KGcfu水位线孔隙率1×1×2(m3)1×1×2200030MPa10MPa10kPa150.351m0.5config fluidgen zon bri size 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 = 142s sz = -40e3s sx = -24e3单元数较少产生的误差单元数较少产生的误差55水下建筑的初始应力生成模型尺寸单元数量饱和密度KGcfu水位线1×1×2(m3)1×1×2200030MPa10MPa10kPa150.353mgen 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 ran 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 = 0s sz = -50e3s sx = -40e356前后处理57前后处理●基本前后处理○命令操作○菜单操作●dd & dip●attach & merge●外界模型的导入●复杂模型的网格检查58FLAC3D的前后处理●命令驱动(推荐)○程序控制○图形界面接口○计算模型输出○指定本构模型及参数○指定初始条件及边界条件，指定结构单元○指定接触面○指定自定义变量及函数(FISH)○求解过程的变量跟踪○进行求解○模型输出59菜单驱动(计算模式)命令栏60菜单驱动(Plot)61FLAC3DFLAC3D的前处理的前处理FLAC3D网格生成的关键特征：网格生成的关键特征：•FLAC3D 是命令驱动.•使用FLAC3D内置基元进行形状组合可形成复杂网格.•用户自定义FISH 函数可以用来修改基元网格以创建更为复杂的网格.•第三方软件导入.62FLAC3D网格基元块体体退化退化块体体楔体楔体金字塔金字塔四面体四面体柱体柱体径向径向块体体径向隧道径向隧道径向柱体径向柱体柱状壳体柱状壳体柱状交叉柱状交叉隧道交叉隧道交叉6364块体体径向柱体径向柱体柱状交叉柱状交叉65Tecplot的后处理6667Tecplot后处理68697071Tecplot的云图和等值线TecplotTecplot中云中云图的效果与的效果与flac3dflac3d的的shade onshade on的效果差的效果差不多不多, ,但是其出三但是其出三维等等值线的功能是卓越的的功能是卓越的. .7273Tecplot的切片功能TecplotTecplot的的sliceslice切片切片功能功能. .和和CADCAD一一样, ,可以任意可以任意切剖面出切剖面出图. .最大最大优点是可以点是可以几个剖面几个剖面同同时出出图. .74整体和截面的整体和截面的数据以数据以单元形元形式存在列表中式存在列表中, ,可自由可自由选择和和组合多个合多个单元元出出图. .75dd 和 dip●Simwe.Com上有近百的讨论贴●地质上的概念，倾向和倾角●建议用ori和norm代替●后处理的切片功能○plo set plane ori (*,*,*) norm (*,*,*)○plo con zdis plane○概念清晰，使用快捷●其他应用○reflect网格○建立水面xyzorinorm76attach 和 merge●attach○用于连接次节点 (sub-grid)○允许网格存在一定随意性○要求成比例(1:2, 1:3等)○可用于模型的检查○慎用●merge○用于节点的合并○外来模型导入的精度差异○tolerance的设置sub-grid 4:25:2tol121merg77外界模型的导入●复杂网格的生成难度大●接口编写○不同软件之间的精度差异○Group的定义○采用.flac3d的文件格式●.flac3d文件的格式○G 1 1.0e+00 1.0e+00 1.0e+00○Z B8 1 2 3 4 5 6 7 8○ZGROUP Soil1 2 3●impgrid & expgrid○仅限于网格3w zones 3sec78复杂网格的检查●FLAC3D生成的复杂网格○attach face◦无接触面时检查整体模型是否存在无接触面时检查整体模型是否存在sub-grid◦有接触面时给定范围进行检查有接触面时给定范围进行检查●其他软件导入的复杂模型○网格划分的检查○gen merge○弹性模型model elastic求解◦“独立独立”节点节点◦“畸形畸形”单元单元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速度约束条件○plo fap red sk体力○plo hist (unbal)不平衡力○plo interface nstress (sstress)接触面单元81第二讲FLAC3D接触面、FISH语言、流固耦合分析82接触面（interface）83接触面单元●应用范围●原理●建模方法○推荐方法○复杂内部接触面的设置方法●参数选择●单桩承载力分析●挡土墙的接触面设置○思路：未知问题的分析方法84接触面单元的用途●网格不连续●岩体介质中的解理、断层、岩层面●地基与土体的接触●箱、槽及其内充填物的接触●空间中无变形的固定“障碍”85接触面的原理●三角形单元 (无厚度!)●参数较多●三种工作模式○粘结界面○粘接滑移○库伦滑动86接触单元模型的建立 (1)●关键要形成同一位置的两个节点(面)●“移来移去”(推荐)○建两个分开的模型○建立接触单元○通过INI * add使模型接触○注意dist的含义○NO merge, NO attach!!!接触面dist123487接触单元模型的建立 (2)●“导来导去”○利用expgrid, impgrid命令进行网格导出与导入○配合DELETE命令○适于内部接触面的建立，或○其他前处理工具建立的网格88接触单元模型的建立 (3)●GEN separate ●INTERFACE wrap○指定正确的group89接触面参数的确定●虚构的为了合并节点而设置的接触面○Kn=ks=10*max[(K+4/3G)/Dzmin]●真实的刚性接触面○如料仓下料○c,D,Tension重要，kn,ks不重要●真实的柔性接触面○断层；水力劈裂材料○试验得到参数○对于kn,ks：岩石断层10~100MPa/m(粘土); 100GPa(岩石)○反分析方法：通过断层中岩石的变形与原岩的变形90单桩承载力分析●软土地基○bulk 1.6878E6 ○shear 3.6167E5 ○coh 15E3 ○fric 12○dens 1.73E3●桩体○bulk 5e9○shear 3.75e9○dens 2.5e30.5m 8m 10m 20m 91计算过程施加桩顶荷载施加桩顶荷载计算结果计算结果92影响因素7.50E+03fricgk3coh/0.70.7fric10g10k23.00E+0420100e100e1cohfrickskn水平因素93计算工况设计1233(9)3123(8)2313(7)2132(6)1322(5)3212(4)3331(3)2221(2)1111(1)cohfrickskn 水平因素工况94计算结果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水平因素水平因素95最优方案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水平因水平因素素96合理步骤单桩分析单桩分析简单网格简单网格接触面参数接触面参数多次试算多次试算理论、实测理论、实测加密网格加密网格接触参数接触参数理想结果理想结果群桩分析群桩分析Pile结构单元结构单元单元参数单元参数理想结果理想结果Pile结构单元结构单元单元参数单元参数理想结果理想结果97挡土墙的接触面设置●对于未知问题的分析思路○3个独立的接触面○3个同ID的接触面○2个独立的接触面○2个独立的接触面并进行底部merge挡墙挡墙 土体土体 123WallSoil98挡土墙的接触面设置(2)Z = 5.16cmZ = 1.14cmZ = ∞不能初始平衡不能初始平衡 Z = 0.35cm？？99挡土墙的接触面设置(3)X = 1.95cmX = 1.91cmX = 1.47cm？？100挡土墙的接触面设置(3)●前两种方法的差别实质○3个ID的独立接触面在相同位置产生互不影响的两个节点○共同ID的接触面在相同位置自动设置为1个节点●最终的结论需要您自己去判断！3 interfaces, 2 IDs 1 interface, 1 ID 101FISH语言102FISH语言简介●软件自带的编程语言●是否一定要学？○视情况而定，需要时查询FISH变量即可●语法简单xxx……end_xxx●注意事项○与FLAC本身的关键字冲突○保留字不可缩写○变量可不定义，因此注意检查程序◦print fish○table, extra等命令使用103一个最简单的FISH程序def abc abc = 1 + 2 * 3 abcd = 1.0 / 2.0endabcprint fish数据格式104函数与变量1.都可以在FISH函数中进行赋值，赋值操作与常规的编程语言类似，按照运算符的优先级先后顺序来执行。

2.函数和变量的赋值遵守数据类型的规则，即整型的计算结果为整型，浮点型的计算结果为浮点型，因此读者在进行除法运算、开方运算时都需要将数据类型设置为浮点型，数字尽量使用小数点以保证运算正确3.变量和函数名的命名规则是不能以数字开头，不能含有中文，并且不能包含如下的字符 , * / + - ˆ = < > # ( ) [ ] @ ; ’ "4.变量和函数名不能与FLAC3D、FISH的保留字相冲突，不要采用过于简单的单词，比如a，hist等，这些都与保留字相冲突5.即使程序中存在与保留字相冲突的变量，FLAC3D也不会提供任何提示，所以提醒读者在编制FISH程序时尽量使用较长的、复杂的变量和函数名105函数与变量5.对变量进行赋值时，不能使用当前函数的函数名放在“=”的右边，比如采用下面的定义abcd = abc + 1.0在FISH程序执行时会提出错误，因为这样会形成递归调用，这种调用方式在FISH程序中是不允许的6.变量和函数的作用是全局的，在命令中的任何地方修改变量的值都会立即生效，因此在实际应用中尽量避免不同的函数中含有相同的变量，因为这样可能会造成赋值错误，并难以检查7.在FLAC3D中可以用如下的命令来引用FISH函数和变量○PRINT用于查看函数和变量的数值；○HISTORY命令可以对函数和变量的数值进行记录；○SET命令用于变量的赋值。

106主要语句●选择语句CASEOF 表达式…默认语句CASE n1…表达式的值为n1时的语句CASE n2…表达式的值为n2时的语句ENDCASE107主要语句●条件语句IF 条件表达式 [THEN]…[ELSE]…ENDIFFISH中条件运算符没有“并”、“或”、“否”这样的符号 if aa > 1.0 if aa < 2.0 执行语句 endifendif表达“1 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_gpid112FISH的编写习惯第一步第一步第二步第二步def abcendabcdef 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) endloopendabc113FISH与建模两个圆形隧道的连接部分两个圆形隧道的连接部分 变直径的隧道部分变直径的隧道部分 114FISH检查●主要采用PRINT fish命令●查看变量的赋值是否合理，主要检查值为0的函数和变量，因为FISH程序中一般定义的变量都有实际的意义，输出为0的变量很可能是与保留字相冲突的变量○（如a就是apply的保留字）●或者由于编写笔误产生的变量○（如数字0与大写字母O）。

115流固耦合分析很难！116流-固耦合分析(单相流)●基本功能●理论框架●计算模式●渗流边界条件，初始条件●单渗流计算及渗流耦合计算 117基本功能●渗流各向同性、各向异性●不同的渗流模型和属性●流体压力，涌入量，渗漏量和不渗水边界●抽水井、点源、体积源●饱和渗流可采用显式差分法、隐式差分法 ●非饱和渗流采用显式差分法●渗流-固体-热的耦合●流体和固体的耦合程度依赖于土体颗粒(骨架)的压缩程度，用Biot系数表示颗粒的可压缩程度 ●循环荷载引起的动水压力变化和土体液化118理论框架●准静态Biot理论 ●多孔介质中遵循Darcy定律的单相渗流●描述多孔介质中流体渗流的变量○孔隙水压力，饱和度，特定排水向量的三个分量 ●质量守恒定律●达西定律●本构定律○考虑流体响应孔隙水压力改变，饱和度改变，体积应变改变和温度改变 119有效应力计算●无渗流模式○不设置CONFIG Fluid○孔压不改变○设置孔压分布◦INITIAL pp◦WATER table◦WATER density◦WATER table face◦SET gravity○手动设置干湿密度F渗流模式§设置CONFIG fluid§设置土体干密度 §渗流模型•MODEL fl_isotropic •MODEL fl_anisotropic •MODEL fl_null §SET fluid off§set WATER bulk = 0120渗流边界条件，初始条件●默认的边界条件是不透水边界 ●孔隙压力自由(不透水边界 )●固定孔隙水压力(透水边界) ○如：井 ●孔隙压力，孔隙率，饱和度和流体属性的初始分布可以用INITIAL命令或者PROPERTY命令定义。

121单渗流计算及渗流耦合计算 ●时间比例 ●完全耦合分析方法 ●孔压固定分析(有效应力分析) ●单渗流得到孔压分布 ●无渗流计算——孔压的力学响应 ●流-固耦合计算 122单渗流得到孔压分布●用途：排水沟；抽水井；耦合计算●计算步骤○CONFIG fluid SET mech off○SET fluid implicit on/off○MODEL fl_; PROP ○STEP; SOLVE age; SET fluid ratio○SET fluid off mech on◦PROP biot_c 0 (or INI fmod 0)123无渗流计算——孔压的力学响应●不排水短期响应●两种分析方法：干法和湿法○干法：Ku=K+a2M◦两种破坏形式两种破坏形式◦WATER或INI获得常孔压，不排水的c,φ (孔压改变较小)◦φ=0,c=cu (M>>K+4/3G)○湿法：耦合体系的短期行为◦使用排水的使用排水的K, c,φ◦若若SET fluid off, Biot_mod(fmod)真实真实124力学过程的特征时间力学过程的特征时间 流体扩散过程的特征时间流体扩散过程的特征时间 完全耦合分析方法•时间比例○短期行为 (不排水)◦ ts(分析时间分析时间)<>tc•施加扰动的属性○流体扰动：渗流可不与力学过程耦合○力学扰动：耦合等级取决于流固刚度比•流固刚度比125流-固耦合计算●CONFIG fluid; M(Kf); K(渗透系数) 真实，则FLAC3D默认耦合计算○Δp→Δev→s○Δev→Δp●预估流/力特征时间●耦合计算前先达到一个平衡状态○SET fluid on mech off; SET fluid off mech on; STEP○SET mech force; SET mech substep n auto; SET fluid substep m (=1)○STEP:渗流步足够小126流固耦合的计算方法●手动调整的STEP求解◦SET fluid on mech off◦STEP◦SET fluid off mech on◦STEP●主从进程的SOLVE求解◦SET mech force◦SET mech substep n auto (从进程)◦SET mech substep m(主进程)◦SOLVE age●自动STEP求解◦STEP127渗流问题(CONFIG fluid)分析步骤●时间比例(ts, tc)○稳态○不排水状态○相当●扰动类型○力学扰动○孔压扰动●流固刚度比Rk○是否>>>1●完全耦合模式○时间比例相当；力学扰动128心墙土坝的渗流 (1)newconfig fluid set fluid 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 am 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网格模型初始孔压129心墙土坝的渗流 (2)rest elastic.savini 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竖向应力沉降130心墙土坝的渗流 (3)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 215solveRatio = 1Ratio = 1E-5131荷载引起的地基土体的超孔隙水压力 10 m3 mK = 500 MPa, G = 300 MPa, c = 10 kPa f = 1540 kPa20 m132计算文件; --- 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 surface ---fix pp 0 range z 9.9 10.1; --- settings ---set fl off; --- test ---step 750[1][1]因为本例中没有设置初始应力，这里只进行了750步的求解。

133计算结果134对主从进程法的讨论SET mech force 设置一个不平衡力的大小，达到这个不平衡力系统认为暂时达到平衡状态；SET 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的收敛准则既可以满足计算流固耦合过程中的计算精度要求，同时又具有较高的计算效率 135对主从进程法的讨论(2)讨论讨论2：子步数的影响：子步数的影响 set mech sub 100 fluid sub 10set mech sub 10 fluid sub 10set mech sub 1 fluid sub 1设置合理的子步数也很重要设置过大，则会导致计算时间大大增加，过小又会造成计算结果的误差 136真空预压的简单模拟•孔压边界条件•ts>>tc○长期分析(排水)•Rk>>1○骨架很软•孔压扰动○进行biot_mod调整砂层软土层粘土层PVD2m8m10mData file: Data file: 137数值分析过程(movie)138课程安排●第一讲：基本介绍、静力分析、前后处理●第二讲：接触面、FISH语言、流固耦合分析●第三讲：动力分析、自定义本构、结构单元●第四讲：FLAC（2D）基本介绍与应用实例●讨论139FLAC3D非线性动力分析非常复杂！Said by Prof. Peter Cundall140为什么要用FLAC做动力分析？●FLAC 可以模拟体系（土，岩石，结构，流体）受到的外部动力荷载（比如地震）或内部动力荷载（比如基础振动、爆炸）。

●可以计算塑性引起的永久变形以及孔隙水压力的消散●土动力学中常用的等效线性方法无法直接处理上述问题 141动力模拟的3个重要问题1.动力荷载与边界条件2.材料响应与阻尼3.土体液化142动力荷载●动力输入的类型○加速度时程○速度时程○应力(压力)时程○力时程●APPLY ○INTERIOR (内部)○TABLE○FISH143Quiet边界●静态(quiet,粘性)边界○Lysmer and Kuhlemeyer(1969) ○模型边界法向和切向设置独立的阻尼器●性能○对于法向p波和s波能很好的吸收○对于倾斜入射的波和Rayleigh波也有所吸收，但存在反射○人工边界仍应当足够远144Quiet边界应用●内部振动(如隧道中的列车振动问题)○动力荷载直接施加在节点上○使用Quiet边界减小人工边界上的反射○不需要FF边界●外部荷载的底部边界○软土地基上的地震荷载不适合用加速度或速度边界条件○使用应力条件t = -2Csrvs●地震底部输入的侧向边界○扭曲了入射波quietquietquiet145Free-field边界•Cundall et al. (1980)•自由场网格与主体网格的耦合粘性阻尼器，自由场网格的不平衡力施加到主体网格边界上•设置条件○底部水平，重力方向为z向○侧面垂直，法向分别为x, y向○其他边界条件在APPLY ff之前相当于一个阻尼器146Free-field边界●APPLY ff将边界上单元的属性、条件和变量全部转移ff单元上；●设置以后主体网格上的改动将不会被FF边界所响应●可存在任意的本构模型以及流体耦合(仅竖向)●FF边界进行小变形计算，主体网格可大变形，FF边界上的变形要相对较小●存在attach的边界将不能设置FF边界●边界上的Interface将不能连续●动力边界设置需在FF边界设置之前147Free-field边界与动力荷载•模型底部边界○fix——施加速度或加速度荷载——刚性边界○Free——施加应力时程荷载——柔性边界•对于软弱的地基不适合施加速度（加速度荷载），而应当施加应力荷载Note that there is a factor of 2 because the input energy divides into a downward- & upward-propagating wave.1482. 材料响应与阻尼1.连续的非的非线性，表性，表观模量随着模量随着应变的增大而降低的增大而降低2.对所有循所有循环应变等等级均存在滞回均存在滞回特性，因此特性，因此导致随着循致随着循环应变的的增加阻尼比增大。

阻尼是率相关增加阻尼比增大阻尼是率相关的3.对于复于复杂波形的各个成分都波形的各个成分都产生生阻尼4.剪切剪切应变会会产生的体生的体积应变，相，相应的，随着剪的，随着剪应变循循环次数的增次数的增加体加体积应变逐逐渐积累149材料响应土体在循土体在循环荷荷载作用下呈作用下呈现出模量衰减和能量消散的特点，出模量衰减和能量消散的特点，那么如何用非那么如何用非线性数性数值方法方法对其其进行模行模拟呢？呢？Nonlinear characteristics of soils (Martin and Seed, 1979)150试验得到的阻尼比、割线模量随循环剪应变的曲线试验得到的阻尼比、割线模量随循环剪应变的曲线 0.00010.0010.010.11Shear Strain Amplitude (%)01020304050Damping Ratio (%)0.00010.0010.010.11Shear Strain Amplitude (%)0.00.10.20.30.40.50.60.70.80.91.0Normalized Shear Modulus, G/GmaxMid-Range Sand Curve(Seed & Idriss, 1970)Sand Fill Inland: Friction =32, hr=0.47, Go=440Sand Fill under Rock Dike: Friction=30, hr=0.43, Go=440151等效线性方法•等效等效线性方法是岩土地震工程中模性方法是岩土地震工程中模拟波的波的传播的最常用的方法。

播的最常用的方法•假定土体是粘假定土体是粘弹性体，参照性体，参照实验室得到的切室得到的切线模量及阻尼比与剪模量及阻尼比与剪应变幅幅值的关的关系曲系曲线，，对地震中每一地震中每一单元的阻尼和模量重新元的阻尼和模量重新赋值Iteration toward strain-compatible shear modulus and damping ratio (after Kramer, 1996)152等效线性方法的特点1.使用振动荷载的平均水平来估算每个单元的线性属性，并在振动过程中保持不变在弱震阶段，单元会变得阻尼过大而刚度太小；在强震阶段，单元将会变得阻尼太小而刚度太大对于不同部位不同运动水平的特性存在空间变异性2.不能计算永久变形等效线性方法模型在加荷与卸荷时模量相同，不能计算土体在周期荷载作用下发生的剩余应变或位移3.塑形屈服模拟不合理在塑性流动阶段，普遍认为应变增量张量是应力张量的函数，称之为“流动法则”然而，等效线性方法使用的塑性理论认为应变张量（而不是应变增量张量）是应力张量的函数因此，塑性屈服的模拟不合理4.大应变时误差大等效线性方法所用割线模量在小应变时与非线性的切线模量很相近，但在大应变时二者相差很大，偏于不安全。

5.本构模型单一等效线性方法本身的材料本构模型包括了应力应变的椭圆形方程，这种预设的方程形式减少了使用者的选择性，但却失去了选择其它形状的适用性方法中使用迭代程序虽然部分考虑了不同的试验曲线形状，但是由于预先设定了模型形式，所以不能反映与频率无关的滞回圈另外，模形是率无关的，因此不能考虑率相关性153完全非线性分析方法 FLAC3D采用完全非线性分析方法，基于显式差分方法，使用由周围区域真实密度得出的网格节点集中质量，求解全部运动方程 154完全非线性分析方法的特点1.可以遵循任何指定的非线性本构模型如果模型本身能够反映土体在动力作用下的滞回特性，则程序不需要另外提供阻尼参数如果采用Rayleigh阻尼或局部（local）阻尼，则在动力计算中阻尼参数将保持不变2.采用非线性的材料定律，不同频率的波之间可以自然地出现干涉和混合，而等效线性方法做不到这一点3.由于采用了弹塑性模型，因此程序可以自动计算永久变形4.采用合理的塑性方程，使得塑性应变增量与应力相联系5.可以方便地进行不同本构模型的比较6.可以同时模拟压缩波和剪切波的传播及两者耦合作用时对材料的影响在强震作用下，这种耦合作用的影响很重要，比如在摩擦型材料中，法向应力可能会动态地减小从而降低土体的抗剪强度。

155使用弹塑性模型●附加考虑的因素：○阻尼，对于屈服面以下应力的循环○体积应变积累，是循环周数与幅值的函数○模量衰减，基于平均应变水平的表格最简单的弹塑性模型往往在描述累计塑性应变方面具有很好的效果，但是对于加速度放大系数的估算上效果不好156弹塑性模型简单的理想弹塑性本构模型仅仅在发生屈服时才会出现滞回特性strainstress注意：    即使这样粗糙的模型也能够作出连续的阻尼比和模量衰减曲线在屈服条件下会产生体积改变，但通常都是剪胀157FLAC3D中的阻尼比1.弹塑性本构模型使用瑞利（粘性）阻尼2.弹塑性本构模型使用滞后（HD）阻尼.3.复杂本构模型拥有连续的屈服应力应变关系和对应的加卸载响应（Wang，UBCSand）158瑞利阻尼瑞利阻尼瑞利阻尼最初应用于结构和弹性体的动力计算中，以减弱系统的瑞利阻尼最初应用于结构和弹性体的动力计算中，以减弱系统的自然振动模式的振幅在计算时，假设动力方程中的阻尼矩阵自然振动模式的振幅在计算时，假设动力方程中的阻尼矩阵C C与与刚度矩阵刚度矩阵K K和质量矩阵和质量矩阵M M有关：有关： 瑞利阻尼中的质量分量相当于连接每个节点和地面的阻尼器，而刚瑞利阻尼中的质量分量相当于连接每个节点和地面的阻尼器，而刚度分量则相当于连接单元之间的阻尼器。

虽然两个阻尼器本身是与度分量则相当于连接单元之间的阻尼器虽然两个阻尼器本身是与频率有关的，但是通过选取合适的系数，可以在有限的频率范围内频率有关的，但是通过选取合适的系数，可以在有限的频率范围内近似获得频率无关的响应近似获得频率无关的响应 159瑞利阻尼瑞利阻尼frequencyratio of damping to critical混合仅有刚度分量仅有质量分量采用叠加的方法得到的阻尼比在较大的频率范围内保持定值 （3:1）Combined curve reachesminimum at:160中心频率的选择 假假设弹性模型性模型计算土石算土石坝不同材料的功率不同材料的功率谱曲曲线 (assuming elastic material)frequencyfrequencyfrequencyfrequencyfrequency161阻尼比的选择 根据根据65%最大最大应变来来选择阻尼比和模量衰阻尼比和模量衰减参参数减参参数162瑞利阻尼的缺点●必须指定中心频率，有时须在一些充满矛盾的数据当中进行选择（场地响应或地震平均频率）●刚度比例项会导致时间步减小163滞后阻尼FLAC Version 5.0 and FLAC3D Version 3.0 提供了滞后阻尼功能，提供了滞后阻尼功能，该项阻尼独立于材料模型之外，阻尼独立于材料模型之外，实质是包含了一个与是包含了一个与应变幅幅值相关相关的切的切线剪切模量乘子剪切模量乘子如果割如果割线模量是由衰减曲模量是由衰减曲线来确定来确定的，那么切的，那么切线模量可以表示模量可以表示为：：From Seed & Idriss (1970)Go = small-strain shear modulus164滞后阻尼表表观应变是偏是偏应变，根据先前，根据先前反反转的点的点进行累行累计得到的。

使得到的使用堆用堆栈的数据的数据结构可以描述构可以描述应变反反转现象因此，在“迷你滞回圈”中的能量损失与主滞回圈的计算是一样的elastic model with hysteresis damping特定的曲特定的曲线可以使用滞后阻尼可以使用滞后阻尼来描述，并在本构模型来描述，并在本构模型计算中算中根据根据应变的大小的大小获得不同的切得不同的切线模量165滞后阻尼的特点1.可以直接采用动力试验中的模量衰减曲线；2.相对于瑞利阻尼而言，滞后阻尼不影响动力计算的时间步；3.可以应用于任意的材料模型，且可以与其它阻尼格式同时使用一个缺点：模量衰减曲线一致，但阻尼比曲线存在差别“Good” fit to Seed & Idris data for G/Gmax (sigmoidal 3-parameter function) – note inconsistent damping result.G/GmaxD - % of critical166滞后阻尼●低循环应变下得到的阻尼比要小于试验结果，这会导致低级的噪声，尤其在高频情况下可以在中心频率上增加一个小量的Rayleigh阻尼(~0.2%刚度比例)，这样也不会降低时步；●若初始应力不为0，剪应力-剪应变曲线可能不匹配。

因此在生成初始应力时就要调用Hyst阻尼；●Hyst阻尼不仅会增加能量损失，还会导致在大循环应变下的平均剪切模量的降低，在输入波的基频接近共振频率的时候，由于可能会导致动力反应幅值的增大；●Hyst阻尼之前要做一次弹性无阻尼求解，以获得发生循环应变的最大水平，若循环应变过大导致剪切模量过多的降低，那么用Hyst阻尼是有问题的；●即使应变较小，使用屈服模型也会增大应变，因此若有广泛屈服的现象，则使用屈服模型，不用Hyst阻尼167动孔压的生成——液化●干沙剪应变循环加载试验○初始加载阶段，沙土通常先压实再膨胀卸载时，沙土遵循与加载相似的路径，但在零应变时，有些残余体积应变存在取决于初始孔隙率，这可能代表纯粹的压实●假定孔隙中充满水○对于常体积测试，有效应力降低，孔隙水压保持不变○对于常荷载测试，(例如，盒子上法向荷载固定)，孔隙水压增加，有效应力减小●有效应力为零时发生液化168动孔压的生成——液化●因此孔隙水压增加不是液化的基本原因●由于颗粒间 (重组以后) 的低接触力导致有效应力的减小●描述液化的模型○高级模型：BSHP (边界面低塑性本构模型, Wang et al. 1990)○简单模型：MC + 体积应变增量模型◦Finn模型：模型：◦Byrne模型：模型：169UBCTOT: Post-liquefactionAssign post-liquefaction propertiestgt170岩土介质爆炸液化模型动孔压模式应力力/应变孔孔压模型模型有效有效应力力/瞬瞬态孔孔压模型模型内内时/能量孔能量孔压模型模型 基于基于Byrne 有效有效应力力弹塑性模型，提出塑性模型，提出一个能一个能够考考虑主主应力力轴旋旋转、、饱和砂土含和砂土含有少量气体、有少量气体、饱和砂和砂土液化后的土液化后的应变软化化和和应力重分布特性的力重分布特性的弹塑性模型。

塑性模型模型模型应用用——大大坝斜坡斜坡场地爆炸液化地爆炸液化图1.有效有效应力等力等值线图2.大大坝震后液化区震后液化区预测171PL-Finn开始初始液化判断?零有效应力计算YES非零有效应力计算零有效应力判断?NO普通的Finn模型YESNO动力时间？结束NOYES0：零有效应力状态1：非零有效应力状态0有效应力非0有效应力172地震波的调整●基线校正○对于地震分析的加速度时程，其积分得到的速度和位移应归0○美国地质调查研究所◦Basic Strong-Motion Accelerogram Processing Software (BAP)○对网格施加一个固定速度从而使残余的位移变为0●动力荷载的频率与单元尺寸的双向调整○高频的输入要求单元尺寸很小○一定的单元尺寸对应输入的最大频率○一般进行滤波处理◦滤掉低能量的高频滤掉低能量的高频◦FFT.FIS◦Origin◦SeismoSignal173地震波的调整El-Centro波FFT修正后的时程修正后FFT5Hz174PL-FinnPL-Finn模型的应用模型的应用●阪神地震的码头分析沉箱顶部水平残余位移最大达5m，平均为3.5m，残余沉降为1 ~ 2m，海侧倾斜角3 ~ 5 o。

sea Siltreplaced sandstoneland SandcaissonclayFLAC3D grid8050 zones10386 gps可液化的砂土可液化的砂土可液化的砂土可液化的砂土175分析参数Group nameConstitutive modelrd(kg/m3)E(MPa)mc(kPa)j(o)clayMC1350500.333020sea siltMC1250200.33030replaced sandMC1350150.33037land sandMC135013.70.33036stoneMC15501000.33040caissonElastic350020000.17‐‐Group nameFluid modelK (cm/s)PorositynDamping ratio Dliquefied parametersPL-Finn parametersclayfl_iso1.0E-60.450.05‐‐sea siltfl_iso1.0E-050.450.05‐‐replaced sandfl_iso1.0E-030.450.05Byrne ModelDr = 40 %C1 = 0.751C2 = 0.533C3 = 0pprc = 0.99k0 = 3105.4n0 = 0.3225k1 = 5503.1n1 = 0.1739land sandfl_iso1.0E-030.450.05Byrne ModelDr = 25 %C1 = 2.432C2 = 0.164C3 = 0pprc = 0.99k0 = 3105.4n0 = 0.3225k1 = 5503.1n1 = 0.1739stonefl_iso1.0E-010.450.05‐‐caissonfl_null‐‐0.05‐‐基本力学参数 流体参数动力参数PL-Finn模型参数176输入加速度南北向地震波 最大值 = 0.6 g竖直向地震波 最大值 = 0.2 g177接触面与自由边界条件沉箱与填石之间的接触面模型周围自由场边界条件178计算结果●位移震后水平位移等值线沉箱顶点处的水平位移、沉降时程曲线3.44 m1.99 m179计算结果（续）●超孔压比震后的超孔压比云图超孔压比时程曲线 A: 置换砂 B: 海砂（底部）C: 海砂（中部）180自定义本构模型181自定义本构模型的基本方法●必要性○试验总结的本构模型○特定条件下的本构模型○交叉学科的本构模型●二次开发环境●自定义本构模型的功能●自定义本构模型的基本方法182二次开发环境●FLAC3D采用面向对象的语言标准C++编写●本构模型都是以动态连接库文件(.DLL文件)的形式提供●VC++2005或更高版本的开发环境●优点○自定义的本构模型和软件自带的本构模型的执行效率处在同一个水平 ○自定义本构模型(.DLL文件)适用于高版本的FLAC(2D)、3DEC、UDEC等其他Itasca软件中183自定义本构模型的功能●主要功能：对给出的应变增量得到新的应力●辅助功能：○模型名称、版本○读写操作●模型文件的编写○基类(class Constitutive Model)的描述○成员函数的描述○模型的注册○模型与FLAC3D之间的信息交换○模型状态指示器的描述 184自定义本构模型的基本方法●头文件(usermodel.h)中进行新的本构模型派生类的声明○修改模型的ID(>100)、名称和版本○修改派生类的私有成员●C++文件(usermodel.cpp)中修改模型结构○(UserModel::UserModel(bool bRegister): ConstitutiveModel)●const char **UserModel::Properties()函数○模型的参数名称字符串●const char **UserModel::States()函数○计算过程中的状态指示器185自定义本构模型的基本方法●double UserModel::GetProperty()和void UserModel:: SetProperty()函数●const char * UserModel::Initialize()函数○参数和状态指示器的初始化，并对派生类声明中定义的私有变量进行赋值●const char * UserModel::Run() 函数○由应变增量计算得到应力增量，从而获得新的应力●const char * UserModel::SaveRestore()函数○对计算结果进行保存。

●程序的调试○在VC++的工程设置中将FLAC3D软件中的EXE文件路径加入到程序的调试范围中，并将FLAC3D自带的DLL文件加入到附加动态链接库(Additional DLLs)中，然后在Initialize()或Run()函数中设置断点，进行调试；○在程序文件中加入return()语句，这样可以将希望得到的变量值以错误提示的形式在FLAC3D窗口中得到186相关文件mohr.hmohr.cpp187一个例子(Duncan-Chang)188其他成功的例子●南京水科院双屈服面模型双屈服面中心截面沉降云图水平位移云图小主应力云图189主要结论●FLAC的动力反应分析非常复杂，用户在进行动力分析之前要具有熟练的静力分析基础●FLAC的动力计算需要花费很多时间，主要原因是由于有限差分法本身需要较小的时间步，因此建议在进行动力分析前首先进行简单模型的分析●“小例子”虽然会额外花费一定时间，但对于熟悉边界条件、荷载施加、阻尼选择等内容来说非常必要190第四讲FLAC（2D）基本介绍与应用实例FLAC/FLAC3D的学习方法及建议191FLAC●GIIC基本操作●Interface的建模●流体计算192GIIC基本操作●差分网格的建模思路●成层土建模●养成Mark的好习惯●修剪模型●撤销的最好方法模型草图模型草图计算网格数量计算网格数量确定模型生成确定模型生成方法方法不规则形状生成不规则形状生成Attach, interface参数化参数化运行FLAC193Build选项卡194Alter195Material196In-situ197StructureBeam，梁单元——可以用于各种类型的支护模拟，包括开挖支护、隧道中的支架等。

可以在其两侧连接Interface单元来模拟岩土介质中的挡土墙，还可以通过Interface单元与FLAC网格相连以模拟土工格栅Liner，衬砌单元——主要用于隧道衬砌，包括混凝土初衬或喷射混凝土初衬Cable，锚索单元——不能承受弯矩，可以施加预应力，常用于模拟受拉构件，包括岩石中的锚杆Pile，桩单元——常用于模拟地基中的桩Rockbolt，岩石锚杆单元——常用于模拟岩石中的锚杆Strip，条形锚单元——用来模拟加固提防或土坝中的多层条带型结构Support，支撑单元——用来模拟液压支柱、木质支撑等198Utility199Setting200Plot201Run202Interface单元的建模●移来移去的方法●Interface建模的“方向”○aside, bside的方向一致203含Interface的回填过程●Interface并不随着单元的“回填”而持续作用●必须每回填步Remove和Interface i未创建新接触面的位移云图未创建新接触面的位移云图删除并创建的位移云图删除并创建的位移云图开始状态的应力云图开始状态的应力云图204流体计算模式●单相流○基本类似于3D的计算○快速渗流●两相流(Two phase flow)○气◦空气：非饱和流动空气：非饱和流动◦瓦斯：采矿工程瓦斯：采矿工程◦有毒气体：环境岩土有毒气体：环境岩土●其它流体○混凝土：预制桩的施工模拟○水泥浆：注浆；浆固205应用实例20 m1:1.58 m12 m粘土，r=1800, E=4MPa, u=0.35,c=5 kPa, f=10砂土，r=1800, E=15MPa, u=0.35, c=0, f=25填土，r=1800, E=8MPa, u=0.35,c=10 kPa, f=104 m206网格划分思路●对称性●分块10 m8 m12 m4 m①②⑤⑥③④64 m245812207几个建议1. 了解FLAC和FLAC3D的适用范围、优点和局限性○任何一种方法都是有一定的适用范围，并不能不能解决所有问题，这学要读者对所使用工具的优点和局限性有清醒的认识。

数值模拟的最终目的是为工程问题的诊断和解决提供服务的，需根据问题的本质选择合适的方法和工具；而非“膜拜”和迷信某种方法，机械地用它去套工程，本末倒置○2D永远要比3D简单、实用208几个建议2. 由简到繁，循序渐进○遵循“由简到繁，循序渐进”的学习方法，切忌盲目求大求全，期望一口气吃成胖子学习时，可进行少量单元的简单数值试验来理解软件的特点和功能，积累一定的经验后再进行复杂的数值模拟试验如果注定要错如果注定要错100次，我宁愿错次，我宁愿错1次的时间次的时间是是10分钟，而不是分钟，而不是10小时从从800,000到到2,000从从2个月到个月到5天天209几个建议3. 充分利用手册○手册是最权威的软件说明书，一定要充分利用尽管FLAC和FLAC3D的手册编制顺序不一定适合中国读者的思维习惯，但应尽量养成查阅手册的习惯，做到常翻常新手册中的例子大多都是为了说明某个特定的问题而设定的，因此在讲述该问题时往往会忽略与该问题无关的一些细节，比如参数选择等，因此读者在学习手册时不要“迷信”某个特定的例子，也不要“纠缠”于某些无关的细节，而是要从这些例子中掌握分析问题的基本方法○也不要“迷信”我讲的这个PPT210几个建议4. 了解计算中每一条语句的含义○初学者由于对FLAC和FLAC3D软件了解的不多，在计算时往往会直接套用软件手册或教科书中的例子，而对例子中某些语句的含义并不是真正的了解，这些“不明其意”的语句往往是造成计算结果不合理的原因。

这里建议读者在使用FLAC和FLAC3D程序时，要对自己编写的命令文件中的每一条语句都有清晰的认识和了解，这就要求读者要勤查手册、注重平时的积累set small是什么意思？211几个建议5. 多做“数值试验”○FLAC和FLAC3D程序功能强大，内容众多，在分析具体问题时，读者往往会遇到如法解决的新问题，这些问题在软件手册或教科书中都很难找到答案，这时读者应该多做一些小的算例，开展数值试验，从而了解程序的功能，达到解决问题的目的212几个建议6. 使用“?”○FLAC3D的命令很多，在初学者看来，记住数量可观的各种命令及语句格式是一件很困难的事情，事实也的确如此幸运的是，FLAC3D在命令窗口中提供了“?”功能，无论在命令的什么位置都可以插入“?”字符，让系统告诉你接下来可以应该输入的是哪些关键字或变量RUN FLAC3D213几个建议7. 夯实知识基础○FLAC和FLAC3D的计算结果和中间时步表现出一些不合实际的结果，需要读者具有足够的专业和数学知识进行判断与解释因此，决定FLAC和FLAC3D使用水平高低的决定性因素取决于使用者的专业素养、工程经验和数理知识因此加强专业知识、数学和力学的学习，夯实知识基础十分重要。

214几个建议8. 相互交流，取长补短○FLAC和FLAC3D命令、关键词和变量繁多，个人学习难免顾此失彼，因此加强交流，与他人共享学习经验是提高FLAC和FLAC3D应用水平的一个捷径互联网的出现，为大家提供了一个讨论和共享的平台，读者可以在相互间的交流、争论中取长补短，共同提高215FLAC/FLAC3D并非无所不能●长期动力问题模拟○动三轴试验、动扭剪试验的模拟○列车荷载、车辆荷载的长期作用●复杂网格模型的建立（前处理）不方便●流固耦合问题需要简化处理，简单的固结问题用FLAC来分析时会异常复杂216谢谢！敬请批评指正！217。