非协调六面体有限元分析方法研究

韦建平，韦宏法，韦志安，黄振之

(1.柳州孔辉汽车科技有限公司，广西柳州 545006；2.柳州五菱柳机动力有限公司，广西柳州 545005)



非协调六面体有限元分析方法研究

韦建平1，韦宏法1，韦志安2，黄振之1

(1.柳州孔辉汽车科技有限公司，广西柳州 545006；2.柳州五菱柳机动力有限公司，广西柳州 545005)

基于平面等参协调单元给出空间六面体有限元模型表达式。针对完备的协调单元由于刚度系数值比精确值大、导致在给定的载荷之下计算模型的变形比实际结构小的问题，建立了一套可以满足分片检验条件的六面体非协调单元的有限元分析方法，有效解决了三维协调元为提高计算精度而导致的计算效率低、计算内存占用大的问题。并通过实例验证了非协调元的计算精度。

等参元；非协调元；有限元；六面体网格

0 引言

协调等参元虽有良好的适应性和表达格式的简明性，也因此得到广泛的应用，但是从严格的意义上说，它的精度和效率仍是不够高的。在三维坐标中，一次和二次完全多项式分别是4项和10项。而三维线性单元和二次单元却分别具有8个和20个节点，也即三维等参元中有二分之一的节点自由度对计算精度是无贡献的。因此，E L WILSON提出了非协调等参单元，对改进等参单元计算精度和提高计算效率是很有意义的[1-2]。

作者结合二维非协调有限元分析方法，并在文献[3]提出的协调等参单元的基础上，建立了一套可以满足分片检验条件的三维六面体非协调单元的有限元解法，有效解决了三维非协调元的计算问题。

1 有限元力学基础

有限元法中，需要用到弹性力学的基本方程和与之等效的变分原理。弹性体内任一点的位移可由沿直角坐标轴方向的3个位移分量u、v和w来表示，它的矩阵形式是[2]：

(1)

上式称作位移列阵或位移向量。对于三维问题，弹性力学基本方程包括平衡方程、几何方程(即应变和位移变化关系)、物理方程(即应力和应变关系)，还有力和几何的边界条件，在V域内，其方程一般形式如下：

平衡方程： Lσ+pv=0 (在V内)

几何方程： ε=LTu (在V内)

物理方程： σ=Dε (在V内)

(2)

nσ=ps(在边界V上)

2 协调六面体等参有限元

2.1 六面体单元概述

图1(a)是一个物理空间中的任意四边形单元。在该单元上建立一个参考坐标系ξη。四边形单元的边被横坐标ξ和纵坐标η平分，它们所表示的方程是ξ=±1和η=±1，如图1(a)所示。在该参考坐标系ξη下，该四边形变换成一个边长等于2的正方形，如图1(b)所示。

图1 四边形等参数单元

同理，空间等参数单元可由平面问题等参单元推广得到。空间八节点单元是直棱的，可以描述形状复杂的三维结构体。八节点等参数单元及其参考坐标系ξηζ如图2所示。在参考坐标系下，直棱的六面体是一个边长等于2的立方体，坐标系ξηζ的原点位于形心处[1]。

图2 八节点六面体等参单元

2.2 六面体单元特性

根据等参数单元的含义，它的位移和坐标都采用相同的形函数表示，即

位移函数：

(3)

坐标函数：

式中：ui、vi、wi(i=1,2,3,4,5,6,7,8)是整体坐标系xyz下节点的位移分量，xi、yi、zi是节点的整体坐标， Ni(ξ,η,ζ)是由参考坐标ξηζ表示的形函数。根据矩形单元的形函数公式，可以得出八节点等参单元的形函数Ni的表达式：

(4)

式中：ξi，ηi和ζi是节点i的局部坐标，分别位于六面体的8个角处，其对应坐标值分别为1或者-1，即：

(5)

图2所示的单元节点中，节点编号1-2-3-4-5-6-7-8是逆时针并且由外向内顺序的，而其参考坐标系ξηζ在整体坐标系xyz下的方向是由单元的节点编号决定的。

根据式(2)的几何方程，等参六面体的应变ε公式为：

(6)

式(6)中，由于位移函数uvw和坐标函数xyz不在同一坐标系下，根据复合函数的求导法则，必须经过必要的转换，记Ni,x、Ni,y、Ni,z分别表示形函数Ni对x、y和z的偏导数，它们与参考坐标系下对应的Ni,ξ、Ni,η、Ni,ζ关系式表示如下：

(7)

式中：矩阵J称为雅克比(Jacobin)矩阵。其表达式如下所示：

根据式(6)、式(7)得出等参六面体的应变ε计算公式为：

(8)

式中：B称为单元应变系数，从而单元应力的计算公式为：

(9)

式中：

(10)

其中：常数A1、A2和A3由式(11)定义：

(11)

综上所述，并根据单元节点的刚度公式，三维等参单元的刚度矩阵表达式如下：

(12)

单元刚度矩阵Ke可划分成n×n个子矩阵形式，其计算公式如下：

(i,j=1,2,…,8)

(13)

由于形函数对整体坐标x和y的导数中包含雅克比矩阵的逆矩阵，式(13)在一般情况下很难得到显示方程，必须采用数值积分来进行求解，该方法可参考文献[1]。

3 非协调元的凝聚

3.1 Wilson非协调元

对于三维线性单元来讲，由于单元的插值函数中包含有非完全的ξη、ξζ和ηζ项，在用它表示纯弯曲应力容易出现误差。

而附加项的位移在三维线性单元的8个节点上都取零值时，它对节点位移没有任何影响，而只对单元内部的位移起了调整作用。待定系数α1、α2和α3与单元边界上的节点无关。

包含附加的无节点位移项的单元位移插值表示如下：

(14)

式中：Ni为式(4)所示的形函数。

等式(14)右边由两部分组成，前者为第2节介绍的具备协调性的等参单元位移形函数，后者为非协调部分的形函数，两者矩阵维数不同，计算过程和计算方法是一致的。

因此，为了改善三维线性单元的性质，解决完备协调元中刚度系数值过大，提高其解精度，需将单元的非协调位移插值函数附加到内部无节点的位移项。

3.2 非协调元刚度的凝聚

根据上文所述，结合几何关系和位移函数，加入非协调部分后，单元的应变矩阵公式变为如下式子：

(15)

将式(14)和式(15)等代入位移能泛函并按照通常步骤，此时可以得到位移、刚度与受力的矩阵式[3]：

(16)

其中：

(17)

式中：Pe为外载荷列阵。

对上述关系式进行展开求解，则单元内部位移、刚度矩阵以及载荷列阵之间的关系式

Keδe=Pe

(18)

式中：

此式即为包含附加内部位移项的单元刚度矩阵和载荷列阵。它是在原单元刚度矩阵和载荷列阵内增加了修正项而得到的。经过凝聚后，单元的自由度仍是原六面体单元的自由度，之后的分析计算步骤也跟之前的标准解题方法和步骤一样。

3.3 非协调元的分片试验

为了检验采用非协调元的任意网格划分时是否能达到连续性的要求，需要对其进行分片试验。若能够通过分片试验，则说明解的收敛性能够得到保证。

(19)

4 有限元实例分析

以单片变截面钢板弹簧为例，对其进行有限元模型的建立与分析，并将分析结果与其他分析方法比较，从而得出文中分析方法的可行之处。

单片簧的参数包括板簧的展开长度L=1 250 mm、根部截面长度L1=120 mm、端部截面长度L3=220 mm、根部厚度h2=17 mm、端部厚度h1=9 mm，板簧的宽度w=60 mm，板簧自由半径值R=1 400 mm，以及作用力F=4 500 N、泊松比μ=0.3，弹性模量E=206 000 MPa。以上尺寸参数决定了板簧的截面形状，由于板簧的对称性，下面以其半长作为研究对象，并根据其结构特点将它等效为根部固定端部加载作用力进行有限元分析，如图3所示。

图3 单片变截面钢板弹簧半展长

(1)板簧结构分析及节点编排

依据变截面簧的结构特点，对六面体进行有限元网格的划分、坐标系的确定及网格节点的编排等，见图4。

图4 变截面簧网格单元及节点编号示意图

(2)网格单元定义

针对比较规则的钢板弹簧，只需要八节点的六面体就能很好地描述其三维实体模型，故这里将对空间六面体单元进行网格单元的定义。网格单元如图4所示。

(3)赋予材料属性

弹簧的材料属性为各向同性材料，当材料有弹性对称轴时，弹性矩阵中的独立变量数目将减少，只有弹性模量和泊松比，用于集成弹性矩阵。

(4)添加边界约束

所谓边界约束，就是在某个方向上限制物体的位移或者受力。对于钢板弹簧，根据其安装位置及结构原理，取半边作为分析对象，可看作是简支梁，约束根部节点3个方向的位移，使根部节点1、2、7、8处的X、Y、Z方向位移均为0。

(5)集中力加载

集中力的加载位置为端部沿厚度方向的节点，将载荷力及作用方向按照节点数平均分配到各节点处，比如图5所示的5、11节点处的y方向。

图5 文中与ABAQUS有限元网格模型

(6)有限元模型的求解

利用文献[2]所示的高斯积分法或者Iron积分方法，并基于MATLAB编写有限元模型的计算程序，根据结构参数建立三维六面体网格的有限元模型，最终计算得出有限元模型的刚度以及节点应力值，如图6、7所示。

由表1、2可见：采用协调单元计算结果与ABAQUS计算结果存在较大误差，而采用非协调单元计算结果与ABAQUS计算结果吻合，验证了非协调单元的计算精度。

图6 文中有限元分析结果

图7 ABAQUS分析结果

载荷/N端部位移/mm最大主应力/MPaABAQUS225069.08599.20协调单元[3]225060.16499.36非协调单元(文中)225068.53605.29

表2 与ABAQUS结果的误差分析 %

5 结论

针对空间六面体协调单元进行推导，并通过分析提出该有限元模型的解刚度值偏大的问题。通过引入了非协调项，使得三维八节点单元中的插值函数二次项完全了，这种三维八节点的非协调元在计算中可以达到与三维二十节点协调元同级的计算精度，但是前者的节点数仅仅是后者的2/5，从而使得在有限元的计算分析中，占用最大内存的平衡方程求解运算效率大大提高，而其计算精度也是非常显著的。

【1】王勖成.有限单元法[M].北京：清华大学出版社,2003(1):209-216.

【2】王勖成.有限单元法基本原理和数值方法[M].北京：清华大学出版社,1997.

【3】徐荣娇.结构分析的有限元法与MATLAB程序设计[D].杭州：浙江大学建筑工程学院,2005.

【4】焦兆平.非协调四节点平面等参位移元新列式方法[J].计算结构力学及其应用,1996,13(2):147-156.

【5】龙驭球,辛克贵.广义协调元[J].清华大学土木工程学报,1987,20(1):1-14.

【6】鹿晓阳,刘玉文,许焕然，等.Wilson非协调元的研究与改进[J].力学学报,1989,21(3):379-384.

【7】陈耀明.汽车悬架论文集[M].苏州：苏州大学出版社,2009：32-35.

Study on Finite Element Analysis Method of Incompatible Hexahedral

WEI Jianping1,WEI Hongfa1,WEI Zhian2, HUANG Zhenzhi1

(1.KH Automotive Technologies (Liuzhou) Co.,Ltd., Liuzhou Guangxi 545006,China;2.Liuzhou Wuling Liuji Power Co.,Ltd., Liuzhou Guangxi 545005,China)

Based on the plane iso-parametric incompatible element, hexahedral finite element model expression was shown. Because the stiffness value of complete-compatible element was more large than normal, finally small deformation was lead to. For this purpose, a hexahedral incompatible finite element analysis method was established which could satisfy the fragmentation test condition. And then, low efficiency and big amount of computations for the element of iso-parametric-compatible was solved effectively.Finally, computational accuracy of the incompatible element was verified through example.

Iso-parametric element; Incompatible element; Finite element; Hexahedral mesh

2015-10-26

韦建平(1983—)，男，硕士研究生，从事车辆动力学与控制研究。E-mail:weihf15@163.com。

U464.11

A

1674-1986(2016)02-010-05