8节点等参单元 有限元程序.docx
《8节点等参单元 有限元程序.docx》由会员分享，可在线阅读，更多相关《8节点等参单元 有限元程序.docx（17页珍藏版）》请在冰豆网上搜索。
8节点等参单元有限元程序
programstress
!
character*64fname1,fname2
dimensionake（16,16）
dimensionld（1000）,is（16）,a（30000）
common/sol1/nn,ne,nm,nc
common/sol2/nf,nd,n,nfd
common/shu1/mea（200,9）
common/shu2/x（1000）,y（1000）
common/shu3/aeu（10,3）
common/shu4/p（1000）
common/shu5/jz（1000,2）
!
write（*,'（a22）'）'inputdatafile:
'
!
read（*,'（a）'）fname1
!
write（*,'（a22）'）'outputdatafile:
'
!
read（*,'（a）'）fname2
open（1,file='fname1_re1.txt'）
open（2,file='fname2_re2.txt'）
read（1,*）nn,ne,nm,nc
nf=2
nd=8
nfd=nf*nd
n=nn*nf
do50i=1,nn
50read（1,*）k,x（i）,y（i）,p（2*i-1）,p（2*i）
read（1,*）（（jz（i,j）,j=1,2）,i=1,nc）
do100i=1,ne
read（1,*）ie,（mea（i,j）,j=1,9）
100continue
read（1,*）（（aeu（i,j）,j=1,3）,i=1,nm）
!
400format（4i5）
!
405format（i5,4f15.6）
!
410format（5x,i5,5x,i5）
!
415format（5x,i5,5x,9i5）
!
420format（3f15.6）
write（2,460）nn,ne,nm,nc
write（2,465）（i,x（i）,y（i）,p（2*i-1）,p（2*i）,i=1,nn）
write（2,470）（（jz（i,j）,j=1,2）,i=1,nc）
write（2,475）（i,（mea（i,j）,j=1,9）,i=1,ne）
write（2,480）（i,（aeu（i,j）,j=1,3）,i=1,nm）
460format（/5x,'总的节点数，单元数，材料数，约束数'//（10x,4i5））
465format（/5x,'节点号',5x,'横坐标',5x,'纵坐标',5x,'节点载荷'//（5x,i5,4f15.6））
470format（/5x,'约束节点号及节点类型'//（5x,i5,5x,i5））
475format（/'节点编号矩阵：
单元号',15x,'八个节点的整体编号',12x,'材料号'//（10x,i5,8i5,5x,i5））
480format（/3x,'材料的材料号',5x,'杨氏模量',10x,'泊松比',10x,'厚度'//（5x,i5,5x,3e15.6））
callfld（nt,ld）
do500i=1,nt
500a（i）=0
do600ie=1,ne
callke（ie,ake）
callfis（ie,is）
do560i=1,nfd
do560j=1,nfd
if（is（i）-is（j））560,520,520
520ni=is（i）
ij=ld（ni）-ni+is（j）
a（ij）=a（ij）+ake（i,j）
560continue
600continue
callfcc（nt,ld,a）
calldecom（nt,a,ld）
!
write（*,*）'kakka'
write（2,850）（i,p（2*i-1）,p（2*i）,i=1,nn）
!
write（*,*）'kakka'
callstr
!
write（*,*）'kakka'
850format（//25x,'节点位移'//5x,'nodeno.',12x,'u',12x,'v',/（7x,i5,2e15.7））
end
!
c应力输出子程序
subroutinestr
common/sol1/nn,ne,nm,nc
common/shu3/aeu（10,3）
common/shu4/p（1000）
common/bmdmp/bmatx（3,16）,dmatx（3,3）
common/shu1/mea（200,9）
dimensions（3,16）,st（3）,ve（16）,is（16）
dimensionss（8）,tt（8）
datass/-1,1,1,-1,0,1,0,-1/
datatt/-1,-1,1,1.,-1,0,1,0/
callzero（3,3,dmatx）
callzero（3,16,bmatx）
do800ie=1,ne
write（2,850）ie
nm1=mea（ie,9）
yong=aeu（nm1,1）
poiss=aeu（nm1,2）
thick=aeu（nm1,3）
callmodps（yong,poiss）
callfis（ie,is）
do750i=1,16
ni=is（i）
ve（i）=p（ni）
750continue
do300i=1,8
exisp=ss（i）
etasp=tt（i）
k=mea（ie,i）
callsfr2（exisp,etasp）
calljacob2（ie,djacb）
callbmatps
calldot（3,3,16,dmatx,bmatx,s）
calldot（3,16,1,s,ve,st）
write（2,880）i,k,st
（1）,st
（2）,st（3）
c1=（st
（1）+st
（2））/2.0
c2=sqrt（（st
（1）-st
（2））**2/4.0+st（3）**2）
psi=c1+c2
psj=c1-c2
s1=amax1（psi,psj,0.0）
s3=amin1（psi,psj,0.0）
s2=st
（1）+st
（2）-s1-s3
!
write（*,*）s1,s2,s3
c3=（s1-s2）**2+（s2-s3）**2+（s1-s3）**2
seq4=sqrt（c3/2.0）
!
write（2,*）seq4
300continue
800continue
880format（i5,5x,i5,5x,f15.5,5x,f15.5,5x,f15.5）
850format（//10x,'第',i3,'个单元的节点应力'/2x,'局部编号',2x,'整体编号',5x,'u方向应力',12x,'v方向应力',12x,'剪应力'）
return
end
!
c这是一个形成局部编号与整体编号对应关系is数组的子程序
SUBROUTINEFIS（IE,IS）
common/shu1/mea（200,9）
common/sol2/nf,nd,n,nfd
dimensionis（nfd）
DO150ID=1,ND
DO100JF=1,NF
I1=（ID-1）*NF+JF
IS（I1）=（MEA（IE,ID）-1）*NF+JF
100CONTINUE
150CONTINUE
RETURN
END
!
c这是一个形成变带宽存储指示数组ld的子程序
SUBROUTINEFLD（NT,LD）
common/sol1/nn,ne,nm,nc
common/sol2/nf,nd,n,nfd
common/shu1/mea（200,9）
dimensionld（n）
LD
（1）=1
DO300K=1,NN
NMIN=1000
DO200I=1,NE
DO150J=1,ND
IF（MEA（I,J）.NE.K）GOTO150
DO100L=1,ND
IF（MEA（I,L）.LT.NMIN）NMIN=MEA（I,L）
100CONTINUE
150CONTINUE
200CONTINUE
DO250I=1,NF
J=（K-1）*NF+I
IF（J.NE.1）LD（J）=LD（J-1）+（K-NMIN）*NF+I
250CONTINUE
300CONTINUE
NT=LD（N）
RETURN
END
!
c这是一个处理约束的子程序
SUBROUTINEFCC（NT,LD,A）
common/sol1/nn,ne,nm,nc
common/sol2/nf,nd,n,nfd
common/shu5/jz（1000,2）
dimensionld（n）,a（nt）
DO350K=1,NC
I=JZ（K,1）
J=JZ（K,2）
NI=（I-1）*NF+J
NJ=LD（NI）
A（NJ）=1.0E25
350CONTINUE
RETURN
END
!
c这是一个解线性方程组的子程序
SUBROUTINEDECOM（NT,A,LD）
common/shu4/p（1000）
common/sol2/nf,nd,n,nfd
dimensiona（nt）,ld（n）
DO20I=1,N
LDI=LD（I）
IF（I.EQ.1）GOTO10
IO=LDI-I
IM1=I-1
MI=1-IO+LD（IM1）
IF（MI.GT.IM1）GOTO10
DO8J=MI,IM1
JO=LD（J）-J
IJ=IO+J
IF（J.EQ.1）GOTO6
JM1=J-1
MJ=1-JO+LD（JM1）
MIJ=MAX0（MI,MJ）
IF（MIJ.GT.JM1）GOTO6
S=0.0
DO2K=MIJ,JM1
IK=IO+K
JK=JO+K
S=S+A（IK）*A（JK）
2CONTINUE
A（IJ）=A（IJ）-S
6p（I）=p（I）-A（IJ）*p（J）
8CONTINUE
10ALDI=A（LDI）
IF（I.EQ.1.OR.MI.GT.IM1）GOTO13
DO12J=MI,IM1
IJ=IO+J
LDJ=LD（J）
S=A（IJ）
A（IJ）=S/A（LDJ）
ALDI=ALDI-S*A（IJ）
12CONTINUE
13A（LDI）=ALDI
p（I）=p（I）/ALDI
20CONTINUE
DO40LDI=2,N
I=N-LDI+2
IO=LD（I）-I
MI=1-IO+LD（I-1）
IM1=I-1
IF（MI.GT.IM1）GOTO40
DO30J=MI,IM1
IJ=IO+J
p（J）=p（J）-A（IJ）*p（I）
30CONTINUE
40CONTINUE
RETURN
END
!
c这是一个形成单元刚度矩阵的子程序
subroutineke（ie,ake）
dimensionposgp
（2）
common/sol1/nn,ne,nm,nc
common/bmdmp/bmatx（3,16）,dmatx（3,3）
common/shu3/aeu（10,3）
common/shu1/mea（200,9）
common/shu2/x（1000）,y（1000）
dimensionbt（16,3）,bd（16,3）
dimensionake（16,16）,ake1（16,16）
posgp
（1）=-0.5773502692
posgp
（2）=0.5773502692
nm1=mea（ie,9）
young=aeu（nm1,1）
poiss=aeu（nm1,2）
thick=aeu（nm1,3）
callmodps（young,poiss）
callzero（16,16,ake）
do80iagus=1,2
do80jagus=1,2
exisp=posgp（iagus）
etasp=posgp（jagus）
callsfr2（exisp,etasp）
calljacob2（ie,djacb）
callbmatps
calltran（3,16,bmatx,bt）
calldot（16,3,3,bt,dmatx,bd）
calldot（16,3,16,bd,bmatx,ake1）
do70i=1,16
do70j=1,16
70ake（i,j）=ake（i,j）+djacb*thick*ake1（i,j）
80continue
!
write（2,*）ie
!
do90i=1,16
!
write（*,*）ie
!
write（2,10）（ake（i,j）,j=1,16）
!
90continue
!
10format（16e15.3）
return
end
!
c这是一个形成弹性矩阵的子程序
subroutinemodps（y,p）
common/bmdmp/bmatx（3,16）,dmatx（3,3）
callzero（3,3,dmatx）
!
write（*,*）y,p
const=y/（1.0-p*p）
dmatx（1,1）=const
dmatx（2,2）=const
dmatx（1,2）=const*p
dmatx（2,1）=dmatx（1,2）
dmatx（3,3）=const*（1.0-p）/2.0
!
do100i=1,3
!
write（*,*）（dmatx（i,j）,j=1,3）
!
100continue
return
end
!
c这是一个计算形函数当前积分点及形函数对局部坐标的导数值的子程序
subroutinesfr2（s,t）
common/sd/sh（8）,deriv（2,8）
s2=s*2.0
t2=t*2.0
ss=s*s
tt=t*t
st=s*t
stt=s*t*t
sst=s*s*t
st2=s*t*2.0
sh
（1）=（-1.0+st+ss+tt-sst-stt）/4.0
sh
（2）=（-1.0-st+ss+tt-sst+stt）/4.0
sh（3）=（-1.0+st+ss+tt+sst+stt）/4.0
sh（4）=（-1.0-st+ss+tt+sst-stt）/4.0
sh（5）=（1.0-t-ss+sst）/2.0
sh（6）=（1.0+s-tt-stt）/2.0
sh（7）=（1.0+t-ss-sst）/2.0
sh（8）=（1.0-s-tt+stt）/2.0
deriv（1,1）=（t+s2-st2-tt）/4.0
deriv（1,2）=（-t+s2-st2+tt）/4.0
deriv（1,3）=（t+s2+st2+tt）/4.0
deriv（1,4）=（-t+s2+st2-tt）/4.0
deriv（1,5）=-s+st
deriv（1,6）=（1.0-tt）/2.0
deriv（1,7）=-s-st
deriv（1,8）=（-1.0+tt）/2.0
deriv（2,1）=（s+t2-ss-st2）/4.0
deriv（2,2）=（-s+t2-ss+st2）/4.0
deriv（2,3）=（s+t2+ss+st2）/4.0
deriv（2,4）=（-s+t2+ss-st2）/4.0
deriv（2,5）=（-1.0+ss）/2.0
deriv（2,6）=-t-st
deriv（2,7）=（1.0-ss）/2.0
deriv（2,8）=-t+st
return
end
!
c这是一个形成jacobi矩阵的子程序
subroutinejacob2（ie,djacb）
dimensionxjacm（2,2）,xjaci（2,2）
common/sd/sh（8）,deriv（2,8）
common/car/cartd（2,8）
common/shu1/mea（200,9）
common/shu2/x（1000）,y（1000）
do200j=1,2
xjacm（1,j）=0.0
xjacm（2,j）=0.0
do200k=1,8
if（j.eq.1）then
xjacm（1,j）=xjacm（1,j）+deriv（1,k）*x（mea（ie,k））
xjacm（2,j）=xjacm（2,j）+deriv（2,k）*x（mea（ie,k））
else
xjacm（1,j）=xjacm（1,j）+deriv（1,k）*y（mea（ie,k））
xjacm（2,j）=xjacm（2,j）+deriv（2,k）*y（mea（ie,k））
endif
200continue
djacb=xjacm（1,1）*xjacm（2,2）-xjacm（1,2）*xjacm（2,1）
!
write（*,*）xjacm（1,1）,xjacm（1,2）,xjacm（2,1）,xjacm（2,2）
!
write（*,*）ie
!
write（*,*）djacb
if（djacb.lt.1.e-6）then
write（6,6100）ie
6100format（'单元编号为:
',i4,'的jcobi行列式的值小于或等于零!
'）
stop'程序运行被中断于子程序jacob2'
endif
xjaci（1,1）=xjacm（2,2）/djacb
xjaci（2,2）=xjacm（1,1）/djacb
xjaci（1,2）=-xjacm（1,2）/djacb
xjaci（2,1）=-xjacm（2,1）/djacb
do300i=1,2
do300k=1,8
cartd（i,k）=0.0
do300j=1,2
cartd（i,k）=cartd（i,k）+xjaci（i,j）*deriv（j,k）
300continue
return
end
!
c这是一个形成应变矩阵的子程序
subroutinebmatps
common/car/cartd（2,8）
common/bmdmp/bmatx（3,16）,dmatx（3,3）
callzero（3,16,bmatx）
ngash=0
do200i=1,8
mgash=ngash+1
ngash=mgash+1
bmatx（1,mgash）=cartd（1,i）
bmatx（1,ngash）=0.0
bmatx（2,mgash）=0.0
bmatx（2,ngash）=cartd（2,i）
bmatx（3,mgash）=cartd（2,i）
bmatx（3,ngash）=cartd（1,i）
200continue
return
end
SUBROUTINEZERO（M,N,A）
DIMENSIONA（M,N）
DO10I=1,M
DO10J=1,N
A（I,J）=0.0
10CONTINUE
RETURN
END
SUBROUTINETRAN（M,N,A,B）
DIMENSIONA（M,N）,B（N,M）
DO20I=1,M
DO20J=1,N
B（J,I）=A（I,J）
20CONTINUE
RETURN
END
SUBROUTINEDOT（M,N,L,A,B,C）
DIMENSIONA（M,N）,B（N,L）,C（M,L）
DO50I=1,M
DO50J=1,L
C（I,J）=0.0
DO40K=1,N
C（I,J）=C（I,J）+A（I,K）*B（K,J）
40CONTINUE
50CONTINUE
RETURN
END
升级会员 