数值模拟

《工程地质数值模拟的理论与方法》课程学习报告解决岩土力学问题的方法主要有实验方法、理论分析方法和数值模拟方法三大类。由于经典力学对解析解的设定非常严格，使其用于岩体工程的求解极为有限。数值模拟是解决岩土工程问题较有效的手段，已被学术界和工程界广泛接受。近30年来，随着岩土工程领域的不断扩展与延伸，数值分析方法得到了长足的发展，分析模式不断改进，分析精度不断提高，出现了有限差分法(FDM)、有限单元法(FEM)、边界元法(BEM)、无限元法(IEM)、拉格朗日元法(LEM)、刚体弹簧模型或刚性有限元法(RBSM或RFEM)、离散元法(DEM)、非连续变形分析方法(DDA)、无单元法(EFM)、流形元法(MM)及其耦合的数值计算方法和以数值模拟为主的渐进破坏模型等各种数值分析技术。但岩土工程材料是一种复杂的地质材料，具有高度非连续性、非均匀性以及各向异性等鲜明的地质特点，在力学性质上表现出强烈的非线性、非弹性等力学行为。因此，只有不断改进数值分析的思路与方法，使其能更好的模拟实际岩土材料的工程力学性质，才能更好的指异岩土工程实践。一、对有限单元法与拉格朗日法的认识有限单元法和拉格朗日法都是广泛应用到实际工程中行之有效的数值模拟方法。两者都被业内专家所推崇，在求解工程技术领域问题上得到了很大的收获。1．有限单元法有限元法主要分析岩土介质的连续小变形和小位移特性。有限元法在连续性分析方面取得了很大的成功，但在解决处理问题、应力与应变解答不连续问题和进行任意路径开裂计算等方面还存在着一些局限。为了充分考虑岩土介质的非连续性、非均匀性和多相性等物理特性，必须对这些连续变形分析方法进行深入的改进和发展。1.1有限单元法的基本原理弹性理论需要解一系列的偏微分方程组，并满足边界条件。常用的解法是变分法，即泛函求极值问题，这里泛函表述的是功(能量)，变分法是求功(能量)的极值。数学上表述，有限单元法就是求结构的泛函极值问题。工程界习惯用熟知的结构力学位移法建立有限元解法，简明清晰。理解清楚对结构的离散化、形函数以及通过平衡条件建立力与位移的关系。1.2有限单元法的分析方法（1）结构的离散化，使连续介质无限多个自由度化为有限个自由度；（2）选用合适的形函数位移模式；（3）建立结点力与结点位移之间的关系，形成单元刚度矩阵；（4）通过结点力的平衡建立总体刚度矩阵方程组，并依约束条件进行处理(消除矩阵的奇异性)；（5）形成荷载矩阵；（6）解方程组，求位移场、应力场。在求解过程中建立刚度矩阵(几何条件明确)的关键依赖于本构方程，它决定于解题的性质、水平与可信度。1.3有限单元法的应用应用有限单元法能较好地反映复杂的岩体介质，可以解弹性、塑性和粘性问题，可以解动载、静载问题和流体问题，也便于施工程序的模拟。目前对岩体工程至关重要的结构面模拟问题——节理元(固体或流体)，在理论上尚欠成熟，无厚度元的质量为零，使动力问题求解困难。计算中多组结构面相交时，剖分处理要非常谨慎，它可能使岩体相互嵌入或分离，甚至求解失败。解弹塑性问题时采用增量理论，是变刚度法。用迭代法求解时，可以用常刚度法或变刚度法，后者较前者收敛快，由于解题中源程序处理的差异，常刚度法省机时，且稳定性好。用有限元法解题前处理工作量很大，特别是复杂岩基的三维问题剖分相当繁杂，很难实现自动剖分。应事先根据结构特点、工作状态、精度要求周密筹划单元的剖分。打印成果的阅读、分析处理等，需通过学习、实践、体验，才能提高工作质量和效率。有限单元法是近似解法，单元剖分的疏密程度与质量、效益密切相关，在理论上要把握好这个度且保证收敛。2．拉格朗日分析法以连续介质大变形分析为目的的拉格朗日元法在实际工程中也得到了较好的应用。基于拉格朗日元理论，美国的ITASCA咨询集团于20世纪80年代编写的专用程序FLAC，具有界面友好，计算能力强，使用简单等特点，并附有功能强大的前后处理器，使计算模拟过程直观明了。拉格朗日元法可以同时考虑岩土体的材料非线性和几何非线性，并能跟踪物体变形的全过程，适于分析岩土力学中的大变形问题。这种方法避免了有限元法进行大型矩阵的复杂计算，但时间步长的选择成了一个新的突出问题，时步过大会导致解答的不稳定，时步太小则会使计算时间过长。2.1拉格朗日的基本原理在确定研究区域后将其离散划分为一系列的网格；当节点受到作用后，每一时步中，外力或速度的变化和时间步长求得节点的应变率，通过本构方程求得应力张量，从而求得节点的不平衡力；再由节点的不平衡力求得节点的速度和位移，并开始下一步的迭代，直到不平衡力足够小或位移趋近平衡为止。2.2拉格朗日的分析方法拉格朗日元法运用流体力学中跟踪质点运动的物质描述方法，即拉格朗日拖带坐标系方法，利用差分格式，按显示时步积分方法进行迭代求解，根据构形的变化不断更新坐标系，以此模拟岩土介质的有限变形和大位移行为。（1）模型的建立，建立网格，定义材料的本构模型及物理特性，指定边界与初始条件等。（2）确定模型的平衡状态，在给定的条件下，通过检测最大不平衡力，节点速度和位移等，确定模型在何种情况下达到平衡。同时检查模型的反应，当模型不能达到平衡时，返回上一步重新建立模型。（3）改变模型的设置，当进行计算，改变节点的应力，约束条件以及材料的特性等，当结果满意时，结束计算过程，并对结果进行必要的分析和解释。2.3拉格朗日的应用拉格朗日法能像有限元法那样适用于多种材料模式和边界条件的非规则区域的连续问题求解，能够处理一般的大变性问题，能模拟岩体沿某一弱面或结构面产生的滑动变形，还能针对不同的材料特性使用相应的本构方程来反映实际材料的动态行为。采用这种方法可以模拟弹性模型、摩尔-库仑模型、遍布节理模型、应变硬化和软化模型等。现已广泛应用于边坡、基础、坝体、隧道、混凝土衬砌、锚杆和锚索设置、地下采场和酮室等岩土工程分析中。二、对数值模拟课程的几点建议数值分析法在岩体力学中的应用越来越广，各种方法互有长短，且本身也处在不断的发展之中。随着岩土工程领域的不断扩展与延伸，数值分析法也得到了长足的发展，分析模式不断改进，分析精度不断提高。岩土工程材料是一种复杂的地质材料，具有高度非连续性、非均匀性以及各向异性等地质特点，在力学性质上表现出强烈的非线性、非弹性等力学行为。对实际的工程问题，精确的分析对实际问题的性质、程度、时效、经济合理等诸多因素的选定提供了很好的理论依据。因此，在这个电信时代学习数值模拟课程是很有必要的，也很重要的。我是勘查技术与工程的学生，由于专业的原因，在这方面了解得很少，只是听土木的同学说过这些方法，觉得他的功能很强大，可以直接计算成图给人一个直观的图形，通过专业的分析就可以解决很多的实际问题。但是对于其基本原理及具体应用没有了解。由于现在计算机的广泛应用，任何领域中的很多问题都可以通过计算机软件来高效地解决。为了能更好的在专业领域里学习，本学期学习了这门课程。为期三个月的学习结束后，对几种比较常用的方法如有限元法、离散法、快速拉格朗日法有了一个大致的了解；对其各自的基本原理及相关的实际工程应用也有了一些了解；对基本的上机模拟过程中有一个浅显的认识；对实际应用中应该注意的问题，相关问题的解决方案等也有点了解……但是这些认识仅仅是基于理论上的，对于种种应用软件，要很好地掌握并应用到自己的专业和科研项目中去，还要不断的在实际应用中获取知识，通过不断的使用，发现问题，分析问题，解决问题来达到对软件的掌握。这门课程的开设是很有必要的，它是一门紧密结合专业知识，结合实践问题的应用学课程，但是我个人认为在教学研究生中有些可以做一些适当的调整可得到更好的效果。室内教学中可以直入主题，少一些方法的发展的介绍，先将一种模拟方法的原理介绍清楚后，结合工程实例介绍该种方法在其中的具体应用过程，最后所分析得到的结果（这一点很重要，会使用软件但是不知道所得出的数据、图形所说明的问题，只能是做无用功）。上机实习课我认为是很有必要的，而且应该根据课时安排安排好上机时间。在每次教授完一种数值模拟的方法和介绍相关实例后就应该安排一次上机实习，让同学在进行实际的动手操作中自己亲身体验一下如何将实际问题通过软件来呈现，遇到问题可以适时地与老师沟通得到解决，这样能更有效的、更快速的掌握数值模拟方法。尽管有部分同学已经掌握了这种方法，但是大部分同学还是比较生疏或不熟练的。本学期上机的课时很少，实习内容也相当简单，但对于没有接触过这种软件的同学来说，是一个很好的学习机会，使他们能够简单了解这种程序。此外，老师无需将同一个实习内容上机几次，对于宝贵的实习课程，希望每一次实习课程都有不同之处，这样才能够较全面的了解这类方法；同时适当将实习报告的难度加大，如结合各自专业领域内的工程实例来进行模拟，并分析结果，交上一份完整的报告。尽管这样老师教学难度加大，同学的学习负担加重，但是学生能够更深入地了解这种方法，也为今后学习工作打下坚实的基础。三、上机实习报告计算模型分别如图1、2、3所示，边坡倾角分别为30°、45°、60°，岩土体参数为：密度ρ＝25kN/m3,弹性模量E＝1×108Pa，泊松比μ＝0.32，抗拉强度σt＝0.8×106Pa，内聚力C＝4.2×104Pa，摩擦角φ＝17°，膨胀角Δ＝20°采用FLAC3D软件建立单位厚度的计算模型，并进行网格剖分，参数赋值，设定合理的边界条件，利用FLAC3D软件分别计算不同坡角情况下边坡的稳定性，并进行结果分析。附换算公式：1kN/m3=100kg/m3剪切弹性模量：)1(2EG体积弹性模量：)21(3EK（一）倾角为30°的边坡其模型如图1所示。首先根据如下程序，对模型进行网格剖分：genzonebrick&p0000p110000p2020p30040&size50110genzonebrick&p040040p1100040p240240p374.64060&p4100240p574.64260p6100060p7100260&size3011图1倾角为30°的边坡(单位：m)网格剖分图如下：FLAC3D2.10ItascaConsultingGroup,Inc.Minneapolis,MNUSAStep39370ModelPerspective22:06:35TueMay272008Center:X:5.000e+001Y:1.000e+000Z:3.000e+001Rotation:X:0.000Y:0.000Z:0.000Dist:2.775e+002Mag.:1Ang.:22.500SurfaceMagfac=0.000e+000图1-a网格剖分图然后对稳定性计算等进行编程，整个程序如下：newsetlogonsetlogfilestability.loggenzonebrick&p0000p110000p2020p30040&size50110genzonebrick&p040040p1100040p240240p374.64060&p4100240p574.64260p6100060p7100260&size30110modelmohrpropdensity=2500.0bulk=1.0E8shear=3.0E7coh42000.0tens1.0E10&friction=17dilation=20fixxyzrangez-0.10.1fixxrangex99.9100.1fixxrangex-0.10.1fixy30°100406025.3640setgravity=10.0plotaddaxeredplotconzdis;defcalfosminf=0.6maxf=2.0loopwhilemaxf-minf0.01fs=(maxf+minf)/2.0refric=atan(0.31/fs)*180/3.14recoh=42000/fscommandinisxx0.0syy0.0szz0.0sxy0.0sxz0.0szz0.0inixvel0.0yvel0.0zvel0.0inixdis0.0ydis0.0zdis0.0profri