平面三角形单元有限元程序设计
平面单元的有限元法
u
1
5
3
2
y
2x
3
5
2
y
则单元刚体位移为
v
4
5
2
3
x
6
y
3
2
5
u
1
5
3
2
y
v
4
5
2
3
x
记为
u v
1 4
0 y 0x
显然，位移函数包含 了单元的刚体位移 （平动和转动）
u v
j j
um
vm
[I]是单位矩阵，
[N]称为形函数矩阵，
Ni只与单元节点坐标有关，称为 单元的形状函数
4-2 平面问题的常应变(三角形)单元
据弹性力学几何方程得单元的应变分量
u
x y
xy
x
4-1 有限单元法的计算步骤
弹性力学平面问题的有限单元法包括五个主要步骤： 1、所分析问题的数学建模 2、离散化 3、单元分
析 4、整体分析与求解 5、结果分析
图 3-1
4-2 平面问题的常应变(三角形)单元
有限单元法的基础是用所谓有限个单元的集合 体来代替原来的连续体，因而必须将连续体简化为 由有限个单元组成的离散体。对于平面问题，最简 单，因而最常用的单元是三角形单元。因平面问题 的变形主要为平面变形，故平面上所有的节点都可 视为平面铰，即每个节点有两个自由度。单元与单 元在节点处用铰相连，作用在连续体荷载也移置到 节点上，成为节点荷载。如节点位移或其某一分量 可以不计之处，就在该节点上安置一个铰支座或相 应的连杆支座。如图3-1
平面三角形单元有限元程序设计
平面三角形单元有限元程序设计P9 m 9 m一、题目如图1所示,一个厚度均匀的三角形薄板,在顶点作用沿板厚方向均匀分布的竖向载荷。
已知:P=150N/m,E=200GPa,=0、25,t=0、1m,忽略自重。
试计算薄板的位移及应力分布。
要求:1.编写有限元计算机程序,计算节点位移及单元应力。
(划分三角形单元,单元数不得少于30个);2.采用有限元软件分析该问题(有限元软件网格与程序设计网格必须一致),详细给出有限元软件每一步的操作过程,并将结果与程序计算结果进行对比(任选取三个点,对比位移值);3.提交程序编写过程的详细报告及计算机程序;4.所有同学参加答辩,并演示有限元计算程序。
有限元法中三节点三角形分析结构的步骤如下:1)整理原始数据,如材料性质、荷载条件、约束条件等,离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。
2)单元分析,建立单元刚度矩阵。
3)整体分析,建立总刚矩阵。
4)建立整体结构的等效节点荷载与总荷载矩阵5)边界条件处理。
6)解方程,求出节点位移。
7)求出各单元的单元应力。
8)计算结果整理。
一、程序设计网格划分如图,将薄板如图划分为6行,并建立坐标系,则刚度矩阵的集成建立与总刚度矩阵等维数的空矩阵,已变单元刚度矩阵的集成。
由单元分析已知节点、单元的排布规律,继而通过循环计算求得每个单元对应的节点序号。
通过循环逐个计算:(1)每个单元对应2种单元刚度矩阵中的哪一种; (2)该单元对应总刚度矩阵的那几行哪几列 (3)将该单元的单元刚度矩阵加入总刚度矩阵的对应行列循环又分为3层循环:(1)最外层:逐行计算(2)中间层:该行逐个计算(3)最里层:区分为第 奇/偶 数个计算XYPXYP节点编号单元编号单元刚度的集成:[][][][][][]''''''215656665656266256561661eZeeeZeZeeeekkkKkkkkkk+⋯++=⇓=⇒==⇒==⇒=⨯⨯⨯⨯⨯⨯M边界约束的处理:划0置1法适用:这种方法适用于边界节点位移分量为已知(含为0)的各种约束。
matlab有限元三角形单元编程
matlab有限元三角形单元编程
在MATLAB中进行有限元分析，可以使用其自带的有限元分析工具箱（FEATool）进行编程。
以下是一个简单的例子，演示如何使用三角形单元进行有限元分析：
1. 打开MATLAB，进入FEATool环境。
2. 创建新的有限元模型。
选择“Model”菜单下的“New Model”选项，进入“Model Builder”界面。
3. 在“Model Builder”界面中，选择“2D Triangle”单元类型，并在绘图区域中绘制出三角形网格。
4. 在“Model Builder”界面中，设置材料属性、边界条件和载荷等参数。
5. 运行有限元分析。
选择“Model”菜单下的“Solve”选项，进行有限元求解。
6. 查看结果。
选择“Model”菜单下的“Results”选项，可以查看位移、应力、应变等结果。
以上是一个简单的例子，演示了如何使用三角形单元进行有限元分析。
在实际应用中，还需要根据具体问题进行详细的建模和计算。
有限元2-弹性力学平面问题有限单元法(2.1三角形单元,2.2几个问题的讨论)分析
第2章弹性力学平面问题有限单元法2.1 三角形单元(triangular Element)三角形单元是有限元分析中的常见单元形式之一，它的优点是:①对边界形状的适应性较好，②单刚形式及其推导比较简单,故首先介绍之。
一、结点位移和结点力列阵设右图为从某一结构中取出的一典型三角形单元。
在平面应力问题中，单元的每个结点上有沿x、y两个方向的力和位移，单元的结点位移列阵规定为：相应结点力列阵为: (式2-1-1)二、单元位移函数和形状函数前已述及，有限单元法是一种近似方法，在单元分析中，首先要求假定（构造）一组在单元内有定义的位移函数作为近似计算的基础。
即以结点位移为已知量，假定一个能表示单元内部（包括边界）任意点位移变化规律的函数。
构造位移函数的方法是：以结点(i,j,m)为定点。
以位移(u i ,v i ,…u m v m )为定点上的函数值，利用普通的函数插值法构造出一个单元位移函数。
在平面应力问题中，有u,v两个方向的位移,若假定单元位移函数是线性的,则可表示成:(,)123u u x y x yααα==++546(,)v v x y x yααα==++(2-1-2)a式中的6个待定常数α1 ,…, α6 可由已知的6个结点位移分量(3个结点的坐标) {}⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧=⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=mjimeddddmjjivuvuvui{}iijjmXYX(2-1-1)YXYiejmmFF FF⎧⎫⎪⎪⎪⎪⎧⎫⎪⎪⎪⎪⎪⎪==⎨⎬⎨⎬⎪⎪⎪⎪⎩⎭⎪⎪⎪⎪⎪⎪⎩⎭确定。
将3个结点坐标(x i,y i ),(x j,y j ),(x m,y m )代入上式得如下两组线性方程:123i i i u x y ααα=++ 123j j j u x y ααα=++ (a)123m m m u x y ααα=++和546i i i v x y ααα=++ 546j j j v x y ααα=++ (b)546m m m v x y ααα=++利用线性代数中解方程组的克来姆法则,由(a)可解出待定常数1α 、2α 、3α :11A Aα=22A Aα=33A Aα=式中行列式:1i i i j j j m m m u x y A u x y u x y =2111i i j j m mu y A u y u y =3111i i j j m m x u A x u x u = 2111i i j j m mAx y A x y x y ==A 为△ijm 的面积,只要A 不为0,则可由上式解出:11()2m m i ij j a u a u a u A α=++ 21()2m m i ij j bu b u b u A α=++ （C ）31()2m mi i j j c u c u c u A α=++式中:m m i j j a x y x y =- m m j i i a x y x y =- m i j j i a x y x y =-m i j b y y =- m j i b y y =- m i j b y y =- （d ）m i j c x x =- m j i c x x =- m j i c x x =-为了书写方便，可将上式记为: m m i j ia x y x y =-m ij by y =- (,,)i j m u u u u ruu u u r m i jc x x =-(,,)i j m u u u u ru u u u r表示按顺序调换下标,即代表采用i,j,m 作轮换的方式便可得到(d)式。
第六章三角形单元的有限元法程序设计
从第j列的主对角线元素起到该列上方第一个非零 元素为止，所含元素的个数称为第j列的列高，记为hj ； 如果把第j列上方第1个非零元素的行号记为mj，则第j 列的列高为 hj = j - mj + 1 其实，hj就是第j行的左带宽，因而必有 UBW= max(hj)
j=1,2, …,N
利用节点位移信息数组 ID （去约束后节点位移自 由度编码），可容易地确定刚度矩阵 [K] 任何一列的 列高。
{ } [S ]{ }
例1：对角受压的正方形薄板，载 荷沿厚度均匀分布，为2N/m。由 于对称性，取1/4部分作为计算对
2N/m y x
象，试用有限元程序进行计算。
2N/m 2m 2m
h 1.0m, E 1.0, 0.0
例2：简支梁，梁高3m,跨度18m，厚度1m，承 受均布荷载10N/m2。已知 E 2 1010 N / m2 , 0.167
y为挤压应力，属次要应 力。材料力学不考虑 , FEM与弹性力学误差较大。
6.2 提高计算精度的方法
（1） 计算结果的整理 计算结果包括位移和应力两个方面。在位移方 面，一般无须进行整理工作。应力结果则需要 整理。通常认为计算出的应力是三角形单元形 心处的应力。而相邻单元之间的应力存在突变， 甚至正、负符号都不相同。为了由计算结果推 算出结构内某一点的接接实际的应力，必须通 过某种平均计算。通常可采用两单元平均法或 绕结点平均法。
平均法整理单元应力
两单元平均法：把两个相邻单元中的常应力加 以平均，用来表示公共边界中点处的应力。 绕结点平均法：把环绕某一结点的各单元常应 力加以平均，用以表示该结点的应力。在内结 点效果较好，而在边界结点可能很差，一般改 为应由内结点的应力外推计算出来。 （2）网格的细分 通过网格的细分，使每个单元的面积缩小，那 么尽管每个单元是应变、常应力单元，仍可较 好地反映结构中的应力变化，使得到的解答收 敛于问题的精确解。
有限元方法课件 第四章 平面三角形单元
第四章 平面三角形单元
§4–1 有限元法的基本思想 §4–2 三角形常应变单元 §4–3 形函数的性质 §4–4 刚度矩阵 §4–5 等效节点力载荷列阵 §4–6 有限元分析的实施步骤 §4–7 计算实例
§4-1 有限元法的基本思想
一、有限元法的基本思想 假想的把一连续体分割成数目有限的小体（单元），
vi (Vi )
i ui (Ui )
m
um (Um )
o
x
图4-2 平面三角形单元
将 (d) 式代入 (b) 式的第一式，经整理后得到
u 1 2ai源自bi x ci yuiaj
bjx cj y
uj
am bm x cm yum
(e)
其中 同理可得
ai
xj xm
yj ym
x j ym xm y j
这样，位移模式 (e) 和 (f) 就可以写为
u Ni ui N j u j N mum v Nivi N jv j Nmvm
(4-11)
也可写成矩阵形式
f
u v
Ni I
NjI
NmI e N e
(4-12)
式中 I是二阶单位矩阵；Ni 、Nj 、Nm 是坐标的函数， 它们反映了单元的位移状态，所以一般称之为形状函数，简 称形函数。矩阵 [N] 叫做形函数矩阵。三节点三角形单元的 形函数是坐标的线性函数。单元中任一条直线发生位移后仍 为一条直线，即只要两单元在公共节点处保持位移相等。则 公共边线变形后仍为密合。
f N e
（4-1）
f ——单元内任一点的位移列阵； e——单元的结点位移列阵；
N ——单元的形函数矩阵；（它的元素是任一点位置坐
有限元第五讲 平面问题(二)——离散化、三角形单元分析
该式建立了用单元节点位移表 达单元上应变分布的关系。
B 称为应变矩阵，其一个子块的计算式为：
(i l , m, n)
•
对简单三角形单元，应变矩阵为：
上面求出的待定系数 a1
~ a6 代回位移多项式，得到：
al 1 y bl 2 cl am bm cm a n ul bn um u cn n
u 1 x
ul 1 al bl x cl y am bm x cm y an bn x cn y um 2 u n ul N l N m N n um N l ul N mum N nun N i ui i l , m , n u n
xl xn yl ym yn
am bm cm
an ul bn um u cn n
2 1 xm
为三角形面积
节点坐标行列式
ak , bk , ck 分别是节点坐标行列式 的第k (k l，m，n)行第1， 2， 3个 元素的代数余子式，均为常数。
~ a3 : yl a1 ym a2 a yn 3
1
a1 1 xl a2 1 xm a 1 x n 3
其中：
yl ym yn
1 1
ul al 1 bl um u 2 c n l
yl ym yn
1
vl al 1 bl vm v 2 c n l
平面三角形3节点有限元程序
平面三角形３结点有限元程序1、程序名：，2、程序功能本程序能计算弹性力学的平面应力问题和平面应变问题；能考虑自重和结点集中力两种荷载的作用，在计算自重时y轴取垂直向上为正；能处理非零已知位移，如支座沉降的作用。
主要输出的内容包括：结点位移、单元应力、主应力、第一主应力与x轴的夹角以及约束结点的支座反力。
程序采用Visual Fortran 编制而成，输入数据全部采用自由格式。
3、程序流程及框图图１-1 程序流程图图１-2 程序框图其中，各子程序的功能如下：INPUT——输入结点坐标、单元信息和材料参数；MR——形成结点自由度序号矩阵；FORMMA——形成指标矩阵MA（N）并调用其他功能子程序，相当于主控程序；DIV——取出单元的3个结点号码和该单元的材料号并计算单元的b i，c i等；MGK——形成整体劲度矩阵并按一维存放在SK（NH）中；LOAD——形成整体结点荷载列阵F；OUTPUT——输出结点位移或结点荷载；TREAT——由于有非零已知位移，对K和F进行处理；DECOMP——整体劲度矩阵的分解运算；FOBA——前代、回代求出未知结点位移δ；ERFAC——计算约束结点的支座反力；KRS——计算单元劲度矩阵中的子块K rs。
4、输入数据及变量说明当程序开始运行时，按屏幕提示，键入数据文件的名字。
在运行程序之前，必须根据程序中输入要求建立一个存放原始数据的文件，这个文件的名字由少于8个字符或数字组成。
数据文件包括如下内容：⑴总控信息，共一条，9个数据NP，NE，NM，NR，NI，NL，NG，ND，NCNP——结点总数；NE——单元总数；NM——材料类型总数；NR——约束结点总数；NI ——问题类型标识，0为平面应力问题，1为平面应变问题；NL ——受荷载作用的结点的数目；NG ——考虑自重作用为1，不计自重为0；ND ——非零已知位移结点的数目；NC ——要计算支座约束反力的结点数目。
⑵ 材料信息，共NM 条，每条依次输入EO ，VO ，W ，tEO ——弹性模量（kN/m 2）；VO ——泊松比；W ——材料容重（kN/m 3）；t ——单元厚度（m ）。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
. 一、题目如图1所示，一个厚度均匀的三角形薄板，在顶点作用沿板厚方向均匀分布的竖向载荷。
已知：P=150N/m ，E=200GPa ，=0.25，t=0.1m ，忽略自重。
试计算薄板的位移及应力分布。
要求： 1. 编写有限元计算机程序，计算节点位移及单元应力。
（划分三角形单元，单元数不得少于30个）；2. 采用有限元软件分析该问题（有限元软件网格与程序设计网格必须一致），详细给出有限元软件每一步的操作过程，并将结果与程序计算结果进行对比（任选取三个点，对比位移值）；3. 提交程序编写过程的详细报告及计算机程序；4. 所有同学参加答辩，并演示有限元计算程序。
有限元法中三节点三角形分析结构的步骤如下：1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。
2）单元分析，建立单元刚度矩阵。
3）整体分析，建立总刚矩阵。
4）建立整体结构的等效节点荷载和总荷载矩阵5）边界条件处理。
6）解方程，求出节点位移。
7）求出各单元的单元应力。
8）计算结果整理。
一、程序设计网格划分如图，将薄板如图划分为6行，并建立坐标系，则刚度矩阵的集成建立与总刚度矩阵等维数的空矩阵，已变单元刚度矩阵的集成。
由单元分析已知节点、单元的排布规律，继而通过循环计算求得每个单元对应的节点序号。
通过循环逐个计算：（1）每个单元对应2种单元刚度矩阵中的哪一种；（2）该单元对应总刚度矩阵的那几行哪几列（3）将该单元的单元刚度矩阵加入总刚度矩阵的对应行列循环又分为3层循环：（1）最外层：逐行计算（2）中间层：该行逐个计算（3）最里层：区分为第 奇/偶 数个计算单元刚度的集成：[][][][][][]''''''215656665656266256561661e Z e e e Z e Z e e e e k k k K k k k k k k +⋯++=⇓=⇒==⇒==⇒=⨯⨯⨯⨯⨯⨯边界约束的处理：划0置1法X Y P X Y P适用：这种方法适用于边界节点位移分量为已知(含为0)的各种约束。
做法：（1）将总刚矩阵〔K〕中相应于已知位移行主对角线元素置1，其他元素改为零；同时将载荷列阵｛R｝中相应元素用已知位移置换。
◎这样，由该方程求得的此位移值一定等于已知量。
（2）将〔K〕中已知位移相应的列的非主对角成元素也置0，以保持〔K〕的对称性。
◎当然，在已知位移分量不为零的情况下，这样做就改变了方程左端的数值，为保证方程成立，须在方程右端减去已知位移对该方程的贡献——已知位移和相应总刚元素的乘积。
◎若约束为零位移约束时，此步则可省去。
特点：（1）经以上处理同样可以消除刚性位移(约束足够的前提下)，去掉未知约束反力。
（2）但这种方法不改变方程阶数，利于存贮。
（3）不过，若是要求出约束反力，仍要重新计算各个划去的总刚元素。
程序如下：变量说明NNODE 单元节点数NPION 总结点数NELEM 单元数NVFIX 受约束边界点数FIXED 约束信息数组NFORCE 节点力数FORCE 节点力数组COORD 结构节点坐标数组LNODS 单元定义数组YOUNG 弹性模量POISS 泊松比THICK 厚度B 单元应变矩阵(3*6)D 单元弹性矩阵(3*3)S 单元应力矩阵(3*6)A 单元面积ESTIF 单元刚度矩阵ASTIF 总体刚度矩阵ASLOD 总体荷载向量ASDISP 节点位移向量ELEDISP 单元节点位移向量STRESS 单元应力%********************************************************** %初始化clearformat short e %设定输出类型clear %清除内存变量NELEM=36 %单元个数（单元编码总数）NPION=28 %结点个数（结点编码总数）NVFIX=2 %受约束边界点数NFORCE=1 %结点荷载个数YOUNG=2e11 %弹性模量POISS=0.25 %泊松比THICK=0.1 %厚度LNODS=[1 2 3;2 4 5;2 5 3;3 5 6;4 7 8;4 8 5;5 8 9;5 9 6;6 9 10;7 11 12;7 12 8;8 12 13;8 13 9;9 13 14;9 14 10;10 14 15;11 16 17;11 17 12; 12 17 18; 12 18 13;13 18 19; 13 19 14;14 19 20;14 20 15;15 20 21;16 22 23;16 23 17;17 23 24;17 24 18;18 24 25;18 25 19;25 19 26;19 26 20;20 26 27;20 27 21;21 27 28] %单元定义数组（单元结点号）%相应为单元结点号（编码）、按逆时针顺序输入COORD=[0 0;-0.75 1.5;0.75 1.5;-1.5 3;0 3;1.5 3;-2.25 4.5;-0.75 4.5;0.75 4.5;2.25 4.5;-3 6;-1.5 6;0 6;1.5 6;3 6;-3.75 7.5;-2.25 7.5; -0.75 7.5;0.75 7.5;2.25 7.5;3.75 7.5;-4.5 9;-3 9;-1.5 9;0 9;1.5 9;3 9;4.5 9] %结点坐标数组%坐标：x，y 坐标（共NPOIN 组）FORCE=[1 0 -15] %结点力数组（受力结点编号, x 方向,y 方向）FIXED=[22 1 1;28 1 1] %约束信息（约束点，x 约束，y 约束）%有约束为1，无约束为0%********************************************************** %生成单元刚度矩阵并组成总体刚度矩阵ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0%********************************************************** for i=1:NELEM%生成弹性矩阵 DD= [1 POISS 0;POISS 1 0;0 0 (1-POISS)/2]*YOUNG/(1-POISS^2)%********************************************************** %计算当前单元的面积A=-det([1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2);1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2);1 COORD(LNODS(i,3),1) COORD(LNODS(i,3),2)])/2%********************************************************** %生成应变矩阵 Bfor j=0:2b(j+1)=COORD(LNODS(i,(rem((j+1),3))+1),2)-COORD(LNODS(i,(re m((j+2),3))+1),2);c(j+1)=-COORD(LNODS(i,(rem((j+1),3))+1),1)+COORD(LNODS(i,(r em((j+2),3))+1),1);endB=[b(1) 0 b(2) 0 b(3) 0;0 c(1) 0 c(2) 0 c(3);c(1) b(1) c(2) b(2) c(3) b(3)]/(2*A);B1( :,:,i)=B;%********************************************************** %求应力矩阵S=D*BS=D*B;ESTIF=B'*S*THICK*A; %求解单元刚度矩阵a=LNODS(i,:); %临时向量,用来记录当前单元的节点编号for j=1:3for k=1:3ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=ASTIF((a(j)*2-1) :a(j)*2,(a(k)*2-1):a(k)*2)+ESTIF(j*2-1:j*2,k*2-1:k*2);%根据节点编号对应关系将单元刚度分块叠加到总刚%度矩阵中endendend%********************************************************** %将约束信息加入总体刚度矩阵（对角元素改一法）for i=1:NVFIXif FIXED(i,2)==1ASTIF(:,(FIXED(i,1)*2-1))=0; %一列为零ASTIF((FIXED(i,1)*2-1),:)=0; %一行为零ASTIF((FIXED(i,1)*2-1),(FIXED(i,1)*2-1))=1; %对角元素为 1end%********************************************************** %生成单元刚度矩阵并组成总体刚度矩阵%**********************************************************if FIXED(i,3)==1ASTIF( :,FIXED(i,1)*2)=0; %一列为零ASTIF(FIXED(i,1)*2,:)=0; %一行为零ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1; %对角元素为 1endend%********************************************************** %生成荷载向量ASLOD(1:2*NPION)=0; %总体荷载向量置零for i=1:NFORCEASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);end%********************************************************** %求解内力ASDISP=ASTIF\ASLOD' %计算节点位移向量ELEDISP(1:6)=0; %当前单元节点位移向量for i=1:NELEMfor j=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);%取出当前单元的节点位移向量endiSTRESS=D*B1(:, :, i)*ELEDISP' %求内力end（程序计算结果和有限元软件得出的结果稍有偏差，可能是程序某些地方数据输入时出了问题，还在寻找具体原因）二、有限元软件分析设置材料参数建模网格划分。