内容发布更新时间 : 2024/12/23 21:52:45星期一 下面是文章的全部内容请认真阅读。

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 \ 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=I-K+1

C=AK(K,L)/AK(K,1) IW=N_BAND-L+1 DO 10 J=1,IW