四边形八节点等参元matlab程序.docx
《四边形八节点等参元matlab程序.docx》由会员分享，可在线阅读，更多相关《四边形八节点等参元matlab程序.docx（14页珍藏版）》请在冰豆网上搜索。
四边形八节点等参元matlab程序
悬臂钢梁，尺寸如图一所示；v=0.3。
h=1，E=2.1e11.
图一悬臂钢梁
图二单元划分与结点编号
Matlab输出结果
附录Ⅰ：
有限元ANSYS分析结果
采用PLANE183单元（四边形八节点）单元得出的结构Y向最大位移为-0.216E-04。
约等于matlab平面四边形八节点等参元结点Y向最大位移-2.4024E-5。
附录Ⅱ：
%---------------四边形八节点等参元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）;%张成总刚矩阵并清零
%调用子程序生成单元刚度矩阵
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）;%张成总荷载向量并清零
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]=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];
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）;
bjd.txt文件数据
2.1E110.31528813
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
3.01.0
3.51.0
4.01.0
4.51.0
5.01.0
170-10000
111
1211
1811
升级会员 