matlab-四节点四边形等参元的刚度矩阵计算程序.doc

课程（论文类）试卷2 011 /2 012 学年第 1 学期课程名称： 弹塑性力学中的有限元分析 课程代码： 论文题目： 平面四节点四边形等参元的刚度矩阵程序 学生姓名： 专业学号： 学院： 机械工程学院 课程（论文）成绩：课程（论文）评分依据（必填）：任课教师签字： 日期： 年 月 日2D梁单元以及2D框架结构计算研究报告姓名： 学号： 日期：2011年11月11日1、 课题要求在matlab平台上编写四节点四边形等参元的刚度矩阵计算程序，要求输入单元的节点坐标及材料参数，输出单元刚度矩阵。2、 单元划分（图）及输入数据列表（1）单元划分：将题目所给的坐标转到自然坐标下（2）输入数据列表：单元节点1、2、3、4的坐标：x1=1;y1=0;x2=2;y2=0;x3=2.25;y3=1.5;x4=1.25;y4=1;材料参数：E=30e12;NU=0.3;h=1;ID=1;Ai=1Aj=1;3、 程序（由附件给出）4、 计算结果及讨论（需有图表来说明）计算结果如下：k = 1.0e+012 *1.4619 0.2211 -1.2150 0.0616 -0.3716 -0.4236 0.1248 0.1409 0.2211 1.1491 0.1440 -0.0419 -0.4236 -0.2582 0.0585 -0.8489 -1.2150 0.1440 2.1968 -0.8933 0.5212 -0.0151 -1.5030 0.7645 0.0616 -0.0419 -0.8933 1.8399 0.0673 -0.5250 0.7645 -1.2729 -0.3716 -0.4236 0.5212 0.0673 1.1645 0.2763 -1.3141 0.0800 -0.4236 -0.2582 -0.0151 -0.5250 0.2763 0.9061 0.1624 -0.1229 0.1248 0.0585 -1.5030 0.7645 -1.3141 0.1624 2.6923 -0.9854 0.1409 -0.8489 0.7645 -1.2729 0.0800 -0.1229 -0.9854 2.2448讨论：利用等参元编写程序的计算结果与参考文献提供的例题对比。参考文献计算单元刚度矩阵采用int函数直接积分，计算结果是精确解。本题单元刚度矩阵的计算采用高斯积分，根据对比发现，两者的计算结果在小数点后四位没有差别，因此，在对精度要求不高的情况下高斯积分满足要求。而且，由于高斯积分的方法是将积分运算简化为了加法运算，使得运算变快，节省了计算时间。需要注意的是，编写MATLAB程序时要尽量减少符号变量的使用，以减少计算机的运算时间。5、 参考文献P.L.Kattan著.MATLAB有限元分析与应用.韩来彬译.北京:清华大学出版社，20046、 附录%Quad2D4NodeGuass_Stiffness%begin%function k=Quad2D4NodeGuass_Stiffness(E,NU,h,x1,y1,x2,y2,x3,y3,x4,y4,ID)%该程序计算平面四节点四边形等参元的刚度矩阵%输入弹性模量E、泊松比NU和厚度h%输入四个节点1、2、3、4的坐标x1,x2,x3,x3,y1,y2,y3,y4%输入平面问题性质参数ID（1为平面应力，2为平面应变）%输出单元刚度矩阵%-syms s t;a=-(1-t)*x1+(1-t)*x2+(1+t)*x3-(1+t)*x4/4;b=-(1-s)*x1-(1+s)*x2+(1+s)*x3+(1-s)*x4/4;c=-(1-t)*y1+(1-t)*y2+(1+t)*y3-(1+t)*y4/4;d=-(1-s)*y1-(1+s)*y2+(1+s)*y3+(1-s)*y4/4;J=a c;b d;L=det(J);H=d -c 0 0 ; 0 0 -b a; -b a d -c/L;Q=t-1 0 1-t 0 1+t 0 -(1+t) 0; s-1 0 -(1+s) 0 1+s 0 1-s 0; 0 t-1 0 1-t 0 1+t 0 -(1+t); 0 s-1 0 -(1+s) 0 1+s 0 1-s/4;if ID=1 D=(E/(1-NU*NU)*1 NU 0;NU 1 0;0 0 (1-NU)/2;else if ID=2 D=(E/1+NU)/(1-2*NU)*1-NU NU 0;NU 1-NU 0;0 0 (1-2*NU)/2; endendT=h*Ai*Aj*(Q*H*D*H*Q)*L;T1=subs(T,s,t,-1/3*sqrt(3),-1/3*sqrt(3);T2=subs(T,s,t,1/3*sqrt(3),-1/3*sqrt(