（机械设计及理论专业论文）刚塑性数值流形方法研究及应用.pdf

中文摘要 摘要 为了研究超出线性弹性变形达到屈服点应力后的应力及应变状态， 特别是在实际加工当中经常遇到的金属塑性成形这种大变形问题，先后 出现了弹塑性有限单元法以及刚塑性有限单元法。然而，在处理大变形 问题时，弹塑性有限单元法需要非常小的时间增量，计算费用非常昂贵 而且累积的误差很大；在可以忽略弹性变形量大变形场合，刚塑性有限 单元法虽然解决了在每一时间增量进行对应力的直接求解从而消除了应 力累积误差的问题，但是其计算精度的提高还是需要极大地依赖于大量 的计算网格数据等方面。 现代的数值流形方法是一种具有一般意义的数值方法，它与传统的 有限单元法单一的覆盖不同，它采用了两套可分离的数学网格和物理网 格，其数学网格定义了求解的精度，从而有效地避开了为了提高计算精 度需要增加大量节点数据的问题。而且，数值流形方法也可以有效的解 决非连续变形等方面问题。该种方法将有限单元法和非连续变形分析统 一了起来。由于数值流形方法的诸多优点，因而受到了越来越广泛的关 注。 本论文就是主要针对传统的刚塑性有限单元法的优缺点，并在此基 础上把数值流形理论方法引入到刚塑性问题的分析当中，从变分原理出 发，逐步建立起刚塑性数值流形方法的理论体系及相关计算公式。考虑 到四边形等参单元的广泛的应用性，因而本论文着重推导了等参四边形 的刚塑性数值流形计算表达式。 关键词；刚塑性，数值流形，等参四边形 a b s n a c t a b s t r a c t i no r d e rt o s t u d y t h es t a t e so fs t r a i na n ds t r e s sw h i c he x c e s st h e b e n d i n gp o i n ta f t e rt h e1 i n e a re l a s t i c i t y ，s p e c i a l l yt h eb i gd i s t o r t i o nw h i c h b eo f t e nm e td u r i n gt h em a n u f a c t u r e ，t h ee l a s t i ca n dp l a s t i cf i n i t ee l e m e n t m e t h o da n dt h er i g i d p l a s t i cf i n i t ee l e m e n tm e t h o do c c u r ss u c c e s s i v e l y b u t ， w h e nd e a l i n gw i t ht h eb i gd i s t o r t i o n ，t h ee l a s t i ca n dp l a s t i cf i n i t ee l e m e n t m e t h o dn e e d st h ev e r ys m a l li n c r e m e n to ft h et i m e ，w a s t e sm u c ht i m eo f c a l c u l a t i n g a n dc p u ，a tt h es a m et i m e ，t h ec a l c u l a t i n gt o l e r a n c ew i l lb e v e r ym u c h ；a tt h eo c c a s i o nw h e r et h ee l a s t i cd i s t o n i o nc a nb en e 9 1 e c t e d ， a l t h o u g ht h er i g i d p l a s t i cf i n i t ee l e m e n tm e t h o dc a nr e s o l v et h ep r o b l e mo f t h ec a l c u l a t i n gt o l e r a n c eb yc o m p u t i n gt h es t r e s sd i r e c t l yd u r i n ge v e r yt i m e s t e p ，i th a st od e p e n do nt h ea b u n d a td a t ao fg r i d st oe n h a n c et h ep r e c i s i o n o fc a l c u l a t i o n t h em o d e r nn u m e r i c a lm a n i f 0 1 dm e t h o d ( n m m ) i san e wm e t h o do f n u m e r i c a iv a l u ew i t hg e n e r a lm e a n i n g w h a ti sd i f 南r e n tf r o mt h et r a d i t o n a l f i n i t ee l e m e n tm e t h o dw h i c ho n l yh a st h es i n g l ec o v e r i n gi st h a tt h en m m a d o p tt w os e t so fs e p a r a t eg r i d so fm a t h sg r i da n dp h y s i c sg r i d ，t h e r e i n t o ， t h em a t h sg r i dd e f i n e st h e p r e c i s i o no fc a i c u l a t i o n ，w h i c hc a na v o i dt h e p r o b l e mt h a tn e e d st oa d dt h ed a t ao fg r i df o re n h a n c i n gt h ep r e c i s i o n f u r t h e r m o r e ， n m ma l s oc a nr e s o l v et h e p r o b l e m o fd i s c o n t i n u o u s d i s t o r t i o n n m mu n i f i e st h en n i t ee l e m e n tm e t h o da n dt h ed i s c o n t i n u o u s d i s t o r t i o na n a l y s i s d u et ot h em a n ym e r i t so fn m m ，n m mh a sb e e n p a i d m o r ea n dm o r ea t t e n t i o n i nt h e d i s s e r t a t i o n ，b a s i n g o nt h em e r i t sa n dt h e d e f 色c t so ft h e t r a d “i o n a lr i g i d - p l a s t i cf i n i t ee l e m e n tm e t h o d ，t h ea c a d e m i cs y s t e ma dt h e c a l c u l a t i n gf o r m u l a sh a v eb e e ne s t a b l i s h e db yi m p o r t i n gt h en m m t ot h e r i g i d p l a s t i c f i n i t ee l e m e n tm e t h o d c o n s i d e r i n gt h ee x t e n s i v ea p p l i c a t i o n o ft h ei s o p a r a m e t r i cq u a d r a n g l e ，t h ef o r m u l a so f r i g i d p i a s t i cf i n i t ee l e m e n t m e t h o dw i t ht h ei s o p a r a m e t r i cq u a d r a n g l ea r ee m p h a s i z e dt od e d u c e k e yw o r d s ：r i g i d p l a s t i c ，n m m ，i s o p a r a m e t r i cq u a d r a n g l e 2t t 第一章绪论 第一章绪论 1 1 学术背景 1 1 1 刚塑性成形的发展及背景 金属塑性成形过程是一个很复杂的过程，材料特性、变形速度、温 度、摩擦条件、坯料形状及尺寸、模具形状对成形过程都有影响，而这 些因素及作用是我们研究的主要对象。 以计算固体应力、应变为目的的有限元法是从小变形理论的弹性有 限元发开始的。6 0 年代，为了研究超出弹性变形达到屈服点的应力和应 变状态，小变形弹塑性有限元法得到了发展。然而发现，在将速率形 式的本构方程向前积分时，需要非常小的时间增量，如果采用这种方法 跟踪大变形，计算费用则为非常昂贵，误差的累积也很大，还有得不到 正确解的可能。7 0 年代，基于有限变形理论的大变形弹塑性有限元法得 到了发展。与此同时，考虑到在很多的实际应用中，塑性变形量远大于 弹性变形量，也就是说弹性变形量在很多实际的问题上面是可以忽略不 计的，在这个基础之上出现了刚塑性的有限单元法。实际上，直到7 3 年，c h l e e 和小林史郎提出了刚塑性有限元法，金属塑性成形有 限元法的发展才取得了根本性转变，向前迈出了一大步。1 9 7 3 年，l e e 和k o b a y a s h i 以矩阵分析的名义导出刚塑性有限元中l a g r a n g e 法的矩 阵方程，剐塑性有限元从塑性成形工艺的数值模拟中显示了巨大的潜力， 成为塑性成形工艺的重要分析工具。k o b a y a s h i 及其合作者用刚塑性有限 元和刚粘塑性有限元分析了墩粗、挤压、扎制、锻造等体积成形以及胀 型、冲孔、方盒拉延等板料成形。1 9 7 5 年，z i e n k i e w i c z 提出了刚塑性有 限元中的罚函数法，并利用这些方法模拟了热态挤压、轧制及板料成形。 1 9 8 2 年m o r i 和o s a k a d a 提出了刚塑性有限元中的材料可压缩法对多孔 类可压缩材料的成形工艺进行了模拟分析。a l t a n 等利用大型刚塑性有限 元分析程序a l p i d ( a n a l v s io f p l a t i ei n c r e a m e n t a ld e f o r m a t i o n ) 对塑性 成形的问题进行了模拟分析。近年来，一些学者对三维成形问题进行了 刚塑性有限元分析。k i k u c h i ，s h e p h a r d 等对有限元分析中的网格自动主 成方法进行了研究，进一步提高了有限元分析的自动化。c h e n g ，y a n g ， 广东工业大学工学硕士学位论文 b a d a w y ，g e l t o n 等研究了网格再划分问题。k o b a y a s h i 等人还提出刚塑 性有限元的反向模拟技术，并用此技术对一些简单的成形问题进行了预 成形设计，这些工作都丰富和发展了塑性有限元的具体处理技术和应用 范围。另外，p i l l i n g e r ，c h e r o t ，s h i a u 和k o b a y a s h i 考虑了塑性成形时 工件的摩擦边界情况，使刚塑性有限元分析塑性加工间题时更符合实际 2 【5 1 。 刚塑性有限单元法不需要像弹塑性那样对应力增量进行求解，而是 在每一个时间增量都直接求出应力，因而就不会出现应力的累积误差方 面的问题，并且，由于它是一种流动型的列式，所以可以取较大的增量 步长，从而有效地减少计算量，在保证足够的精度下提高计算的效率。 需要注意的是，在刚塑性有限元分析中，它必须随时间向前携带的历史 变量仅仅是与材料结构变化有关的变量；它通常采用率方程，也即列式 本身是根据小应变增量建立的，这样，变形后的构件可以通过在离散空 问上对速度进行积分而获得，从而避开了几何非线性问题，这些特点都 使得刚塑性有限元列式比较简单，易于实现编程计算”。 但是，在另一方面，刚塑性有限单元法也有明显的局限性”，譬如： 由于忽略了弹性变形，这种方法仅仅适合于塑性变形区的分析，不能直 接分析弹性区的变形和应力状态；弹性变形区的解仅仅在整体意义上是 正确的，即维持平衡方程；它不能处理卸载问题和计算由此带来的残余 应力和残余应变。因而与大变形的弹塑性有限元法相比，载变形量很小 的场合，刚塑性有限元分析的精度较差。虽然这样，在工程实践当中， 并不总是需要考虑卸载、残余应变以及残余应力，所以，刚塑性有限元 法仍不失为一种有效的分析手段。特别是在金属体积成形方面，刚塑性 有限单元法是主要的数值模拟方法。 但是，不管怎么样，受制于有限元根本思想的束缚，在有限单元法 里面，无论是将位移作为未知量还是将速度作为未知量，这些未知量的 构造基本都只是一场常函数的模式，因而要提高计算精度的话就需要增 加节点数量，但是这样就会引起计算量大大增大等问题。现代的数值流 形方法能有效的解决这个问题，它不再简单地将这些根本未知量作为常 函数的模式来进行处理，而是将数学流形的概念引入到位移场或是速度 2 第一章绪论 场的构造当中，也即将多项式的概念引入到位移场或是速度场构造当中， 从而有效的解决节点数目与计算量、计算精度的矛盾等方面的问题，在 不增加节点数目的情况下就可提高计算精度。另一方面，由于现代数值 流形方法将传统的有限元、块体计算等问题统起来，在理论研究或工 程的实际应用当中发挥着越来越大的作用。 1 1 2 数值流形的数值方法 数值流形方法是石根华创立的一种更新、层次也更高的现代数值方 法，它是以数值流形为核心，在非连续变形分析( d d a ) 的块体系统非 连续运动学理论的基础上融入了有限元和解析法的连续分析方法。从而 形成了一种至少包括有限元、d d a 、解析法在内的一种新的统一计算形 式，有限元以及d d a 都是数值流形方法中的两个特例【。 在物理上，材料对象通常会具有不同的形状，也会有不同的解析区 域：断裂的、连续的等。特别是在对分析由多种材料构成的不连续的分 析区域的时候，通常的数值分析方法会遇到很大的困难，对具有这些特 征的材料进行分析的时候其处理的过程也非常的复杂和麻烦，因为常规 分析的近似方法只在代表一小部分的局部连续区域内是可行和有用的。 在数值流形方法当中，其“流形”的概念是把许多个别的重叠的局部区 域连接在一起去覆盖全部的材料体，因此，总体形状函数可用局部覆盖 所定义的函数来计算。总的来说，数值流形方法的基本的思想就是在求 解区域上构造一组覆盖函数，使得这些覆盖函数具有如下的两个基本特 征：1 ) 局部非零性：覆盖函数只在局部区域内不为零，也即在局部区域 内有解：2 ) 在求解的区域内，覆盖函数之和为l ，也即局部区域内的总 体位移函数由所有的覆盖函数所组成 1 。 在数值流形方法中，其与数学流形的区别在于：数学流形在整个求 解区域内总体位移函数是处处连续可光滑高度微分的，它可完全定义与 覆盖无关：而数值流形的总体位移函数只在覆盖的基础上来进行定义， 而且是可分片微分，在覆盖的接触面上是完全非连续的。传统的数值计 算方法通常只是只有一种网格，而数值流形方法则需要定义两套分开而 且完全独立的网格：数学网格以及物理网格。其中，数学网格只是定义 广东工业大学工学硕士学位论文 了近似解的精度，可由用户来进行选择，它由占整个材料体的许多的有 限重叠覆盖所组成，譬如规则的网格、有限元网格或是级数的收敛域都 能转换为有限数学网格；物理网格则代表了材料条件，不能由用户人为 的选择。作为实际的材料边界，物理网格定义了积分区域。实质上，物 理网格包括了材料体的边界、裂缝、块体以及不同材料区域的交接面， 譬如不变化的水面也是物理网格的一部分。有了数学网格和物理网格， 就可以形成物理覆盖系统了。所谓物理覆盖，就是裂缝或是块体边界等 把一个数学覆盖分成两个或更多的完全不连续的区域，这些区域就定义 为物理覆盖。因此，物理覆盖是不连续缝对数学网格的再次剖分 ”。 数值流形方法有可能为偏微分方程的求解创造出更高的效率和更好 的算法。例如经典的解析法虽然有很高的收敛速度和编程的优点，特别 是在处理无穷域和奇点问题上常常优于有限元法，但是，由于不能常常 满足求解区域上的全部条件从而限制了它的应用发展。采用数值流形， 则可在不同的区域采用不同的级数进行展开，这一优点就很好的解决了 上述的问题，而且还可在不同的连续区域把解析法和有限元法结合起来 进行协调计算。另一方面，鉴于目前非常成熟的有限单元法的计算程序 【8 “”，结合数值流形的计算特点，也容易实现刚塑性数值流形方法计算 程序的编制。 数值流形方法可根据实际需要对物理覆盖系统上的每一个物理覆盖 定义各自独立的覆盖位移函数，然后在几个覆盖的公共区域那将其所有 覆盖上的独立覆盖加权求和就能形成适用于该域的总体位移函数。因此， 数值流形方法的求解具有适应性，传统的有限元方法只是数值流形方法 的一个特例，即数学覆盖与物理覆盖完全重合时，数值流形方法便变为 有限元方法。 相对于有限元，数值流形有以下三个主要的差别： l 、数值流形方法中覆盖位移函数可以是任意级数形式，而不是形函 数。 2 、在数值流形中，接触理论基本上采用不连续变形分析方法的理论。 3 、数值流形方法提供了惯性力、速度等的计算公式，在静力计算过 程中，若有逐渐转化为时间变化的运动趁势或运动状态，仍可继续分析 4 第一章绪论 和模拟。 基于上述的数值流形方法相对于有限元方法的优点，数值流形方法 是一种比有限元更为普遍的方法。在解决不连续变形( d d a ) 或大位移 方面得到了突出的应用。 目前，数值流形方法的理论研究相对比较集中于弹塑性材料的领域， 而对材料成形具有非常意义的刚塑性材料的数值流形的理论研究还是比 较薄弱。本论文所研究的主要课题就是将刚塑性理论引入数值流形方法 中，建立起基于刚塑性理论的数值流形理论体系和数值分析方法。该课 题的研究对制造与成形过程的数值模拟具有较重要的理论意义和应用价 值。 总之，作为一种新的数值计算方法，数值流形方法具有现代的有限 单元法和经典的解析法所没有的优点，它可以把有限元、d d a 和解析法 融合到统一的计算形式当中，所以它可以同时把连续和非连续、数值法 和解析法放在一起来进行计算，是当前最有发展前景的一种前沿的数值 方法之一。 1 2 国内外文献综述 数值流形方法是石根华9 0 年代初继创立块体理论和d d a f 非连续变 形分析) 之后，并在此基础上首创的一种更新的、层次更高的现代数值方 法。流形方法以拓扑流形和微分流形为基础，它有分开的且独立的两套 网格：数学覆盖和物理网格。数学覆盖定义近似解的精度，由用户选择 且独立于物理网格划分，可以是任何规则或不规则的格子，物理网格则 包括材料的边界、裂缝、块体和不同材料区域的交接面，它不能人为选 择。物理网格对数学覆盖再剖分就形成了覆盖材料全域的流形方法求解 的物理覆盖系统。对物理覆盖系统上的每一个物理覆盖可以定义各自独 立的覆盖位移函数，然后在几个覆盖的公共区域( 即流形单元) 内将其 所有覆盖上的独立覆盖位移函数加权求和就能形成适应于该域的总体位 移函数。 借助于建立的物理覆盖系统，流形方法能够很方便地分析有裂缝的 结构或块体以及有很大柔性的和明显可见的大变形或大位移。而且流形 广东工业大学工学硕士学位论文 方法覆盖位移函数的采用也能够为偏微分方程的求解创造出更高的效 率，提供更好的算法。例如经典的解析法有很高的收敛速度和编程简单 等优点，尤其对处理无穷区域问题和奇点问题，常优于有限单元法。但 是它不能满足求解域上的全部条件，应用受到限制。采用流形方法的覆 盖函数，则可在分析结构的不同区域采用不同的级数展开式，或在求解 域的不同地方采用不同的方法，如有的地方采用解析法，有的地方采用 有限元法等，而且在不同区域的联接处，整体求解的位移函数是自然协 调的，软件实现时与通常的协调方法没有多大的区别。这样，流形方法 很容易地就解决了人们多年来一直研究而至今尚未有好的解决办法的局 部区域解析法以及有限元与解析法相结合的办法，同时对恢复经典解析 方法的应用和研究也将起到积极的推动作用【l ”。 流形方法是一种高度统一的数值方法，它用连续的和非连续的有限 覆盖系统把连续的和非连续的变形问题分析融为一体，统一解决了有限 元、非连续变形分析( d d a ) 和解析法的计算问题。国内外的不少学者 在此方面也展开了积极的工作。1 9 9 6 年，已有学者将流形方法应用于工 程实践及模拟裂纹扩展中的实现过程，目前在岩土工程分析中已经得到 了广泛的应用。c h e ng 和o h n i s h i 等用完备的高阶覆盖函数来提高流形 方法的求解精度，并给出了其具体实现过程【2 0 】。1 9 9 9 年，有研究者将增 量流形方法推广到岩石大变形问题，并模拟了岩石中裂缝的发展和大变 形2 1 【22 1 。2 0 0 0 年，还有学者应用流形方法思想，通过引入广义节点的 概念，对传统有限元法进行改进，提出了广义节点有限元法 2 3 2 “。还有 学者对裂隙扩展、多裂隙结构体的破坏及跟踪、移动中悬梁的自重弯曲 破坏以及格流形方法与局部子域b e m 法相结合进行裂隙扩展的模拟等 方面做了大量研究 2 7 】【28 1 另外，还有一些研究者介绍和探讨了数值流形 方法及其变分原理及应用前景 2 9 【3 ，在数学上给出了严格的证明推导， 建立起数值流形方法的数学基础在流形方法中考虑侧向影响，以改进 平面模型的模拟”】，以及将数值流形方法用于求解弹性地基板问题等方 面均有研究文献报道【”】近期的研究报道主要有对物理覆盖系统的构 成、自动剖分、网格重构技术及程序设计等方面【3 3 3 7 】；在对总体位移函 数的构成方面也有更具体的进一步研究同时，还有利用流形概念对传 第一章绪论 统有限元方法进行改进、建立可具有任意高阶多项式插值函数的广义自 由度有限元方法王芝银利用数值流形方法中的惯性力概念对传统流变 分析的有限元格式或是数值流形方法进行了改进” 3 ”也有相关学者通 过使用高阶覆盖位移函数来提高数值流形方法的求解精度4 0 “2 1 ，部分学 者将数值流形理论引入到自适应分析中进行研究及应用 4 ”。这些研究成 果的出现完善和发展了流形方法的内容。但是，也可以看出，将数值流形 方法与刚塑性力学结合起来的研究还是没有引起足够多研究人员的重 视，虽然已有学者将数值流形方法导入弹塑性分析中 4 ”，甚至是充分地 将数值流形方法引入到分析不连续问题、裂纹扩展等复杂问题上去【45 1 ， 但是弹塑性数值流形方法还是比较多的应用于岩土力学分析中。对刚塑 性数值流形方法的研究还是处于开始阶段。由于流形方法是一种新的数 值方法，它的发展才刚刚起步，研究对象还比较狭窄，目前大多数是针 对岩土工程的块体或不连续的裂缝数值模拟，应该说，流形方法还应该 有更大的发展空间和更多的工作要做。 1 3 研究的主要内容 本论文针对目前国内外的数值流形方法的研究动向，将数值流形理 论引入到刚塑性材料的理论研究当中，主要是通过刚塑性材料的基本力 学方程，结合数值流形理论的一般思想，从而建立起刚塑性材料的数值 流形方法的求解体系，并较为详细地推导了刚塑性数值流形方法的相关 计算公式。采用c c + + 语言编制了相关的计算程序，最后利用本论文的 刚塑性数值流形方法对两个实例进行了计算分析。计算的结果表明该方 法的有效性和正确性。 广东工业大学工学硕士学位论文 第二章数值流形方法的基本理论 2 1 基本理论 2 1 1 有限覆盖的数学覆盖和物理覆盖 数值流形方法是把有限覆盖系统连接在一起，用以覆盖全部的材料体，材 料体的总体形状可用局部覆盖所定义的函数来计算。这种计算方法使用了两套相 互独立的数学网格和物理网格。数学网格定义数值解的精度，由用户自行选择； 物理网格是由材料的物理及几何性质所决定，用户不能随意选择。在数值流形方 法中，这两种网格都需要进行单位分解并转换成覆盖。数学覆盖可以是任意形状 的，而物理覆盖则完全是由材料的物理及几何性质所确定。一般地，物理覆盖具 有如下的规定1 ，3 6 ： ( 1 ) 覆盖的区域是包含在数学覆盖中的材料，或用数学语言来说是数学覆 盖和材料场的交集： ( 2 ) 如果材料体的边界、块体边界或是裂缝把数学覆盖划分成完全的隔离 区域，每一区就是一个物理覆盖。因此，物理覆盖是数学覆盖的再剖分。 图2 1 有一条裂缝的一般覆盖 f i g 2 1t h e g e n e r a l c o v e r sw i t ho n ec r a c k 如图2 1 所示，三角形表示的是材料体的形状，中间的一条粗黑线条表示 的是这个材料体的裂缝。这个材料体由一个圆表示的数学覆盖形和两个矩形表示 的数学覆盖k 及吒所覆盖。图2 1 中，数学覆盖k 被材料裂缝分成两个物理覆 盖2 ，及2 ：，数学覆盖巧则被分成3 ，及3 ：两个物理覆盖，而数学覆盖k 并没有被 第二章数值流形方法的基本理论 裂缝进一步划分，实质上1 ，就是1 。而图2 1 中的物理覆盖1 ，2 2 3 2 则表示该小部分 是由3 个数学覆盖所覆盖起来的。 数值流形方法将这两种覆盖重叠在一起，并用物理覆盖把数学覆盖进行重 新的剖分裁剪，把数学覆盖叠入物理覆盖，形成新的连续和非连续的供计算的有 限覆盖系统。在这个系统中，物理覆盖代替单元的节点，覆盖的并集交线代替单 元的边界，覆盖重叠的交集代替单元，这样，原始单元就有可能被块体的边界、 裂缝等划分为几个流形意义上的单元。 其实，数值流形方法能够进行连续材料的有限元法的计算。在将有限元网 格转换为数值流形方法的有限覆盖之后，裂缝可作为新的物理网格输入，这样同 一个有限元网格就可计算同一材料体中的裂缝。有限元网格可用来定义数值流形 方法的有限覆盖。如考虑任一节点，含有这一节点的所有单元形成一个数学覆盖。 2 1 2 有限覆盖系统的覆盖函数以及权函数 对数值流形方法，覆盖位移函数是对各个物理覆盖独立定义的。局部位移 函数可连接起来形成整个材料体上的总体位移函数。总体位移函数是一般性的， 它足以适合位于移动边界内的连续或不连续的材料的广泛变化。位移函数与材料 边界无关，如果材料只是占有单元的一部分，位移函数仍是相同的。 一般地，二维物理覆盖u 上的覆盖函数“，( x ，y ) 可以表示成如下形式： ，( 工，y )工，_ y v ( 2 1 ) 这个函数可以是常量、线性、高阶多项式或是局域级数。这些覆盖函数用下面的 权函数w ，( 工，力连续起来： ”r ( w ) 0 训5 v( 2 2 ) 【( x ，力= ox ，y 萑 并且 ，( x ，y ) = 1 j ，y e u ， 权函数w 。( j ，y ) 的含义就是加权平均 函数“。“y ) 的百分数。 ( 2 3 ) 它对所有含有( x ，力的物理覆盖取每个覆盖 用权函数w ，( x ，y ) 就可定义整个物理覆盖系统的总体函数。如下 广东j = 业大学j _ 学硕士学位论文 ”“力= 宝w “力m “力 r l ( 2 4 ) 上式中，”的含义是表示组成整个物理覆盖系统的物理覆盖的个数。实质上， 当数值流形方法退化为有限单元法的时候，这里的权函数就相当于有限元法里面 的形函数。这一点在后续的论述中会得到体现。 2 1 3 数值流形单元 数值流形方法中，单元定义为物理覆盖重叠的部分。 如果是对二维数值流形方法的计算，定义在物理覆盖己，上的覆盖位移函数 “，) 和p ，) ，则： 峨，一 叼! ， ( 25 ) v 。( x ，y )( # ， u ， 利用权函数嵋y ) ，整个物理覆盖系统上的总体位移函数可根据覆盖位移函数 确定。在整个材料体上的总体位移函数“，( t y ) 和v ( x ，y ) 足两个总体函数： 隧辫喜吣一，能瓣；| ；善k 阮 】l d f ) = 喜删旺s ， 其中： k c 圳= 乏：譬暑：譬暑 慨 = 缸。一。一。f k “力j = 啊，) 协) = d 。口。叱r 以上各式中下标f 表示第j 个物理覆盖。 物理覆盖的每一个交集定义为单元，并且： ：! 。，y ! 。18 5 兰 ( 2 7 ) z ( x ，) = 0p 而单元是覆盖玑t 。) ，以( 2 j ，玑( 的交集，则： 糕嵩老： 泣s ， 单元p 的整体位移函数就是 第二章数值流形方法的基本理论 熙辨轨心朋溉篙斟 = 喜善 w “n x ，y 訾p x y 。：。，。，y 是。，。，y ， 2 ：i ：j 1 ) q r，、 = k ，亿力，) = 窆阮( 训) 缸) ( x ，力e 。 ( 2 9 ) 上述各式中， d ) 也称为物理覆盖的自由度。 2 1 4 一般覆盖的位移函数模式及整体位移模式 在覆盖u ，上，( “，( x ，y ) ，v 。( x ，y ) ) 是在z 和y 任何方向上的点( x ，y ) 的位移。覆 盖位移函数通常可取如下形式： 覆盖u ，上的常量函数，此时m = l ： 髅瓣鼢 眨 髅瓣：础i i t ， 巨 ( 2 1 1 ) 明显的可以看出来，如果采用常量函数的表达方式，在二维有限元的情况下， 形函数正好是权函数。对四边形的流形单元( 下面会对其进行论述) 而言，当采 用胛= 3 的完全一阶近似时，可以知道物理覆盖的自由度数将会是6 ，而传统的 有限单元法对二维问题而言其各节点的自由度数却是2 ，这一点与有限单元法的 区别在以后各小节的论述中将会得到充分的体现，在此可以粗略的看出，数值流 形方法在求解的精度方面将会大大的优于传统的有限单元法，因为其“节点”的 位移函数构造形式在数学上就是可以多阶的，而非简单的常量表达。 如果采用的是等参四边形流形单元，在覆盖玑上其位移函数采用完全一阶 近似的形式的情况下，即= 3 时，将引入( s ，卵) 为自变量，是因为在这里引入的 东工业大学工学硕士学位论文 是有限元等参单元里的母单元表达方式，也即将自变量( x ，y ) 经过等参变换为 ( s 川) ，这样就可以有效的解决复杂网格划分的问题，就可以采用有限元成熟的 积分方案，为以后的积分运算带来方便。 单元e 是g = 4 个覆盖己，u 形成的交集四边形流形单元 是四个物理覆盖的交集。如果使用形函数作为流形方法中的权函数，即： ( s ，) = m ( s ，叩) 则流形单元e 的整体位移函数可采用如下形式： 书譬券j = 霎c 岛叩叶i 譬舅) = 喜m c 功氍芝宝) ：至 l _ l f 如 i 一： = 雾， ：硎窑 l 如 【氐 n!o n 。n # ， o 0 ，。 o o ，。 o 4 n ? n j 。 h = 1 = 喜甜紫m z 捌划= 莩雾k 功舰 智智l o m 兀( 五力| i 一副f 鲁鲁”“ = 阮( s ，叩) 】 d j 其中： ，。( x ，力= 1 ，：( t y ) = x ，3 ( 五y ) = y 揍 ( 2 1 2 ) 磊如丸如如以 ，ll刊iiill【 o 嘎o o m x毗o o m m o g 己目 第二章数值流形方法的基本理论 k c s ，叩， = 1 毪。，y m z 。，y ， ， k c 功】= 亿z z ) = ：芸2 乏是： 对四边形流形单元而言，也即： d ，) 剖 ，( 功 乞， 叩) 。幢咖 f 2 4 ( 们 ，幢功 f 25 幢功 f 。( 功 f ze 慨功j ) 】- 亿。l ：叫m m 墨，功气。z 。掣z 眨 由上述可阻看出，流形单元位移函数模式与有限元位移函数模式的差别： 1 ) ：当在流形单元的位移函数模式中令= 1 ( 组成单元的某个覆盖配上的 位移模式不再为多项式的形式而变为常函数的模式) 即变为有限元位移函数模 式。 2 ) ：也正是在数值流形中任一覆盖以采用了多项式的形式取代有限元的常量 模式，使得计算的精度得到了提高，也即可以采用较少的结点数也即单元数就可 以获取较大的计算精度了。 3 ) ：也正因为流形单元的任一组成覆盖采取了多项式形式，因而其未知量也 增加了( 由两个未知量( “。h ) 变为了六个未知量，在流形单元计算中，由其单元 的位移函数模式可以看出：这六个未知量实际上就作为某个覆盖的多项式列式中 的系数日l 2 f 砖( ，= 1 ，2 ，3 ) ) 。 4 ) ：从以后的推导过程中将可以发现，流形单元的各个矩阵表达式之所以与 有限元有差别都是由于流形单元的覆盖系统采用了多项式列式而引起的，从这点 可以看出，有限元其实就是数值流形的一种特例，也即埘= l 情况下的特例。 需要注意的是，如文献【4 4 】所论述的那样，“在数值流形方法中，常用的覆 盖函数并非总是最佳选择。随着刚度矩阵的形成，选用常用的覆盖函数时，在非 对角元上出现绝对值很大的元素”，并且这“会增加刚度矩阵的条件数，尤其是 在用刚性弹簧约束位移的情况下”。建议在必要的时候“采用局部化较好的覆盖 函数，取代常用的关联于全局坐标的覆盖函数”。或是可以采用文献【3 7 】的方法 来进行修正：也即“将覆盖位移函数对应的求解物理覆盖自由度全部置为零，具 体编程时用乘大数法引入此边界条件”。 2 1 5 一般覆盖的平衡方程式的系数矩阵 如果物理覆盖的个数为n ，每个物理覆盖有2 坍个未知数，即： 广东工业大学工学硕士学位论文 d ，= p 。d ，：d ，2 。一lt 2 。 7f = 1 ，2 ，”( 2 1 4 ) 总势能具有如下的形式： 肚扣叫怪 f e ，l e + 研谚蟛珥e l = 圭 d ) 医】 d ) + d r p 卜c + c 丘1 2置i3 k 2 2足2 3 墨：氏 足。2 岸n ( 2 1 5 ) 由于每个物理覆盖有2 个未知数，因而每个覆盖都有2 m 个自由度。由 d r 是一个l ( 2 m x 而矩阵， d 是一个( 2 m n ) 1 的矩阵，由此可以知道上式中的 子矩阵k 。是一个2 聊2 的矩阵。根据总势能的表达式： k 1k 。 7( 2 1 6 ) 利用最小总势能推导，下式中的第f 行由2 m 个线性方程组成： 要 0 ，州厶2 m ( 2 - 1 7 ) 出t ， 子矩阵心的第r 行第j 列元素是： 矗n 刚幺。晰 眩 e 的第，行元素是e 为零时也即扫 为零时的导数： a e f 0 1 面_ b ”。唧， 扫1 ，2 ，2 “ ( 2 1 9 ) 式中，d ，是覆盖f 的位移。而 罢0 ，罢：o ( 2 2 0 ) 翻卯 则分别表示作用在覆盖上的所有的载荷或是接触力在x 和y 方向上的平衡 式。 q b b 珥 n 月 h k b k k 第二章数值流形方法的基本理论 平衡方程式如下： 薰篓嵯 足2 lk 2 2足2 3一k 2 。| | d 2 k 3 1足3 2”k 3 。l d 3 k 。lk 。2五一k 。j l d 。 曩 只 _ e c ( 2 2 1 ) 上式中，f 是覆盖f 上的载荷，它被分配到2 m 个位移变量上，子矩阵丘“与 覆盖f 的材料特征和k 。有关，此处f j 是由覆盖f 和覆盖，之间重叠和接触来规 定的。 2 1 6 流形单元的刚度矩阵 对一般覆盖的流形方法而言，流形单元是物理覆盖相交的公共区域。刚度 矩阵的积分域就是整个流形单元，正如上述所说，积分域实质上就是数个物理覆 盖u ，虬的交集。 流形单元整体位移函数： 骶辨鼽砖 蒜掰 在数值流形方法中，应力与应变的关系和有限单元法相同。 针 抛( x ，y ) 玉 加( z ，y ) 砂 丝坚盟塑! 塑 砂 瓠 又对某一物理覆盖u ： 防( 五y ) 】_ k 。z ：一乙j l ，( x ，力：( x ，y ) ，( x ，力 lf 2 1 ( x ，y )f 2 2 ( 工，y )屯，3 ( 】力 所以： f 1 。g 一1 ( x ，y ) 乞2 ( x ，y ) h 翱( 囊y ) t 翔( x ，_ y ) j ( 2 2 2 ) 对覆盖矿， 陋，】- 则： 广东工业大学工学硕士学位论文 西，( z ，y ) 僦 ，( 置力 洲 堡坐：塑垦兰塑 砂 出 可表示为 h ， 勺 = 协州) 占啦) 川 西l3 积 西23 西l3a 23 却苏 = 陋。王d 。】 对弹性材料，由单元的弹性应力所做的应变能e 为 e 。= l l 。+ s 芦y + y 矿。出母 式中，s 为单元的材料面积。 所以： 耻f 扣叽+ 髟q u 蚴 = 导 d 。) 7 瞳】7 陋k 汹砂 = 吾娩) 1 小b 。r 陋p 。蚴i d 。) = 要 d 。y 0 陋。r 陋p 。她 式中，陋】是弹性矩阵。 1 6 西1 2 。 础 a r 22 。 口v 西l ，h研2 ，2 ， 却a o ( 2 2 3 ) ( 2 2 4 ) ( 2 2 5 ) ( 2 - 2 6 ) h 川 ， = 等等辩卜琶可维 肋肋一蹦 rl，k=_【 日 第二章数值流形方法的基本理论 因此，流形单元的刚度矩阵为： k 刚， js 瞳j r 陋拉】 ( 2 2 7 ) 2 1 7 流形单元一般覆盖上的点载荷矩阵 在数值流形方法中，点载荷可加载在单元内的任一点上。假设点载荷 忙，作用在单元内的某一点k ，) 上，力作用点上的位移是： 翟 = 麓 黧：羔；) c ：z s ， 那么，点载荷引起的势能是： e ，= 一k “( ，) + c v ( ，) j 蚶圳乏 眨z 。， 所以，载荷矩阵是： 阮( 枞) r ： l 1yj 五( 五( 五( l ( 并且： 也挑，】7 k ) = k j r 乏 = 苫苍 甜 o y n 。 删q 洲哥 = k x n j f f n t x f 。n j x f y n j y f 。n i y f y n 飞1 ( 2 3 0 ) ( 2 3 1 ) ( 2 3 2 ) 2 1 8 固定点矩阵 为引入位移边界条件，在数值流形方法当中，可假定固定点是在单元的 ( ，儿) 上用两个十分刚性的弹簧( 弹簧刚度为n 以及n ) 约束住，其位移应该 是： llhiniij 舱辨协 ， 假定的弹簧所产生的弹簧力是： ：叩一。粤3 l ( 2 3 4 ) 协f 2 1 啊v 翩，儿) f 。3 4 弹簧的应变能是： 影= 圭o ，( ，) 2 + 以v o ( ，蜘) 2 ) = 三埘洲7 陪冲( 枞胞 旺s s ， y ) 啦 珈圳= 舡( 粕，) 】k 。：，( ，虬) 】( ，蜘) b 爿 并且： 帕 7 l 号k 川 j 忉h 乩z 3 埘 眨s e ， 由z 的表达式，可得g 个物理覆盖的固定点矩阵k 町) j ： k k 6 ) - e ( i ) 。( 1 ) e ( 2 ) f ( 1 ) e g e ( 1 ) “2 ) e ( 2 ) k 。( 2 ) = 眯m ，】 1 8 ( 2 3 7 ) 以o 1。riooor brb，b， i( k k k 帅 帅 蝴 肠肠一岛 第二章数值流形方法的基本理论 oi 。1 曼岭詈 瞄爱等品警另 0l y mi 弘，描h p ：n 。n o p ；鼬q 。n o p 。y n 。n i o o p _ n 。n i o x p y n 。n io ) 节y n 。n i p 。x n i ni o p ；f n 。nr o p ；) 睁n 。n j o o p y x n 。n i o p y # n l nj o p y 弓剁t n i p ：y n 。n l o p 。秘n 。n i o p ；y 。n 。n | o o p y y n 。n i q p y ) 0 刑。n i o p y y | n i n i ( 2 - 3 8 ) ( 1 ) 当只有x 方向固定，y 方向的p ，= o ； ( 2 ) 当只有_ ) ，方向固定，z 方向的n = o ； ( 3 ) 如果x 和y 方向都固定，以o ，p ，o 。 实际上，在对固定点施加约束的时候，约束弹簧的刚度是非常大的，因而在 通常情况下我们可取以= p 。 2 1 9 初应力矩阵 在数值流形方法当中，其是按时序来进行计算的。前一步计算得到的应力将 作为下一步的初始应力，因此，初始应力对数值流形方法来说是基本的。 对流形单元，其初始常应力为： b 。 = 仃；r 吕) ( 2 3 9 ) 其势能是： 驴= j f + 盯一+ r 吕胁 = 见) 7 陋。r p 。扭方 m。川。埘。 矾。彬。 广东工业大学工学硕士学位论文 = s d 。) 7 陋。r p 。) 因此，载荷矩阵是： 一s 陋。r p 。) = 一s 并且记为： ( 2 4 0 ) ( 2 4 1 ) ( 2 4 2 ) 2 1 1 0 一股覆孟髓运发矩阵 在每一时间步中，都应使得位移足够小，从而使得最后的结果与具体的时 间步的大小无关。考虑当前的时间步时，流形单元上的任一点( x ，y ) 与时间相关 的位移是0 ( x ，弘f ) ，