单元刚度矩阵MATLAB编程.docx

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

2、on 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=*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=z

3、eros(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,; %*输入材料参数矩阵(弹性模量，泊松比，壁厚)*B3,S3,K3=ele_mat_tri3(xy3,mat)三、算例分析 算例1：如图1所示三角形单元，结点坐标为1(0,0)，2(5,2)，3(1,4)，弹性模量为200GPa，泊松比为、厚度为。试求应变矩阵，应力矩阵和单元刚度矩阵。图1 算例1三角形单元解：根据如图1所示三角形单元及其几何和材料参数，编制主程序如下：clear;clc;%*输入结点坐标数组*xy3=0,0;5,2;1,4; mat=2e11,; %*输入材料参数矩阵(弹性模量，泊松比，壁厚)*B3,S

5、3,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，泊松比为、厚度为。

6、试求应变矩阵，应力矩阵和单元刚度矩阵。图2 算例2三角形单元解：根据如图2所示三角形单元及其几何和材料参数，编制主程序如下：clear;clc;%*输入结点坐标数组*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+1

7、0+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,，弹性模量为200GPa，泊松比为、厚度为。试求应变矩阵，应力矩阵和单元刚度矩阵。图3 算例3三角形单元解：根据如图3所示三角形单元及其几何和材料参数，编制主程序如下：clear;clc;%*输入结点坐标数组*xy3=0,0;3,0;,*sqrt(3); mat=2e11,; %*输入材料参数矩阵(弹性模量，泊松比，壁厚)*B3,S3,K3=ele_mat_tri3(xy

8、3,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)，其中：输入变

9、量xy4为结点坐标数组，mat为材料参数矩阵；输出变量K4为单元刚度矩阵。（要求给出3个不同算例进行验证，并绘制出单元形状和结点号）二、程序代码 通用函数function K4 = 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;%*数

10、值积分(Guass, n=4)*C(1)=;C(2)=;C(3)=;C(4)=;A(1)=;A(2)=;A(3)=;A(4)=;sum=0;for i = 1:4 for j = 1:4 y1=C;y2=C; k=1:4;%* PN1(1)=*(y2(j)-1);PN2(1)=*(y1(i)-1); PN1(2)=*(y2(j)-1);PN2(2)=*(y1(i)+1); PN1(3)=*(y2(j)+1);PN2(3)=*(y1(i)+1); PN1(4)=*(y2(j)+1);PN2(4)=*(y1(i)-1); for k=1:4 PN(:,k)=PN1(k),PN2(k); end J

11、=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); su

12、m = sum+A(i)*A(j)*m; endendK4=vpa(sum,9); 主程序clc;clear;%*输入结点坐标数组*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所示四边形等参

13、单元及其几何和材料参数，编制主程序如下：clc;clear;%*输入结点坐标数组*xy4=3,2;8,3;7,8;4,7; %*输入材料参数矩阵(弹性模量，泊松比，壁厚)*mat=2e11,; K4=ele_mat_quad4(xy4,mat)运行程序，得到单元刚度矩阵K4（Pa）如下：+09+09+09+08+09+09+09+08+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

14、+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）如下：+0

15、9+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四边形等参单元解：

16、根据如图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+

17、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结点四边形等参元刚度矩阵的功能函数%*变量说明*%xy8-结点坐标数组%mat-材料参数矩阵（弹性模量，泊松比

18、，壁厚）%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)*(-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

19、)*(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);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*P

20、N(:,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;%

21、*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:4 for j = 1:4 sum = 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; %*输入材料参数矩阵(弹性模量，泊松比，壁厚)*mat=3e6,; K8 = ele_mat_qua

22、d8(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)运行程序，得到单元刚度矩阵K8（GPa）如下： 算例2：如图10所示曲四边形等参单元，已知8个结点整体坐标系内的坐标为1(-5,-4); 2(6,-12); 3(10,8); 4(-8,12); 5(0,-6); 6(5,-1); 7(0