基于Matlab语言的按平面三角形单元划分的结构有限元程序设计.docx
《基于Matlab语言的按平面三角形单元划分的结构有限元程序设计.docx》由会员分享，可在线阅读，更多相关《基于Matlab语言的按平面三角形单元划分的结构有限元程序设计.docx（24页珍藏版）》请在冰豆网上搜索。
基于Matlab语言的按平面三角形单元划分的结构有限元程序设计
基于语言的按平面三角形单元划分的结构有限元程序设计
专业：
建筑与土木工程
班级：
建工研12-2
姓名：
韩志强
学号：
471220580
基于语言的按平面三角形单元划分结构有限元程序设计
一、有限单元发及语言概述
1.有限单元法
随着现代工业、生产技术的发展，不断要求设计高质量、高水平的大型、复杂和精密的机械及工程结构。
为此目的，人们必须预先通过有效的计算手段，确切的预测即将诞生的机械和工程结构，在未来工作时所发生的应力、应变和位移因此，需要寻求一种简单而又精确的数值分析方法。
有限单元法正是适应这种要求而产生和发展起来的一种十分有效的数值计算方法。
有限元法把一个复杂的结构分解成相对简单的“单元”，各单元之间通过结点相互连接。
单元内的物理量由单元结点上的物理量按一定的假设内插得到，这样就把一个复杂结构从无限多个自由度简化为有限个单元组成的结构。
我们只要分析每个单元的力学特性，然后按照有限元法的规则把这些单元“拼装”成整体，就能够得到整体结构的力学特性。
有限单元法基本步骤如下：
（1）结构离散：
结构离散就是建立结构的有限元模型，又称为网格划分或单元划分，即将结构离散为由有限个单元组成的有限元模型。
在该步骤中，需要根据结构的几何特性、载荷情况等确定单元体内任意一点的位移插值函数。
（2）单元分析：
根据弹性力学的几何方程以及物理方程确定单元的刚度矩阵。
（3）整体分析:
把各个单元按原来的结构重新连接起来，并在单元刚度矩阵的基础上确定结构的总刚度矩阵，形成如下式所示的整体有限元线性方程：
①
式中，
是载荷矩阵,
是整体结构的刚度矩阵,
是节点位移矩阵。
（4）载荷移置:
根据静力等效原理，将载荷移置到相应的节点上，形成节点载荷矩阵。
（5）边界条件处理:
对式①所示的有限元线性方程进行边界条件处理。
（6）求解线性方程:
求解式①所示的有限元线性方程，得到节点的位移。
在该步骤中，若有限元模型的节点越多，则线性方程的数量就越多，随之有限元分析的计算量也将越大。
（7）求解单元应力及应变根据求出的节点位移求解单元的应力和应变。
（8）结果处理与显示进入有限元分析的后处理部分，对计算出来的结果进行加工处理，并以各种形式将计算结果显示出。
2.简介
在用有限元法进行结构分析时，将会遇到大量的数值计算，因而在实用上是一定要借助于计算机和有限元程序，才能完成这些复杂而繁重的数值计算工作。
而是当今国际科学界最具影响力和活力的软件。
它起源于矩阵运算，并已经发展成一种高度集成的计算机语言。
它提供了强大的科学计算，灵活的程序设计流程，高质量的图形可视化与界面设计，便捷的与其他程序和语言接口的功能。
在各国高校与研究单位起着重大的作用。
“工欲善其事，必先利其器”。
如果有一种十分有效的工具能解决在教学与研究中遇到的问题，那么语言正是这样的一种工具。
它可以将使用者从繁琐、无谓的底层编程中解放出来，把有限的宝贵时间更多地花在解决问题中，这样无疑会提高工作效率。
目前，已经成为国际上最流行的科学与工程计算的软件工具，现在的已经不仅仅是一个“矩阵实验室”了，它已经成为了一种具有广泛应用前景的全新的计算机高级编程语言了，有人称它为“第四代”计算机语言，它在国内外高校和研究部门正扮演着重要的角色。
语言的功能也越来越强大，不断适应新的要求提出新的解决方法。
可以预见，在科学运算、自动控制与科学绘图领域语言将长期保持其独一无二的地位。
为此，本例采用语言编程，以利用其快捷强大的矩阵数值计算功能。
二、问题描述
一矩形薄板，一边固定，承受150集中荷载，结构简图及按平面三角形单元划分的有限元模型图如下所示。
材料参数：
弹性模量
；泊松比：
；薄板厚度
。
在本例中，所取结构模型及数据主要用于程序设计理论分析，与工程实际无关。
三、参数输入：
单元个数4
节点个数6
受约束边界点数2
节点荷载个数1
弹性模量2e8
泊松比0.2
厚度0.002
单元节点编码数组=
节点坐标数组=
节点力数组=[40-150]
约束信息数组=
以上数值数据为程序运行前输入的初始数据，存为“471220580”文本格式，同时必须放在工作目录下，路径不对程序不能自动读取指定初始文件，运行出错。
初始数据文本文件输入格式如下图：
四、语言程序源代码：
1.程序中变量说明
单元节点数
总结点数
单元数
受约束边界点数
约束信息数组
节点力数
节点力数组
结构节点坐标数组
单元定义数组
弹性模量
泊松比
厚度
B单元应变矩阵（3*6）
D单元弹性矩阵（3*3）
S单元应力矩阵（3*6）
A单元面积
单元刚度矩阵
总体刚度矩阵
总体荷载向量
节点位移向量
单元节点位移向量
单元应力
1数据文件指针
2.语言程序代码
%******************************************************************************
%初始化及数据调用
e%设定输出类型
%清除内存变量
1（'471220580',''）;%打开输入数据文件，读入控制数据
（1,'',1）,%单元个数（单元编码总数）
（1,'',1）,%结点个数（结点编码总数）
（1,'',1）%受约束边界点数
（1,'',1）,%结点荷载个数
（1,'',1）,%弹性模量
（1,'',1）,%泊松比
（1,'',1）%厚度
（1,'',[3]）'%单元定义数组（单元结点号）
%相应为单元结点号（编码）、按逆时针顺序输入
（1,'',[2]）'%结点坐标数组
%坐标：
x，y坐标（共组）
（1,'',[3]）'%结点力数组（受力结点编号,x方向方向）
（1,'',[3]）'%约束信息（约束点，x约束，y约束）
%有约束为1，无约束为0
%*****************************************************************************
%生成单元刚度矩阵并组成总体刚度矩阵
（2*,2*）;%生成特定大小总体刚度矩阵并置0
%*****************************************************************************
1
%生成弹性矩阵D
[10;
10;
00
（1）/2]*（1^2）
%*****************************************************************************
%计算当前单元的面积A
（[1（（i,1）,1）（（i,1）,2）;
1（（i,2）,1）（（i,2）,2）;
1（（i,3）,1）（（i,3）,2）]）/2
%*****************************************************************************
%生成应变矩阵B
0:
2
b
（1）=（（i,（（
（1）,3））+1）,2）（（i,（（
（2）,3））+1）,2）;
c
（1）（（i,（（
（1）,3））+1）,1）（（i,（（
（2）,3））+1）,1）;
[b
（1）0b
（2）0b（3）0;...
0c
（1）0c
（2）0c（3）;...
c
（1）b
（1）c
（2）b
（2）c（3）b（3）]/（2*A）;
B1（）;
%*****************************************************************************
%求应力矩阵*B
*B;
'*S**A;%求解单元刚度矩阵
（i,:
）;%临时向量,用来记录当前单元的节点编号
1:
3
1:
3
（（a（j）*2-1）（j）*2,（a（k）*2-1）（k）*2）…
（（a（j）*2-1）（j）*2,（a（k）*2-1）（k）*2）（j*2-1*2*2-1*2）;
%根据节点编号对应关系将单元刚度分块叠加到总刚
%度矩阵中
%*****************************************************************************
%将约束信息加入总体刚度矩阵（对角元素改一法）
1
（i,2）1
（:
（（i,1）*2-1））=0;%一列为零
（（（i,1）*2-1）,:
）=0;%一行为零
（（（i,1）*2-1）,（（i,1）*2-1））=1;%对角元素为1
（i,3）1
（（i,1）*2）=0;%一列为零
（（i,1）*2,:
）=0;%一行为零
（（i,1）*2（i,1）*2）=1;%对角元素为1
%*****************************************************************************
%生成荷载向量
（1:
2*）=0;%总体荷载向量置零
1
（（（i,1）*2-1）（i,1）*2）（i,2:
3）;
%*****************************************************************************
%求解内力
\'%计算节点位移向量
（1:
6）=0;%当前单元节点位移向量
1
1:
3
（j*2-1*2）（（）*2-1（）*2）;
%取出当前单元的节点位移向量
i
*B1（:
:
i）*'%求内力
（1）;%关闭数据文件
五、程序运行结果
=
4
=
6
=
2
=
1
=
200000000
=
2.0000001
=
2.0000003
=
126
234
245
256
=
00
10
20
21
11
01
=
40-150
=
111
611
D=
2.08330084.16670070
4.16670072.08330080
008.3333007
A=
5.0000001
D=
2.08330084.16670070
4.16670072.08330080
008.3333007
A=
5.0000001
D=
2.08330084.16670070
4.16670072.08330080
008.3333007
A=
5.0000001
D=
2.08330084.16670070
4.16670072.08330080
008.3333007
A=
5.0000001
=
0
0
-8.0607004
-1.5848003
-1.0281003
-4.4727003
1.1937003
-4.6947003
8.6670004
-1.8880003
0
0
（说明：
以上12个输出结果数据表示节点位移向量，即各节点分别在X，Y方向上产生的位移。
）
i=
1
=
-1.6793005
-3.3586004
-1.3207005
i=
2
=
-5.5503004
-5.5503004
-5.5503004
i=
3
=
5.5503004
-4.9526004
-9.4497004
i=
4
=
1.6793005
-2.7040004
-1.7932004
（说明：
以上四组输出结果数据表示单元应力，i为单元号。
）
六、建模比较分析
利用有限元分析软件，完全按照前面程序设计的各项条件参数以及单元划分方式建立模型，其求解结果与以上程序计算结果比较，验证程序是否正确。
模型节点单元如下图所示：
模型求解变形图如下所示：
求解节点位移结果列表显示如下：
（单位：
m）
U
*****1*****
10.00000.00000.00000.0000
2-0.8060703-0.15848020.00000.1778002
3-0.1028102-0.44727020.00000.4589302
40.1193702-0.46947020.00000.4844102
50.8667003-0.18880020.00000.2077402
60.00000.00000.00000.0000
4404
0.1193702-0.46947020.00000.4844102
求解单元应力结果列表显示如下：
（单位：
2）
S
*****1*****
1182
1-0.1679306-33586.0.0000-0.13207060.0000
2-0.1679306-33586.0.0000-0.13207060.0000
6-0.1679306-33586.0.0000-0.13207060.0000
2182
2-55503.-55503.0.0000-55503.0.0000
3-55503.-55503.0.0000-55503.0.0000
4-55503.-55503.0.0000-55503.0.0000
3182
255503.-49526.0.0000-94497.0.0000
455503.-49526.0.0000-94497.0.0000
555503.-49526.0.0000-94497.0.0000
4182
20.1679306-27040.0.0000-17932.0.0000
50.1679306-27040.0.0000-17932.0.0000
60.1679306-27040.0.0000-17932.0.0000
结论
通过比较语言设计程序运行结果和建模分析结果可知，两种方式算出的结果完全一致，说程序设计正确。
所以本程序适用于按三角形单元划分的平面结构有限元分析。
e
1（'',''）;
（1,'',1）
（1,'',1）
（1,'',1）
（1,'',1）
（1,'',1）
（1,'',1）
（1,'',1）
（1,'',[3]）
（1,'',[2]）
（1,'',[3]）
（1,'',[3]）
（2*,2*）;%生成特定大小总体刚度矩阵并置0
1
%生成弹性矩阵D
（1^2）*[10;
10;
00
（1）/2]
%计算当前单元的面积A
（[1（（i,1）,1）（（i,1）,2）;
1（（i,2）,1）（（i,2）,2）;
1（（i,3）,1）（（i,3）,2）]）/2
%生成应变矩阵B
0:
2
b
（1）=（（i,（（
（1）,3））+1）,2）（（i,（（
（2）,3））+1）,2）;
c
（1）（（i,（（
（1）,3））+1）,1）（（i,（（
（2）,3））+1）,1）;
[b
（1）0b
（2）0b（3）0;...
0c
（1）0c
（2）0c（3）;...
c
（1）b
（1）c
（2）b
（2）c（3）b（3）]/（2*A）;
B1（）;
%求应力矩阵*B
*B;
'*S**A;
（i,:
）;
1:
3
1:
3
（（a（j）*2-1）（j）*2,（a（k）*2-1）（k）*2）…
（（a（j）*2-1）（j）*2,（a（k）*2-1）（k）*2）（j*2-1*2*2-1*2）;
%将约束信息加入总体刚度矩阵
1
（i,2）1
（:
（（i,1）*2-1））=0;
（（（i,1）*2-1）,:
）=0;
（（（i,1）*2-1）,（（i,1）*2-1））=1;
（i,3）1
（（i,1）*2）=0;
（（i,1）*2,:
）=0;
（（i,1）*2（i,1）*2）=1;
%生成荷载向量
（1:
2*）=0;
1
（（（i,1）*2-1）（i,1）*2）（i,2:
3）
%求解内力
\'
（1:
6）=0;
1
1:
3
（j*2-1*2）（（）*2-1（）*2）;
i
*B1（:
:
i）*'
（1）;
e
1（'',''）;
（1,'',1）
（1,'',1）
（1,'',1）
（1,'',1）
（1,'',1）
（1,'',1）
（1,'',1）
（1,'',[,3]）
（1,'',[,2]）
（1,'',[,3]）
（1,'',[,3]）
（2*,2*）;%生成特定大小总体刚度矩阵并置0
1
%生成弹性矩阵D
（1^2）*[10;
10;
00
（1）/2]
%计算当前单元的面积A
（[1（（i,1）,1）（（i,1）,2）;
1（（i,2）,1）（（i,2）,2）;
1（（i,3）,1）（（i,3）,2）]）/2
%生成应变矩阵B
0:
2
b
（1）=（（i,（（
（1）,3））+1）,2）（（i,（（
（2）,3））+1）,2）;
c
（1）（（i,（（
（1）,3））+1）,1）（（i,（（
（2）,3））+1）,1）;
[b
（1）0b
（2）0b（3）0;...
0c
（1）0c
（2）0c（3）;...
c
（1）b
（1）c
（2）b
（2）c（3）b（3）]/（2*A）;
B1（）;
%求应力矩阵*B
*B;
'*S**A;
（i,:
）;
1:
3
1:
3
（（a（j）*2-1）（j）*2,（a（k）*2-1）（k）*2）…
（（a（j）*2-1）（j）*2,（a（k）*2-1）（k）*2）（j*2-1*2*2-1*2）;
%将约束信息加入总体刚度矩阵
1
（i,2）1
（:
（（i,1）*2-1））=0;
（（（i,1）*2-1）,:
）=0;
（（（i,1）*2-1）,（（i,1）*2-1））=1;
（i,3）1
（（i,1）*2）=0;
（（i,1）*2,:
）=0;
（（i,1）*2（i,1）*2）=1;
%生成荷载向量
（1:
2*）=0;
1
（（（i,1）*2-1）（i,1）*2）（i,2:
3）
%求解内力
\'
（1:
6）=0;
1
1:
3
（j*2-1*2）（（）*2-1（）*2）;
i
*B1（:
:
i）*'
（1）;
升级会员 