单元刚度矩阵(等参元)MATLAB编程
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
《有限元法》实验报告
专业班级力学（实验）1601
姓名田诗豪
学号 10
提交日期
实验一（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=*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 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;
%*********输入结点坐标数组********
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,S3,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，泊松比为、厚度为。试求应变矩阵，应力矩阵和单元刚度矩阵。
图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+10+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,,];