06-边坡变形-flac3d前后处理.ppt

边坡岩石力学与工程 flac3d前后处理 河海大学岩土工程科学研究所石崇 2 85 Flac3d前后处理和三维模型构建方法 主要认识和结论 FLAC3D计算结果的后处理 FLAC3D计算中的应力场 FLAC3D三维模型的快速构建 Sufer的三维出图和数据导出 目录 CAD的前处理和数据导出 3 85 在CAD中选择建模区域命令 Trim对等高线截断Move平移 减小坐标起始值导出为R12格式的dxf文件 CAD的数据处理 4 85 模型区域等高线的导出 借助第三方软件GID导入CAD生成的dxf文件 将CAD中的等高线信息读出 GID可以作为CAD和ANSYS的过渡平台 其可转换DXF IGES格式 5 85 在GID软件中显示等高线上的点信息 将等高线的点列为文本输出 作为sufer三维可视化出图的高程点数据 6 85 sufer出图只需要三维的高程点数据 格式写成如下 X1 Y1 Z1X2 Y2 Z2 也可利用其它软件或编程读出CAD中的等高线数据信息 7 85 Flac3d前后处理和三维模型构建方法 主要认识和结论 FLAC3D计算结果的后处理 FLAC3D计算中的应力场 FLAC3D三维模型的快速构建 Sufer的三维出图和数据导出 目录 CAD的前处理和数据导出 8 85 9 85 在滤波器中过滤CAD中导出的局部失真地形点 设置z方向上某一高程区域内的地形点为有效地形点 如模型范围实际地形高程在500 2000以内 则输入z2000 10 85 在插值方法中选择合适的插值方法 注意sufer只能提供栅格化插值方法 如要插值不规则网格 需自己编程实现 设置插值区域 X方向MinMaxY方向MinMax 11 85 Sufer可以生成多种效果图形 对于三维地形的显示 主要采用线框图和表面图 其中表面图的效果最好 其插值效果受插值点个数和插值方法制约 12 85 地图 三维旋转可以旋转查看地形 设置光照 网格可以显示不同的效果 注意在sufer插值后 生成一个后缀名为grd的文件 此为存贮栅格化插值点三维信息的数据文件 13 85 用sufer打开生成的dgx grd文件 选择菜单 文件 文件另存为在保存类型中的下拉菜单中悬着ASCIIXYZ dat格式 即为sufer生成的插值点数据信息 选择 选项 网格信息 可以查看插值的精度和网格间距 14 85 flac3d执行命令流command txt 即可读入sufer生成的地表数据文件 投影生成flac3d ansys和3dec的模型的命令流 15 85 公路 1875m 16 85 AutoCAD在DWG文件中以实体名存储各种点 线 面 实体 实体名是AutoCAD中专门设置的一种特殊数据类型 是赋予图形实体的数字标号 通过AutoLisp语言开发的程序可以访问各种点 线 面 实体信息 例如边坡中的地形线 开挖面 地下厂房 引水洞 平硐等现阶段常用的插值方法有 Kriging法 Delaunay法 Polymesh法 Shepard法 Multiquadric法等等 滑动最小二乘法 MovingLeastSquaresMethod 简称MLSM 是一种精度较高的方法 其克服了经典最小二乘法在拟合中的不足 实质为方差泛函极小化序列与最小二乘法 Sufer中都是栅格插值方法 所以如果要按给定任意网格插值的话 需要自己编程 17 85 18 85 Flac3d前后处理和三维模型构建方法 主要认识和结论 FLAC3D计算结果的后处理 FLAC3D计算中的应力场 FLAC3D三维模型的快速构建 Sufer的三维出图和数据导出 目录 CAD的前处理和数据导出 19 85 导入sufer插值的栅格点 在cad中按精度需要建立多条多义线polyline 类似与ansys中的askin命令的做法 在CAD也可以生成由多条多义线过渡圆滑的nurbs曲面 CAD操作命令Loft 按顺序选择几条多义线 回车生成面域 20 85 在CAD中建立底面域 由封闭的多义线围成面 利用REGION生成面域 采用EXTRUDE命令按Z方向拉伸成体 21 85 EXTRUDE命令按Z方向拉伸成体 让其延伸至实际地面以上一定距离 为后面的CAD中布尔切割操作做准备 22 85 利用slice命令中的曲面切割功能 将底面拉伸的块体沿实际地貌切割为两个块体 同理断层和地层的切割也采用slice中的三点法或曲面法 采用slice切割出的地层与地貌交线 23 85 利用CAD中的EXPORT命令将生成的实体SOLID模型输出为SAT格式的数据文件 ANSYS可以直接读入CAD导出的SAT模型 24 85 在ANSYS进行网格剖分 建议全部剖分为四面体 可以设置线的剖分段数来控制局部和整体的网格质量 由于FLAC3D中最终计算的单元均为四面体单元 因此采用四面体剖分精度影响不大 25 85 FLAC3D数值模型分解后的最小单元类型为Brick 六面体 Wedge 三棱柱体 Pyramid 金字塔体 和Tetrahedral 四面体 单元 26 85 在ansys选择菜单 FILE ReadInputfrom选择APDL命令流文件 在仿真论坛可以下载此程序 执行后在电脑的D 下生成ansys单元节点信息的两个文件01 node dat02 ele dat 27 85 将D 下生成ansys单元节点信息的两个文件01 node dat02 ele dat拷至程序目录下 运行flac 6weishu exe后 程序生成了flac3d可以直接读入的模型数据 28 85 打开flac3d 读入生成的模型数据 Flacmodel haitang Flac3D读入后 保存为sav文件 29 85 30 85 结合CAD的extrude slice union subtract intersect等三维布尔操作功能 编制VBA AUTOLISP程序可以快速构建更为复杂的实体模型 保存为sat格式后 可以直接被ansys读入 31 85 Flac3d前后处理和三维模型构建方法 主要认识和结论 FLAC3D计算结果的后处理 FLAC3D计算中的应力场 FLAC3D三维模型的快速构建 Sufer的三维出图和数据导出 目录 CAD的前处理和数据导出 32 85 岩体地应力是岩体工程最基本也是最重要的工程荷载之一 因此 初始应力分布应该作为设计和分析的基础 相应地 也应当作为数值计算的初始条件 初始应力场的重要性 33 85 inisxx 5e6syy 1e7szz 2e6Inisxxs0gradgxgygz Thecomponents 11 or xx 22 or yy and 33 or zz aresettocompressivestressesof 5 0e6 1 0e7and 2 0e6 respectively throughoutthegrid TheINIcommandsetsallstressestothegivenvalues respectively throughoutthegrid butthereisnoguaranteethatthestresseswillbeinequilibrium genzonebricksize666modelelasinisxx 5e6syy 1e7szz 2e7applysxx 5e6rangex 0 10 1applysxx 5e6rangex5 96 1applysyy 1e7rangey 0 10 1applysyy 1e7rangey5 96 1applyszz 2e7rangez 0 10 1applyszz 2e7rangez5 96 1 Example3 20Initialandboundarystressesinequilibrium 34 85 赋予初始应力 genzonebricksize666modelelaspropbulk5e9shear3e9inisxx 5e6pausesolve 只用ini 所有单元将有sxx应力 但solve变为0 因为没有约束 力引起位移最后平衡 所以必须加上应力边界条件或者位移边界条件 35 85 36 85 加上应力边界条件 genzonebricksize666modelelaspropbulk5e9shear3e9inisxx 5e6applysxx 5e6rangex 0 10 1applysxx 5e6rangex5 96 1pauseplconsxxsolveplconsxx注意是约束x方向 变形参数并不影响应力场 它只影响位移场 如果考虑重力的时候 变形参数影响泊松比 从而影响侧压力系数 才会影响水平的应力场分布 37 85 表明所有单元只加上SXX 5e6Solve是平衡的 即应力不发生改变了 改变bulk shear和边界条件 solve后应力不发生变化 38 85 由于施加的初始应力场和应力边界条件产生了平衡 所以solve不产生cyc应力循环 只有不平衡力大于容许值时 才会进行cyc 设置不同的变形参数并不影响位移结果 因为系统平衡了 如果人为设置执行计算流程 如cyc10000 由于没有位移边界条件约束 模型则会慢慢随时间产生位移 但应力场是不会发生变化的 因此 intial应力得到的位移解相比较于无intial应力的位移解更符合实际地层条件 即认为是0 Cyc20000或step20000 39 85 genzonebricksize666modelelaspropbulk5e2shear3e2inisxx 5 0e6fixxrangex 0 10 1fixxrangex5 96 1pauseplconsxxsolveplconsxx注意是约束x方向 加上位移边界条件 40 85 表明所有单元只加上SXX 5e6Solve是平衡的 即应力不发生改变了 改变bulk shear和边界条件 solve后应力不发生变化 41 85 施加应力场后 如果没有应力和位移约束 模型会沿着临空面产生变形 如果施加了位移边界条件约束了模型的变形 在边界上则产生了反向的平衡力 从而产生了平衡 所以solve不产生cyc应力循环 只有不平衡力大于容许值时 才会进行cyc 设置不同的变形参数并不影响位移结果 因为系统平衡了 如果人为设置执行计算流程 如cyc10000 由于位移边界条件约束 边界上位移为0 内部则会慢慢随时间产生位移 但应力场是不会发生变化的 42 85 new Example3 24Initialstressstateforanirregularfreesurfacegenzonebricksize151510p00 0 0edge 100 0modelelasticpropshear3e8bulk5e8defmountaingp gp headloopwhilegp nullzz sqrt gp xpos gp 2 gp ypos gp 2 dz 0 06 sin 0 2 zz 100 0 SumFouriertermsfordz dz 0 06 sin 0 22 zz 20 3 quasi randomsurfacedz dz 0 04 sin 0 33 zz 33 3 topology gp zpos gp 0 5 gp zpos gp 1 0 dz gp gp next gp end loopendmountainfixxrangex 1 1fixxrangex99 9100 1fixyrangey 1 1fixyrangey99 9100 1fixzrangez 1 1 模型高度约为50m 模型长宽约为100m 43 85 仅考虑自重 算至平衡 模型高度约为50m 泊松比 0 25 0 33333模型最大高程50m左右自重来看 基底szz 2000 10 50 1 0e6sxx syy 0 333 1 0e6 0 333e6 modelelasticpropbulk5e8shear3e8setgrav0 0 10inidensity 2000solve 采用自重来计算初始应力场 44 85 仅考虑自重 算至平衡基底szz 2000 10 50 1 0e6sxx syy 0 333 1 0e6 0 333e6 45 85 propbulk5e8shear3e8 propbulk5e11shear3e11此时泊松比不变侧压力系数也不变 46 85 按照自重求初始地应力 形成的应力场表现为竖直向 自重应力 水平向 自重应力 侧压力系数k0 k0 v 1 v K0在这个方法中显然是和泊松比有关 然而土的泊松比变化不大 在0 25 0 35左右 由此求得的k0在0 3 0 6之间 和phi 在45 27度之间 在上面 泊松比是通过语句propbulk5e9shear3e9间接来定义的 这种方法会产生很大的位移 当然FLAC有一个投机取巧的方法来减小位移 就是在求初始应力的时候 把材料设为弹性 并且用很大的弹性模量确 比如上面的modelelas propbulk5e11shear3e11 修改变形参数不影响初始应力场的分布 只影响位移场的大小 在此之后 初始化位移为0 并把材料的参数恢复正常 47 85 但是这种方法存在两个问题 第一是 由此求得初始侧压力系并不和实际一致 因此初始地应力也和实际不一致 在这一点上用第二种方法可以得到比较正确的初始应力场 或者也可以通过正确的k0反算一个不切实际的泊松比v代入来求 第二也是由于设置材料为弹性的 所以是不会有破坏产生 不管初始产生的剪应力有多大 然而在材料设为弹塑性时 如果剪应力过大 材料就会在初始应力 没有其它外力作用 自动破坏 从而产生很大的变形 因此有时候第一步计算时就会产生很大的变形 很有可能就是这个原因 采用弹性材料计算的应力场在地形比较陡峭的地方有比较大的Mises应力 导致当采用弹塑性材料时破坏而不收敛 特别是有软弱夹层的地方 那么 土体 土层 在原位应力下 土体基本上处于弹性状态 因此用弹性的方法估计初始应力是否可以满足计算的要求 并不全然 在现实中土也有处于屈服状态的 象边坡的潜在滑动面应力状态 就很有可能处于殘余应力状态 为证明这一点 可以在FLAC建一个边坡模型 在线弹性模型下求得平衡 然后再把弹塑性模型换回 在不加任何外载的情况下再求解一次 我相信会有很大的位移产生 并且会出性塑性屈服点 并可以预测 边坡的安全系数越小 出现的塑性屈服点越多 那么 采用塑性的本构模型来算自重应力场 是否更为贴切呢 48 85 iniszz 1 00e6 grad0 0 0 020e6 inisxx 0 36e6 grad0 0 0 072e6 inisyy 0 36e6 grad0 0 0 072e6 solve 手动构建初始应力场 利用ini设置相同于自重形成的初始应力场自重应力场szz 1 0e6 sxx syy 0 36e6从公式上看地表szz 1 0e6 50 0 02e6 0 0e6sxx 0 36e6 50 0 072e6 0 0e6 TheINIcommandsetsallstressestothegivenvalues respectively throughoutthegrid butthereisnoguaranteethatthestresseswillbeinequilibrium 49 85 利用自重应力计算的初始应力场 szz 1 0e6 sxx syy 0 36e6 利用ini设置的初始应力场 未计算solve 系统未平衡时iniszz 1 00e6 grad0 0 0 020e6 整体规律与自重应力场相同inisxx 0 36e6 grad0 0 0 0072e6 inisyy 0 36e6 grad0 0 0 0072e6 50 85 利用ini设置的初始应力场 计算solve至系统平衡时 发现应力场降低了一个数量级 说明位移边界条件并不能使系统平衡 如何设置才能平衡呢 利用ini设置的初始应力场 未计算solve 系统未平衡时iniszz 1 00e6 grad0 0 0 020e6 整体规律与自重应力场相同inisxx 0 36e6 grad0 0 0 0072e6 inisyy 0 36e6 grad0 0 0 0072e6 51 85 以上分析可知 单纯用自重来计算应力场 耗时较长 如果仅通过ini手动设置应力场 有时并不能使系统平衡 这说明加入人为得到的初始地应initial 在弹性情况下是不会对初始应力造成影响的 那如何设置才能平衡呢 可以想像 ini是人为设置的一个单元应力 其可以加快平衡迭代的过程 在计算初始应力场的时候 如果将岩体设定为弹性的话 在给定的重力与约束条件下 初始应力是不会对最终的应力结果造成影响的 ini是模型的内力 setgrav是外力施加对模型的作用力 没有外力怎么有内力 这样做其实只是减少程序的计算步骤 这样设置以后整个模型就已经是平衡了 同时设置重力的原因不是为了平衡减少程序计算步骤 而是因为如果未初始化应力的话 在setgrav的时候 模型会在自重应力下发生应力变形 所以 初始化应力的原因是防止模型在自重应力下发生变形 一 施加自重 施加位移边界条件modelelasticpropbulk5e11shear3e11setgrav0 0 10inidensity 2000iniszz 1 00e6 grad0 0 0 020e6 inisxx 0 36e6 grad0 0 0 0072e6 inisyy 0 36e6 grad0 0 0 0072e6 fixxrangex 1 1fixxrangex99 9100 1fixyrangey 1 1fixyrangey99 9100 1fixxyzrangez 1 1Solve需527步平衡 一 施加自重 施加应力边界条件modelelasticpropbulk5e11shear3e11setgrav0 0 10inidensity 2000iniszz 1 00e6 grad0 0 0 020e6 inisxx 0 36e6 grad0 0 0 0072e6 inisyy 0 36e6 grad0 0 0 0072e6 applysxx 0 36e6 grad0 0 0 0072e6 applysyy 0 36e6 grad0 0 0 0072e6 fixxyzrangez 1 1Solve需875步平衡 52 85 一 施加自重 施加应力边界条件modelelasticpropbulk5e11shear3e11setgrav0 0 10inidensity 2000iniszz 1 00e6 grad0 0 0 020e6 inisxx 0 36e6 grad0 0 0 0072e6 inisyy 0 36e6 grad0 0 0 0072e6 applysxx 0 36e6 grad0 0 0 0072e6 applysyy 0 36e6 grad0 0 0 0072e6 fixxyzrangez 1 1Solve需875步平衡 一 施加自重 施加位移边界条件modelelasticpropbulk5e11shear3e11setgrav0 0 10inidensity 2000iniszz 1 00e6 grad0 0 0 020e6 inisxx 0 36e6 grad0 0 0 0072e6 inisyy 0 36e6 grad0 0 0 0072e6 fixxrangex 1 1fixxrangex99 9100 1fixyrangey 1 1fixyrangey99 9100 1fixxyzrangez 1 1Solve需527步平衡 53 85 在某些构造区域 很可能水平应力远大于竖直应力 此时不能用自重来计算初始应力场 此时采用ini来影响水平应力场的分布 Thereisnosimplewaytodeduceanequilibriumstressdistributionforthisgrid itmustbemodeledasaboundary valueproblem However wemayinsertinitialstressesinordertospeeduptheconvergenceandinfluencethefinalstressdistribution Forexample ifweknowthatthereisahighhorizontalin situstress withonlyasmalldecreaseinstressnearthesurface wecaninitialize xxand yyto 0 2e6attheapproximatelocationofthesurface increasingto 4 0e6atthebottom Theverticalstress zz canbesettocorrespondtotheaverageoverburden setgrav0 0 10inidensity 2000iniszz 1 00e6 grad0 0 0 020e6 inisxx 4 00e6 grad0 0 0 076e6 inisyy 4 00e6 grad0 0 0 076e6 fixxrangex 1 1fixxrangex99 9100 1fixyrangey 1 1fixyrangey99 9100 1fixxyzrangez 1 1Solve需402步平衡 inisxx 4 00e6 grad0 0 0 076e6 inisyy 4 00e6 grad0 0 0 076e6 地表sxx 4 0e6 50 0 2e6所以 0 076e6 54 85 由此可见 ini是人为设置的一个单元应力 不仅可以加快平衡迭代的过程 而且可以采用ini来影响水平应力场的分布 同时设置的ini水平初始应力是不会对竖直应力场造成影响的 一般而言这样是为了计算比较快 因为这样设置以后 此时几乎已经平衡 了 计算步较少了 变形也就很少了嘛 当然最好再在以后计算时初始化位移为0 当然还有另一种方法就是开始时设成弹性的 modelelas 或者还是mohr的但参数给的要大一些 目的是防止破坏 直接运行到平衡达来获取初始应力值 在真正计算时再改成mohr并将参数改为真实的数值 这种方法更适合于单元不太多 因为比较简单 但如果单元数目很多的话 那么此种方法速度会很慢的 55 85 1 Donotinitializestresses allowgravitytocompactthelayer 2 Initializehorizontalstressonly notverticalstress 3 Imposeconstantstressatthelateralboundariesratherthanzerohorizontaldisplacement 4 Removeirregularoverburdenfrominitialgridofuniformthickness 5 Allowplasticflowtooccur thusremovingstressconcentrations 6 Builduptheprofilelayerbylayer equilibrateeachlayer 7 初始地应力的问题 应该说是一个困难的问题 并没有一个完整的解决方案 在有限元工具中 比如ansys或者abaqus 都可以先固定边界 然后施加重力得到一个弹性应力场 然后再将这个应力场作为初始应力场导入和重力平衡 的确可以不产生位移 但关键的问题是 采用弹性材料计算的应力场在地形比较陡峭的地方有比较大的米赛斯应力 导致当采用弹塑性材料时破坏而不收敛 特别是有软弱夹层的地方 往往这些区域都不大 在靠近地表比较陡峭的地方 虽然存在这些缺点 但仍不失为一个好方法 采用ini条件加线性的初始应力 这个方法flac有 abaqus也有 但这个方法的缺点不言而喻 不合适地形起伏的情况 但最起码他可以控制Ko 这点又是前面这方法所欠缺的 所以 这两招都有其合适的场合 关键是看工程的实际作用范围了 在工程的影响范围之外的地应力影响比较小 是否可以比较粗糙一些 56 85 1 利用重力来计算竖直向应力场2 人为施加水平向应力场 加快收敛速度 同时可以得到实际的水平应力场分布3 采用塑性本构计算初始应力场 允许塑性流动发生 这样可以减小局部应力集中的现象 4 等比例提高bulk和shear参数 是不会改变泊松比的 这样在不改变应力场的前提下 也减小了自重产生的位移场 这也是自重产生的应力场在后面的计算中是要清零的原因 5 如有原位测点的地应力资料 分段计算出不同区位的梯度施加力和梯度值 由于flac以原点来计算梯度力 inisxxs0 grad0 0 gz rangzz1z2 测点应力场h1处为s1 h2处为s2从公式上看sxx h1 s1 s0 0 x1 0 y1 gz z1sxx h2 s2 s0 0 x2 0 y2 gz z2 57 85 河谷边坡中的地应力场是一个形如岩体开挖以后的二次应力场 不适宜用构造应力分量与自重应力分量的叠加 理论上 河谷发育是在构造运动格局基本确定以后的地表地质作用 模拟河谷侵蚀下切来获得局部地应力场与真实地应力场更接近 河谷地应力场模拟思想 58 85 3 河床 距离 1 驼峰应力分布 拉裂区 压致拉裂区 张剪型破裂区 Morh columb准则 3 2 2 2 1 1 3 2 1 3 应力集中区 应力 过渡区 原岩应力区 应力松弛区 河流地质作用控制区 河谷地应力场的基本特征 河谷地构造应力是一个非常不均匀的 受河谷发育改造的局部地应力场 59 85 河谷地应力场模拟的基本考虑 x 坝址区两岸为单斜构造 右岸地形呈现5级阶面 因此 将河谷侵蚀下切假设为5期 1300m 高程 60 85 假设五期河谷侵蚀下切 61 85 地应力场模拟结果正确性验证 谷底高应力区 岩芯饼化 应力松弛带 62 85 Flac3d前后处理和三维模型构建方法 主要认识和结论 FLAC3D计算结果的后处理 FLAC3D计算中的应力场 FLAC3D三维模型的快速构建 Sufer的三维出图和数据导出 目录 CAD的前处理和数据导出 63 85 Flac3d只能出云图 不能出等值线图 而工程中大多要求出等值线图 Flac3d三维截剖面出图效果不好 不能截几个剖面一起出图所以 有必要导出其计算的应力位移信息 利用sufer或tecplot出图 Sufer 优点 要求数据结构简单 可以采用白化功能进行局部区域的等值线出图 缺点 需要给定边界 bln文件 不能出三维等值线 TECPLOT 优点 不需边界bln文件 可以出三维等值线图 截剖面方便 可同时多剖面出图缺点 数据结构复杂 不能隐藏部分单元出图 64 85 DEFINEGetGPDispXP1 p gp headloopn 1 P total node id gp id p1 p message 1 string gp xpos P1 p message 1 message 1 string gp ypos P1 p message 1 message 1 string gp xdisp P1 p 1000 message 1 message 1 string gp ydisp P1 p 1000 message 1 message 1 string gp zdisp P1 p 1000 status write message 1 P1 p gp next P1 p endloopENDGetGPDispX 采用flac内嵌的fish语言进行所有节点的位移值提取 65 85 definitializationIO READ 0IO WRITE 1IO FISH 0IO ASCII 1Sufer dispX string Sufer dispX txt arraymessage 1 endinitialization status open Sufer dispX IO WRITE IO ASCII ifstatus 0thenGetGPDispXstatus closeii out SuccessfullyWriteDataIntoFile Sufer dispX elseii out OpenFileError Status string status endif 采用flac内嵌的fish语言进行所有节点的位移值提取 也可采用setlogon的方式记录位移 但是log文件中也包含了其它信息 建议采用文本通道加格式输出的方式 比较方便 66 85 67 85 originlbydynamax simwe modifiedbyHaitang Hohai 18 12 2006byHaitang Hohai Modifications 1 adaptedtoalltypes