岩土工程数值方法

刘兴隆 厉梅

摘要：逐渐发展起来的一些岩土分析手段与数学理论，如信息量法、层次分析法、随机模拟法、无网络法、数值流形法、离散元法、分形理论、可靠度分析、人工神经元网络和智能岩石力学等，已经呈现出综合应用的趋势，对于岩体力学研究而言，岩石破坏过程的渐进性、岩体内部初始损伤的存在及块体之间的不连续特征是必须考虑的因素，因此建立在连续介质力学基础上的传统有限单元法具有明显的局限性。各种新方法的涌现从不同方面推动了岩石力学数值计算方法的进步。

关键词：岩土数值模拟 有限元法 无网络伽辽金法扩展有限元法数值流形法离散元法

Abstract: gradually developed some geotechnical analysis method and mathematical theory, such as information method, the analytic hierarchy process (ahp), random simulation method, the numerical manifold method, no network, discrete element method, fractal theory, reliability analysis, artificial neural network and intelligent rock mechanics etc, has presented a comprehensive application trend, for research in rock mechanics, rock failure process of rock mass progressive, the existence of the internal initial damage and block the discontinuous characteristics between is must consider factors so based on continuum mechanics on the basis of the traditional finite element method has obvious limitation. All kinds of the emerging of the new method from different aspects promote the rock mechanics numerical calculation method of progress.

Keywords: geotechnical numerical simulation finite element method without network petro-galerkin method was expanded numerical manifold method finite element method of discrete element method

中图分类号：O241 文献标识码：A文章编号：

岩土数值模拟是否正确，其解决问题的重要基础仍然是地质工作，“地质体运动真实行为的理解比精确计算更为重要”。基础工作主要是研究地球内外动力作用的规模、机制和发展演化规律及所产生的不良地质问题，对之进行分析、评价以及提出有效的防治对策和措施。传统的岩土工程地质评价方法有自然历史分析法、工程地质对比法、模型模拟试验法、图解法等。

纯粹数值模拟方法的研究，近年来重新成为国际岩土工程界的热点。有限单元法、有限差分法等传统数值方法的使用仍然是解决实际应用性岩土工程问题的主要方法。近年来，逐渐发展起来的一些分析手段与数学理论，如信息量法、层次分析法、随机模拟法、无网络法、数值流形法、离散元法、分形理论、可靠度分析、人工神经元网络和智能岩石力学等，已经呈现出综合应用的趋势。

目前，有限元法显然成为计算力学中解决工程问题的主要数值手段，然而随着其应用范围的扩展，其固有的一些缺陷也日益突出。在金属冲压成形、高速碰撞、流固耦合等涉及特大变形的领域中，基于拉格朗日法的有限元划分网格可能产生严重的扭曲，甚至使得单元的雅可比行列式为负值，不仅在计算中需要网格重构，而且严重地影响解的精度；对高速冲击等动态问题，显式时间积分的步长取决于有限元网格的最小尺寸，因而网格的扭曲将使得时间积分步长过小，这将大幅度地增加计算工作量，降低运算效率；对裂纹的动态扩展问题，由于裂纹的扩展方向不能提前确定，在计算过程中不断地重新划分网格以模拟裂纹的动态扩展过程成为难题。由于有限元近似基于网格，因此必然难于处理与原始网格线不一致的不连续性和大变形。网格重构不仅计算费用昂贵，而且会损害计算精度。鉴于此，近些年来国际上许多著名的计算力学学者，如 T. Belytschko， O.C. Zienkiewicz，S.N. Atluri，J.T. Oden，W.K. Liu 等都对无网格方法表现出了极大的兴趣，并进行了大量的研究工作。无网格方法采用基于点的近似，可以彻底或部分地消除网格，不需要网格的初始划分和重构，不仅可以保证计算的精度，而且可以大大减小计算的难度。然而，由于目前的无网格近似一般没有解析表达式，且大都基于伽辽金原理，因此计算量很大，要超出传统的有限元法；另外，无网格近似大都是拟合，因此对于位移边界的处理比较困难，多采用拉格朗日乘子法处理。 目前已提出了十余种无网格法，其主要区别在于离散微分方程的方法（如伽辽金法、配点法、最小二乘法、彼得洛夫-伽辽金法等）和建立近似函数的方法（移动最小二乘近似、核近似、重构核质点近似、单位分解法、hp云团法、径向基函数法、点插值法等）。

众多岩土力学数值方法中，无网络伽辽金法(EFGM) 。EFGM 的特点是采用移动最小二乘法构造位移函数，从弱变分原理形成控制方程，利用罚函数满足边界条件。因此，无网络法只需结点信息，勿需大量的划割网格的前处理工作，尤其适用于岩体的裂隙扩展、大变形和动力学分析等。由于是任意点的加权拟合面，导致边界条件的处理仍然得不到很好的解决。目前，较之于有限元方法，无网格方法的问题一点不比它的优点少。首先无网格方法确实提高了插值的连续性，但同时也极大地提高了计算量，求解域内的每个点的插值系数都需要对相应节点耦合矩阵求逆，光滑性越高耦合矩阵的阶数就越高，相应的求逆就越耗费计算量，较之于有限元方法，这一部分计算负担实在太重。对于被奉为经典的二维裂缝扩展问题，仔细考察就会发现裂缝扩展问题的关键技术在于加强函数（enrichment）的采用，其实用不用无网格插值都可以，用有限元插值也没有任何问题，这就是后来发展的，也是现在最热门的扩展有限元方法（XFEM）。

1999年，以美国西北大学Belytschko教授为代表的研究组首先提出扩展有限元思想，2000年，他们正式使用扩展有限元法(XFEM)这一术语。XFEM是迄今为止求解不连续力学问题最有效的数值方法，它在标准有限元框架内研究问题，保留常规有限元的所有优点，但并不需要对结构内存在的几何或物理界面进行网格剖分。XFEM与常规有限元的最根本区别在于所使用的网格与结构内部的几何或物理界面无关，从而克服了在诸如裂纹尖端等高应力和变形集中区进行高密度网格剖分所带来的困难，当模拟裂纹扩一展时也无需对网格进行重新剖分。例如在处理裂纹问题时，XFEM包括以下三方面内容：第l不考虑结构的任何內部细节(例如材料特性的变化和/或内部儿何的跳跃)按照结构的少卧丁外形尺寸生成有限元网格;第2，采用其它方法(如水平集法)确定裂纹的实际位置，跟踪裂纹的生长;第3，借助于对所研究问题解的已有知识(不必知道封闭形式解)。改进影响区内单几的形状函数，以反映裂纹的存在和生长。由于改进的形状函数在单元内部具有“单位分解”特性，扩展有限单元的刚度知阵具有与常规有限单元样的优点，即对称、稀疏且带状。

1996年Melenk和Babuska及Duarte和Oden先后提出了单位分解法(PUM)，其基本思想是任意函数都可以用域内组局部函数表示，即

其中，为有限单元形状函数，它形成个单位分解

基于此，可以对有限元形状函数根据需要进行改进。

1997年，Babuska、Menlenk证明了PUM的收敛性并将之应用于求解高波数Helmholtz方程。PUM容许在相容的试探空间中增加用户定义的局部特性，因而对CFEM无法求得或求解代价太大的问题，可体现出PUM的独特优势。其特征为:(1)PUM容许在试探空间中包含对微分方程的先验知识;(2)利用PUM能很容易地构造出任何期望的试探空间，因而可以获得适用于高阶微分方程变分形式的试探空间。

当前最为热门的另外一种数值方法，数值流形法(MM) 采用有限覆盖体系，系统控制方程的建立依据最小势能原理，特别适合模拟断续介质材料的变形，如岩体开挖导致的大位移。该方法采用连续和非连续的覆盖函数统一解决有限元、DDA和解析方法的计算问题, 最早由石根华博士在 1991 年提出，并率先应用在块体与节理岩体的变形模拟中。经过近十年的发展，二维数值流形方法已经拥有一套比较完善的理论，但在工程方面的应用则不多见。主要的困难来自对双重覆盖的理解，覆盖函数的选取以及程序的编制。

数值流形方法采用数学覆盖和物理覆盖两套覆盖系统来定义计算区域：数学覆盖定义近似求解精度，可由用户自行选择，如规则的格子、有限元的网格或级数的收敛域，都可以转化为有限数学覆盖；物理覆盖系统是由数学覆盖和物理网格共同形成，物理网格由材料本身决定，定义其积分区域，包括材料体的边界、裂隙、块体和不同材料区域的交界面等。如果物理网格将一个数学覆盖分成两个或更多的完全不连续的区域，这些区域定义为物理覆盖。所以说物理网格将数学覆盖细分成物理覆盖系统。如图 1 所示, 分析一个包含一条不贯穿的裂隙的多边形块体, 采用三角形的数学覆盖网格将分析区域全部覆盖, 共有 6个数学覆盖 V1, V2, V3,V4,V5, V6, 每个数学覆盖由含有该节点的所有三角形单元形成。例如，节点 1的数学覆盖由单元 123组成，节点 2的数学覆盖由单元 123，253，245 共同组成。物理网格（图 1 中的粗线条）将数学覆盖V1 分成一个物理覆盖11，将数学覆盖 V2 分成两个物理覆盖 21和 22，21和 22可有不同的材料性质，沿交界边（即图中裂隙）可以有不连续的力学行为，可以拉伸，扭曲甚至完全裂开。也就是说，原始的有限元单元 245被裂隙分割成 2个流形单元21425 和 224151，它们在接触边上的节点是不同的，在接触处可以有不连续的位移，这是与有限元分析的最大的不同，也是该方法的优势所在。而在边 23 上，两边的物理覆盖都是相同的，如物理覆盖 21和 21，物理覆盖 31和 31。这表明，沿着 边23，两边的材料体是相同的，所有的力学行为是连续的，这又比不连续变形分析（DDA）方法中采用的所有的边都是不连续的假设要更符合材料本身的状况。可以说数值流形方法是利用覆盖技术将有限元，DDA和解析方法统一到内部的新的数值方法，而有限元和DDA只是数值流形方法的两种特殊情况.

图1 有限覆盖系统中流形单元的节点编号

流形元法和无单元法的出发点在于计算步骤的简化和较粗的积分,但精度并不降低。不用单元,只用结点使得它对开裂问题极为有效。

离散元法 (Distinct Element Method)是Cundall于20世纪70年代初所提出的,最初它的研究对象主要是岩石等非连续介质的力学行为。Cundall提出了第一个实用的离散元模型,并用它来模拟岩石块体的渐进运动过程。后来,Cundall和Strack提出用于模拟颗粒体的二维程序,得到与动光弹实验极为吻合的结果。至今为止,经过许多学者的共同努力,离散元法已经得到了长足的发展。早期的离散元法只能处理离散刚度块体系统,后来该方法被扩充了,可用于模拟变形块体。1980年美国ITASCA咨询集团开发离散元法程序UDEC并投放到市场。Lorig和Brady开发出离散元—边界元耦合计算程序。Cundall等开发了用于模拟节理岩体的三维离散元程序(3DEC)。离散元在我国起步较晚,但是发展迅速。

离散元法是专门用来解决不连续介质问题的数值模拟方法。该方法把节理岩体视为由离散的岩块和岩块间的节理面所组成,允许岩块平移、转动和变形,而节理面可被压缩、分离或滑动。因此,岩体被看作一种不连续的离散介质。其内部可存在大位移、旋转和滑动乃至块体的分离,从而可以较真实地模拟节理岩体中的非线性大变形特征。离散元法的一般求解过程为:将求解空间离散为离散元单元阵,并根据实际问题用合理的连接元件将相邻两单元连接起来;单元间相对位移是基本变量,由力与相对位移的关系可得到两单元间法向和切向的作用力;对单元在各个方向上与其它单元间的作用力以及其它物理场对单元作用所引起的外力求合力和合力矩,根据牛顿运动第二定律可以求得单元的加速度;对其进行时间积分,进而得到单元的速度和位移。从而得到所有单元在任意时刻的速度、加速度、角速度、线位移和转角等物理量。

综上所述：对于岩体力学研究而言，岩石破坏过程的渐进性、岩体内部初始损伤的存在及块体之间的不连续特征是必须考虑的因素，因此建立在连续介质力学基础上的传统有限单元法具有明显的局限性。以上各種新方法的涌现从不同方面推动了岩石力学数值计算方法的进步，具有相当重要的意义，但任何一种数值方法均不能解决所有的岩石力学问题，需要从地质分析的角度出发，根据岩土工程失效的内在机制选择合适的计算方法。

注：文章内所有公式及图表请以PDF形式查看。