《有限单元法》编程作业

0
y 0
Ni
0
Nj
0
Nm
y x
1 2A
bi
0
ci
0 ci bi
bj 0
cj
0 cj bj
bm 0
cm
0
cm
bm
(2-1-7)
应力矩阵[S]：
Biblioteka Baidu
S
D
B
bi
ci
bj
cj
bm
cm
E 21
2
A
bi
1 2
ci
ci
1 2
bi
bj
1 2
c
j
cj
1 2
bj
bm
1 2
cm
cm
输入其它数据 形成整体刚度阵 形成节点荷载向量 引入支承条件
1 2
6 单元刚度矩阵
7 单元面积
3 8
解方程，输出位移
9
10 求应力，输出应力
结束
图 1 整个程序流程图
3、程序主要标示符及变量说明
1、变量说明： Node ------- 节点定义 gElement ---- 单元定义 gMaterial --- 材料定义，包括弹性模量，泊松比和厚度 gBC1 -------- 约束条件 gNF --------- 集中力 gk------------总刚 gDelta-------结点位移
E
1
E
2
1
5.程序的验证
采用一悬臂梁作为考题，悬臂梁的尺寸为 5m 1m 1m ，端部承受 F 100 KN 集中
力的作用。取 E 200 GPa ， 0.3 。将计算结果和理论计算结果进行对比。对 比截面为悬臂端截面和距离支座 1/2 跨度的截面下端的位移和截面上下缘的应 力。为验证划分单元数多少会对计算的精度产生影响，现将该悬臂梁划分成 22 单元和 40 单元两种情况下进行计算。
1、程序编制总说明
a.该程序采用平面三角形等参单元，能解决弹性力学的平面应力、平 面应变问题。
b.能计算单元受集中力的作用。 c.能计算结点的位移和单元应力。 d.考题计算结果与理论计算结果比较，并给出误差分析。 e.程序采用 MATLAB R2008a 编制而成。
2、Matlab 程序编制流程图
开始 输入结构控制参数
将三角形单元的位移函数用矩阵表示：
f
(x,
y)
u
v
Ni
0
0
Ni
Nj 0 0 Nj
u i vi
Nm
0
N0 u j mv j
u
m
vm
由节点位移求应变──几何矩阵[B] 将 {f}=[N]{d}e 代入:
H N de [B]de
其中：
x
0
B
H
N
0
Ni
0 Nj
0
Nm
4 求弹性矩阵
位移-应变矩阵 5
2、包含的子程序说明： PlaneStructualModel ———定义有限元模型 SolveModel ———————求解有限元模型 DisplayResults ——————显示计算结果 k = StiffnessMatrix( ie )———计算单元刚度 AssembleStiffnessMatrix( ie, k )—形成总刚 es = ElementStress( ie )————计算单元应力
-1.13998604e+006 -1.25961706e+005 -8.37330822e+005 -9.08521263e+004 -5.57259307e+005 -5.60915016e+004 -2.75500197e+005 -3.21284714e+004 -6.10434315e+004 1.49745101e+006 1.29688331e+005 1.13999174e+006 1.25936919e+005 8.37334870e+005 9.08600989e+004 5.58019942e+005 5.29239673e+004 2.83418910e+005 2.40288921e+003 1.09482831e+005
1 2
bm
单元刚度矩阵为： K t BT D B dxdy
将上式积分后可得：
子块：
Kii
Kij
Kim
K
K
ji
K jj
K
jm
Kmi Kmj Kmm
Krs
Et 4(1 2)
A
brbs
12
crcs
crbs
12
brcs
brcs
12
crbs
crcs 12 brbs
上述是平面应力问题的单刚，对于平面应变问题，只要使用以下变换：
4、理论基础和求解过程
4.1、构造插值函数：
uuvNvui(((xxx,,,yyy)))ui 41Nj(x52,xxy)uj 63 yyNm(x, y)um
Ni
1 2A
(ai
bi
x
ci
y)
(i, j,m)
v Ni(x, y)vi N j(x, y)vj Nm(x, y)vm
同理可得： 4.2 位移插值函数及应变应力求解：
②第二次将该悬臂梁划分成 40 个单元，33 个节点。如下图所示。
图 3 悬臂梁 40 单元划分示意图
计算结果如下：
1.输入节点坐标： [x y] [0,0;0,0.5;0,1;0.5,0;0.5,0.5;0.5,1;1,0;1,0.5;1,1;1.5,0;1.5,0.5;1.5,1;2,0;2 ,0.5;2,1;2.5,0;2.5,0.5;2.5,1;3,0;3,0.5;3,1;3.5,0;3.5,0.5;3.5,1;4,0;4,0.5;4 ,1;4.5,0;4.5,0.5;4.5,1;5,0;5,0.5;5,1]
4.1、构造插值函数 ..................................... 5 4.2 位移插值函数及应变应力求解......................... 5 5.程序的验证 ............................................. 6 附录：程序代码 .......................................... 15
;3.5,1;4,0.5;4.5,0;4.5,1;5,0;5,0.5;5,1]
2.输入单元信息
各个单元包括的节点号 材料号： [1,4,2,1;4,6,2,1;4,7,6,1;7,9,6,1;7,10,9,1;10,12,9,1;10,13,12,1;13,15,12, 1;13,16,15,1;16,19,15,1;16,18,19,1;2,5,3,1;2,6,5,1;6,8,5,1;6,9,8,1;9,11, 8,1;9,12,11,1;12,14,11,1;12,15,14,1;15,17,14,1;15,19,17,1;19,20,17,1]
湖南大学 《有限单元法》编程大作业
专业：土木工程 姓名： 学号： 2013 年 12 月
目录
程序作业题目： ........................................... 3 1、程序编制总说明 ........................................ 3 2、Matlab 程序编制流程图.................................. 3 3、程序主要标示符及变量说明 .............................. 4 4、理论基础和求解过程 .................................... 5
3.材料参数 弹性模量 泊松比 厚度：
[2e11,0.3,1]
4.约束 节点号，自由度号，约束值: [1,1,0;1,2,0;2,1,0;2,2,0;3,1,0;3,2,0]
5.荷 载 ： 节点号 自由度号 集中力值
[20,2,-1e5]
6.程序计算输出的节点位移，各个单元应力：
节点位移
节点号
x 方向位移
1 -3.37036821e-021
2 -2.84383452e-027
3
3.37037253e-021
4 -3.40670099e-006
5
3.40670104e-006
6
7.52240584e-012
7 -9.79666516e-006
8
9.79665910e-006
9 -8.58390708e-011
程序作业题目：
完成一个包含以下所列部分的完整的有限元程序( Project) 须提供如下内容的文字材料（1500 字以上）：
① 程序编制说明； ② 方法的基本理论和基本公式； ③ 程序功能说明； ④ 程序所用主要标识符说明及主要流程框图； ⑤ 1～3 个考题：考题来源、输出结果、与他人成果的对比结果（误差 百分比）； ⑥ 对程序的评价和结论（包括正确性、适用范围、优缺点及其他心得 等）。 须提供源程序、可执行程序和算例的电子文档或文字材料。选题可根据各自的论 文选题等决定。
①第一次将该悬臂梁划分成 22 个单元，共 20 个节点。划分示意图如下所示。
计算结果如下：
图 2 悬臂梁 22 单元划分示意图
1.输入节点坐标 [x y] [0,0;0,0.5;0,1;0.5,0;0.5,1;1,0.5;1.5,0;1.5,1;2,0.5;2.5,0;2.5,1;3,0.5;3.5,0
4.60022636e+005 -9.79789911e+004 -4.19810113e+005 -6.18504870e+004 3.62593854e+005 -8.09752548e+004 -3.03024973e+005 -5.92052900e+004 2.43166517e+005 -8.06856809e+004 -1.80611027e+005 -5.98498412e+004 1.18743602e+005 -8.13737464e+004 -6.95369585e+004 -6.65415096e+004 -6.10434315e+004 -6.10434315e+004 4.49235302e+005 -5.02550275e+005 4.32289421e+005 -1.29744125e+005 -4.59999584e+005 -9.79710172e+004 4.19851971e+005 -6.18393692e+004 -3.62740453e+005 -8.10071118e+004 3.02682445e+005 -5.91854535e+004 -2.41554234e+005 -7.99369095e+004 1.82773870e+005 -5.81917136e+004 -1.22337290e+005 -7.82806950e+004 4.55675769e+004 -4.78276920e+004 -2.90517169e+005 -1.09482831e+005
19 -6.56868986e-008
20 1.97720706e-005
y 方向位移 -1.60393066e-021 1.02108017e-021 -1.60393495e-021 -3.26656027e-006 -3.26657679e-006 -8.50006819e-006 -1.77436683e-005 -1.77435451e-005 -2.88973765e-005 -4.31200502e-005 -4.31209125e-005 -5.87215705e-005 -7.63748124e-005 -7.63670297e-005 -9.48469722e-005 -1.14326100e-004 -1.14363730e-004 -1.33997894e-004 -1.34104720e-004 -1.34913125e-004
单元应力
单元号 1 单元号 2
X-STR -1.49745098e+006 -1.29687202e+005
Y-STR
XY-STR
-4.49235295e+005 -5.02547734e+005
-4.32295689e+005 -1.29742976e+005
单元号 3 单元号 4 单元号 5 单元号 6 单元号 7 单元号 8 单元号 9 单元号 10 单元号 11 单元号 12 单元号 13 单元号 14 单元号 15 单元号 16 单元号 17 单元号 18 单元号 19 单元号 20 单元号 21 单元号 22
10 -1.45272100e-005
11 1.45274441e-005
12 1.90988361e-010
13 -1.76782564e-005
14 1.76798752e-005
15 -9.34997929e-009
16 -1.92338727e-005
17 1.92804757e-005
18 -1.93406988e-005