有限元 第8讲等参单元2

I b f ( )d b( )d
a
a
问题是 : 如何构造多项式( ),使其对f ( )有最好的逼近?
如果n个结点 i (i 1, 2,..., n) 等距分布，则前面 的插值型求积公式称为Newton-cotes求积公式。
Newton-cotes求积公式具有n-1次代数精度 几个常用求积公式
{R}e 1 1 1 [N ]T {p}t J ddd 1 1 1
在ξ=1的面上受到面力作用，面力移置的公式为:
{R}e 1 1
1 1
[N ]T1{P}t J 1 dd
其中在点 ξ0,η0,ζ0集中力移置的公式为:
{R}e [N ]T(0 ,0 , 0 ){P}
2011/4/19
有限单元法及软件 应用
山东建筑大学 2011
进度安排
10 材料非线性问题 11 几何非线性问题 12 热传导问题 13 有限元Fortran程序设计 14 ANSYS有限元软件 期末考试
1进有度限安元方排法概述
2 数理力学基础 3 简单杆系结构有限元法 4 弹性力学平面有限元方法 5 等参元和高斯积分 6 空间问题有限元法 7 梁结构单元 8 板壳问题有限元法 9 结构动力问题有限元法
– 梯形公式，n=1 – Simpson公式，n=2
b
a
f
( x)dx

ba 2

f
(a)

f
(b)
b a
f
( x)dx

ba 6

f
(a) 4
f
(a
b) 2
f
(b)
k
k 1
7
2011/4/19
写成统一的形式
1 1
f
( )d

n k 1
Hk
f
(k )
u(x, y) u( ,) N( ,)qe 8 3
坐标变换
u(x, y) u( ,) N( ,)qe
x y z





J


x
y
z


x y z

B( ,, ) N( ,, )
示为形函数对局部坐标的偏微分，

N i x N i y


[
J
]1

N i N i

Ni
Ni
z

雅可比矩阵 J
只有在雅可比矩阵可逆 的情况下，可以求解出 Ni,x Ni,y Ni,z！
基本思路是：在单元上选择某些特征点（积分点），求 出被积函数在这些积分点上的数值，然后用一些权函数 乘这些函数值，最后求和就可得到近似积分值。有限元 分析中，最常用的高斯数值积分法。
基于区间内n个结点i (i 1, 2,..., n), 将多项式( )取为Lagrange插值多项式,
n
即
2. 坐标变换式和位移模式
形函数计算
4
2011/4/19
按照上节介绍的等参单元分析的基本步骤可以得到 三维单元的单元刚度矩阵。雅可比矩阵为，
x y z
[J
]


x

y

z



x
y
z
ቤተ መጻሕፍቲ ባይዱ


形函数对整体坐标的偏微分可以用雅可比矩阵表
[K e ] 1 1 [B]T [D][B]t | J | d d 1 1 f ( ,)d d
1 1
1 1
6
2011/4/19
数值积分方法是一种近似的方法。
b a
f
( )d

n i 1
Ai
f
(i )
一个函数的定积分可以通过n个结
点的函数值的加权组合来表示
f
( )d

n i 1
Ai
f
(i )
1 f ( )d 1
n
Hk f (k )
k 1
n 2 则右边四个参数 H1 1 H 2 2
如果左边的 f ( ) 也只含有四个参数，则它们之间的关系就能完全确定。
f ( ) c0 c1 c2 2 c3 3
H k 为求积系数或称为积分系数。
牛顿-柯斯特积分公式精度 n 1
高斯求积法
如果不事先规定积分点的位置，而是允许这一些点位于能得到精度最好的
积分值之处，在给定积分点数目的条件下，这样做可以提高所构造的求积
公式的精度。
1 1
f
( )d

n k 1
Hk
f
(k )
如果规定可以取n个分点，我们必须求出2n个未知量
(k f (k ) k 1, 2,, n)
对它进行精确积分，并用积分结果代替原函数的积分，其误差是 O(2n )
这种求积公式具有代数精度 2n 1
这样求出的数值积分公式称为高斯积分公式。
31
Gauss-Legendre求积公式
n个插值结点非等距分布 结点和积分权系数可以查表
1 1
ANSYS提 供的
Solid95 单元
Ni

1 8
(1


i
)(1

i)(1


i
)(
i
i
i
2)
(i 1,2,...,8)
Ni

1 (1 4
2 )(1 i)(1 i
)
Ni

1 4
(1


2
)(1


i
)(1 i )
Ni

1 (1 4
[K ]e 1 1 1 [B]T [D][B] J ddd 1 1 1
5.6 数值积分
等参单元刚度矩阵的每个元素都是局部 坐标的函数，等参数变换后具有非常复 杂的形式，在有限元程序中不用解析的 办法来计算局部坐标系中的积分，而采 用数值积分方法。
通常采用高斯积分方法计算单元刚度矩 阵中的元素及等效节点载荷列阵的元素。
n
v Ni ( ,, )vi i 1
n
w Ni ( ,, )wi i 1
n
x Ni ( ,, )xi i 1
n
y Ni ( ,, ) yi i 1
n
z Ni ( ,, )zi i 1
三维问题的单元数目大，节点自由度多，导致 计算规模大，对计算机硬件的要求很高。为缩 短计算时间，有许多问题需要采用巨型计算机， 如CRAY，或并行计算机。
雅可比矩阵计算公式
n
i1
[J ] n i1 n i1
N i
xi
N i
xi
N i
xi

N1
N 2


N1
N 2

N1
N 2

n
i 1
N i
yi
n
i 1
N i
yi
n
i 1
N i
yi
n
i 1
N i
zi

n
i1
n
i1
N i
z
i

N
zi

... ...
N n
N n


x1

x2
.
y1 y2 .
z1
z
2

.
...
N
n


x
n
yn
z
n

三维问题有限元法有以下两个主要难点： (1) 单元划分比较复杂无法采用人工方法完成复杂三维实体的 单元划分，需要有功能强大的单元划分程序，从CAD模型直接 生成离散的单元网格。现在的有限元软件可以读入IGES、STL 等格式的图形交换文件。六面体单元的计算精度比较高，但是 对于复杂三维实体无法实现六面体单元的自动划分。采用四面 体单元能够实现单元自动划分，但是四面体单元的计算精度比 较低。

5
2011/4/19
利用雅可比矩阵的行列式，将整体坐标系下的积 分转换为在局部坐标系下的积分。在整体坐标系中的 体积微元为，
dV dx (dy dz) dxdydz
微矢量在局部坐标系中表示为，
dx

x
d
1

x
d2
x
d
3



dy

左侧代入积分： I 2c0 2c2 / 3
I

H1 f
(1)
H2
f
(2 )

i
)(1


i
)
(i 1,...,8)
其中，i ,i , i 为结点的局部坐标。
Solid95的基本单元
二十结点基本单元
与六面体八结点等参单元相比，六面体二十结点等参单 元能更好地适应不规则的形状，计算误差比较小，基本 单元如图所示，其形函数为
ANSYS提供的Solid95单元是六面体二十节点等参单元， 每个节点有代表x、y、z三个方向位移的三个自由度，可 以退化为五面体棱柱、五面体金字塔形和四面体单元。 Solid95单元的基本单元如图所示。
• ANSYS提供的Solid45单元就是六面体八节点等 参单元，每个节点有代表x、y、z三个方向位 移的三个自由度（DOF，Degree of Freedom）， 可以退化为五面体棱柱和四面体单元，单元局 部坐标为r、s、t，六面体八结点等参单元的 基本单元如图所示。
(2）计算规模大
n
u Ni ( ,, )ui i 1
等参数单元
5.1 概述 5.2 等参单元定义的给出 5.3 四节点四边形等参数单元 5.4 平面八节点等参数单元 5.5 六面体等参单元 5.6 数值积分
1
2011/4/19
5.5 六面体等参单元
多数弹性力学问题需要按照三维空间问题来求解。三维弹 性力学问题的有限元法的基本步骤与平面问题的步骤一样，在 分析三维问题时，所选择的单元主要为四面体单元和六面体单 元。每个单元节点上定义有三个位移分量u、v、w。
常用的三维等参单元有六面体八节点等参单元 和六面体二十结点等参单元。等参单元的位移 模式和坐标变化式采用相同的形函数，如上
ANSYS提供的Solid45单元
2
2011/4/19
Solid45的基本单元
八结点基本单元
六面体八节点等参单元的基本单元如图所示，其形函数 为，
Ni

1 8
(1


i
)(1
dy dz
y
z
dd3
y
z
dd
2

y
z
dd3
y
z
dd
1





y
z
dd
2

y
z
dd
1
x y z dxdydz x y z ddd J ddd x y z
y
d
1

y
d
2

y
d
3



dz z d 1 z d2 z d 3



其中
1,
2,
3
为局部坐标系中ξ,η,ζ方向上的单位向量。
最后，用高斯积分计算出单元刚度矩阵。同样， 用上节中类似的公式就可以在局部坐标下完成单元的 载荷移置。体力移置的公式为
( )=
l n1
i
(
)
f
(i
)

a0

a1

an1
n1
i 1
其中lin1( )为(n 1)阶Lagrange插值函数
n
lin1(
)


j 1
ji
j i j
,显然lin1( j )

1 0
,i ,i
j j
(i ) f (i )
数值积分的基本思想
对于一个定积分 I b f ( )d ,构造一个多项式( ), a
使得在区间内n个点(i ,i 1, 2,..., n)上( )与f ( )相同, 即 (i ) f (i ) (i 1, 2,..., n)
则用( )来近似代替f ( ),积分式变为
2 )(1 i )(1 i)
(i 9,11,17,19) (i 10,12,18,20) (i 13,14,15,16)
其中ξi,ηi,ζi为单元结点在局部坐标系中的坐标 单元刚度矩阵为，
K e [B]T [D][B]dxdydz
3
2011/4/19
x( ,) N( ,)xe
则用( )来近似代替f ( ),积分式变为
I
b f ( )d
a
b( )d
a
b a
n i 1
l n1
i
(
)
f
(i
)d
n

i 1
b a
l n1
i
(
)d
n
f (i ) Ai f (i )
i 1
至少具有n-1次 代数精度