三角形高阶

单元位移模式和形函数矩阵N e 如下：
单元几何矩阵Be如下：
单元应力矩阵Se如下： 本构矩阵 D 对平面应力和平面应变略有不同：
E; 平面应力
E1
=
{ 1
E − μ2
；
平面应变
������；平面应力 ������1 = { ������ ； 平面应变
1 − ������ 单位刚度矩阵Ke如下： Ke = ∫ BeT DBe������V
… …
N4 0
0 N4
]
δe
=
Neδe
图 6：5 节点单元
N1 = L1 N2 = L2(1 − 3L3) N3 = 3(1 − 3L2) N4 = 2L2L3 在有限元分析中，经常会碰到一些比较 特殊的 几何 区域， 需要较 高的计 算精 度， 例如凹角、台阶和孔洞等突变区域；有多种材料连接的区域；边界条件比较复杂 的区域；需要特别关注的区域。这些区域有可能只在一个方向或很小的范围对精度 要求高，这时，这一类 4 节点三角形单元和 5 节点单元就可以很好的满足这一要 求。 3 评价基于 ANSYS: Ansys 中选择三角形plane182 单元和三角形plane183 单元，分别代表三角 形三节点单元和六节点单元，用来分别对简支梁的受力情况进行计算:
开卷试题二
一 类 三 角形高次元的 设想、实现与评价
1、 设 想 :
作为一种较为普及的有限元单元，三角形单元非常适合与求解二维问题。并且能够
处理二维边界问题，此外它的形函数相对便于构建。
线性三角形单元在有限元分析中有着计算简单,计算量小,几何适应性好,单元易分割等
优点。特别是当对包含锐角的、几何形状复杂的二维模型进行网格划分时，往往会采用三
这些单元有着下面几个共同特点：有中节点的边是高阶的，没有中节点的边是
线性的，若所有中节点都不存在则该单元退化成为线性单元。这样的节点分布可以
满足部分边的精度要求又可以减少单元整体的自由度。
以图 5 中的 4 节点三角形单元为 Nhomakorabea，其位移差值和形函数为：
ue
=
{UV}
=
[N01
0 N1
N2 0 0 N2
5%
Ｅ１８３ ８０
5.29944 5.8753 5.625 6%
4%
表２
从上表中可以看出，ＰＬＡＮＥ１８２单元（线性三角形单元）当单元总数增大使得拉
应力和压应力的计算精度逐渐增大，在单元总数均大于４０的情况下，二次六节点计算出来
的误差在５％左右，而线性三角形单元所计算的误差超过了６０％，当其单元总数达到８０
，
该单元应变、应力随坐标完全呈线性变化，属于高精度单元。 单元位移模式和形函数可以通过三角形单元的面积坐标来描述。
图 2 面积坐标
如图 2 所示，所示，三角形任意一点 P 的面积坐标为(������������，������������，������������)， 定义为������������ =
图 12 长度方向 20 份高度方向 2 份
图 13 长度为 40 份高度为 2 份 当采用 PLANE183 单元时计算结果如下：
图 14 单元密度为 15
图 17 单元密度为 7.5 现在回到模型上（见图 9）： 应为按照平面问题处理，所以在求解集中力点处的挠度可以使用位移法。
UY = F4������×2������×12×2×4������×23，其中I = bh3
采用线性三角形单元PLANE182的压力云图如下：
图 18 单元总数为 40 的压力分布云图 图 19 单元总数为 80 压力分布云图
图 20 单元总数为 160 的压力云图 图 21 单元总数为 800 时的压力云图
采用高阶三角形单元 PLANE182 的压力云图如下： 图 22 单元总数为 40 的压力分布云图 图 23 单元总数为 80 的压力分布云图
Ve
将单元几何矩阵Be和本构矩阵 D 代入上式，可得单元刚度矩阵Ke的表达式如下：
通过上述理论推导,我们得到了三角形二阶单元的位移模式,并分析了其完备性 和协调性给出了形函数,几何矩阵,应力矩 阵和单 元刚度 的具 体表达 式。 2.2 三角形 10 节点(三角形的每条边上插入 2 个两个节点,另加一个形心) 其位移模式为：
位移模式为：
图 1 三角形（2 阶）6 节点单元
很显然该模式为完全二次多项式，能够较好地模拟刚体位移和常应变，因此其满足完备性 条件。此外由于该位移模式决定了单元边界上位移呈二次抛物线分布，相邻单元公共边界 上有三个公共节点，保证了相邻单元在边界上位移的连续性，因而是协调元，单元满足收 敛条件。其应变可以写为：
但是，高次三角形单元仍有一些问题没有解决，如高次三角形形函数的普遍性还不理 想，不能用原有的线性三角形单元的形函数的方法来推导，和一元函数的插值方法不相 同，故有线性三角形到高次三角形单元的转化是一个具有挑战的问题。因此这里主要讨论 和构造二次三角形单元。
2、 三 角 形二阶单 元和三阶单元 2.1 三角形 6 节点（在三边中点处分别插入一个节点）
33%
800
0.004996
0.0050242
1%
PLANE183
40
80
0.005245 0.005307
0.0050242
4%
0.0050242
6%
表1 同样是单元数为40，线性三角形单元的误差达到了50%，而PLANE183 的二次六节点插值三 角形单元的 误差只有4%，可以看出PLANE182 的三角形单元即线性三角形单元的精度要差 很多，而PLANE182 要想将精度控制在5%之内，要将单元数量增加至800.而PLANE183 即二 次六节点二次插值三角形单元计算出位移的相对误差随着网格数量的增加 没有明显变化， 而且有增大的趋势，所以6 节点三角形单元不需要将单元数量设置大来获得高精度。
０时期拉应力的误差仍然在１５％，压应力的误差控制到了５％左右，可见二次六节点三角
形单元的计算精度要明显的高于线性三角形单元。
图 4 三角形 6 阶（10 节点）单元
其面积坐标的形式函数为：
ue
=
{UV}
=
[N01
0 N2 N1 0
0 N2
… …
N10 0
0 N10
]
δe
=
Neδe
1 Ni = 2 (3Li − 1)(3Li − 2)Li, i = 1,2,3
9 N4 = 2 L1 L2(3L1 − 1)
9 N5 = 2 L1 L2(3L2 − 1)
������������ = (2������������ − 1)������������ ，������ = 1，2，3 ������4 = 4������1������2 ������5 = 4������2������3 ������6 = 4������3������1
单元节点位移向两以及单元节点载荷向量定义如下：
10 节点三角形的位移模式是完全三次多项式，因此该模式 满足完备性条件。每条边上都是三次函数。同时每条边都有 4 个节点，可以确定一个三次函数，因此该单元的位移模式 也满足协调性条件。
三角形三阶 10 节点单元的位移模式和形函数和 6 节点 三角形不同，但几何矩阵，应力矩阵和单元刚度矩阵的推导 方法与之前一样。
EI
12
解得： UY=0.0050223mm
纵向位移方面，采用 PLANE182 和 PLANE183 单元的对比如下
单元类型
单元总数
竖直方向位移 理论解
相对误差
PLANE182 40 80
0.002481 0.003131
0.0050242
51%
0.0050242
38%
160
0.003362
0.0050242
图 7 建立模型
首先建立高度为20mm长度为 150mm的 梁模型 ，厚度 设置为 1，杨式 模量设 置 为2.1e5MPa。泊松比设置 为0. 3， 该梁两 端简支 中间受 垂直向 下的力 10N。
然后对模型划分网格，长度方向上的单 元先后 设置 为15， 7. 5， 3. 75将高度 方向 上的单元设置为10以下是网格划分情况：
图 8 长度方向分为 10 份高度方向分为 2 份
图 9 长度方向 20 份高度方向 2 份
图 10 长度方向 40 份高度方向 20 份 添加约束
在座节点添加 XY 方向的位移约束，由节点只添加 Y 方向上的约束，并在重点加载 10N 的集中力。 当采用 PLANE182 单元时计算结果如下：
图 11 长度方向 10 份高度 2 份
ａ）
ａ）
ＰＬＡＮ ４０
1.56219 1.59952 5.625 72%
72%
Ｅ１８２ ８０
2.08817 2.15279 5.625 63%
62%
１６０ 2.32886 2.37878 5.625 59%
58%
８００ 4.78208 5.28041 5.625 15%
6%
ＰＬＡＮ ４０
5.2703 5.32207 5.625 6%
������������������，������������
=
������������，������
������
������
=
������������。面积坐标与直角坐标的转换关 系如下 ：
������
图 3 节点编号 三角形二阶单元的节点顺序如图3 所示，节点1,2,3 逆时针排列，节点4,5,6 分别位于 12 边，23 边，31 边的中点。利用面积坐标，形函数可写为：
角形单元，而且其简单的拓扑性质使得研制网格生成器变得相对容易。但是,由于其位移差
值函数是线性的,这就不可避免的产生一些缺点。线性三角形单元的差值函数为：