基于水平集坐标的二维节理动态扩展过程无网格法模拟.doc

第XX卷 第X期 庄晓莹等 7 基于水平集坐标的二维节理动态扩展过程无网格法模拟庄晓莹1,2，黄润秋3，朱合华1,2(1. 同济大学 岩土及地下工程教育部重点实验室, 上海 200092；2. 同济大学 地下建筑与工程系, 上海 200092；3. 成都理工大学 地质灾害防治与地质环境保护国家重点实验室, 成都）摘要：本文发展了模拟岩体节理张开和闭合力学行为的无网格法近似函数构造方法，采用水平集坐标描述二维节理几何信息，沿任意几何形状节理两侧布置跳跃结点，节理扩展时更新水平集方程捕捉节理几何形状。在节理张开时引入断裂线，采用衍射准则构造无网格近似差值函数。当节理闭合时，发展了类似于Goodman接触摩擦单元的无网格函数构造方法，考虑节理面上下两侧的相对滑动和法向接触。与Goodman单元不同之处在于，本文构造方法基于材料点而非结点，对于几何形状复杂的节理在模型生成和计算上更灵活简便。最后，应用本文方法对二维节理受拉、节理受压的静态和动态破裂扩展过程做了分析，并通过对照解验证了本文方法的正确性。关键词：岩石力学；节理扩展；无网格法；节理闭合；水平集；接触摩擦；中图分类号：O 319.56 文献标识码：A 文章编号：10006915(2005)01000102MODELLING 2D JOINT PROPAGTION IN ROCK USING THE MESHLESS METHODS AND LEVEL SETSZHUANG Xiao-ying1,2 HUANG Run-qiu3 ZHU He-hua1,2(1. Key Laboratory of Geotechnical Engineering of China Ministry of Education, Tongji University, Shanghai 200092, China;2. Department of Geotechnical Engineering, Tongji University, Shanghai 200092, China;3. State Key Laboratory of Geohazard Prevention and Geoenvironment Protection, Chengdu University of Technology, Chengdu, China)Abstract：In this paper, the meshless method, element-free Galerkin method is used to model 2D crack propagation in jointed rock mass where the level set method is used to describe the crack geometry. To capture the displacement jump across the joint when the join is open, the diffraction rule is used to modify the weight function, and hence construction of the discontinuous meshless approximation. The joint closure is modeled by a new meshless approximation, which is similar to the well known Goodman element in the family of finite elements to model the sliding and contact between joint surfaces. Different from Goodman element, the construction of stiffness matrix is based on material points instead of nodes and thereby is more flexible to deal with joint of complex geometry. Test examples are tested to model the crack propagation in rock samples subjected to compression and tension showing the correctness of the method.Key words：rock mechanics; joint propagation; meshless methods; joint closure; level sets; contact and friction;1 引 言岩体稳定性分析一直是岩土工程领域的重要研究课题，而其中最关键的问题之一就是节理面扩展贯通的机理和模拟方法的研究，因为节理面的存在和发展控制着岩体的整体强度、变形和稳定性，准确地预测和分析节理面的力学行为是岩体稳定性评估和工程治理的重要依据。针对节理岩体破坏机理和等效力学模型，国内外学者从实验和解析角度出发做了大量研究，如Lajtai 1根据剪切实验提出了岩体剪切破坏、张拉破坏和挤压破坏模式和节理岩体强度准则，Oslen2讨论了节理岩体脆性断裂模式与边界条件关系，李术才等3研究了岩石断裂面压剪脆性破坏条件, 周小平等4采用 Dugdale-Barenblatt 模型分析了断续节理的贯通机制，并提出了节理扩展准则。Li等5通过实验方法对次生节理和翼节理的生成和扩展规律做了研究，认为节理扩展路径通常按折线形式发展。刘明远6等认为节理扩展贯通需考虑节理的粗糙度和起伏度。上述这些节理岩体破坏模型与数值模拟方法相结合，可实现破坏过程的动态分析7；同时室内试验结果也对验证理论和数值模拟方法的正确性提供了依据。当岩体受到大量不同走向和分布节理切割呈碎屑，可视为连续松散介质，采用连续介质分析方法8-10，如有限元法、有限差分法、边界元法、无网格法和扩展有限单元法等；若岩体呈块体分布，需考虑机动平衡条件，则采用非连介质分析方法，如离散单元法（DEM）11、流形方法（NMM）12和非连续变形分析（DDA）13。对介于两者之间，岩体主体部分较完整，存在若干个与对分析对象尺度相近的控制结构面，此时节理分布和扩展对岩体的稳定性起控制作用，可采用引入不连续界面的连续介质力学分析方法如无网格法，许多等岩体工程如水电坝坡等正属于此类问题。这些节理面在开挖卸荷、降雨等外界因素作用下会扩展乃至贯通。现有节理岩体稳定性分析数值方法，主要集中在块体分析13、节理扩展9,14、强度折减法15、塑性极限分析16等。节理扩展模拟研究主要集中在非连续分析方法、无网格法和扩展有限元法这三种方法中，如焦玉勇等9研究了模拟节理岩石试件动态扩展的二维数值流形方法，Zhu等17应用无网格法Shepard方法分析了节理岩石边坡稳定性。采用基于连续介质力学的数值方法模拟节理开裂扩展过程最困难的问题就是裂缝前沿计算网格更新和对含裂纹奇异问题的求解精度和稳定性。起源上世纪70 年代，快速发展于90 年代的无网格法，提供了克服此类问题的可能。它摆脱了限元法在处理各种大变形、动边界问题的由网格存在引起的模拟困难，尤其适用于断裂模拟8-9。本文将结合无网格法模拟动边界问题的优势，采用水平集坐标描述和捕捉任意二维节理几何形状，研究二维节理动态扩展过程的无网格模拟方法。2 二维断裂问题的无网格法模拟2.1无网格近似函数构造无网格法是指对无需事先定义计算网格，通过离散方式求解偏微分方程的一类数值方法的统称8-9。其中无网格伽辽金法(EFGM)因其数值稳定性高、求解精度高、收敛快成为目前使用最为广泛的无网格法。本文采用此方法来分析二维节理动态扩展问题，但文中所提出算法也可同样用于其他无网格法，如局部彼得罗夫-伽辽金法(MLPG)和无网格最小二乘Shepard法(MSLS)。在无网格伽辽金法中位移插值函数采用滑动最小二乘法进行构造，对于分析域内任一点，通常积分点或材料点，其位移表达如下8： (1)其中为多项式基函数组成的列向量，m的项数。对二维问题，取线性基函数时时，m=2，其它高次基函数见8。 (2) (3) (4) 为nn的对角阵，其对角元素依次为各结点权函数值，本文采用如下径向权函数： (5)其中，为到结点I的距离，为结点I的影响域半径，为结点坐标向量。 是mn矩阵且每列元素为8。为长度为n的形函数向量，其中I为结点编号，n为包含在影响域内的结点总数。 2.2 捕捉节理位移跳跃的权函数在无网格法中，节理、空洞、裂隙等不连续面可通过修改权函数影响域来确定，如可视准则和衍射准则。可视准则计算具有计算简便、代价小的特点，但缺点是存在裂纹尖端虚假扩展，因此计算精度降。相比之下，衍射准则不存在此问题18，衍射准则计算较为复杂，尤其在三维断裂模拟中，衍射路径的计算代价较大。但在二维分析中，其计算代价可以接受，本文中采用此方法。图1 节理两侧影响域采用衍射准则Fig1. The diffraction rule for modified nodal support along either sides of the joint衍射准则原理如图1所示，若某结点影响域被节理线切割，则在计算该结点到某材料的距离时采用如下方程： (6)其中，为衍射前结点I和材料点的距离，和分别为结点I和材料点到节理尖端的最短距离，如图1虚线所示。需注意随变化。衍射准则如同模拟光线绕过障碍物，由此反映节理或裂缝存在引起的两侧变形及应力不连续性。3 二维节理的水平集坐标描述3.1 节理几何形状的水平集坐标描述在水平集法中，目标界面 (interface) 的几何形状由距离方程(signed distance function) 隐式地描述19，即不显式地离散对象本身，而采用在记录域内离散点到界面的最短距离。水平集法使用两组正交方程描述断裂面，例如图1 所示折线节理面（黑），水平集函数（红色虚线）和（蓝色）分别捕捉空间任意一点到节理面前端的法向和切向最短距离，节理面本身是和的交集18。分析域内任意一点到节理的法向最短距离为，切向最短距离为，节理端部的最短距离为： (7)在节理扩展过程中，不断更新求解如下的水平集的 Hamilton-Jacobi方程18： (8)其中为节理沿梯度方向的扩展速度，由材料破坏准则决定，如金属材料疲劳破坏时，可根据Paris疲劳准则计算，对于类似岩石、混凝土等的脆性材料破坏时，可设成与单位计算时间步内的扩展长度，需要通过迭代试算后确定其合适大小18。和求解方法类似。图2 二维折线型节理的水平集描述Fig. 2 The level sets representation of the 2D kinked crack3.2基于水平集坐标的节理两侧结点布置在无网格计算中需要采用水平集坐标在节理上下两侧布置结点，可更好地模拟由节理引起的位移不连续性，提高计算精度19。图3 节理两侧布置结点Fig. 3 Offset nodes along either side of the joint材料点如图1中的所示。需要对于无厚节理可在节理线上下两侧偏移一较小距离，本文计算取 = 0.01Lk，其中Lk为节理迹长或特征长度。对于直线型节理Lk为节理长度，对于折线型节理则取水平集坐标点邻近段折线长度。对于有厚度的填充型节理，可在无网格计算模型中先采用水平集坐标生成无厚节理线，再分别向断裂线两侧偏移水平集坐标均匀地布置结点，表示节理厚度。3.3 节理两侧位移跳跃项的水平集判断在计算式（6）衍射准则前，需先判定材料点、结点和节理之间的空间位置关系。如图4所示，虽然材料点位和结点分别位于节理上下两侧，且在I的影响域内，但由于点落在I与节理端部连线的阴影区内，因此他们之间不存在变形跳跃，无需通过衍射准则计算，反之若落在阴影区外则应按式(6)计算。由图4，通过几何运算可知当 (9)时，点和I需考虑变形不连续。其中，和分别为结点I的法向和切向水平集坐标。图4 采用水平集坐标判定节理变形不连续性Fig.4 The displacement jump based on level sets3.4 节理开裂及扩展准则本文采用岩石抗拉破坏准则和压剪破坏准则，当某点主应力值超过岩石的抗拉强度则破坏，即式认为裂纹扩展。节理扩展方向按断裂力学分析，因此需要在节理端部计算应力强度因子，这里采用J积分方法节理，计算方法详见20。节理扩展可能是多种断裂模式的组合，扩展方向角由如下准则确定21 (10)由上式可推得 (11)其中和分为模式I和II应力强度因子，为节理扩展方向角，需要注意的是 的方向和大小均相对于节理局部坐标系定义，如图5所示。4 基于材料点的节理闭合模拟节理受压时采用建立如图5所示局部坐标构造无网格近似函数。图5 二维节理闭合行为的无网格计算模型Fig. 5 The meshless model for 2D joints设节理上下侧相对位移向量为，可表达为 (12)其中，和表示节理上下侧材料点坐标， (13)为节理上侧材料点位移形函数 (14)其中M为包含材料点在影响域内的结点数目。类似地可导出节理下侧位移为， (15)其中 (16)N表示含在影响域内的结点数目。将式（13）和（14）代入（12），经过重写和整理，可得节理两侧相对位移表达式为 (17)其中 (18) (19)节理切向和法向应力可表示为： (20)其中，Dc是节理材料矩阵， (21)和分别为节理法向和切向刚度。于是节理能量泛函表达为 (22)如图所示由于节理的相对位移定义在局部坐标系和上，而结点列向量和位移插值函数式（17）定义在全局坐标上，因此需要将映射到局部坐标系上。定义转换矩阵： (23)节理局部坐标可通过水平集坐标的梯度确定：， (24)则节理上下两侧相对位移为 (25)对变分，并将式（25）代入（22），可得 (26)于是可得包含节理于影响域内的结点刚度矩阵为 (27)其中，是 (M+N)(M+N)的方阵且包含了所有含材料点和的在影响域内结点。5 测试算例5.1单一节理张开扩展首先，验证本文方法在模拟节理受压的正确性，考虑如图所示的二维裂纹试件，高度h=18，宽度w=9，裂缝长度为2a = 2，节理倾角45，试件两边受均匀拉伸力初始值为1KPa，模拟节理扩展时通过荷载控制法增加拉力。计算参数见表1，不考虑节理抗拉强度。无网格模型采用3719规则结点排布，节理两侧按2.3节方法布置结点。根据脆性断裂力学解析解，当试件尺寸为无穷大，即考虑为远端应力场时，裂纹两尖端应力强度因子为21： (28) (29)其中，a为节理半长。本算例试件尺寸相对节理较大，可近似认为接近于解析解。节理的水平集坐标如图7(a)所示，根据还可方便地计算出计算域内任意点到节理端部的最短距离和夹角，可用于计算应力强度因子如图7(b)。图6 中心斜节理受拉试件几何形状Fig. 6 The geometry of slanted centre joint under tension(a) 法向切向水平集 (b) 径向水平集图7 节理几何信息的水平集坐标描述Fig. 7 The level sets description of joint geometry在本文算例中，KI解析解0.8862，KII解析解为0.8862，本文方法所得结果为KI = 0.8987，KII =0.9013，误差分别为1.41%和1.70%，计算结果精度良好。按脆性断裂假定计算，扩展步长为节理初始长度的10%。第一步节理上端扩展方向角相对水平方向夹角为-8.1度，下端节理为171.1度，如图8所示。8。到了第二步，扩展方向角相对水平方向夹角为，第三步后节理基本接近水平。节理在扩展10步后的路径如图8所示，符合文献22的参考结果。(a) 上端翼裂纹 (a) 下端翼裂纹图8 节理扩展第1步时出现的翼裂纹Fig.8 The wing crack occurred in the first step propagation表1节理试件材料参数Table 1 Material parameters of jointed rock sample弹性模量（GPa）泊松比密度（kg/m3）节理抗拉强度(GPa)16.50.2526500.07图9 受拉节理扩展10步后计算结果Fig.9 The crack geometry of joint under tension after propagation for 10 steps5.2单一节理受压扩展图10 中心节理受压试件几何形状Fig. 10 The geometry of slanted centre joint under compression考虑如图10所示受压斜节理，其倾角为30度。为验证对比计算结果，计算模型尺寸与文献12相同。节理初始长度2a=2，计算材料参数见表1。节理受压时，为防止接触贯入，设置成较大数值，这里采用20倍的E，设置为0.1E。计算时通过荷载控制法迭代，当节理张开时按衍射准则处理，闭合时按式（27）计算。初始状态时，在初始步时，KII = -0.8661 , 节理扩展方向角为35.2，相对于水平方向夹角为9.8，扩展16步后节理路径如图12所示，与文献12扩展形状相吻合良好。 Fig. 12受拉节理扩展16步后计算结果Fig.12 The crack geometry after propagation for 10 steps5 结 语本文发展了基于无网格法和水平集坐标的节理扩展模拟方法。节理受拉时通过修正衍射权函数引入节理两侧位移不连续变形项；节理受压时通过材料点建立上下两侧节理面的接触摩擦关系。上述两种无网格近似函数构造方法，在节理附近区域均采用了水平集坐标，对处理任意几何形状的节理具有通用性。测试算例结果表明，本文方法对受拉压两种状态的节理试件的模拟的准确性和有效性。水平集法的内存占用和计算时间较长。下一步，作者希望通过建立局部水平集法减少计算代价。如文献所指出5的节理面起伏度对节理扩展有影响，拟在下一步研究中考虑。 参考文献(References)：1 Lajtai E Z. Shear strength of weakness planes in rock J. International Journal of Rock Mechanics and Mining Sciences and Geomechanics 1969, 6(5): 499-515.2 Olson J E. Joint pattern development: effects of subcritical crack growth and mechanical crack interaction J. Journal of Geophysical Research, 1993, 98 (B7): 12251-12265.3 李术才, 陈卫忠, 朱维申, 王书法.加锚节理岩体裂纹扩展失稳的突变模型研究J. 岩石力学与工程学报, 2003, 22(10): 1661-1666. (LI Shucai, CHEN Weizhong, ZHU Weishen, WANG Shufa. Catastrophe model of crack propagating in bolt-supported jointed rock mass J. Chinese Journal of Rock Mechanics and Engineering, 2003, 22(10):1661-1666. (in Chinese)4 周小平, 王建华, 哈秋舲.压剪应力作用下断续节理岩体的破坏分析J. 岩石力学与工程学报, 2003, 22(): 1437-1440. (ZHOU Xiaoying, WANG Jianhua, HA Qiuling. Failure analysis 0f intermittent jointed rocks under compression-shear stress J. Chinese Journal of Rock Mechanics and Engineering, 2003, 22(9): 1437-1440 (in Chinese)5 刘远明, 夏才初, 李宏哲. 节理研究进展及在非贯通节理岩体研究的应用J. 地下空间, 2007, 3: 682-693. (LIU Mingyuan, XIA Caichu, LI Hongzhe. Development in joints research and its application to rock mass containing discontinuous joints J. Chinese Journal of Underground Space and Engineering, 2007, 3(4):682-693. (in Chinese).6 LI Y P CHEN L Z WANG Y H. Experimental research on pre-cracked marble under compression J. International Journal of Solids and Structures 2005, 42(9-10) 2 5052 516.7 JIANG L. A review of techniques advances and outstanding issues in numerical modelling for rock mechanics and rock engineering J. International Journal of Rock Mechanics and Mining Sciences & Geomechanics, 2003, 40(): 283-353.8 BELYTSCHKO T, BLACK T. Elastic crack growth in finite elements with minimal remeshing J. International Journal for Numerical Methods in Engineering, 1999, 45(), 601-620.9 NGUYEN V, RABCZUK T, BORDAS S et al. Meshless methods: A review and computer implementation aspects. Mathematics and Computers in Simulation, 2008, 79(3), 763-813.10 庄晓莹，蔡永昌，朱合华. 锚固边坡稳定性分析有限元模型及锚固效应的探讨J. 岩土工程学报, 2008, 30(7): 1099-1004. (ZHUANG Xiaoying, CAI Yongchang, ZHU Hehua FEM modeling method for anchored slope stability analysis and discussions on anchoring effect. Chinese Journal of Geotechnical Engineering, 2008, 30(7): 1099-1004.11 雷远见, 王水林. 基于离散元的强度折减法分析岩体边坡稳定性J. 岩土力学 2006, 27(10): 1693-1698. (LEI Yuanjian WANG Shuilin. Slope stability analysis of jointed rock slope by strength reduction method based on UDEC J. Chinese Journal of Rock and Soil Mechanics 2006 27(10) 1693-169. (in Chinese)12 焦玉勇，张秀丽，刘泉声，陈卫忠.用非连续变形分析方法模拟岩石裂纹扩展J. 岩石力学与工程学报, 2007, 26(4):682-691. JIAO Yuyong, ZHANG Xiuli, LIU Quansheng, CHEN Weizhong. Simulation of rock crack propagation using discontinuous deformation analysis method J. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(4):682-691. (in Chinese)13 甯尤军，杨军，陈鹏万. 节理岩体爆破的DDA方法模拟J. 岩土力学, 2010, 31(7):2259-2263. (NING You-jun, YANG Jun, CHEN Peng-wan. Numerical simulation of rock blasting in jointed rock mass by DDA method J. Chinese Journal of Rock and Soil Mechanics, 2010,31(7):2259-2263. (in Chinese)14 庄晓莹, 朱合华, 蔡永昌, 周丁恒. 基于无网格-图论方法的岩体边坡稳定性分析研究J. 力学季刊, 2008，29(4): 537-543. (ZHUANG Xiaoying, ZHU Hehua, CAI Yongchang, ZHOU Dingheng. Rock slope stability analysis based on meshless method and shortest path algorithm J. Chinese Quarterly of Mechanics, 2008, 29(4):537-543. (in Chinese) 15 陈炜, 王均星. 节理岩质边坡的块体元塑性极限分析下限法. 岩土工程学报, 2008, 30(2): 272-277. (CHEN Wei, WANG Junxing. Plastic lower bound limit analysis using block element method for jointed rock slope J. Chinese Journal of Geotechnical Engineering, 2008, 30(2): 27