四边形八节点等参元matlab程序Word文档下载推荐.doc
《四边形八节点等参元matlab程序Word文档下载推荐.doc》由会员分享，可在线阅读，更多相关《四边形八节点等参元matlab程序Word文档下载推荐.doc（10页珍藏版）》请在冰豆网上搜索。
%只输出了节点位移、单元中心点的应力
%*******************************************************************%***************
%变量说明
%Evh
%弹性模量泊松比厚度
%NPOINNELEMNVFIXNNODENFPOIN
%总结点数,单元数,约束结点个数,单元节点数,受力结点数
%COORDLNODS
%结构节点整体坐标数组,单元定义数组,
%FPOINFORCEFIXED
%结点力数组，总体荷载向量,约束信息数组
%HKDISP
%总体刚度矩阵,结点位移向量
%******************************
clearall
formatshorte
FP1=fopen（'
bjd.txt'
'
rt'
）;
%打开数据文件
%%读入控制数据
E=fscanf（FP1,'
%f'
1）;
%弹性模量
v=fscanf（FP1,'
%泊松比
h=fscanf（FP1,'
%厚度
NELEM=fscanf（FP1,'
%d'
%单元数
NPOIN=fscanf（FP1,'
%总结点数
NNODE=fscanf（FP1,'
%单元节点数
NFPOIN=fscanf（FP1,'
%受力结点数
NVFIX=fscanf（FP1,'
%约束结点个数
LNODS=fscanf（FP1,'
[NNODE,NELEM]）'
;
%单元定义：
单元结点号（逆时针）
COORD=fscanf（FP1,'
[2,NPOIN]）'
%结点号x,y坐标（整体坐标下）
FPOIN=fscanf（FP1,'
[3,NFPOIN]）'
%节点力：
结点号、X方向力（向右正），Y方向力（向上正）
FIXED=fscanf（FP1,'
[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:
HK（（a（j）*2-1）:
a（j）*2,（a（k）*2-1）:
a（k）*2）=HK（（a（j）*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方向的节点力
%将约束信息加入总刚,总荷载
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,:
c2）=0;
HK（c2,c2）=1;
FORCE（c2）=0;
end
%===========================================================
DISP=HK\FORCE%计算节点位移向量
%———————————求解单元应力————————————————
stress=zeros（3,NELEM）;
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）;
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）;
stressm=D*BM*u'
stress（:
m）=stressm;
stress%输出应力
functionKe=K（E,v,h,x1,y1,x3,y3,x5,y5,x7,y7）
%=========单元刚度矩阵===============
%E弹性模量
%v泊松比
%h厚度
%x1,y1,x3,
升级会员 