单元刚度矩阵等参元MATLAB编程.docx

1、单元刚度矩阵等参元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

2、,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=zer

3、os(1,3);%*for i=1:3 if i=1 j=2; m=3; elseif i=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;%*输

4、入结点坐标数组*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; %*输入材料参数矩

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+105.06E+10-4.43E+09-2.53E+102.22E+10-8.86E+09-5.06E+101.77E+10-1.27E+10-8.

6、86E+096.33E+10-1.65E+10-8.23E+09-4.12E+091.65E+102.06E+10-8.23E+09得到单元刚度矩阵K3（Pa）如下：2.91E+101.71E+10-2.12E+10-1.42E+10-7.91E+09-2.85E+091.71E+105.48E+10-1.57E+104.43E+09-1.42E+09-5.92E+10-2.12E+10-1.57E+105.17E+10-8.55E+09-3.05E+102.42E+10-1.42E+104.43E+09-8.55E+091.96E+102.28E+10-2.41E+10-7.91E+09-1

7、.42E+09-3.05E+102.28E+103.84E+10-2.14E+10-2.85E+09-5.92E+102.42E+10-2.41E+10-2.14E+108.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; %*输入材料参数矩阵(弹性模量，泊松

8、比，壁厚)*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+107.60E+100.00E+000.00E+001.60E+10-2.66E+10-4.56E+102.66E+100.00E+000.00E+004.56E+10-1.

9、48E+10-2.47E+100.00E+002.47E+101.48E+100.00E+00得到单元刚度矩阵K3（Pa）如下：1.06E+113.85E+10-9.50E+10-1.85E+10-1.11E+10-1.99E+103.85E+106.51E+10-1.99E+10-3.09E+10-1.85E+10-3.42E+10-9.50E+10-1.99E+109.50E+100.00E+000.00E+001.99E+10-1.85E+10-3.09E+100.00E+003.09E+101.85E+100.00E+00-1.11E+10-1.85E+100.00E+001.85E+

10、101.11E+100.00E+00-1.99E+10-3.42E+101.99E+100.00E+000.00E+003.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,

11、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+107.60E+10-1.54E+100.00E+003.07E+10-2.66E+10-4.39E+102.66E+10-4.39E+100.00E+008.77E+10-1.43E+1

12、0-2.47E+10-1.43E+102.47E+102.85E+100.00E+00得到单元刚度矩阵K3（Pa）如下：5.47E+101.92E+10-4.40E+107.12E+08-1.07E+10-1.99E+101.92E+103.25E+10-7.12E+084.11E+08-1.85E+10-3.29E+10-4.40E+10-7.12E+085.47E+10-1.92E+10-1.07E+101.99E+107.12E+084.11E+08-1.92E+103.25E+101.85E+10-3.29E+10-1.07E+10-1.85E+10-1.07E+101.85E+102

13、.14E+100.00E+00-1.99E+10-3.29E+101.99E+10-3.29E+100.00E+006.58E+10实 验 二（30分）一、实验内容编写一个计算平面4结点四边形等参元的刚度矩阵的MATLAB函数文件K4 = ele_mat_quad4(xy4,mat)，其中：输入变量xy4为结点坐标数组，mat为材料参数矩阵；输出变量K4为单元刚度矩阵。（要求给出3个不同算例进行验证，并绘制出单元形状和结点号）二、程序代码通用函数function K4 = ele_mat_quad4(xy4,mat)%生成平面4结点四边形等参元的刚度矩阵的功能函数%*变量说明*%xy4-结点坐

14、标数组%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.652145154862

15、546;A(3)=0.652145154862546;A(4)=0.347854845137454;sum=0;for i = 1:4 for j = 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); for k=1:4 PN(:,k

16、)=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=m

17、at(3)*B*D*B*det(J); sum = sum+A(i)*A(j)*m; endendK4=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

18、局部坐标系规则单元图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+091.45E+09-3.54E+09-2.37E+08-1.54E+09-1.67E+091.08E+094.56E+081.45E+093.79E+09-3.79E+084.98E+08-1.67E+09-1

19、.87E+095.98E+08-2.41E+09-3.54E+09-3.79E+086.62E+09-2.37E+091.38E+095.71E+08-4.47E+092.18E+09-2.37E+084.98E+08-2.37E+094.53E+094.28E+08-2.17E+092.18E+09-2.86E+09-1.54E+09-1.67E+091.38E+094.28E+085.24E+091.34E+09-5.08E+09-9.97E+07-1.67E+09-1.87E+095.71E+08-2.17E+091.34E+094.74E+09-2.42E+08-6.98E+081.0

20、8E+095.98E+08-4.47E+092.18E+09-5.08E+09-2.42E+088.47E+09-2.54E+094.56E+08-2.41E+092.18E+09-2.86E+09-9.97E+07-6.98E+08-2.54E+095.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;%*输入结点坐标

21、数组*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+094.50E+09-3.83E+091.50E+09-4.17E+09-4.50E+09-3.33E+08-1.50E+094.50E+091.33E+10-1.50E+094.67E+09-4.50E+09-6.67E+091.50E+09-1.13E+10-3.83E+09-1.50E+098.33E+09-4.50E+09-3.33E+081.50

22、E+09-4.17E+094.50E+091.50E+094.67E+09-4.50E+091.33E+10-1.50E+09-1.13E+104.50E+09-6.67E+09-4.17E+09-4.50E+09-3.33E+08-1.50E+098.33E+094.50E+09-3.83E+091.50E+09-4.50E+09-6.67E+091.50E+09-1.13E+104.50E+091.33E+10-1.50E+094.67E+09-3.33E+081.50E+09-4.17E+094.50E+09-3.83E+09-1.50E+098.33E+09-4.50E+09-1.50

23、E+09-1.13E+104.50E+09-6.67E+091.50E+094.67E+09-4.50E+091.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=

24、ele_mat_quad4(xy4,mat)运行程序，得到单元刚度矩阵K4（Pa）如下：7.17E+092.50E+09-4.17E+093.50E+09-2.83E+09-3.50E+09-1.67E+08-2.50E+092.50E+098.67E+095.00E+083.33E+09-3.50E+09-1.33E+095.00E+08-1.07E+10-4.17E+095.00E+081.02E+10-6.50E+09-1.67E+085.00E+08-5.83E+095.50E+093.50E+093.33E+09-6.50E+092.07E+10-2.50E+09-1.07E+105

25、.50E+09-1.33E+10-2.83E+09-3.50E+09-1.67E+08-2.50E+097.17E+092.50E+09-4.17E+093.50E+09-3.50E+09-1.33E+095.00E+08-1.07E+102.50E+098.67E+095.00E+083.33E+09-1.67E+085.00E+08-5.83E+095.50E+09-4.17E+095.00E+081.02E+10-6.50E+09-2.50E+09-1.07E+105.50E+09-1.33E+103.50E+093.33E+09-6.50E+092.07E+10实 验 三（40分）1、实验内容编写一个计算平面8结点四边形等参元刚度矩阵的MATLAB函数文件K8 = ele_mat_quad8(xy8,mat)，其中：输入变量xy8为结点坐标数组，mat为材料参数矩阵；输出变量K8为单元刚度矩阵。（要求给出3个不同算例进行验证，并绘制出单元形状和结点号）2、程序代码通用函数function K8 = ele_mat_quad8(xy8,mat)%生成平面8结点四边形等参元刚度矩阵的功能函数%*变量说明*%xy8-结点坐标数组%mat-材料参数矩阵（弹性模量