《岩石力学的数值模拟(讲义).docx》由会员分享，可在线阅读，更多相关《岩石力学的数值模拟(讲义).docx（14页珍藏版）》请在课桌文档上搜索。
1、第10章岩石力学的数(ft模拟随着计算机软硬件技术的快速发展，使岩石力学有了长足的进步，特殊在岩石力学的数值计算和模拟方面发展尤为快速,使得很Z着看力学解析方法难于解决的问题得以重新相识.正如钱学森在给中国力学学会“力学一一iffi接21世纪新的挑战”的一封信中M力%发展均势总结的那样“今日力学是一门用计用机计算去I可答一切宏观的实际科学技术问8S，计豫方法特别象要,o岩石力学和其他力学学科一样，须要数优计豫方法并推动岩石力学的发展。岩石介质不同干金屈材料,在数值计算方面具有其独特的特点2O5:(I)岩石介质是赋存于地壳中的各向异性自然介质.(2)岩石介质被众多的节理、裂缝等弱面所切割而呈现高
2、度的非均质性，而其物理、化学及力学性侦具有时机性特点。(3)岩石介质赋存时以受J卡为主，而且抗压强度远大于抗拉强度.(4)岩石力学及JC程同超在时空分布上较广，从本质上讲都是三年同SS.(5)岩石工程一般无法进行原型试验，而试抬空测得的数据不能干脆应用于工程设计和计算.(6)岩石力学及工程具有数据有限问胞.数值计。方法羟过几十年的发展，目前已形成很多种岩石力学计W方法，主要有有限元法、边界元法、彳i限差分法、离散元法、流形元法、拉格朗口元法、不连续变形法及无单元法等。它们各有优缺点，有限元的理论基础和应用比较成熟，在金属材料和构件的计算中应用特别胜利.但它是以连续介质为基础,好像及岩体的非连续
3、性有确定差距,源形元等数值方法虽然考虑了岩体中节理效应,但其理论基础还不完全成熟.信任在不久的将来，确定会出现完全适合于*体材料和工程的数侑计算方法2O62O81。10.1 岩石力学的有限元分析2O92I3)行限元法(IicntmctiKK1.F1.iM)是岩石力学数(IIi计算方法中最为广泛应用的一种.自20世纪50年头发展至今，布取元己胜利地求解了很多困碓的岩石力学及工程同腮.被广阔界石力学探讨及工程技术人员嗡为斛决岩石工程何阳的有效工具.彳i限元法是依据变分原理求解数学物理方程的一种数值方法,有限无法把连续体离散成有限个单元，旬个单元的场函数只包含有限个节点参Ift的简洁场函数.这些有限
4、个总元的场函数集合构成整个结构连续体场函数.依据能成方程和加权函数方程可建立有限个求解参数的方程组,求解这些离敌方程粗，就是有限元法的精饰所在，虽然求解时把连续函数转化为求解有限个离散点处的函数值，但只要单元划分得充分小时，足可以酒意计算要求。有限元法求耨问题时一般遵循以下步骤：(I)有限元计分模型的建立，包括模型总元的划分、确定边界条件.(2)对单元体进行力学分析，包括求解节点位移、单元应变和单元应力。(3)对计算模型进行分析。(4)进行计算分析.10.1.1 线弹性有K1.元法的基本方程线弹性有限元是弹塑性有限元、扭伤有限元、流变有限元等非线性有限元的基础，线弹性有跟元假定岩石介质连鼻均施
5、、小变形和完全弹性.有限元法求解弹性力学同胞时通常以位移作为基本未知用单元位移是以单元节点位移为基本未知疑.选择合埋的位移插色函数，将单元位移表达为节点坐标的连续函效，捅依函数也可称为形函数，不同形态的单元具有不I可的形函数.图10-1为三种最常见单元形式.即三角形、四边形及四面体单元。它们的形南敢分别为:图IO-I有限元的三种茶本单元形式（八）三角形单元(b)四边形单元(C)四面体明元三角形的形函数式中，S为三角形面积；ai=xiym；bi=yi-yt；ci=xfxj四边形的形函数M=1+49(1.+7)0=1.2.3.4)4式中，位移球为“=EN,%:v=NiviiX=EN用；y=Niy1
6、.(i=i.j.k.m)四面体的形函数Ni=-(,+b,x+ciy+diz)OV式中，V为四面体的体枳.单元在出角坐标轴中位移重St分别为V.v,因此第元的位移矩阵为we=(M,v,ie=()(10.1)式中，,.为单元位移矩阵：NJ为形函数电际：13/为单元节点位移列阵.依据几何方程，对位移矩阵求体导数，可以知到应变矩阵=网但(10.2)式中.小I为连续单元节点位移和单元应变的电阵.也称为应变足际.对于三角形的元,BJ为常数矩阵，元素值取决于单元节点坐标差.依据本构方程，可以得到单元节点位移及单元应力矩阵之间的关系(j=DtfJ=ID1.IWJ103)式中，D1.为弹性矩阵.应用虚功惊理和最
7、小势能原理可以推导出IR元附度矩阵的在达式AfJ=7D1.kv(10.4)各单元的体积力和面力依据静力的等效原则移置到各型元的节点上.其等效节点力为(10.5)式中，&为作用于单元体积力IP)的等效节点荷我；Qt为作用于单元面枳力Q的等效节点荷载.设环绕某节点i共有4个单元，则i节点上的外加苛或怵,)为dF,=XjBrD1.,ddV(10.18)由于）随位移改变而改变，因此计尊时必需进行迭代求解。初应力法求解依据卜述计算步骤实现：（I）把全部荷我划分成若干个增量，在每一级增限段内.依据增玳弹塑性平衡方程进行求解,（2）计算各单元的应力增量及当前的应用dii=Bdiidii=Ddii10.19）
8、,b=,b.+M1.式中.下标i表示第i级荷载培崎：表示第j次迭代.（3）依据岩石的屈服准则，由各单元应力推断单元是否屈服.对于型性雎元，计算应力修正项井修正应力d%1.=DM1.1.dGn叫=矶-m*(10.20)（4）塑性IR元通过修正项m,j计算等效节点力，全部塑性总元的等效节点力会加构成总的修正荷毂矢量mj=（dp）jdv10.21）（5）在修正荷载作用下进行下次迭代运算，此时基本方程为Kd,j=dF,j（10.22）由系进行（2）（5）步计算直至全部的理性单元达到收敛精度要求。再进行下一步的荷就增量计算。（6）曳新施加下一鞭荷我增1.tK+,王复计算（1）（5）步，直至计算完毕，通过
9、累加各级荷载作用的计算结果,求得总位移“和总应力b一般初应力法的收敛速度比较缓慢.因此遹常采纳常榭度和变刚度法相结合的方法加速收敛.初应变法依据弹型性增及理论，总的应变呆可表示为i=drdp（10.23）式中,de,.为弹性应变增盘：de。为型性应变增S1.类似初应力法可以把电性应变if1.依旬0看做初应变,因为在计算过程中而0的计算格式和界性系统中的初应变“与）一样.图1.（K6初应变法从图106中可以看出，荷我增玳W所对应的位移为小。.按线弹性计算为由以,两界的位移增疗之差称为初始位移，它是由材料电性引起的附加位移。及模型系统初始位移对应的单元位移为“%,那么单元的初应变为(J)=IBI(
10、JJ0)(10.24)10.1.4岩石节理单元的行限元方法岩体中常存在大量节理，而节埋的变形性质和岩体力学性质有特别明显差别.从某种程度上讲，节理的存在确定希*体的力学性旗。因此分析节理的变形性质，对于探讨岩体工程的变形态况至关垂要，在进行有限元分析时，这种节理一般来熟特地的节理电元来处理，节理单元首先IhGtXK1.man提出“小？川,GOOdman认为节理不是一个实体.它只在若干点处接触.因此节埋单元也应满意这一特点，图167为无厚度节理单元，节点I及%2及3具有相同的坐标，沿节理面的应力分别为法向应力b“和剪切应力口.把节理单元两(M对应的位移定义为应变,相对的法向位移】，称为法向应变，
11、相对的切向位移称为切向应变J，那么节理单元的本恂关系为=0(10.25)图10-7无厚度节理单元式中，/)为单元的刚度矩阵,对于节埋单元长度上任何一点S处的应变可定义为界面两侧相应位移之差J=(I-I)(V4-%)+-(V3-v2)emmethod.BEM)是和彳i限元法同步发展的一种数值计算方法。及有限元相比有以卜特点，边界元法把一个均质区域看做一个大单元，只把它的边界离散化,区域内不划分单元.场变量到处连续：对于无限区域.场变用自动满意无穷边界条件及自然衣面状态.对于半无限区域.场变珏也自动满意无穷远边界条件及自然灰而状态.有限元法是全区域离散化，而边界元是把基本方程转化为边界积分方程，只
12、对边界盛敢化建立相应的方程纲进行求斜，这样边界元使三维向胞降为二维问题求解，使二维问题转化为一维问遨求解.当物体的表面积和体积之比较小时,边界元的划分胞元数要比有限元少数倍和十几倍，这样也使待解的方程数目、处埋和存储的数据依降低同样的倍数.大大节约了机时和肾遨商用。当仅需求斛物体内部几个点的应力时，有限元仍不得不划分整个区域，才能确定这几个点的应力值,而边界元当知道边界的应力解时，就可以依据须要去求物体内部个别点的解,化应力梯度较高处,有限元法的剖分密度常常受到限制.而边界元可以便利地确定应力梯度的分布.但K1.蕾计算机硬件的飞速发展，边界元的优势渐渐显得不很明显得边界元法矩阵方程中系数阵的元
13、素结构比有限元法刚度阵中的元素困就。有照元刚度阵属带状稀疏阵，而边界无法的系数阵为湎阵,因此时于面枳和体根之比较大的薄壁结构而言.边界元的优越性就不明1.边界元比较适合求解无限区域和半无限区域问应，如深埋普道是一个典型的例子，但边界元在计算非均质介质问题、非雄性何期以及模拟工程开挖过程等方面不如有限元便利有效.有限元及边界元划分单元的比较如图10-8.rao8有限元法及边界元法比较(a力学模型和边界条件b)有限元单元划分(C边界元单元划分求解边界方程组土要仔Massnnct和Rizzc分别提出的干脆和间接两种数值解法.干脆法是干脆建立关于边界的枳分方程，通过离散化求解边界未知参数，并进而求解计
14、能区域内场函数瓶.间接法是设立个在求解区域内含有若干系数的满意基本方程加的解,代入边界条件求出系数.进而求得边界及区域内各点的场函数值.10.2.1边界元分析原理a111.在无限的弹性体内作用一单位力引起四周的应力和位移的解析解是边界元求解弹性体问题的基础，如图10-9在平面无限体内i点作用一X方向的单位力，其基本的弹性解在1848年田Ke1.vin解出%=-7Ti1.T固1吗+畀刃一停为/3D411(1.-v)rcnxfxkJx1.xkJJ(10.32)uU=1.:(3_4V)如+白舁8v(1-v)rex,vxkJ式中.为i点/方向单位力引起J点8方向的应力和位移：A为Kroncckcr符号
15、.当/=*时，传=1，当火时，a=0；r为。/两点之间的距离(矢径)；，“，,为pr，点作用面外法向对于k和/的方向余弦;F为矢径r对外法向的方向倒数：;ocn当不考虑体力时，依据功的互等定理和以上的KCMn基本解.可以建立弹性体边界上人，两点之间的应力和位移关系(如图10-10)c,tt,+u=t(10.33)图10-9开尔文基本解条件图10/0边界积分方程的建立式中.J的大小取扶于边界状况,当边界为光滑曲线时.C等于0.5：当边界曲线不光滑时.c的值依据刚体位移原则求解，当边界作用分布力时，/点绕行一周，按小加原理积分得c*ut+uiV=jf1.(10.34)式(10.34)就是边界元中的
16、边界积分方程。10.2.2边界单元及边界的离散c1.cmcntncth,加上块体B的IR力，它们对块体B的汆心产生合力F和合力矩M,即七七%z1.=1.z-iz1.-1.-FFv-M(10.64)XC)F(为-yj图1047块体的集合及作用于个别块体上的力式中，勺，为块体受力作用点坐标；“为块体无心坐标，假如上式中的合力和合力矩不等于零，那么块体B将依据牛顿其次定律的规律运动，令块体B的旗W;为,“，转动惯最为1，块体绕其重心转动角度为.那么块体的运动方程为(10.65)式中，“为块体位移：g为重力加速度.对上式可果纳差分方法进行求解.对X方向的运动方程采纳向前分格式进行数frt积分，这样可以
17、得到岩块质心沿方向的速度和位移11t(0+)=zix0)+M1r10.66)ur(f0+)=Mt0)+zi,式中，川为初始时间：Ar为计算时步。对于y方向的运动和行动，都有类似的算式.虽然离散单元可被视为不变形的刚体，但在单元接触的边界仍存在变形，如图10-18设块体之间的法向力H1.两块体之间的“弟合”位移为如两者之间成正比关系F1.I=KM0-67)式中，Kn为接触面的法向刚度系数.自己的固定位置.处于平衡状态.“外界外力或边界位移约束条件发生改变时,某些块体在史力和外力作用下,依抵运动方程产生确定的加速侬并产生位移，因而使块体在空间的位徨发生改变。产生位移后块体及相邻块体在空间位置上发生
18、IRf1.S依据块体力和位移关系，使更多的块体由于块体间的重登而殳力，于是又产生新的运动和位移，依次迭代下去，遍及整个块体集合,计算模拟各个块体位移和转动的整个过程.直到以终收敛达到平询状态或所10.5岩石力学的流形元分析实际的工程岩体常被一些节理等结构面所切割,形成一种处于连续及非连续之间的拟连续介质，囚此连续介质力学数值方法，如有限元、边界元、有限差分方法等，处理岩体节理时采纳节理单元方式，而非连续介质力学数值方法，离微元、I1.连续变形分析(disconinousdefo11natitmana1.ysis.DDA)处理事先划分好的块体模型比较相宜。而新近兴起的流形元对于处理连埃及非连续介
19、质朋合同跑比较相宜,对于岩石材料尤为适合.流形元方法(manifo1.de1.ementmethod.MEM)是我国8名留美学者石根华、林德琼博士于20世纪90年头初提出并开发的种全新数值计算方法.波形元法以拓扑流形学为基础.应用有限网盖技术,通过在计算区域内各物理覆盖上隹.立通用的濯盖函数和以加权求和形成总体位移函数.从而把连续和非连续变形力学问即统一起来考虑.它以数值流形为核心.在DDA基础上，联合有限元和解析法的连续分析方法，从而形成在一切空间(包括彳!限元、DDA和解析法统一的计优形式I?”?切，有限殂薇技术是在数学流形分析中常常采纳的一种方法。石根华采纳这一方法,统一斛决了连续和非连
20、续变形的力学问即.有限覆盅技术构造的覆i函数具有以下两个基本性质：同部非零性，装瓶函数只在局部区域内不为零，即在局部区域内有解：在求解区域内，焉能用数之和为1,即局部区域内全部海盐团数组成的总体位移函数，它及数学流形的区分在于数学流形在整个求湃区域内总体位移南散到处连续并光滑可导，它UJ完全定义而及网盖无关.而数(fi/形的总体位移函数是在榄前基础上定义的,且可分片微分,在覆彘的接触向上是完全非连续的.总之，数值流形的基本特征是覆忐函数在同部区域内连续，各分片区域之间班族函数不连续，通过采纳连续和非连续粗靛函数的方法，把连续和非连续变形的计算统一到数值流形中.10.5.1数低流形方法的行限覆盅
21、流形元有两套有限覆盅技术，数学覆盅和物理覆虚.其实覆盖的概念和书限元的网格基本一样，也就是说流形元有两层独立的网格，数学网格和物理冏格，数学网格侬定义数值解的精度,它由用户依据须要确定；而物理网格由材料的形态、裂隙、边界以及分区截面等材料的物理和集合性质确定,用户不能随意改动和选择.因此物理置盅”,就是计算边界、材料分区和裂隙等所组成的刈格：数学置维V1是数学网格节点所构成的单元,数学覆荔可以是任何形态的，而物理班战完全由材料的物理和几何性侦确定。流形元方法就是把物理和数学两层网前重登在一起，并用物理覆盖把数学陵盖重新划分,把数学用盖重登入物丹法盖.形成新的连城和非连续耦合的有限覆盅系统,在这
22、个系统中,物理覆盖代替单元的节点.用益的并集交战代替单元的边界,覆养曲光的交集代替垠元.图10-21有限圈盖和流形单元的描述219（八）柄个块体的有限强盅(b)块体内有条裂隙的有限覆盅图10-21分别采纳四个行限单元网格翟燕个四边形的材料区域，四个有限单元的节点辑(I.2.3).(2.4.5),(2,5.3).3,5.6)。图10-21（八）表示在有4个原始单元的数学隈盅细线表示上描i有两个块体的物理覆盖(机线表示)的行限没盅系统.图10-21为同数学桁靛上有条不连通裂隙的块体状况.图中大数字表示节点相鼓码.小数字下标表示被不连续的物理覆靛划分的分区码，图中I”1“2”2?等分别表示节点1，2
23、的物理覆盖。10.5.2流形无网格的覆盖函数和权函数位移函数及材料边界无关,假如材料只占有单元的一部分，位移函数仍旧是相同的.对三角形单元，相应三个节点有三个用蔽包含这个单元，因此斑个单元是它三个节戊的三个题盖公共区域，行人中，!要计算权函数.5.F)表示为节点元1,2.3的坐标.则相应节点e(I)e(2),的坐标和位移可表示如下：坐标位移e(1.):(1.)ve)e(2)：(X心M2)-(%2%2)”3):(Xwyd3)T(%3),03)流形元中单元随意一点/(X，,V)的权函数为式中1.2A1/12f22fi1.XqX一X心匕匕2)一匕(j-Xm工匕-4)e(1.)*f(2)XgKtI)M
24、-M22jIf22.X(10.72)/32hyX3)X2)Xzd-XZD(10.73)Xgfr1.UI/11W2(工#=%MU)；hi它取节点位移的百分数,(10.74)式中,、!(x,y)为节点j没盖的权函数，其意义为加权的平均，对三节点单元，三个节点的三个权函数之和为1,权函数相当于有限元中的形函数,物理用鼓。上的攫龙函数为3,(x.y),匕(,y).可以是找性的，也可以是非线性的.每一个节点用装生成的位移为覆忐函数及权函数的乘积，总元各点的最终位移则为整个物理粮法系统上的总体函数，如覆萩函数用标准：阶级数，JE么总体函数为心3=二(10,75)vr1.)1y)Wrt0(-.V)/-27-
25、I建立了每个流形电元上的位移函数以后,就可以依据有限IR元方法形成单元的单刚及总刚度矩阵进行求解，10.53数值流形方法的平衡方程式对二维三角形的元力有三个物理攫养或节点力,力,力.每个物理覆盖或节点i有两个未知数(Ui,Vi)(10.76)将全部的势能相加，总势能有以下形式q1K*D2K6+DFi+C(10.77)KnnD11Fn式中，O=(OIK=Kjr总势能包括应变势能、初应力势能、点荷载势能、体背或势能、惯性力势能、接触弹簧的应变势能和摩擦力势能.对鬼蒜，具有以下方程式,(10.78)上式表示作用在泡盖r上的所用荷载或接触力在X和)方向上的平衡方程式，10.5.4单元IS或接触的进入理
26、论在有限用差中有连续介质部分和不连续边界，但不连续边界必箭连接成一个系统.不连续边界的位移必需满意在接触面之间无拉力和无嵌入，以覆靛接触为拓础的运动学理论，用刚性接触弹簧计算不连续变形，它符合以卜两个原则；有限陵盖上全部点的位移小于规定的极限值?：转动角小于3因此计算时时间步选得足筋小。这样各点位移(u(,y),v(,y)转动r(x,y)和变形%，g，%n.能精确地依据海孟未知数IQ1.的线性函数来描述。考虑单元e和随意一点(.%)，)e届么它的应变为=方约(x,y)M210.79)r=1.式中，8为系数：。为覆盖位移：P为置蒜(U,),U2),。“凉的交集：180,180OEOE2a180.
27、180OEOE4为了防止接触的两条边的嵌入.在进入线上设置倒性弹费.洋细设置有以下三种状况：（1）的对用的接触.若两个角都小于180,在角顶点和它首先越过或梢有越过的线之间加一个弹簧，（2）角对角的接触,有一个角大于180.如两条进入城中任何一条被其相对应的另一角的顶所越过,在沿这条边的垂出方向加一个弹簧.如两条进入线都被越过，则在两条进入线的边上各加一个刚性弹费.（3）角及边的接触.在进入我的边上加一个弹黄.由于在有限覆盖之间加入了弗箫，在计舞也盖接触他阵时，应考虑弹贲的存在。如图10-12,从角点P1.到PiP,的进入线的距离d为j1x+,+jd=-=-X,+U2y,2+V2（10.80）
28、1勺+“33+式中.I=（X2+42-4-二）2+（立+丫2-g-二尸法向弹簧的势能为W”=犷=d+G门+率）ra10-22进入线的刚性并筑式中，”用坐标及位移矩阵T(x,.Y)形式在示时,此公式中各项详细如下假如点(.r1.),.y0)在进入我P2Py上，它也是用点P1的假定接触点,则明切弹彼在PzP,方向上及小和几点接触,剪切弹簧的势能足叱=#/I1。”针33十(10.85)式中(10.86)(10.87)G卜10.88)当摩擦角不为零时.法向弹费力就会产生摩擦力.在修一边摩擦力F的势能为VV1-=FDr(10.89)式中向=沁一);:在外一边摩擦力户的势能为Wf=-FDiiir(C10.
29、90)式中G=M(n,)b)一，,1.1.1 MfJ对以上各项势能取微小值,可以得到各项系数，加到总矩阵中去。10.5.5 数值流形的单纯形积分流形方法解积分方程时及有限元不同,有限元处理的是连续介质，划分的博格般采纳简沽的标准形态，可以进行复合形枳分，流形方法一般是解非连续体，它的有限徵族形态很困难，不能随意简化。因此它须要在更合形态上进行单纯形积分，其处理方法为将随意网雉的形态分解为最简洁的中.纯形.在单纯形上进行单纯形积分然后求其矢量和.图10-23为复合形上的单纯形积分.已知多边形A%R凡RRI.同时为=尸6,对于的意一点人，：维单纯形的面积AP俨2,PP：P”PgP、，凡孔PS和八/
30、5P1.的代数和是多边形的面枳4(其中单纯形岛的面枳是负的，其余都是正的)=八+C1.B玛+昂玛B+为鸟鸟+&居A(io.9)图10-23双合形上的单纯形枳分对于被枳函数f(x.y),多边形的标准枳分是单纯形枳分的代数和A鸟孱6=/(x.y)D(x,y)=Xv(J)fiPiPif(x,y)(ix1.y10.92)I-I式中的三角形面积具有正负之分，则Sign(J)中的)为行列式(10.93)对于随意形态的三维史合体可用三维单纯膨积分.同时单纯形积分可扩展到”维.在数伯流形的积分方程中除计算体枳的枳格外，还要计算一次和二次矩，它们在计分时还要进行坐标转换。10.5.6 流形方法的应用岩石工程常遇
31、到不连续结构而和困难边界问起.在正常的数值分析方法中,一般实行特殊的处埋方法或简化.而流形元采纳数学网格和物埋网格两套覆盖系统对计算域内不连续部分和连续介质部分曲亍统一分析，成为全新统,的数例方法，因此它具有特别深远的发展潜力和广泛的应用前景，(1)自由表面和柔性边界,分析时不受边界条件的阻碍.(2)单元形态随意，可以同时诳行连续和非连续分析.(3)能M守恒，并符合库仑定律，(4)可进行动力学和龄力学分析，也可进行小变形和大变形分析.作为应用实便,图10-24为层状结构岩石受集中力破坏开裂破坏状况.图10-25为岩石块体边坡滑动结果，滑动时中间块体分别成两相邻块体，其结果及试验室试验结果样.图10-24层状结构岩石受集中力破坏开裂破坏状况224)（八）图IO25岩石块体边坡滑动结果224
