介老师第一课有限元计算探讨

1、PAGE PAGE 6有限元编程和计算中的若干问题探讨介玉新（清华大学水利水电工程系，北京，100084）【摘要】针对有限元编程和计算中的若干问题进行探讨，指出了一些应当注意的概念和细节，比如坐标轴和符号约定、流量边界和流速边界名称的辨析、固结分析中流量边界的处理等。也指出了数值积分中的一些错误之处。结合编程和计算经验，给出了程序不能正常运算时的一些解决方法，以及人们对有限元认识的误区与原因分析。【关键词】有限元；固结分析；数值积分；程序1 前言笔者近年来主要从事土工数值计算等工作，进行有限元编程和数值分析。在科研和教学中有一些简单的经历和体会。在此谈一下，希望能够抛砖引玉，与从事相同工作的同

2、仁共勉。文中谈到的问题相信大多数人都曾经发现过，只是没有说出来而已。2 坐标轴及符号约定xyz图1 坐标轴一般弹塑性力学分析中约定坐标轴为右手螺旋方式，如图1。三个坐标轴方向满足，。在编程中应当遵循这种隐含约定，否则可能会使面积计算出现负值以及其它一系列的问题。xyz图1 坐标轴在堤坝二维渗流分析中，对xy坐标系常常隐含约定y坐标为竖直方向。在一般地下水二维渗流分析的文献中，xy坐标则常常定义在水平面内。注意到这些细微差别，对阅读相关文献会有好处。土力学中应力应变的符号规定与一般弹塑性力学符号规定相反。在编程时，建议采用一般弹塑性力学的符号规定。虽然这种规定使固结分析时土力学中习惯为正的孔隙水

3、压力变为负号，但能够使之在整个有限元计算中与一般弹塑性力学中的定义保持一致。如果将孔隙水压力按土力学符号定义，则容易造成公式推导和编程中的很大困扰。原因是在连续介质力学中坐标轴、边界荷载、应力分量以及位移、应变等等是在统一的正负号约定下定义的，牵一发而动全身。要尽量保持符号约定的一致性。3 流量边界还是流速边界在渗流计算中，第二类边界条件为流量边界，常写为 ， 或 （1）式中q（或）为单位面积上的渗流量。、为边界外法线的方向余弦。一些人觉得既然是“单位面积”上的流量，为何不称作“流速”边界。原因在于“流速”是矢量，“流量”是标量。边界条件中“单位面积上的渗流量”其实是流速在边界外法线方向上的分

4、量，等于，即与外法线方向单位矢量的点积。所以将此边界条件称为“流量”边界比叫做“流速”边界更贴切些。 很多貌似合理的直观感觉其实是经不起推敲的。比如人们总以为制造衣服或地毯的纤维丝比土软，其实其杨氏模量是土的几百至上千倍。又如为了说明自然界存在很大的吸力，人们常拿上百米高的大树举例，但如果树木等将水份输运到顶部只是依靠力学上的吸力，那么树上的枯枝为什么没有得到水分？还有一旦植物死去，为什么哪怕一棵小草也会干枯呢？4 数值积分中的错误笔者在编程时，发现二维三角形的Hammer积分当积分点数为4时，权系数应分别为-27/96和25/96，而不应当是书中所写的-27/48和25/48。两者相差2倍。

5、笔者曾就此询问文献1的第一作者，作者的回答是在写书时其实注意到了这个问题，但是看到国内外很多书中都是如此，所以就照直采用。这个事情让人费解。笔者核对了许多书，发现几乎所有的书本都在这里犯了错误。尊重原作者固然是一种美德，但容忍其错误却并不值得提倡。其实这个错误很容易就可以发现。不需推导公式，只需将被积函数取为1，就能轻易地发现问题。相信第一个得到这个公式的人应当是正确的，但在辗转反抄的过程中何时引用错了就不得而知。相信这许多年来也有不少人发现了这个错误，但不知为何大家都没有在书中将它改过来。文献2中作者用“*”做注，指出了正确的权系数计算值，但仍没有直接在原表格中进行更正。文献中的错误在所难免

6、。ASCE的Journal of Geotechnical and Geoenvironmental Engineering 中也有不少打印错误。在程序编写时，将用到的公式亲自推导一遍，就能避免很多麻烦。另外，一些简单方法也有助于方便地发现错误，比如，a 量纲分析：检查量纲是否匹配；b 检查对称性：公式是否对称，x、y、z方向三个分量轮换是否改变公式含义；c 用特殊数值进行检验：比如取变量为0、1、或等进行分析，就可以检验公式是否有误。 另外，多找几本书比较着学习，也有利于发现问题，避免错误。5 固结计算中的问题比奥固结理论的平衡方程为： （2）其中p为孔隙水压力。上述及以下的符号规定均按一般

7、弹塑性力学的符号规定。连续性方程， （3）其中vx、vy、vz为流速，v为体应变。根据达西定律，有， ， 其中h为总水头，（式中负号是由于孔隙水压力p按一般弹塑性力学定义所引起）。由此代入式（3）可以得到以孔压表示的连续性方程： （4）以式（2）、（4）为基础可以推导出相应的有限元格式： （5a） （5b）注意连续方程的流量边界条件，用孔压表示变为 （6）有些书中漏掉了上式中的项，将有限元格式（5b）的右边项写为 （7）其中，为形函数。一些程序中也犯了上述错误。这个错误虽然幸运地尚不影响对超静孔压问题的计算，但将导致在考虑总孔压的计算中在不透水边界处有流量不为0的错误。事实上，不但在不透水边界

8、上，而且在不同的内部单元边界上也会出错。相信很多人都曾发现了这个问题。但在修改上一些人可能走了弯路。考虑项后，式（7）应当改写为 （8）于是，除了边界流量外，需要针对每个单元边界计算。除了在整个计算域的边界上要对该式进行计算外，在不同的内部单元边界上也要进行这项操作（在相同的内部单元边界上倒不必处理）。这种计算是相当繁琐的，且极易在编程中出现失误。其实只要利用高斯公式将边界积分改写为体积分，即将式（8）进一步写为 （9）其中 ， 就能简便地解决这个问题。因为体积分计算要比边界积分计算方便许多。6 程序不能正常运算的解决办法笔者也常遇到有人询问程序算不下去的情况。此时往往是输入文件有误。在确认输

9、入文件正确的情况下，应当检查单元和节点数目是否超过相关数组所定义的数值。还有一种情况就是被0除。遇到分母为0的情况，程序往往中断运行。编程中消除被0除的情况可以采用如下方法：对于这种算式，如果分母b 0时，同时有a 0和y 0，此时可以在上述算式前加上一句“时，取的语句”（其中为一小量，式中正负号根据b本身的正负号确定）。这样，将b取为一小量，可以避免被0除的情况；如果b 0时，y趋于一常数或，程序中可以按和两种情况用“ifelseendif”等选择语句分别加以处理；如果时处理比较繁琐，且b=0的情况非常罕见，也可以不予处理，只在算式前加上诸如“如果，输出警告信息”等语句，以便于一旦程序算不下

10、去时排查出错原因。对于等算式，在计算前加上“时，取”及“时，取”的语句。将x限制在其定义域内，也能避免程序算不下去的困扰。同时将x超出其定义域的情况输出出来，便于计算完成后分析其对整个计算的影响。对于等算式，也有类似对定义域要求的情况需要处理。对于应力计算中有时遇到的，表面上是一个简单算式，只要分母不为0即可正常计算，但考虑到该式其实等于，其中为Lode角，所以也应当限制其值域在-11之间，以避免导致其它意外错误。总之，编程时只要注意处理好分母为0、定义域、值域、数组越界等情况，相信能够解决绝大多数程序算不下去的困扰。养成良好的编程和应用程序的习惯，按部就班，一般就不会出现疏漏。7 对有限元认

11、识的误区及分析高等数学中微积分的本质是由连续体中一个点或单元体的信息去推求二维（面）或三维（体）的特征，其纽带是偏微分方程。求解偏微分方程的方法可以是解析法，也可以是数值解法。在计算机出现之前，解析方法占有主导地位，被数学家和力学家发挥得淋漓尽致。计算机的发展带动了数值方法的发展。应当说有限元的发明是数值方法的最大成就。对绝大多数力学问题，有限元法是当前惟一有效的方法。微积分是解决由点到面（体）推广的最有效的数学解析工具，有限元则是解决这一问题的最成功的数值解算手段。有限元的本质是求解偏微分方程问题，即用数值方法解偏微分方程。对土工数值计算来说，除了好的程序之外，正确的计算还需要把握好以下各个

12、环节：对问题本身的理解是否清楚？对程序的特点是否真正掌握？试验和参数是否恰当，可信度有多少？对问题的简化是否抓住了主要矛盾且符合实际？有限元网格的划分是否合理？对边界条件的处理是否恰当？数据文件是否严格按程序要求准备，并正确输入？输出成果是否违背自己的预期？任何一个环节出错，都有可能导致求解错误。目前对有限元有一些误区，一种是事事都希望用有限元算一下，另一种则是对有限元加以排斥，说它是糊弄凑数。但总的来讲后一种情况可能更多一些。分析其原因，可能在于：（1）因有限元公式符号复杂，涉及内容太多，令人望而生畏。由畏生怨，因此予以拒绝排斥。这在有限元发展之初可能较多，随着有限元在本科和研究生教学中的普

13、及，人们大都受过这方面的良好教育，在接受有限元方面已不成问题；（2）不恰当的使用。由于使用者对有限元以及待解决的问题认识不清，或参数问题，或对计算结果不能正确解释，因此导致错误。这可能是有限元招致批评的主要原因；（3）过分期待。理论上有限元可以解决岩土力学中的所有问题，但事实上并不能事事都能解决或解决得很好。土不是钢材那样的弹性材料，单单勘察取样和土的本构关系就会制约数值计算的准确性。过分要求，让有限元做力所不能及的事情，就容易使有限元因出错而受到批评；（4）直觉错误。计算结果与预期或直觉不符，因此产生对有限元的不信任感。这种事情更多发生在生产单位。他们实践经验丰富，所以很多情况下他们的直觉或

14、经验是对的。有限元计算不准往往与参数选择有关，不能怪有限元本身。有的时候，人的直觉以及经验其实也是错误的。另外，工程师们已习惯于对强度问题的计算，对与变形相关的有限元计算，总习惯上要与根据强度计算得到变形结果进行对照。土工有限元计算中施工顺序对变形的影响很大，强度问题计算却不必考虑施工顺序。概念和理解以及对直观现象解释上的错误常常使工程师们非常困惑，这也导致有限元受到一些误解。8 结束语随着计算机以及以有限元为代表的各种计算手段的发展，数值计算已深入人心。人们对有限元计算也不再认为是单纯的数字游戏。事实上，有限元分析从边界条件简化、模型参数选择、输出结果判读等各个环节，也同工程设计施工一样，需要深厚的理论基础和敏锐的经验判断，决不是简单的数据输入和结果输出。把数值计算理想化或简单化都是不对的。近年来，随着商业软件在国内的传播，使用商业软件渐渐成为一种趋势。但目前的商业