《单元刚度矩阵(等参元)MATLAB编程(共16页).docx》由会员分享，可在线阅读，更多相关《单元刚度矩阵(等参元)MATLAB编程(共16页).docx（16页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上有限元法实验报告 专业班级 力学（实验）1601 姓 名 田诗豪 学 号 提交日期 2019.4.24 实验编号实验一实验二实验三总分得 分实 验 一（30分）一、实验内容编写一个计算平面3结点三角形单元的应变矩阵、应力矩阵和单元刚度矩阵的MATLAB函数文件B3,S3,K3 = ele_mat_tri3(xy3,mat)，其中：输入变量xy3为结点坐标数组，mat为材料参数矩阵；输出变量B3为应变矩阵，S3为应力矩阵，K3为单元刚度矩阵。（要求给出3个不同算例进行验证，并绘制出单元形状和结点号）二、程序代码l 通用函数function B3,S3,K3 = ele_
2、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);%*for
3、 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;l 主程序clear;clc;%*输入结点坐标数组*xy3
4、=0,0;5,1;1,4; mat=3e6,0.5,1.0; %*输入材料参数矩阵(弹性模量，泊松比，壁厚)*B3,S3,K3=ele_mat_tri3(xy3,mat)三、算例分析l 算例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.86E+096.3
6、3E+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.42E+09-3
7、.05E+102.28E+103.84E+10-2.14E+10-2.85E+09-5.92E+102.42E+10-2.41E+10-2.14E+108.33E+10l 算例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; %*输入材料参数矩阵(弹性模量，泊松比，壁厚)*B
8、3,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.48E+10-
9、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+101.11E
10、+100.00E+00-1.99E+10-3.42E+101.99E+100.00E+000.00E+003.42E+10l 算例3：如图3所示三角形单元，结点坐标为1(0,0)，2(3,0)，3(1.5,1.53)，弹性模量为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,K
11、3=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+10-2.
12、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.14E
13、+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个不同算例进行验证，并绘制出单元形状和结点号）二、程序代码l 通用函数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.4055;C(2)=0.4886;C(3)=-0.4886;C(4)=-0.4055;A(1)=0.7454;A(2)=0.2546;A(3)=0.2546;A(4)=0.7454;sum=0;for i = 1:4 for j = 1:4 y1=C;y2=C;
15、 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)=PN1(k),PN2(k); end J=PN*xy4; JN=inv(J); J1=JN(1,:);J2=JN(2,:);%*应变矩阵* B1=J1*PN(:,1),0;0,
16、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; endendK4=vpa(sum,9);l 主程序clc;clear;%*输入结点坐标数组*x
17、y4=3,2;8,3;7,8;4,7; %*输入材料参数矩阵(弹性模量，泊松比，壁厚)*mat=3e6,0.5,1.0; K4=ele_mat_quad4(xy4,mat)三、算例分析l 算例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
18、,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.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.
19、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.08E+095.98E+08-4.47E+092.18E+09-5.08E+09-2.42E+088.47E+09-2.54E+094.56E+08-2.41E+092.18
20、E+09-2.86E+09-9.97E+07-6.98E+08-2.54E+095.97E+09l 算例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
21、(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.50E+09-4.17E+094.50E+091.50E+094.67E+09-4.50E+091.33E+10-1.50E+09-1.13E+104.50E+09-6.6
22、7E+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.50E+09-1.13E+104.50E+09-6.67E+091.50E+094.67E+09-4.50E+091.33E+10l 算例3：如图7所示四边形等参单元，已知
23、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+092.50E+09-4.17E+093.50E+09-2.83
24、E+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.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
25、+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为结点坐标数组，
26、mat为材料参数矩阵；输出变量K8为单元刚度矩阵。（要求给出3个不同算例进行验证，并绘制出单元形状和结点号）2、程序代码l 通用函数function K8 = ele_mat_quad8(xy8,mat)%生成平面8结点四边形等参元刚度矩阵的功能函数%*变量说明*%xy8-结点坐标数组%mat-材料参数矩阵（弹性模量，泊松比，壁厚）%K8-单元刚度矩阵%D-弹性矩阵%*PN = sym(zeros(2,8);D=mat(1)/(1-mat(2)2)*1,mat(2),0;mat(2),1,0;0,0,(1-mat(2)/2;syms y1 y2 realN(1)=(1-y1)*(1-y2)*(
27、-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-y12)*(1-y2)/2;N(6)=(1-y22)*(1+y1)/2;N(7)=(1-y12)*(1+y2)/2;N(8)=(1-y1)*(1-y22)/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=inv(J);J
28、1=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
29、,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)=0.4055; C(2)=0.4886;C(3)=-0.4886; C(4)=-0.4055;A(1)=0.7454; A(2)=0.2546;A(3)=0.2546; A(4)=
30、0.7454;sum=0;for i = 1:4 for j = 1:4 sum = sum+A(i)*A(j)*subs(m,y1,y2,C(i),C(j); endendK8=vpa(sum,9);l 主程序clc;clear;%*输入结点坐标数组*xy8=2,2;7,2;6,7;3,6;5,3;6,4;4,6;3,4; %*输入材料参数矩阵(弹性模量，泊松比，壁厚)*mat=3e6,0.5,1.0; K8 = ele_mat_quad8(xy8,mat)3、算例分析l 算例1：如图9所示曲四边形等参单元（图8为局部坐标系规则单元），已知8个结点整体坐标系内的坐标为1(2,2); 2(7,
31、2); 3(6,7); 4(3,6); 5(5,3); 6(6,4); 7(4,6); 8(3,4)。弹性模量为100GPa，泊松比为0.35、厚度为0.05m。试求单元刚度矩阵。图8 局部坐标系规则单元图9 算例1曲四边形等参单元解：根据如图9所示曲四边形等参单元及其几何和材料参数，编制主程序如下：clc;clear;%*输入结点坐标数组*xy8=2,2;7,2;6,7;3,6;5,3;6,4;4,6;3,4; %*输入材料参数矩阵(弹性模量，泊松比，壁厚)*mat=1.8e11,0.5,0.1; K8 = ele_mat_quad8(xy8,mat)运行程序，得到单元刚度矩阵K8（GPa）
32、如下：3.79 1.47 1.00 -0.10 0.91 0.16 1.92 -0.12 -7.51 -3.49 4.99 2.55 -2.87 -0.50 -2.22 0.03 1.47 5.37 -0.08 0.68 0.16 1.42 -0.14 2.40 -3.59 -3.08 2.55 2.15 -0.50 -1.92 0.12 -7.03 1.00 -0.08 2.05 -0.55 1.87 0.32 2.00 -0.45 -1.78 1.38 -2.52 -0.49 -2.63 -0.13 0.02 0.00 -0.10 0.68 -0.55 2.07 0.34 2.40 -0.
33、45 1.50 1.48 2.24 -0.58 -6.13 -0.13 -1.23 0.00 -1.55 0.91 0.16 1.87 0.34 5.87 2.38 3.38 0.19 2.76 2.16 -6.15 -3.32 -7.41 -1.51 -1.21 -0.40 0.16 1.42 0.32 2.40 2.38 5.90 0.21 2.42 2.16 4.79 -3.22 -10.36 -1.61 -3.15 -0.40 -3.43 1.92 -0.14 2.00 -0.45 3.38 0.21 6.68 -1.15 -2.68 0.08 -1.51 0.36 -9.18 -0.
34、21 -0.61 1.30 -0.12 2.40 -0.45 1.50 0.19 2.42 -1.15 3.98 0.08 0.39 0.36 -3.02 -0.12 -3.51 1.21 -4.15 -7.51 -3.59 -1.78 1.48 2.76 2.16 -2.68 0.08 36.97 11.70 -29.09 -15.32 4.65 0.41 -3.32 3.08 -3.49 -3.08 1.38 2.24 2.16 4.79 0.08 0.39 11.70 31.58 -15.32 -29.16 0.41 -2.12 3.08 -4.65 4.99 2.55 -2.52 -0
35、.58 -6.15 -3.22 -1.51 0.36 -29.09 -15.32 37.14 14.78 -0.72 2.09 -2.14 -0.66 2.55 2.15 -0.49 -6.13 -3.32 -10.36 0.36 -3.02 -15.32 -29.16 14.78 40.11 2.09 1.91 -0.66 4.50 -2.87 -0.50 -2.63 -0.13 -7.41 -1.61 -9.18 -0.12 4.65 0.41 -0.72 2.09 18.97 1.45 -0.81 -1.60 -0.50 -1.92 -0.13 -1.23 -1.51 -3.15 -0.
36、21 -3.51 0.41 -2.12 2.09 1.91 1.45 9.72 -1.60 0.30 -2.22 0.12 0.02 0.00 -1.21 -0.40 -0.61 1.21 -3.32 3.08 -2.14 -0.66 -0.81 -1.60 10.29 -1.74 0.03 -7.03 0.00 -1.55 -0.40 -3.43 1.30 -4.15 3.08 -4.65 -0.66 4.50 -1.60 0.30 -1.74 16.00 l 算例2：如图10所示曲四边形等参单元，已知8个结点整体坐标系内的坐标为1(-5,-4); 2(6,-12); 3(10,8); 4(
37、-8,12); 5(0,-6); 6(5,-1); 7(0,4); 8(-4,2)。弹性模量为200GPa，泊松比为0.35、厚度为0.1m。试求单元刚度矩阵。图10 算例2曲四边形等参单元解：根据如图10所示曲四边形等参单元及其几何和材料参数，编制主程序如下：clc;clear;%*输入结点坐标数组*xy8=-5,-4;6,-12;10,8;-8,12;0,-6;5,-1;0,4;-4,2; %*输入材料参数矩阵(弹性模量，泊松比，壁厚)*mat=2e11,0.35,0.1; K8 = ele_mat_quad8(xy8,mat)运行程序，得到单元刚度矩阵K8（GPa）如下：40.8 11.0 16.5 1.2 10.6 2.9 11.6 2.9 -33.1 -0.6 -12.2 -4.7 -3.7 3.2 -30.6 -15.9 11.0 22.7 1.3 9.6 2.9 7.2 2.8 9.7 -0.9 -9.3 -4.7 -13.5 3.2 3.6 -15.5 -29.8 16.5 1.3 16.4 -3.9 9.4 1.1 5.5 -0.6 -26