有限元四边形八节点程序解析.docx
- 文档编号:27731991
- 上传时间:2023-07-04
- 格式:DOCX
- 页数:15
- 大小:23.08KB
有限元四边形八节点程序解析.docx
《有限元四边形八节点程序解析.docx》由会员分享,可在线阅读,更多相关《有限元四边形八节点程序解析.docx(15页珍藏版)》请在冰豆网上搜索。
有限元四边形八节点程序解析
_______土木工程______学院___15_____级__结构工程____专业姓名__杨尚荣______学号__2111516045____
………………………………(密)………………………………(封)………………………………(线)………………………………
广州大学
作业(论文)题目:
平面四边形八节点有限元程序
所修课程名称:
有限元法及其应用
任课教师姓名:
张永山
上课时间:
星期五5-8节
布置作业(论文)日期:
2016年6月28日
任课教师打分:
2016年6月28日
悬臂钢梁,尺寸如图所示;v=0.3。
h=1,E=2.1e11.
悬臂钢梁
单元划分与结点编号
原始输入数据(bjd.txt)
2.1E110.31528824
1231320191812
3451422212013
5671524232214
7891626252415
910111728272616
0.00.0
0.50.0
1.00.0
1.50.0
2.00.0
2.50.0
3.00.0
3.50.0
4.00.0
4.50.0
5.00.0
0.00.5
1.00.5
2.00.5
3.00.5
4.00.5
5.00.5
0.01.0
0.51.0
1.01.0
1.51.0
2.01.0
2.51.0
%******************************************************************
%四边形八节点等参元matlab计算程序
%主程序
%******************************************************************
%变量说明
%Evh
%弹性模量泊松比厚度
%NPOINNELEMNVFIXNNODENFPOIN
%总结点数,单元数,约束结点个数,单元节点数,受力结点数
%COORDLNODS
%结构节点整体坐标数组,单元定义数组,
%FPOINFORCEFIXED
%结点力数组,总体荷载向量,约束信息数组
%HKDISPFP1
%总体刚度矩阵结点位移向量数据指针文件
%BDstress
%单元应变矩阵单元弹性矩阵单元应力
%******************************************************************
clc%清屏
clearall%清除内存变量
formatshorte%设定输出类型
FP1=fopen('bjd.txt','rt');%打开数据文件
%%读入控制数据
E=fscanf(FP1,'%f',1);%弹性模量
v=fscanf(FP1,'%f',1);%泊松比
h=fscanf(FP1,'%f',1);%厚度
NELEM=fscanf(FP1,'%d',1);%单元数
NPOIN=fscanf(FP1,'%d',1);%总结点数
NNODE=fscanf(FP1,'%d',1);%单元节点数
NFPOIN=fscanf(FP1,'%d',1);%受力结点数
NVFIX=fscanf(FP1,'%d',1);%约束结点个数
LNODS=fscanf(FP1,'%f',[NNODE,NELEM])';%单元定义数组:
单元结点号(逆时针)
COORD=fscanf(FP1,'%f',[2,NPOIN])';%结构节点坐标数组(整体坐标下)
FPOIN=fscanf(FP1,'%f',[3,NFPOIN])';
%节点力数组:
结点号、X方向力(向右正),Y方向力(向上正)
FIXED=fscanf(FP1,'%d',[3,NVFIX])';
%约束信息数组(n,3)n:
受约束节点数目,(n,1):
约束点号
%(n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0
[HK,FORCE]=H(E,v,h,NELEM,NPOIN,NFPOIN,NVFIX,LNODS,COORD,FPOIN,FIXED);
%形成总纲和总荷载
DISP=HK\FORCE%求位移
stress=STRESS(NELEM,COORD,LNODS,NNODE,DISP,E,v)%求应力
%********************************************************************
%生成总刚矩阵和总荷载
%********************************************************************
function[HK,FORCE]=H(E,v,h,NELEM,NPOIN,NFPION,NVFIX,LNODS,COORD,FPOIN,FIXED)
HK=zeros(2*NPOIN,2*NPOIN);
form=1:
NELEM%m为单元号
Ke=K(E,v,h,...
COORD(LNODS(m,1),1),COORD(LNODS(m,1),2),...
COORD(LNODS(m,3),1),COORD(LNODS(m,3),2),...
COORD(LNODS(m,5),1),COORD(LNODS(m,5),2),...
COORD(LNODS(m,7),1),COORD(LNODS(m,7),2));%调用单元刚度矩阵(读取四个角节点)
a=LNODS(m,:
);%临时向量,用来记录当前单元的节点编号
%对总刚度矩阵的处理
forj=1:
8%对行进行循环(按照节点号循环)
fork=1:
8%对列进行循环(按照节点号循环)
HK((a(j)*2-1):
a(j)*2,(a(k)*2-1):
a(k)*2)=HK((a(j)*2-1):
a(j)*2,(a(k)*2-1):
a(k)*2)+...
Ke(j*2-1:
j*2,k*2-1:
k*2);
end%把单刚2乘2子块加入到总纲矩阵中(其中a(j)*2,a(k)*2为求取j节点和k节点第二个位移的整体编码)
end
end
%—————————————————————————————————
%对荷载向量进行处理
FORCE=zeros(2*NPOIN,1);%张成总荷载向量并清零
fori=1:
NFPION%对节点荷载个数进行循环
b1=FPOIN(i,1)*2-1;
b2=FPOIN(i,1)*2;%FPION(i,1)为作用点
%荷载对应的节点号
FORCE(b1)=FPOIN(i,2);%FPION(i,2)为x方向的节点力
FORCE(b2)=FPOIN(i,3);%FPION(i,3)为y方向的节点力
end
%—————————————————————————————————
%将约束信息加入总刚---划零置一(即对总纲和总荷载进行边界条件的处理)
fori=1:
NVFIX%对约束个数进行循环
ifFIXED(i,2)==1%如果x方向存在约束
c1=2*FIXED(i,1)-1;%零位移所在位置
HK(c1,:
)=0;%将一约束序号处的总刚列向量清0
HK(:
c1)=0;%将一约束序号处的总刚行向量清0
HK(c1,c1)=1;%将行列交叉处的元素置为1
FORCE(c1)=0;
end
ifFIXED(i,3)==1%如果y方向存在约束
c2=2*FIXED(i,1);
HK(c2,:
)=0;
HK(:
c2)=0;
HK(c2,c2)=1;
FORCE(c2)=0;
end
End
%********************************************************************
%求单元应力
%********************************************************************
functionstress=STRESS(NELEM,COORD,LNODS,NNODE,DISP,E,v)
stress=zeros(3,NELEM);
form=1:
NELEM%对单元个数进行循环
u(1:
16)=0;%单元位移列向量清零;
d=LNODS(m,:
);%临时向量,用来记录当前单元的节点编号
fori=1:
NNODE%对节点个数进行循环
u(i*2-1:
i*2)=DISP(d(i)*2-1:
d(i)*2);
%从总位移向量中取出当前单元的节点位移
end
D=(E/(1-v*v))*[1v0;v10;00(1-v)/2];%弹性矩阵
%形成应变矩阵BM
BM=zeros(3,16);
fori=1:
NNODE
J=Jacobi(COORD(LNODS(m,1),1),COORD(LNODS(m,1),2),...
COORD(LNODS(m,3),1),COORD(LNODS(m,3),2),...
COORD(LNODS(m,5),1),COORD(LNODS(m,5),2),...
COORD(LNODS(m,7),1),COORD(LNODS(m,7),2),0,0);
[N_s,N_t]=DHS(0,0);
B1i=J(2,2)*N_s(i)-J(1,2)*N_t(i);
B2i=-J(2,1)*N_s(i)+J(1,1)*N_t(i);
BM(1:
3,2*i-1:
2*i)=[B1i0;0B2i;B2iB1i]/det(J);
end
stressm=D*BM*u';
stress(:
m)=stressm;
end%输出应力
%********************************************************************
%求单元刚度矩阵
%********************************************************************
functionKe=K(E,v,h,x1,y1,x3,y3,x5,y5,x7,y7)
%E弹性模量v泊松比h厚度
%x1,y1,x3,y3,x5,y5,x7,y7为4个角结点的坐标
%矩阵尺寸为16乘16
Ke=zeros(16,16);%即(2*NNODE,2*NNODE)
D=(E/(1-v*v))*[1v0;v10;00(1-v)/2];%弹性矩阵
%高斯积分采用3x3个积分点
W1=5/9;W2=8/9;W3=5/9;%加权系数
W=[W1W2W3];
r=15^(1/2)/5;
x=[-r0r];%积分点
fori=1:
3
forj=1:
3
B=eleB(x1,y1,x3,y3,x5,y5,x7,y7,x(i),x(j));
J=Jacobi(x1,y1,x3,y3,x5,y5,x7,y7,x(i),x(j));
Ke=Ke+W(i)*W(j)*B'*D*B*det(J)*h;
end%求解单元刚度矩阵(根据虚功原理)
End
%********************************************************************
%求雅克比矩阵
%********************************************************************
functionJ=Jacobi(x1,y1,x3,y3,x5,y5,x7,y7,s,t)
%单元坐标
%用1,3,5,7点的坐标表示2,4,6,8点的坐标
x2=(x1+x3)/2;y2=(y1+y3)/2;
x4=(x3+x5)/2;y4=(y3+y5)/2;
x6=(x5+x7)/2;y6=(y5+y7)/2;
x8=(x7+x1)/2;y8=(y7+y1)/2;
x=[x1x2x3x4x5x6x7x8];
y=[y1y2y3y4y5y6y7y8];
%调用形函数对局部坐标的导数
[N_s,N_t]=DHS(s,t);
%求Jacobi矩阵的行列式的值
x_s=0;y_s=0;
x_t=0;y_t=0;
fori=1:
8
x_s=x_s+N_s(i)*x(i);y_s=y_s+N_s(i)*y(i);
x_t=x_t+N_t(i)*x(i);y_t=y_t+N_t(i)*y(i);
end
J=[x_sy_s;x_ty_t];
%********************************************************************
%求应变矩阵B
%********************************************************************
functionB=eleB(x1,y1,x3,y3,x5,y5,x7,y7,s,t)
[N_s,N_t]=DHS(s,t);%调用导函数
J=Jacobi(x1,y1,x3,y3,x5,y5,x7,y7,s,t);%求Jacobi矩阵
B=zeros(3,16);%求应变矩阵B
fori=1:
8
B1i=J(2,2)*N_s(i)-J(1,2)*N_t(i);
B2i=-J(2,1)*N_s(i)+J(1,1)*N_t(i);
B(1:
3,2*i-1:
2*i)=[B1i0;0B2i;B2iB1i];
end
B=B/det(J);
%********************************************************************%设置形函数
%********************************************************************
functionN=shape(s,t)
%ξ,η
N
(1)=(1-s)*(1-t)*(-s-t-1)/4;
N(3)=(1+s)*(1-t)*(s-t-1)/4;
N(5)=(1+s)*(1+t)*(s+t-1)/4;
N(7)=(1-s)*(1+t)*(-s+t-1)/4;
N
(2)=(1-t)*(1+s)*(1-s)/2;
N(4)=(1+s)*(1+t)*(1-t)/2;
N(6)=(1+t)*(1+s)*(1-s)/2;
N(8)=(1-s)*(1+t)*(1-t)/2;
%********************************************************************
%型函数的导函数
%********************************************************************
function[N_s,N_t]=DHS(s,t)
%形函数求导
%ξ,η
N_s
(1)=-1/4*(1-t)*(-s-t-1)-1/4*(1-s)*(1-t);
N_s(3)=1/4*(1-t)*(s-t-1)+1/4*(1+s)*(1-t);
N_s(5)=1/4*(1+t)*(s+t-1)+1/4*(1+s)*(1+t);
N_s(7)=-1/4*(1+t)*(-s+t-1)-1/4*(1-s)*(1+t);
N_s
(2)=1/2*(1-s)*(1-t)-1/2*(1+s)*(1-t);
N_s(4)=1/2*(1+t)*(1-t);
N_s(6)=1/2*(1-s)*(1+t)-1/2*(1+s)*(1+t);
N_s(8)=-1/2*(1+t)*(1-t);
N_t
(1)=-1/4*(1-s)*(-s-t-1)-1/4*(1-s)*(1-t);
N_t(3)=-1/4*(1+s)*(s-t-1)-1/4*(1+s)*(1-t);
N_t(5)=1/4*(1+s)*(s+t-1)+1/4*(1+s)*(1+t);
N_t(7)=1/4*(1-s)*(-s+t-1)+1/4*(1-s)*(1+t);
N_t
(2)=-1/2*(1+s)*(1-s);
N_t(4)=1/2*(1+s)*(1-t)-1/2*(1+s)*(1+t);
N_t(6)=1/2*(1+s)*(1-s);
N_t(8)=1/2*(1-s)*(1-t)-1/2*(1-s)*(1+t);
结果输出
DISP=
0
0
-8.5685e-10
-2.2935e-08
1.7374e-08
-6.6907e-08
4.7180e-08
-1.2425e-07
9.5082e-08
-1.1254e-07
9.6672e-08
-7.9875e-08
0
0
-2.3101e-07
-4.5306e-07
-4.1913e-07
-1.0324e-06
-5.0975e-07
-1.7562e-06
-5.2554e-07
-2.4610e-06
0
0
4.8027e-08
-8.1458e-08
1.0689e-07
-1.6005e-07
1.6537e-07
-8.0480e-08
1.8846e-07
-1.0125e-06
1.7417e-07
-2.5257e-06
0
0
5.1705e-08
-3.6646e-08
8.7374e-08
-8.7254e-08
1.0378e-07
-1.4436e-07
1.1930e-07
-2.4763e-07
1.9278e-07
-9.9302e-08
3.6992e-07
-1.7104e-07
6.0259e-07
-4.6703e-07
7.6873e-07
-1.0457e-06
8.9733e-07
-1.7467e-06
9.7325e-07
-2.6221e-06
stress=
1.0134e+041.2190e+041.2150e+044.3631e+03-2.6402e+03
1.6096e+02-5.6745e+02-4.3455e+02-1.6249e+031.2050e+03
-2.3339e+03-1.7762e+031.4189e+04-7.9471e+03-8.5723e+03
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 有限元 四边形 八节 程序 解析