fortran.zip_Fortran

PROGRAM FEM2D DIMENSION IJK_ELE(500,3),X(500),Y(500),IJK_U(50,3),P_IJK(50,3), &RESULT_N(500),AK(500,100) DIMENSION STS_ELE(500,3),STS_ND(500,3) OPEN(4,FILE='BASIC.IN') OPEN(5,FILE='NODE.IN') OPEN(6,FILE='ELEMENT.IN') OPEN(8,FILE='DATA.OUT') OPEN(9,FILE='FOR_POST.DAT') READ(4,*)ID,N_ELE,N_NODE,N_BC,N_LOAD IF(ID.EQ.1)WRITE(8,20) IF(ID.EQ.2)WRITE(8,25) 20 FORMAT(/5X,'=========PLANE STRESS PROBLEM========') 25 FORMAT(/5X,'=========PLANE STRAIN PROBLEM========') CALL READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR,PT, & IJK_ELE,X,Y,IJK_U,P_IJK) CALL BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE, & IJK_ELE,X,Y,PE,PR,PT,AK) CALL FORM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK, & RESULT_N) CALL DO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N) CALL SOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N) CALL CAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N, & STS_ELE,STS_ND) c to putout a data file WRITE(9,70)REAL(N_NODE),REAL(N_ELE) 70 FORMAT(2f9.4) WRITE(9,71)(X(I),Y(I),RESULT_N(2*I-1),RESULT_N(2*I), & STS_ND(I,1),STS_ND(I,2),STS_ND(I,3),I=1,N_NODE) 71 FORMAT(7F9.4) WRITE(9,72)(REAL(IJK_ELE(I,1)),REAL(IJK_ELE(I,2)), &REAL(IJK_ELE(I,3)),REAL(IJK_ELE(I,3)), &STS_ELE(I,1),STS_ELE(I,2),STS_ELE(I,3),I=1, N_ELE) 72 FORMAT(7f9.4) c CLOSE(4) CLOSE(5) CLOSE(6) CLOSE(8) CLOSE(9) END c c to get the original data in order to model the problem SUBROUTINE READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR, &PT,IJK_ELE,X,Y,IJK_U,P_IJK) DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),IJK_U(N_BC,3), & P_IJK(N_LOAD,3),NE_ANSYS(N_ELE,14) REAL ND_ANSYS(N_NODE,3) READ(4,*)PE,PR,PT READ(4,*)((IJK_U(I,J),J=1,3),I=1,N_BC) READ(4,*)((P_IJK(I,J),J=1,3),I=1,N_LOAD) READ(5,*)((ND_ANSYS(I,J),J=1,3),I=1,N_NODE) READ(6,*)((NE_ANSYS(I,J),J=1,14),I=1,N_ELE) DO 10 I=1,N_NODE X(I)=ND_ANSYS(I,2) Y(I)=ND_ANSYS(I,3) 10 CONTINUE DO 11 I=1,N_ELE DO 11 J=1,3 IJK_ELE(I,J)=NE_ANSYS(I,J) 11 CONTINUE N_BAND=0 DO 20 IE=1,N_ELE DO 20 I=1,3 DO 20 J=1,3 IW=IABS(IJK_ELE(IE,I)-IJK_ELE(IE,J)) IF(N_BAND.LT.IW)N_BAND=IW 20 CONTINUE N_BAND=(N_BAND+1)*2 IF(ID.EQ.1) THEN ELSE PE=PE/(1.0-PR*PR) PR=PR/(1.0-PR) END IF RETURN END c C to form the stiffness matrix of element SUBROUTINE FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE) DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),BB(3,6),DD(3,3), & AKE(6,6), SS(6,6) CALL CAL_DD(PE,PR,DD) CALL CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB) DO 10 I=1,3 DO 10 J=1,6 SS(I,J)=0.0 DO 10 K=1,3 10 SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J) DO 20 I=1,6 DO 20 J=1,6 AKE(I,J)=0.0 DO 20 K=1,3 20 AKE(I,J)=AKE(I,J)+SS(K,I)*BB(K,J)*AE*PT RETURN END c c to form banded global stiffness matrix SUBROUTINE BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,IJK_ELE,X,Y,PE, & PR,PT,AK) DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),AKE(6,6),AK(500,100) N_DOF=2*N_NODE DO 40 I=1,N_DOF DO 40 J=1,N_BAND 40 AK(I,J)=0 DO 50 IE=1,N_ELE CALL FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE) DO 50 I=1,3 DO 50 II=1,2 IH=2*(I-1)+II IDH=2*(IJK_ELE(IE,I)-1)+II DO 50 J=1,3 DO 50 JJ=1,2 IL=2*(J-1)+JJ IZL=2*(IJK_ELE(IE,J)-1)+JJ IDL=IZL-IDH+1 IF(IDL.LE.0) THEN ELSE AK(IDH,IDL)=AK(IDH,IDL)+AKE(IH,IL) END IF 50 CONTINUE RETURN END c c to calculate the area of element SUBROUTINE CAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE) DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE) I=IJK_ELE(IE,1) J=IJK_ELE(IE,2) K=IJK_ELE(IE,3) XIJ=X(J)-X(I) YIJ=Y(J)-Y(I) XIK=X(K)-X(I) YIK=Y(K)-Y(I) AE=(XIJ*YIK-XIK*YIJ)/2.0 RETURN END c c to calculate the elastic matrix of element SUBROUTINE CAL_DD(PE,PR,DD) DIMENSION DD(3,3) DO 10 I=1,3 DO 10 J=1,3 10 DD(I,J)=0.0 DD(1,1)=PE/(1.0-PR*PR) DD(1,2)=PE*PR/(1.0-PR*PR) DD(2,1)=DD(1,2) DD(2,2)=DD(1,1) DD(3,3)=PE/((1.0+PR)*2.0) RETURN END c c to calculate the strain-displacement matrix of element SUBROUTINE CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB) DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),BB(3,6) I=IJK_ELE(IE,1) J=IJK_ELE(IE,2) K=IJK_ELE(IE,3) DO 10 II=1,3 DO 10 JJ=1,3 10 BB(II,JJ)=0.0 BB(1,1)=Y(J)Y(K) BB(1,3)=Y(K)-Y(I) BB(1,5)=Y(I)-Y(J) BB(2,2)=X(K)-X(J) BB(2,4)=X(I)- X(K) BB(2,6)=X(J)-X(I) BB(3,1)=BB(2,2) BB(3,2)=BB(1,1) BB(3,3)=BB(2,4) BB(3,4)=BB(1,3) BB(3,5)=BB(2,6) BB(3,6)=BB(1,5) CALL CAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE) DO 20 I1=1,3 DO 20 J1=1,6 20 BB(I1,J1)=BB(I1,J1)/(2.0*AE) RETURN END c c to form the global load matrix SUBROUTINE FORM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK, & RESULT_N) DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),P_IJK(N_LOAD,3), & RESULT_N(N_DOF) DO 10 I=1,N_DOF 10 RESULT_N(I)=0.0 DO 20 I=1,N_LOAD II=P_IJK(I,1) RESULT_N(2*II-1)=P_IJK(I,2) 20 RESULT_N(2*II)=P_IJK(I,3) RETURN END c c to deal with BC(u) (here only for fixed displacement) using "1-0" method SUBROUTINE DO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N) DIMENSION RESULT_N(N_DOF),IJK_U(N_BC,3),AK(500,100) DO 30 I=1,N_BC IR=IJK_U(I,1) DO 30 J=2,3 IF(IJK_U(I,J).EQ.0)THEN ELSE II=2*IR+J-3 AK(II,1)=1.0 RESULT_N(II)=0.0 DO 10 JJ=2,N_BAND 10 AK(II,JJ)=0.0 DO 20 JJ=2,II 20 AK(II-JJ+1,JJ)=0.0 END IF 30 CONTINUE RETURN END c c to solve the banded FEM equation by GAUSS elimination SUBROUTINE SOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N) DIMENSION RESULT_N(N_DOF),AK(500,100) DO 20 K=1,N_DOF-1 IF(N_DOF.GT.K+N_BAND-1)IM=K+N_BAND-1 IF(N_DOF.LE.K+N_BAND-1)IM=N_DOF DO 20 I=K+1,IM L=IK+1 C=AK(K,L)/AK(K,1) IW=N_BAND-L+1 DO 10 J=1,IW M=J+I-K 10 AK(I,J)=AK(I,J)-C*AK(K,M) 20 RESULT_N(I)=RESULT_N(I)-C*RESULT_N(K) RESULT_N(N_DOF)=RESULT_N(N_DOF)/AK(N_DOF,1) DO 40 I1=1,N_DOF-1 I=N_DOF-I1 IF(N_BAND.GT.N_DOF-I-1)JQ=N_DOF-I+1 IF(N_BAND.LE.N_DOF-I-1)JQ=N_BAND DO 30 J=2,JQ K=J+I-1 30 RESULT_N(I)=RESULT_N(I)-AK(I,J)*RESULT_N(K) 40 RESULT_N(I)=RESULT_N(I)/AK(I,1) WRITE(8,50) 50 FORMAT(/12X,'* * * * * RESULTS BY FEM2D * * * * *',//8X, &'--DISPLACEMENT OF NODE--'//5X,'NODE NO',8X,'X-DISP',8X,'Y-DISP') DO 60 I=1,N_NODE 60 WRITE(8,70) I,RESULT_N(2*I-1),RESULT_N(2*I) 70 FORMAT(8X,I5,7X,2E15.6) RETURN END c c calculate the stress components of element and node SUBROUTINE CAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N, &STS_ELE,STS_ND) DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),DD(3,3),BB(3,6), &SS(3,6),RESULT_N(N_DOF),DISP_E(6) DIMENSION STS_ELE(500,3),STS_ND(500,3) WRITE(8,10) 10 FORMAT(//8X,'--STRESSES OF ELEMENT--') CALL CAL_DD(PE,PR,DD) DO 50 IE=1,N_ELE CALL CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB) DO 20 I=1,3 DO 20 J=1,6 SS(I,J)=0.0 DO 20 K=1,3 20 SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J) DO 30 I=1,3 DO 30 J=1,2 IH=2*(I-1)+J IW=2*(IJK_ELE(IE,I)-1)+J 30 DISP_E(IH)=RESULT_N(IW) STX=0 STY=0 TXY=0 DO 40 J=1,6 STX=STX+SS(1,J)*DISP_E(J) STY=STY+SS(