岩土工程有限元大作业

岩土工程有限元大作业
题目：平面问题3结点三角形有限元法
的Matlab程序实现
专业：建筑与土木工程（岩土方向）
姓名：
学号：_________
2014年12 月31 日
1.概述
对于二维平面问题，运用有限单元法求解时，单元的类型有3结点三角形单元、4结点矩形单元以及高精度的三角形单元（6结点、10结点）等。

三角形单元对复杂边界有较强的适应能力，因此很容易将一个二维域离散成有限个三角形单元，在边界上以若干段直线近似原来的曲线边界，随着单元增多，这种拟合将越精确。

本文即基于相对比较简单且运用广泛的3结点三角形单元的有限元法编制Matlab 程序，解决一个悬臂梁端点承受竖向集中荷载作用时的位移计算问题，并绘出变形前后的图形作对比。

在分析过程中，将此问题简化为平面应力问题，取单位厚度进行分析。

以下一共分为两部分介绍，第一部分是理论部分，介绍三角形有限元法的基本求解过程；第二部分是程序部分，介绍具体的Matlab 程序实现过程。

2.三角形有限元法的基本求解过程 2.1 假设单元的位移函数
如图2.1所示为一端部承受竖向集中荷载作用的悬臂梁，将此问题简化为二维平面问题，梁取单位厚度，固定端简化为两个不动铰支座，已知悬臂梁的长度为8米，高度为1米，材料为钢材，它的弹性模量取2.06e11Pa ，泊松比取0.25，右端处作用的竖向荷载为100kN 。

二维域Ω被离散成如图所示的九个三角形单元，后续编程中的结点号和单元号都以此为基础，任取一三角形单元，设其结点编码为m j i ,,，以逆时针编码为正向(见图2.2)。

该三角形单元在x 、y 方向的位移函数表达式为：
m m j j i i u N u N u N u ++= (2.1) m m j j i i v N v N v N v ++=
(2.2)
其中()m j i i N i ,, =称为单元的插值函数或形函数
()y c x b a A
N i i i i ++=
21
(2.3)
图2.1 任意区域的三角形单元离散 j (x j , y j )
x
y
i (x i , y i )
m (x m , y m )
图2.2 3结点三角形单元
),,(m j i x
x c y y b y x y x a m j i
m j i j m m j i ⎪⎩⎪
⎨⎧+-=-=-= (2.4)
上式),,(m j i 表示下标轮换，即i m m j j i →→→,,。

(2.1)和(2.2)式简记为
()∑==i
i i m j i i N ,, u u
(2.5)
(2.5)式的矩阵形式是
[
]
e m j i m j
i
v u Na a a a N N N u =⎪⎭
⎪
⎬⎫⎪⎩⎪⎨⎧=⎭
⎬⎫
⎩⎨⎧=
(2.6)
其中⎥⎦⎤⎢
⎣⎡=i i
i N N 00N ，⎭
⎬⎫⎩⎨⎧=i i i v u a 。

2.2 单元应变矩阵和弹性矩阵
将位移函数(2.6)式代入平面问题的几何方程得到单元的应变表达式
[
]
e e m j i
Ba a B B B ε==
(2.7)
其中B 称为应变矩阵，
⎥⎥⎥⎥⎥
⎥⎥⎦
⎤⎢⎢⎢⎢⎢
⎢⎢⎣⎡∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂=x N y
N y N x N x N y
N y
N x
N x N y
N y
N x N m m m m j j j
j i i i
i i 0000
B (2.8)
将(2.3)式代入上式得
⎥⎥⎥⎦
⎤⎢⎢
⎢⎣
⎡=m m
m m j j
j j i i i i
b c c b b c c b b c c b A 00000
021B (2.9)
对于平面应力问题，单元的弹性矩阵如下：
⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢
⎢⎢⎣
⎡
-=2100001
1ννν
D B (2.10)
其中2
1ν-=
E
D 。

2.3 整体平衡方程
根据最小势能原理或变分原理，对总势能取驻值，即0=δπ得
()()∑∑=e
e
e
e
e
F a K
(2.11)
其中，单元刚度矩阵
tdxdy e
e B D B K ⎰Ω
=T
(2.12)
单元载荷列阵
⎰⎰⎰ΩΩ++-=e e
e
s e tds tdxdy tdxdy σ
T N f N σB F T T 0T
(2.13)
本文分析的问题只是在悬臂梁右端某一结点处作用一竖向荷载，荷载列阵相对比较简单，如果涉及到分布荷载，需要局部坐标到整体坐标的转换，形成整体荷载列阵。

2.4 组集单元矩阵到整体方程
按对号入座的原则对单元平衡方程进行组集得到整体求解的结构控制方程
F Ka = (2.14)
2.5 引入边界条件
边界条件包括外力边界条件和位移边界条件，对于本文分析的问题，外力边界条件主要是在右端点处有一个y 方向的集中荷载。

位移边界条件主要是左端的两个结点处在x 、y 两个方向位移为零。

2.6 求解整体离散方程
在求出单元刚度矩阵组集而成的整体刚度矩阵和单元等效荷载列阵以后，求解整体离散方程就变成求解线性方程组的问题了。

经过组集单元矩阵到整体矩阵形成2.14式中的整体刚度矩阵K 以及整体荷载列阵F ，求解位移矩阵就只需要用到Matlab 中的左除符号“\”即可。

3.Matlab程序实现过程
Matlab的求解程序应该包括一个主程序和若干个子程序，若干个子程序分别实现输入数据、计算单元弹性矩阵、应变矩阵、刚度矩阵以及单元刚度矩阵的组集、边界条件的处理、绘图等功能。

另外在主程序中定义全局变量，当在各个子程序中需要使用全局变量时，只需在各子程序中用命令global声明即可。

以下是程序中使用的全局变量的说明（首字母g代表global，表示全局变量Global Variables）
gMP------材料属性（Material Properties）
gNC------结点坐标（Node Coordinates）
gENN------单元结点编号（Element Node Number）
gBC------边界条件（Boundary Conditions）
gKA-------整体刚度矩阵
gFA-------整体荷载列阵
gND------计算得到的结点位移值（Nodal Displacement）
3.1求解主程序
3.2输入数据
3.3计算单元弹性矩阵子程序
3.5计算单元刚度矩阵子程序
3.7计算整体刚度矩阵子程序
3.8处理边界条件子程序3.8.1位移边界条件
3.8.2外力边界条件
3.9绘制变形前后的对比图子程序
说明：trisurf(Tri,x,y,z,c)三角表面图函数，用于生成三角形表面图。

Tri矩阵的每一行表示一个三角形的三个结点的结点编号，x,y,z表示三个结点的结点坐标。

C用来控制颜色。

3.10程序运行结果
3.10.1位移计算结果
3.10.1变形前后对比图（变形放大了1000倍）
4.结论
有限单元法在工程分析中应用广泛，它的基本思想是设法将实际上是无穷多自由度的连续介质问题近似的简化为由有限个“结点”构成的有限个自由度问题，并以这些结点的“自由度”为未知量，设法将控制方程近似的化为一组线性代数方程，然后用计算机求解。

本文通过平面问题三结点三角形有限单元法编制Maltab程序，实现悬臂梁承受集中荷载时的变形计算。

通过编制程序，认识到了有限元求解问题的一般流程，加深了对课程内容的理解与掌握。

在编程的过程中也遇到很多问题，比如这个程序中的核心部分——单元刚度的组集问题，各个单元刚度矩阵按照对号入座的原则组集成整体刚度矩阵，这里要使用两个循环语句才能实现。

另外结点固定的位移边界条件如何处理也不太清楚，在查阅相关有限元书籍后，知道了有划零置一法和乘大数法等，本文中采用的是划零置一法，即在整体矩阵中将与零节点位移相对应的行列中，将主对角元素改为以1，其他元素改为0，在荷载列阵中将与零节点位移相对应的元素改为0。

这种方法引入强制边界条件比较简单，不改变原来方程的阶数和未知量的顺序编号。

最后需要说明的是本文的荷载边界条件是点荷载，编程比较容易处理，进一步的研究中可以考虑均匀分布荷载以及变化的荷载作用下，程序如何编写。

10。