刘晓宇-CDEM的流固耦合及块体破裂模型

1、CDEM的流固耦合及块体破裂模型地质体的非连续性 岩土体作为一种地质体，在漫长的地质过程中会形成不同尺度的结构面（断层、软弱夹层、节理裂隙）。结构面的存在控制着岩体的变形、强度。地质体的流固耦合特性 由于不同尺度、方向的结构面所构成的网络系统成为地下水储存、运移的通道，这使得岩土体又赋存于一定的渗流场环境当中，渗流场与应力场之间的耦合作用(相互影响作用)是岩体力学的重要特性之一。地质体的渐进破坏 岩土体在各种外部因素作用下，会在不同部位产生局部应力集中，并在某些变形加剧区域产生破裂面。随着外部荷载、材料强度的变化以及破裂面释放的应力在岩土体内部的转移和调整，新的破裂面继续萌生，已有的破裂面不断

2、扩展、汇合，直至整个岩土体整体破坏。一、不连续地质体建模技术一、不连续地质体建模技术二、裂隙岩体渗流应力耦合模型二、裂隙岩体渗流应力耦合模型三、裂纹扩展模拟技术三、裂纹扩展模拟技术主要内容主要内容结构面（不连续面）网格划分块体切割Incomplete fracture geometry information in field surveys 不连续地质体建模技术预设大尺度结构面软弱结构面软弱结构面手工预设贯穿岩土体的大尺度结构面Incomplete fracture geometry information in field surveys 不连续地质体建模技术规则裂隙分布手工预设结构面In

3、complete fracture geometry information in field surveys 不连续地质体建模技术不规则裂隙分布随机结构面Incomplete fracture geometry information in field surveys 不连续地质体建模技术裂隙统计分布 根据勘测所得数据，运用统计学方法建立包括裂隙面的密度、迹长、产状和开度在内的各组裂隙几何参数的概率模型 裂隙系统取向迹长Incomplete fracture geometry information in field surveys 不连续地质体建模技术裂隙统计分布 假定裂隙面在空间上为圆盘

4、状，利用每组优势裂隙面的数目，然后利用Monte-Carlo 模拟原理生成该组裂隙中各裂隙面的几何参数。其中，利用均匀分布随机数生成器，随机生成服从均匀分布的圆盘中心点坐标；利用迹长（或半迹长）的概率分布类型和分布参数，利用随机生成圆盘的直径；利用裂隙产状的分布特征和规律随机生成圆盘的倾向和倾角；利用裂隙开度的统计规律随机生成圆盘的开度。 RzZyYxXzZyYxXcccccc222)()()(0cos)(sinsin)(sincos)(Incomplete fracture geometry information in field surveys 不连续地质体建模技术裂隙统计分布 利用随机

5、生成的各个裂隙圆盘的中心点坐标、倾向、倾角和圆盘直径或半径，生成一系列的平面方程和柱面方程，逐个切割完整的岩体的模型。 第一条裂隙切割第二条裂隙切割第三条裂隙切割Incomplete fracture geometry information in field surveys 不连续地质体建模技术裂隙统计分布 利用随机生成的各个裂隙圆盘的中心点坐标、倾向、倾角和圆盘直径或半径，生成一系列的平面方程和柱面方程，逐个切割完整的岩体的模型。 平面切割圆柱切割Incomplete fracture geometry information in field surveys 不连续地质体建模技术裂隙统计

6、分布 利用随机生成的各个裂隙圆盘的中心点坐标、倾向、倾角和圆盘直径或半径，生成一系列的平面方程和柱面方程，逐个切割完整的岩体的模型。Successive space subdivision approach 不连续地质体建模技术组号条数倾向角开度迹长均值方差均值方差均值110136.218.60.0050.0013.4221056.611.50.0050.0013.440)()()(ccczZCyYBxXA22222)()()()()()(cccccczZCyYBxXARzZyYxXcossinsinsincosCBA，Successive space subdivision approach

7、 不连续地质体建模技术 根据块体单元的切割面自动生成三维裂隙网络系统，并进行三角形单元的网格划分。最后将利用Monte-Carlo方法生成的具有统计分布的开度赋予各个裂隙单元，产生服从开度统计分布规律的裂隙面渗透系数。 一、不连续地质体建模技术一、不连续地质体建模技术二、裂隙岩体渗流应力耦合模型二、裂隙岩体渗流应力耦合模型三、裂纹扩展模拟技术三、裂纹扩展模拟技术主要内容主要内容+裂隙裂隙岩体岩体渗流渗流应力应力耦合耦合作用作用 模拟岩体应力场的模拟岩体应力场的不规则可变形块体离散不规则可变形块体离散元模型（元模型（CDEM）模拟裂隙渗流场的模拟裂隙渗流场的裂隙网络渗流模型裂隙网络渗流模型裂隙岩

8、体渗流应力耦合模型三维不规则、可变形块体 将块体看作内部没有任何不连续面的连续体。基于连续介质力学思想，块体的力学行为可采用三维动力学基本微分方程来描述 btbtbbbbSfiiiiiiijijiiiidStudStudVbudVdVuudVuu 为求解块体域内的动力学问题，采用四面体网格划分，将块体离散为四面体单元 Face CenterVertexBlock Center )()()()(tFtuKtuCtuMexts 接触单元连接单元 当相邻两块体之间存在结构面时，邻居节点为接触单元，遵循非线性本构关系 当相邻两块体之间完全连续时，邻居节点为连接单元，之间采用位移约束 IbIbuu或连接

9、弹簧，遵循线弹性本构关系 裂隙网络渗流方程 非稳定渗流的基本微分方程 fqyyxxfhfyxonqnyhKnxhKonhhinthSyhKxhK 2222裂隙渗透系数123gtKKyxxyztqq对于裂隙流而言，认为渗流只在沿平行于裂隙面的方向发生 达西定律和质量守恒方程thSvhvK裂隙网络渗流方程 VertexSegment Center)3 , 2 , 1()(01nedsqNdxdyNthSdxdyyNyhkxNxhkMelMeeleeeleyelex 采用三角形单元进行有限元数值计算 裂隙渗流模型计算流程)()()()(tQtKutuCtuM 水头：0.001s 水头：1s 水头：5

10、s 水头：10s 水头：20s 水头：40s 水头：80s 水头：160s 水头：240s 节理裂隙岩体渗流应力耦合作用的数值模拟节理裂隙岩体渗流应力耦合作用的数值模拟 试件试件 长长宽宽高高: 40m40m40m渗流单元渗流单元 长长宽宽: 2.5m2.5m初始裂隙隙宽：初始裂隙隙宽：10mm重力重力: 9.8m/s2; 水密度水密度: 1.0e3; 水流粘滞系数水流粘滞系数: 1.0e-6边界线边界线(x=0, z=20): 流量流量q1=0.8m2/s边界线边界线(x=40, z=20): 水头水头H1=50m三组相三组相互正交互正交的裂隙的裂隙组成的组成的规则裂规则裂隙网络隙网络 渗渗

11、流流区区域域内内的的水水头头分分布布 正视图正视图 侧视图侧视图 裂隙岩体渗流应力耦合模型H: 水头水头q: 流量流量Kx, Ky, Kz: 渗透系数渗透系数ij: 应力应力 ij: 应变应变渗流场渗流场模型：模型： ),( ),(),( ),( ),(),(),( 0)()()(22111zyxzyxqzyxqzyxzyxHzyxHzyxzHKzyHKyxHKxzyx应力场应力场模型：模型： ),( ),( )(),( 2 ),( )(21 ),( 022,2,uiiijijijijVijijjiijijjiSzyxuuSzyxHtnzyxGzyxuuzyxf各条裂隙的隙宽各条裂隙的隙宽空间

12、位置空间位置 切向拖曳压力法向渗透压力: 2/: )(0wfwtJbtpzHp裂隙岩体渗流应力耦合模型数值算例 水力耦合Young modulusPoisson ratiodensity8.751070.21.81030 xu0 xu0zyxuuu不考虑水力耦合考虑水力耦合水平位移杨氏模量密度内聚力摩擦系数泊松比30210010500.3粘性系数重力加速度缝宽b密度贮水系数1.01069.81.010-31.01035.010-3算例：裂隙渗流作用下坡体失稳过程模拟The study domain: 20m1m20m水平位移：20000步 水平位移：30000步 水平位移：50000步 水平位

13、移：80000步 水平位移：90000步 一、不连续地质体建模技术一、不连续地质体建模技术二、裂隙岩体渗流应力耦合模型二、裂隙岩体渗流应力耦合模型三、裂纹扩展模拟技术三、裂纹扩展模拟技术 单元界面破裂单元界面破裂 块体内部破裂块体内部破裂主要内容主要内容Incomplete fracture geometry information in field surveys 裂纹扩展模拟技术界面破裂Ks： ：剪切刚度剪切刚度Kn： ：法向刚度法向刚度 当相邻两块体之间完全连续时，邻居节点采用连接型弹簧，遵循线弹性本构关系 当连接型弹簧满足张拉-压剪复合破坏准则时，连接型弹簧将发生张拉或压剪模式的弹脆性

14、破坏，转变为接触型弹簧其中，张拉模式的弹脆性破坏采用最大张拉力准则模拟 ctnAF压剪模式的弹脆性破坏采用莫尔-库仑准则模拟 cnsAcFF0tan 当相邻两块体之间存在结构面时，邻居节点采用接触型弹簧，遵循非线性本构关系 块体块体3块体块体1块体块体2Incomplete fracture geometry information in field surveys 界面破裂算例顺层岩体反倾岩体含软弱带岩体软弱土层一、不连续地质体建模技术一、不连续地质体建模技术二、裂隙岩体渗流应力耦合模型二、裂隙岩体渗流应力耦合模型三、裂纹扩展模拟技术三、裂纹扩展模拟技术 单元界面破裂单元界面破裂 块体内部破

15、裂块体内部破裂 弹脆性材料（岩体）弹脆性材料（岩体） 弹塑性材料（土体）弹塑性材料（土体）主要内容主要内容Incomplete fracture geometry information in field surveys 裂纹扩展模拟技术块体内部破裂自适应网格模拟技术扩展有限元技术块体的破裂模型 C 3 拉破坏面3 3 拉伸破坏通过相对位移选择滑动面：3331111t10tan0nc 剪切破坏常规单元切割切割点与角点重合裂纹扩展切割边界单元切割块体一条边与切割面关系不相交被剖分为两部分该边至少有一个角点在切割面上块体一个面与切割面关系不相交该边两角点P1、P2将继承父节点与其它块体节点的混合节

16、点关系该面面中心点将继承父面中心点与其它块体节点的混合节点关系被剖分为两部分边被完全切割边被不完全切割边界边边被完全切割边被不完全切割边界边面被完全切割面被不完全切割边界面Key point Block cutting technique 0)( 0)( 0)(inoninxxx对于给定的切割平面方程，给出切割面、切割面上方+和切割面下方-的level set函数 0)(min)(min :0)(min)(min :iIiiIiiIiiIielelelel不切割块体切割块体)(122111xxxx块体边与切割面的交点坐标： 利用level-set函数遍历块体所有面，检查各面与切割面之间的关系

17、常规单元切割插值关系混合节点混合节点接触节点1。寻找混合节点；2。寻找接触节点；3。寻找插值节点；边界单元切割接触节点混合节点混合节点接触节点接触节点自由节点1。寻找混合节点；2。寻找接触节点；3。寻找插值节点；裂纹扩展切割混合节点接触节点或混合节点1。寻找混合节点；2。寻找接触节点；3。寻找插值节点；混合节点块体一条边被切割面剖分为两部分1。建立新节点S1、S2与角点P1、P2之间的插值关系；2。将位于切割点的所有新节点定义为克隆关系；Y该边两角点P1、P2将继承父节点与其它块体节点的克隆关系该边是否为边界边Y将新产生的两个节点S1、S2定义为接触关系；N该边是否为完全被切割NY1。消除原有

18、新节点的插值关系；2。将每一对新节点（S1,S2）都定义为接触关系；3。将位于接触面上方的所有新节点定义为克隆关系；4。将位于接触面下方的所有新节点定义为克隆关系；块体一个面被切割面剖分为两部分将面心节点看作自由节点Y建立面心节点与相邻面内节点之间的插值关系该面是否为边界面YN该边是否为完全被切割NY1。消除原有面心节点的插值关系；2。将两个新面心节点定义为克隆关系；切割点与角点重合接触节点混合节点内部块体边界块体1。寻找混合节点；2。寻找接触节点；3。寻找插值节点；混合节点切割点与角点重合1。寻找混合节点；2。寻找接触节点；3。寻找插值节点；该边至少有一个角点在切割面上Y该边是否为边界边Y该

19、边两角点P1、P2将继承父节点与其它块体节点的克隆关系N1。将子块体新产生的、位于接触面上的两个节点S1、S2定义为接触关系；2。将节点S1、S2继承的父节点与其它块体节点的克隆关系进行重新筛选：位于切割面之上的子块体将继承的父节点中与位于切割面之上的其它块体节点的克隆关系；位于切割面之下的子块体将继承的父节点中与位于切割面之下的其它块体节点的克隆关系。该点是否为完全被切割Y1。将子块体新产生的、位于接触面上的两个节点S1、S2定义为接触关系；2。将节点S1、S2继承的父节点与其它块体节点的克隆关系进行重新筛选：位于切割面之上的子块体将继承的父节点中与位于切割面之上的其它块体节点的克隆关系；位

20、于切割面之下的子块体将继承的父节点中与位于切割面之下的其它块体节点的克隆关系。Incomplete fracture geometry information in field surveys 裂纹扩展模拟技术单裂纹扩展(1). 每个迭代时步结束后，计算裂纹前端单元的应力场；(2). 当该单元应力场满足破坏准则时，裂纹在该单元内扩展至其边界；(3). 裂纹扩展方向沿破坏准则给定方向；Incomplete fracture geometry information in field surveys 裂纹扩展模拟技术单裂纹扩展sP/2P/2Pbas=12; b=4; Incomplete frac

21、ture geometry information in field surveys 裂纹扩展模技术多裂纹扩展Incomplete fracture geometry information in field surveys 裂纹扩展模技术多裂纹扩展底面三个方向约束P=3.0107(Pa)一、不连续地质体建模技术一、不连续地质体建模技术二、裂隙岩体渗流应力耦合模型二、裂隙岩体渗流应力耦合模型三、裂纹扩展模拟技术三、裂纹扩展模拟技术 单元界面破裂单元界面破裂 块体内部破裂块体内部破裂 弹脆性材料（岩体）弹脆性材料（岩体） 弹塑性材料（土体）弹塑性材料（土体）主要内容主要内容弹塑性材料本构模型 弹

22、塑性材料本构方程klijklijD )( )( 塑性阶段弹性阶段kleijklijemnklmnsteijsteijkleijklijklFDFHDFFDDDDE00YF 弹性阶段：peeE 塑性阶段：)(pE屈服函数：0FF一致性条件FEFFEFEFEFFEEFp)sgn(流动法则弹塑性材料本构模型 弹塑性材料本构方程klijklijD )( )( 塑性阶段弹性阶段kleijklijemnklmnsteijsteijkleijklijklFDFHDFFDDDD弹塑性材料本构模型 )sin3(3cos6)sin3(3sin2 12ckakIJFDrucker-Prager 模型：tannFMo

23、hr-Coulumb 模型：强度折减方法计算边坡安全系数40000步80000步120000步160000步200000步240000步280000步320000步 问题由来 其一是缺乏统一的边坡极限状态的判定标准，亦即边坡失稳判据 数值解非收敛判据：认为当边坡处于极限平衡状态时，非线性有限元方程组的迭代过程将不收敛。由于收敛的标准通常是由用户指定的失衡力的大小来控制的，且不收敛可由很多其他因素而引起，并不一定意味着结构已经崩溃。 塑性区贯通判据：认为某一幅值广义剪应变或广义塑性应变等值线由坡底贯通到坡项，则意味边坡已经失稳破坏。研究者必须有相当的经验才能选取合理的临界应变幅值。 位移突变判据：认为当边坡达到极限状态时，某些点的位移会趋于流动状态。由于在某些复杂情况下选取不同的特征点，可能会给出不同的安全系数。并且，选取哪一点作为极限状态也尚未取得共识。 由于失稳判据或难以进行量化，或在分析过程中受人为因素的影响，或存在诸多不确定性，其将直接影响到安全系数的计算结果。 问题由来