单元刚度矩阵MATLAB编程.docx
《单元刚度矩阵MATLAB编程.docx》由会员分享,可在线阅读,更多相关《单元刚度矩阵MATLAB编程.docx(16页珍藏版)》请在冰豆网上搜索。
![单元刚度矩阵MATLAB编程.docx](https://file1.bdocx.com/fileroot1/2023-4/16/4aa87854-09f2-49bb-bbb3-135f05642633/4aa87854-09f2-49bb-bbb3-135f056426331.gif)
单元刚度矩阵MATLAB编程
《有限元法》实验报告
专业班级力学(实验)1601
姓名田诗豪
学号
提交日期
实验编号
实验一
实验二
实验三
总分
得分
实验一(30分)
一、实验内容
编写一个计算平面3结点三角形单元的应变矩阵、应力矩阵和单元刚度矩阵的MATLAB
函数文件[B3,S3,K3]=ele_mat_tri3(xy3,mat),其中:
输入变量xy3为结点坐标数组,mat为材料参数矩阵;输出变量B3为应变矩阵,S3为应力矩阵,K3为单元刚度矩阵。
(要求给
出3个不同算例进行验证,并绘制出单元形状和结点号)
、程序代码
通用函数
function[B3,S3,K3]=ele_mat_tri3(xy3,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);
°%*********************************
fori=1:
3
ifi==1
j=2;
m=3;
elseifi==2
j=3;
m=1;
else
j=1;
m=2;
end
b(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,,];
%****输入材料参数矩阵(弹性模量,泊松比,壁厚)****
[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
+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所示三角形单元及其几何和材料参数,编制主程序如下:
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
+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
+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为材料参数矩阵;输出变量K4为单元刚度矩阵。
(要求给出3个不同算例进行验证,并绘制出单元形状和结点号)
二、程序代码
通用函数
functionK4=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泊松比
为、厚度为。
试求单元刚度矩阵。
图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
+09
+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_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)。
弹性模量为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
+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、程序代码
通用函数
functionK8=ele_mat_quad8(xy8,mat)
毗成平面8结点四边形等参元刚度矩阵的功能函数
%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];symsy1y2real
N
(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)*(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);
fori=1:
8
PN(:
i)=[PN1(i),PN2(i)]';end
J=PN*xy8;JN=inv(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*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(:
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;
fori=1:
4
forj=1:
4
sum=sum+A(i)*A(j)*subs(m,[y1,y2],[C(i),C(j)]);end
end
K8=vpa(sum,9);
主程序
clc;clear;
xy8=[2,2;7,2;6,7;3,6;5,3;6,4;4,6;3,4];
%****输入材料参数矩阵(弹性模量,泊松比,壁厚)****
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,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];
%****输入材料参数矩阵(弹性模量,泊松比,壁厚)****
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排版的能力,最重要的是加深了对有限单元法的理解,加深了自己的力学功底。
本次实验的所有算例,不仅都通过我的程序得出结果,还经过宋戎妆同学的程序验证,故可以保证程序与算例的正确性。
在这里衷心地感谢周博教授的指导,同时感谢几位帮助过我的同学。