四边形八节点等参元matlab程序Word格式文档下载.docx
《四边形八节点等参元matlab程序Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《四边形八节点等参元matlab程序Word格式文档下载.docx（12页珍藏版）》请在冰豆网上搜索。
%弹性模量泊松比厚度
%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）;
%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:
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
%—————————————————————————————————
%对荷载向量进展处理
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方向的节点力
%将约束信息参加总刚,总荷载
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;
ifFIXED（i,3）==1
c2=2*FIXED（i,1）;
HK（c2,:
c2）=0;
HK（c2,c2）=1;
FORCE（c2）=0;
%===========================================================
DISP=HK\FORCE%计算节点位移向量
%———————————求解单元应力————————————————
stress=zeros（3,NELEM）;
NELEM
u（1:
16）=0;
d=LNODS（m,:
NNODE
u（i*2-1:
i*2）=DISP（d（i）*2-1:
d（i）*2）;
%从总位移向量中取出当前单元的节点位移
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,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,y3,x5,y5,x7,y7为4个角结点的坐标
%矩阵尺寸为16x16
Ke=zeros（16,16）;
D=（E/（1-v*v））*[1v0;
%高斯积分采用3x3个积分点书74页
W1=5/9;
W2=8/9;
W3=5/9;
%加权系数
W=[W1W2W3];
r=15^（1/2）/5;
x=[-r0r];
%积分点
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;
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）;
B（1:
B2iB1i];
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
升级会员 