1 C 这是一个采用平面四边形八节点等参单元的平面有限元分析程序。 2 PROGRAM PLANEFEM 3 IMPLICIT REAL*8(A-H,O-Z) 4 CHARACTER*80 LINECHAR,NEWLINECHAR 5 COMMON A(30000),L(4000) 6 COMMON /SOL/NPOIN,NELEM,NTYPE,NMATS 7 OPEN(5,FILE='FEMDATA',STATUS='OLD') 8 OPEN(6,FILE='FEMOUT',STATUS='UNKNOWN') 9 C 以下进行的是从数据文件FEMDATA中读入数据文件主标题信息。 10 READ(5,5000) LINECHAR 11 LOCATECHAR=INDEX(LINECHAR,'输入') 12 IF(LOCATECHAR.NE.0) THEN 13 LINECHAR(LOCATECHAR:LOCATECHAR+3)='输出' 14 ENDIF 15 WRITE(6,5100) LINECHAR 16 5000 FORMAT(A) 17 5100 FORMAT(80('*')/A/80('*')/) 18 C 以下进行的是先读入一行字符，然后从字符中读入网格单元总数NELEM。 19 READ(5,5000) LINECHAR 20 WRITE(6,5000) LINECHAR 21 LOCATECHAR=INDEX(LINECHAR,'=') 22 LOCATECHAR=LOCATECHAR+1 23 NEWLINECHAR=LINECHAR(LOCATECHAR:80) 24 READ(NEWLINECHAR,5200) NELEM 25 5200 FORMAT(I5) 26 C 以下进行的是先读入一行字符，然后从字符中读入单元节点总数NPOIN。 27 READ(5,5000) LINECHAR 28 WRITE(6,5000) LINECHAR 29 LOCATECHAR=INDEX(LINECHAR,'=') 30 LOCATECHAR=LOCATECHAR+1 31 NEWLINECHAR=LINECHAR(LOCATECHAR:80) 32 READ(NEWLINECHAR,5200) NPOIN 33 C 以下进行的是先读入一行字符，然后从字符中读入问题类型编号NTYPE。 34 READ(5,5000) LINECHAR 35 WRITE(6,5000) LINECHAR 36 LOCATECHAR=INDEX(LINECHAR,'=') 37 LOCATECHAR=LOCATECHAR+1 38 NEWLINECHAR=LINECHAR(LOCATECHAR:80) 39 READ(NEWLINECHAR,5200) NTYPE 40 C 以下进行的是先读入一行字符，然后从字符中读入材料类型总数NMATS。 41 READ(5,5000) LINECHAR 42 WRITE(6,5000) LINECHAR 43 LOCATECHAR=INDEX(LINECHAR,'=') 44 LOCATECHAR=LOCATECHAR+1 45 NEWLINECHAR=LINECHAR(LOCATECHAR:80) 46 READ(NEWLINECHAR,5200) NMATS 47 C 以下进行的是根据以输入的控制参数确定虚拟数组在实数组中的起始位置。 48 N1=1 49 N2=N1+NPOIN*2 50 N3=N2+NMATS*3 51 N4=N3+NPOIN*2 52 NN2=N1+NELEM*8 53 NN3=NN2+NPOIN 54 NN4=NN3+NELEM 55 NN5=NN4+NELEM 56 C 以下进行的是调用子程序进行有限元分析计算。 57 CALL INPUT(A(N1),A(N2),L(N1),L(NN2),L(NN3)) 58 CALL STIF(A(N1),A(N2),L(N1),L(NN3)) 59 CALL LOAD(A(N1),L(N1),L(NN4)) 60 CALL ASEM(A(N3),A(N4),L(N1),L(NN2),L(NN4),L(NN5)) 61 CALL SOLVE(A(N3),A(N4),L(NN5)) 62 CALL STRE(A(N1),A(N2),A(N3),L(N1),L(NN3)) 63 WRITE(*,5300) 64 5300 FORMAT(////////////////) 65 STOP ' *** 程序正常运行结束 ***' 66 END 67 C***************************************************************************** 68 C 这是一个进行数据输入并进行数据检验的子程序。 * 69 C***************************************************************************** 70 SUBROUTINE INPUT(COORD,PROPS,LNODS,IFPRE,MATNO) 71 IMPLICIT REAL*8(A-H,O-Z) 72 CHARACTER*80 LINECHAR,NEWLINECHAR,CHARCOMP 73 LOGICAL NOTREADCHAR 74 DIMENSION COORD(2,1),PROPS(3,1),LNODS(8,1),IFPRE(1),MATNO(1) 75 DIMENSION NEROR(6) 76 COMMON/SOL/NPOIN,NELEM,NTYPE,NMATS 77 CHARCOMP='aa(这个字符常量在ASCII序列中介于数字和字符之间)' 78 NNODE=8 79 C 以下进行的是对材料类型编号赋初值。 80 I=1 81 DO WHILE (I.LE.NELEM) 82 MATNO(I)=1 83 I=I+1 84 ENDDO 85 C 以下进行的是输入材料类型编号信息标题。 86 READ(5,5000) LINECHAR 87 WRITE(6,5000) LINECHAR 88 5000 FORMAT(A) 89 C 以下进行的是输入材料类型编号。 90 DO WHILE (NMATS.GT.1) 91 NUMBERREAD=1 92 DO WHILE (NUMBERREAD.LE.NMATS) 93 READ(5,*) I,MATNO(I) 94 NUMBERREAD=NUMBERREAD+1 95 ENDDO 96 ENDDO 97 C 以下进行的是输出材料类型编号。 98 I=1 99 DO WHILE (I.LE.NELEM) 100 WRITE(6,6000) I,MATNO(I) 101 I=I+1 102 ENDDO 103 6000 FORMAT(2X,'单元编号为：',I4,10X,'材料类型编号为：'I6) 104 C 以下进行的是输入单元节点总体编号标题。 105 READ(5,5000) LINECHAR 106 LOCATECHAR=INDEX(LINECHAR,'输入') 107 IF(LOCATECHAR.NE.0) THEN 108 LINECHAR(LOCATECHAR:LOCATECHAR+3)='输出' 109 ENDIF 110 WRITE(6,5000) LINECHAR 111 C 以下进行的是输入单元节点总体编号。 112 NOTREADCHAR=.TRUE. 113 DO WHILE (NOTREADCHAR) 114 READ(5,5000) LINECHAR 115 IF (LINECHAR(1:2).LT.CHARCOMP(1:2)) THEN 116 READ(LINECHAR,6100) IELEM,(LNODS(I,IELEM),I=1,NNODE) 117 J=1 118 DO WHILE (J.LE.10) 119 KELEM=IELEM+J 120 IF (KELEM.GT.NELEM) GOTO 100 121 K=1 122 DO WHILE (K.LE.8) 123 LNODS(K,KELEM)=LNODS(K,IELEM)+J 124 K=K+1 125 ENDDO 126 100 J=J+1 127 ENDDO 128 ELSE 129 NOTREADCHAR=.FALSE. 130 END IF 131 ENDDO 132 6100 FORMAT(9I5) 133 C 以下进行的是输出单元节点总体编号。 134 IELEM=1 135 DO WHILE (IELEM.LE.NELEM) 136 WRITE(6,6200) IELEM,(LNODS(I,IELEM),I=1,NNODE) 137 IELEM=IELEM+1 138 ENDDO 139 6200 FORMAT(2X,'单元编号为：',I4,' 节点总体编号为：',8I5) 140 C 以下进行的是输出单元节点的直角坐标信息标题。 141 LOCATECHAR=INDEX(LINECHAR,'输入') 142 IF(LOCATECHAR.NE.0) THEN 143 LINECHAR(LOCATECHAR:LOCATECHAR+3)='输出' 144 ENDIF 145 WRITE(6,5000) LINECHAR 146 C 以下进行的是输入单元节点的直角坐标信息。 147 NOTREADCHAR=.TRUE. 148 DO WHILE (NOTREADCHAR) 149 READ(5,5000) LINECHAR 150 IF (LINECHAR(1:2).LT.CHARCOMP(1:2)) THEN 151 READ(LINECHAR,6300) IPOIN,(COORD(J,IPOIN),J=1,2) 152 ELSE 153 NOTREADCHAR=.FALSE. 154 END IF 155 ENDDO 156 6300 FORMAT(I5,2F15.5) 157 C 以下进行的是输出单元节点的极坐标信息标题。 158 LOCATECHAR=INDEX(LINECHAR,'输入') 159 IF(LOCATECHAR.NE.0) THEN 160 LINECHAR(LOCATECHAR:LOCATECHAR+3)='输出' 161 ENDIF 162 WRITE(6,5000) LINECHAR 163 C 以下进行的是输入单元节点的极坐标信息。 164 PI=DATAN(1.D0)*4.D0 165 NOTREADCHAR=.TRUE. 166 DO WHILE (NOTREADCHAR) 167 READ(5,5000) LINECHAR 168 IF (LINECHAR(1:2).LT.CHARCOMP(1:2)) THEN 169 READ(LINECHAR,6300) IPOIN,R,ANGLE 170 ANGLE=ANGLE/180.D0*PI 171 COORD(1,IPOIN)=R*DCOS(ANGLE) 172 COORD(2,IPOIN)=R*DSIN(ANGLE) 173 NUMBERREAD=NUMBERREAD+1 174 ELSE 175 NOTREADCHAR=.FALSE. 176 END IF 177 ENDDO 178 C 以下进行的是根据输入信息计算出单元节点的坐标。 179 I=1 180 DO WHILE (I.LE.NELEM) 181 INODE=1 182 DO WHILE (INODE.LE.NNODE) 183 NODST =LNODS(INODE,I) 184 IGASH=INODE+2 185 IF(IGASH.GT.NNODE) IGASH=1 186 NDOFN=LNODS(IGASH,I) 187 MIDPT=INODE+1 188 NODMD=LNODS(MIDPT,I) 189 J=1 190 DO WHILE (J.LE.2) 191 IF(DABS(COORD(J,NODMD)).LT.1.E-6) THEN 192 COORD(J,NODMD)=(COORD(J,NDOFN)+COORD(J,NODST))/2.0 193 ENDIF 194 J=J+1 195 ENDDO 196 INODE=INODE+2 197 ENDDO 198 I=I+1 199 ENDDO 200 C 以下进行的是输出单元节点的坐标信息。 201 WRITE(6,6400) 202 I=1 203 DO WHILE (I.LE.NPOIN) 204 WRITE(6,6500) I,(COORD(J,I),J=1,2) 205 I=I+1 206 ENDDO 207 6400 FORMAT(4X,'节点编号',10X,'Ｘ坐标',16X,'Ｙ坐标') 208 6500 FORMAT(2X,I10,F16.5,6X,F16.5) 209 C 以下进行的是输出单元节点约束代码信息标题。 210 LOCATECHAR=INDEX(LINECHAR,'输入') 211 IF(LOCATECHAR.NE.0) THEN 212 LINECHAR(LOCATECHAR:LOCATECHAR+3)='输出' 213 ENDIF 214 WRITE(6,5000) LINECHAR 215 C 以下进行的是输入单元节点约束代码信息。 216 NOTREADCHAR=.TRUE. 217 DO WHILE (NOTREADCHAR) 218 READ(5,5000) LINECHAR 219 IF (LINECHAR(1:2).LT.CHARCOMP(1:2)) THEN 220 READ(LINECHAR,6600) I,IFPRE(I) 221 ELSE 222 NOTREADCHAR=.FALSE. 223 END IF 224 ENDDO 225 6600 FORMAT(2I5) 226 C 以下进行的是输出单元节点约束代码信息。 227 LINECHAR='节点约束代码信息为一个两位数的阿拉伯数字，其中，' 228 WRITE(6,5000) LINECHAR 229 LINECHAR='十位数表示Ｘ方向的约束状况，' 230 NEWLINECHAR='个位数表示Ｙ方向的约束状况，' 231 LOCATECHAR=LEN_TRIM(LINECHAR) 232 LOCATECHAR=LOCATECHAR+1 233 LINECHAR(LOCATECHAR:)=NEWLINECHAR(1:LOCATECHAR) 234 WRITE(6,5000) LINECHAR 235 LINECHAR='具体的，0表示自由，1表示固定，2表示给定非零位移。' 236 WRITE(6,5000) LINECHAR 237 I=1 238 DO WHILE (I.LE.NPOIN) 239 WRITE(6,6700) I,IFPRE(I) 240 I=I+1 241 ENDDO 242 6700 FORMAT(2X,'节点编号为：',I4,13X,'约束代码为：',I4) 243 C 以下进行的是输入材料的特性参数信息标题。 244 LINECHAR='以下输出的是材料特性参数（格式为：' 245 NEWLINECHAR='材料类型编号NMATS,弹性模量Ｅ,泊松比μ,厚度t）' 246 LOCATECHAR=LEN_TRIM(LINECHAR) 247 LOCATECHAR=LOCATECHAR+1 248 LINECHAR(LOCATECHAR:)=NEWLINECHAR(1:46) 249 WRITE(6,5000) LINECHAR 250 C 以下进行的是输入材料的特性参数。 251 NUMBERREAD=1 252 DO WHILE (NUMBERREAD.LE.NMATS) 253 READ(5,*) I,(PROPS(J,I),J=1,3) 254 WRITE(6,6800) I 255 WRITE(6,6900) (PROPS(J,I),J=1,3) 256 NUMBERREAD=NUMBERREAD+1 257 ENDDO 258 6800 FORMAT(2X,'材料类型编号为：',I2,' 的特性参数为：') 259 6900 FORMAT(2X,1PE13.3,1PE13.3,1PE13.3) 260 C 以下进行的是对输入和计算出的数据进行检验。 261 IERROR=1 262 DO WHILE (IERROR.LE.6) 263 NEROR(IERROR)=0 264 IERROR=IERROR+1 265 ENDDO 266 C 以下进行的是对输入和计算出的单元节点坐标进行检验。 267 IPOIN=1 268 DO WHILE (IPOIN.LE.NPOIN) 269 KPOIN=IPOIN-1 270 JPOIN=1 271 DO WHILE (JPOIN.LE.KPOIN) 272 J=1 273 DO WHILE (J.LE.2) 274 IF(COORD(J,IPOIN).NE.COORD(J,JPOIN)) GOTO 200 275 J=J+1 276 ENDDO 277 NEROR(1)=NEROR(1)+1 278 WRITE(6,7000) IPOIN,JPOIN 279 JPOIN=JPOIN+1 280 ENDDO 281 200 CONTINUE 282 IPOIN=IPOIN+1 283 ENDDO 284 7000 FORMAT(2X,'坐标错误！节点编号为',I4,' 的坐标等于',I4,' 的坐标') 285 C 以下进行的是对输入和计算出的材料类型编号进行检验。 286 I=1 287 DO WHILE (I.LE.NELEM) 288 IF(MATNO(I).LE.0.OR.MATNO(I).GT.NMATS) THEN 289 NEROR(2)=NEROR(2)+1 290 END IF 291 I=I+1 292 ENDDO 293 IF(NEROR(2).NE.0) THEN 294 WRITE(6,7100) NEROR(2) 295 END IF 296 7100 FORMAT(2X,' 材料类型编号错误,错误个数为',I2,' 个') 297 C 以下进行的是对输入和计算出的单元节点总体编号进行检验。 298 I=1 299 DO WHILE (I.LE.NELEM) 300 J=1 301 DO WHILE (J.LE.8) 302 IF(LNODS(J,I).LE.0.OR.LNODS(J,I).GT.NPOIN) THEN 303 NEROR(3)=NEROR(3)+1 304 END IF 305 J=J+1 306 ENDDO 307 I=I+1 308 ENDDO 309 IF(NEROR(3).NE.0) THEN 310 WRITE(6,7200) NEROR(3) 311 END IF 312 7200 FORMAT(2X,'单元节点总体编号错误,错误个数为',I5,' 个') 313 C 以下进行的是对输入和计算出的单元节点编号进行检验。 314 I=1 315 DO WHILE (I.LE.NPOIN) 316 KSTAR=0 317 IELEM=1 318 DO WHILE (IELEM.LE.NELEM) 319 KZERO=0 320 J=1 321 DO WHILE(J.LE.8) 322 IF(LNODS(J,IELEM).EQ.I) THEN 323 KZERO=KZERO+1 324 KSTAR=IELEM 325 END IF 326 J=J+1 327 ENDDO 328 IF(KZERO.GT.1) THEN 329 NEROR(4)=NEROR(4)+1 330 WRITE(6,7300) I,IELEM 331 END IF 332 IELEM=IELEM+1 333 ENDDO 334 IF(KSTAR.EQ.0)THEN 335 NEROR(5)=NEROR(5)+1 336 WRITE(6,7400) I 337 END IF 338 I=I+1 339 ENDDO 340 7300 FORMAT(2X,'节点编号',I5,' 在单元',I5,' 中重复。') 341 7400 FORMAT(2X,'总体节点编号',I5,' 在单元中从未出现') 342 C 以下进行的是对输入的单元节点约束代码进行检验。 343 I=1 344 DO WHILE (I.LE.NPOIN) 345 IF(IFPRE(I).LT.0) THEN 346 NEROR(6)=NEROR(6)+1 347 END IF 348 IF(IFPRE(I).GT.22) THEN 349 NEROR(6)=NEROR(6)+1 350 END IF 351 I=I+1 352 ENDDO 353 IF(NEROR(6).NE.0) THEN 354 WRITE(6,7500) NEROR(6) 355 END IF 356 7500 FORMAT(2X,'总体节点编号为',I5,' 的节点约束代码出现错误。') 357 C 以下进行的是判断出错类型。 358 KEROR=0 359 IERROR=1 360 DO WHILE (IERROR.LE.6) 361 IF(NEROR(IERROR).NE.0) THEN 362 WRITE(6,7600) IERROR 363 KEROR=KEROR+1 364 END IF 365 IERROR=IERROR+1 366 ENDDO 367 7600 FORMAT(2X,'出现第 ',I1,' 类错误。') 368 C 以下进行的是判断是否出错,如出错,则中断程序运行,否则返回调用主程序。 369 7700 FORMAT(2X,'*****在输入数据文件中发现错误。程序运行中断*****') 370 IF(KEROR.NE.0) THEN 371 WRITE(6,7700) 372 STOP '*** 程序运行被中断,这是由于输入数据出现错误引起的。***' 373 END IF 374 RETURN 375 END 376 C***************************************************************************** 377 C 这是一个形成单元刚度矩阵的子程序。 * 378 C***************************************************************************** 379 SUBROUTINE STIF(COORD,PROPS,LNODS,MATNO) 380 IMPLICIT REAL*8(A-H,O-Z) 381 DIMENSION COORD(2,1),PROPS(3,1),LNODS(8,1),MATNO(1) 382 DIMENSION POSGP(3),WEIGP(3),ESTIF(136),ELCOD(2,8),DBMAT(3) 383 COMMON/SOL/NPOIN,NELEM,NTYPE,NMATS 384 COMMON/GCOD/GPCOD(2,9) 385 COMMON/BMDMP/BMATX(3,16),DMATX(3,3) 386 OPEN(7,FILE='ESTIF',FORM='UNFORMATTED') 387 C 以下进行的是高斯积分点位置赋值。 388 POSGP(1)=-DSQRT(0.6D0) 389 POSGP(2)=0.0 390 POSGP(3)=DSQRT(0.6D0) 391 C 以下进行的是加权系数赋值。 392 WEIGP(1)=5.0/9.0 393 WEIGP(2)=8.0/9.0 394 WEIGP(3)=5.0/9.0 395 C 以下进行的是应力分量个数赋值。 396 NSTRE=3 397 C 以下进行的是形成弹性矩阵[D]。 398 DO 90 IELEM=1,NELEM 399 WRITE(6,615) IELEM 400 615 FORMAT(2X,'网格单元编号为：',I5) 401 LPROP=MATNO(IELEM) 402 DO 10 I=1,8 403 LNODE=LNODS(I,IELEM) 404 DO 10 J=1,2 405 10 ELCOD(J,I)=COORD(J,LNODE) 406 YOUNG=PROPS(1,LPROP) 407 POISS=PROPS(2,LPROP) 408 THICK=PROPS(3,LPROP) 409 CALL MODPS(YOUNG,POISS) 410 C 以下进行的是存放[K]的一维数组赋初值零。 411 KEVAB=0 412 DO 20 I=1,16 413 DO 20 J=1,I 414 KEVAB=KEVAB+1 415 ESTIF(KEVAB)=0.0 416 20 CONTINUE 417 C 以下进行的是计算积分点处形函数及其导数之值以及对整体坐标的偏导数之值。 418 KGASP=0 419 DO 80 IGAUS=1,3 420 DO 80 JGAUS=1,3 421 KGASP=KGASP+1 422 EXISP=POSGP(IGAUS) 423 ETASP=POSGP(JGAUS) 424 WRITE(6,705) KGASP 425 705 FORMAT(6X,'高斯积分点编号为：',I5) 426 CALL SFR2(EXISP,ETASP) 427 CALL JACOB2(IELEM,DJACB,KGASP,ELCOD) 428 DVOLU=DJACB*WEIGP(IGAUS)*WEIGP(JGAUS)*THICK 429 CALL BMATPS 430 C 以下进行的是弹性矩阵[D]赋初值零。 431 KEVAB=0 432 DO 70 IEVAB=1,16 433 DO 40 ISTRE=1,NSTRE 434 DBMAT(ISTRE)=0.0 435 DO 30 JSTRE=1,NSTRE 436 DBMAT(ISTRE)=DBMAT(ISTRE)+BMATX(JSTRE,IEVAB)*DMATX(JSTRE,ISTRE) 437 30 CONTINUE 438 40 CONTINUE 439 DO 60 JEVAB=1,IEVAB 440 KEVAB=KEVAB+1 441 BTDBM=0.0 442 DO 50 ISTRE=1,NSTRE 443 BTDBM=BTDBM+DBMAT(ISTRE)*BMATX(ISTRE,JEVAB) 444 50 CONTINUE 445 ESTIF(KEVAB)=ESTIF(KEVAB)+BTDBM*DVOLU 446 60 CONTINUE 447 70 CONTINUE 448 80 CONTINUE 449 WRITE(7) ESTIF 450 90 CONTINUE 451 REWIND 7 452 RETURN 453 END 454 C***************************************************************************** 455 C 这是一个形成弹性矩阵的子程序。 * 456 C***************************************************************************** 457 SUBROUTINE MODPS(Y,P) 458 IMPLICIT REAL*8(A-H,O-Z) 459 COMMON/SOL/NPOIN,NELEM,NTYPE,NMATS 460 COMMON/BMDMP/BMATX(3,16),DMATX(3,3) 461 DO 100 I=1,3 462 DO 100 J=1,3 463 DMATX(I,J)=0.0 464 100 CONTINUE 465 C 如果NTYPE=1,表示为平面应力问题。 466 IF(NTYPE.EQ.1) THEN 467 CONST=Y/(1.0-P*P) 468 DMATX(1,1)=CONST 469 DMATX(2,2)=CONST 470 DMATX(1,2)=CONST*P 471 DMATX(2,1)=CONST*P 472 DMATX(3,3)=CONST*(1.0-P)/2.0 473 ENDIF 474 C 如果NTYPE=2,表示为平面应变问题。 475 IF(NTYPE.EQ.2) THEN 476 CONST=Y*(1.0-P)/((1.0+P)*(1.0-2.0*P)) 477 DMATX(1,1)=CONST 478 DMATX(2,2)=CONST 479 DMATX(1,2)=CONST*P/(1.0-P) 480 DMATX(2,1)=CONST*P/(1.0-P) 481 DMATX(3,3)=Y/(2.0*(1.0+P)) 482 ENDIF 483 RETURN 484 END 485 C***************************************************************************** 486 C 这是一个计算形函数当前积分点及形函数对局部坐标的导数值的子程序。 * 487 C***************************************************************************** 488 SUBROUTINE SFR2(S,T) 489 IMPLICIT REAL*8(A-H,O-Z) 490 COMMON/SD/SHAPE(8),DERIV(2,8) 491 C 以下进行的是为简化表达式定义一些变量。 492 S2=S*2.0 493 T2=T*2.0 494 SS=S*S 495 TT=T*T 496 ST=S*T 497 STT=S*T*T 498 SST=S*S*T 499 ST2=S*T*2.0 500 C 以下进行的是给形函数赋值。 501 SHAPE(1)=(-1.0+ST+SS+TT-SST-STT)/4.0 502 SHAPE(2)=(1.0-T-SS+SST)/2.0 503 SHAPE(3)=(-1.0-ST+SS+TT-SST+STT)/4.0 504 SHAPE(4)=(1.0+S-TT-STT)/2.0 505 SHAPE(5)=(-1.0+ST+SS+TT+SST+STT)/4.0 506 SHAPE(6)=(1.0+T-SS-SST)/2.0 507 SHAPE(7)=(-1.0-ST+SS+TT+SST-STT)/4.0 508 SHAPE(8)=(1.0-S-TT+STT)/2.0 509 C 以下进行的是给形函数对局部坐标的导数赋值。 510 DERIV(1,1)=(T+S2-ST2-TT)/4.0 511 DERIV(1,2)=-S+ST 512 DERIV(1,3)=(-T+S2-ST2+TT)/4.0 513 DERIV(1,4)=(1.0-TT)/2.0 514 DERIV(1,5)=(T+S2+ST2+TT)/4.0 515 DERIV(1,6)=-S-ST 516 DERIV(1,7)=(-T+S2+ST2-TT)/4.0 517 DERIV(1,8)=(-1.0+TT)/2.0 518 DERIV(2,1)=(S+T2-SS-ST2)/4.0 519 DERIV(2,2)=(-1.0+SS)/2.0 520 DERIV(2,3)=(-S+T2-SS+ST2)/4.0 521 DERIV(2,4)=-T-ST 522 DERIV(2,5)=(S+T2+SS+ST2)/4.0 523 DERIV(2,6)=(1.0-SS)/2.0 524 DERIV(2,7)=(-S+T2+SS-ST2)/4.0 525 DERIV(2,8)=-T+ST 526 RETURN 527 END 528 C***************************************************************************** 529 C 这是一个形成雅可比矩阵的子程序。 * 530 C***************************************************************************** 531 SUBROUTINE JACOB2(IELEM,DJACB,KGASP,ELCOD) 532 IMPLICIT REAL*8(A-H,O-Z) 533 DIMENSION XJACM(2,2),XJACI(2,2),ELCOD(2,8) 534 COMMON/SD/SHAPE(8),DERIV(2,8) 535 COMMON/CAR/CARTD(2,8) 536 COMMON/GCOD/GPCOD(2,9) 537 C 以下进行的是计算当前积分点的整体坐标。 538 DO 100 I=1,2 539 GPCOD(I,KGASP)=0.0 540 DO 100 J=1,8 541 GPCOD(I,KGASP)=GPCOD(I,KGASP)+ELCOD(I,J)*SHAPE(J) 542 100 CONTINUE 543 C 以下进行的是形成雅可比矩阵[J]的各元素。 544 DO 200 I=1,2 545 DO 200 J=1,2 546 XJACM(I,J)=0.0 547 DO 200 K=1,8 548 XJACM(I,J)=XJACM(I,J)+DERIV(I,K)*ELCOD(J,K) 549 200 CONTINUE 550 C 以下进行的是计算雅可比行列式┃Ｊ┃的值。 551 DJACB=XJACM(1,1)*XJACM(2,2)-XJACM(1,2)*XJACM(2,1) 552 WRITE(6,6000) DJACB 553 6000 FORMAT(14X,'雅可比行列式┃Ｊ┃的值为：',F12.5) 554 IF(DJACB.LT.1.E-6)THEN 555 WRITE(6,6100) IELEM 556 6100 FORMAT('单元编号为：',I4,' 的雅可比行列式的值小于或等于零！') 557 STOP '****** 程序运行被中断于子程序JACOB2 ******' 558 END IF 559 C 以下进行的是形成雅可比矩阵[J]的逆矩阵的各元素。 560 XJACI(1,1)=XJACM(2,2)/DJACB 561 XJACI(2,2)=XJACM(1,1)/DJACB 562 XJACI(1,2)=-XJACM(1,2)/DJACB 563 XJACI(2,1)=-XJACM(2,1)/DJACB 564 C 以下进行的是计算形函数对整体坐标的导数。 565 DO 300 I=1,2 566 DO 300 K=1,8 567 CARTD(I,K)=0.0 568 DO 300 J=1,2 569 CARTD(I,K)=CARTD(I,K)+XJACI(I,J)*DERIV(J,K) 570 300 CONTINUE 571 RETURN 572 END 573 C***************************************************************************** 574 C 这是一个形成应变矩阵的子程序。 * 575 C***************************************************************************** 576 SUBROUTINE BMATPS 577 IMPLICIT REAL*8(A-H,O-Z) 578 COMMON/CAR/CARTD(2,8) 579 COMMON/SD/SHAPE(8),DERIV(2,8) 580 COMMON/BMDMP/BMATX(3,16),DMATX(3,3) 581 COMMON/SOL/NPOIN,NELEM,NTYPE,NMATS 582 COMMON/GCOD/GPCOD(2,9) 583 C 以下进行的是应变矩阵[B]赋初值零。 584 DO 100 I=1,16 585 DO 100 J=1,3 586 100 BMATX(J,I)=0.0 587 C 以下进行的是形成应变矩阵[B]。 588 NGASH=0 589 DO 200 I=1,8 590 MGASH=NGASH+1 591 NGASH=MGASH+1 592 BMATX(1,MGASH)=CARTD(1,I) 593 BMATX(1,NGASH)=0.0 594 BMATX(2,MGASH)=0.0 595 BMATX(2,NGASH)=CARTD(2,I) 596 BMATX(3,MGASH)=CARTD(2,I) 597 BMATX(3,NGASH)=CARTD(1,I) 598 200 CONTINUE 599 RETURN 600 END 601 C***************************************************************************** 602 C 这是一个计算等效节点载荷的子程序。 * 603 C***************************************************************************** 604 SUBROUTINE LOAD(COORD,LNODS,NLOAD) 605 IMPLICIT REAL*8(A-H,O-Z) 606 CHARACTER*80 LINECHAR,NEWLINECHAR 607 DIMENSION COORD(2,1),LNODS(8,1),NLOAD(1) 608 DIMENSION ELCOD(2,3),ELOAD(16),POSGP(3),WEIGP(3),NOPRS(3) 609 DIMENSION PRESS(3,2),PGASH(2),DGASH(2) 610 COMMON/SOL/NPOIN,NELEM,NTYPE,NMATS 611 COMMON/SD/SHAPE(8),DERIV(2,8) 612 OPEN(8,FILE='ELOAD',FORM='UNFORMATTED') 613 C 以下进行的是高斯积分点位置赋值。 614 POSGP(1)=-DSQRT(0.6D0) 615 POSGP(2)=0.0 616 POSGP(3)=DSQRT(0.6D0) 617 C 以下进行的是加权系数赋值。 618 WEIGP(1)=5.0/9.0 619 WEIGP(2)=8.0/9.0 620 WEIGP(3)=5.0/9.0 621 C 以下进行的是输入两行单元受载信息标题。 622 READ(5,5000) LINECHAR 623 LOCATECHAR=INDEX(LINECHAR,'输入') 624 IF(LOCATECHAR.NE.0) THEN 625 LINECHAR(LOCATECHAR:LOCATECHAR+3)='输出' 626 ENDIF 627 WRITE(6,5000) LINECHAR 628 READ(5,5000) LINECHAR 629 LOCATECHAR=INDEX(LINECHAR,'输入') 630 IF(LOCATECHAR.NE.0) THEN 631 LINECHAR(LOCATECHAR:LOCATECHAR+3)='输出' 632 ENDIF 633 WRITE(6,5000) LINECHAR 634 LOCATECHAR=INDEX(LINECHAR,'=') 635 LOCATECHAR=LOCATECHAR+1 636 NEWLINECHAR=LINECHAR(LOCATECHAR:80) 637 READ(NEWLINECHAR,5200) NUMBER 638 5000 FORMAT(A) 639 5200 FORMAT(I5) 640 C 以下进行的是单元受载边数赋初值零。 641 DO 10 I=1,NELEM 642 10 NLOAD(I)=0 643 C 以下进行的是输入受载网格单元的编号及其受载边数。 644 755 FORMAT(2X,'网格单元编号为：',I4,', 单元受载边数为：',I3) 645 DO 20 I=1,NUMBER 646 READ(5,*) KEL,NEDGE 647 NLOAD(KEL)=NEDGE 648 WRITE(6,755) KEL,NLOAD(KEL) 649 20 CONTINUE 650 DO 200 IELEM=1,NELEM 651 NEDGE=NLOAD(IELEM) 652 C 如果当前单元受载边数为零,则转移到循环的终端语句进行下一个网格单元。 653 IF(NEDGE.EQ.0) GO TO 200 654 C 以下进行的是输入单元受载信息标题。 655 READ(5,5000) LINECHAR 656 LOCATECHAR=INDEX(LINECHAR,'输入') 657 IF(LOCATECHAR.NE.0) THEN 658 LINECHAR(LOCATECHAR:LOCATECHAR+3)='输出' 659 ENDIF 660 WRITE(6,5000) LINECHAR 661 C 以下进行的是单元等效节点载荷列阵赋初值零。 662 DO 50 I=1,16 663 ELOAD(I)=0.0 664 50 CONTINUE 665 C 以下进行的是输入受载边三个节点的总体编号及当前受载边均布载荷σ和τ。 666 DO 160 IEDGE=1,NEDGE 667 READ(5,*) (NOPRS(I),I=1,3) 668 READ(5,*) SGMAR,TAU 669 C 如果均布载荷为零,则输入三个节点的法向载荷σ(1,2,3)及切向载荷τ(1,2,3)。 670 IF(ABS(SGMAR).LT.1.D-8.AND.ABS(TAU).LT.1.D-8) THEN 671 READ(5,*) ((PRESS(J,I),J=1,3),I=1,2) 672 ELSE 673 C 如果均布载荷不为零,则给三个节点赋值。 674 DO 60 I=1,3 675 PRESS(I,1)=SGMAR 676 PRESS(I,2)=TAU 677 60 CONTINUE 678 END IF 679 C 给当前受载边三个节点坐标赋值。 680 DO 100 I=1,3 681 LNODE=NOPRS(I) 682 DO 100 J=1,2 683 ELCOD(J,I)=COORD(J,LNODE) 684 100 CONTINUE 685 T=-1.0 686 DO 150 IGAUS=1,3 687 S=POSGP(IGAUS) 688 CALL SFR2(S,T) 689 DO 110 J=1,2 690 PGASH(J)=0.0 691 DGASH(J)=0.0 692 C 以下进行的是计算当前积分点的σ和τ及偏导数在当前积分点的值。 693 DO 110 I=1,3 694 PGASH(J)=PGASH(J)+PRESS(I,J)*SHAPE(I) 695 DGASH(J)=DGASH(J)+ELCOD(J,I)*DERIV(1,I) 696 110 CONTINUE 697 DVOLU=WEIGP(IGAUS) 698 PXCOM=DGASH(1)*PGASH(2)-DGASH(2)*PGASH(1) 699 PYCOM=DGASH(1)*PGASH(1)+DGASH(2)*PGASH(2) 700 C 查找当前受载边的第I个节点是当前单元的第几个节点,查找结果为第I个节点。 701 DO 120 I=1,8 702 NLOCA=LNODS(I,IELEM) 703 IF(NLOCA.EQ.NOPRS(1)) GO TO 130 704 120 CONTINUE 705 130 J=I+2 706 KOUNT=0 707 DO 140 K=I,J 708 KOUNT=KOUNT+1 709 N=(K-1)*2+1 710 M=N+1 711 C 如果单元节点局部编号大于8,则应等于1。 712 IF(K.GT.8)THEN 713 N=1 714 M=2 715 END IF 716 C 以下进行的是累加得到当前单元的等效节点载荷。 717 ELOAD(N)=ELOAD(N)+SHAPE(KOUNT)*PXCOM*DVOLU 718 ELOAD(M)=ELOAD(M)+SHAPE(KOUNT)*PYCOM*DVOLU 719 140 CONTINUE 720 150 CONTINUE 721 160 CONTINUE 722 WRITE(6,795) ELOAD 723 795 FORMAT(2X,'单元等效载荷列阵为：'/4(1P4E15.5/)) 724 WRITE(8)ELOAD 725 200 CONTINUE 726 RETURN 727 END 728 C***************************************************************************** 729 C 这是一个形成整体刚度矩阵和整体载荷列阵的单元组集子程序。 * 730 C***************************************************************************** 731 SUBROUTINE ASEM(V,A,LNODS,IFPRE,NLOAD,MAXA) 732 IMPLICIT REAL*8(A-H,O-Z) 733 DIMENSION V(1),A(1) 734 DIMENSION LNODS(8,1),IFPRE(1),NLOAD(1),MAXA(1) 735 DIMENSION NCODF(2),ESTIF(136),ELOAD(16) 736 COMMON/SOL/NPOIN,NELEM,NTYPE,NMATS 737 MAXA(1)=1 738 DO 100 IPOIN=1,NPOIN 739 KMIN=NPOIN+1 740 DO 40 IELEM=1,NELEM 741 DO 10 I=1,8 742 LNODE=LNODS(I,IELEM) 743 IF(LNODE.EQ.IPOIN)GO TO 20 744 10 CONTINUE 745 GO TO 40 746 20 CONTINUE 747 DO 30 I=1,8 748 LNODE=LNODS(I,IELEM) 749 IF(LNODE.GE.KMIN)GO TO 30 750 KMIN=LNODE 751 30 CONTINUE 752 40 CONTINUE 753 KHIGH=(IPOIN-KMIN)*2 754 DO 60 J=1,2 755 KDOFN=(IPOIN-1)*2+J+1 756 MAXA(KDOFN)=MAXA(KDOFN-1)+KHIGH+J 757 60 CONTINUE 758 100 CONTINUE 759 NN=NPOIN*2 760 NNM=NN+1 761 NWK=MAXA(NNM)-1 762 WRITE(6,735) NWK 763 735 FORMAT(2X,'剖面元素之和=',I7) 764 DO 120 IWK=1,NWK 765 A(IWK)=0.0 766 120 CONTINUE 767 DO 130 IN=1,NN 768 130 V(IN)=0.0 769 DLARG=1.0D+15 770 REWIND 8 771 DO 300 IELEM=1,NELEM 772 WRITE(6,135) IELEM 773 135 FORMAT(2X,'网格单元编号为：',I5) 774 READ(7) ESTIF 775 IF(NLOAD(IELEM).GT.0)THEN 776 READ(8) ELOAD 777 END IF 778 DO 280 INODE=1,8 779 LNODI=LNODS(INODE,IELEM) 780 ICODF=IFPRE(LNODI) 781 NCODF(1)=ICODF/10 782 NCODF(2)=ICODF-NCODF(1)*10 783 DO 260 IDOFN=1,2 784 IEVAB=(INODE-1)*2+IDOFN 785 ICOLN=(LNODI-1)*2+IDOFN 786 DO 240 JNODE=1,8 787 LNODJ=LNODS(JNODE,IELEM) 788 DO 220 JDOFN=1,2 789 JEVAB=(JNODE-1)*2+JDOFN 790 JROW=(LNODJ-1)*2+JDOFN 791 IF(JROW.GT.ICOLN)GO TO 220 792 IF(JEVAB.GT.IEVAB)THEN 793 KEVAB=JEVAB*(JEVAB-1)/2+IEVAB 794 ELSE 795 KEVAB=IEVAB*(IEVAB-1)/2+JEVAB 796 END IF 797 LDIS=MAXA(ICOLN)+(ICOLN-JROW) 798 A(LDIS)=A(LDIS)+ESTIF(KEVAB) 799 220 CONTINUE 800 240 CONTINUE 801 IF(NLOAD(IELEM).GT.0)THEN 802 V(ICOLN)=V(ICOLN)+ELOAD(IEVAB) 803 END IF 804 IF(NCODF(IDOFN).EQ.0)GO TO 260 805 KEVAB=IEVAB*(IEVAB+1)/2 806 LDIS=MAXA(ICOLN) 807 A(LDIS)=A(LDIS)+DLARG*ESTIF(KEVAB) 808 260 CONTINUE 809 280 CONTINUE 810 300 CONTINUE 811 C******* 812 C ENTER GIVED DISPLACEMENTS 813 C******* 814 DO 400 I=1,NPOIN 815 ICODF=IFPRE(I) 816 NCODF(1)=ICODF/10 817 NCODF(2)=ICODF-NCODF(1)*10 818 DO 320 J=1,2 819 IF(NCODF(J).EQ.2)THEN 820 C******* 821 READ(5,*) X 822 C******* 823 ICOLN=(I-1)*2+J 824 LDIS=MAXA(ICOLN) 825 V(ICOLN)=A(LDIS)*X 826 END IF 827 320 CONTINUE 828 400 CONTINUE 829 RETURN 830 END 831 C***************************************************************************** 832 C 这是一个求解结构平衡方程组的子程序。 * 833 C***************************************************************************** 834 SUBROUTINE SOLVE(V,A,MAXA) 835 IMPLICIT REAL*8(A-H,O-Z) 836 DIMENSION V(1),A(1),MAXA(1) 837 COMMON/SOL/NPOIN,NELEM,NTYPE,NMATS 838 K=0 839 NN=NPOIN*2 840 NNM=NN+1 841 NWK=MAXA(NNM)-1 842 DO 640 N=1,NN 843 KN=MAXA(N) 844 KL=KN+1 845 KU=MAXA(N+1)-1 846 KH=KU-KL 847 IF(KH)610,590,550 848 550 K=N-KH 849 IC=0 850 KLT=KU 851 DO 580 J=1,KH 852 IC=IC+1 853 KLT=KLT-1 854 KI=MAXA(K) 855 ND=MAXA(K+1)-KI-1 856 IF(ND)580,580,560 857 560 KK=MIN(IC,ND) 858 C=0 859 DO 570 L=1,KK 860 M1=KI+L 861 M2=KLT+L 862 570 C=C+A(M1)*A(M2) 863 A(KLT)=A(KLT)-C 864 580 K=K+1 865 590 K=N 866 B=0.0 867 DO 600 KK=KL,KU 868 K=K-1 869 KI=MAXA(K) 870 C=A(KK)/A(KI) 871 B=B+C*A(KK) 872 600 A(KK)=C 873 A(KN)=A(KN)-B 874 610 IF(A(KN)) 620,620,640 875 620 WRITE(6,2000) N,A(KN) 876 STOP 877 640 CONTINUE 878 C******** 879 C REDUCE RIGHT--HAND--SIDE LOAD VECTOR 880 C******** 881 650 CONTINUE 882 DO 680 N=1,NN 883 KL=MAXA(N)+1 884 KU=MAXA(N+1)-1 885 KUL=KU-KL 886 IF(KUL)680,660,660 887 660 K=N 888 C=0.0 889 DO 670 KK=KL,KU 890 K=K-1 891 670 C=C+A(KK)*V(K) 892 V(N)=V(N)-C 893 680 CONTINUE 894 C******* 895 C BACK-SUBSTITUTE 896 C******** 897 DO 700 N=1,NN 898 K=MAXA(N) 899 V(N)=V(N)/A(K) 900 700 CONTINUE 901 IF(NN.EQ.1)GO TO 740 902 N=NN 903 DO 730 L=2,NN 904 KL=MAXA(N)+1 905 KU=MAXA(N+1)-1 906 KUL=KU-KL 907 IF(KUL)730,710,710 908 710 K=N 909 DO 720 KK=KL,KU 910 K=K-1 911 720 V(K)=V(K)-A(KK)*V(N) 912 730 N=N-1 913 740 CONTINUE 914 2000 FORMAT(//10X,'***** 停止运行！系数矩阵非正定！*****'//'NON 915 + POSITIVE PIVOT FOR EQUATION'I4//' PIVOT=',E20.12) 916 C******** 917 C OUTPUT DISPLACEMENT 918 C******** 919 WRITE(6,903) 920 903 FORMAT(//80('=')/36X,'移动输出'/80('=')//) 921 WRITE(6,905) 922 905 FORMAT(32X,'Ｘ方向',10X,'Ｙ方向') 923 M2=0 924 DO 1000 I=1,NPOIN 925 M1=M2+1 926 M2=M1+1 927 1000 WRITE(6,915)I,(V(J),J=M1,M2) 928 915 FORMAT(2X,'高斯积分点编号为：',I5,1PE16.5,1PE16.5) 929 RETURN 930 END 931 C***************************************************************************** 932 C 这是一个计算单元应力的子程序。 * 933 C***************************************************************************** 934 SUBROUTINE STRE(COORD,PROPS,V,LNODS,MATNO) 935 IMPLICIT REAL*8(A-H,O-Z) 936 CHARACTER*80 LINECHAR 937 DIMENSION COORD(2,1),PROPS(3,1),V(1),LNODS(8,1),MATNO(1) 938 DIMENSION POSGP(3),ELCOD(2,8) 939 DIMENSION ELDIS(16),BE(3),S(3),SP(2) 940 COMMON/SOL/NPOIN,NELEM,NTYPE,NMATS 941 COMMON/BMDMP/BMATX(3,16),DMATX(3,3) 942 COMMON/GCOD/GPCOD(2,9) 943 POSGP(1)=-DSQRT(0.6D0) 944 POSGP(2)=0.0 945 POSGP(3)=DSQRT(0.6D0) 946 WRITE(6,616) 947 616 FORMAT(80('=')/26X,'以下是高斯积分点的应力输出'/80('=')/) 948 DO 90 IELEM=1,NELEM 949 WRITE(6,636) IELEM 950 636 FORMAT(2X,'网格单元编号为：',I5) 951 LPROP=MATNO(IELEM) 952 DO 20 I=1,8 953 LNODE=LNODS(I,IELEM) 954 LDOFN=LNODE*2-2 955 DO 10 J=1,2 956 IEVAB=I*2-2+J 957 MDOFN=LDOFN+J 958 ELDIS(IEVAB)=V(MDOFN) 959 ELCOD(J,I)=COORD(J,LNODE) 960 10 CONTINUE 961 20 CONTINUE 962 YOUNG=PROPS(1,LPROP) 963 POISS=PROPS(2,LPROP) 964 THICK=PROPS(3,LPROP) 965 CALL MODPS(YOUNG,POISS) 966 KGASP=0 967 DO 80 IGAUS=1,3 968 DO 80 JGAUS=1,3 969 KGASP=KGASP+1 970 EXISP=POSGP(IGAUS) 971 ETASP=POSGP(JGAUS) 972 CALL SFR2(EXISP,ETASP) 973 CALL JACOB2(IELEM,DJACB,KGASP,ELCOD) 974 CALL BMATPS 975 DO 40 I=1,3 976 BE(I)=0.0 977 DO 30 J=1,16 978 BE(I)=BE(I)+BMATX(I,J)*ELDIS(J) 979 30 CONTINUE 980 40 CONTINUE 981 DO 60 L=1,3 982 S(L)=0.0 983 DO 50 J=1,3 984 S(L)=S(L)+DMATX(L,J)*BE(J) 985 50 CONTINUE 986 60 CONTINUE 987 WRITE(6,755) KGASP 988 755 FORMAT(/2X,'高斯积分点的编号为：',I2) 989 LINECHAR='高斯积分点的坐标为：' 990 WRITE(6,775) LINECHAR,GPCOD(1,KGASP),GPCOD(2,KGASP) 991 775 FORMAT(2X,A26,' X=',F9.3,5X,'Y=',F9.3) 992 C ***** 993 XA=(S(1)+S(2))*0.5 994 XI=(S(1)-S(2))*0.5 995 XE=S(3) 996 X0=DSQRT(XI*XI+XE*XE) 997 SP(1)=XA+X0 998 SP(2)=XA-X0 999 C ****** 1000 WRITE(6,815) 1001 815 FORMAT(14X,'正应力σx',6X,'正应力σy',7X,'剪应力τ') 1002 WRITE(6,835) S 1003 835 FORMAT(10X,1PE13.4,2X,1PE13.4,2X,1PE13.4) 1004 WRITE(6,855) 1005 855 FORMAT(14X,'主应力σ1',6X,'主应力σ2') 1006 WRITE(6,875) SP 1007 875 FORMAT(10X,1PE13.4,2X,1PE13.4) 1008 80 CONTINUE 1009 90 CONTINUE 1010 RETURN 1011 END 1012