有限元分析程序与算例
WRITE (*,5)
5 FORMAT(/10X,'The Joint Coordinates'//1X,'joint',11X,'X',17X,'Y') WRITE (*,6)(M,X(M),Y(M),M=1,NJ) 6 FORMAT(1X,I4,2F18.6) READ (3,*)(IS(M),VS(M),M=1,NS) WRITE (*,7)
7 FORMAT(/10X,'The Information of Supports'//4X,'IS',9X,'VS') WRITE (*,8)(IS(M),VS(M),M=1,NS) 8 FORMAT(1X,I5,F16.6) RETURN END
! Determin location vector of element SUBROUTINE LCVCT (NM,IST,IEN,LV,NJ,N) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION IST(NM),IEN(NM),LV(6,NM) DO 100 M=1,NM I=IST(M)*3 J=IEN(M)*3 LV(1,M)=I-2 LV(2,M)=I-1 LV(3,M)=I LV(4,M)=J-2 LV(5,M)=J-1 LV(6,M)=J 100 CONTINUE N=NJ*3 RETURN END
! Determine location of diagonal elements of global stiffness SUBROUTINE LCDIA (NM,N,LV,MIN,IBDW,LD,MAXBDW,NA) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION LV(6,NM),MIN(N),IBDW(N),LD(N) DO 100 I=1,N MIN(I)=I 100 CONTINUE DO 400 M=1,NM MINLV=LV(1,M) DO 200 I=2,6 IF (LV(I,M).LT.MINLV) MINLV=LV(I,M) 200 CONTINUE DO 300 I=1,6 IF (MINLV.LT.MIN(LV(I,M))) MIN(LV(I,M))=MINLV
- 3 -
有限元分析程序与算例
300 CONTINUE 400 CONTINUE MAXBDW=0 DO 500 I=1,N IBDW(I)=I-MIN(I)+1 IF (IBDW(I).GT.MAXBDW) MAXBDW=IBDW(I) 500 CONTINUE LD(1)=IBDW(1) DO 600 I=2,N LD(I)=LD(I-1)+IBDW(I) 600 CONTINUE NA=LD(N) RETURN END
! Calculate length, cosine & sine of member SUBROUTINE RLCS (M,NM,NJ,IST,IEN,X,Y,RL,C,S) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION IST(NM),IEN(NM),X(NJ),Y(NJ) I=IST(M) J=IEN(M) X1=X(J)-X(I) Y1=Y(J)-Y(I) RL=SQRT(X1*X1+Y1*Y1) C=X1/RL S=Y1/RL RETURN END
! Calculate element stiffness matrix along local axes SUBROUTINE KEBAR (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,C,S,E1,E2,E3,E4) IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION IST(NM),IEN(NM),X(NJ),Y(NJ),AR(NM),RI(NM) CALL RLCS (M,NM,NJ,IST,IEN,X,Y,RL,C,S) E1=E*AR(M)/RL E2=12.0*E*RI(M)/(RL*RL*RL) E3=0.5*E2*RL E4=0.6666667*E3*RL RETURN END
! Calculate element stiffness matrix along global axes SUBROUTINE KE (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,AE) IMPLICIT DOUBLE PRECISION (A-H,O-Z)
- 4 -
有限元分析程序与算例
DIMENSION IST(NM),IEN(NM),AR(NM),RI(NM),X(NJ),Y(NJ),AE(6,6) CALL KEBAR (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,C,S,E1,E2,E3,E4) A1=E1*C*C+E2*S*S A2=(E1-E2)*C*S A3=E1*S*S+E2*C*C A4=E3*S A5=E3*C A6=E4 AE(1,1)=A1 AE(2,1)=A2 AE(2,2)=A3 AE(3,1)=-A4 AE(3,2)=A5 AE(3,3)=A6 AE(4,1)=-A1 AE(4,2)=-A2 AE(4,3)=A4 AE(4,4)=A1 AE(5,1)=-A2 AE(5,2)=-A3 AE(5,3)=-A5 AE(5,4)=A2 AE(5,5)=A3 AE(6,1)=-A4 AE(6,2)=A5 AE(6,3)=0.5*A6 AE(6,4)=A4 AE(6,5)=-A5 AE(6,6)=A6 RETURN END
! Form global stiffness matrix A SUBROUTINE FORMA (E,NM,NJ,N,NA,IST,IEN,AR,RI,X,Y,LV,AE,LD,A) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION
IST(NM),IEN(NM),AR(NM),RI(NM),X(NJ),Y(NJ),LV(6,NM),AE(6,6),LD(N) DOUBLE PRECISION A(NA) DO 300 M=1,NM CALL KE (M,E,NM,NJ,IST,IEN,AR,RI,X,Y,AE) DO 200 I=1,6 DO 100 J=1,I IF (LV(I,M).GE.LV(J,M)) THEN A(LD(LV(I,M))-LV(I,M)+LV(J,M))=A(LD(LV(I,M))-LV(I,M)+LV(J,M))+AE(I,J)
- 5 -
有限元分析程序与算例
ELSE A(LD(LV(J,M))-LV(J,M)+LV(I,M))=A(LD(LV(J,M))-LV(J,M)+LV(I,M))+AE(I,J) ENDIF
100 CONTINUE 200 CONTINUE 300 CONTINUE RETURN END
! Introduce support conditions into global stiffness matrix A SUBROUTINE AS(NS,N,NA,IS,LD,A)
IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION IS(NS),LD(N),A(NA) DO 100 M=1,NS I=3*(IS(M)/10)-3+MOD(IS(M),10) A(LD(I))=1D22 100 CONTINUE RETURN END
! Solve equations(1)-decomposition of matrix A SUBROUTINE LDLT(N,NA,LD,A,T)
IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION LD(N),A(NA),T(N) DOUBLE PRECISION SUM DO 300 I=2,N LDI=LD(I) I1=I-LDI+LD(I-1)+1 DO 200 J=I1,I-1 LDJ=LD(J) J1=J-LDJ+LD(J-1)+1 IF(I1.GT.J1) J1=I1 SUM=0.0D0 DO 100 K=J1,J-1 SUM=SUM+T(K)*A(LDJ-J+K) 100 CONTINUE T(J)=A(LDI-I+J)-SUM A(LDI-I+J)=T(J)/A(LDJ) A(LDI)=A(LDI)-T(J)*A(LDI-I+J) 200 CONTINUE 300 CONTINUE RETURN END
! Solve equations(2)-forward & back substitution SUBROUTINE SLVEQ (N,NA,MAXBDW,LD,A,B)
- 6 -
有限元分析程序与算例
IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION LD(N),A(NA),B(N) DO 200 I=2,N LDI=LD(I) I1=I-LDI+LD(I-1)+1 DO 100 J=I1,I-1 B(I)=B(I)-A(LDI-I+J)*B(J) 100 CONTINUE 200 CONTINUE DO 300 I=1,N B(I)=B(I)/A(LD(I)) 300 CONTINUE DO 500 I=N-1,1,-1 IMIN=I+MAXBDW IF(IMIN.GT.N) IMIN=N DO 400 J=I+1,IMIN LDJ=LD(J) J1=J-LDJ+LD(J-1)+1 IF(I.GE.J1) B(I)=B(I)-A(LDJ-J+I)*B(J) 400 CONTINUE 500 CONTINUE RETURN END
! Initialize joint load vector B SUBROUTINE B0 (LC,N,NLJ,B)
IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION B(N) WRITE (*,1)LC,NLJ
1 FORMAT(/15X,'Loading Case',I3 //10X,'The Loadings at Joints'//17X,'NLJ=',I4) DO 100 I=1,N B(I)=0.0D0 100 CONTINUE RETURN END
! Read data of loading at joint & form joint load vector B SUBROUTINE IOLJB (N,NLJ,ILJ,PX,PY,PM,B) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION ILJ(NLJ),PX(NLJ),PY(NLJ),PM(NLJ),B(N) READ (3,*)(ILJ(M),PX(M),PY(M),PM(M),M=1,NLJ) WRITE (*,1)
1 FORMAT(/2X,'ILJ',11X,'PX',14X,'PY',15X,'PM') WRITE (*,2)(ILJ(M),PX(M),PY(M),PM(M),M=1,NLJ) 2 FORMAT(1X,I4,2F16.4,F18.5) DO 100 M=1,NLJ
- 7 -