北航有限元第4讲 等参元和高斯积分
高斯积分ppt课件
结论:插值型求积公式的代数精度d满足：n-1 d2n-1
5
定理: 若f(2n)(x)在[a,b]上连续，则高斯求积公式的余项为
Rn
f (2n) ()
(2n)!
b a
Hale Waihona Puke (x)wn2(
x)dx
其中(a,b),w(x)=(x-x1)(x-x2)…..(x-xn)。
高斯求积公式的系数Ak恒为正,故高斯求积公式是稳定的.
1
n
f (x)
-1
wi f (xi )
i0
的代数精度最高不
证明：分别取 f(x)=1, x，x2，...xn 时代入公式，并让其成为等式得
A1 + A2 + …… + An =∫ab1dx.= b-a x1 A1 + x2 A2+ …… +xn An =∫abxdx.= （b2-a 2)/2
......
1
3
t(5)
2
A(5) 5
1
3
t(5)
5
0.69314719
积分精确值为
I=ln2=0.69314718… 由此可见，高斯公式精确度是很高的
9
2.Gauss - Chebyshev 求积公式
1 1
f (x) 1 x2
dx
n
Ak
k 1
f (xk )
其中高斯点为Chebyshev 多项式Tn(x)的零点
• Guass求积公式有多种,他们的Guass点xk, Guass系数 Ak都有表可以查询.
6
Gauss - Legendre 求积公式
1
n
f (x)dx
1
第四章 等参数有限元方法
第四章 等参数有限元方法§4.1 引言在有限元法中，提高计算精度的有效方法是提高单元的计算精度。
一种高精度的单元不仅要求采用高阶插值函数的位移模式，而且应能较好地适应物体的边界几何形状。
§4.2 等参数单元的概念本节首先从平面任意四边形单元着手，介绍等参元的一些基本概念。
平面问题的有限元法中，最简单而又常用的是常应变三角形单元，其次是具有4结点的矩形单元。
矩形单元，由于它的位移是坐标的二次函数，因而能比常应变三角形单元较好地反映实际应力的变化情况，但是矩形单元不能适应曲线边界和非正交的直线边界，也不能随意地改变大小，如果能改用任意的四边形单元，如图4-1所示，上述缺点就能克服，但是对于任意四边形单元，如采用矩形单元的双线性位移模式xy y x u 4321ββββ+++=xy y x v 8765ββββ+++= （4.2.1）则相邻单元的公共边界上位移将是不连续的。
这是由于在单元不平行于坐标轴的任意一个边界上，上述位移模式是二次函数C Bx Ax u ++=2而不是线性变化的，因而在边界上的插值函数值不能由同一边界上两个结点的位移唯一的确定。
双线性位移模式不能直接地用于任意四边形单元。
但是对于矩形单元，由于在边界上位移模式是线性的，即在矩形单元中的边界上位移完全由同一边界上的两个端点的位移唯一确定，因而单元间位移是连续的。
也就是说单元是协调的。
为解决任意四边形单元的协调性问题，我们通过坐标变换，首先将xy 平面上的任意四边形单元变换为ξη平面上以原点为中心，边长为2的正方形单元，而xy 平面上的结点1，2，3，4分别映射为ξη平面上的结点1，2，3，4。
这种变换不是对应于整个求解域进行，而是针对每个单元进行的，ηξ,是局部坐标，它只应用于单元范围内，而y x ,为整体坐标，它适用于所有单元。
下面我们考虑局部坐标系下的单元。
设位移模式为ξηβηβξββ4321+++=uξηβηβξββ8765+++=v （4.2.2）每个结点的位移可用位移矢量i α表示，即⎥⎦⎤⎢⎣⎡=i i i v u α )4,3,2,1(=i 每个单元有8个结点位移分量，于是单元结点的位移向量可表示为[]Te v u v u v u v u 443322114321=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=ααααα e α为单元结点位移列阵。
有限元方法--第4章
+
∫∫∫ [N ] {p }dV
T V
dV = J d ξ d η d ζ
{ε} =[εx εy εz γxy γ yz γzx] T=[A] {d}
{ε } = [B1
B2 B3 B4 {δ } = [ B]{δ }
e
[A ]微分算子矩阵 [B ]应变矩阵
]
e
六、应力矩阵
{σ} =[σx σy σz τxy τyz τzx] T=[D] {ε} =[D][B] {δ}e =[S] {δ}e [S ]应力矩阵，[D ]弹性矩阵
二、坐标变换
ζ
5
6
1
8
7
o
x = ∑ Ni xi , y = ∑ Ni yi , z = ∑ Ni zi
i =1 i =1 i =1
8
8
8
η
4
三、位移模式
u = ∑ Niui , v = ∑ Ni vi , w = ∑ Ni wi
i =1 i =1 i =1
8
8
8
ξ
2
3
1 Ni = (1 + ξ0 )(1+η0 )(1+ ζ 0 )(ξ0 +η0 + ζ 0 − 2) 角结点： 8 1 ξi = 0，ηi = ±1, ζ i = ±1：Ni = (1− ξ 2 )(1 +η0 )(1+ ζ 0 ) 4 1 ξi = ±1, ηi = 0，ζ i = ±1：Ni = (1−η 2 )(1 + ξ0 )(1+ ζ 0 ) 4 1 ξi = ±1, ηi = ±1, ζ i = 0：Ni = (1− ζ 2 )(1 + ξ0 )(1+η0 ) 4 其中： ξ0 = ξiξ , η0 = ηiη, ζ 0 = ζ iζ , i = 1,2,Λ ,8
有限元讲义弹性力学平面问题有限单元法(四结点四边形等参元,八结点曲线四边形等参元,问题补充)分析
2.6 四结点四边形单元(The four-node quadrilateral element)前面介绍了四结点的矩形单元其位移函数:xy y x U 4321αααα+++=xy y x V8765αααα+++=为双线性函数，应力，应变在单元内呈线性变化，比常应力三角形单元精度高。
但它对边界要求严格。
本节介绍的四结点四边形等参元，它不但具有较高的精度，而且其网格划分也不受边界的影响。
对任意四边形单元（图见下面）若仍直接采用前面矩形单元的位移函数，在边界上它便不再是线性的（因边界不与x,y 轴一致），这样会使得相邻两单元在公共边界上的位移可能会出现不连续现象（非协调元），而使收敛性受到影响。
可以验证，利用坐标变换就能解决这个问题，即可以通过坐标变换将整体坐标中的四边形（图a ）变换成在局部坐标系中与四边形方向无关的边长为2的正方形。
正方形四个结点i,j,m,p 按反时钟顺序对应四边形的四个结点i j m p 。
正方形的 1-=η 和 1=η 二条边界，分别对应四边形的i ，j 边界和p,m 边界；ξ=-1和ξ=+1分别对应四边形的i ，p 边界和j ，m 边界。
如果用二组直线等分四边形的四个边界线段，使四边形绘成一个非正交网格，那么该非正交网格在正方形上对应着一个等距离的规则网格（见图a, b ）。
当然, 局部坐标上的A 点与整体坐标的A 点对应。
一、四结点四边形等参单元的形函数及坐标变换由于可以将整体坐标下的四边形单元变换成局部坐标下的正方形单元，对于这种正方形单元，自然仍取形函数为： ξηαηαξαα2321+++=U ξηαηαξαα8765+++=V引入边界条件，即可得位移函数：∑=ijmpi i U N Ui ijmpi V N V ∑==写成矩阵形式：{}{}[]{}ee p i p i ed N d N N N N V U f =⎥⎦⎤⎢⎣⎡=⎭⎬⎫⎩⎨⎧=000 式中形函数: ()()()ηηξξηξi i i N ++=1141, ()p m j i ,,, 按照等参元的定义,我们将坐标变换式亦取为: p p m m j j i i i ijmpi x N x N x N x N x N x +++==∑p p m m j j i i i ijmpi y N y N y N y N y N y +++==∑ ()162-- 式中形函数N 与位移函数中的完全一致。
有限元基本知识归纳
有限元知识点归纳1.、有限元解的特点、原因？答：有限元解一般偏小，即位移解下限性原因：单元原是连续体的一部分，具有无限多个自由度。
在假定了单元的位移函数后，自由度限制为只有以节点位移表示的有限自由度，即位移函数对单元的变形进行了约束和限制，使单元的刚度较实际连续体加强了，因此，连续体的整体刚度随之增加，离散后的刚度较实际的刚度K为大，因此求得的位移近似解总体上将小于精确解。
2、形函数收敛准则（写出某种单元的形函数，并讨论收敛性）P49(1)在节点i处N i=1,其它节点N i=0;(2)在单元之间，必须使由其定义的未知量连续；(3)应包含完全一次多项式；(4)应满足∑Ni=1以上条件是使单元满足收敛条件所必须得。
可以推证，由满足以上条件的形函数所建单元是完备协调的单元，所以一定是收敛的。
4、等参元的概念、特点、用时注意什么？（王勖成P131）答：等参元—为了将局部坐标中几何形状规则的单元转换成总体（笛卡尔）坐标中的几何形状扭曲的单元，以满足对一般形状求解域进行离散化的需要，必须建立一个坐标变换。
即：为建立上述的变换，最方便的方法是将上式表示成插值函数的形式，即：其中m是用以进行坐标变换的单元节点数，xi,yi,zi是这些结点在总体（笛卡尔）坐标内的坐标值，Ni’称为形状函数，实际上它也是局部坐标表示的插值函数。
称前者为母单元，后者为子单元。
还可以看到坐标变换关系式和函数插值表示式：在形式上是相同的。
如果坐标变换和函数插值采用相同的结点，并且采用相同的插值函数，即m=n，Ni’=Ni,则称这种变换为等参变换。
5、单元离散？P42答：离散化既是将连续体用假想的线或面分割成有限个部分，各部分之间用有限个点相连。
每个部分称为一个单元，连接点称为结点。
对于平面问题，最简单、最常用的离散方式是将其分解成有限个三角形单元，单元之间在三角形顶点上相连。
这种单元称为常应变三角形单元。
常用的单元离散有三节点三角形单元、六节点三角形单元、四节点四边形单元、八节点四边形单元以及等参元。
有限元法基础5等参元与数值积分
有限元法基础
5. 3 弹性力学中等参单元的一般列式方法 有限元方程为 单元刚度矩阵为
34
有限元法基础
5. 3 弹性力学中等参单元的一般列式方法
1）母单元为
自然坐标系列
坐标变换
位移插值
Jacobi矩阵
应变的计算
求B时需建立
35
有限元法基础
5. 3 弹性力学中等参单元的一般列式方法 单元矩阵计算时
等参变换单元矩阵的变化：
等参变换
单元矩阵的变化：B、K、dΩ换的概念 由于插值函数使用自然坐标，涉及到求导和积分的变 换，如B矩阵的偏微分计算，K矩阵的积分计算。
13
有限元法基础
5.1等参变换的概念 1）导数之间的变换
由复合函数求导规则有
写成矩阵形式
J 称为Jacobi 矩阵
30
有限元法基础
5.2 等参变换的条件与收敛性 只要
Ni 满足形函数性质，完备性就得到满足， 插值函数能够反映刚体位移和常应变。
31
有限元法基础
5.2 等参变换的条件与收敛性
协调性 单元间边界上的位移场：
➢具有相同的节点和相同的节点数 ➢插值函数相同，有连续的位移场 ➢插值函数满足
32
有限元法基础
38
有限元法基础
5. 3 弹性力学中等参单元的一般列式方法 例：无限元 1）一维问题：2节点单元
通常u2是已知的。 39
有限元法基础
5. 3 弹性力学中等参单元的一般列式方法 例：无限元 2）二维问题：4节点单元
40
有限元法基础
5. 3 弹性力学中等参单元的一般列式方法 坐标变换
反映了1-2边的变化率。 位移插值函数依然与传统单元一样。 通常节点2和节点3的量是已知的。
有限元讲义等参单元
第五页，共24页
§6.2 平面四节点等参单元
• 该局部坐标系使得在x-y平面上的任意四边形与ξ-η平面上的正方形之间形成 了1-1对应的映射。正方形的4个顶点对应任意四边形单元的四个节点； 4条边 对应任意四边形单元的4条边；正方形内任一点p(ξ,η)对应于任意四边形内一 点p(x,y)。
• 上述映射是利用母单元描述实际单元力学特性的桥梁。由于该坐标变换式中 采用了与位移插值相同的节点和参数（插值函数），因此称为等参变换。而 所有采用等参变换的单元称为等参单元。等参单元是一个单元家族，目前在 通用程序中广泛采用。
4
x N i x i
i1
4
y N i y i
i1
第九页，共24页
函数对单元的节点坐标和节点位移在单元上进行插值。这种单元称为等参单元。
• 等参单元的提出对于有限元法在工程实践中的应用具有重要意义。
第三页，共24页
§6.2 平面四节点等参单元
1、局部坐标系与位移模式
下图为一个4节点任意四边形单元，单元有8个自由度。将矩形单元 放松为4节点任意四边形单元将带来许多好处。
4) 高阶等参元精度高，描述复杂边界和形状的能力强，所需单元少，在结构应 力分析中应用最广泛。
第二十四页，共24页
• 三维等参单元的力学分析和数学分析原理和方法与平面等参元相同。
第二十二页，共24页
§6.5 等参变换的条件 等参变换的条件
• 等参变换要保证母单元与实际单元之间形成1-1对应的映射， 数学上的条件是变换的Jacobi行列式大于零。要保证这个条 件，单元的几何必须满足一定要求，主要是： ①单元形状不能过度畸变；
平面有限元分析-等参单元
等参元变换的条件为 J ≠ 0 ，因此在有限元网格划分时，要 特别注意这一点。
等参单元等效节点力（4节点）
（1）集中力引起的单元节点载荷
单元内某点受到集中载荷P=[Px Py]T，移置到单元节 点上的等效节点力为：
j
同理 得
dη = ∂x dηi + ∂y dη j
∂η
∂η
∂x dξ
dA =
dξ × dη
=
∂ξ
∂x
dη
∂µ
∂y dξ ∂x
∂ξ
∂y
dη
=
∂ξ
∂x
∂η ∂µ
∂y
∂ξ
∂y
dξdη
∂η
等参单元刚度（4节点）
因为
∂x ∂y
J
=
∂ξ
∂x
∂ξ
∂y
∂η ∂η
雅可比行列式 Jacobi
曲边面积元dA:
dA = J dξdη
8 平面问题有限元分析 等参单元
8.1等参曹单国元华刚度（4节点） 8.2等参单元等效节点力（4节点） 8.3矩形单元（8节点） 8.4等参单元（8节点） 8.5高斯积分法
等参单元刚度（4节点）
4
4
=u ∑= Ni (ξ ,η )ui ,v ∑ Ni (ξ ,η )vi
=i 1 =i 1
4
4
=x ∑= Ni (ξ ,η ) xi , y ∑ Ni (ξ ,η ) yi
i =1
4
y = ∑ Ni (ξ ,η ) yi
i =1
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
节点条件： ui = u (ξi ,ηi ) vi = v (ξi ,ηi )
(ξ1 ,η1 ) = (−1, −1) (ξ 2 ,η 2 ) = (1, −1)
(ξ3 ,η3 ) = (1,1) (ξ 4 ,η4 ) = (−1,1)
北京航空航天大学
u1 = α1 + α 2ξ1 + α 3η1 + α 4ξ1η1 u = α + α ξ + α η + α ξ η 2 1 2 2 3 2 4 2 2 u3 = α1 + α 2ξ3 + α 3η3 + α 4ξ3η3 u4 = α1 + α 2ξ 4 + α 3η4 + α 4ξ 4η4
1 N1 = 4 (1 − ξ )(1 − η ) N = 1 (1 + ξ )(1 − η ) 2 4 1 N 3 = 4 (1 + ξ )(1 + η ) N = 1 (1 − ξ )(1 + η ) 4 4
北京航空航天大学
x N1 = y 0
β1 1 1 1 1 y1 β 2 1 −1 1 1 − 1 y2 = β3 4 −1 − 1 1 1 y3 β4 1 − 1 1 − 1 y4
x = α1 + α 2ξ + α 3η + α 4ξη = N1 x1 + N 2 x2 + N 3 x3 + N 4 x4 y = β1 + β 2ξ + β3η + β 4ξη = N1 y1 + N 2 y2 + N 3 y3 + N 4 y4
x = x(ξ )
坐标插值函数：
x(ξ ) = a0 + a1ξ
a0 = ( xi + x j ) / 2 a1 = ( x j − xi ) / 2
节点条件： x(−1) = a0 − a1 = xi x(1) = a0 + a1 = x j
北京航空航天大学
x(ξ ) = ( xi + x j ) / 2 + ξ ( x j − xi ) / 2
北京航空航天大学
当点的运动轨迹已知时，通常采用自然法 确定点的运动规律、速度、加速度。 在自然坐标系中表示质点速度，是非常简 单的，因为无论质点处在什么位置上速度 都只有切向分量，而没有法向分量。
北京航空航天大学
新途径：建立局部自然坐标系进行单元分析
xi
xj
局部自然坐标和整体直角坐标可以建立一种映射关系
N(ξ ) = [(1 − ξ ) / 2
N(ξ ) = [(1 − ξ ) / 2
(1 + ξ ) / 2]
(1 + ξ ) / 2]
北京航空航天大学
du ( x(ξ )) d (N(ξ )q e ) dN(ξ ) e ε ( x(ξ )) = = = q = B(ξ )q e dx dx dx
Ωe Sp
(1 + ξ ) / 2]
Sp
P (1) = ∫ NT (ξ )bx dV + ∫ e NT (ξ ) px dA = ∫ e NT (−1) px dA 1 1 R = ∫ e px dA = R1 = 1 S p 0 0 0 P (2) = ∫ NT (ξ )bx dV + ∫ e NT (ξ ) px dA = ∫ e NT (1) px dA
u ( x(ξ )) = N(ξ )q e
N(ξ ) = [(1 − ξ ) / 2 (1 + ξ ) / 2]
北京航空航天大学
观察： 单元的几何坐标与位移用同样的节点和相同的 形状函数通过插值的方式表示。 形状函数是用自然坐标给出的，表达式很简单
x(ξ ) = N(ξ )xe
u ( x(ξ )) = N(ξ )q e
北京航空航天大学
平面问题四边形等参单元的推导
坐标映射
( x4 , y4 ) ( x3 , y3 )
P ( x, y )
( x1 , y1 )
P (ξ ,η )
( x2 , y2 )
整体直角坐标
（一般四边形）
P ( x, y )
单元局部自然坐标
（规格化的矩形）
映射
P (ξ ,η )
北京航空航天大学
P ( x, y )
e 1 e
T
e
e
1 − 1 −1 1
北京航空航天大学
单元外力功 W = q P e
e eT
等效节点力
P e = ∫ NT (ξ )bx dV + ∫ e NT (ξ ) px dA
Ωe Sp
N(ξ ) = [(1 − ξ ) / 2
对于本例
自然坐标系下的分 析结果与整体直角 坐标系下的分析结 果完全相同。
a0 = (ui + u j ) / 2 a1 = (u j − ui ) / 2
ui (1 + ξ ) / 2] u j
u ( x(ξ )) = (ui + u j ) / 2 + ξ (u j − ui ) / 2 = [(1 − ξ ) / 2
e
单元刚度矩阵
1 1 e U = ∫ E [B(ξ )q e ] B(ξ )q e Ae (l e / 2)d ξ 2 −1 T 1 eT 1 e e U = q [ ∫ (l / 2)B (ξ ) E e Ae B(ξ )dξ ]q e −1 2 1 eT e e e U = q Kq 2
e
E e Ae K = ∫ (l / 2)B (ξ ) E A B(ξ )d ξ = e −1 l
节点条件： xi = x(ξi ,ηi ) yi = y (ξi ,ηi )
x1 = α1 + α 2ξ1 + α 3η1 + α 4ξ1η1 x = α + α ξ + α η + α ξ η 2 1 2 2 3 2 4 2 2 x3 = α1 + α 2ξ3 + α 3η3 + α 4ξ3η3 x4 = α1 + α 2ξ 4 + α 3η4 + α 4ξ 4η4
第4讲 等参单元和数值积分
金朝海 jch666@
北京航空航天大学
实际问题常常需要使用一些几何形状不太规整 的单元来逼近原问题。直接研究这些不规整单 元的表达式比较困难（在整体坐标系下构造位 移插值函数，则计算形状函数矩阵、单元刚度 矩阵及等效节点载荷列阵时十分冗繁）。事实 上，形状不规整的单元和形状规整的单元（矩 形单元、正六面体单元）可以建立一种映射关 系，使得物理坐标系中的整体坐标和自然坐标 系中的局部坐标一一对应。 等参单元的提出为有限元法成为现代工程实际 等参单元 领域最有效的数值分析方法迈出了极为重要的 一步。
x(ξ ) = [(1 − ξ ) / 2 xi (1 + ξ ) / 2] x j
x(ξ ) = N(ξ )x e
N(ξ ) = [(1 − ξ ) / 2
局部坐标到物理坐标的变换
(1 + ξ ) / 2]
单元内坐标由节点坐标插值表示
北京航空航天大学
单元位移函数： u ( x(ξ )) = a0 + a1ξ 节点条件： u ( x(−1)) = a − a = u 0 1 i u ( x(1)) = a0 + a1 = u j
0 N1
N2 0
0 N2
N3 0
0 N3
N4 0
x1 y 1 x2 0 y2 N 4 x3 y3 x4 y 4
x(ξ ,η ) = N(ξ ,η )x e
北京航空航天大学
位移函数
( x3 , y3 )
( x4 , y4 )
(−1,1)
(1,1)
( x2 , y2 ) ( x1 , y1 )
(−1, −1)
(1, −1)
u ( x(ξ ,η ), y (ξ ,η )) = u (ξ ,η ) = α1 + α 2ξ + α 3η + α 4ξη v( x(ξ ,η ), y (ξ ,η )) = v(ξ ,η ) = β1 + β 2ξ + β3η + β 4ξη
u1 = α1 − α 2 − α 3 + α 4 u = α + α − α − α 2 1 2 3 4 u3 = α1 + α 2 + α 3 + α 4 u4 = α1 − α 2 + α 3 − α 4
α1 1 1 1 1 u1 α 2 1 −1 1 1 − 1 u2 = α 3 4 −1 − 1 1 1 u3 α 4 1 − 1 1 − 1 u4
北京航空航天大学
关于坐标系
直角坐标系( 直角坐标系 x , y , z) 极坐标(r, 极坐标 ϕ) ，2维 维 球坐标系(r,θ, ϕ) 球坐标系 柱坐标系 (ρ, ϕ, z)
自然坐标系
北京航空航天大学
自然坐标系: 自然坐标系
选轨迹上任一点O为原点 选轨迹上任一点 为原点 用轨迹长度 描写质点位置 轨迹长度S
x1 = α1 − α 2 − α 3 + α 4 x = α + α −α −α 2 1 2 3 4 x3 = α1 + α 2 + α 3 + α 4 x4 = α1 − α 2 + α 3 − α 4