单元刚度矩阵MATLAB编程.docx

1、单元刚度矩阵MATLAB编程有限元法实验报告专业班级 力学（实验）1601姓 名 田诗豪 学 号 提交日期 实验编号实验一实验二实验三总分得分实验一（30分）一、实验内容编写一个计算平面 3结点三角形单元的应变矩阵、应力矩阵和单元刚度矩阵的 MATLAB函数文件B3,S3,K3 = ele_mat_tri3（xy3,mat） ，其中：输入变量xy3为结点坐标数组，mat 为材料参数矩阵；输出变量 B3为应变矩阵，S3为应力矩阵，K3为单元刚度矩阵。（要求给出3个不同算例进行验证，并绘制出单元形状和结点号）、程序代码通用函数fun ction B3,S3,K3 = ele_mat_tri3（xy

2、3,mat）%fe成平面3结点三角形单元的应变矩阵、应力矩阵和单元刚度矩阵的功能函数%xy3 结点坐标数组%mat 材料参数矩阵（弹性模量，泊松比，壁厚）%B3 应变矩阵%S3 应力矩阵%K3 单元刚度矩阵%* xyh=1,xy3(1,1),xy3(1,2);1,xy3(2,1),xy3(2,2);1,xy3(3,1),xy3(3,2);A=*det(xyh);A=abs(A);D=mat(1)/(1-mat (2)A2)*1,mat (2) ,0;mat (2) ,1,0;0,0,(1-mat (2)/2;b=zeros(1,3);c=zeros(1,3);%*for i=1:3if i=1

3、j=2;m=3;elseif i=2j=3;m=1;elsej=1;m=2;endb(i)=xy3(j,2)-xy3(m,2);c(i)=xy3(m,1)-xy3(j,1);end%*B31=1/(2*A)*b(1),0;0,c(1);c(1),b(1);B32=1/(2*A)*b(2),0;0,c(2);c(2),b(2);B33=1/(2*A)*b(3),0;0,c(3);c(3),b(3);B3=B31,B32,B33;%*S3=D*B3;%*K3=A*mat (3)*B3*D*B3;主程序clear;clc;xy3=0,0;5,1;1,4;mat=3e6,;%*输入材料参数矩阵(弹性模

4、量，泊松比，壁厚)*B3,S3,K3=ele_mat_tri3(xy3,mat)三、算例分析算例1:如图1所示三角形单元，结点坐标为 1(0,0) ，2(5,2)，3(1,4)，弹性模量为200GPa泊松比为、厚度为。试求应变矩阵，应力矩阵和单元刚度矩阵。图1算例1三角形单元解：根据如图1所示三角形单元及其几何和材料参数，编制主程序如下:xy3=0,0;5,2;1,4;mat=2e11,;%*输入材料参数矩阵(弹性模量，泊松比，壁厚)*B3,S3,K3=ele_mat_tri3(xy3,mat)运行程序，得到应变矩阵 B3如下:得到应力矩阵S3( Pa)如下:+10+ 10+ 10+09+10

5、+10+09+ 10+ 10+ 10+09+10+ 10+09+09+ 10+10+09得到单元刚度矩阵 K3 ( Pa)如下:+10+ 10+10+ 10+09+09+10+ 10+10+09+09+10+10+ 10+10+09+ 10+10+10+09+09+ 10+ 10+10+09+09+10+ 10+ 10+10+09+ 10+10+ 10+ 10+10算例2 :如图2所示三角形单元，结点坐标为 1(0,0)，2(3,0)，3(0,5)，弹性模量为200GPa泊松比为、厚度为。试求应变矩阵，应力矩阵和单元刚度矩阵。图2算例2三角形单元解：根据如图2所示三角形单元及其几何和材料参数，

6、编制主程序如下:xy3=0,0;3,0;5,0;mat=2e11,;%*输入材料参数矩阵(弹性模量，泊松比，壁厚)*B3,S3,K3=ele_mat_tri3(xy3,mat)运行程序，得到应变矩阵 B3如下:得到应力矩阵S3( Pa)如下:+ 10+ 10+ 10+00+00+10+ 10+ 10+ 10+00+00+10+ 10+ 10+00+ 10+10+00得到单元刚度矩阵K3( Pa)如下：+11+10+ 10+ 10+10+10+ 10+10+ 10+ 10+10+10+ 10+10+ 10+00+00+10+ 10+ 10+00+ 10+10+00+ 10+ 10+00+ 10+

7、10+00+ 10+ 10+ 10+00+00+10算例3:如图3所示三角形单元，结点坐标为1(0,0) ,2(3,0) ,3, 丁 3,弹性模量为200GPa 泊松比为、厚度为。试求应变矩阵，应力矩阵和单元刚度矩阵。图3算例3三角形单元解：根据如图3所示三角形单元及其几何和材料参数，编制主程序如下:clear;clc;xy3=0,0;3,0;,*sqrt(3);mat=2e11,;%*输入材料参数矩阵(弹性模量，泊松比，壁厚)*B3,S3,K3=ele_mat_tri3(xy3,mat)运行程序，得到应变矩阵 B3如下:得到应力矩阵S3（ Pa）如下:+ 10+ 10+ 10+ 10+00+

8、10+ 10+ 10+ 10+ 10+00+10+ 10+ 10+ 10+ 10+10+00得到单元刚度矩阵K3 (Pa)如下：+ 10+10+ 10+08+10+10+ 10+ 10+08+08+10+10+ 10+08+ 10+ 10+10+10+08+08+ 10+ 10+10+10+ 10+ 10+ 10+ 10+10+00+ 10+ 10+ 10+ 10+00+10实验二（30分）一、 实验内容编写一个计算平面4结点四边形等参元的刚度矩阵的 MATLAB数文件K4 = ele_mat_quad4（xy4,mat），其中：输入变量xy4为结点坐标数组，mat为材料参 数矩阵；输出变量K

9、4为单元刚度矩阵。（要求给出3个不同算例进行验证，并绘 制出单元形状和结点号）二、 程序代码通用函数function K4 = ele_mat_quad4(xy4,mat)毗成平面4结点四边形等参元的刚度矩阵的功能函数主程序xy4=3,2;8,3;7,8;4,7;%*输入材料参数矩阵(弹性模量，泊松比，壁厚)*mat=3e6,;K4=ele_mat_quad4(xy4,mat)三、算例分析算例1 :如图5所示四边形等参单元(图 4为局部坐标系规则单元)，已知4个结点整体坐标系内的坐标为 1(3,2); 2(8,3); 3(7,8); 4(4,7) 。弹性模量为200GPa泊松比为、厚度为。试求

10、单元刚度矩阵。图4局部坐标系规则单元图5算例1四边形等参单元解：根据如图5所示四边形等参单元及其几何和材料参数，编制主程序如下:clc;clear;%*输入结点坐标数组*xy4=3,2;8,3;7,8;4,7;%*输入材料参数矩阵(弹性模量，泊松比，壁厚)* mat=2e11,;K4=ele_mat_quad4(xy4,mat)运行程序，得到单元刚度矩阵K4( Pa)如下+08+09+09+09+08+09+09+09+09+09+08+08+09+09+08+09+09+08+09+09+09+08+09+09+08+08+09+09+08+09+09+09+09+09+09+08+09+0

11、9+09+07+09+09+08+09+09+09+08+08+09+08+09+09+09+08+09+09+08+09+09+09+07+08+09+09算例2 :如图6所示四边形等参单元，已知4个结点整体坐标系内的坐标为1(1,1);2(4,1); 3(4,3); 4(1,4) 。弹性模量为 180GPa泊松比为、厚度为。试求单元刚度矩阵。图6算例2四边形等参单元解：根据如图6所示四边形等参单元及其几何和材料参数，编制主程序如下:clc;clear;%*输入结点坐标数组*xy4=1,1;4,1;4,3;1,3;%*输入材料参数矩阵(弹性模量，泊松比，壁厚 广* mat=,;K4=ele_

12、mat_quad4(xy4,mat)运行程序，得到单元刚度矩阵 K4( Pa)如下:+09+09+09+09+09+09+08+09+09+10+09+09+09+09+09+10+09+09+09+09+08+09+09+09+09+09+09+10+09+10+09+09+09+09+08+09+09+09+09+09+09+09+09+10+09+10+09+09+08+09+09+09+09+09+09+09+09+10+09+09+09+09+09+10算例3:如图7所示四边形等参单兀，已知4个结点整体坐标系内的坐标为1(1,1);2(4,1); 3(5,3); 4(2,3) 。弹性

13、模量为 180GPa泊松比为、厚度为。试求单元刚度矩阵。图7算例3四边形等参单元解：根据如图7所示四边形等参单元及其几何和材料参数，编制主程序如下:clc;clear;%*输入结点坐标数组*xy4=1,1;4,1;5,3;2,3;%*输入材料参数矩阵(弹性模量，泊松比，壁厚)* mat=,;K4=ele_mat_quad4(xy4,mat)运行程序，得到单元刚度矩阵K4 (Pa)如下：+09+09+09+09+09+09+08+09+09+09+08+09+09+09+08+10+09+08+10+09+08+08+09+09+09+09+09+10+09+10+09+10+09+09+08+

14、09+09+09+09+09+09+09+08+10+09+09+08+09+08+08+09+09+09+08+10+09+09+10+09+10+09+09+09+10实验三（40分）1、实验内容编写一个计算平面8结点四边形等参元刚度矩阵的 MATLAB函数文件K8 = ele_mat_quad8（xy8,mat），其中：输入变量xy8为结点坐标数组，mat为材料参 数矩阵；输出变量K8为单元刚度矩阵。（要求给出3个不同算例进行验证，并绘 制出单元形状和结点号）2、程序代码通用函数function K8 = ele_mat_quad8（xy8,mat）毗成平面8结点四边形等参元刚度矩阵的功

15、能函数%xy8 结点坐标数组%mat 材料参数矩阵（弹性模量，泊松比，壁厚）%K8 单元刚度矩阵%D 弹性矩阵%*PN = sym(zeros(2,8);D=mat(1)/(1-mat (2)A2)*1,mat (2) ,0;mat (2) ,1,0;0,0,(1-mat (2)/2; syms y1 y2 realN(1)=(1-y1)*(1-y2)*(-1-y1-y2)/4;N(2)=(1+y1)*(1-y2)*(-1-y2+y1)/4;N(3)=(1+y1)*(1+y2)*(-1+y2+y1)/4;N(4)=(1-y1)*(1+y2)*(-1+y2-y1)/4;N(5)=(1-y1A2)

16、*(1-y2)/2;N(6)=(1-y2A2)*(1+y1)/2;N(7)=(1-y1A2)*(1+y2)/2;N(8)=(1-y1)*(1-y2A2)/2;%* i=1:8;PN1(i)=diff(N(i),y1);PN2(i)=diff(N(i),y2);for i=1:8PN(:,i)=PN1(i),PN2(i); endJ=PN*xy8;JN=i nv(J);J1=JN(1,:);J2=JN(2,:);B1=J1*PN(:,1),0;0,J2*PN(:,1);J2*PN(:,1),J1*PN(:,1);B2=J1*PN(:,2),0;0,J2*PN(:,2);J2*PN(:,2),J1

17、*PN(:,2);B3=J1*PN(:,3),0;0,J2*PN(:,3);J2*PN(:,3),J1*PN(:,3);B4=J1*PN(:,4),0;0,J2*PN(:,4);J2*PN(:,4),J1*PN(:,4);B5=J1*PN(:,5),0;0,J2*PN(:,5);J2*PN(:,5),J1*PN(:,5);B6=J1*PN(:,6),0;0,J2*PN(:,6);J2*PN(:,6),J1*PN(:,6);B7=J1*PN(:,7),0;0,J2*PN(:,7);J2*PN(:,7),J1*PN(:,7);B8=J1*PN(:,8),0;0,J2*PN(:,8);J2*PN(:

18、,8),J1*PN(:,8);B=B1,B2,B3,B4,B5,B6,B7,B8;%*m=mat(3)*B*D*B*det(J);%* 数值积分(Guass, n=4)*C(1)=; C(2)=;C(3)=; C(4)=;A(1)=; A(2)=;A(3)=; A(4)=;sum=0;for i = 1:4for j = 1:4sum = sum+A(i)*A(j)*subs(m,y1,y2,C(i),C(j); endendK8=vpa(sum,9);主程序clc;clear;xy8=2,2;7,2;6,7;3,6;5,3;6,4;4,6;3,4;%*输入材料参数矩阵(弹性模量，泊松比，壁厚

19、)*mat=3e6,;K8 = ele_mat_quad8(xy8,mat)3、算例分析算例1 :如图9所示曲四边形等参单元(图 8为局部坐标系规则单元)，已知8个结点 整体坐标系内的坐标为 1(2,2); 2(7,2); 3(6,7); 4(3,6); 5(5,3); 6(6,4); 7(4,6);8(3,4)。弹性模量为100GPa泊松比为、厚度为。试求单元刚度矩阵。图8局部坐标系规则单元图9算例1曲四边形等参单元解：根据如图9所示曲四边形等参单元及其几何和材料参数，编制主程序如下:clc;clear;%*输入结点坐标数组*xy8=2,2;7,2;6,7;3,6;5,3;6,4;4,6;3

20、,4;%*输入材料参数矩阵(弹性模量，泊松比，壁厚)*mat=,;K8 = ele_mat_quad8(xy8,mat)1(-5,-4);弹性模量为算例2:如图10所示曲四边形等参单元，已知8个结点整体坐标系内的坐标为2(6,-12); 3(10,8); 4(-8,12); 5(0,-6); 6(5,-1); 7(0,4); 8(-4,2)200GPa泊松比为、厚度为。试求单元刚度矩阵。图10算例2曲四边形等参单元解：根据如图10所示曲四边形等参单元及其几何和材料参数，编制主程序如下:clc;clear;xy8=-5,-4;6,-12;10,8;-8,12;0,-6;5,-1;0,4;-4,2

21、;%*输入材料参数矩阵(弹性模量，泊松比，壁厚)*mat=2e11,;K8 = ele_mat_quad8(xy8,mat)算例3:如图11所示曲四边形等参单元，已知8个结点整体坐标系内的坐标为 1(-6,-4);2(6,-4); 3(6,4); 4(-6,4); 5(0,-3); 6(3,0); 7(0,3); 8(-3,0)。弹性模量为 200GPa泊松比为、厚度为。试求单元刚度矩阵。图11算例3曲四边形等参单元解：根据如图11所示曲四边形等参单元及其几何和材料参数，编制主程序如下:clc;clear;xy8=-6,-4;6,-4;6,4;-6,4;0,-3;3,0;0,3;-3,0;%*输入材料参数矩阵(弹性模量，泊松比，壁厚)*mat=2e11,;K8 = ele_mat_quad8(xy8,mat)总结与备注这次有限元编程实验给了我很多收获，不仅锻炼了 Matlab 编程的能力，还 加强了 Word排版的能力，最重要的是加深了对有限单元法的理解，加深了自己 的力学功底。本次实验的所有算例， 不仅都通过我的程序得出结果， 还经过宋戎妆同学的 程序验证，故可以保证程序与算例的正确性。在这里衷心地感谢周博教授的指导，同时感谢几位帮助过我的同学。