有限元入门
好的入门资料：
1.朱伯芳的《有限单元法原理与应用》
2.刘尔烈 《有限单元法及程序设计-刘尔烈》
3. 单元刚度矩阵组装及整体分析.doc
4.矿大的ppt ：第五章-三角形单元的有限元法
5.下面的例子是刘尔烈书中例子：
读入文件：
1 EXAMPLE 2 6,4,5,3,1 3 1,0,1,0 4 3,1,2,5,2,4,2,5,3,6,3,5 5 0,2,0,1,1,1,0,0,1,0,2,0 6 1,1,0,2,1,0,4,1,1,5,0,1,6,0,1 7 1,-0.5,-1.5,3,-1.,-1.,6,-0.5,-0.5
程序代码：
1 PROGRAM PSSPAP 2 C PLAN STRESS/STRAIN PROBLAM ANALYSIS PROGRAM 3 DIMENSION LND(500,3),X(200),Y(200),JZ(50,3), 4 & PJ(50,3),P(500),TL(20),AK(500,100),AKE(6,6) 5 OPEN (6,FILE='PSSPAPIN.TXT') 6 OPEN (8,FILE='PSSPAPOUT.TXT') 7 READ(6,10) TL 8 10 FORMAT(20A4) 9 READ(6,*)NJ,NE,NZ,NPJ,IPS 10 WRITE(8,10)TL 11 IF(IPS.EQ.1) WRITE(8,20) 12 IF(IPS.EQ.2) WRITE(8,30) 13 20 FORMAT(/1X,'PLATE STRESS PROBLEM') 14 30 FORMAT(/1X,'PLATE STRAIN PROBLEM') 15 NPJ0=NPJ 16 IF(NPJ.EQ.0) NPJ=1 17 CALL READ(NJ,NE,NZ,ND,NPJ,NPJ0,IPS,E,PR, 18 & T,V,LND,X,Y,JZ,PJ) 19 N=2*NJ 20 DO I=1,N 21 DO J=1,ND 22 AK(I,J)=0.0 23 END DO 24 END DO 25 DO IE=1,NE 26 CALL MKE(IE,NJ,NE,LND,X,Y,E,PR,T,AKE) 27 DO I=1,3 28 DO II=1,2 29 IH=2*(I-1)+II 30 IDH=2*(LND(IE,I)-1)+II 31 DO J=1,3 32 DO JJ=1,2 33 L=2*(J-1)+JJ 34 IL=2*(LND(IE,J)-1)+JJ 35 IDL=IL-IDH+1 36 IF(IDL.GT.0) THEN 37 AK(IDH,IDL)=AK(IDH,IDL)+AKE(IH,L) 38 END IF 39 END DO 40 END DO 41 END DO 42 END DO 43 END DO 44 CALL MF(NJ,NE,NPJ,NPJ0,N,T,V,LND,X,Y,PJ,P) 45 CALL RKR(NZ,ND,N,JZ,AK,P) 46 CALL SLOV(NJ,N,ND,AK,P) 47 CALL MADE(NE,NJ,N,E,PR,LND,X,Y,P) 48 CLOSE (6) 49 CLOSE (8) 50 STOP 51 END 52 SUBROUTINE READ(NJ,NE,NZ,ND,NPJ,NPJ0,IPS,E,PR,T,V, 53 & LND,X,Y,JZ,PJ) 54 DIMENSION LND(500,3),X(NJ),Y(NJ),JZ(NZ,3),PJ(NPJ,3) 55 READ(6,*) E,PR,T,V 56 READ(6,*)((LND(I,J),J=1,3),I=1,NE) 57 READ(6,*)(X(I),Y(I),I=1,NJ) 58 READ(6,*)((JZ(I,J),J=1,3),I=1,NZ) 59 IF(NPJ0.NE.0) THEN 60 READ(6,*)((PJ(I,J),J=1,3),I=1,NPJ) 61 END IF 62 ND=0 63 DO IE=1,NE 64 DO I=1,3 65 DO J=1,3 66 IW=IABS(LND(IE,I)-LND(IE,J)) 67 IF(ND.LT.IW) THEN 68 ND=IW 69 END IF 70 END DO 71 END DO 72 END DO 73 ND=(ND+1)*2 74 IF(IPS.NE.1) THEN 75 E=E/(1.0-PR*PR) 76 PR=PR/(1.0-PR) 77 END IF 78 END 79 SUBROUTINE MKE(IE,NJ,NE,LND,X,Y,E,PR,T,AKE) 80 DIMENSION LND(NE,3),X(NJ),Y(NJ),B(3,6),D(3,3),S(3,6), 81 & AKE(6,6) 82 CALL MA(IE,NJ,NE,LND,X,Y,AE) 83 CALL MD(E,PR,D) 84 CALL MB(IE,NJ,NE,LND,X,Y,AE,B) 85 DO I=1,3 86 DO J=1,6 87 S(I,J)=0.0 88 DO K=1,3 89 S(I,J)=S(I,J)+D(I,K)*B(K,J) 90 END DO 91 END DO 92 END DO 93 DO I=1,6 94 DO J=1,6 95 AKE(I,J)=0.0 96 DO K=1,3 97 AKE(I,J)=AKE(I,J)+S(K,I)*B(K,J)*AE*T 98 END DO 99 END DO 100 END DO 101 END 102 SUBROUTINE MA(IE,NJ,NE,LND,X,Y,AE) 103 DIMENSION LND(500,3),X(NJ),Y(NJ) 104 I=LND(IE,1) 105 J=LND(IE,2) 106 K=LND(IE,3) 107 XIJ=X(J)-X(I) 108 YIJ=Y(J)-Y(I) 109 XIK=X(K)-X(I) 110 YIK=Y(K)-Y(I) 111 AE=0.5*(XIJ*YIK-XIK*YIJ) 112 RETURN 113 END 114 SUBROUTINE MD(E,PR,D) 115 DIMENSION D(3,3) 116 DO I=1,3 117 DO J=1,3 118 D(I,j)=0.0 119 END DO 120 END DO 121 D(1,1)= E/(1.0-PR*PR) 122 D(1,2)=E*PR/(1.0-PR*PR) 123 D(2,1)=D(1,2) 124 D(2,2)=D(1,1) 125 D(3,3)=0.5*E/(1.0+PR) 126 RETURN 127 END 128 SUBROUTINE MB(IE,NJ,NE,LND,X,Y,AE,B) 129 DIMENSION LND(500,3),X(NJ),Y(NJ),B(3,6) 130 I=LND(IE,1) 131 J=LND(IE,2) 132 K=LND(IE,3) 133 DO II=1,3 134 DO JJ=1,6 135 B(II,JJ)=0.0 136 END DO 137 END DO 138 B(1,1)=Y(J)-Y(K) 139 B(1,3)=Y(K)-Y(I) 140 B(1,5)=Y(I)-Y(J) 141 B(2,2)=X(K)-X(J) 142 B(2,4)=X(I)-X(K) 143 B(2,6)=X(J)-X(I) 144 B(3,1)=B(2,2) 145 B(3,2)=B(1,1) 146 B(3,3)=B(2,4) 147 B(3,4)=B(1,3) 148 B(3,5)=B(2,6) 149 B(3,6)=B(1,5) 150 DO I1=1,3 151 DO J1=1,6 152 B(I1,J1)=0.5/AE*B(I1,J1) 153 END DO 154 END DO 155 RETURN 156 END 157 SUBROUTINE MF(NJ,NE,NPJ,NPJ0,N,T,V,LND,X,Y,PJ,P) 158 DIMENSION LND(500,3),X(NJ),Y(NJ),PJ(NPJ,3),P(N) 159 DO I=1,N 160 P(I)=0.0 161 END DO 162 IF(NPJ0.NE.0) THEN 163 DO I=1,NPJ 164 II=PJ(I,1) 165 P(2*II-1)=PJ(I,2) 166 P(2*II)=PJ(I,3) 167 END DO 168 END IF 169 IF (V.NE.0) THEN 170 DO IE=1,NE 171 CALL MA(IE,NJ,NE,LND,X,Y,AE) 172 PE=-V*AE*T/3.0 173 DO I=1,3 174 II=LND(IE,I) 175 P(2*II)=P(2*II)+PE 176 END DO 177 END DO 178 END IF 179 RETURN 180 END 181 SUBROUTINE RKR(NZ,ND,N,JZ,AK,P) 182 DIMENSION P(N),JZ(NZ,3),AK(500,100) 183 DO I=1,NZ 184 IR=JZ(I,1) 185 DO J=2,3 186 IF (JZ(I,J).NE.0) THEN 187 II=2*IR+J-3 188 AK(II,1)=1.0 189 DO JJ=2,ND 190 AK(II,JJ)=0.0 191 END DO 192 IF (II.GT.ND) JO=ND 193 IF (II.LE.ND) JO=II 194 DO JJ=2,JO 195 AK(II-JJ+1,JJ)=0.0 196 END DO 197 P(II)=0.0 198 END IF 199 END DO 200 END DO 201 RETURN 202 END 203 SUBROUTINE SLOV(NJ,N,ND,AK,P) 204 DIMENSION P(N),AK(500,100) 205 NJ1=N-1 206 DO K=1,NJ1 207 IF (N.GT.K+ND-1) IM=K+ND-1 208 IF (N.LE.K+ND-1) IM=N 209 K1=K+1 210 DO I=K1,IM 211 L=I-K+1 212 C=AK(K,L)/AK(K,1) 213 IW=ND-L+1 214 DO J=1,IW 215 M=J+I-K 216 AK(I,J)=AK(I,J)-C*AK(K,M) 217 END DO 218 P(I)=P(I)-C*P(K) 219 END DO 220 END DO 221 P(N)=P(N)/AK(N,1) 222 DO I1=1,NJ1 223 I=N-I1 224 IF(ND.GT.N-I+1) JO=N-I+1 225 IF(ND.LT.N-I+1) JO=ND 226 DO J=2,JO 227 K=J+I-1 228 P(I)=P(I)-AK(I,J)*P(K) 229 END DO 230 P(I)=P(I)/AK(I,1) 231 END DO 232 WRITE(8,50) 233 50 FORMAT (/8X,5('**'),' RESULTE OF CALCULATION ', 234 & 5('**')//1X,'NODAL DISPLACEMENTS' //3X,'NODE', 235 & 5X,'X-DISP.',8X,'Y-DISP.') 236 DO I=1,NJ 237 WRITE(8,70) I,P(2*I-1),P(2*I) 238 70 FORMAT(1X,I5,2E15.6) 239 END DO 240 RETURN 241 END 242 SUBROUTINE MADE(NE,NJ,N,E,PR,LND,X,Y,P) 243 DIMENSION LND(500,3),X(NJ),Y(NJ),D(3,3), 244 & B(3,6),S(3,6),ST(3),P(N),DE(6) 245 WRITE(8,10) 246 10 FORMAT(/1X,'ELEMENT STRESSES'/) 247 CALL MD(E,PR,D) 248 DO IE=1,NE 249 CALL MA(IE,NJ,NE,LND,X,Y,AE) 250 CALL MB(IE,NJ,NE,LND,X,Y,AE,B) 251 DO I=1,3 252 DO J=1,6 253 S(I,J)=0.0 254 DO K=1,3 255 S(I,J)=S(I,J)+D(I,K)*B(K,J) 256 END DO 257 END DO 258 END DO 259 DO I=1,3 260 DO J=1,2 261 IH=2*(I-1)+J 262 IW=2*(LND(IE,I)-1)+J 263 DE(IH)=P(IW) 264 END DO 265 END DO 266 DO I=1,3 267 ST(I)=0.0 268 DO J=1,6 269 ST(I)=ST(I)+S(I,J)*DE(J) 270 END DO 271 END DO 272 STX=ST(1) 273 STY=ST(2) 274 TXY=ST(3) 275 AST=(STX+STY)*.5 276 RST=SQRT(0.25*(STX-STY)**2+TXY*TXY) 277 STMA=AST+RST 278 STMI=AST-RST 279 IF (STY.EQ.STMI) CETA=0.0 280 IF (STY.NE.STMI) CETA=90.-57.29578*ATAN 281 & (TXY/(STY-STMI)) 282 WRITE(8,60) IE,STX,STY,TXY,STMA,STMI,CETA 283 60 FORMAT (1X,'ELEMENT NO.=',I5/3X,' STX=',E15.6, 284 & 2X,' STY=',E15.6,2X,' TXY=',E15.6/3X,'STMA=', 285 & E15.6,2X,'STMI=',E15.6,2X,'CETA=',E15.6) 286 END DO 287 RETURN 288 END
输出结果：
1 EXAMPLE 2 3 PLATE STRESS PROBLEM 4 5 ********** RESULTE OF CALCULATION ********** 6 7 NODAL DISPLACEMENTS 8 9 NODE X-DISP. Y-DISP. 10 1 0.000000E+00 -0.525275E+01 11 2 0.000000E+00 -0.225275E+01 12 3 -0.108791E+01 -0.137363E+01 13 4 0.000000E+00 0.000000E+00 14 5 -0.824176E+00 0.000000E+00 15 6 -0.182418E+01 0.000000E+00 16 17 ELEMENT STRESSES 18 19 ELEMENT NO.= 1 20 STX= -0.108791E+01 STY= -0.300000E+01 TXY= 0.439560E+00 21 STMA= -0.991704E+00 STMI= -0.309621E+01 CETA= 0.123458E+02 22 ELEMENT NO.= 2 23 STX= -0.824176E+00 STY= -0.225275E+01 TXY= 0.000000E+00 24 STMA= -0.824176E+00 STMI= -0.225275E+01 CETA= 0.000000E+00 25 ELEMENT NO.= 3 26 STX= -0.108791E+01 STY= -0.137363E+01 TXY= 0.307692E+00 27 STMA= -0.891531E+00 STMI= -0.157001E+01 CETA= 0.325476E+02 28 ELEMENT NO.= 4 29 STX= -0.100000E+01 STY= -0.137363E+01 TXY= -0.131868E+00 30 STMA= -0.958147E+00 STMI= -0.141548E+01 CETA= 0.162391E+03
