有限元法剪切锁闭现象的研究论文答辩

图3.1 具有剪切变形影响的梁变形
9
设剪切变形为 ，则
（3.1） 设单元的总挠度函数 和截面转角函数 的单元插值模式为（线性插值 ）： ( x) N1 ( x)1 N2 ( x)2 （3.2） ( x) N1 ( x)1 N2 ( x)2 （3.3）
( 其中 ( x)，x) 分别为单元节点1和2的挠则
曲率为：
Bb ( x ) q e
1 1 B b 0 0 l l
B b 为曲率几何矩阵
(3.11)
11
将上式代入得到单元势能泛函：
e 1 eT b e 1 eT s e q K q q K q P eT q e 2 2


1 eT e e q K q P eT q e 2
16
梁为150mm长， 2.5mm宽， 5mm高；一端固定；自由端承受5N的 末端荷载。材料的杨氏模量E为70GPa，泊松比为0.0。采用梁的理论 ，在给定载荷P作用下，梁末端的挠度为 (3.22) 其中I = bd3/12， 是长度，b是宽度，d是梁的高度。P = 5N时末端挠 度是3.09mm。
4
2.1 求解精度的估计
在这里我们考察平面单元的求解精度与收敛速度，单元的位移场 可 以展开为以下级数： u u (2.1) u ui ( )i x ( )i y
x
y
对于满足完备性（completeness）和协调性（compatibility）要求的 协调元，当 时，有限元分析的结果是单调收敛的如图(2.1)
8
3.1 Timoshenko梁单元
当梁不是细长梁时，梁变形后的横截面垂直于中性层的假设不再成立（ Kirchhoff假设）这时需要考虑梁的剪切变形，考虑剪切变形的几何描 述如图3.1所示。
变形前的截面 变形后的截 面 变形前的截面 中性层 垂直于中性 层的截面
γ
中性层
β dω/dx (a)变形前 (b)变形前
l 1 2 l2 l l GA GA 2 3 K s ( B s )T B s dx 0 l l 1 2 l l2 6 2
0 0 0 1 0 0 0 1
1 l 2 l 2 l2 6 0 l2 3
(4.4)
24

图4.1
25

将(4.6)代入（4.5）的到
Biblioteka Baidu
K1s qe 0
(4.7)
26
27
图４.2 使用缩减积分的二维单元中的积分点
28
单元 CPS4R CPS8R C3D8R
表4.1 使用缩减积分单元的梁挠度的归一化结果
网格尺寸 (高度 x 长度)
1x6 20.3* 1.000 70.1* 2x12 1.308 1.000 1.323 4x12 1.051 1.000 1.063 8x24 1.012 1.000 1.015
=F, x
F F , y x y
解：
(4.15)
y 0 y

xy
x 0 x
3
有限元分析的误差
结构体由于本身存在有自然的连接关系即自然节点,所以它们的离散化 叫自然离散 人为的在连续体内部与边界上划分节点,以单元连续的形式来逼近原来 复杂的几何形状,这种过程叫逼近性离散(approximated discretization) 控制误差的h方法（h-version ,h-method）和p方法(p-version ,pmethod)。

图4.4 变形梯度的变化
32
图4.5使用位移场增强而非变形梯度增强非协调模式单元间潜在的运动非 协调性，ABAQUS/Standard使用后一表达式为其非协调模式单元

图4.6 非协调模式单元的扭曲网格
33
图4.7 平行和交错扭曲对非协调元的影响
34
4.2 厚板单元解析试函数法
在有限元的离散法中吸取解析法的优点并加以利用，很早就有文献论 述。利用基本解析解作为试函数而构造出的广义协调元
31
4.1.3 非协调元
非协调元是为克服完全积分的，一阶线性单元的剪切闭锁问题的一种 尝试。由于剪切闭锁是单元的位移场不能模拟与弯曲相关的变形而引 起的，所以在一阶单元中引入一个增强单元变形梯度的附加自由度。 这种对变形梯度的增强可以允许变形梯度在一阶单元的单元域上有一 个线性变化，如图4.4(a) 所示，而标准的单元数学公式在单元中只能 得到一个常数变形梯度，如图4.4（b）所示，这导致与剪力闭锁相关 的非零剪切应力。 在弯曲问题中，非协调元可能产生与二次单元相当精度的结果，但计 算成本却明显的下降。然而，它们对单元的扭曲很是敏感。
(4.11)

,
x l
y l
(4.12)
36
代入得到:
l 如果 h 0
D h2 2 Cl 5(1 )l 2
2 2 F 2 2 F
(4.13)
(4.14)
0 ,则厚板理论的解析解退化为薄板理论的解析
方法II (, ) 从假设 入手 方法Ⅲ ( , )从假设 入手
2.块体问题 将在缩减积分详细介绍
22
4.1 缩减积分法、选择性缩减积分法、非协 调元

0 0 l EI 0 1 K b ( B b )T EI B bdx 0 l 0 0 0 1

d dx
N1 ( x ) (1
x x ), N 2 ( x ) l l
（3.4） （3.5）
由（3.1）可得到剪应变为：

d( x) 1 1 x x ( x ) ( )1 ( )2 (1 ) 1 ( ) 2 dx l l l l
(4.8)
35
设F(x,y)是双调和函数，满足双调和方程：
其中： 则下列为上式微分方程的解：
F
22 F 0
2 2 2 2 2 x y
D 2 F F F , x , y C x y
(4.9)
(4.10)
设 l 是一个特征长度。引入无量纲坐标 ， ：
论文题目：
有限元法剪切锁闭现象的研究
专 业：工程力学 答辩人：徐 鹏 导 师：李欣宇
武汉科技大学理学院 2007年6月16日
答辩提纲
课题背景及任务 有限元分析的误差 有限单元特征及剪切闭锁产生 剪切闭锁的对策 结论及下一步工作 致谢
2
课题背景及任务
当梁不是细长梁时，梁变形后的横截面垂直于中面的假设不再成 （Kirchhoff假设）立，需考虑剪切变形的影响，于是Timoshenko提 出了剪切变形的梁理论，如果对扰度函数与转角函数进行独立插值， 且考虑剪切变形的影响，构造出的单元为Timoshenko梁单元 （Timoshenko beam element）。此单元在退化为细长梁单元时，会 导致剪切闭锁现象。经典薄板理论以Kirchhof假设为基础，忽略了横 向剪切的影响，故计算结果与实验值相比，总是低估了挠度，高估了 自然频率。由于横向剪切变形对厚板及复合材料的板单元影响较大， 故忽略剪切变形的薄板理论已不再适用，于是Mindlin等人提出了板的 剪切变形理论,以此构造出的为Mindlin板单元（Mindlin plate element）。与Timoshenko梁单元类似，当板很薄时，会发生剪切闭 锁与零能模式。对于块体，由于形函数的阶数过低，而单元受到复杂 应力作用，如产生弯曲扭转等，导致形函数不能真实表达。本文讨论 剪切闭锁产生的原因，与对策。
6
2.2提高精度的ｈ方法和ｐ方法
（１）ｈ方法（h-version）:不改变各单元 上的形状函数，只通过渐渐加密有限元的 网格使结果趋于准确解 （２）ｐ方法（p-version）:同ｈ方法相反 ，它是不改变单元网格，而采用较高阶的 多项式来进行插值。
7
3.有限元的单元特征及剪切闭锁的产生
Timoshenko梁单元的剪切闭锁现象 Mindlin板单元的剪切闭锁现象 块体单元的闭锁现象
(3.12)
其中单元的刚度阵由剪切和弯曲变形刚度阵组成

Ke Ks Kb
(3.13)
12
13
3.2Mindlin板单元
G 2 A r 2 dx
14
15
3.3 块体单元
以如图3.2悬臂梁的静力分析为例来说明块体的剪切闭锁现象
图3.2 在自由端受到点载荷P的悬臂梁
C3D20R
0.999**
1.000
1.000
1.000
注：* 没有刚度抵抗所加载荷, ** 在宽度方面使用了两个单元。
29

图4.3 承受弯矩M的使用缩减积分的线性单元的变形
30
4.1.2 选择性缩减积分
实质上同缩减积分没有区别，只是在处理实际问题时加以改善了，对 于可能会发生锁闭的局部采用缩减积分，减少计算量。
tip
Pl 3 , 3EI
17
如图3.4所示，使用几种不同的单元网格划分对悬臂梁问题进行模拟 。模拟既采用线性又采用二次的完全积分单元，并说明了单元阶数( 一阶与二阶)和网格密度对结果精度的影响。

图3.4 用于悬臂梁模拟的网格划分
18
由上知道四节点的线性单元（CPS4）与8节点的六面体单元预测的挠 度值基本上是不可以用的。
在厚板平衡微分方程中，如果将内力改用位移表示，就得到厚板位移 法的基本微分方程。考虑载荷为零时的齐次问题，则有：

2 x 1 2 x 1 2 y D( 2 ) C( x ) 0 2 x 2 y 2 xy x 1 2 1 2 y 2 y x ) C( y) 0 D( 2 2 2 xy 2 x y y 2 2 C( 2 2 )0 x y x y
19

图3.6 8节点的六面体单元
20
4 剪切闭锁的对策
1.板（梁）问题
x 厚板理论挠度和转角是独立位移 y
x y 0 而剪应变 0 y x

21
方法I (, ) 从假设 入手
ui(1) ui O(h 2 ) 4 (2) 2 ui ui O((h / 2) )
ui(1) ui O(h s ) (2) ui ui O((h / 2) s )
(2.2)
（2.3）
(2.4)
可由此估计出准确解：
1 ui (4ui(2) ui ) 3
由有限个单元的试函数来逼近整体域的场函数所引起的误差，即离散 误差。
10
曲率为：

d ( x) 1 1 ( ) 1 ( ) 2 dx l l
(3.6)
将节点的位移列阵记为：
qe 1 1 2
2
T
(3.7) (3.8) (3.9) (3.10)
则（3.5）化为：
B s ( x ) qe
1 1 B b 0 0 l l
(4.1)

(4.2)
1 l 2
23
剪切应变能的表达式:
具体的有一点积分得到的
Ur
1 eT s e q K q 2
(4.3)

l 1 2 l2 l GA 2 4 Ks l l 1 2 l l2 4 2
1 l 2
1 l 2
l 2 2 l 4 0 l2 4
图2.1 有限元分析结果的收敛情况（纵坐标是总势能）
5
还可以就两次网格划分所计算的结果进行外推以估计结果的准确值如 第一次网格划分的结果是 ,然后进一步将各单元尺寸减半进行网格划 分,得到结果为 .假设该单元的收敛速度是 ,则其准确解可以按如下方 法估计:
具体对平面3节点三角形单元有s=2，上式可化为：