構造力学 構造工学特論
FORTRANプログラム
平面応力 平面ひずみ（剛性マトリックス法）
FEM.FOR
C ******** STRESS ANALYSIS OF 2-DIMENSIONAL PROBLEM ********
C KIKUTI,SUGIDATE AND MIYAMOTO 1983,7,7
C
DIMENSION F(80),INDEX(80),DIS(80),REAC(80),NOKX(20),NOKY(20)
1 ,WD(20),WS(20),W(20,4),NONGL(20),ANGL(20)
COMMON KAKOM(50,3),X(40),Y(40),T(50),SM(6,6),TSM(80,80),STR(3,6)
2 ,BI(80,80)
C
C ******** DATA INPUT ********
C
1000 READ(5,10,END=3000) NSTRE
IF(NSTRE) 4,5,4
4 WRITE(6,6)
6 FORMAT(1H0,30H******** PLANE STRESS ********)
GO TO 7
5 WRITE(6,8)
8 FORMAT(1H0,30H******** PLANE STRAIN ********)
7 READ(5,10) NODT,NELT,KOX,KOY,NF
10 FORMAT(12I5)
DO 15 I=1,NELT
15 READ(5,11)NO,(KAKOM(I,J),J=1,3),T(I)
11 FORMAT(4I5,F10.0)
DO 20 I=1,NODT
20 READ(5,25) NO,X(I),Y(I)
25 FORMAT(I5,2F10.0)
READ(5,26) E,PO
26 FORMAT(2F10.0)
READ(5,10) (NOKX(I),I=1,KOX)
READ(5,10) (NOKY(I),I=1,KOY)
READ(5,10) NANGL
IF(NANGL) 27,90,27
27 READ(5,28) (NONGL(K),ANGL(K),K=1,NANGL)
28 FORMAT(I5,F10.0)
90 NT=NODT*2
DO 30 I=1,NT
30 F(I)=0.
DO 35 I=1,NF
READ(5,29) NO,FX,FY
29 FORMAT(I5,2F10.3)
F(2*NO-1)=FX
35 F(2*NO)=FY
C
C ******** DATA PRINT ********
C
WRITE(6,45) NODT,NELT,KOX,KOY,NF
45 FORMAT(1H0,3X,5HNODT=,I3,3X,5HNELT=,I3,3X,4HKOX=,I3,3X,4HKOY=,I3
& ,3X,3HNF=,I3)
WRITE(6,50)
50 FORMAT(1H0,1X,48H********* NODE NO. AND PLANE THICKNESS *********)
DO 55 I=1,NELT
55 WRITE(6,60) I,(KAKOM(I,J),J=1,3),T(I)
60 FORMAT(1H0,I5,5X,3I5,F10.3)
WRITE(6,65)
65 FORMAT(1H0,1X,35H********* ELASTIC MODULUS *********)
WRITE(6,70) E,PO
70 FORMAT(1H0,5X,2HE=,F10.0,5X,3HPO=,F8.3)
WRITE(6,72)
72 FORMAT(1H0,1X,48H********* FIXED NODAL POINT IN X-DIRECTION *****)
WRITE(6,74) (NOKX(I),I=1,KOX)
74 FORMAT(1H0,20I5)
WRITE(6,76)
76 FORMAT(1H0,1X,48H********* FIXED NODAL POINT IN Y-DIRECTION *****)
WRITE(6,74) (NOKY(K),K=1,KOY)
IF(NANGL) 700,700,201
201 WRITE(6,202)
202 FORMAT(1H0,43H********* DIAGONAL FIXED NODAL POINT ******)
WRITE(6,204) (NONGL(K),ANGL(K),K=1,NANGL)
204 FORMAT(1H0,I5,5X,F10.3)
700 WRITE(6,78)
78 FORMAT(1H0,1X,48H********* DIMENSION(X,Y) AND FORCE(FX,FY) ******)
WRITE(6,80)
80 FORMAT(1H0,14X,4HX(I),6X,4HY(I),6X,5HFX(I),5X,5HFY(I))
WRITE(6,82) (I,X(I),Y(I),F(2*I-1),F(2*I),I=1,NODT)
82 FORMAT(1H0,I5,5X,4F10.3)
C
C ******** SSENBLY OF ELEMENT STIFFNESS MATRIX *******
C
DO 100 I=1,NT
DO 100 J=1,NT
100 TSM(I,J)=0.
DO 110 NE=1,NELT
CALL PLANS(NSTRE,NE,E,PO,NODT,NELT,AREA)
WRITE(6,1100) NE,AREA
1100 FORMAT(1H0,5X,I5,5HAREA=,E15.5)
WRITE(6,2500)
2500 FORMAT(1H0,1X,40H****** ELEMENT STIFFNESS MATRIX ********)
DO 1200 I=1,6
1200 WRITE(6,1300) NE,(SM(I,J),J=1,6)
1300 FORMAT(1H0,I5,5X,6E15.5)
110 CALL ASMAT(NE,3,2,NT,NELT,6)
WRITE(6,1400)
1400 FORMAT(1H0,1X,48H***** ASSENBLY OF ELEMENT STIFFNESS MATRIX *****)
DO 1500 I=1,NT
1500 WRITE(6,1600) (TSM(I,J),J=1,NT)
1600 FORMAT(1H0,5X,8E15.5)
C
IF(NANGL) 93,92,93
93 DO 500 N=1,NANGL
ANGL(N)=ANGL(N)*3.1416/180.
W(N,1)=COS(ANGL(N))
W(N,2)=-SIN(ANGL(N))
W(N,3)=SIN(ANGL(N))
W(N,4)=COS(ANGL(N))
NI=NONGL(N)
DO 500 NJ=1,NODT
DO 520 I=1,2
DO 520 J=1,2
JA=(NT-1)*2+J
SM(I,J)=0.
DO 520 K=1,2
KA=(NI-1)*2+K
KI=(K-1)*2+I
520 SM(I,J)=SM(I,J)+W(N,KI)*TSM(KA,JA)
DO 500 I=1,2
IA=(NI-1)*2+I
DO 500 J=1,2
JA=(NJ-1)*2+J
500 TSM(IA,JA)=SM(I,J)
DO 530 N=1,NANGL
NJ=NONGL(N)
DO 530 NI=1,NODT
DO 540 I=1,2
IA=(NI-1)*2+I
DO 540 J=1,2
SM(I,J)=0.
DO 540 K=1,2
KA=(NJ-1)*2+K
KJ=(K-1)*2+J
540 SM (I,J)=SM(I,J)+TSM(IA,KA)*W(N,KJ)
DO 530 I=1,2
IA=(NI-1)*2+I
DO 530 J=1,2
JA=(NJ-1)*2+J
530 TSM(IA,JA)=SM(I,J)
C
C ******** REARRANGEMENT OF STIFFNESS MATRIX AND FORCE VECTER *******
C
92 DO 115 I=1,NT
115 INDEX(I)=I
DO 120 I=1,KOX
N=2*NOKX(I)-1
120 INDEX(N)=0
DO 125 I=1,KOY
N=2*NOKY(I)
125 INDEX(N)=0
MM=0
DO 130 I=1,NT
IF(INDEX(I)-0) 131,130,131
131 MM=MM+1
INDEX(MM)=I
130 CONTINUE
DO 135 I=1,MM
IA=INDEX(I)
F(I)=F(IA)
DO 135 J=1,MM
JA=INDEX(J)
135 TSM(I,J)=TSM(IA,JA)
WRITE(6,1700)
1700 FORMAT(1H0,1X,48H******* REARRANGEMENT OF STIFFNESS MATRIX ******)
DO 1800 I=1,IA
1800 WRITE(6,1900) (TSM(I,J),J=1,JA)
1900 FORMAT(1H0,5X,8E15.5)
C
C ***** CALCULATION OF DISPLACEMENT *****
C
DO 140 I=1,NT
140 DIS(I)=0.
CALL MATIN(MM)
WRITE(6,2000)
2000 FORMAT(1H0,1X,35H******** INVERSE OF MATRIX ********)
DO 2100 I=1,MM
2100 WRITE(6,2200) (BI(I,J),J=1,MM)
2200 FORMAT(1H0,5X,8E15 .5)
DO 141 I=1,MM
REAC(I)=0.
DO 141 K=1,MM
141 REAC(I)=REAC(I)+TSM(I,K)*F(K)
DO 143 I=1,MM
IA=INDEX(I)
143 DIS(IA)=REAC(I)
C
IF(NANGL) 95,94,95
C
95 DO 600 N=1,NANGL
NI=NONGL(N)
DO 620 I=1,2
WD(I)=0
DO 620 K=1,2
KA=(NI-1)*2+K
IK=(I-1)*2+K
620 WD(I)=WD(I)+W(N,IK)*DIS(KA)
DO 600 I=1,2
IA=(NI-1)*2+I
600 DIS(IA)=WD(I)
C
C ********* CALCULATION OF STRESS AND REACTION *********
94 DO 150 I=1,NT
150 REAC(I)=0.
DO 155 NE=1,NELT
CALL PLANS(NSTRE,NE,E,PO,NODT,NELT,AREA)
DO 160 I=1,3
IA=KAKOM(NE,I)
WD(2*I-1)=DIS(2*IA-1)
160 WD(2*I)=DIS(2*IA)
DO 165 I=1,3
WS(I)=0.
DO 165 K=1,6
165 WS(I)=WS(I)+STR(I,K)*WD(K)
DO 170 I=1,3
170 TSM(NE,I)=WS(I)
WS1=WS(1)
WS2=WS(2)
WS3=WS(3)
CALL PRST(WS1,WS2,WS3,SA,SB,A)
TSM(NE,4)=SA
TSM(NE,5)=SB
TSM(NE,6)=A
DO 180 I=1,6
WS(I)=0.
DO 180 K=1,6
180 WS(I)=WS(I)+SM(I,K)*WD(K)
DO 185 I=1,3
IA=KAKOM(NE,I)
DO 185 J=1,2
IR=2*(IA-1)+J
IW=2*(I-1)+J
185 REAC(IR)=REAC(IR)+WS(IW)
155 CONTINUE
C
IF(NANGL) 97,96,97
C KYOKUSYO ZAHYO E HENKAN SURU
97 DO 814 N=1,NANGL
NI=NONGLE(N)
DO 813 I=1,2
WD(I)=0.
DO 813 K=1,2
KA=(NI-1)*2+K
KI=(K-1)*2+I
813 WD(I)=WD(I)+W(N,KI)*REAC(KA)
DO 814 I=1,2
IA=(NI-1)*2+I
814 REAC(IA)=WD(I)
C
C ******** PRINT OF CALCULATEED RESULTS
C
96 WRITE(6,200)
200 FORMAT(1H0,1X,43H******** DISPLACEMENT AND REACTION ********)
WRITE(6,205)
205 FORMAT(1H0,1X,8HNODE NO.,10X,1HU,14X,1HV,13X,2HRX,13X,2HRY)
WRITE(6,210) (I,DIS(2*I-1),DIS(2*I),REAC(2*I-1),
1 REAC(2*I),I=1,NODT)
210 FORMAT(1H0,I5,5X,4E15.5)
WRITE(6,215)
215 FORMAT(1H0,1X,47H********* STRESS AND PRINCEPAL STRESS *********)
WRITE(6,220)
220 FORMAT(1H0,1X,8HELEM.NO.,7X,2HSX,11X,2HSY,11X,3HSXY,15X,2HS1,11X
1 ,2HS2,10X,5HANGLE)
DO 225 I=1,NELT
225 WRITE(6,230) I,(TSM(I,J),J=1,6)
230 FORMAT(1H0,I5,5X,3E13.5,5X,2E13.5,F10.1)
GO TO 1000
3000 STOP
END
C
C ******** SUBROUTINE OF STIFFNESS MATRIX FOR 2-DIMENSIONAL PROBLEM
C
SUBROUTINE PLANS(NSTRE,NE,E,PO,NODT,NELT,AREA)
C
DIMENSION B(3,6),D(3,3)
COMMON KAKOM(50,3),X(40),Y(40),T(50),SM(6,6),TSM(80,80),STR(3,6)
C
I=KAKOM(NE,1)
J=KAKOM(NE,2)
K=KAKOM(NE,3)
XI=X(I)
XJ=X(J)
XK=X(K)
YI=Y(I)
YJ=Y(J)
YK=Y(K)
AREA=0.5*((XK-XJ)*YI+(XI-XK)*YJ+(XJ-XI)*YK)
IF(AREA-0.) 100,100,5
5 DO 10 I=1,3
DO 10 J=1,6
10 B(I,J)=0.
B(1,1)=YJ-YK
B(1,3)=YK-YI
B(1,5)=YI-YJ
B(2,2)=XK-XJ
B(2,4)=XI-XK
B(2,6)=XJ-XI
B(3,1)=B(2,2)
B(3,2)=B(1,1)
B(3,3)=B(2,4)
B(3,4)=B(1,3)
B(3,5)=B(2,6)
B(3,6)=B(1,5)
DO 15 I=1,3
DO 15 J=1,6
15 B(I,J)=B(I,J)*0.5/AREA
DO 20 I=1,3
DO 20 J=1,3
20 D(I,J)=0.
IF(NSTRE) 25,30,25
25 D(1,1)=E/(1.-PO**2)
D(1,2)=PO*D(1,1)
D(2,1)=D(1,2)
D(2,2)=D(1,1)
D(3,3)=E*0.5/(1.+PO)
GO TO 35
30 D(1,1)=E*(1.-PO)/((1.+PO)*(1.-2.*PO))
D(1,2)=E*PO/((1.+PO)*(1.-2.*PO))
D(2,1)=D(1,2)
D(2,2)=D(1,1)
D(3,3)=E*0.5/(1.+PO)
C
C ******** STRESS MATRIX ********
C
35 DO 40 I=1,3
DO 40 J=1,6
STR(I,J)=0.
DO 40 K=1,3
40 STR(I,J)=STR(I,J)+D(I,K)*B(K,J)
C
C ********* STIFFNESS MATRIX *********
C
DO 45 I=1,6
DO 45 J=1,6
SM(I,J)=0.
DO 45 K=1,3
45 SM(I,J)=SM(I,J)+B(K,I)*STR(K,J)
IF(NSTRE) 46,50,46
46 TH=T(NE)
GO TO 55
50 TH=1.
55 DO 60 I=1,6
DO 60 J=1,6
60 SM(I,J)=SM(I,J)*AREA*TH
GO TO 70
C **********************************
100 WRITE(6,110)
110 FORMAT(1H0,1X,29H********* DATA MISS *********)
WRITE(6,115) NE, I,J,K
115 FORMAT(1H0,5X,12HELEMENT NO.=,I3,6X,9HNODE NO.=,I3,2I10)
WRITE(6,120) XI,XJ,XK
120 FORMAT(1H0,23X,5H--X--,2X,3F10.3)
WRITE(6,125) YI,YJ,YK
125 FORMAT(1H0,23X,5H--Y--,2X,3F10.3)
70 CONTINUE
RETURN
END
C
C
C
SUBROUTINE PRST(SGMAX,SGMAY,TAUXY,PS1,PS2,ANGLE)
SXPY=0.5*(SGMAX+SGMAY)
SXMY=0.5*(SGMAX-SGMAY)
TEMP=SQRT(SXMY**2+TAUXY**2)
PS1=SXPY+TEMP
PS2=SXPY-TEMP
IF(SXTY) 100,101,100
100 TMP=28.64789*ATAN(TAUXY/SXMY)
IF(SXMY) 102,101,103
102 ANGLE=TMP+90.0
GO TO 110
101 IF(TAUXY) 104,105,106
104 ANGLE=135.0
GO TO 110
105 ANGLE=360.0
GO TO 110
106 ANGLE=45.0
GO TO 110
103 IF(TAUXY) 107,108,108
107 ANGLE=TMP+180.0
GO TO 110
108 ANGLE=TMP
110 RETURN
END
C
SUBROUTINE MATIN(NN)
DIMENSION IND(90)
COMMON KAKOM(50,3),X(40),Y(40),T(50),SM(6,6),BI(80,80),STR(3,6)
DO 102 KO=1,NN
102 IND(KO)=KO
DO 103 KO=1,NN
W=0.0
DO 104 IO=KO,NN
IF(ABS(BI(IO,1))-W) 104,104,200
200 W=ABS(BI(IO,1))
IR=IO
104 CONTINUE
IF(IR-KO) 201,106,201
201 DO 107 JO=1,NN
W=BI(KO,JO)
BI(KO,JO)=BI(IR,JO)
107 BI(IR,JO)=W
M=IND(KO)
IND(KO)=IND(IR)
IND(IR)=M
106 W=BI(KO,1)
NNA=NN-1
DO 108 JO=1,NNA
108 BI(KO,JO)=BI(KO,JO+1)/W
BI(KO,NN)=1./W
DO 109 IO=1,NN
IF(IO-KO) 202,109,202
202 W=BI(IO,1)
DO 110 JO=1,NNA
110 BI(IO,JO)=BI(IO,JO+1)-W*BI(KO,JO)
BI(IO,NN)=-W*BI(KO,NN)
109 CONTINUE
103 CONTINUE
DO 111 KO=1,NNA
IF(KO-IND(KO)) 203,111,203
203 KA=KO+1
DO 112 IO=KA,NN
IF(KO-IND(IO)) 112,204,112
204 IR=IO
GO TO 114
112 CONTINUE
114 DO 115 JO=1,NN
W=BI(JO,KO)
BI(JO,KO)=BI(JO,IR)
115 BI(JO,IR)=W
IND(IR)=IND(KO)
IND(KO)=KO
111 CONTINUE
RETURN
END
C
SUBROUTINE ASMAT(NE,NOD,NHEN,NT,NELT,NS)
C
COMMON KAKOM(50,3),X(40),Y(40),T(50),SM(6,6),TSM(80,80),STR(3,6)
C
DO 10 I=1,NOD
DO 10 J=1,NOD
KI=(KAKOM(NE,I)-1)*NHEN
KJ=(KAKOM(NE,J)-1)*NHEN
IS=(I-1)*NHEN
JS=(J-1)*NHEN
DO 10 K=1,NHEN
DO 10 L=1,NHEN
KIK=KI+K
KJL=KJ+L
ISK=IS+K
JSL=JS+L
10 TSM(KIK,KJL)=TSM(KIK,KJL)+SM(ISK,JSL)
RETURN
END