Fortran平面四边形八结点等参元有限元程序.docx
《Fortran平面四边形八结点等参元有限元程序.docx》由会员分享，可在线阅读，更多相关《Fortran平面四边形八结点等参元有限元程序.docx（19页珍藏版）》请在冰豆网上搜索。
Fortran平面四边形八结点等参元有限元程序
%-------------------------------------------------------------------------
%这是一个平面四边形八结点等参元有限元主程序，本程序具有以下功能：
%1,本程序能计算平面应力问题。
%2，自动画出变形图和荷载图。
%3，输出节点位移，节点力,节点应力，和节点主应力。
%时间：
2004年12月
%--------------------------------------------------------------------------
functionN8plane
%变量说明列表:
%NfreeNum节点自由度数（本程序为2）;EleNodeNum单元节点数（本程序为8）
%---------------------------单变量---------------------------------------
%FP1数据文件指针（输入）;NPOIN节点总数;NELEM单元总数;NVFIX约束自由度数;
%E弹性模量;POISSON泊松比;h单元厚度;method问题类型（1.单元应力2.单元应变）
%----------------------------数组----------------------------------------
%LNODS单元定义;COORD坐标;FPOIN点荷载定义;k单刚矩阵;K总刚矩阵
%Oxy单元中点应力;O12单元中点主应力
%----------------------------矢量----------------------------------------
%FIXED约束信息;NFORCE荷载信息;ASLOD等效节点荷载;Ecoord单元坐标
%DISP总体位移;EDisp单元位移
%------------------------------------------------------------------------
%主调函数:
%N8ek生成单刚N8Assemble组装总刚;N8stress求解应力;
%------------------------------------------------------------------------
%
%
%
clear
%-----输入数据存储的文件名------------
file_in=input（'数据输入的文件:
','s'）;
file_out=input（'数据输出的文件:
','s'）;
%--------------------------------------------------------------------------
%打开文件读取数据
%--------------------------------------------------------------------------
FP1=fopen（file_in,'r'）;
NfreeNum=2;EleNodeNum=8;
NPOIN=fscanf（FP1,'%d',1）;NELEM=fscanf（FP1,'%d',1）;NVFIX=fscanf（FP1,'%d',1）;
E=fscanf（FP1,'%f',1）;POISSON=fscanf（FP1,'%f',1）;h=fscanf（FP1,'%f',1）;method=fscanf（FP1,'%d',1）;
%-------读取结构信息--------
LNODS=fscanf（FP1,'%f',[8,NELEM]）';%单元定义
COORD=fscanf（FP1,'%f',[2,NPOIN]）';%坐标
FIXED=fscanf（FP1,'%f',NVFIX）';%约束信息
NFORCE=fscanf（FP1,'%d',2）';%荷载信息
ifNFORCE
（1）~=0
FPOIN=fscanf（FP1,'%f',[3,NFORCE
（1）]）';%节点力
end
%--------------定义非结点力-------------------------------------------------
ifNFORCE
（2）~=0
NF=NFORCE
（1）;
np_number=fscanf（FP1,'%d',1）;%读入均布侧压的单元边数
PBorder=zeros（np_number,6）;
fori=1:
np_number
PBorder（i,1）=fscanf（FP1,'%d',1）;%单元编号
PBorder（i,2）=fscanf（FP1,'%d',1）;%侧压边左结点号
PBorder（i,3）=fscanf（FP1,'%d',1）;%侧压边中结点号
PBorder（i,4）=fscanf（FP1,'%d',1）;%侧压边右结点号
PBorder（i,5）=fscanf（FP1,'%f',1）;%侧压力的大小
PBorder（i,6）=fscanf（FP1,'%d',1）;%侧压力的方向：
1--x方向，2--y方向
%--------------------------------------------------------------------------
%以下语句将均布荷载转化为集中力
%--------------------------------------------------------------------------
ifPBorder（i,6）==1%判断均布力的方向1--x方向，2--y方向
xi=COORD（PBorder（i,2）,1）;%提取侧压边两端点坐标
yi=COORD（PBorder（i,2）,2）;
xj=COORD（PBorder（i,4）,1）;
yj=COORD（PBorder（i,4）,2）;
Blength=sqrt（（xj-xi）^2+（yj-yi）^2）;%边界长度
FPOIN（NF+1,1）=PBorder（i,2）;%给侧压力下端点分配力，并输入到节点力矩阵中
FPOIN（NF+1,2）=PBorder（i,5）*Blength*h/6;
FPOIN（NF+1,3）=0;
FPOIN（NF+2,1）=PBorder（i,3）;%给侧压力中端点分配力，并输入到节点力矩阵中
FPOIN（NF+2,2）=PBorder（i,5）*Blength*h*2/3;
FPOIN（NF+2,3）=0;
FPOIN（NF+3,1）=PBorder（i,4）;%给侧压力上端点分配力，并输入到节点力矩阵中
FPOIN（NF+3,2）=PBorder（i,5）*Blength*h/6;
FPOIN（NF+3,3）=0;
NF=NF+3;
else
xi=COORD（PBorder（i,2）,1）;%提取侧压边两端点坐标，为计算边界长度作准备
yi=COORD（PBorder（i,2）,2）;
xj=COORD（PBorder（i,4）,1）;
yj=COORD（PBorder（i,4）,2）;
Blength=sqrt（（xj-xi）^2+（yj-yi）^2）;%边界长度
FPOIN（NF+1,1）=PBorder（i,2）;%给侧压力左端点分配力，并输入到节点力矩阵中
FPOIN（NF+1,2）=0;
FPOIN（NF+1,3）=PBorder（i,5）*Blength*h/6;
FPOIN（NF+2,1）=PBorder（i,3）;%给侧压力中端点分配力，并输入到节点力矩阵中
FPOIN（NF+2,2）=0;
FPOIN（NF+2,3）=PBorder（i,5）*Blength*h*2/3;
FPOIN（NF+3,1）=PBorder（i,3）;%给侧压力右端点分配力，并输入到节点力矩阵中
FPOIN（NF+3,2）=0;
FPOIN（NF+3,3）=PBorder（i,5）*Blength*h/6;
NF=NF+3;%NF为转化后集中力的个数
end
end
end
%--------------------------------------------------------------------------
%均布荷载转化完毕
%--------------------------------------------------------------------------
fclose（FP1）;
K=zeros（2*NPOIN,2*NPOIN）;
%-----生成单刚,组装总刚-----
fori=1:
NELEM
forj=1:
EleNodeNum
Ecoord（（j-1）*NfreeNum+1:
j*NfreeNum）=COORD（LNODS（i,j）,1:
NfreeNum）;%提取i单元各节点的坐标
end
k=N8ek（E,POISSON,h,Ecoord,method）;%高斯积分求解单刚
K=N8Assemble（K,k,LNODS,NfreeNum,EleNodeNum,i）;%叠加总刚
ifk~=k'
error（'单刚矩阵不对称，请检查'）%判断刚度对称性
end
end
ifK~=K'
error（'总刚矩阵不对称，请检查'）
end
%-------生成荷载向量---------
ASLOD（1:
2*NPOIN）=0;
ASLOD1（1:
2*NPOIN）=0;
fori=1:
NF%NF为转化后集中力的个数
ASLOD1（（FPOIN（i,1）*2-1）:
FPOIN（i,1）*2）=FPOIN（i,2:
3）;
ASLOD=ASLOD+ASLOD1;%实现相同节点力的叠加
ASLOD1（1:
2*NPOIN）=0;
%NPOIN（i,1）为作用点,NPOIN（i,2:
3）为x,y方向的节点力
end
%----------------------------------------------
K1=K;
%-----------加约束-----------
forj=1:
NVFIX
a=FIXED（j）;
K1（1:
NfreeNum*NPOIN,a）=0;K1（a,1:
NfreeNum*NPOIN）=0;K1（a,a）=1;
ASLOD（a）=0;
end
DISP=K1\ASLOD';%求出节点位移
%-----------------
Fnode=K*DISP;%求出节点力
%------求单元节点应力--------
EDisp=zeros（16,1）;
sigma=zeros（NPOIN,3）;
nodetimes=zeros（NPOIN）;%统计节点求解的次数,因不同单元可能共用某个节点
fori=1:
NELEM
forj=1:
EleNodeNum
EDisp（j*2-1:
j*2）=DISP（2*LNODS（i,j）-1:
2*LNODS（i,j））;
end
forj=1:
EleNodeNum
Ecoord（（j-1）*NfreeNum+1:
j*NfreeNum）=COORD（LNODS（i,j）,1:
NfreeNum）;
end
Enodesigma=N8stress（E,POISSON,Ecoord,EDisp,method）;%求解i单元的节点应力
forj=1:
EleNodeNum
num=LNODS（i,j）;
nodetimes（num）=nodetimes（num）+1;
ifnodetimes（num）==0
sigma（num,1:
3）=Enodesigma（j,1:
3）;
else
sigma（num,1:
3）=sigma（num,1:
3）+Enodesigma（j,1:
3）;%公共节点处应力叠加
end
end
end
fori=1:
NPOIN
sigma（i,1:
3）=sigma（i,1:
3）/nodetimes（i）;%i节点应力
mainsigma（i,1:
3）=mainstress2d（sigma（i,1:
3））;%i节点主应力
end
%-------------------------------------------------------------------------
%
%
%--------------------------------------------------------------------------
%程序输出部分
%--------------------------------------------------------------------------
%
FP2=fopen（file_out,'w'）;%打开输出文件
fprintf（FP2,'----------------------------------------------------------------------------------------\n'）;
fprintf（FP2,'这是一个平面四边形八结点等参元有限元主程序，本程序具有以下功能：
\n'）;
fprintf（FP2,'1,本程序能计算平面应力问题。
\n'）;
fprintf（FP2,'2，自动画出变形图和荷载图。
\n'）;
fprintf（FP2,'3，输出节点位移，节点力,节点应力，和主应力。
\n'）;
fprintf（FP2,'编者：
李春林指导老师：
张永山老师时间：
2004年12月\n'）;
fprintf（FP2,'-----------------------------------------------------------------------------------------\n'）
%---------结果数据后处理-----------------------------------------------
fprintf（FP2,'这个算例有%3d个单元共%3d个结点，以下为结果：
\n',NELEM,NPOIN）;
fprintf（FP2,'-----------------------------------------------------------------------------------------\n'）;
%---------输出结点位移------------------------------
fprintf（FP2,'节点位移：
\n'）;
fori=1:
1:
NPOIN
fprintf（FP2,'第%3d号结点,X轴方向位移=%.3e,Y轴方向位移=%.3e\n',...
i,DISP（（i-1）*2+1）,DISP（（i-1）*2+2））;
end
fprintf（FP2,'-----------------------------------------------------------------------------------------\n'）;
%--------结点位移输出完毕----------------------------
%----------输出节点力--------------------------------
fprintf（FP2,'节点力：
\n'）;
fori=1:
1:
NPOIN
fprintf（FP2,'第%3d号结点,X轴方向力=%.3e,Y轴方向力=%.3e\n',...
i,Fnode（（i-1）*2+1）,Fnode（（i-1）*2+2））;
end
fprintf（FP2,'-----------------------------------------------------------------------------------------\n'）;
%--------单元力输出完毕----------------------------
fprintf（FP2,'-----------------------------------------------------------------------------------------\n'）;
%-------输出节点应力--------------------------------
fprintf（FP2,'节点应力：
\n'）;
fori=1:
1:
NPOIN
fprintf（FP2,'第%3d节点,X轴方向=%.3e,Y轴方向=%.3eXY剪应力=%.3d\n',...
i,sigma（i,1）,sigma（i,2）,sigma（i,3））;
end
fprintf（FP2,'-----------------------------------------------------------------------------------------\n'）;
%--------节点应力输出完毕----------------------------
%
fprintf（FP2,'-----------------------------------------------------------------------------------------\n'）;
%-------输出节点主应力--------------------------------
fprintf（FP2,'节点主应力：
\n'）;
fori=1:
1:
NPOIN
fprintf（FP2,'第%3d节点,主应力1=%.3e,主应力2=%.3e方向角（弧度）=%.3d\n',...
i,mainsigma（i,1）,mainsigma（i,2）,mainsigma（i,3））;
end
fprintf（FP2,'----------------------------------------------------------------------------------------\n'）;
%--------单元主应力输出完毕----------------------------
%
%
%-------程序算例图形输出------------------------------
%-------------------------------------------------------------------------
clf
forie=1:
NELEM%绘出计算原图及变形图
plotne（ie,COORD,LNODS）;
ploted（ie,LNODS,DISP,COORD）;
end
%--------------------------------------------------------------------------
%
%
%--------------------------------------------------------------------------
fori=1:
NPOIN%绘制结点编号，即在图中显示各结点编号
nx=COORD（i,1）;
ny=COORD（i,2）;
plot（nx,ny,'.'）;
aa=num2str（i）;
text（nx,ny,aa）;
end
%--------------------------------------------------------------------------
%
%
%--------------------------------------------------------------------------
foribc=1:
NVFIX%绘出边界条件
plotbc（ibc,COORD,FIXED）;
end
%-------------------------------------------------------------------------
forik=1:
NFORCE
（1）%绘出集中荷载图
plotfxy（ik,COORD,FPOIN）;
end
ifNFORCE
（2）~=0
forip=1:
np_number%绘出均布荷载图
plotmp（ip,PBorder,COORD,h）;
end
end
%------------------------------------------------------------------------
fclose（FP2）;%关闭输出文件
return;
functionD=Dmtrx（E,NU,method）
%各向同性单元弹性矩阵
ifNU<0||NU>=1
error（'thePoissonratioshouldonlybegreaterthan0andlessthan1.'）
end
ifmethod==1
%平面应力
D=E/（1-NU^2）*[1NU0;NU10;00（1-NU）/2];
elseifmethod==2
%平面应变
D=E*（1-NU）/（（1+NU）*（1-2*NU））*[1NU/（1-NU）0;NU/（1-NU）10;00（1-2*NU）/2/（1-NU）];
else
error（'variable"method"shouldonlybe1:
planestressor2:
planestrain.'）
end
functionmainsigma=mainstress2d（sigma）
%平面应力下的主应力及方向（弧度）
mainsigma
（1）=1/2*（sigma
（1）+sigma
（2））+1/2*sqrt（（sigma
（1）-sigma
（2））^2+4*sigma（3）^2）;
mainsigma
（2）=1/2*（sigma
（1）+sigma
（2））-1/2*sqrt（（sigma
（1）-sigma
（2））^2+4*sigma（3）^2）;
ifabs（sigma
（1）-sigma
（2））<1e-30
ifsigma（3）>1e-5
mainsigma（3）=pi/4;
elseifsigma（3）<-1e-5
mainsigma（3）=-pi/4;
else
mainsigma（3）=0;
end
else
mainsigma（3）=0.5*atan（2.0*sigma（
升级会员 