BEAMGN平面梁大挠度分析程序.doc

BEAMGN平面梁大挠度分析程序 BEAMGN平面梁大挠度分析程序 * * THE NONLINEAR FINITE ELEMENT ANALYSIS * * PROGRAM FOR BEAM T.L. AND U.L. METHOD * * PROGRAM MAIN * * 变量说明 * FILE1: 输入文件名 * FILE2: 输出文件名 * FILE3: 图形结果数据文件名 * NTU: 计算方法指示数，1 为 T.L.，2 为 U.L. * LOAD：载荷因子数 * EF: 拉伸刚度 * EI: 抗弯刚度 * NE: 单元总数 * IME: 增量步总数 * ACC: 计算精度 * MM: 迭代信息数，1 为等刚度，2 为变刚度 * NRM: 广义外力作用个数 * NU: 初始位移约束个数 * XX1,YY1,XX2,YY2: 梁两端点坐标 * PNM: 广义力信息，PNM(I,1),PNM(I,2),PNM(I,3)分别为作用点号，方向和数值 * UXY: 结点位移约束信息，UXY(I,1),UXY(I,2),UXY(I,3)为作用点号，方向和数值 * M: 结构自由度数 * NE: 单元总数 * BO：单元两结点连线与坐标方向夹角 * PP: 总结点载荷数组 * U: 总位移数组 * P: 结点增量载荷数组 * IMPLICIT REAL*8 (A-H,O-Z) DIMENSION WK(63,63),U(63),DU(63),RS(63,1),PP(63,1),RN(20), & XY(21,2),DLL(20),BO(20),XYO(21,2),TXY(20,2),TL1(20), & PNM(10,4),UXY(10,3),FP(63,1),P(63,1),UI(63),TL2(20), & PP1(63,1),UU(63) * INPUT DATA * CHARACTER FILE1*12,FILE2*12 WRITE(*,*)PLEASE INPUT THE INPUT FILE NAME READ (*,(A) FILE1 OPEN (1, FILE=FILE1, STATUS=OLD) OUT=2 IF(OUT.EQ.2) THEN WRITE(*,*)PLEASE INPUT THE OUTPUT FILE NAME READ (*,(A) FILE2 OPEN (2, FILE=FILE2, STATUS=NEW) WRITE(2,*)- THE RESULTS OF CALCULATION - END IF READ(1,*)NTU,LOAD WRITE(*,*)- THE RESULTS OF CALCULATION - READ(1,*) EF,EI READ(1,*) NE,IME,ACC,MM,NRM,NU NP=NE+1 READ(1,*) XX1,YY1 READ(1,*)XX2,YY2 XXX=(XX2-XX1)/NE YYY=(YY2-YY1)/NE DO 10 I=1,NP XY(I,1)=XXX*(I-1) XY(I,2)=YYY*(I-1) 10 CONTINUE DO 11 I=1,NRM READ(1,*) PNM(I,1),PNM(I,2),PNM(I,3) 11 CONTINUE DO 20 I=1,NU READ(1,*) UXY(I,1),UXY(I,2),UXY(I,3) 20 CONTINUE M=3*NP M1=M+1 PI=3.1415926 DO 50 I=1,NE 50 BO(I)=0.0 DO 60 I=1,M PP(I,1)=0.0 U(I)=0.0 60 P(I,1)=0.0 * FORM TOTAL LOAD COLUMN MATRIX * DO 70 K=1,NRM II=3*(PNM(K,1)-1)+PNM(K,2) 70 P(II,1)=PNM(K,3)/IME DO 75 I=1,NP DO 75 J=1,2 75 XYO(I,J)=XY(I,J) NO=0 * FORM LOCAL COORDINATE * 110 CALL FLC(NE,M,RS,UI,DU,RN,XY,BO,DLL) 1500 NO=NO+1 * FORM INCREMENT LOAD PP * IF(LOAD.GT.0) THEN IF(NTU.EQ.2) THEN DO 90 I=1,M 90 PP(I,1)=P(I,1) ELSE DO 95 I=1,M 95 PP(I,1)=P(I,1)*NO END IF ELSE E=M-1 F=M-2 PP1(E,1)=0. PP1(F,1)=0.0 UU(M)=DABS(U(M) IF(NTU.EQ.2) THEN DO 15 I=1,M 15 PP(I,1)=0.0 PP(E,1)=P(E,1)*COS(UU(M) PP(F,1)=P(E,1)*SIN(UU(M) ELSE PP(E,1)=P(E,1)*COS(UU(M)+PP1(E,1) PP(F,1)=P(E,1)*SIN(UU(M)+PP1(F,1) PP1(E,1)=PP(E,1) PP1(F,1)=PP(F,1) END IF END IF NUM=0 1600 CONTINUE * FORM THE ELEMENT STIFFNESS MATRIX WK AND * * ASSEMBLE THE TOTAL STIFFNESS MATRIX TWK * CALL WKK(NE,DLL,RN,BO,WK,EF,EI) 1800 NUM=NUM+1 * CALCULATE THE INEQUILIBRIUM FORCE FP * DO 105 I=1,M 105 FP(I,1)=PP(I,1)-RS(I,1) * SOLVE INCREMENT DISPLACEMENT DU AND TOTAL * * DISPLACEMENT UI WITHIN INCREMENT SETP * CALL JFC(WK,FP,DU,M,M1,NU,UXY) DO 130 I=1,M UI(I)=UI(I)+DU(I) 130 CONTINUE * CALCULATE THE REACTION COLUMN MATRIX RS * CALL RSS(DLL,RS,RN,BO,NE,UI) * JUDGE CONVERGENCY |DU/UI|< 0.001? * AXM=0.0 DO 140 I=1,NP II=3*(I-1)+2 AX=DABS(DU(II)/UI(II) IF(AX.GT.AXM) THEN AXM=AX END IF 140 CONTINUE IF(AXM.GT.ACC) THEN IF(MM.EQ.1) THEN GOTO 1800 ELSE GOTO 1600 END IF END IF * CALCULATE TOTAL DISPLACEMENT U AND * * REVISE NODES COORDINATE XY * IF(NTU.EQ.1) THEN DO 145 I=1,M 145 U(I)=UI(I) ELSE DO 148 I=1,M 148 U(I)=U(I)+UI(I) END IF DO 180 I=1,NP I1=3*(I-1)+1 I2=I1+2 DO 180 J=I1,I2 XY(I,1)=XYO(I,1)+U(I1) XY(I,2)=XYO(I,2)+U(I1+1) 180 CONTINUE III=NO DO 450 J=1,2 450 TXY(III,J)=XY(NP,J) * OUTPUT CALCULATION RESULTS U * IF(OUT.EQ.1) THEN WRITE(*,1000)NO,NUM WRITE(*,1001) DO 200 I=1,NP II=3*(I-1)+1 WRITE(*,1002)I,U(II),U(II+1),U(II+2) 200 CONTINUE ELSE WRITE(2,1000)NO,NUM WRITE(2,1001) DO 300 I=1,NP II=3*(I-1)+1 WRITE(2,1002)I,U(II),U(II+1),U(II+2) 300 CONTINUE WRITE(2,1007)(XY(I,J),J=1,2),I=1,NP) 1007 FORMAT(1X,2E14.2) END IF 1000 FORMAT(/1X,INCREMENT STEP:,I2,/1X,ITERATIVE TIME:,I2) 1001 FORMAT(/1X,NO,11X,U,14X,W,12X,DW/DX) 1002 FORMAT(1X,I2,3F15.6) IF(NO.NE.IME) THEN IF(NTU.EQ.2) THEN GO TO 110 ELSE GOTO 1500 END IF ELSE END IF WRITE(2,*) WRITE(2,*)-Graph- DO 500 I=1,NO TL1(I)=TXY(I,1)/300.0 500 TL2(I)=DABS(TXY(I,2)/300.0) WRITE(2,1005) WRITE(2,1003)(I,TL1(I),I=1,NO) WRITE(2,1006) WRITE(2,1003)(I,TL2(I),I=1,NO) 1003 FORMAT(1X,I2,F14.2) 1005 FORMAT(/1X,NO,5X,(L-U)/L) 1006 FORMAT(/1X,NO,5X,W/L) END * * THIS SUBROUTINE IS USED TO CALCULATE * * NODES REACTION * * SUBROUTINE RSS(DLL,PP,E FA,B0,NE,U) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 B(2,6),DLL(20),E(2,1),P(6,1),BL(2,6),BLT(6,2),B0(20), & PG(6,1),EFA(20),PP(63,1),U(63),P1(6,1),B00(2,6),P2(1,6), & BLL(3),CC(6,2),BT(6,2),C(2,6),T(6,6),B00T(6,2),PGT(1,6) COMMON/A1/D(2,2)/A2/X(3),W(3) COMMON/A3/X1(20),X2(20),Y1(20),Y2(20),Z1(20),Z2(20) DO 22 I=1,3*(NE+1) 22 PP(I,1)=0.0 DO 33 N=1,NE CALL TT(T,B0,N,CCOS,SSIN) NI=3*(N-1)+1 U1=U(NI)*CCOS+U(NI+1)*SSIN W1=-U(NI)*SSIN+U(NI+1)*CCOS U2=U(NI+3) *CCOS+U(NI+4)*SSIN W2=-U(NI+3)*SSIN+U(NI+4)*CCOS O1=U(NI+2) O2=U(NI+5) X1(N)=U1 X2(N)=U2 Y1(N)=W1 Y2(N)=W2 Z1(N)=O1 Z2(N)=O2 DL=DLL(N) DO 1 I=1,2 DO 1 J=1,6 1 B(I,J)=0.0 DO 11 I=1,6 P1(I,1)=0.0 11 P(I,1)=0.0 DO 10 I=1,3 B00(1,1)=-1/DL B00(1,2)=0.0 B00(1,3)=0.0 B00(1,4)=1/DL B00(1,5)=0.0 B00(1,6)=0.0 B00(2,1)=0.0 B00(2,2)=12.*X(I)/(DL*3)-6./DL*2 B00(2,3)=-4./DL+6.*X(I)/(DL*2) B00(2,4)=0.0 B00(2,5)=6./(DL*DL)-12.*X(I)/(DL*3) B00(2,6)=-2./DL+6.*X(I)/(DL*2) G2=6.*X(I)*X(I)/DL*3-6.*X(I)/DL/dl G3=1.-4.*X(I)/dl+3.*X(I)*X(I)/DL*2 G5=6.*X(I)/(DL*DL)-6.*X(I)*X(I)/DL*3 G6=-2.*X(I)/DL+3.*X(I)*X(I)/DL*2 BLL(I)=G2*W1+G3*O1+G5*W2+G6*O2 BL(1,1)=0.0 BL(1,2)=BLL(I)*G2 BL(1,3)=BLL(I)*G3 BL(1,4)=0.0 BL(1,5)=BLL(I)*G5 BL(1,6)=BLL(I)*G6 DO 88 II=1,6 88 BL(2,II)=0.0 DO 20 K=1,2 DO 20 J=1,6 20 B(K,J)=B00(K,J)+BL(K,J) E(1,1)=-1/DL*U1+1/DL*U2+BLL(I)*BLL(I)/2 E(2,1)=B00(2,2)*W1+B00(2,3)*O1+B00(2,5)*W2+B00(2,6)*O2 CALL TTT(2,6,B,BT) CALL MTMULT(6,2,2,BT,D,CC) CALL MTMULT(6,2,1,CC,E,P1) DO 30 K=1,6 30 P(K,1)=P(K,1)+P1(K,1)*W(I)*DL 10 CONTINUE EFA(N)=-P(1,1) CALL TTT(6,1,P,P2) CALL MTMULT(1,6,6,P2,T,PGT) CALL TTT (1,6,PGT,PG) DO 40 I=1,6 II=NI+I-1 40 PP(II,1)=PP(II,1)+PG(I,1) 33 CONTINUE RETURN END * * THIS SUBROUTINE IS USED TO SOLVE * * INCREAMENT DISPLACEMENT BY AN * * CHOLESKYS METHOD * * SUBROUTINE JFC(WK,P,X,N,M,NU,UXY) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(63,64),WK(63,63),X(63),P(63,1),UXY(10,3) DO 20 I=1,N DO 20 J=1,N A(I,J)=WK(I,J) 20 A(I,M)=P(I,1) DO 30 I=1,NU KI=3*(UXY(I,1)-1)+UXY(I,2) A(KI,KI)=A(KI,KI)*10E15 A(KI,M)=UXY(I,3)*A(KI,KI) 30 CONTINUE * CALCULATION FIRST ROW OF UPPER UNIT * * TRIANGULAR MATRIX * DO 40 J=2,M 40 A(1,J)=A(1,J)/A(1,1) * CALCULATION OTHER ELEMENTS OF U AND L MATRICES * DO 100 I=2,N J=I DO 70 II=J,N SUM=0. JM1=J-1 DO 60 K=1,JM1 60 SUM=SUM+A(II,K)*A(K,J) 70 A(II,J)=A(II,J)-SUM IP1=I+1 DO 90 JJ=IP1,M SUM=0. IM1=I-1 DO 80 K=1,IM1 80 SUM=SUM+A(I,K)*A(K,JJ) 90 A(I,JJ)=(A(I,JJ)-SUM)/A(I,I) 100 CONTINUE * SOLVE FOR X(I) BY BACK SUBSTITUTION * X(N)=A(N,N+1) L=N-1 DO 120 NN=1,L SUM=0. I=N-NN IP1=I+1 DO 110 J=IP1,N 110 SUM=SUM+A(I,J)*X(J) 120 X(I)=A(I,M)-SUM RETURN END * * THIS SUBROUTINE IS USED TO FORM * * TRANSFORMANTION MATRIX T * * SUBROUTINE TT(T,BO,N,CCOS,SSIN) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION T(6,6),BO(20) CCOS=DCOS(BO(N) SSIN=DSIN(BO(N) DO 20 I=1,6 DO 20 J=1,6 20 T(I,J)=0.0 T(1,1)=CCOS T(1,2)=SSIN T(2,1)=-SSIN T(2,2)=CCOS T(3,3)=1.0 T(4,4)=CCOS T(4,5)=SSIN T(5,4)=-SSIN T(5,5)=CCOS T(6,6)=1.0 RET URN END * * THIS SUBROUTINE IS USED TO FORM LOCAL * * COORDINATE * * SUBROUTINE FLC(NE,M,RS,UI,DU,RN,XY,BO,DLL) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION RS(63,1),UI(63),DU(63),RN(20),BO(20), & XY(21,2),DLL(20) PI=3.1415926 DO 10 I=1,M RS(I,1)=0.0 UI(I)=0.0 10 DU(I)=0.0 DO 20 I=1,NE 20 RN(I)=0.0 DO 30 N=1,NE I=N J=N+1 XO=XY(J,1)-XY(I,1) YO=XY(J,2)-XY(I,2) IF(XO.GT.0.0) THEN IF(YO.GT.0.0) THEN BO(N)=DATAN(YO/XO) ELSE BO(N)=PI*2.-DABS(DATAN(YO/XO) END IF ELSE IF(YO.GT.0.0) THEN BO(N)=PI-DABS(DATAN(YO/XO) ELSE BO(N)=PI+DABS(DATAN(YO/XO) END IF END IF DLL(N)=DSQRT(YO*YO+XO*XO) 30 CONTINUE RETURN END * * THIS SUBROUTINE IS USED TO FORM THE * * ELEMENT STIFFNESS MATRIX WK AND ASSEMBLE * * THE TOTAL MATRIX TWK * * SUBROUTINE WKK(NE,DLL,RN,BO,TWK,EF,EI) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 K(6,6),K1(6,6),KW(6,6),K0(6,6),GT(6,2),B00(2,6), & BL(2,6),C(2,6),BLL(3),BLT(6,2),KK(6,6),G(2,6), & K3(6,6),K4(6,6),RN(20),DLL(20),BO(20),T(6,6),CC(6,2), & TWK(63,63),WT(6,6),K2(6,6),B00T(6,2),T0(6,6) COMMON/A1/D(2,2)/A2/X(3),W(3) COMMON/A3/X1(20),X2(20),Y1(20),Y2(20),Z1(20),Z2(20) NI=3*(NE+1) DO 91 I=1,NI DO 91 J=1,NI 91 TWK(I,J)=0.0 DO 200 N=1,NE DL=DLL(N) U1=X1(N) U2=X2(N) W1=Y1(N) W2=Y2(N) O1=Z1(N) O2=Z2(N) X(1)=0.1127*DL X(2)=0.5*DL X (3)=0.8873*DL W(1)=0.555556/2. W(2)=0.8888889/2. W(3)=0.555556/2. DO 211 I=1,6 DO 211 J=1,6 K0(I,J)=0.0 K1(I,J)=0.0 K2(I,J)=0.0 K3(I,J)=0.0 K4(I,J)=0.0 K(I,J)=0.0 KW(I,J)=0.0 211 CONTINUE D(1,1)=EF D(1,2)=0.0 D(2,1)=0.0 D(2,2)=EI DO 10 I=1,3 B00(1,1)=-1/DL B00(1,2)=0.0 B00(1,3)=0.0 B00(1,4)=1/DL B00(1,5)=0.0 B00(1,6)=0.0 B00(2,1)=0.0 B00(2,2)=12.*X(I)/(DL*3)-6./DL*2 B00(2,3)=-4./DL+6.*X(I)/(DL*2) B00(2,4)=0.0 B00(2,5)=6./(DL*DL)-12.*X(I)/(DL*3) B00(2,6)=-2./DL+6.*X(I)/(DL*2) G(1,1)=0.0 G(1,2)=6.*X(I)*X(I)/DL*3-6.*X(I)/DL/dl G(1,3)=1.-4.*x(i)/dl+3.*x(i)*x(i)/dl*2 G(1,4)=0.0 G(1,5)=6.*X(I)/(DL*DL)-6.*X(I)*X(I)/DL*3 G(1,6)=-2.*X(I)/DL+3.*X(I)*X(I)/DL*2 DO 99 II=1,6 99 G(2,II)=0.0 BLL(I)=G(1,2)*W1+G(1,3)*O1+G(1,5)*W2+G(1,6)*O2 BL(1,1)=0.0 BL(1,2)=BLL(I)*G(1,2) BL(1,3)=BLL(I)*G(1,3) BL(1,4)=0.0 BL(1,5)=BLL(I)*G(1,5) BL(1,6)=BLL(I)*G(1,6) DO 88 II=1,6 BL(2,II)=0.0 88 CONTINUE CALL TTT(2,6,B00,B00T) CALL MTMULT(6,2,2,B00T,D,CC) CALL MTMULT(6,2,6,CC,B00,KW) DO 20 M=1,6 DO 20 J=1,6 20 K0(M,J)=K0(M,J)+KW(M,J)*W