有限元钢架结构分析手算+matlab+ansys模拟.docx
《有限元钢架结构分析手算+matlab+ansys模拟.docx》由会员分享，可在线阅读，更多相关《有限元钢架结构分析手算+matlab+ansys模拟.docx（29页珍藏版）》请在冰豆网上搜索。
有限元钢架结构分析手算+matlab+ansys模拟
有限元大作业——钢架结构分析
选题人：
日期：
2016年6月2日
目录：
第一章：
问题重述1
一、题目容：
1
二、题目要求：
1
第二章：
有限元法手工求解2
一、平面两单元离散化2
二、单元分析2
三、单元组装5
四、边界条件引入及组装总体方程5
五、求解整体刚度方程，计算节点2的位移和转角6
六、求节点1、3支撑反力6
七、设定数据，求解结果7
八、绘制轴力图、弯矩图、剪力图8
第三章、matlab编程求解：
9
一、总体流程图绘制：
9
二、输入数据：
9
三、计算单元刚度矩阵：
10
四、建立总体刚度矩阵：
10
五、计算未约束点位移：
10
六、计算支反力：
10
七、输出数据：
10
八、编程：
10
第四章有限元求解11
一、预处理11
二、模型建立：
12
二、分析计算14
三、求解结果15
四、绘制图像16
第五章结果比较19
第六章心得体会19
一、王小灿：
19
二、明哲：
20
三、国威20
第七章附录22
一、matlab程序22
第一章：
问题重述
一、题目容：
图示平面钢架结构
图1.1题目容
二、题目要求：
（1）采用平面梁单元进行有限元法手工求解，要求写出完整的求解步骤，包括：
a）离散化：
单元编号、节点编号；
b）单元分析：
单元刚度矩阵，单元节点等效载荷向量；
c）单元组长：
总体刚度矩阵，总体位移向量，总体节点等效载荷；
d）边界条件的引入及总体刚度方程的求解；
e）B点的位移，A、C处支撑反力，并绘制该结构的弯矩图、剪力图和轴力图。
（2）编制通用平面钢架分析有限元Matlab程序，并计算盖提，与手工结果进行比较；
（3）利用Ansys求解，表格列出B点的位移，A、C处支反力，绘制弯矩图、剪力图和轴力图，并与手算和Matlab程序计算结果比较。
（4）攥写报告，利用A4纸打印；
（5）心得体会，并简要说明各成员主要负责完成的工作。
第二章：
有限元法手工求解
一、平面两单元离散化
将平面梁离散为两个单元，单元编号分别为①和②，节点号分别为1、2、3；
如图2-1所示：
图2-1单元离散化示意图
二、单元分析
首先建立整体坐标系与局部坐标系如图所示；
1、求单元刚度矩阵
对于单元①，求局部坐标系的单元刚度矩阵：
由于单元①局部坐标系与整体坐标系的夹角为：
，则单元①的局部坐标变换矩阵为：
可以得到在总体坐标系下的单元①的刚度矩阵：
对于单元②，求局部坐标系的单元刚度矩阵：
由于单元②局部坐标系与整体坐标系的夹角为
则
。
2、求单元节点等效载荷向量
将P等效在单元①两侧节点1,2上：
将均布载荷等效在单元②两侧的节点2，3上：
与作用在节点上的力叠加为整体坐标系下的节点载荷：
三、单元组装
将两个整体坐标系下的单元刚度矩阵组装为整体刚度矩阵：
四、边界条件引入及组装总体方程
由于节点1、3为固定约束，所以节点1和3的x、y方向的位移以及转角均为0，节点2无位移约束，不存在支反力，所以力约束即为外力约束。
五、求解整体刚度方程，计算节点2的位移和转角
提取节点2位移的相关要素：
求得：
六、求节点1、3支撑反力
根据总体方程，提取求解节点1支撑反力所需方程：
根据总体方程，提取求解节点2支撑反力所需方程：
七、设定数据，求解结果
设定各个数据：
氏模量：
泊松比：
力：
截面面积：
惯性矩：
将数据代入结果。
节点2的位移和转角：
节点1支撑反力：
节点3支撑反力：
八、绘制轴力图、弯矩图、剪力图
应用材料力学的分析方法，对梁单元进行分析。
轴力图：
图2-2轴力图
剪力图：
图2-3剪力图
弯矩图
图2-4弯矩图
第三章、matlab编程求解：
一、总体流程图绘制：
图3.1总体流程图
二、输入数据：
考虑到后续计算和以下参数相关：
节点个数，单元数，氏模量，惯性矩，单元长度，单元截面积，单元的旋转角度，节点与单元的对应关系，力与转矩的约束以及结构约束。
考虑到钢架结构，每个单元的氏模量，惯性矩，单元长度，单元截面积以及单元的旋转角度都可能不一样，所以采用矩阵的形式进行输入。
（注：
由于本题除长度外一样，故将其余几项改为常量进行计算）
单元与节点对应关系为：
一个单元对应2个节点，且按顺序连接。
力与转矩的约束以及结构约束：
应包括约束值，作用节点，作用类型，3种，并以作用节点与作用类型来反推此约束在完整的约束矩阵中的位置。
三、计算单元刚度矩阵：
图3.2单元刚度矩阵生成流程图
考虑到每个单元的刚度矩阵与坐标变换的矩阵形式相同，只是数据不同，故采取建立模板，利用eval（），函数来带入不同单元的值，生成一系列单元刚度矩阵，并用一个三维数组存储这些矩阵。
四、建立总体刚度矩阵：
考虑到每个单元刚度矩阵都是6×6的形式，表述了2个节点间的相互关系；故建立元胞数组，并使元胞数组的阶数与节点个数相同，利用元胞数组存储节点间关系。
首先建立与节点个数相同阶数的空元胞数组，之后检索每个单元刚度矩阵对应的2个节点间的关系，将其分离成4个3×3的矩阵，按节点与单元对应关系，存储到元胞数组中。
最后将元胞数组展开形成的大矩阵即为总体刚度矩阵。
五、计算未约束点位移：
利用总体位移与外力间的关系，采用矩阵求解，求取非约束点的位移。
并针对结果进行对应处理，使结果与作用点、作用形式对应。
六、计算支反力：
利用约束点位移皆零的特点，简化总体刚度矩阵，同时由于部分节点的部分方向上为力而非支反力，再度简化总体刚度矩阵。
利用两次简化后的刚度矩阵与计算出的位移结果相乘，求得不计直接作用在节点约束方向上时的支反力，将结果加上由于直接作用在节点约束方向上时产生的支反力，即为最后的支反力结果。
七、输出数据：
将计算所得的未约束点位移与支反力，采用与输入方式相似的方式进行处理并进行输出。
八、编程：
见附录一
第四章有限元求解
一、预处理
1、选择单元类型：
ANSYSMainMenu:
Preprocessor→ElementType→Add/Edit/Delete…→Add…→beam：
2Delastic3→OK（返回到ElementTypes窗口）→Close
图4.1选择单元类型
2、定义材料参数：
ANSYSMainMenu:
Preprocessor→MaterialProps→MaterialModels→Structural→Linear→Elastic→Isotropic:
EX:
3e10（弹性模量），PRXY:
0.3（泊松比）→OK
图4.2定义材料参数
3、定义单元截面积和惯性矩：
ANSYSMainMenu:
Preprocessor→Realconstant→Add→Typebeam3→Ok→Cross-sectionalareaAREA:
0.05（横截面积）AreamomentofinteiaIZZ:
1（惯性矩）→OK
图4.3定义单元截面积和惯性矩
二、模型建立：
1、画出关键点：
ANSYSMainMenu:
Preprocessor→Modeling→Creat→Keypoint→InActiveCS→Nodenumber1→X:
0,Y:
0,Z:
0→Apply→Nodenumber2→X:
0,Y:
1,Z:
→Apply→Nodenumber3→X:
2,Y:
1,Z:
0→OK
2、构造连线：
ANSYSMainMenu:
Preprocessor→Modeling→Creat→Line→lines→straightline→依次连接特征点→Ok
图4.4模型建立
3、划分网格：
ANSYSMainMenu:
Preprocessor→Meshing→Meshtool→Set→选择1，2节点之间部分→Apply→选择2，3节点之间部分→单元长度分别为0.1和0.2→OK
Meshing→Meshtool→Mesh→分别选择1和2，2和3节点之间部分→OK
图4.4划分网格
4、添加约束和载荷：
左下角和右上角添加约束：
ANSYSMainMenu:
Preprocessor→Solution→Defineloads→Apply→Structural→Displacement→Onnodes→选择1节点→ALLDOF→Apply→Onnodes→选择1节点→ALLDOF→OK
添加顶部均布载荷：
ANSYSMainMenu:
Preprocessor→Solution→Defineloads→Apply→Structural→Pressure→Onbeams→选择顶部所有的单元→VALIpressurevaluenodeI:
1000VALJpressurevaluenodeJ:
1000→OK
添加力矩和力：
ANSYSMainMenu:
Preprocessor→Solution→Defineloads→Apply→Structural→Force/Monment→Onnodes→选择2节点→Apply→LABMZ
VALUE100.（输入力矩）→Onnodes→选择8节点→Apply→LABFXVALUE1000（输入力）
图4.5添加约束和载荷
二、分析计算
ANSYSMainMenu:
Solution→Solve→CurrentLS→OK→ShouldtheSolveCommandbeExecuted?
Y→Close（Solutionisdone!
）→关闭窗口
图4.6求解模型
三、求解结果
1、位移
ANSYSMainMenu:
GeneralPostproc→Listresult→Nodalsolution→DOFsolution→X-componentofdisplacement→Apply→Y-componentofdisplacement→OK
图4.7x方向位移解
图4.8y方向位移解
2、支反力：
ANSYSMainMenu:
GeneralPostproc→Listresult→ReactionSolu→Allitems→OK
图4.9支反力结果
四、绘制图像
1、设置参数
ANSYSMainMenu:
GeneralPostproc→ElementTable→DifineTable→Add在userlabelforitem中输入FX-I,在Resultsdataitem中选择Bysequencenum并输入smisc,1→Apply
在userlabelforitem中输入FX-J,在Resultsdataitem中选择Bysequencenum,,并输入smisc,7→Apply
在userlabelforitem中输入FY-I,在Resultsdataitem中选择Bysequencenum,,并输入smisc,2→Apply
在userlabelforitem中输入FY-J,在Resultsdataitem中选择Bysequencenum,,并输入smisc,8→Apply
在userlabelforitem中输入MZ-I,在Resultsdataitem中选择Bysequencenum,,并输入smisc,6→Apply
在userlabelforitem中输入MZ-J,在Resultsdataitem中选择Bysequencenum,,并输入smisc,12→OK
图4.10图表参数设置
2、图像输出
a）轴力图：
ANSYSMainMenu:
GeneralPostproc→Plotresult→Contourplot→LineElemRes→选择FX_IFX_J→Apply
图4.11轴力图
b）剪力图：
ANSYSMainMenu:
GeneralPostproc→Plotresult→Contourplot→LineElemRes→选择FY_IFY_J→Apply
图4.12剪力图
c）弯矩图：
ANSYSMainMenu:
GeneralPostproc→Plotresult→Contourplot→LineElemRes→选择MZ_IZ_J→OK
图4.13弯矩图
第五章结果比较
结果
手算
MATLAB
ANSYS
-0.41422×10-8
-0.41422×10-8
-0.41422×10-8
-0.33023×10-7
-0.33023×10-7
-0.33023×10-7
-1033.1
-1033.1
-1033.1
49.535
49.535
49.535
64.501
64.501
64.501
3.1067
3.1067
3.1067
1950.5
1950.5
1950.5
-1462.3
-1462.3
-1462.3
通过对比知道，三种方式的结果完全一样，显示了结果的正确性。
第六章心得体会
在有限元课程结束之际，通过这份大作业，我整理总结了这个学期学到的有限元思想。
深化我对有限元计算流程的理解，提升了我的matlab的编程水平。
编程时首先要对整体的流程有一个清晰地构想，从数据的获取开始到计算输出进行分步处理；通过手算例题确定计算过程中使用的参数。
再通过模板获得每个单元的矩阵。
难点在于如何将这些单元正确的拼接到一起，一开始用的是biadiag对角拼接指令，但由于每个单元对应的节点并非按顺序排列，不能实现所需效果；之后我重新计算例题，探究到组装的本质是将每个单元的矩阵按照矩阵所对应节点的信息进行分块，并累加到总体矩阵中，由于矩阵实现较为麻烦，所以我想利用元胞数组来存储节点信息，并成功实现了所需功能，又利用cat拼接指令，将元胞数组展开成一个大型矩阵完成拼接。
再获得总体刚度矩阵后，只要找到对应的位移与外力的参数行，再进行乘除计算就相对简单了。
我认为在学习过程中，提高编程能力是很重要的一个方面，有着良好的编程能力，可以让很多工程问题得以用计算机解决，更容易获得结果。
而想要提高编程能力，首先要有一定的编程思想与数学建模能力。
我认为编程思想指的是，要对所做的事情有一个过程性与结构性的认识，并根据使用对象进行相关的调整。
过程性指的是在编程要有一个完整的流程图，将复杂问题转化成若干简单的小问题，针对小问题进行求解；结构性指的是在编程中要利用条件if，循环for等结构简化程序，同时也要对每部分程序的输入、输出以及执行作用有一个明确的认知；最后根据使用对象，对输入、输出，进行调整，也可以用GUI编制界面方便使用。
以上是我对这次大作业的心得体会，希望老师辅导校正
第七章附录
一、matlab程序
clc
clear
formatcompact
formatshortG
jd=input（'请输入节点数：
'）;
dy=input（'请输入单元数：
'）;
E=input（'请输入氏模量E：
'）;
I=input（'请输入惯性矩I：
'）;
L=input（'请输入单元长度L：
'）;
A=input（'请输入单元截面积：
'）;
FAI=input（'请输入单元相对旋转角度：
'）;
%输入对应关系时，小节点放前面[单元节点1节点2]
dy_jd=input（'请输入单元与节点对应关系：
'）;
%输入力与扭矩约束[值作用节点作用类型]（转矩为3x方向为1y方向为2）
lys=input（'力与转矩约束矩阵：
'）;
%输入结构约束[作用节点作用类型]（转角为3x方向为1y方向为2）
wys=input（'结构约束矩阵：
'）;
%原始数据
%L=1;
%E=3*10^10;
%P=1000;
%A=0.05;
%dy=2;jd=3;LL=[L2*L];I=20*A;
%dy_jd=[112;223];
%FAI=[pi/20];
%q=P/L;M=P*L/10;
%lys=[44/125*P11;-12*P*L/12513;81/125*P21;-P22;-67/750*P*L23;-P32;P*L/333];
%wys=[11;12;13;31;32;33];
%对力约束与位移约束式子分别进行编号处理
wys（:
3）=（wys（:
1）-1）*3+wys（:
2）;
lys（:
4）=（lys（:
2）-1）*3+lys（:
3）;
%对力约束与位移约束式子进行排序
lys=sortrows（lys,4）;
wys=sortrows（wys,3）;
%单元刚度矩阵
symsfaieailreal
k=[e*a/l00-e*a/l00;
012*e*i/l^36*e*i/l^20-12*e*i/l^36*e*i/l^2;
06*e*i/l^24*e*i/l0-6*e*i/l^22*e*i/l;
-e*a/l00e*a/l00;
0-12*e*i/l^3-6*e*i/l^2012*e*i/l^3-6*e*i/l^2;
06*e*i/l^22*e*i/l0-6*e*i/l^24*e*i/l];
t=[cos（fai）,sin（fai）,0;
-sin（fai）,cos（fai）,0;
0,0,1];
%坐标变换矩阵
T=blkdiag（t,t）;
%总体坐标系下的单元刚度矩阵
K=T'*k*T;
%带入每个单元的数，生成单元刚度矩阵kk，其每一页对应相应页数的单元的刚度矩阵
forj=1:
dy;
e=E;
i=I;
l=LL（j）;
a=A;
fai=FAI（j）;
kk（:
:
j）=eval（K）;
end
%生成总体刚度矩阵KK
%采用元胞数组的方式对各项进行保存
%生成空元胞数组，元胞数组的行列大小与节点数相同
forj=1:
jd;
forjj=1:
jd;
ling1{j,jj}=zeros（3）;
end
end
ling2=ling1;
%将对单元刚度矩阵部分分成4分加入元胞数组中
forj=1:
dy;
kk1=kk（1:
3,1:
3,j）;
kk2=kk（1:
3,4:
6,j）;
kk3=kk（4:
6,1:
3,j）;
kk4=kk（4:
6,4:
6,j）;
ling2{dy_jd（j,2）,dy_jd（j,2）}=kk1+ling2{dy_jd（j,2）,dy_jd（j,2）};
ling2{dy_jd（j,2）,dy_jd（j,3）}=kk2+ling2{dy_jd（j,2）,dy_jd（j,3）};
ling2{dy_jd（j,3）,dy_jd（j,2）}=kk3+ling2{dy_jd（j,3）,dy_jd（j,2）};
ling2{dy_jd（j,3）,dy_jd（j,3）}=kk4+ling2{dy_jd（j,3）,dy_jd（j,3）};
end
%将元胞数组进行拼接，形成总体刚度矩阵
forj=1:
jd;
ling3（:
:
j）=cat（2,ling2{j,:
}）;
end
KK=ling3（:
:
1）;
forj=2:
jd;
KK=[KK;ling3（:
:
j）];
end
%消去有已知位移的行与列
b=KK;
b（:
wys（:
3））=[];
b（wys（:
3）,:
）=[];
kjiejuzhen=inv（b）;
%提取对应外力
lyss=lys;
forj=1:
size（wys,1）;
forjj=1:
size（lys,1）;
iflyss（jj,4）==wys（j,3）;
lyss（jj,:
）=0;
end
ifjj==size（lyss,1）;
break
end
end
end
lyss（all（lyss==0,2）,:
）=[];
%求解weiyijie=[作用值作用节点作用类型（转角为3x方向为1y方向为2）序列]
weiyijie=kjiejuzhen*lyss（:
1）;
weiyijie（:
1）=weiyijie;
weiyijie（:
2）=lyss（:
2）;
weiyijie（:
3）=lyss（:
3）;
weiyijie（:
4）=lyss（:
4）;
%计算不计作用在约束方向上时的支反力lysjiee=[作用值作用节点作用类型（转角为3x方向为1y方向为2）序列]
lysjie（:
1）=KK（wys（:
3）,lyss（:
4））*weiyijie（:
1）;
lysjie（:
2:
4）=wys（:
1:
3）;
%将作用在约束方向上时的支反力加在上面的求解结果上
forj=1:
size（lysjie,1）
forjj=1:
size（lys,1）;
iflysjie（j,4）==lys（jj,4）;
lysjie（j,1）=lysjie（j,1）-lys（jj,1）;
end
end
end
%答案
weiyijie
lysjie
升级会员 