单元刚度矩阵等参元MATLAB编程.docx
《单元刚度矩阵等参元MATLAB编程.docx》由会员分享，可在线阅读，更多相关《单元刚度矩阵等参元MATLAB编程.docx（23页珍藏版）》请在冰豆网上搜索。
单元刚度矩阵等参元MATLAB编程
《有限元法》实验报告
专业班级力学（实验）1601
姓名田诗豪
学号1603020210
提交日期2019.4.24
实验编号
实验一
实验二
实验三
总分
得分
实验一（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）
%生成平面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=0.5*det（xyh）;
A=abs（A）;
D=mat
（1）/（1-mat
（2）^2）*[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,0.5,1.0];
%****输入材料参数矩阵（弹性模量，泊松比，壁厚）****
[B3,S3,K3]=ele_mat_tri3（xy3,mat）
三、算例分析
●算例1：
如图1所示三角形单元，结点坐标为1（0,0），2（5,2），3（1,4），弹性模量为200GPa，泊松比为0.35、厚度为0.5m。
试求应变矩阵，应力矩阵和单元刚度矩阵。
图1算例1三角形单元
解：
根据如图1所示三角形单元及其几何和材料参数，编制主程序如下：
clear;clc;
%*********输入结点坐标数组********
xy3=[0,0;5,2;1,4];
mat=[2e11,0.35,0.5];
%****输入材料参数矩阵（弹性模量，泊松比，壁厚）****
[B3,S3,K3]=ele_mat_tri3（xy3,mat）
运行程序，得到应变矩阵B3如下：
-0.1111
0.0000
0.2222
0.0000
-0.1111
0.0000
0.0000
-0.2222
0.0000
-0.0556
0.0000
0.2778
-0.2222
-0.1111
-0.0556
0.2222
0.2778
-0.1111
得到应力矩阵S3（Pa）如下：
-2.53E+10
-1.77E+10
5.06E+10
-4.43E+09
-2.53E+10
2.22E+10
-8.86E+09
-5.06E+10
1.77E+10
-1.27E+10
-8.86E+09
6.33E+10
-1.65E+10
-8.23E+09
-4.12E+09
1.65E+10
2.06E+10
-8.23E+09
得到单元刚度矩阵K3（Pa）如下：
2.91E+10
1.71E+10
-2.12E+10
-1.42E+10
-7.91E+09
-2.85E+09
1.71E+10
5.48E+10
-1.57E+10
4.43E+09
-1.42E+09
-5.92E+10
-2.12E+10
-1.57E+10
5.17E+10
-8.55E+09
-3.05E+10
2.42E+10
-1.42E+10
4.43E+09
-8.55E+09
1.96E+10
2.28E+10
-2.41E+10
-7.91E+09
-1.42E+09
-3.05E+10
2.28E+10
3.84E+10
-2.14E+10
-2.85E+09
-5.92E+10
2.42E+10
-2.41E+10
-2.14E+10
8.33E+10
●算例2：
如图2所示三角形单元，结点坐标为1（0,0），2（3,0），3（0,5），弹性模量为200GPa，泊松比为0.35、厚度为0.5m。
试求应变矩阵，应力矩阵和单元刚度矩阵。
图2算例2三角形单元
解：
根据如图2所示三角形单元及其几何和材料参数，编制主程序如下：
clear;clc;
%*********输入结点坐标数组********
xy3=[0,0;3,0;5,0];
mat=[2e11,0.35,0.5];
%****输入材料参数矩阵（弹性模量，泊松比，壁厚）****
[B3,S3,K3]=ele_mat_tri3（xy3,mat）
运行程序，得到应变矩阵B3如下：
-0.3333
0.0000
0.3333
0.0000
0.0000
0.0000
0.0000
-0.2000
0.0000
0.0000
0.0000
0.2000
-0.2000
-0.3333
0.0000
0.3333
0.2000
0.0000
得到应力矩阵S3（Pa）如下：
-7.60E+10
-1.60E+10
7.60E+10
0.00E+00
0.00E+00
1.60E+10
-2.66E+10
-4.56E+10
2.66E+10
0.00E+00
0.00E+00
4.56E+10
-1.48E+10
-2.47E+10
0.00E+00
2.47E+10
1.48E+10
0.00E+00
得到单元刚度矩阵K3（Pa）如下：
1.06E+11
3.85E+10
-9.50E+10
-1.85E+10
-1.11E+10
-1.99E+10
3.85E+10
6.51E+10
-1.99E+10
-3.09E+10
-1.85E+10
-3.42E+10
-9.50E+10
-1.99E+10
9.50E+10
0.00E+00
0.00E+00
1.99E+10
-1.85E+10
-3.09E+10
0.00E+00
3.09E+10
1.85E+10
0.00E+00
-1.11E+10
-1.85E+10
0.00E+00
1.85E+10
1.11E+10
0.00E+00
-1.99E+10
-3.42E+10
1.99E+10
0.00E+00
0.00E+00
3.42E+10
●算例3：
如图3所示三角形单元，结点坐标为1（0,0），2（3,0），3（1.5,1.5
），弹性模量为200GPa，泊松比为0.35、厚度为0.5m。
试求应变矩阵，应力矩阵和单元刚度矩阵。
图3算例3三角形单元
解：
根据如图3所示三角形单元及其几何和材料参数，编制主程序如下：
clear;clc;
%*********输入结点坐标数组********
xy3=[0,0;3,0;1.5,1.5*sqrt（3）];
mat=[2e11,0.35,0.5];
%****输入材料参数矩阵（弹性模量，泊松比，壁厚）****
[B3,S3,K3]=ele_mat_tri3（xy3,mat）
运行程序，得到应变矩阵B3如下：
-0.3333
0.0000
0.3333
0.0000
0.0000
0.0000
0.0000
-0.1925
0.0000
-0.1925
0.0000
0.3849
-0.1925
-0.3333
-0.1925
0.3333
0.3849
0.0000
得到应力矩阵S3（Pa）如下：
-7.60E+10
-1.54E+10
7.60E+10
-1.54E+10
0.00E+00
3.07E+10
-2.66E+10
-4.39E+10
2.66E+10
-4.39E+10
0.00E+00
8.77E+10
-1.43E+10
-2.47E+10
-1.43E+10
2.47E+10
2.85E+10
0.00E+00
得到单元刚度矩阵K3（Pa）如下：
5.47E+10
1.92E+10
-4.40E+10
7.12E+08
-1.07E+10
-1.99E+10
1.92E+10
3.25E+10
-7.12E+08
4.11E+08
-1.85E+10
-3.29E+10
-4.40E+10
-7.12E+08
5.47E+10
-1.92E+10
-1.07E+10
1.99E+10
7.12E+08
4.11E+08
-1.92E+10
3.25E+10
1.85E+10
-3.29E+10
-1.07E+10
-1.85E+10
-1.07E+10
1.85E+10
2.14E+10
0.00E+00
-1.99E+10
-3.29E+10
1.99E+10
-3.29E+10
0.00E+00
6.58E+10
实验二（30分）
一、实验内容
编写一个计算平面4结点四边形等参元的刚度矩阵的MATLAB函数文件K4=ele_mat_quad4（xy4,mat），其中：
输入变量xy4为结点坐标数组，mat为材料参数矩阵；输出变量K4为单元刚度矩阵。
（要求给出3个不同算例进行验证，并绘制出单元形状和结点号）
二、程序代码
●通用函数
functionK4=ele_mat_quad4（xy4,mat）
%生成平面4结点四边形等参元的刚度矩阵的功能函数
%*********变量说明****************
%xy4------------------结点坐标数组
%mat------------------材料参数矩阵（弹性模量，泊松比，壁厚）
%K4-------------------单元刚度矩阵
%y1y2-----------------局部坐标系结点坐标
%D--------------------弹性矩阵
%*********************************
y1y2=[-1,-1;1,-1;1,1;-1,1];
D=mat
（1）/（1-mat
（2）^2）*[1,mat
（2）,0;mat
（2）,1,0;0,0,（1-mat
（2））/2];
%*********数值积分（Guass,n=4）****
C
（1）=0.861136311594055;C
（2）=0.339981043584886;
C（3）=-0.339981043584886;C（4）=-0.861136311594055;
A
（1）=0.347854845137454;A
（2）=0.652145154862546;
A（3）=0.652145154862546;A（4）=0.347854845137454;
sum=0;
fori=1:
4
forj=1:
4
y1=C;y2=C;
k=1:
4;
%*********************************
PN1
（1）=0.25*（y2（j）-1）;PN2
（1）=0.25*（y1（i）-1）;
PN1
（2）=-0.25*（y2（j）-1）;PN2
（2）=-0.25*（y1（i）+1）;
PN1（3）=0.25*（y2（j）+1）;PN2（3）=0.25*（y1（i）+1）;
PN1（4）=-0.25*（y2（j）+1）;PN2（4）=-0.25*（y1（i）-1）;
fork=1:
4
PN（:
k）=[PN1（k）,PN2（k）]';
end
J=PN*xy4;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）];
B=[B1,B2,B3,B4];
%*********************************
m=mat（3）*B'*D*B*det（J）;
sum=sum+A（i）*A（j）*m;
end
end
K4=vpa（sum,9）;
●主程序
clc;clear;
%*********输入结点坐标数组********
xy4=[3,2;8,3;7,8;4,7];
%****输入材料参数矩阵（弹性模量，泊松比，壁厚）****
mat=[3e6,0.5,1.0];
K4=ele_mat_quad4（xy4,mat）
三、算例分析
●算例1：
如图5所示四边形等参单元（图4为局部坐标系规则单元），已知4个结点整体坐标系内的坐标为1（3,2）;2（8,3）;3（7,8）;4（4,7）。
弹性模量为200GPa，泊松比为0.35、厚度为0.05m。
试求单元刚度矩阵。
图4局部坐标系规则单元
图5算例1四边形等参单元
解：
根据如图5所示四边形等参单元及其几何和材料参数，编制主程序如下：
clc;clear;
%*********输入结点坐标数组********
xy4=[3,2;8,3;7,8;4,7];
%****输入材料参数矩阵（弹性模量，泊松比，壁厚）****
mat=[2e11,0.35,0.5];
K4=ele_mat_quad4（xy4,mat）
运行程序，得到单元刚度矩阵K4（Pa）如下：
4.01E+09
1.45E+09
-3.54E+09
-2.37E+08
-1.54E+09
-1.67E+09
1.08E+09
4.56E+08
1.45E+09
3.79E+09
-3.79E+08
4.98E+08
-1.67E+09
-1.87E+09
5.98E+08
-2.41E+09
-3.54E+09
-3.79E+08
6.62E+09
-2.37E+09
1.38E+09
5.71E+08
-4.47E+09
2.18E+09
-2.37E+08
4.98E+08
-2.37E+09
4.53E+09
4.28E+08
-2.17E+09
2.18E+09
-2.86E+09
-1.54E+09
-1.67E+09
1.38E+09
4.28E+08
5.24E+09
1.34E+09
-5.08E+09
-9.97E+07
-1.67E+09
-1.87E+09
5.71E+08
-2.17E+09
1.34E+09
4.74E+09
-2.42E+08
-6.98E+08
1.08E+09
5.98E+08
-4.47E+09
2.18E+09
-5.08E+09
-2.42E+08
8.47E+09
-2.54E+09
4.56E+08
-2.41E+09
2.18E+09
-2.86E+09
-9.97E+07
-6.98E+08
-2.54E+09
5.97E+09
●算例2：
如图6所示四边形等参单元，已知4个结点整体坐标系内的坐标为1（1,1）;2（4,1）;3（4,3）;4（1,4）。
弹性模量为180GPa，泊松比为0.5、厚度为0.1m。
试求单元刚度矩阵。
图6算例2四边形等参单元
解：
根据如图6所示四边形等参单元及其几何和材料参数，编制主程序如下：
clc;clear;
%*********输入结点坐标数组********
xy4=[1,1;4,1;4,3;1,3];
%****输入材料参数矩阵（弹性模量，泊松比，壁厚）****
mat=[1.8e11,0.5,0.1];
K4=ele_mat_quad4（xy4,mat）
运行程序，得到单元刚度矩阵K4（Pa）如下：
8.33E+09
4.50E+09
-3.83E+09
1.50E+09
-4.17E+09
-4.50E+09
-3.33E+08
-1.50E+09
4.50E+09
1.33E+10
-1.50E+09
4.67E+09
-4.50E+09
-6.67E+09
1.50E+09
-1.13E+10
-3.83E+09
-1.50E+09
8.33E+09
-4.50E+09
-3.33E+08
1.50E+09
-4.17E+09
4.50E+09
1.50E+09
4.67E+09
-4.50E+09
1.33E+10
-1.50E+09
-1.13E+10
4.50E+09
-6.67E+09
-4.17E+09
-4.50E+09
-3.33E+08
-1.50E+09
8.33E+09
4.50E+09
-3.83E+09
1.50E+09
-4.50E+09
-6.67E+09
1.50E+09
-1.13E+10
4.50E+09
1.33E+10
-1.50E+09
4.67E+09
-3.33E+08
1.50E+09
-4.17E+09
4.50E+09
-3.83E+09
-1.50E+09
8.33E+09
-4.50E+09
-1.50E+09
-1.13E+10
4.50E+09
-6.67E+09
1.50E+09
4.67E+09
-4.50E+09
1.33E+10
●算例3：
如图7所示四边形等参单元，已知4个结点整体坐标系内的坐标为1（1,1）;2（4,1）;3（5,3）;4（2,3）。
弹性模量为180GPa，泊松比为0.5、厚度为0.1m。
试求单元刚度矩阵。
图7算例3四边形等参单元
解：
根据如图7所示四边形等参单元及其几何和材料参数，编制主程序如下：
clc;clear;
%*********输入结点坐标数组********
xy4=[1,1;4,1;5,3;2,3];
%****输入材料参数矩阵（弹性模量，泊松比，壁厚）****
mat=[1.8e11,0.5,0.1];
K4=ele_mat_quad4（xy4,mat）
运行程序，得到单元刚度矩阵K4（Pa）如下：
7.17E+09
2.50E+09
-4.17E+09
3.50E+09
-2.83E+09
-3.50E+09
-1.67E+08
-2.50E+09
2.50E+09
8.67E+09
5.00E+08
3.33E+09
-3.50E+09
-1.33E+09
5.00E+08
-1.07E+10
-4.17E+09
5.00E+08
1.02E+10
-6.50E+09
-1.67E+08
5.00E+08
-5.83E+09
5.50E+09
3.50E+09
3.33E+09
-6.50E+09
2.07E+10
-2.50E+09
-1.07E+10
5.50E+09
-1.33E+10
-2.83E+09
-3.50E+09
-1.67E+08
-2.50E+09
7.17E+09
2.50E+09
-4.17E+09
3.50E+09
-3.50E+09
-1.33E+09
5.00E+08
-1.07E+10
2.50E+09
8.67E+09
5.00E+08
3.33E+09
-1.67E+08
5.00E+08
-5.83E+09
5.50E+09
-4.17E+09
5.00E+08
1.02E+10
-6.50E+09
-2.50E+09
-1.07E+10
5.50E+09
-1.33E+10
3.50E+09
3.33E+09
-6.50E+09
2.07E+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------------------材料参数矩阵（弹性模量
升级会员 