单元刚度矩阵(等参元)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:3ifi==1j=2;m=3;elseifi==2j=3;m=1;elsej=1;m=2;endb(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.11110.00000.22220.0000-0.11110.00000.0000-0.22220.0000-0.05560.00000.2778-0.2222-0.1111-0.05560.22220.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.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.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];%****输入材料参数矩阵(弹性模量，泊松比，壁厚)****[B3,S3,K3]=ele_mat_tri3(xy3,mat)运行程序，得到应变矩阵B3如下：-0.33330.00000.33330.00000.00000.00000.0000-0.20000.00000.00000.00000.2000-0.2000-0.33330.00000.33330.20000.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-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+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√3)，弹性模量为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.33330.00000.33330.00000.00000.00000.0000-0.19250.0000-0.19250.00000.3849-0.1925-0.3333-0.19250.33330.38490.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.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+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个不同算例进行验证，并绘制出单元形状和结点号）二、程序代码通用函数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:4forj=1:4y1=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:4PN(:,k)=[PN1(k),PN2(k)]';endJ=PN*xy4;JN=inv(J