结构温度应力仿真计算中的应用研究.pdf

长江科学院 硕士学位论文 数值流形法在大体积混凝土结构温度应力仿真计算中的应用研 究 姓名：明峥嵘 申请学位级别：硕士 专业：岩土工程 指导教师：林绍忠;苏海东 20070601 长江科学院硕士学位论文 i 摘 要 水工混凝土结构往往体积庞大、结构形式复杂，用有限元法对其进行温度场 和温度应力仿真计算时，面临着计算工作量大和前处理困难等问题。数值流形法 具有自适应分析方便、网格切分容易等优点，为大体积混凝土结构的温度应力仿 真计算开避了新的途径。 本文将数值流形法应用于大体积混凝土结构温度场和温 度应力仿真计算,在已有研究成果的基础上开展了如下研究： (1) 根据有限元法的混凝土徐变递推公式，推导了适合于数值流形法的混凝 土递推公式。 (2) 参考有限元法的等效热传导方程和残留比法，推导了适合于数值流形法 的通水冷却计算公式。 (3) 开发了基于高阶数值流形法的大体积混凝土结构温度场及温度应力仿 真计算程序，并通过多个算例验证了公式和程序的正确性。开发的程序可以分析 二维和三维问题。 以上研究工作为数值流形法在大体积混凝土结构温度应力仿真计算中的应 用奠定了坚实的基础。 关键词： 数值流形法 仿真分析 大体积混凝土 温度场 温度应力 徐变 通水冷却 长江科学院硕士学位论文 ii abstract because of huge body and complicated profile of structures much computation is required and difficulty in mesh generation is encountered in the thermal stress analysis simulating the construction process of the hydraulic concrete structures by using fem. being convenient in adaptive analysis and mesh generation, numerical manifold method (nmm) provides a new approach for the thermal stress analysis. in this paper nmm is applied to the thermal stress analysis and the following works are carried out based on the existing researches: (1) according to the recursive formulae of concrete creep based on fem, the recursive formulae suitable for nmm are deduced. (2) referring to the equivalent equation of heat conduction and the remainder ratio method based on fem, the formulae suitable for nmm to model the pipe cooling is derived. (3) upon the derived fomulae, programs of high-order nmm for thermal field analysis and thermal stress analysis simulating the construction process of 2-d and 3-d mass concrete structures are developed and verified by examples. above works provide a sound base for the application of nmm in the thermal stress analysis simulating the construction process of the mass concrete structures. key words: nmm, simulation analysis, mass concrete, thermal field, thermal stress, creep, pipe cooling 长江科学院硕士学位论文 1 第一章 绪论 1.1 引言 大坝、船闸、泄洪建筑物、电站厂房、大型机器设备的基础等大体积混凝土 结构，在现代化建设中发挥了不可替代的作用，尤其是在水利水电建设事业中。 在大体积混凝土结构中，温度应力往往是其他所有外部荷载产生的应力的总和， 甚至还要多。大体积混凝土结构属于高次超静定结构，水泥水化热及环境温度波 动引起结构内部温度场的变化， 在受到地基及结构自身的约束时就会产生温度应 力。混凝土是一种脆性材料，抗拉强度低，如果产生的拉应力超过混凝土的抗拉 强度，结构就会开裂，甚至发展成深层裂缝或贯穿裂缝。因此，为了保证结构的 安全和稳定，必须采取措施控制温度和防止裂缝，而温度应力仿真计算以及结合 具体条件的施工过程的仿真模拟能为温控措施的制定提供依据， 这也是必不可少 的。 1.2 温度应力计算的常用方法 早在上世纪五十年代起，我国的科技工作者就开始对温度应力计算方法的 研究， 并取得了很多重要的成果， 成功地应用到很多大型水利工程的仿真计算中， 为大体积混凝土结构温度应力仿真计算的进一步发展打下了坚实的基础。 常用的 温度应力分析方法包括解析解法、差分法、有限元法、边界单元法等。 1.2.1 解析解法 对于比较简单的情况， 比如无限大的平板温度场的求解等， 可以求出理论解。 张子明、傅作新1提出了模拟碾压混凝土坝成层浇筑过程的一维温度场的解析 解。张德兴、李梦佳2考虑地基对板的双向约束以及非均匀受力的计算模型，竖 向约束采用文克尔假定，水平约束采用折线假定，推导出有关地基上大体积混凝 土温度应力计算的解析方法。王振波3等将复杂的三维热传导问题转化为简单的 一维热传导问题，得到瞬态温度场的解析解，为工程实际应用提供了一种简捷有 效的计算方法。 长江科学院硕士学位论文 2 但对于复杂问题， 求解解析解是非常困难的， 通常需要用其他数值方法求解。 目前求解混凝土温度场和温度应力的数值方法主要有基于对热传导方程进行差 分求解的差分方法系列和基于泛函变分原理的有限元方法系列。 1.2.2 差分方法 差分求解系列包括:普通差分系列和破开算子法等。 普通差分方法求解温度场，需要在时间上划分步长，在空间上划分网格。空 间上的网格，通常要求结点之间按维等距。此外，不同的差分格式，对时间域的 步长有不同的限制，在有限元方法获得广泛应用之前，差分方法是求解温度场问 题的主要手段。 破开算子法是一种数值计算中的分步解法，通过引入一个或若干个中间变 量，将偏微分方程中的时间微商破开成为两个或更多的部分，从而得到相应多个 一维空间（或若干维）的偏微分方程，可利用不同的差分格式求解，使计算大大 简化，适用于求解多维空间的不恒定场问题。这种方法最初由苏联学者提出。董 福品、朱伯芳等4将这种方法运用于混凝土温度场和温度徐变应力方面的研究， 并给出了相应的公式。计算的实例表明，用破开算子法求解温度场和温度应力， 其精度和速度均不亚于有限元法等数值方法。 在基本的差分方法的基础上， 还衍生出了各种基于差分原理的温度场温度徐 变应力求解方法。但应指出，由于差分方法对空间结点距离的特殊要求，难以适 应实际结构形状复杂多变的情况， 它所要求的空间形状通常是比较规则的。 此外， 某些差分格式对时间步长也有特殊的要求， 难以适应仿真分析中为了提高效率进 行变步长的要求，从而限制了各种优化方法的采用。另一方面，采用差分方法求 解温度场、温度徐变应力不利于编制集成化的计算分析软件。 1.2.3 有限元法、边界元法 基于泛函变分原理的温度场温度应力数值求解方法主要包括:有限元法、边 界元法等。 有限元法的基本思想最早出现于 20 世纪 40 年代初期，但是直到 1960 年美 长江科学院硕士学位论文 3 国的克拉夫在一篇论文中首次使用“有限元法”这个名词。在 20 世纪 60 年代末 70 年代初，有限元法在理论上已经基本成熟，并开始陆续出现了商业化的有限 元分析软件。有限元法的基本思想是将连续的结构离散成有限个单元的集合休， 同时选定场函数的结点值作为基本未知量， 并在每一单元中假设一近似插值函数 以表示单元中场函数的分布， 进而利用力学中的某些变分原理去建立用以求解结 点未知量的有限元方程， 从而将一连续的无限自由度问题转化为离散域中有限个 自由度的问题， 一经求解就可以利用解得的结点值和假设的插值函数确定单元上 以至整个集合体上的场函数。有限元法作为一种发展相当成熟的数值分析方法， 在数值计算的各个领域都有着成功的应用。 基于泛函变分原理推导出的有限元法 计算格式，可以对各种复杂的实际情况，在不同的计算网格条件下求解得到不同 精度的位移解答，当网格划分较为理想时，其解答能较好地符合或逼近真实解。 它所求得的温度场和位移场都具有一阶连续的特性， 能够满足一般工程问题的需 要。对于瞬态温度场问题，只需与隐式的差分格式相结合，有限元法即可不受时 间步长的限制。由于易于编制大规模集成化的有限元分析软件，有限元法在结构 位移、应力分析等方面有着相当成熟和成功的应用，也是目前在大体积混凝土结 构仿真分析中被广泛采用的数值方法。 边界元方法是一种很有前途的数值计算方法，其特点是仅对边界进行离散 化，在计算域内并不离散化，之后利用无限介质中点荷载的解答，对所有边界上 的单元求和，得到代数方程组，这样就降低了问题的维数，二维问题可以采用一 维单元，三维问题可以采用二维单元，大大减少了单元的数目，方便了单元的剖 分，而且数据的处理工作与代数方程的阶数显著减少。此外，由于利用了无限介 质中的点荷载理论的解答，可以直接将此方法用于无限及半无限介质。边界元法 的缺点在于，它所形成的刚度矩阵往往是一个满秩矩阵，无法利用其稀疏性。边 界元法在求解温度场及温度应力已经有所应用，实例表明效果良好。但是作为一 种有待于进一步深入发展的数值分析方法，其单独大规模的应用实例还不多，往 往是和有限单元法耦合求解一些特殊的边界问题。 目前有限元法是应用最为广泛的数值方法。虽然经过几十年的发展，有限 元分析技术已相当成熟，并且广泛应用于各个领域中，取得了很好的效果，但在 用于大体积混凝土结构温度应力仿真时， 其计算规模大、 网格划分困难等问题就 显得非常棘手，一直以来也未能得到很好的解决。为获得形态满意的网格，必须 长江科学院硕士学位论文 4 耗费大量的人力物力进行网格剖分， 时间消耗相当大， 对计算机硬件要求也比较 高。 为克服这些困难， 国内外科技工作者做了大量的工作， 并取得了丰硕的成果， 其中包括：在方程组求解方面，林绍忠、苏海东5, 6, 7采取刚度矩阵局部分解技 术和迭代法，提高方程组求解速度；在减少计算步数方面，朱伯芳8, 9, 10提出了 徐变分析的隐式解法和分区异步长算法；在合并单元方面，朱伯芳11, 12提出了 并层算法，王建江13及李克量14等提出非均质单元方法进行碾压混凝土坝 （rccd）仿真，林绍忠15等提出了超级有限元法，chen yaolong16，张建斌等 17提出了浮动网格法，凌道盛等18提出了层合单元概念。这些合并方法的提出 大大减少了计算单元并显著降低了计算规模， 但有的要配以满足协调性要求的单 元衔接措施后才能用于分块浇筑情况。 因此，在保证计算精度的前提下，如何减少工作量，提高计算效率，是当前 大体积混凝土结构温度应力仿真计算面临的主要问题。 1.3 通水冷却的模拟方法 1931年美国垦务局在欧瓦希（owyhee）拱坝进行了混凝土水管冷却的现场试 验，结果比较满意。两年后胡佛(hoover)拱坝开始施工，全面采用水管冷却，效 果良好，此后在全世界得到广泛应用。我国于1956年兴建第一座混凝土拱坝 响洪甸拱坝时，开始采用水管冷却，以后成为混凝土坝温度控制的重要手段。 在长期经验积累及理论研究过程中， 人们对人工冷却的施工工艺及计算方法 都有了更深刻的认识。 以采用钢管或铝管等金属冷却水管为主要特点的人工冷却 近年来已有被高强聚乙烯水管逐渐取代的趋势， 计算理论也逐渐从最初的理论模 型建立而转向探求便于工程实际应用的可行的计算手段。 美国垦务局19研究了二期冷却的计算方法， 用分离变量法得到了无热源平面 问题的严格解答和空间问题的近似解答。 朱伯芳19, 20, 21, 22, 23, 24, 25, 26曾先后得到了 有热源平面问题的严格解答和空间问题的近似解答， 并提出了水管冷却的有限元 分析方法、非金属水管冷却计算方法，及水管冷却效果的近似解法等效负热 源法。麦家煊27将水管冷却理论解与有限元方法相结合，两者取长补短，在提高 计算精度的同时， 节省了大量机时， 为通水冷却的求解提供了一条新思路。 刘宁、 刘光廷28提出了水管冷却效应的有限元子结构模拟技术， 将冷却水管所在单元视 长江科学院硕士学位论文 5 为子结构，给出了相应的有限元法计算公式和计算步骤，解决了计算机储存量、 单元网格的处理等一系列问题， 不失为一种有效的求解途径。 蔡建波29用杂交元 实现了有冷却水管的平面不稳定温度场的求解， 赵代深30采用全过程仿真粘弹性 空间的有限单元法对接缝灌浆水管冷却问题进行了计算研究，对冷却水温、冷却 延续天数进行了敏感性分析，提出了只进行一期水管冷却，并在适当“超冷”情 况下进行大坝接缝灌浆的建议。陈里红、傅作新31也做了类似的工作。stucky 和derron32研究过水管布置方式对冷却效果的影响。朱伯芳33对聚乙烯水管与 金属水管冷却两者之间的关系进行了研究， 提出了聚乙烯水管等效间距的计算方 法，为两种水管冷却的效果和经济比较提供了依据。 在上述方法中，因等效负热源法用已有的温度场的有限元程序就可计算，十 分方便，因而在工程中被广泛采用。在该方法中，将冷却效果的降温作用当作混 凝土自身的一种特殊的吸热特性，即产生负的热源，这样实际上是将与位置相关 的不均匀降温的温度场描述成与水管埋设位置无关的平均意义上的温度场， 使问 题得到了极大地简化，不失为一种有效的方法，但进一步的精度需求要求科研工 作者们致力于提出既方便于实际工程计算， 又能较精确反应实际冷却效果的计算 方法,不过大体积混凝土温度应力计算的复杂性决定了解决此问题是一个漫长的 过程。文献34基于等效负热源的基本方法，考虑了表面散热对水管冷却效果的 影响，据计算，其影响可达23 oc 。在文献25的基础上，朱伯芳26又提出了考 虑外界温度影响的水管冷却等效热传导方程，使等效负热源法趋于完善。 1.4 数值流形法及其研究现状 1.4.1 数值流形法简介 留美学者石根华 35 提出一种新型力学计算方法数值流形方法 （nmmnumerical manifold method） ，该方法采用两套相互独立的网格，即 反映数值解精度的数学网格和表示物体几何边界和材料分区的物理网格， 将整个 研究区域划分成有限个相互重叠的集合（物理覆盖） 。在各个物理覆盖上独立定 义局部覆盖函数， 通过权函数加权平均得到整个求解域上的总体函数。 两套网格 的交集称为流形元，和有限元一样，是基本计算单元，但可具有更复杂的几何形 状。根据数值流形法的思想，石根华给出了一整套二维问题的计算公式，并发布 长江科学院硕士学位论文 6 了相应的计算程序， 其中采用三角形网格做为数学网格， 局部位移覆盖函数取常 数，在流形元上对单项式函数（即多项式函数的基底）的积分采用其提出的单纯 形积分法36进行精确积分。 1.4.2 数值流形法的优点 数值流形法除了在非连续变形分析和裂缝扩展模拟中无需重新划分网格外， 相对于有限元法在仿真分析中的诸多难点， 还具有以下几个优势， 可为温度应力 仿真计算提供强有力的手段： (1) 数学网格与物理网格相互独立， 两者边界可以不重合， 只要求前者覆盖 后者。这种特性在仿真分析中有着巨大的优势：首先，数学网格采用规则的四面 体或长方体单元，物理网格反映结构边界、材料分区及分层分块，流形元网格的 形成过程是数学网格被物理网格切分的过程， 只涉及相对简单的几何运算， 因此 流形元网格形成的难度要比有限元网格小很多， 速度快很多,且完全是自动划分， 人工只需输入简单的边界信息。 其次,可以考虑仿真分析过程中物理网格的变化， 比如建造和开挖造成的结构边界的改变及早期、 中期、 后期的通水冷却混凝土区 域变化， 可以通过局部修改物理网格来适应这些变化， 使人工输入的数据量减到 最小。 (2) 覆盖函数可以是多项式函数或其他级数，随其阶数的提高或级数项的增 多，计算精度得以提高。采用有限元网格作为数学网格，有限元的形函数就是 流形元的权函数。有限元形函数的协调性保证了整体位移函数的协调性，因此 不同区域、不同坐标方向可以方便地采用不同阶次的覆盖函数，自适应分析方 便。在不同的计算时段也可以调整覆盖函数的阶次，从而达到事半功倍的效果， 可在保证计算精度的前提下有效减小计算规模。 1.4.3 数值流形方法的国内外研究现状 由于所具有的优点，数值流形法提出后就立即吸引了国内外学者的广泛兴 趣。为在我国推广、发展这一方法，国内学者裴觉民翻译出版了石根华37的专 著数值流形方法与非连续变形分析 ，并著文38对该方法给予了高度评价。 长江科学院硕士学位论文 7 目前， 国内外关于数值流形法的研究十分活跃， 针对这些原有公式和程序功能的 改进、发展和应用开展了大量研究工作。 (1) 在单纯形积分方面，林绍忠39应用kronecker kronecker 乘积、hadamardhadamard 乘积 和拉直等矩阵特殊运算， 提出了单纯形积分的递推公式， 使得积分计算量显著减 少， 而且在计算高阶单项式函数积分的同时， 还附带获得所有低阶单项式函数的 积分。 (2) 在数学网格方面，jeen-shang lin jeen-shang lin 40，蔡永昌、骆少明4143、王水 林44、张大林45采用四边形网格，提出二维流形元自动剖分、重分技术；张国 新46、彭自强47分别采用6 结点三角形网格和8 结点六面体网格等二维高精度 网格。 (3) 在高阶覆盖函数方面，彭自强47、田荣48等采用一阶覆盖函数以及lu ming lu ming 49、苏海东50、kourepinis.d kourepinis.d 51等采用任意阶覆盖函数；林绍忠52等采 用矩阵特殊运算技术推导单元矩阵， 简化了高阶流形元公式的推导过程； 居炎飞 53等还探讨了数值流形法的p型自适应分析。 (4) 在应用研究方面，张国新54, 55, 56、yaw-jeng chiou yaw-jeng chiou 57、章光58、王 水林59等模拟裂缝扩展及开展地震破坏分析；takeshi sasakitakeshi sasaki60、王书法61, 62、 曹文贵63进行了基于数值流形法的岩体弹塑性分析和加锚岩体变形分析；骆少 明43等将其应用到金属成型过程分析中；苏海东64, 65等还采用高阶覆盖函数的 形式运用到大变形求解中。 (5) 在三维数值流形法研究方面，刚刚起步，已取得一些成果，如姜冬茹66 和林绍忠67提出的基于两凸多面体求交的流形元网格剖分技术，苏海东65等的 基于ansys的流形单元网格生成技术、骆少明69的面向对象设计的三维弹性分析 程序编程技术、苏海东50的高阶流形元矩阵的子程序代码自动生成技术和郑榕 明70的非连续变形分析探讨，以及林绍忠52等的基于矩阵特殊运算技术的高阶 单元矩阵推导技术及其推导的温度场和温度应力的基本计算公式71，为数值流 形法应用于温度应力的仿真模拟奠定了基础。 (6) 在覆盖函数多项式的形式方面，彭自强72采用局部化覆盖函数，如在 二维情况下用 n i m i yyxx)()( 代替全局覆盖函数的基底 nm yx， 以改善刚度矩 阵的性态。林绍忠73等提出了一种改进的局部覆盖函数 长江科学院硕士学位论文 8 () t i i s i i r i i ntsr tsr rsti c zz b yy a xx uzyxu = + 0, ,，其基底在覆盖区域内均不大于 1，广义位移不会因为数学单元尺寸的大小而在数值上相差很大，进一步改善了 刚度矩阵的性态。 1.5 本文研究目的与研究内容 目前数值流形法的研究工作主要局限于二维问题，主要针对常规荷载作用 下的岩体非连续变形分析和结构应力分析，文71中虽然推导了基于数值流形法 的瞬态温度场和温度应力的计算公式， 为数值流形法在温度应力仿真模拟中的应 用奠定了基础， 但其没有考虑混凝土徐变和水管冷却的影响， 也没有开发模拟施 工过程的仿真计算程序。 数值流形法具有很多优点， 特别是其在适时调整单元精 度和网格划分方面的方便性， 为大体积混凝土结构的温度应力仿真计算开避了新 的途径，但还有一些实现技术问题需要深入研究。为此，本文从以下几方面开展 研究， 为数值流形法在大体积混凝土结构温度应力仿真计算中的应用奠定坚实的 基础： (1) 徐变递推公式的推导 在混凝土徐变的模拟方面，目前有限元法中采用的徐变递推公式74是基于 数值积分点的应力状态， 而数值流形法多采用单纯形积分进行精确积分， 其单元 应力呈多项式函数分布， 这种递推公式已不再适用。 本文在此基础上推导了适合 于流形法的递推公式。 (2) 考虑水管冷却时瞬态温度场计算公式的推导 在通水冷却的模拟方面，目前有限元法中常采用残留比法和等效热传导方 程。本文在此两种公式的基础上推导了适合于流形法的通水冷却计算公式。 (3) 基于数值流形法的模拟施工过程的仿真计算程序编制 基于推导的公式， 编制了温度场和应力场高阶流形法仿真分析程序， 并通过 算例进行验证。 长江科学院硕士学位论文 9 第二章 数值流形法的基本概念和基本公式 2.1 引言 数值流形方法（简称流形法）是石根华博士于1991年提出的一种新的数值方 法。 流形法是基于覆盖技术建立的，它的基木思想是将许多只占局部区域的、相 互重叠的有限覆盖连接成一个覆盖系统，用这个覆盖系统去覆盖求解的区域，在 各个有限覆盖上构造覆盖函数并用覆盖上的权函数连接成总体函数去逼近求解 域的真实场函数。 对于任意形状的区域，都可以将其离散为局部区域的集合，这为有限元法所 采用。同样，对于任意形状的区域，都可以使用有限覆盖组成覆盖系统对其完全 覆盖，这就是流形法的思想。从有限元法对求解域的离散到数值流形法对求解域 进行覆盖，虽然在具体实施中有很多相同或相似的地方，但它在数值方法的概念 上却是一个很大的转变。 流形法采用有限覆盖的思想相对于有限元法的离散方法 有较多的优点。对求解域的离散，有限元法必须依赖于物理区域的形状、边界和 材料分区而进行。而在流形法中，有限覆盖和求解域相对独立，它并不要求有限 覆盖完全符合求解区域的物理边界，只要求其完全覆盖求解域，对覆盖的形状与 范围没有限制，这为处理求解域的复杂构形提供了更大方便。 流形法能统一解决有限元法、dda和解析解法的计算问题，它比传统的数 值分析方法更具一般性，并对边值问题提供了一种更贴近自然的求解方法。有限 元法和dda是流形法的两种特殊情况。 本章介绍流形法的基本概念，并根据文52，71，列出三维流形法的基本公 式。 2.2 流形法的覆盖系统 “流形”一词来源于拓扑流形和微分流形，它是把许多个别重叠的区域连接 在一起去覆盖全部的材料区域， 所以材料的总体位移函数可以用局部覆盖所定义 的函数来近似。流形法采用与物理网格相对独立的数学覆盖，数学覆盖只定义近 似解的精度，而物理网格定义其积分区域。 长江科学院硕士学位论文 10 数学覆盖由用户选择，它的形状任意并相互重叠，规则的格子、圆、级数的 收敛域等都可以作为数学覆盖。物理网格则是材料的固有特征，它主要是材料的 边界、裂隙、不同材料之间的交界面等。数学覆盖和物理网格两者相对独立，但 也有一定的依赖关系，数学覆盖必须覆盖求解对象的整个区域，但其可以超出求 解区域而不必完全与求解区域的边界相符合。 物理网格与数学覆盖共同构成物理覆盖， 如果材料边界或裂缝把数学覆盖划 分成完全隔离的小连续区域，则每一个区域就是一个物理覆盖。也就是说物理覆 盖是物理网格对数学覆盖的剖分。 物理覆盖的公共区域定义为流形元，与有限元一样，是基本的计算单元，但 可具有任意的几何形状。目前，流形法一般采用有限元网格作为数学网格（以后 将构成数学网格的有限元称为数学单元，物理覆盖称为数学网格节点） 。对任一 结点，与该结点相连的所有有限元构成一个数学覆盖（在几何拓扑中称之为 “星”） ，流形元即是有限元网格与物理区域的交集。 举个例子来说明以上问题。如图2.1所示，阴影部分为一连续实体，网格 123541构成了整个数学覆盖系统，五个数学覆盖u1、u2、u3、u4、u5见图2.2所 示，物理覆盖即为数学覆盖与物理区域的交集，如图2.2中阴影区域c1、c2、c3、 c4、 c5所示。 采用有限元网格做为数学网格时， 每一个有限元就是一个数学单元， 图2.3中m1、m2、m3即为数学单元，有限元网格节点1、2、3、4、5即是数学网格 结点。物理区域与数学单元的交集即为流形元，见图2.4所示。 图图2.1 数学网格与物理区域 数学网格与物理区域 fig.2.1 mathematical mesh and physical area 长江科学院硕士学位论文 11 图图 2.2 数学覆盖与物理覆盖 数学覆盖与物理覆盖 fig.2.2 mathematical and physical covers 图图 2.3 数学单元 数学单元 fig.2.3 mathematical elements 长江科学院硕士学位论文 12 图图 2.4 流形元 流形元 fig.2.4 manifold elements 2.3 覆盖函数与单元位移函数 流形法在每一个物理覆盖上必须定义一个覆盖函数，覆盖函数可以是常数、 多项式函数或其它形式级数。覆盖函数与物理边界无关，如果物体只占数学网格 的一部分，覆盖函数仍然是相同的。设定义在物理覆盖ci上的覆盖函数为: = = m j ijcj i i i dzyxt zyxw zyxv zyxu 1 ),( ),( ),( ),( （21） 式中， t , ijijijij wvud =为覆盖函数的待定系数，又称为广义位移； cj t),(zyx是 多项式基底， 如zyx, 1等。 传统有限元法采用0阶覆盖函数， 其 iii wvu,是常量， 即结点位移。对于完全n阶覆盖函数，二维情况下其项数 /22)1)(+=nnm；三 维情况下 /63)2)(1)(+=nnnm。 设流形元e由n个物理覆盖的交集构成，则单元位移函数由n个物理覆盖上 的覆盖函数通过权函数加权平均得到，具体表示如下: = = = n i m j ijij n i i i i i d)z , y, x(n )z , y, x(w )z , y, x(v )z , y, x(u w )z , y, x(w )z , y, x( v )z , y, x(u 111 (22) 其中， = = n i i zyxw 1 1),(； wi n wii zyxtfzyxwtf t 1 ),(),( = = ； t , 21n iiii fffl=f， t 21 , wnwww tttl=t； jiwicjicjij twtntftf tt =； t , 21n jjjwcjj ttttl=tt，nttt wcjj 1, = (23) 长江科学院硕士学位论文 13 以上公式中 i w为权函数， i f 是其系数列阵， w t 是多项式基底集。采用有限 元网格作为数学网格时， i w即是有限元的形函数，n为有限元的结点数。流形法 中数学网格与物理网格不要求重合， 因此可以采用形态简单的有限元来构造显式 权函数，如平面的三角形单元、矩形单元和三维的四面体单元和长方体单元。 2.4 位移的偏导数 从式(22)可见，求位移的导数即要求对 ij n求导数。用)32 , 1(，=rxr分别表示 zyx,，根据偏导原则有 jjijji j i ij ,ij )()( xx n ntcftcf t f ttt oo= = = (24) 式中，符号“o”表示矩阵的hadamard乘积运算，设 rs a=a和 rs b=b是 同阶，则 rsrsb a=bao； j c为列阵，是对基底集 j t 偏导产生的常数； t , 21 g jjjj tttl=t是偏导后的基底集，其中包括 “零”基底， j c中对应 常数为零。如果在 j t中删去零基底，则 j c中也不含相应元素，为减少计算量 和方便编程，可予以保留，但对有关计算采用“跳过”的办法。 2.5 整体平衡方程的建立 2.5.1 流形法的平衡方程 数值流形方法较全面地考虑了参与体系平衡的各种平衡力项。 对于不同的问 题，参与平衡方程的各势能项不同，总的来说主要有: (1)应变能势能；(2)初应力 势能；(3)点荷载势能；(4)体荷载势能；(5)惯性力势能；(6)用于不连续变形分析 的接触弹簧力势能和摩擦力势能等。 在得到各项势能的表达式后，根据最小势能原理，写出系统总势能泛函的极 值条件即为整体平衡方程。其一般形式为 rkx = (25) 其中，k 总体刚度矩阵； x 总体位移列向量； r 总体荷载列向量。 长江科学院硕士学位论文 14 由分析以上过程可知， 流形法的整体平衡方程的建立与有限元法基本上是相 同的，单元刚度矩阵向整体刚度阵组装过程也与有限元法完全相同。以下给出单 元刚度矩阵及等效荷载矩阵的计算公式。 2.5.2 单元刚度矩阵 流形元刚度矩阵 文52，71应用矩阵特殊运算技术进行单元分析，推导了高阶流形法的单元 矩阵（包括单元刚度矩阵、弹簧支撑刚度矩阵、等效荷载矩阵等）公式，将被积 函数表示为便于单纯形积分的矩阵相乘形式。 与有限元类似，流形元刚度矩阵的子块为： = v klijijkl vd tdb bk，mljnki1,1,=； (26) 式中， pq d=d为弹性矩阵；b为几何矩阵。如三维实体元，其 ij b如下： = = = 3 1 13 23 12 3 2 1 0 0 0 0 0 0 0 0 0 ,ij ,ij,ij ,ij,ij ,ij,ij ,ij ,ij ,ij ij n nn nn nn n n n a (27) = = = 0 0 1 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 , 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 , 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 321 aaa (28) 将式(27)代入式(26)得 = = v klijijkl vnnd)( )( 3 1 , t 3 1 , adak = = 3 1 3 1 t d v ,kl,ij vnndaa (29) 321 ,aaa起到选择矩阵的作用，因此 ijkl k的(r,s)元为： + = = )s()r( v ,kl,ijpqrsijkl kq,kp, vnnd)( 1313 3 1 3 1 dk (210) 长江科学院硕士学位论文 15 式中，下标集),()k,k ,k(356524641 921 =l。 同理，对于二维实体元,下标集合 ),()k ,k ,k ,k(2331 4321 =， + = )s()r( kq ,kp 1212 ，21=，21=。 如果在积分域v内， pq d是常数，则 = = 3 1 3 1 kd)( pqrsijkl k (211) = v ,kl,ij vnnkd )(v)( lk v ljji cfttcfoo d tt = (212) 式中， j t是对基底集 j t在方向上的偏导数， t l t类似； 被积函数 （ t lj tt） 是一矩阵，但其中各元素仅是两个基底偏导后的乘积，它们的阶次按文71确定 后，就可以应用文39中的积分公式进行精确积分。以下有关公式都是如此。 若为各向同性体，其拉梅常数为和g，式(211)进一步简化为： = += 3 1 )( kgkgk rssrrsrsijkl k (213) 式中， rs 为kronecker 符号。 )21)(1 ( + = e ， )1 (2+ = e g (214) 弹簧支撑刚度矩阵 传统流形法采用一系列孤立的硬弹簧来模拟边界约束条件。设弹簧刚度为 p k，支撑点位于某单元内点),( 000 zyx，弹簧方向余弦矩阵为l，则其刚度矩阵 的子块为： plkjipklijijkl z ,y,xz ,y,xnnktftfkk)()( 000 t 000 t =； (215) llk = 0 0 0 0 0 0 0 0 t p p k (216) 特别地，若弹簧方向与坐标轴 r x方向一致，则 rsprsp k=)(k (217) 长江科学院硕士学位论文 16 当同一点有多个方向的位移约束时，上式可叠加。 在高阶覆盖函数情况下，对于表面位移约束条件，支撑点点数的不确定性会 给计算精度带来一定影响。为适应高阶流形法分析，可采用分布式硬弹簧来模拟 某单元边界面a上的约束条件，其刚度矩阵的子块为面积分形式: ak a ljia a klijijkl ddnnkfttfkk =a)(a tt ； (218) llk = 0 0 0 0 0 0 0 0 t a a k (219) 式中， a k为单位面积的弹簧刚度。 杆单元刚度矩阵 设杆单元的两端分别位于流形元 1 e和 2 e内的点),( 111 zyx和点),( 222 zyx处， 杆轴线的方向余弦矩阵为l，则其刚度矩阵的子块为 , , , 2121222111 b e kl e ij ee ijklb e kl e ij e ijklb e kl e ij e ijkl nnnnknnkkkkk= (220) llk = t t t 0 0 0 0 0 0 g g e l a b (221) 式中，),(),( 222 t 111 t 21 zyxnzyxn lk e klji e ij tftf=分别为流形元 1 e和 2 e所属数学 单元的 ij n在点),( 111 zyx和点),( 222 zyx的值；la,分别为杆的面积和长度， t ge, 分别为杆的弹模和等效剪切模量，后者通常为零。 2.5.3 单元等效荷载矩阵 体积力等效荷载 设体积力 t , zyxv qqq=g，则其等效荷载矩阵的子块为： v v ji v ijij vvngtfgr d d t v = (222) 面力等效荷载 设分布力 ffzyxs ppptpp= t ,作用于流形元表面a上， f p为系数矩阵， 长江科学院硕士学位论文 17 f t为基底集，等效荷载列阵的子元为： i a jff a ijij aanfttppr = dd t s (223) 形状复杂的表面积分区域a可分为若干个空间三角形分别进行单纯形积分 后求和。 设三角形结点为1，2，3（123指向体外） ， 其坐标为3 , 2 , 1,=izyx iii ， 如果面力为指向体内的法向分布力，则 f p为行向量，等效荷载为： = xyyx zxxz yzzy ijf a fij baba baba baba da a fttpr )( 2 1 t (224) 式中， 323232 313131 , , zzbyybxxb zzayyaxxa zyx zyx = = 222 )()()( 2 1 zxxzyzzyxyyx babababababaa+= 初应力等效荷载 初应力等效荷载 设初应力 st= 0 , s为系数矩阵， t为基底集，等效荷载的子块为： )( dd t 3 1 t 0 t = = v jij v ijij vv cfttsabro (225) 其r方向的子元为： = + = 3 1 13 )()( )r( krij rr； = v jij v)( d t cfttsro (226) 初应变等效荷载 设初应变 et= 0 ,e为系数矩阵， t为基底集； 如果在积分域v内，d是 弹性矩阵，为一常数阵，则等效荷载的子块为： )( dd t 3 1 t 0 t = = v jij v ijij vv cfttdeadbro (227) 其r方向的子元为： = + = 3 1 )1(3 )()( r krij rr 长江科学院硕士学位论文 18 )( d t ji v j vcfttdero = (228) 集中力等效荷载 设集中力 t , zyx ppp=p作用于点),( 000 zyx， 其等效荷载矩阵的子块为： ptfpr),( 000 t zyxn jiijij = (229) 2.5.4 单元应力 计算获得位移 ij d后，单元应力按下式计算： 00) (d+= = = n i m j ijij 11 db (230) 将式(27)代入式(230)得 jji m j n i ij tcfda)( t 11 3 1 o = = (231) 上式中，列阵 j t是基底集，其它均为为系数矩阵。 2.6 本章小结 本章简要介绍了流形法的基本概念，并根据文52，71介绍了用流形法分析 位移场的相关计算公式，作为下文章节温度场及温度徐变应力分析的基础。 长江科学院硕士学位论文 19 第三章 用流形法分析温度场及温度应力 文71中已经推导出用流形法分析温度场和温度应力的计算公式，本文是在 其基础上进一步将流形法应用于大体积混凝土结构温度应力的仿真分析中， 故有 必要在此对有关基本公式做一回顾。 3.1 热传导原理 3.1.1 热传导方程 三维热传导方程： d d z t y t x t a t + + + = 2 2 2 2 2 2 (31) 式中， t 混凝土温度； 时间； a 导温系数, ca/=； 混凝土绝热温升； 混凝土导热系数； c 混凝土比热； 混凝土密度 3.1.2 边界条件 第一类边界条件：混凝土表面温度t是时间的已知函数，即 )()(ft= (32) 式中，)(f为已知函数。 第二类边界条件：混凝土表面的热流量是