（水工结构工程专业论文）数值流形方法及其工程应用研究.pdf

摘要 实际工程中普遍存在连续与非连续变形问题，其数值模拟方法呈现两极分化的态势。 一方面是以有限元方法( f e m ) 为主的、基于完全连续或等效连续假说的连续性分析方 法；另一方面则是以非连续变形分析( d d a ) 与离散单元法( d e m ) 为代表的基于块 体理论、离散介质假说的完全非连续性分析方法。数值流形方法( m 心i ) 采用有限覆 盖技术建立了连续与非连续问题的统一求解格式，得到了学术界的广泛关注。本文深入 研究数值流形方法，在进一步完善和改进的基础上，探讨其在水利工程中的应用技术。 主要成果如下： ( 1 ) 在石根华早期数值流形方法计算程序的基础上，完善了边界约束的处理方法， 并采用面向对象程序设计( 0 0 p ：0 b j e c t e d 一0 i _ i e n t e d p m 盯a l l l m i n g ) 方法，运 用s u a lc + 十研制了w i n d o w sx p 平台上的可视化数值流形方法计算程序， 并用算例予以验证。 ( 2 ) 详细比较了数值流形方法和有限元方法等的相互关系，探讨了数值流形方法 权函数的形成过程和特点，以及对接触问题的处理方法。以加权余量法为基 础，研究了n n m 、f e m 、e f g m 和d d a 等的试函数构造方式。 ( 3 ) 运用所研制的数值流形方法计算程序，分析重力坝应力状态，并与有限元法 计算结果进行比较，证明计算结果是合理可信的。应用数值流形方法计算大 变形的能力，采用降强法深入分析了重力坝的滑动变形规律，建立了确定坝 基面极限内摩擦角的准则，提出了重力坝抗滑稳定分析的数值流形方法。 ( 4 ) 考虑非稳定渗流的影响，运用所研制的数值流形方法计算程序，深入研究了 某水库水位下降产生的非稳定渗流场对水库滑坡体稳定性的影响。 关键词数值流形方法有限覆盖技术连续一非连续重力坝应力分析抗滑 稳定边坡稳定降强法 a b s t r a c t t h ep r o b l e m so fc o m i n u o u sa i l dd i s c o n t i n u o u sd e f 0 咖a t i o ni na c t u a le n g m e e r i n 卫a r e u b i q u i t o u s ，t h ec 时e n tn u m 砸c a lm e t h o d si nd e a l i n g 谢山i t 盯es u i t a b l ef o re i l e ri d e a l c o n t i l l u 哪m a t e r i a l s ，o rc o m p l e t c l yd i s c o m i n u o l i sm e d i a o nt 1 1 eo n eh a n d ，s u c hm e t h o d sa s f b ma n db e ma r eb a s e do nm eh y p o m e s i so f t o t a lc o n t i n u i t yo re q u i v a l e n tc o n t i n u i t 矿o nt l l e o t h e rh a n d ，t h e r ea r ed i s c o n t i n u o u sd e f o 蛐a t i o na i l a l y s i sm e m o d ss u c ha sd d aa i l dd e m w h i c ha r eo n l eb a s i so fm eb l o c k ym e o r ya n dt h ea s s u m p t i o no fd i s c r e t em e d i a t 1 l e n l l m e r i c a lm a n i f b i dm e m o dh a sb e e ne s t a b l i s h e de l e m e n t a r i l va n dr e c e i v e dw o r l d w i d e a _ c t e 而o n s mt h ep a p e r ，m en m mi sl u c u b r a f e d ni n c l u d e st h ef 0 1 l o w i n gc o n t 捌1 t s ： ( 1 )t h ee a r i yn m m p r o 掣锄w r o t eb ys h ii nd o sh a sb e e ni m p m v c di nw i n d o w s o s ，a 1 1 dw h i c hi sp m v e db ys o m ee x a l p l e s ( 2 )s o m en l l i n e r i c a im e t h o d sa r ed i s c u s s e d ，s u c h 船f e m ，e f g m ，d d ae t c o nt h e b a s eo fm e t h o do fw e i 出e dr e s i d l l a l s ，t 1 1 e 仃y 劬c t i o na n dw e i 曲t e df l l i l c t i o no f n m ma r es t u d i e d ( 3 ) t h e a 1 1 a l y s i s o f s 订e s so f 争a v i t y d a m i ss t u d i e d b y n m m ，a n d t l l er e s u l t i s p m v e d b yf e m 眦g 1 ca b i l i t yo f c a c u l a t i n gl 哪ed e f o m a t i o no f n m m ，t h eo r d e ro f s l i p p a g ed e f o m a t i o no f 铲a v i t yd 枷i ss t l l d i e db ys a f e t yr e d u c t i o nm e l o d ，a n d t j l em e t h o dt oo b t a i nt 1 1 el i m i t e dm c t i o na n 甜eo fi m e 商a c eb e t w e e nd a ma l l d f o u n d a t i o ni se s t a b l i s h e d ，弛dm es 诅b i l 姆a g a i n s ts i i d i n go f 田v 姆d 锄i s s t i l d i e db yn m m ( 4 ) t 被i n go n er e s e r v o i rr o c ks l o p sa st h ee x a i i l p l e ，t l l er o c ks l o ps t a b i l i 锣i s d j s c u s s e d ，w h e nt h er e s e r v o i rw a t e rl e v e lf 矾l s k e yw o r d s ： n i l i l l 嘶c a lm a i l i f 0 1 dm e t h o d ；f i “t ec o v e rt e c m q u e ；c o n t i n u o u s d i s c o n t i i l u o u s d e f b 玎n a t i o n ；g r a v i t ) rd 棚；柚a l y s i so fs c r e 蹦a 1 1 a l y s i so fs t d b i l 时a g a i l l s t s l i d i n g ；s i o p es 协b i l i t y ；s a f e t yr e d u c t i o nm e t h o d 学位论文独创性声明： 本人所呈交的学位论文是我个人在导师指导下进行的研究工作及取得 的研究成果。尽我所知，除了文中特别加以标注和致谢的地方外，论文中 不包含其他人已经发表或撰写过的研究成果。与我一同工作的同事对本研 究所做的任何贡献均已在论文中作了明确的说明并表示了谢意。如不实， 本人负全部责任。 论文作者( 签名) ： 学位论文使用授权说明 河海大学、中国科学技术信息研究所、国家图书馆、中国学术期刊( 光 盘版) 电子杂志社有权保留本人所送交学位论文的复印件或电子文档，可 以采用影印、缩印或其他复制手段保存论文。本人电子文档的内容和纸质 论文的内容相一致。除在保密期内的保密论文外，允许论文被查阅和借阅。 论文全部或部分内容的公布( 包括刊登) 授权河海大学研究生院办理。 论文作者( 签名) ：- 埠_ _ 盔l a 卯绰，月 乃 日 河海大学硕士学位论文 1 1 引言 第一章绪论 岩体是经过漫长的时间而形成的，经历过变形、遭受过破坏的地质体【l 】。从总体上 来说，岩体是一种非连续介质，是由许多不同走向不同规模的断层面、节理面、裂隙面、 层面等各种地质界面切割成的、极其复杂的介质体。岩体结构对于岩体性质具有很大的 影响，岩体结构的复杂性导致岩体呈现出强烈的非连续性、非均质性及各向异性等特点。 水利工程中也存在着各式各样的连续或非连续问题，如混凝土结构施工接缝和建筑基础 与地基接触造成的非连续问题。实际工程中所涉及的岩体和水工建筑物结构的变形一般 首先是损伤演化、小变形积累，到宏观裂纹的形成与扩展，最后导致结构体的大变形、 大位移的一个复杂的渐进变形与破坏过程。为了求解这类数据有限( d a t al i m i t c d ) 、了 解程度有限( u n d e r s t a n d 血gh m i t c d ) 的复杂结构体系，往往采取一些简化的措施分析结 构体的变形：当结构体内的应力未超过某一临界值时，结构体不会产生明显的断裂，宏 观上呈现出连续性，此时将其抽象为连续介质模型f 或等效的连续体模型) 进行数值分 析；而当结构体内的应力超过这一临界值时，不连续面的控制作用变得十分明显，呈现 出强烈的不连续性，此时应该将其抽象为离散体系统模型进行数值分析。总的来说，实 际工程中连续与非连续问题是普遍、客观存在的。 解决实际工程问题的方法主要有实验研究、理论分析与数值模拟三大类，实际工程 问题的整体结论依赖于三源知识的交叉、结合与相互作用。各种新的理论分析方法在分 析实际工程问题时晟终归结为与计算技术结合的数值模拟，一些新的实验研究方法也是 最终与数值模拟结合才在工程应用中得以体现。数值模拟技术己经成为种经济有效的 辅助手段，因此对于实际工程中数值模拟方法的研究是必须的和重要的。刘怀恒【2 】、周 维垣【3 】、王泳嘉与冯夏庭f 4 j 、奕茂田【5 1 、李宁叼等对实际工程中各种数值模拟方法进行 了系统的阐述。 1 2 数值模拟方法简介 近几十年来，在众多学者的共同努力下，实际工程中的数值模拟方法得到了极大的 发展，出现了多种数值计算方法，下面对各种方法进行简单介绍。 ( 1 ) 有限差分法 有限差分法( f d m ：f i n i t e d i f f e r e n c e m e m o d ) 【7 】是一种出现较早的数值计算方法。 它的产生是由于多数偏微分方程只有在简单边界条件情况下才能用严格的方法求解，而 对于一般情况无法得出精确的解析解。差分法即是将微分方程用相应的差分方程来代替 的数值方法，其基本思想是将基本方程和表示边界条件的方程近似地用差分方程来代 替，这样，就将求解微分方程的问题转化为求解代数方程的问题，而后者不会遇到数学 上的困难。有限差分法虽然方法简单，但是精度不高，已逐渐被有限元法所代替。 河海大学硕士学位论文 分片的刚体位移模式逼近实际整体位移场，以块体问联结弹簧反映结构内部的弹性，并 用界面应力表征结构内部的应力。刚体有限元法对于连续状态的应力分析可给出较高的 计算精度，对于临界状态能够估算出体系的极限载荷，并可有效地进行少量块体界面问 的摩擦接触分析。但目前所发展的刚体有限元法仅能够对小位移、小变形情况下的结构 进行分析计算，不能充分反映非连续介质的大变形与失稳过程。 ( 7 ) 非连续变形分析 非连续变形分析( d d a ：d i s c o n t i n u o l l sd e f o 衄a t i o na n a l 粥i s ) 【1 3 】是石根华与 g o o d m a t l 继提出岩体稳定性分析的关键块体理论( k b t ：k e yb l o c kt h e o r y ) 1 2 2 】之后， 2 0 世纪8 0 年代末发展起来的一种模拟散体系统大变形的数值分析方法。叻a 方法是 在块体运动学基础上，并部分吸收离散单元法关于接触形式的描述方法和刚块一弹簧模 型中位移函数的构造方法等方面优点发展起来的一种与有限元法并行的数值分析方法。 d d a 方法以被有限的不连续面切割的离散块体为基本单元，通过块体之间的接触和作 用在各块体上的位移约束条件，把若干个单独的块体连接起来并构成一个块体系统。利 用最小势能原理建立单个块体的单元刚度矩阵及块体系统的整体平衡方程来求解块体 系统的变形与应力。由于非连续变形分析方法采用统一的有限元格式求解刚体位移和块 体变形，所以不仅可以计算破坏前的小位移，也可以计算破坏后的大位移，不仅允许块 体自身有位移和变形，而且允许块体间有滑动、转动、张开等运动形式，对求解滑动、 崩塌、爆炸及贯入等问题十分有效。 ( 8 ) 数值流形方法 数值流形方法( n m m ：n l l r n c r j c a lm a i l i f o l dm e t l l o d ) f 1 4 是一种新兴的数值计算方 法，是石根华于1 9 9 1 年提出的。该方法吸收了有限元法和d d a 方法的优点，利用现 代流形分析中的有限覆盖技术，对连续问题和非连续问题建立了种统一的求解格式， 是一种相当有发展潜力的数值方法。 ( 9 ) 无单元法 自2 0 世纪7 0 年代起，国内外许多学者致力于)无单元法自2 0 世纪70年代起，国内外许多学者致力于探索一类不需要生成网格或单元，仅用系 列点对求解域进行离散的无单元法( e l 咖e n t - f r e 刮g r i d l e s s 舢e s h l e s sm e n l o d s ) 1 5 羽。但是直到进入1 9 9 0 年，滑动最小二乘方法( m l s a ：m o 、，i n gl e a s ts q u a r em 舶o d ) 应用于边值问题的求解，才标志着无单元法的真正诞生。 当然还 有许多其它方法，如极限平衡法、半解析元法【2 8 】、无限元法( 埘i n i t ee l 锄吼tm e m od ) 1 2 9 1 、拉格朗日元法( l a 舭n 西趾e l e m e n t metllod)【30】等。所有这些方法的不断融合、 发展都将使我们更好的解决工程实际问题。1 3数值流形方法的研究现状 河海大学硕士学位论文 覆盖权函数与覆盖函数构成了数值流形方法总体逼近函数的两个主要组成要素，下 面分别对其进行讨论。 2 4 1 覆盖权函数 根据权函数的一般定义，数值流形方法的覆盖权函数应满足单位分解性质，即 ( f ) ( x ) 0石也( n ( 2 2 ) 心m ( x ) = ox 仨也( ，) ( 2 3 ) ( x ) = l ( 2 4 ) t e 吒m 覆盖权函数在其拓扑星处取最大值1 ，在其对应的数学网格的边界上或边界外取为 0 。 用三角形有限元网格转化为数学网格的数值流形方法，流形单元的节点总是三个， 其权函数可按照三角形有限单元的形函数进行构造。如图2 7 所示，记流形单元e 的三 个物理覆盖所对应的数学网格的拓扑星为，( r = f ，j ，砷，拓扑星的坐标为 ( ，”) ( ，= f ，后) 。为了与节点相区别，下文称围成流形单元的点为“角点”。 形式 其中 匿图 为流形单元 “露 为流形单元节点 。 为流形单元角点 图2 7 三角形有限单元覆盖流形单元节点和角点 在流形单元上，线性覆盖权函数k 。( ，= f ，七) 可用拓扑星的坐标表示为以下一般 心= q + 印算+ c ，y ( 2 - 5 ) 按照三角形有限单元法构造形函数的方法可导出如下形式的覆盖权函数 心：( q + 4 x + 。y ) ( 2 6 ) = 五( q + 包工+ o 2 - 6 ( 2 7 ) 一般情况下，是的意注要需倍 河海大学硕士学位论文 第二章数值流形方法的基本理论0 3 1 ，4 5 ，6 7 j 本章根据搜集到的国内外资料，阐述了数值流形方法的基本概念，较系统地分析整 理了数值流形方法的基本理论。 2 1 数值流形方法的基本流程 数值流形方法采用数学与物理双重网格将整个研究区域剖分成有限个相互重叠覆 盖的集合，然后在各个覆盖上独立定义局部位移函数，最后通过权函数将各个局部位移 函数联结在一起，并根据最小势能原理构造出整个求解域上的总体控制方程( 平衡方 程) 。数值流形方法的基本流程如图2 1 所示。 图2 1 数值流形方法基本流程 2 2 数值流形方法的一般有限覆盖系统 “流形”词来源于拓扑学。拓扑学就是研究图形在拓扑变换下不变性质的学科 嘟j ，流形是其重要研究对象。一个流形就是一个局部欧氏的拓扑空间。数值流形方法 中的“流形”与传统的流形不同。传统流形中的全局函数是高度可微的并且与覆盖无关， 而数值流形方法的全局函数是定义于覆盖上的、分片可微的，且在接触面上可以是非连 续的。 一个工程问题的求解区域在数学上可视为一个流形，它的响应( 如位移) 就是流形 的一个变换。对于某一流形进行复杂变换，通常可以将其分解为一些简单的几何图形如 三角形或多边形，然后用一些易于分析的图形来覆盖这些简单的图形。这样，对一个复 7 第二章数值流形方法的基本理论 这个面积并不等于流形单元以的面积。 式( 2 6 ) 与三角形有限元的形函数在形式上完全一样。只是对于有限元而言，式 ( 2 6 ) 中的是有限单元面积的2 倍，而在流形单元中，它是由三个覆盖的拓扑星所 组成的三角形的面积的2 倍，而不是流形单元的面积的2 倍。由以上分析可知，对于有 限单元覆盖流形方法，覆盖权函数实质上就是有限元的形函数。那么，如果用四边形有 限元网格作为数学网格的流形方法，流形单元的节点总是四个，则权函数可以选取四边 形有限单元等参变换的插值函数。 2 4 2 覆盖函数 数值流形方法的覆盖函数可表达为任意函数的形式，如定义覆盖函数u 。为 也( n = ( 。) 见( i ， ( 2 8 ) 其中乜为覆盖p ( 力的覆盖函数的未知数向量，e ( f ) 可采用多项式或非多项式等不同形 式。 对于二维流形单元口的物理覆盖“f ) ，其覆盖函数可写成如下形式 恐跚= 种圳躲) 协， 或表示成 ( f ) ( 五y ) = 正( f ) ，( 工，_ y ) 以( f ) 2 h = e ( f ) 7 d 美，) ，= l ( 2 1 0 ) m f n k ( f ) ( 五y ) = 厶哪( 五y ) 吐( f ) ，：，= e 以。) 式中m 9 ) 表示覆盖函数的未知系数个数之半( 三维流形单元为三分之一) ，f ) 、f ) 为 覆盖p ( f ) 的未知数向量( 广义自由度) ，碚、瓴y ) 为基向量a 对于多项式形式的覆盖函 数，基向量瑶、( x ，y ) 的形式决定覆盖函数的阶数。 一阶：只7 = ( 1 ，五力 ( 2 - 1 1 ) 二阶：e “) 7 = ( 1 ，工，弘x 2 ，秽，2 ) ( 2 - 1 2 ) 三阶：e = ( 1 ，工，y ，z 2 ，掣，y 2 ，z 2 y ，叫2 ，y 3 ) ( 2 - 1 3 ) 相应地，未知系数个数也随覆盖函数阶次的增加而增加。 2 5 本章小结 本章阐述了数值流形方法的基本流程，介绍了数值流形方法的有限覆盖系统、覆盖 函数和覆盖权函数。 第三章数值流形方法的基本方程及分析过程 第三章数值流形方法的基本方程及分析过程 3 1 , 4 5 , 6 7 1 本章从数值流形方法的基本理论出发，推导了数值流形方法的基本公式，列出了有 限单元覆盖的单元矩阵和接触矩阵，阐述了数值流形方法的积分方法，并研制开发了相 应的数值流形方法计算程序。 3 1 数值流形方法的总体控制方程 议分研城门伺n7 f 、别埋镘盂，母1 、况彬早兀伺q1 、咧理後盂ln q r 侃彤早兀j ， 每个物理覆盖有2 m 个未知系数( 广义自由度) ，则物理覆盖函数为 髅辨水7 o 川 阮1 ) 净- ， 流形单元上的总体位移函数为 鼢辨轧c 训uo ( x , y ) = 窆i = 1 芝y = ll l k 。予m 。v ，e ( o ：je ( i ) j 之篡2 卜j 1 ” jl “e ( f ) ，j = 喜 c 2 ) = 阢】 皿 ( 3 2 ) 其中 蚋 = 心“乎。”心。，羔ik 乎”比。，羔。 ，o = ，2 ，c s s ， 皿) = 【z ( y ) 】= z ( 。j ( 五y ) 正( ：) ( 工，y ) ( 。) ( 五y ) ( 3 5 ) n - 乜。 q ( 2 1 ( 3 - 6 ) t 、见分别为流形单元e 的插值函数矩阵和未知系数向量。 在得到流形单元上的总体位移函数后，就可建立弹性力学边值问题中的能量泛函表 耻 薹 加“咖；缈枷 河海大学硕士学位论文 达式。系统的总势能为 兀= 艺( 兀。+ 1 i 。+ n p + n 。+ n ，+ n ，) 昙( d 1 7 砬7 乜7 ) + ( d 1 7d 2 7 见7 ) 墨。墨： 如。： e 巧： 弘 置。 。 ： 巧。 q d 2 ： e ( 3 - 7 ) 其中：1 7 。为单元应变能；。为初应力势能；丌。为点荷载势能；。为体荷载势能；兀， 为惯性力势能：兀，为约束形成的势能。因为每个物理覆盖有2 m 个未知系数，由式( 3 7 ) 给出的系数矩阵的子矩阵i 岛l 是一个2 m 2 m 的矩阵。 d f 】和k 】是2 m 1 的子矩阵，d l 表示物理覆盖i 的2 m 个未知系数( 节点的自由度) 。 根据最小势能原理，所有几何可能位移中，真实位移使系统总势能取极小值，则有 o r t 等= o ，( r = 1 ，2 ，2 m ) ( 3 8 ) o a i ， 得到的联立平衡方程式( 总体控制方程) 如下 足i i 墨2 岛。： 芷。l 量。2 正 五 ： e ( 3 - 9 ) 其中 巧 的r 行s 列元素为 孚 ( 删，2 ，2 m ) (310)0 d ，a a t jij 。 鼻的，行元素是n 为零( 即所有变量d i ，见等于零) 时的导数，即 a n ( o ) j 峨】。 p = i ，2 ，2 m ) ( 3 一1 1 ) 数值流形方法总体控制方程的建立、总体刚度矩阵的组装过程与有限单元法完全相 同，它的基础仍然是最小势能原理。不同的是数值流形方法能量泛函的积分求解是在流 形单元 1 进行的，冈此存数值流形方法巾凡涉及面积分求解的均是指在流形单元上的面 积分。流形单元的形状允许很不规则，甚至是凹的。当流形单元为l a g r a n g e 型单纯形 单元时，单元刚度矩阵可运用单纯形积分解析求解。 4 吒：k 第三章数值流形方法的基本方程及分析过程 在具体计算时，数值流形方法的计算是按时间步逐步计算的。对每一时间步，要进 行开合迭代阻使当前系统在所有接触界上满足无嵌入与法向弹簧不受拉两个条件。数值 流形方法使用的总体方程迭代算法是超松驰迭代法。 3 2 有限单元覆盖的单元矩阵 3 2 1 单元刚度矩阵 单元刚度矩阵的建立与有限元类似，首先由弹性力学公式 盯) = 【e 】 s ) = 【e 】【l 】 址 = 占】【三】【e 】 皿) = 陋】阻】 晓 ( 3 1 2 ) 其中 m 【e 】为弹性矩阵，对于平面应力问题 旦o 唧 。旦 泖 aa 砂缸 ：三 。1 1 一l ，2 1 l ， yl 0 0 o o 1 一p 2 ( 3 1 3 ) ( 3 1 4 ) 流形单兀上产生的应变能。为 n 。= 吉( t 吒+ e v f y y + r 。) d x d y = 扣) 7l 孵九州趣 蛐p = 去 见) 【蜒】 见 ( 3 1 5 ) 其中【疋】为单元刚度矩阵。 3 2 2 初应力矩阵 数值流形方法在每一时步开始计算时，初应力为上一时步末的应力状态，在最开始 的时候，其初应力为外力( 包括外荷载、体力、温度应力等) 。 单元e 上的初始应力为 河海大学硕士学位论文 。) = 盯( q + 勺嘭+ 焉协毋 = l 【，卜 i o 堋a 吖7 即0 淳忉 i +i 删a 吲7 即 b l ，oi 可以叠加到荷载矩阵。 加荷点可以在单元内任何地方，设有点荷载力f = 【e r 作用在单元e 的点 p ( x o ，y o ) 处，点p ( x o ，y o ) 的位移为 = u ( x o , y o ) n ，= - ( f ，u ( x o ，) + f ，( x o ，) ) 叫川x o , y o ) 1 鼢 _ _ 吖州甜 z 。， 根据最小势能原理可以将点荷载矩阵 噼删7 鼢 、( 3 - 2 1 ) 6 _ 3( -，fj一畔，一聋l，、l = l_-tf畔。，iji 为能势的 l l jp 元啦则 笙三兰塑堕堡丝塑堡箜量查查墨墨坌堑垫堡一 3 2 4 体荷载矩阵 设作用在单元e 常体应力为 正 7 ，则由体荷载产生的势能为 其中 舡l 粉y jl 。j 妫 ( 3 - 2 2 ) ( 3 - 2 3 ) 为体荷载矩阵，叠加到荷载矩阵中 若体力不为常量，则体荷载矩阵变为 枷，7 协蚴 ( 3 - 2 4 ) 3 2 5 惯性力矩阵 惯性力矩阵等效于有限元的质量矩阵，该矩阵是数值流形方法最重要的矩阵，若给 一小时间步，惯性力矩阵将控制整个材料体的所有点的运动，且不连续接触的计算将是 稳定的。 在当前时间步，设【“( f ) “f ) r 为单元e 的任意点“力与时间有关的位移，m 为单 位面积的质量，则单位面积的惯性力为 件一略吲等 。2 5 ) 单元e 的惯性力势能可表示为 n 。= 一j j pv 】 2 出咖 ：f f 盹唰掣蚴 ( 3 2 6 ) 设 皿( o ) ： o 是时间步开始时的位移， 砬( ) j = 见 是时问步结束时的位移， 为时间步长，忽略3 阶及以上项，则展开式为 见( ) ) ： 见( o ) ) + 型8 t + 譬掣 一 趔+ 尘i ：i 堡亟( 3 - 2 7 ) 魂 2西2 讲而 、l，j 、，j，b鼢帐 栅 啪 删 脚 力汀 o 舭。 河海大学硕士学位论文 驾掣= 素 删) 一云掣= 素 必) 一和。) c 。 式中 御 = 挈为起始时间步的速度所以势能为 毕m 阶n 】蚴旧讣扣o ) ) ( 3 - 2 9 ) 因此 惯性力矩阵为 一1 9 吲7 吲嘞j ( 万2 m 一警o ) ) ( 3 _ 3 。) 其中质量矩阵 群 = 警岫他】蚴j 协， 叠加到刚度矩阵，速度矩阵 f = 警九互】妫扣o ) ( 3 - s z ) 叠加到荷载矩阵。 3 2 6 约束矩阵 设单元上的点p ( x o ，y o ) 受约束，可在此点沿x 、y 两个方向加设刚性弹簧来约束它， 其刚度分别为见、p y ，则约束力可表示为 阱 ：德舄 ， 弹簧应变能可表示为 n = 互1 ( n “( 而，y o ) 2 + p y v ( 而，) 2 ) = 圭 几”( 而，儿) p , v ( x o , y o ) j 1 q u 。( ( x o ，, y o ) ) = 三 见) 7 阢( ，蜘灌l 苫兰l 【z ( ，虬) 见) ( 3 3 4 ) 则约束矩阵为 科h 叫心跏( ) 】 ( 3 _ 。s ) 叠加到刚度矩阵。 当点p ( x o ，y o ) 受x 、y 两个方向约束时，取p ，= 矽、，= p ；当点p ( x o ，y o ) 只受工方向 第三章数值流形方法的基本方程及分析过程 约束时，取n = p ，p y = o ； 点p ( x o ，儿) 只受y 方向约束时，取以= o ，p ，= p 。其 中p 是一个相对于材料的弹性模量来讲很大的数，一般可取1 0 3 1 0 5 e 。 3 3 有限单元覆盖的接触矩阵 3 3 1 接触判断 数值流形方法通过在不连续面两侧采用不同的覆盖函数作为场函数的局部近似，刻 画了不连续面处整体函数的非连续性。对于非连续问题，块体的运动不允许它们之间受 拉伸和相互嵌入，无拉伸和无嵌入的直接数学描述是用组合的不等式给出的。当块体系 统运动或变形时，块体只是沿块体边界接触，而无嵌入的不等式可通过施加刚度很大的 弹簧锁住块体在一个或几个方向的运动而转换为两个块体接触时的等式，这些等式将加 于总体平衡方程上。如果两个块体在他们之间产生拉伸在锁定弹簧撤离后它们将分开。 因此，无拉伸和无嵌入约束可简化为施加或撤离弹簧。在计算时，必须反复求解总体平 衡方程，即可修正系统中有拉伸和嵌入的块体。 为了计算块体系统各个接触面上联结弹簧贮存的变形能，首先必须判断各个块体之 间的接触关系。 由于块体位置、块体形状和块体之间的接触随时间步变化，所以我们必须知道在下 一时间步所有可能接触的成对块体。如果两个块体发生接触，那么它们之间的最小距离 应小于某一极限值。 设点a ( ，m ) 、p 2 ( x 2 ，y 2 ) 分别是块体4 和b 上的任意点，那么接触距离为 r ? d = ( 屯一五) 2 + ( y 2 一y 1 ) 2 ( 3 - 3 6 ) 只有当d 0 时，块体4 和曰是接触或重叠的，而当d 0 时，块体爿和b 之间是分离的。 二维情况下，数值流形方法用一个块体所有的角、边( 面) 和另外一个块体的角、 边的距离来确定块体之间的接触关系。对于三维情况，计算量很大。 3 3 2 法向接触矩阵 如图3 1 所示，设p l 是顶点，p 2 p 3 是进入线，法向弹簧垂直p 2 p 3 。 ( a ) 接触前 p 图3 - 1 法向接触示意图 1 9 ( a ) 接触后 一塑塑茎兰堡主兰些堡苎 若单元碍口，之间的法向接触距离为d ，进入线长为，法向弹簧剐度为以，晶为 进入点与进入边围成的三角形的面积，则法向接触距离为 d = 十斜十孚( 3 - 3 7 ) 其中 h i = 知训) 7y 2 - y 3 l 任s s ， g = ； 巧( ，儿) 7 。y 五 一- 毛y i + ； i ( 而，乃) 7 。y 恐t - 一y 2 c ，s 9 ， 法向接触弹簧的势能为 f i = 譬扩= 譬 日) 7 h ) 日) 7 q ) + b 7 g ) g ) 7 哆 + 她 7 ( 日) g 7 9 + 埘f g 日 7 咖：( 和口) 7 + 2 ( 争埘 g ) 十( 争2 ( 3 - 4 0 ) 其中b 日) 日) 7 、只 日) g ) 7 、只 g 口) 7 和见 g ) g ) 7 叠加到剐度矩阵，而 嘿( 争 日) 和也( g 叠加到荷载矩阵。 3 3 3 切向接触矩阵 如图3 2 所示，设p i 是顶点，p 2 p 3 是进入线，p o 是在p 2 p 3 上与pj 假定的接触点， 剪切弹簧在p 2 p 3 方向上。 、 ：二：，=，：j：，7：， 三二当! ，二二一 f。，二【l一 p 2 p o 。p3 图3 - 2 切向接触示意图 若假定的接触点p 0 和p t 沿p 2 p 3 方向的剪切位移为j ，进入线长为，弹簧切向刚 度为以，则单元f 和_ ，2 _ n n 剪切位移为 d 。：霹e 丽 2 7 1 f 毛+ 川一而一一十h 一一v 。】f 羹：瓮：芰：瓮 c s 4 ， 第三章数值流形方法的基本方程及分析过程 瓯= 一一y 。一儿 y x 3 3 - _ y x 2 z 因p 0 和p l 之间距离很小，所以 d z s ，o ，i 。f 为一屯儿一虬】 麓：翟) = 一s ，o ，i 。1 屯一鼍儿一儿】 i + ；【t 一而 可记 d = 7 d f ) + 斜+ 孚 ( 3 - 4 2 ) ( 3 - 4 3 ) ( 3 4 4 ) = 癣川) 】馁：羔 悖4 s ， g = ； 弓( ，乩) 妻： c ，4 e ， 剪切弹簧的势能为 兀。= p - - ，- “d 2 = 譬i 口 7 日 7 皿 + 哆 7 g g q + 口 7 何 g 7 q + 哆 7 g 日) 7 d j + 2 ( 孚) d f 7 日 + 2 争 q g + ( 孚) 2l ( 3 _ 4 7 ) 其中n h 日 7 、p o ( h g 7 、以 g 日 7 和以 g g 7 叠加到刚度矩阵，而 一岛粤) 日 和一只粤) g 叠加到荷载矩阵。 3 3 4 摩擦力矩阵 如图3 - 2 所示，设p l 是顶点，p 2 p 3 是进入线，p o 是在p 2 p 3 上与p l 假定的接触点。 当单元f 和j 在接触边界上发生相对滑动时，如果摩擦角不为0 ，滑动面上存在摩擦力， 摩擦力是根据法向压力计算的，即 f = p d s t a ng o ( 3 4 8 ) 其中p 为法向弹簧刚度，d 为嵌入距离，t a n t p 为摩擦系数，s 为符号函数，由摩擦力方 向确定。 单元f 上p l 处的摩擦力f 的势能为 n ，= 了f 。h ， 妻：芰) = 孚 口) 7 豫c 五，乃，】7 羔：羔 = ， 口 7 日 c ，4 ， 2 l l，、r】 ，、【 1 _ j 肛 一 耽 河海大学硕士学位论文 式甲 = 知” ( x 3 - x 2 ： b s 。， 其中- f h 叠加到荷载矩阵。 单元_ ，上p o 处的摩擦力f 的势能为 ，= 一手c v o ， 羔：主) = 一芋 b 7 弓c ，蜘， 7 妻：羔 = 一f 哆) 7 g ( 3 5 1 ) 式中 g = ； 弓( ，) y 屯3 - 一屯y ： c s - s z ， 其中f g j 叠加到荷载矩阵。 3 4 积分方法 数值流形方法的算法收敛性和精确性主要依赖于复合形形状的解析法积分，为此石 根华在数值流形方法中引进单纯形积分方法，单纯形积分是弗维一般形状区域的精确 解，积分区域可以是任何n 维多边形。 单纯形在0 ，l ，2 ，3 ，l 维空间有最简单的形状( 图3 3 ) ，与普通积分不同，单纯形 积分只以单纯形作为积分域。单纯形有正或负的方向，正负方向分别规定正负体积。 p op 0p 1p o p ip o p i 二维单纯形p o p l p 2 是一个有向的三角形，它的体积是 j ：躬羔i s ， i i 屯y 2 l 单纯形p | p o p 2 的体积是单纯形p o p l p 2 的负体积。 二维单纯形积分定义为 第三章数值流形方法的基本方程及分析过程 l 躺m ，y ) d ( ) 一i g n ( j ) j p o bm ，y ) d x d y ( 3 - 5 4 ) p 5 p 4 如图3 4 所示，已知多边形p i p 2 p 3 p 4 p 5 ，对任意点p o ，二维单纯形体积( 面积) p o p p 2 ， p o p 2 p 3 ，p o p 3 p 4 ，p o p 4 p 5 及p o p 5 p l 的代数和是多边形p i p 2 p 3 p 4 p 5 的面积4 ，设p 6 _ p l ，则 4 = 三喜i 主。妻= 圭喜窿。妻，i c ，p s s ， 对一般被积函数，( z ，j ，) ，多边形p l p 2 p 3 p a p 5 的标准积分是单纯形积分的代数和，表 l 蝴b m y 以力= 善k 肥y ) d ( 捌 ( 3 砌5 6 3 5 计算程序及算例 3 5 1 程序的完善和改进 本文以石根华早期开发的d o s 环境下的数值流形程序为基础，采用面向对象程序 设计( o o p ：o b j e c t e d o r i e n t e d - p r o g r a m m i n g ) 方法对其进行完善和改进，并运用v i s u a l c 十+ 语言将其改写成为w i n d o w s ) 口环境下的可视化数值流形方法计算程序。 ( 1 ) 对约束条件处理的完善 原数值流形方法计算程序只能处理双向约束的固定边界，不能计算半无限体自重作 用下的应力场。而在实际工程问题中，常有大量受约束的点( 面或节点) 需处理，且约 束的方向、形式各不相同。为此，本文对原固定边界的约束方法进行了改进。 改进后的数值流形方法计算程序可以处理三种约束情况：双向约束；z 向约束，y 向自由；y 向约束，工向自由。 ( 2 ) 计算过程的可视化 、 、 、 ， 、 ，# h p n 塑查兰堡主兰焦丝苎 原数值流形方法计算程序是在d o s 平台上用c 语言开发的，不仅应用不方便，且 不适用于目前广泛应用的w i n d o w s 平台。为此，本文采用面向对象程序设计方法改写 了原数值流形方法计算程序，实现了计算过程的可视化。 本文所完成的程序包含数据文件的存取功能、前处理功能、网格自动划分功能、数 值计算功能、数据文件输出功能。程序界面如图3 - 5 所示。 ( a )前处理 ( b ) 计算过程可视化 图3 5 数值流形方法计算程序的可视化界面 3 s 2 程序流程图 数值流形方法核心程序的计算过程及程序流程图如图3 - 6 所示。 第三章数值流形方法的基本方程及分析过程 图3 - 6 数值流形方法程序流程图 河海大学硕士学位论文 3 5 3 算例 ( 1 ) 岩体边坡 如图3 7 所示为某水库一岩体边坡的概化模型，考虑空库和岩体自重。截取边界约 束情况为：底边固定，左、右边均为x 向约束，y 向自由。n m m 计算模型如图3 - 8 所 示，模型的材料参数如表3 - 1 所示。 律； 图3 7 某水库岩体边坡概化模型( 单位：m ) 图3 - - 8n m m 计算模型 表3 - 1 模型材料参数 该岩体边坡的应力场如图3 9 所示，由图可见，该岩体边坡应力场符合一般力学规 第三章数值流形方法的基本方程及分析过程 ( a ) 水平向应力分布 ( b ) 竖向应力分布 ( c ) 剪应力分布 图3 - 9 自重作用下岩体边坡的应力场( 单位：p a ) ( 2 ) 重力坝 如图3 - 1 0 所示为某1 6 0 m 高混凝土重力坝概化模型，考虑坝体上游的水压力和坝 体自重。截取边界约束情况为：底边固定，左、右边均为石向约束，y 向自由。n m m 计算模型如图3 - 1 1 所示，模型的材料参数如表3 - 2 所示。 河海大学硕士学位论文 图3 1 0 混凝土坝概化模型( 单位：m ) 图3 1 1 数值流形方法计算模型 表3 - 2 基本材料参数 材料参数坝体坝基 密度p ( k g m 3 ) 弹性模量e g p a 泊松比y 2 4 0 0 2 5 o 1 6 7 第三章数值流形方法的基本方程及分析过程 n m m 和f e m 两种方法计算所得坝基和坝体应力场如图3 1 2 、图3 一1 3 所示，两种 方法坝体的竖向应力分布如表3 - 3 所示。由表3 3 可见，数值流形方法的计算成果是合 理可信的。 表3 - 3 竖向应力分布特征比较( 单位：m p a ) 方法应力性质应力范围应力集中部位及大小 拉应力主要分布于坝踵尖角处，但是范围不 拉应力 o 0 7 6 大，最大拉应力为0 7 6 m p a 。 n m m 在坝体和基础的大部分区域里出现压应力， 压应力 一5 2 o 在坝趾处出现很小的应力集中区，压应力最 高可达5 2 m p a 。 坝踵尖角处存在拉应力区，但范围很小，并 拉应力 0 1 1 2 且很快过度到压应力区。 f e m 在坝体和基础的太部分区域里出现压应力， 压应力5 8 加 坝址处的压应力达5 8 m p a 以上。 由数值流形方法与材料力学法计算所得竖向应力分布规律比较可知，地基变形对坝 体应力分布的影响只限于坝体下部，而坝体的中、上部分基本上符合直线变化的假定。 另外，数值流形方法中网格剖分大小不同对应力分布也有影响。总体上说，网格剖 分不同对应力的影响非常小，影响基本可以忽略不计。而在坝踵和坝趾较小的范围内， 应力分布随网格剖分的不同有较大的差异。究其产生的原因还有待进一步研究。 河海大学硕士学位论文 ( a )水平向应力分布 ( b )竖向应力分布 ( c ) 剪应力分布 图3 1