具有求解岩土破裂演化问题能力的有限棱元法.pdf

第3 7 卷第1 期中国矿业大学学报 V 0 1 ．3 7N o ．1 2 0 0 8 年1 月J o u r n a lo fC h i n aU n i v e r s i t yo fM i n i n g T e c h n o l o g yJ a n ．2 0 0 8 具有求解岩土破裂演化问题能力的有限棱元法 宋宏伟，赵高峰 中国矿业大学建筑工程学院深部岩土力学与地下工程国家重点实验室，江苏徐州 2 2 1 0 0 8 摘要为了解决岩土破裂演化过程数值模拟问题，以有限元法三角形单元为例，从其棱边变形导 出了其单元的刚度矩阵，发现有限元法的单元可以进一步划分到棱边一级，提出了一种新的有限 元算法，即有限棱元算法 F E E M ．以棱划分单元，可以方便地实现力学介质从连续状态到非连 续状态的破裂演化分析．用所提出的有限棱元法开发了程序，并用其分析了地下岩洞的应力和变 形，模拟了圆盘破裂演化过程，其结果与利用其它大型商业软件计算而得的结果相近． 关键词数值分析；算法；有限棱元法；破裂演化；连续一非连续；破碎 中图分类号T D3 5 3 文献标识码A 文章编号1 0 0 01 9 6 4 2 0 0 8 0 1 0 0 1 5 - 0 4 F i n i t eE d g eE l e m e n tM e t h o df o rU n i t e dS i m u l a t i o no f C o n t i n u o u s D i s c o n t i n u o u sP r o b l e m s S O N GH o n g w e i ，Z H A OG a o f e n g S c h o o lo fA r c h i t e c t u r ea n dC i v i lE n g i n e e r i n g ，S t a t eK e yL a b o r a t o r yo fG e o m e c h a n i c sa n d D e e pU n d e r g r o u n dE n g i n e e r i n g ，C h i n aU n i v e r s i t yo fM i n i n g T e c h n o l o g y ，X u z h o u ，J i a n g s u2 2 1 0 0 8 ，C h i n a A b s t r a c t I no r d e rt os o l v et h ep r o b l e m se n c o u n t e r e di nm o d e l i n gt h er u p t u r eo fg e o t e c h n i c a l m a t e r i a l s ，t h eu n i ts t i f f n e s sm a t r i x ，t a k i n gt r i a n g l ee l e m e n to ff i n i t ee l e m e n tm e t h o d F E M a s a ne x a m p l e ，i sd e d u c e db a s e do nt h es t r a i n so ft h ee d g e s ．I tw a sd i s c o v e r e dt h a ta nF E Me l e m e n tc a nb ef u r t h e rs u b d i v i d e di n t ot h es m a l l e re l e m e n t so fs e a m e de d g e s ．S oan e wF E Mw a s p r o p o s e dc a l l e dt h ef i n i t ee d g ee l e m e n tm e t h o d F E E M ．B yu s i n ge d g e sa st h ee l e m e n t s ， F E E Mc a nc o n v e n i e n t l ym o d e lt h ef a i l u r ep r o c e s so fg e o t e c h n i c a lm a t e r i a l sa st h e ye v o l v ef r o m b e i n gc o n t i n u o u st ob e c o m i n gd i s c o n t i n u o u s ，f i n a l l yb r e a k i n gu p ．An u m e r i c a ls i m u l a t i o ns y s t e r nb a s e do nF E E Mw a sd e v e l o p e d ．T h es t r e s sa n dd i s p l a c e m e n to fa nu n d e r g r o u n dc a v e r n w e r ea n a l y z e da n dt h er u p t u r ep r o c e s so fad i s cw a ss i m u l a t e du s i n gt h i ss y s t e m ．T h er e s u l t s arec o l s et ot h o s ep r o d u c e db yo t h e rc o m m e r c i a ls o f t w a r e ． K e yw o r d s n u m e r i c a la n a l y s i s ；a r i t h m e t i c ；f i n i t ee d g ee l e m e n tm e t h o d ；f a i l u r ee v o l v e m e n t ； c o n t i n u o u s d i s c o n t i n u o u s ；b r e a k i n gu p 岩石破裂是岩土工程研究中最为关心的问题 之一[ 1 。4 ] ，数值求解破裂演化过程在岩土等力学工 程中具有广泛的应用价值．目前应用最广泛的数值 分析方法有限元法 F E M 和离散元法 D E M 邓。6 ] ，其对破裂问题的求解能力相差很大．有限元法 基于连续介质理论算法，适用于求解连续介质问 题，对破裂问题求解能力较弱；基于离散介质算法 的D E M 方法求解破裂问题能力较强，对连续介质 问题求解能力较弱，理论基础也较薄弱．而且，这两 种方法一般还必须事先知道可能的破裂面位置． 为了解决这些限制，很多学者对数值算法开展 了积极研究．1 9 9 4 年石根华提出的数值流形方法 收稿日期2 0 0 7 0 5 1 1 基金项目国家重点基础研究发展计划 9 7 3 项目 2 0 0 2 C B 4 1 2 7 0 3 作者简介宋宏伟 1 9 5 8 一 ，男，云南省玉溪市人，教授，博士生导师，工学博士，从事岩石力学与地下工程稳定性方面的研究 E m a i l h w s o n g l 6 3 1 6 3 ．c o r f l T e l 0 5 1 6 - 8 3 8 8 3 3 0 6 万方数据 中国矿业大学学报第3 7 卷 和非连续变形分析的D D A 方法口] ，用于分析非连 续介质的应力应变及破坏过程．1 9 9 5 年卓家寿和 章青等人提出了有别于传统有限元的界面元№j ，实 现了界面元与有限元方法的结合，可用于求解不连 续介质力学问题．1 9 9 9 年唐春安等提出了岩石破 裂过程分析系统R F P A 2 D 系统【⋯，可以模拟开采过 程中上覆岩层断裂和移动．2 0 0 4 年P O T Y O N D Y A 等提出了约束的颗粒元数值模型 B P M [ 1 ⋯，开发 了P F C 程序，可以广泛用于模拟岩石的破裂．2 0 0 4 年李世海等人提出了一种基于连续介质模型的离 散元方法 C D E M 口1 - - ，用质点和弹簧的集来代替 有限元的单元，可模拟破坏演化的过程． 本文提出了一种具有求解连续一非连续问题 的新数值分析方法，其主要特点是统一实现了有限 元和离散元计算功能，为岩土破裂演化等问题的求 解提供了一个新手段． 1 棱元算法刚度矩阵的推导 建立单元的刚度矩阵算法是有限元算法的核 心和关键，在建立这一刚度矩阵时，一般需要构造 一个单元插值形函数．在这里，我们不用形函数，而 直接采用棱参数来构造有限元的单元刚度矩阵． 以三角形单元 见图1 为例，其节点的全局坐 标为A z A ，y A ，B z B ，y B ，C x ∞Y c ，以三角形 单元的某条棱边建立线性局部坐标系押，则该条三 角形棱边的应变为 ￡。一 U ，一U ， ／l ， 1 式中“，“，分别为棱边的两结点沿其局部坐标系 以方向的位移；z 为这两点间的边长度． r 图1三角形单元及其棱边的局部坐标系 F i g ．1 T r i a n g l ee l e m e n ta n de d g e s ’ l o c a lc o o r d i n a t es y s t e m 翮一 划 t l lK l E l A l O t 3 lK3 E 3 A3 t 1 2K l E l A l 0 t 3 2K 3 E 3 A3 应变从全局坐标系转换到局部坐标系的变换 矩阵为 D c o s 2 卢s i n 2 卢C O Sp s i n 卢] 一[ 矗。] ， 2 式中c o s 卢一罕，s i n 卢一举． 假设在三角单元在全局坐标系下的应变为 o ，￡。，以。 ’，那么该棱的应变为 卜] [ e 。] 一[ d ，d d 。] l 勺| ． 3 1 只。J 三角形单元有3 条棱边，以式 3 同理，3 条棱 边的应变可以写成下列矩阵形式 睁鞋舞崩．㈤ 对式 4 求逆，得到以棱边应变表示的三角形 单元的应变为 讣雕薯r 卧㈤ 再由图1 有，位移从全局坐标系转换到局部坐 标系的变换矩阵为 l ’一c o s 』9 s i np ] 一[ ￡。] ． 6 那么对三角形单元的某棱边的应变为 e 。一 t l “J t 2 功一t l “。一t 2 口， ， 7 式中 “，让 ， “，，口， 分别为全局坐标系下该 棱2 个结点在z ，Y 方向的位移分量． 将三角形单元的某条棱边的刚度与棱边长度 有如下关系 』e 。一竽净卜扣≥』刚一竽衄净 l F K △“ lF K △“【 F K △“ K 一竽≥一丽K ， 8 则式 7 可以写成 ￡。一丽K f l “， ￡2 t 。- - t l u i 一￡2 u 。 ． 9 三角单元的3 条棱边的应变又可以写成下列 矩阵的形式 0 t 2 2K 2 E 2 A 2 t 3 2K 3 E 3 A 3 t 0 2一1，‘一l o KA“一“ 。 ㈣一聃一跳 KA趣一泌 。 仙一酏龇一跳 。 KA省一租 ， ㈨一弛捣一眦 。 万方数据 第1 期宋宏伟等具有求解岩土破裂演化问题能力的有限棱元法 将其代入式 5 得到 r ] 剧一 啊1 1d 1 2d 1 3 ] _ 1 d 2 1d 2 2d 2 3 L d 3 1d 3 2d 3 3 J 由式 1 1 知，三角形单元的应变矩阵为 同l ld 1 2 B 一1d 2 ld 2 2 L d 3 ld 3 2 一些一业 E l A l O E 1 A l O t l lK 1 E l A l t 2 lK 2 E 2 A 2 O t 2 lK 2 E 2 A 2 t 3 lK 3 E 3 A 3 t 1 2K 1 E l A l t 2 2K 2 E 2 A 2 垫垫00 E 3 A 3E 3 A 3 这样用棱边为单元推导的三角形单元的刚度 矩阵为 艇曲8 一J B D B ‰d 1 3 这里将这种由有限单元的棱边，即由2 个结点 加1 条棱，构成的单元的有限元分析方法称为有限 棱元法 F E E M ．在使用有限棱元法进行力学分析 时，首先用传统的有限元剖分的算法，将介质实体 进行有限元网格划分，再对各个单元进行棱元划 分．将棱单元看成是传统有限元单元的进一步划 分，即先划分成三角形单元、四边形单元，再以它们 的边作为棱元进一步划分． O t 2 2K 2 E 2 A 2 t 3 2K 3 E 3 A 3 O0 t 2 1K 2 E 2 A 2 t 3 1K 3 E 3 A 3 t 2 2K 2 E 2 A 2 t 3 2K 3 E 3 A 3 ． 1 2 2 有限棱元法的特点与验证 有限棱元法融合了有限元和离散元的优点，具 有从连续到非连续的分析能力．其特点如下 1 可以与有限元法一样方便地处理连续介质 的分析． 利用式 1 3 计算单元刚度矩阵的方法，我们开 发了2 维的数值计算系统R o c k B o x ，利用其与大型 商业软件A N S Y S ，F L A C 对同一地下硐室进行分 析，其分析得到的位移及应力分布图都具有相同的 分布规律 见图2 ． 图2 软件计算结果对比 F i g ．2C o m p a r i s o no ft h ea n a l y t i c a ls o l u t i o n s 表1 列出的对比计算数据显示，计算结果基本 一致．这说明棱元法的正确性，这也证明了数值分 析时，实体可以划分为有限元的三角形单元，还可 以进一步划分为棱单元． 表1软件计算结果对比 T a b l e1 C o m p a r i s o no ft h ea n a l y t i c a ls o l u t i o n s 2 可从连续介质分析开始，顺利过渡到非连 续介质问题的分析，即具有统一求解连续一非连续 力学介质问题、进行力学介质破裂演化模拟的能 力． 以棱单元来建立一个数值模型，当开始时棱单 元受力与变形处于弹性状态，介质处于连续状态； 当棱元的受力和变形不断增大，棱元发生断裂，根 据棱断裂情况不同，程序可以自动转到损伤计算、 非连续状态计算和破碎 离散 计算，而不需要进行 重新划分网格．从而统一实现了有限元与离散元算 法的计算功能． 利用R o c k B o x 对平面应变园盘进行模拟 弹 性模量1 0G P a ，泊松比0 ．2 ，抗拉压强度分别为 1 ．3 ，5 0M P a ，其结果表明 见图3 ，得到了从初始 损伤，破裂产生，到最后破碎的破裂演化过程． 们 饥 妣 伽 m 饥 蚝一A越一M 。 汹一聃批一聃 。 晡一A掩一M ， 仙一聃捣一跳 。 K A G b 泓一附 。 批一酏 堕A 鱼k 泓一蹦 。 ㈣一吼 ● ” 嬲 孙 d d d 万方数据 1 8中国矿业大学学报第3 7 卷 3 结论 图3R o e k B o x 模拟结果 F i g ．3M o d e l i n gr e s u l tb yR o c k B o x 1 利用本文方法导出了有限元的单元刚度矩 阵，证明了有限单元可以进一步划分到棱边一级， 以此为基础提出了一种新算法，即有限棱元法．该 方法以有限元单元的棱作为单元，算法上把有限元 与离散元很好地整合在一起，具有统一求解连续一 非连续介质问题的特点，为岩土介质破裂演化规律 的研究提供了一种新的研究手段． 2 基于有限棱元法，开发了程序，其计算实例 表明，模拟破裂演化具有较好效果，证明了棱元法 具有连续到破裂演化的统一分析能力． 3 要使棱元法能很好地用于模拟岩土破裂演 化等问题，还需对棱元法的岩土破裂判据等进一步 深入研究． 参考文献 [ 1 ] [ 2 ] [ 3 ] [ 4 ] S O N GH o n g w e i ，L US h o u m i n g ．R e p a i r i n go fa d e e pm i n ep e r m a n e n ta c c e s st u n n e lu s i n gb o l t ，m e s h a n ds h o t c r e t e l - J ] ．T u n n e l l i n ga n dU n d e r g r o u n dS p a c eT e c h n o l o g y ，2 0 0 1 ，1 6 3 2 3 5 2 4 0 ． S O N GH o n g - w e i ，Z H A OJ i a n ，J I AY i n gx u a n ．N u m e r i e a ls i m u l a t i o no nt r a n s v e r s ee f f e c to fb o l ti n d i s c o n t i n u o u sr o c km a s s [ C ] ／／L i nY u n m e i ．P r o c e e d i n g s o ft h e2 t hI n t e r n a t i o n a lS y m p o s i u mo nN e wD e v e l o p m e n ti nR o c kM e c h a n i c sa n dE n g i n e e r i n g ．P r i n c e t o n R i n t o nP r e s s ，2 0 0 2 2 7 7 2 8 0 ． 段艳燕，宋宏伟，赵坚．岩体非连续结构面对隧道 稳定性影响的数值分析E J ] ．岩土工程界，2 0 0 2 ，7 9 3 4 - 3 5 ，3 9 ． D U A NY a n y a n ，S O N GH o n e w e i ，Z H A OJ i a n ． N u m e r i c a ls i m u l a t i o no nE f f e c to fr o c kd i s c o n t i n u o u s s t r u c t u r eo ns t a b i l i t yo fT u n n e l [ J ] ．G e o t e c h n i c a l E n g i n e e r i n gW o r l d ，2 0 0 2 ，7 9 3 4 3 5 ，3 9 ． 宋宏伟，贾颖绚，段艳燕．开挖中的围岩破裂性质与 支护对象研究I - J ] ．中国矿业大学学报，2 0 0 6 ，3 5 l - s ] [ 6 ] [ 7 ] [ 8 3 [ 9 ] [ 1 0 ] [ 1 1 ] 2 1 9 2 19 6 S N GH o n g w e i ，JI AY i n g x u a n ，D U A NY a n y a n ． S t u d yo nc h a r a c t e r i s t i c so fr o c kb r o k e nb ye x c a v a t i o n a n dr o o fs u p p o r t i n go b j e c t [ J ] ．J o u r n a lo fC h i n aU n i v e r s i t yo fM i n i n g T e c h n o l o g y ，2 0 0 6 ，3 5 2 1 9 2 1 9 6 ． 蔡美峰，何满潮，刘东燕．岩石力学与工程[ M ] ．北 京科学出版社，2 0 0 2 1 1 9 - 1 2 1 ． C U N D A I 。I 。PA 。S T R A C KO DI 。．AD i s c r e t en u m e r i c a lm o d e lf o rg r a n u l a ra s s e m b l i e sl - J ] ．G e o t e c h n i q u e ，1 9 7 9 ，2 9 1 4 7 6 5 ． S H IG e n - h u a ．D i s c o n t i n u o u sd e f o r m a t i o na n a l y s i s a n e wn u m e r i c a lm o d e lf o rt h es t a t i ca n dd y n a m i c so f d e f o r m a b l eb l o c ks t r u c t u r e s [ J ] ．E n g i n e e r i n gC o m p u t a t i o n s S w a ns e a ，W a l e s ，1 9 9 2 ，9 2 1 5 7 - 1 6 8 ． 章青，周资斌，卓家寿．分区界面元一有限元一无限 元混合模型[ J ] ．计算力学学报，2 0 0 5 ，2 2 1 8 - 1 1 ． Z H A N GQ i n g ，Z H O UZ i - h i n g ，Z U OJ i a - s h o u ． M i x e dm o d e lf o rP I S E F E M i n f i n i t eE M l - J ] ．J o u r n a l o fC o m p u t a t i o n a lM e c h a n i c s ，2 0 0 5 ，2 2 1 8 - 1 1 ． 唐春安，徐曾和，徐小荷．岩石破裂过程分析R F P A 2 0 系统在采场上覆岩层移动规律研究中的应用 E J ] ．辽宁工程技术大学学报自然科学版，1 9 9 9 ，1 8 5 4 5 6 4 5 8 ． T A N GC h u n - a n ，X UZ e n g h e ，X UX i a o h e ．A p p l i c a t i o no fa n a l y s i ss y s t e mR F P A 2 0o fr o c kf a i l u r e p r o c e s si nr e s e a r c h i n gm o v i n gr u l eo fc o v e r i n gw o r k f a c e [ J ] ．J o u r n a lo fL a o n i n gT e c h n i c a lU n i v e r s i t y N a t u r a lS c i e n c e ，1 9 9 9 ，1 8 5 4 5 6 - 4 5 8 ． P O T Y O N D Y AD0 。C U N D A I 。I 。PA ．Ab o n d e d p a r t i c l em o d e lf o rr o c k [ J ] ．I n t e r n a t i o n a lJ o u r n a lo f R o c kM e c h a n i c s M i n i n gS c i e n c e s ，2 0 0 4 ，4 1 1 1 3 2 9 1 3 6 4 ． I 。ISH ，Z H A MH ，W A N GYN ，e ta 1 ．An e w c o m p u t a t i o n a lm o d e lf o rD E M b l o c ka n dp a r t i c l e m o d e l [ J ] ．I n t e r n a t i o n a lJ o u r n a lo fR o c kM e c h a n i c s M i n i n gS c i e n c e s ，2 0 0 4 ，4 1 3 4 3 6 ． 责任编辑王继红 万方数据