第２４卷第２期 黑龙江工程学院学报（自然科学版） Ｖｏ１．２４ Ｎｏ．２ ２０１０年６月 Ｊｏｕｒｎａｌ ｏｆ Ｈｅｉｌｏｎｇｊｉａｎｇ Ｉｎｓｔｉｔｕｔｅ ｏｆ Ｔｅｃｈｎｏｌｏｇｙ Ｊｕｎ．，２０１０ 四节点基面力元（ＢＦＥＭ）模型计算性能研究 侯秋喜，孙占青，董占龙 （北京Ｘ－业大学建筑工程学院，北京１００１２４） 摘要：在单元柔度矩阵表达式的基础上，依据已推导出的空间任意多面体单元柔度矩阵的具体表达式和利用Ｉ．ａ— ｇｒａｎｇｅ乘子法推导出的以基面力为基本未知量的线弹性问题的余能原理有限元控制方程的张量表达式及节点位移 表达式，以四节点四边形单元应用为例，通过大型有限元软件ＡＢＡＱＵＳ对比该算法与传统位移有限元法的精确性 和稳定性。 关键词：基面力；余能原理；有限元；线弹性；ＡＢＡＱＵＳ 中图分类号：０３４３；０２４２．２ 文献标志码：Ａ 文章编号：１６７１～４６７９（２０１０）０２—００４４—０４ Ｒｅｓｅａｒｃｈ ｏｎ ｆｏｕｒ—ｎｏｄｅ ｐｌａｎｅ ｅｌｅｍｅｎｔ ｂａｓｅｄ ｏｎ ｂａｓｅ ｆｏｒｃｅ ＨＯＵ Ｑｉｕ—ｘｉ，ＳＵＮ Ｚｈａｎ－ｑｉｎｇ，ＤＯＮＧ Ｚｈａｎ—ｌｏｎｇ （Ｃｏｌｌｅｇｅ ｏｆ Ａｒｃｈｉｔｅｃｔｕｒｅ ａｎｄ Ｃｉｖｉｌ Ｅｎｇｉｎｅｅｒｉｎｇ，Ｂｅｉｊｉｎｇ Ｕｎｉｖｅｒｓｉｔｙ ｏｆ Ｔｅｃｈｎｏｌｏｇｙ，Ｂｅｉｊｉｎｇ １００１２４，Ｃｈｉｎａ） Ａｂｓｔｒａｃｔ：Ｂａｓｅｄ ｏｎ ｔｈｅ ｅｘｐｒｅｓｓｉｏｎ ｏｆ ｅｌｅｍｅｎｔ ｃｏｍｐｌｉａｎｃｅ ｍａｔｒｉｘ ｄｅｒｉｖｅｄ ｂｙ Ｇａｏ ＹＣ（２００３）ａｎｄ ｔｈｅ ｓｐｅｃｉｆｉｃ ｅｘｐｒｅｓｓｉｏｎ ｏｆ ｓｐａｃｉａｌ ａｒｂｉｔｒａｒｉｌｙ ｐｏｌｙｈｅｄｒｏｎ ｅｌｅｍｅｎｔ ｃｏｍｐｌｉａｎｃｅ ｍａｔｒｉｘ，ｕｓｉｎｇ ｔｈｅ ｔｅｎｓｏｒ ａｎｄ ｎｏｄｅ ｄｉｓｐｌａｃｅ— ｍｅｎｔ ｅｘｐｒｅｓｓｉｏｎ ｂｙ Ｉｉｎｅａｒ ｃｏｍｐｌｅｍｅｎｔａｒｙ ｆｉｎｉｔｅ ｅｌｅｍｅｎｔ ｃｏｎｔｒｏｌ ｅｑｕａｔｉｏｎ ｅｘｐｒｅｓｓｅｄ ｂｙ ｔｈｅ ｂａｓｅ ｆｏｒｃｅ ｄｅ— ｒｉｖｅｄ ｂｙ Ｌａｇｒａｎｇｅ ｍｕｌｔｉｐｌｉｅｒ ｍｅｔｈｏｄ，ｔｈｅ ｐａｐｅｒ ｔａｋｉｎｇ ｆｏｕｒ—ｎｏｄｅ ｂｉｌｉｎｅａｒ ｐｌａｎｅ ｓｔｒｅｓｓ ｑｕａｄｒｉｌａｔｅｒａｌ ｆｏｒ ｅｘ— ａｍｐｌｅ，ｃｏｎｔｒａｓｔｓ ｔｈｅ ａｃｃｕｒａｃｙ ａｎｄ ｓｔａｂｉｌｉｔｙ ｏｆ ｔｈｉｓ ｍｅｔｈｏｄ ａｎｄ ｔｒａｄｉｔｉｏｎａｌ ｄｉｓｐｌａｃｅｍｅｎｔ ｆｉｎｉｔｅ ｅｌｅｍｅｎｔ ｂｙ ＡＢＡＱＵＳ． Ｋｅｙ ｗｏｒｄｓ：ｂａｓｅ ｆｏｒｃｅ；ｃｏｍｐｌｅｍｅｎｔａｒｙ ｅｎｅｒｇｙ ｐｒｉｎｃｉｐｌｅ；ｆｉｎｉｔｅ ｅｌｅｍｅｎｔ；ｌｉｎｅａｒ ｅｌａｓｔｉｃ；ＡＢＡＱｕｓ 有限元法由于能精确模拟具有各种复杂的边界 问题；模拟裂纹扩展时的网格重新剖分问题，以及通 条件和变化的材料性能，有效地分析几何非线性问 过位移的偏导数求解应力而带来的精度损失问题等 题和材料非线性问题而得以广泛应用。基于假设位 等。目前，高性能有限元法的研究仍是国际计算力 移场的有限元法是求解各种结构和同体力学问题的 学界研究的热点问题。 主要方法之一，假设位移场的有限元模型及其应用 ２００３年高玉臣【“。ｊ基于基面力的概念，给出了推 研究已取得很大进展。ＴｏｎｇＥｌｌ提出了杂交位移元 导空间任意多面体单元柔度矩阵显式表达式的思 表达格式；ＷｉｌｓｏｎＥ 等提出了非协调位移元表达格 路，为基面力的概念在余能原理有限元领域的应用 式；卞学磺【３ 建立了基于多场变分原理的杂交元（某 奠定了理论基础。在高玉臣工作的基础上，彭一江 些场变量仅在单元交界面定义）表达式。龙驭球ｌ４］ 等［７ ］基于基面力概念详细推导出空间任意多面体 等提出了广义协调元方法，其长处是在平均位移的 单元柔度矩阵的具体表达式，利用Ｌａｇｒａｎｇｅ乘子法 意义上保证单元问的位移协调，广义协调元自由度 推导出线弹性问题的余能原理有限元控制方程，并 少，精度高，能保证收敛；陈绍春和石钟慈ｌ５］提出了 分别给出了四边形单元和任意多边形单元的应用算 双参数法，其构造单元的过程与广义协调元基本一 例。但针对该方法的计算性能还有待于进一步 致，也保证了单元间必要的连续性。但假设位移的 考核。 常规有限元法在分析一些典型问题时存在明显不 本文将在高玉臣、彭一江等工作的基础上，针对 足：如接近不可压缩问题；大变形情况下的网格畸变 二维问题，进一步研究其性能；并与有限元软件 收稿日期：２０１０—０３ ０４ ＡＢＡＱＵＳ中的Ｑ４单元对比该方法的精确性和稳 作者简介：侯秋喜（１９８４ ），男，硕士研究生，研究方向：结构优化设 定性。 计与可靠性分析． 第２期 侯秋喜，等：四节点基面力元（ＢＦＥＭ）模型计算性能研究 １ 基面力的概念 考虑三维弹性体区域，Ｑ表示变形后的矢径， ３２ （ 一１，２，３）表示Ｌａｇｒａｎｇｅ坐标，则变形后的坐标 架为 Ｑｉ一 ｄ Ｃ‘ ． 为了描述Ｑ点附近的应力状态，在向量 Ｑ 、 出 Ｑ２、 。 上做一个平行六面体微元， Ｑ 、 Ｑ２、ｄ 。Ｑ。所对应的面上的力记为ｄＴ 、ｄＴＯ、 ｄＴＯ，并作如下定义： —ｃＴ Ｌｒ斗 ｄｒＨ — ｄ ’， 一０，’ （２） 这里约定３＋１—１，１—１—３。 式中： （ｉ一１，２，３）称为坐标系ｚ 中Ｑ点的基面 力。有关基面力的理论及应用可参见文献［９］～ ［１２］。 ２ 单元余能ＷＤ的表达式 ２．１余能密度Ｗ 表达式 在小变形条件下，只有两个独立应力的不变数 ‘， 及‘，２，其表达式为Ｌ８．］ Ｊ１一 ：Ｕ， （３） Ｊ２： ： ． （４） 式中： 为应力张量，ｕ为单位张量。 对各向同性材料，余能将是Ｉ， 及Ｉ， 的函数。 如果要求余能为０，应力亦为０，于是余能密度应是 应力不变量的二次函数ｌ８］，即 ｗｃ一 （Ｉ，ｚ一带１２）． （５） 式中：Ｅ为弹性模量，ｖ为泊松比。 将Ｉ， 及．，。带入式（５），可导出用应力表示的余 能密度 、一 ［ｆｉ＇Ｏ＂－ （ ） ］． （６） ２．２单元余能ｗ。的表达式 如上所述，对各向同性材料，当 被忽略时，余能 密度可表示为 ｗ＝＝＝游［ ： 一 （ ）。］． （７） 式中：Ｅ为弹性模量， 为泊松比， 为单元平均应力， 为应力偏量，Ｕ为单位张量。 单元Ｄ的余能为 ｗ。一 ［ ： 一 （ ］．（８） 式中： 为单元的体积。 现在，将式 一 丁０￣ｔｏ代人式（８），可得 一 ｛（ １ ）：（ ）一 ［（ ）：ｕｎ Ｗ。一 ｛（ （Ｗ＠ｒｐ一 ［（ ｔｏ）：Ｕ］ ）． （１０） 根据并矢运算规则有 （ ）：（Ｔｏ ）＝＝＝（ｔｏ・ＴＯ）（ ・ｒｎ）． （１１） 根据向量点积交换律规则，可以得到 （ ）：（ＴＯ ｒｐ）一（ ・ＴＯ）（ｔｏ・ｒ口）． （１２） 若令 一 ・ ，则式（１２）可写为 （ ｔｏ）：（ＴＯ ｒｐ）一（ ・ＴＯ）ｔｏ口．（１３） 又因为单位张量 可写为 Ｕ—Ｐ ｏ Ｐ ． （１４） 式中：Ｐ 为协变矢基，Ｐ为逆变矢基。 故（Ｔｏ ｔｏ）：Ｕ一（ ）（ Ｐ ）．（１５） 根据并矢运算规则，可得 （ ）：Ｕ一（ ・Ｐ ）・（ ・Ｐ ）．（１６） 根据向量点积交换律规则得 （ ｒａ）：Ｕ＝＝＝（ ・ｒｎ）（Ｐ ・Ｐ ）． （１７） 由于 （Ｐ ・Ｐ ）一１， 故（ ）：己，一 ・ ． （１８） 因此，单元Ｄ的余能Ｗｏ的表达式可写为 Ｗ。＝＝＝ ［（ ・ＯＴ）ｔｏ 一 （ ・ｔｏ）。］． （１９） 这样，就将单元的余能Ｗ。用单元面力 和 Ｔｐ表示出来。 注意：式（１９）中包含着求和约定， 为单元的体 积。 ３ 基面力表示的支配方程 设弹性体划分为 Ｐ个区Ｖ（或单元Ｖ），则余能 原理有限元的支配方程是下列余能泛函的约束极值 问题。 系统的余能泛函 Ⅱ： 一 （Ｊ ｗｃ 一ｊ ・丁 ｄｓ）・（２ｏ） 式中：Ｔ‘Ｓ 为位移边界Ｓ 上作用的应力， 为相应 给定位移。 黑龙江工程学院学报（自然科学版） 第２４卷 对常应力单元，单元的余能泛函可写成单元面 力Ｔ的显示表达式 列于表１，并绘于图２；径向位移“ 的数值解，以及 与理论解和Ｑ４计算结果的比较列于表２。其中Ｑ４ ⅡＴ—ｗ 死・ ． 约束条件 （２１） 的计算结果由节点值线性内插得到。 ∑Ｔ ０， × 一０，（在 内） （２２） ～Ｆ一０，（在 上） （２３） （２４） Ｔ（ｖＡ’＋Ｔ Ｂ 一０．（在Ｓ＾Ｊ｛上） 式中： 为单元的ａ面上形心点的已知位移， 表 示单元的已知应力边界， 表示单元 与相邻单 元 之间的边界， 和Ｆ分别表示单元应力边界 上的面力和已知面力，丁 和丁‘ 分别表示单元 和相邻单元 在两者公共边界上的面力。 利用Ｌａｇｒａｎｇｅ乘子法，放松平衡条件的约束， 则构建新的单元余能泛函可写成 Ⅱ（丁， ， ）一ＩＩ（Ｔ）＋ Ｃ Ｃ （∑Ｔ。）十 （ ×ｔｏ）． （２５） 式中：丁为单元面力， 一 １， ２， ３］， 为Ｌａｇｒａｎｇｅ 乘子，面对平面问题 一［ ， ］，而系统的新泛函为 １－ｉ一∑［Ⅱ（Ｔ，２，２ ）］． （２６） 由广义余能原理，泛函的驻值条件可写为 Ⅱ一０． （２７） 由式（２７）可以得到结构关于单元面力丁，以及 单元的Ｌａｇｒａｎｇｅ乘子 和 的一组线性代数方程 组，即结构的余能原理有限元支配方程 —ａ丁 ０。 Ⅱ（丁， ， ） ！．：．．．．．．．．． ， ．．．．．．．一 一Ｏ． （２８） ８２ ａⅡ（丁， ， ４） — ＿～一０．‘ ４ 应用算例 ４．１外壁固定厚壁圆筒承受内压作用问题 如图１所示，厚壁圆桶在外壁ｒ—ｂ处被固定， 在内壁ｒ—ｎ处承受内压Ｐ 作用。计算时取ｌ／４结 构按平面应变问题考虑，其中ａ－－－０．５，６—１，Ｐ 一 １，Ｅ一１，ｖ一０．３。计算所得径向应力 、环向应力 的数值解及其与理论解』” 和Ｑ４计算结果的比较 图１ 受内压的厚壁圆筒及有限元网格 表１厚壁圆桶的应力解 Ｏ Ｏ —Ｏ Ｏ 一０ 一ｌ ０．２ ０．１ ０ ０ １ ０．２ ０ ５ ０ ６ ０ ７ ０ ８ ０ ９ ｌ Ｏ （ｂ）ｒ一 图２，一 关系曲线 表２厚壁圆桶的位移解 ４．２悬臂曲梁顶端承受剪力作用问题 如图３所示，悬臂曲梁在上端承受水平剪力作 用。计算时，按平面应力问题考虑，其中取以一１、６一 第２期 １．５、 一１、Ｅ一１、ｖ＝０．３。 侯秋喜，等：四节点基面力元（ＢＦＥＭ）模型计算性能研究 ・４７・ 积分显示形式，不仅表达简洁，而且编程、计算简单 易行，还可避免传统方法因采用数值积分而造成的 精度损失。 ２）与传统的有限元法相比，基于基面力的有限 元计算方法在计算应力时具有更高的精度。 ３）上述基于基面力概念的余能原理有限元法的 基本思想不仅适用于线弹性问题，而且可以推广应 用于结构的几何非线性问题中，具有较好的应用前 图３承受水平力的曲梁及有限元网格 计算所得环向应力 数值解，以及与理论解和 Ｑ４解的比较列于表３，并绘于图４。 表３曲梁的应力解 （ａ）口一４５。的ｒ一 关系 （ｂ）ｒ一１．３２５ ０的 一卯关系 图４ ｒ， 关系曲线 ５ 结束语 研究结果表明，采用本文方法的数值解与理论 解吻合较好，与传统的有限元法相比，本文方法还有 如下优点： １）基于基面力概念推导出的单元柔度矩阵具有 景。 参考文献 ［Ｉ￣ＴＯＮＧ Ｐ．Ｎｅｗ ｄｉｓｐｌａｃｅｍｅｎｔ ｈｙｂｒｉｄ ｆｉｎｉｔｅ ｅｌｅｍｅｎｔ ｍｏｄｅｌｓ ｆｏｒ ｓｏｌｉｄ ｃｏｎｔｉｎｕａ［Ｊ］．Ｉｎｔ Ｊ Ｎｕｍｅｒ Ｍｅｔｈ Ｅｎｇｎｇ，１９７０，２： ７３—８３． ［２￣ＷＩＬＳＯＮ Ｅ Ｌ，ＴＡＹＬ０Ｒ Ｒ Ｌ，ＤＯＨＥＲＴＹ Ｗ Ｐ，ｅｔ ａ１． Ｉｎｃｏｍｐａｔｉｂｌｅ ｄｉｓｐｌａｃｅｍｅｎｔ ｍｏｄｅｌｓ［Ｊ］．Ｎｕｍｅｒｉｃａｌ ａｎｄ Ｃｏｍｐｕｔｅｒ Ｍｅｔｈｏｄｓ ｉｎ Ｓｔｒｕｃｔｕｒａｌ Ｍｅｃｈａｎｉｃｓ，Ａｃａｄｅｍｉｃ Ｐｒｅｓｓ，１９７３：４３—５７． ［３］卞学镄．杂交应力有限元法的研究进展ＥＪ］．力学进展， ２００１，３１（３）：３４４—３４９． ［４Ｊ龙驭球，辛克贵．广义协调元［Ｊ］．土木工程学报，１９８７ （１）：１ １４． ［５］陈绍春，石钟慈．构造单元刚度矩阵的双参数法［Ｊ］．计算 数学，１９９１（３）． ［６］ＧＡＯ Ｙ Ｃ Ａ ｎｅｗ ｄｅｓｃｒｉｐｔｉｏｎ ｏｆ ｔｈｅ ｓｔｒｅｓｓ ｓｔａｔｅ ａｔ ａ ｐｏｉｎｔ ｗｉｔｈ ａｐｐｌｉｃａｔｉｏｎｓ［Ｊ］．Ａｒｃｈｉｖｅ ｏｆ Ａｐｐｌｉｅｄ Ｍｅｃｈａｎｉｃｓ， ２００３，７３（３／４）：１７１—１８３． ［７］彭一江，金明．基于基面力概念的一种新型余能原理有限 元方法［Ｊｌ应用力学学报，２００６，２３（４）：６４９—６５２． ［８］彭一江，金明．基于基面力概念的余能原理任意网格有限 元方法『Ｊ］．工程力学学报，２００７，２４（１０）：４１—４５． ［９］高玉臣．固体力学基础［Ｍ］．北京：中国铁道出版社， １９９９． ［１０］ＧＡＯ Ｙ Ｃ Ｌａｒｇｅ ｓｔｒａｉｎ ａｎａｌｙｓｉｓ ｏｆ ａ ｒｕｂｂｅｒ ｗｅｄｇｅ ｃｏｒｎ－ ｐｒｅｓｓｅｄ ｂｙ ａ ｌｉｎｅ ｌｏａｄ ａｔ ｉｔｓ ｔｉｐ［Ｊ］．Ｉｎｔ Ｊ Ｅｎｇ Ｓｃｉ，１９９８， ３６：８３１—８４２． ［１１］ＧＡＯ Ｙ Ｃ，ＧＡＯ Ｔ．Ｌａｒｇｅ ｄｅｆｏｒｍａｔｉｏｎ ｃｏｎｔａｃｔ ｏｆ ａ ｒｕｂ ｂｅｒ ｎｏｔｃｈ ｗｉｔｈ ａ ｒｉｇｉｄ ｗｅｄｇｅ口］．Ｉｎｔ Ｊ Ｓｏｌｉｄｓ Ｓｔｒｕｃｔ， ２０００，３７：４３１９－４３３４． ［１２］ＧＡ０ Ｙ Ｃ．Ａｓｙｍｐｔｏｔｉｃ ａｎａｌｙｓｉｓ ｏｆ ｔｈｅ ｎｏｎｌｉｎｅａｒ Ｂｏｕｓｓ— ｉｎｅｓｑ ｐｒｏｂｌｅｍ ｆｏｒ ａ ｋｉｎｄ ｏｆ ｉｎｃｏｍｐｒｅｓｓｉｂｌｅ ｒｕｂｂｅｒ ｍａｔｅ— ｒｉａｌ（ｃｏｍｐｒｅｓｓｉｏｎ ｃａｓｅ）［Ｊ］．Ｊ Ｅｌａｓｔｉｃｉｔｙ，２００１，６４：１１ １— １３Ｏ． ［１３］吴家龙．弹性力学［Ｍ］．北京：高等教育出版社，２００１． ［责任编辑：刘文霞］
