平面四边形四节点等参单元Fortran源程序.doc

C ************************************************
C * FINITE ELEMENT PROGRAM *
C * FOR Two DIMENSIONAL ELASticity PROBLEM *
C * WITH 4 NODE *
C ************************************************
PROGRAM ELASTICITY
character*32 dat,cch
DIMENSION SK(80000),COOR(2,300),AE(4,11),MEL(5,200),
& WG(4),JR(2,300),MA(600),R(600),iew(30),STRE(3,200)
MON /CMN1/ NP,NE,NM,NR
MON /CMN2/ N,MX,NH
MON /CMN3/ RF(8),SKE(8,8),NN(8)
WRITE(*,*)'PLEASE ENTER INPUT '
READ(*,'(A)')DAT
OPEN(4,'OLD')
OPEN(7,FILE='OUT',STATUS='UNKNOWN')
READ(4,*)NP,NE,NM,NR
WRITE(7,'(A,I6)')'NUMBER OF NODE---------------------NP=',np
WRITE(7,'(A,I6)')'NUMBER OF ELEMENT------------------NE=',ne
WRITE(7,'(A,I6)')'NUMBER OF MATERIAL-----------------NM=',nm
WRITE(7,'(A,I6)')'NUMBER OF surporting---------------NC=',Nr
CALL INPUT (JR,COOR,AE,MEL)
CALL CBAND (MA,JR,MEL)
DO I=1,NH
SK(I)=0、0
enddo
CALL SK0(SK,MEL,COOR,JR,MA,AE)
do I=1,N
R(I)=0、0
enddo
pause 'aaa'
stop
READ(4,*)NCP,NBE,iz
WRITE(*,'(5i8)')NCP,NBE,iz
WRITE(7,'(5i8)')NCP,NBE,iz
IF(NCP、GT、0)CALL CONCR(NCP,R,JR)
IF(NBE、GT、0) CALL BODYR(NBE,R,MEL,COOR,JR,AE)
IF(iz、GT、0)then
do jj=1,iz
READ (4,*)Js,nse,(WG(I),I=1,4)
read(4,*)(iew(m),m=1,nse)
CALL FACER(iew,NSE,R,MEL,COOR,JR,WG)
enddo
endif
CALL DECOP (SK,MA)
CALL FOBA (SK,MA,R)
CALL OUTDISP(NP,R,JR)
CALL STRESS (COOR,MEL,JR,AE,R,STRE)
WRITE(7,'(A)')' PROGRAM SAFF HAS BEEN ENDED'
WRITE(*,'(A)')' PROGRAM SAFF HAS BEEN ENDED'
STOP
c RETURN
END
C *********************************************
SUBROUTINE INPUT (JR,COOR,AE,MEL)
DIMENSION JR(2,*),COOR(2,*),AE(4,*),MEL(5,*)
MON /CMN1/ NP,NE,NM,NR
MON /CMN2/ N,MX,NH
DO 70 I=1,NP
READ(4,*) IP,X,Y
COOR(1,IP)=X
COOR(2,IP)=Y
70 CONTINUE
DO 11 J=1,NE
READ(4,*)NEE,NME,(MEL(I,NEE),I=1,4)
MEL(5,NEE)=NME
11 CONTINUE
DO 10 I=1,NP
DO 10 J=1,2
10 JR(J,I)=1
DO 20 I=1,NR
READ(4,*) IP,IX,IY
JR(1,IP)=IX
JR(2,IP)=IY
20 CONTINUE
N=0
DO 30 I=1,NP
DO 30 J=1,2
IF (JR(J,I)) 30,30,25
25 N=N+1
JR(J,I)=N
30 CONTINUE
DO 55 J=1,NM
READ (4,*)JJ,(AE(I,JJ),I=1,4)
WRITE(*,910) JJ,(AE(I,JJ),I=1,4)
55 CONTINUE
910 FORMAT (/20X,'MATERIAL PROPERTIES'/(3X,I5,4(1x,E8、3)))
RETURN
END
C **********************************************
SUBROUTINE CBAND (MA,JR,MEL)
DIMENSION MA(*),JR(2,*),MEL(5,*),NN(8)
MON /CMN1/ NP,NE,NM,NR
MON /CMN2/ N,MX,NH
DO 65 I=1,N
65 MA(I)=0
DO 90 IE=1,NE
DO 75 K=1,4
IEK=MEL(K,IE)
DO 95 M=1,2
JJ=2*(K-1)+M
NN(JJ)=JR(M,IEK)
95 CONTINUE
75 CONTINUE
L=N
DO 80 I=1,2*4
NNI=NN(I)
IF(NNI、EQ、0) GO TO 80
IF(NNI、LT、L) L=NNI
80 CONTINUE
DO 85 M=1,2*4
JP=NN(M)
IF(JP、EQ、0) GO TO 85
JPL=JP-L+1
IF(JPL、GT、MA(JP)) MA(JP)=JPL
85 CONTINUE
90 CONTINUE
MX=0
MA(1)=1
DO 10 I=2,N
IF(MA(I)、GT、MX) MX=MA(I)
MA(I)=MA(I)+MA(I-1)
10 CONTINUE
NH=MA(N)
WRITE(7,'(A,I8)')'TOTAL DEGREES OF FREEDOM-----------N= ',N
WRITE(7,'(A,I8)')'MAX-SEMI-BANDWIDTH-----------------MX=',MX
WRITE(7,'(A,I8)')'TOTAL-STORAGE----------------------NH=',NH
500 FORMAT (/5X,'FREEDOM N='
*,I5,3X,'SEMI-BANDWI、 MX=',I5,3X,
* 'STORAGE NH=',I7)
RETURN
END
C **********************************************
SUBROUTINE SK0(SK,MEL,COOR,JR,MA,AE)
DIMENSION SK(*),MEL(5,*),COOR(2,*),JR(2,*),MA(*),
* AE(4,*),XYZ(2,4),iven(4)
MON /CMN1/ NP,NE,NM,NR
MON /CMN2/ N,MX,NH
MON /CMN3/ RF(8),SKE(8,8),NN(8)
MON /CMN4/ NEE,NME
MON /GAUSS/ RSTG(3),H(3)
H(1)=0、5555555555555560
H(2)=0、8888888888888890
H(3)=H(1)
RSTG(1)=-0、7745966692414830
RSTG(2)=0、00
RSTG(3)=-RSTG(1)
DO 10 IE=1,NE
NEE=IE
NME=MEL(5,IE)
DO 75 K=1,4
IEK=MEL(K,IE)
iven(k)=IEK
DO 95 M=1,2
JJ=2*(K-1)+M
NN(JJ)=JR(M,IEK)
95 XYZ(M,K)=COOR(M,IEK)
75 CONTINUE
CALL STIF(XYZ,AE,iven)
DO 60 I=1,8
DO 60 J=1,8
II=NN(I)
JJ=NN(J)
IF ((JJ、EQ、0)、OR、(II、LT、JJ)) GO TO 60
JN=MA(II)-(II-JJ)
SK(JN)=SK(JN)+SKE(I,J)
60 CONTINUE
70 CONTINUE
write(7,1111) ((ske(i,j),j=1,8),i=1,8)
1111 format(2x,8f12、2)
10 CONTINUE
RETURN
END
C *********************************************
SUBROUTINE STIF(XYZ,AE,iven)
DIMENSION AE(4,*),DNX(2,4),XYZ(2,*),iven(*),
* RJAC(2,2)
MON /CMN1/ NP,NE,NM,NR
MON /CMN2/ N,MX,NH
MON /CMN3/ RF(8),SKE(8,8),NN(8)
MON /CMN4/ NEE,NME
MON /GAUSS/ RSTG(3),H(3)
DO 40 I=1,8
RF(I)=0、00
DO 30 J=1,8
SKE(I,J)=0、00
30 CONTINUE
40 CONTINUE
E=AE(1,NME)
U=AE(2,NME)
GAMA=AE(3,NME)
D1=E*(1、00-U)/((1、00+U)*(1、00-2、00*U))
D2=E*U/((1、00+U)*(1、00-2、00*U))
D3=E*0、50/(1、00+U)
DO 120 I=1,4
II=2*(I-1)
I1=II+1
I2=II+2
DO 115 J=1,4
JJ=2*(J-1)
J1=JJ+1
J2=JJ+2
DXX=0
DXY=0
DYX=0
DYY=0
DO 99 IS=1,3
S=RSTG(IS)
SH=H(IS)
DO 98 IR=1,3
R=RSTG(IR)
RH=H(IR)
CALL FDNX (XYZ,DNX,DET,R,S,RJAC,iven,NEE)
DNIX=DNX(1,I)
DNIY=DNX(2,I)
DNJX=DNX(1,J)
DNJY=DNX(2,J)
DXX=DXX+DNIX*DNJX*DET*RH*SH
DXY=DXY+DNIX*DNJY*DET*RH*SH
DYX=DYX+DNIY*DNJX*DET*RH*SH
DYY=DYY+DNIY*DNJY*DET*RH*SH
98 CONTINUE
99 CONTINUE
SKE(I1,J1)=DXX*D1+DYY*D3
SKE(I2,J2)=DYY*D1+DXX*D3
SKE(I1,J2)=DXY*D2+DYX*D3
SKE(I2,J1)=DYX*D2+DXY*D3
115 CONTINUE
120 CONTINUE
RETURN
END
C *********************************************
SUBROUTINE CONCR(NCP,R,JR)
DIMENSION R(*),JR(2,*),XYZ(2)
DO 100 I=1,NCP
READ (4,*) IP,PX,PY
XYZ(1)=PX
XYZ(2)=PY
DO 95 J=1,2
L=JR(J,IP)
IF(L、EQ、0) GO TO 95
R(L)=R(L)+XYZ(J)
95 CONTINUE
100 CONTINUE
RETURN
END
C **********************************************
SUBROUTINE BODYR(NBE,R,MEL,COOR,JR,AE)
DIMENSION R(*),MEL(5,*),COOR(2,*),JR(2,*),
& AE(4,*),XYZ(2,4),iven(4)
MON /CMN1/ NP,NE,NM,NR
MON /CMN2/ N,MX,NH
MON /CMN3/ RF(8),SKE(8,8),NN(8)
MON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)
MON /GAUSS/ RSTG(3),H(3)
H(1)=1、0
H(2)=1、0
RSTG(1)=-0、57735
RSTG(2)=-RSTG(1)
DO 10 IE=1,NBE
DO I=1,8
RF(I)=0、00
ENDDO
c READ(4,*)NEE
NEE=ie
NME=MEL(5,NEE)
GAMA=AE(3,NME)
DO 75 K=1,4
IEK=MEL(K,NEE)
iven(k)=iek
DO 95 M=1,2
JJ=2*(K-1)+M
NN(JJ)=JR(M,IEK)
95 XYZ(M,K)=COOR(M,IEK)
75 CONTINUE
DO 99 IS=1,2
S=RSTG(IS)
SH=H(IS)
DO 98 IR=1,2
RR=RSTG(IR)
RH=H(IR)
CALL FUN8 (XYZ,RR,S,DET)
DO 30 I=1,4
J=2*I
RF(J)=RF(J)-FUN(I)*RH*SH*DET*GAMA
30 CONTINUE
98 CONTINUE
99 CONTINUE
CALL ASLOAD (R)
10 CONTINUE
RETURN
END
C *********************************************
SUBROUTINE FACER(iew,NSE,R,MEL,COOR,JR,WG)
DIMENSION R(*),MEL(5,*),COOR(2,*),JR(2,*),wg(*)
* ,XYZ(2,4),iew(*),PR(2)
MON /CMN1/ NP,NE,NM,NR
MON /CMN2/ N,MX,NH
MON /CMN3/ RF(8),SKE(8,8),NN(8)
MON /CMN4/ NEE,NME
MON /GAUSS/ RSTG(3),H(3)
H(1)=1、0
H(2)=1、0
RSTG(1)=-0、57735
RSTG(2)=-RSTG(1)
nwf=0
nnf=0
ir=wg(1)+0、1
if(ir、eq、1)nwf=1
if(ir、eq、2)nnf=1
DO 510 IE=1,NSE
DO I=1,8
RF(I)=0、00
ENDDO
nee=iew(ie)
DO 575 K=1,4
IEK=MEL(K,NEE)
DO 595 M=1,2
JJ=2*(K-1)+M
NN(JJ)=JR(M,IEK)
595 XYZ(M,K)=COOR(M,IEK)
575 CONTINUE
IF(NWF、EQ、1) then
GAMA=WG(2)
Z0=WG(3)
NSU=WG(4)+0、1
CALL SURLOD (NSU,XYZ,PR,Z0,GAMA,1)
endif
IF(NNF、EQ、1) then
q=WG(2)
NSU=WG(4)+0、1
do j=1,2
PR(J)=q
enddo
CALL SURLOD (NSU,XYZ,PR,Z0,GAMA,2)
endif
CALL ASLOAD (R)
510 CONTINUE
RETURN
END
C *********************************************
SUBROUTINE SURLOD (NSU,XYZ,PR,Z0,GAMA,NSI)
DIMENSION XYZ(2,*),RST(3),PR(2),KCRD(4),KFACE(2,4),
& FVAL(4),NODES(2),FACT(4)
MON /CMN1/ NP,NE,NM,NR
MON /CMN2/ N,MX,NH
MON /CMN3/ RF(8),SKE(8,8),NN(8)
MON /CMN4/ NEE,NME
MON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)
MON /GAUSS/ RSTG(3),H(3)
DATA KCRD/1,1,2,2/
DATA KFACE/1, 4,
* 2, 3,
* 1, 2,
* 4, 3/
DATA FVAL/-1、00,1、00,-1、00,1、00/
FACT(1)=1、0
FACT(2)=-1、0
FACT(3)=-1、0
FACT(4)=1、0
FACTNUS=FACT(NSU)
DO I=1,2
J=KFACE(I,NSU)
NODES(I)=J
ENDDO
IF (NSI、EQ、1) THEN
DO I=1,2
J=NODES(I)
Z=Z0-XYZ(2,J)
PR(I)=0、00
IF (Z、GT、0、00) PR(I)=Z*GAMA
ENDDO
ENDIF
ML=KCRD(NSU)
IF(ML、EQ、1)MM=2
IF(ML、EQ、2)MM=1
RST(ML)=FVAL(NSU)
DO 70 LX=1,2
RST(MM)=RSTG(LX)
CALL FUN8 (XYZ,RST(1),RST(2),DET)
PXYZ=0、00
DO 25 I=1,2
J=NODES(I)
PXYZ=PXYZ+FUN(J)*PR(I)
25 CONTINUE
A1=XJAC(MM,2)
A2=-XJAC(MM,1)
30 DO 60 I=1,2
J=NODES(I)
K2=2*J
K1=K2-1
Q=PXYZ*FUN(J)*H(LX)*FACTNUS
RF(K1)=RF(K1)+Q*A1
RF(K2)=RF(K2)+Q*A2
60 CONTINUE
70 CONTINUE
RETURN
END
C *********************************************
SUBROUTINE ASLOAD (R)
DIMENSION R(*)
MON /CMN1/ NP,NE,NM,NR
MON /CMN3/ RF(8),SKE(8,8),NN(8)
DO 20 I=1,8
L=NN(I)
IF (L、EQ、0) GO TO 20
R(L)=R(L)+RF(I)
20 CONTINUE
RETURN
END
C ***********************************************
SUBROUTINE DECOP (SK,MA)
DIMENSION SK(*),MA(*)
MON /CMN2/ N,MX,NH
DO 50 I=2,N
L=I-MA(I)+MA(I-1)+1
K=I-1
L1=L+1
IF (L1、GT、K) GO TO 30
DO 20 J=L1,K
IJ=MA(I)-I+J
M=J-MA(J)+MA(J-1)+1
IF (L、GT、M) M=L
MP=J-1
IF (M、GT、MP) GO TO 20
DO 10 LP=M,MP
IP=MA(I)-I+LP
JP=MA(J)-J+LP
SK(IJ)=SK(IJ)-SK(IP)*SK(JP)
10 CONTINUE
20 CONTINUE
30 IF (L、GT、K) GO TO 50
DO 40 LP=L,K
IP=MA(I)-I+LP
LPP=MA(LP)
SK(IP)=SK(IP)/SK(LPP)
II=MA(I)
SK(II)=SK(II)-SK(IP)*SK(IP)*SK(LPP)
40 CONTINUE
50 CONTINUE
RETURN
END
C *************************************************************
SUBROUTINE FOBA (SK,MA,R)
DIMENSION SK(*),MA(*),R(*)
MON /CMN2/ N,MX,NH
DO 10 I=2,N
L=I-MA(I)+MA(I-1)+1
K=I-1
IF (L、GT、K) GO TO 10
DO 5 LP=L,K
IP=MA(I)-I+LP
R(I)=R(I)-SK(IP)*R(LP)
5 CONTINUE
10 CONTINUE
DO 20 I=1,N
II=MA(I)
45 R(I)=R(I)/SK(II)
20 CONTINUE
DO 30 J1=2,N
I=2+N-J1
L=I-MA(I)+MA(I-1)+1
K=I-1
IF (L、GT、K) GO TO 30
DO 25 J=L,K
IJ=MA(I)-I+J
55 R(J)=R(J)-SK(IJ)*R(I)
25 CONTINUE
30 CONTINUE
RETURN
END
C *****************************************************************
SUBROUTINE STRESS(COOR,MEL,JR,AE,R,STRE)
DIMENSION XYZ(2,4),DNX(2,4),AE(4,*),STRE(3,*),
& COOR(2,*),MEL(5,*),JR(2,*),RJAC(2,2),SIG(3),
& B(3,8),R(*),iven(4)
MON /CMN1/ NP,NE,NM,NR
MON /CMN3/ RF(8),SKE(8,8),NN(8)
MON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)
DO 106 IE=1,NE
NME=MEL(5,IE)
DO 300 K=1,4
IEK=MEL(K,IE)
DO 310 M=1,2
310 XYZ(M,K)=COOR(M,IEK)
DO 320 M=1,2
JRR=2*(K-1)+M
320 NN(JRR)=JR(M,IEK)
300 CONTINUE
E=AE(1,NME)
U=AE(2,NME)
D1=E*(1、00-U)/((1、00+U)*(1、00-2、00*U))
D2=E*U/((1、00+U)*(1、00-2、00*U))
D3=0、50*E/(1、00+U)
SS=0、0
RR=0、0
CALL FDNX (XYZ,DNX,DET,RR,SS,RJAC,iven,IE)
DO 30 I=1,4
II=2*(I-1)
J1=II+1
J2=II+2
BI=DNX(1,I)
CI=DNX(2,I)
B(1,J1)=BI
B(2,J1)=0、
B(3,J1)=CI
B(1,J2)=0、
B(2,J2)=CI
B(3,J2)=BI
30 CONTINUE
DO 55 II=1,3
SIG(II)=0、00
55 CONTINUE
DO 70 K=1,8
NA=NN(K)
IF (NA、EQ、0) GO TO 70
DO 60 L=1,3
SIG(L)=SIG(L)+B(L,K)*R(NA)
60 CONTINUE
70 CONTINUE
SX=D1*SIG(1)+D2*SIG(2)
SY=D2*SIG(1)+D1*SIG(2)
SXY=D3*SIG(3)
STRE(1,IE)=SX
STRE(2,IE)=SY
STRE(3,IE)=SXY
106 CONTINUE
CALL OUTSTRE(NE,STRE)
RETURN
END
C *********************************************
SUBROUTINE FDNX (XYZ,DNX,DET,R,S,RJAC,iven,NEE)
DIMENSION XYZ(2,*),DNX(2,*),RJAC(2,2),iven(*)
MON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)
CALL FUN8 (XYZ,R,S,DET)
IF (DET、LT、1、0E-5)THEN
WRITE(7,600) NEE,R,S,det
WRITE(7,*) (iven(m),m=1,4)
STOP
ENDIF
REC=1、00/DET
RJAC(1,1)=REC*XJAC(2,2)
RJAC(2,2)=REC*XJAC(1,1)
RJAC(2,1)=-REC*XJAC(2,1)
RJAC(1,2)=-REC*XJAC(1,2)
DO 30 K=1,4
DO 20 I=1,2
DNX(I,K)=0、
DO 25 M=1,2
DNX(I,K)=DNX(I,K)+RJAC(I,M)*PN(M,K)
25 CONTINUE
20 CONTINUE
30 CONTINUE
600 FORMAT (1X,'ERR0R*** NEGTIVE OR ZERO '
* 'JACOBIAN DETERMINANT FOR '
* 'ELEMENT'/'ELE、=',I5,' R=',F10、5,6X,'S=',F10、5,
* 'det=',f12、5)
RETURN
END
C *********************************************
SUBROUTINE FUN8 (XYZ,R,S,DET)
DIMENSION XYZ(2,*),XI(4),ETA(4)
MON /CMN5/ FUN(4),PN(2,4),XJAC(2,2)
DATA XI/-1、0,1、0,1、0,-1、0/
DATA ETA/-1、0,-1、0,1、0,1、0/
DO 10 I=1,4
G1=(1、0+XI(I)*R)
G2=(1、0+ETA(I)*S)
FUN(I)=0、25*G1*G2
PN(1,I)=0、25*XI(I)*G2
PN(2,I)=0、25*ETA(I)*G1
10 CONTINUE
DO 80 I=1,2
DO 75 J=1,2
DET=0、00
DO 70 K=1,4
DET=DET+PN(I,K)*XYZ(J,K)
70 CONTINUE
XJAC(I,J)=DET
75 CONTINUE
80 CONTINUE
DET=XJAC(1,1)*XJAC(2,2)
* -XJAC(2,1)*XJAC(1,2)
RETURN
END
C *********************************************