四边形八节点等参元matlab程序.docx
《四边形八节点等参元matlab程序.docx》由会员分享，可在线阅读，更多相关《四边形八节点等参元matlab程序.docx（12页珍藏版）》请在冰豆网上搜索。
四边形八节点等参元matlab程序
悬臂钢梁，尺寸如图一所示；v=0.3。
h=1，E=2.1e11.
图一悬臂钢梁
图二单元划分与结点编号
Matlab输出结果
附录Ⅰ：
有限元ANSYS分析结果
采用PLANE183单元〔四边形八节点〕单元得出的结构Y向最大位移为6E-04。
约等于matlab平面四边形八节点等参元结点Y向最大位移。
附录Ⅱ：
%---------------四边形八节点等参元matlab计算程序----------------------------
%———————————主程序—————————
%*******************************************************************%************************************
%2012年
%本程序只能处理集中荷载作用下的情况
%只输出了节点位移、单元中心点的应力
%*******************************************************************%***************
%变量说明
%Evh
%弹性模量泊松比厚度
%NPOINNELEMNVFIXNNODENFPOIN
%总结点数,单元数,约束结点个数,单元节点数,受力结点数
%COORDLNODS
%结构节点整体坐标数组,单元定义数组,
%FPOINFORCEFIXED
%结点力数组，总体荷载向量,约束信息数组
%HKDISP
%总体刚度矩阵,结点位移向量
%******************************
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]）';%结点号x,y坐标（整体坐标下）
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
%*******************************************************************
%*******************************************************************
%========平面应力问题的求解==============
%
%*******************************************************************
%*******************************************************************
%—————————————————————
%刚度矩阵的生成
%计算刚度矩阵，并对约束条件进展处理
Ke=zeros（2*NNODE,2*NNODE）;%单元刚度矩阵并清零
HK=zeros（2*NPOIN,2*NPOIN）;%X成总刚矩阵并清零
%调用子程序生成单元刚度矩阵
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
end
end
%—————————————————————————————————
%对荷载向量进展处理
FORCE=zeros（2*NPOIN,1）;%X成总荷载向量并清零
fori=1:
NFPOIN
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
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
c2=2*FIXED（i,1）;
HK（c2,:
）=0;
HK（:
c2）=0;
HK（c2,c2）=1;
FORCE（c2）=0;
end
end
%—————————————————————————————————
%===========================================================
%===========================================================
DISP=HK\FORCE%计算节点位移向量
%===========================================================
%===========================================================
%———————————求解单元应力————————————————
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
stress%输出应力
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个角结点的坐标
%矩阵尺寸为16x16
Ke=zeros（16,16）;
D=（E/（1-v*v））*[1v0;v10;00（1-v）/2];%弹性矩阵
%高斯积分采用3x3个积分点书74页
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
functionB=eleB（x1,y1,x3,y3,x5,y5,x7,y7,s,t）
%调用导函数
[N_s,N_t]=DHS（s,t）;
%求Jacobi矩阵
J=Jacobi（x1,y1,x3,y3,x5,y5,x7,y7,s,t）;
%求应变矩阵B
B=zeros（3,16）;
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）;
functionJ=Jacobi（x1,y1,x3,y3,x5,y5,x7,y7,s,t）
%-------Jacobi-----------
%单元坐标
%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]=DH
升级会员 