PROGRAM BUCKLE1C C====================================================================== C AN FEM SOLVER FOR BUCKLING PROBLEM WITH MOMENT AT BOTH ENDS OF BEAM C EQUATION: D2U/DXDX + ALPHA*U=0; ALPHA = P/(EI) C P=APPLIED FORCE AT ENDS OF BEAM, E=YOUNG MODULUS, I=2ND MOMENT OF C INERTIA. RL=LENGTH OF ELEMENT, IBTYPE(1)=BOUNDARY CONDITION AT THE C LEFT END OF BEAM, IBTYPE(2)=BC AT RIGHT END. IBTYPE(I)=1 FOR FIXED C DISPLACEMENT, IBTYPE(I)=2 FOR PRESCRIBED SLOPE AT THE END. C ****** SYMMETRIC TRI-DIAGONAL MATRIX SOLVER****************** C JUNE 16, 1999 EIJI FUKUMORI C DOWN AND UPPER STREAM SIDE ELEMENTS: NEW-TWO-NODE-PARABOLIC. C====================================================================== IMPLICIT REAL*8 ( A-H , O-Z ) PARAMETER ( ND=2, MXE=100, MXN=MXE+1,NBW=ND ) DIMENSION NODEX(MXE,ND),EI(MXE),X(MXN),A(MXN,NBW),RHS(MXN), * IBTYPE(2), BV(2) C====================================================================== C (1) READING OF DATA CALL INPUT ( MXE,MXN,ND,P,NE,NNODE,NODEX,EI,X,IBTYPE,BV ) C====================================================================== C (2) CONSTRUCTION OF FEM-MATRIX EQUATION CALL MATRIX ( MXE,MXN,ND,NBW,P,NE,NNODE,NODEX,EI,X,A,RHS,IBTYPE ) C====================================================================== C (3) IMPLEMENTATION OF BOUNDARY CONDITION CALL FORM ( MXN, NBW, NNODE, A, RHS, IBTYPE, BV ) C====================================================================== C (4) SOLVE FOR UNKNOWN VARIABLES CALL SYSTEM ( MXN, NBW, NNODE, A, RHS ) C====================================================================== C (5) PRINTING RESULTS OPEN ( 1,FILE='SOLUTION.FEM', STATUS='UNKNOWN') WRITE (*,*) 'SOLUTION IS IN A FILE OF SOLUTION.FEM.' DO I = 1 , NNODE WRITE(1,100) I, X(I), RHS(I) 100 FORMAT ( ' NODAL #=', I3, ' X=',G16.7, 'DISPLACEMENT=',G20.11 ) END DO CLOSE (1) STOP' NORMAL TERMINATION' END C C SUBROUTINE INPUT ( MXE,MXN,ND,P,NE,NNODE,NODEX,EI,X,IBTYPE,BV ) IMPLICIT REAL*8 ( A-H , O-Z ) DIMENSION NODEX(MXE,ND),EI(MXE),X(MXN),IBTYPE(2), BV(2) OPEN ( 1,FILE='BUCKLE.DAT', STATUS='OLD') READ(1,*) P READ(1,*) NE DO I = 1 , NE READ(1,*) IEL, (NODEX(IEL,J),J=1,ND),EI(IEL) END DO NNODE = NE + 1 DO I = 1 , NNODE READ(1,*) NODE, X(NODE) END DO READ(1,*) IBTYPE(1), BV(1) READ(1,*) IBTYPE(2), BV(2) CLOSE (1) RETURN END C C SUBROUTINE MATRIX (MXE,MXN,ND,NBW,P,NE,NNODE,NODEX,EI,X,A,RHS, * IBTYPE) IMPLICIT REAL*8 ( A-H , O-Z ) DIMENSION NODEX(MXE,ND),EI(MXE),X(MXN),A(MXN,NBW),RHS(MXN), * IBTYPE(2) DO I = 1 , NNODE A(I,1) = 0. A(I,2) = 0. RHS(I) = 0. END DO DO IEL = 1 , NE I = NODEX(IEL,1) J = NODEX(IEL,2) RL = X(J) - X(I) ALPHA = P / EI(IEL) C-------- TOP ELEMENT IF ( IEL .EQ. 1 ) THEN IF ( IBTYPE(1) .EQ. 1 ) THEN A(I,1) = A(I,1) - 1./RL + ALPHA*RL/3. A(I,2) = A(I,2) + 1./RL + ALPHA*RL/6. A(J,1) = A(J,1) - 1./RL + ALPHA*RL/3. END IF IF ( IBTYPE(1) .EQ. 2 ) THEN A(I,1) = A(I,1) - 4./(3.*RL) + ALPHA*RL*8./15. A(I,2) = A(I,2) + 4./(3.*RL) + ALPHA*RL*2./15. A(J,1) = A(J,1) - 4./(3.*RL) + ALPHA*RL*3./15. END IF END IF C------- LAST ELEMENT IF ( IEL .EQ. NE ) THEN IF ( IBTYPE(2) .EQ. 1 ) THEN A(I,1) = A(I,1) - 1./RL + ALPHA*RL/3. A(I,2) = A(I,2) + 1./RL + ALPHA*RL/6. A(J,1) = A(J,1) - 1./RL + ALPHA*RL/3. END IF IF ( IBTYPE(2) .EQ. 2 ) THEN A(I,1) = A(I,1) - 4./(3.*RL) + ALPHA*RL*3./15. A(I,2) = A(I,2) + 4./(3.*RL) + ALPHA*RL*2./15. A(J,1) = A(J,1) - 4./(3.*RL) + ALPHA*RL*8./15. END IF END IF C-------- IN BETWEEN IF ( (IEL .GT. 1) .AND. (IEL .LT. NE) ) THEN A(I,1) = A(I,1) - 1./RL + ALPHA*RL/3. A(I,2) = A(I,2) + 1./RL + ALPHA*RL/6. A(J,1) = A(J,1) - 1./RL + ALPHA*RL/3. END IF C END DO RETURN END C C SUBROUTINE FORM ( MXN, NBW, NNODE, A, RHS, IBTYPE, BV ) IMPLICIT REAL*8 ( A-H , O-Z ) DIMENSION IBTYPE(2),BV(2),A(MXN,NBW), RHS(MXN) IF ( IBTYPE(1) .EQ. 1 ) THEN A(1,1) = 1. RHS(1) = BV(1) RHS(2) = RHS(2) - BV(1)*A(1,2) A(1,2) = 0. ELSE RHS(1) = RHS(1) - BV(1) END IF IF ( IBTYPE(2) .EQ. 1 ) THEN A(NNODE,1) = 1. RHS(NNODE) = BV(2) RHS(NNODE-1) = RHS(NNODE-1) - BV(2)*A(NNODE-1,2) A(NNODE-1,2) = 0. ELSE RHS(NNODE) = RHS(NNODE) - BV(2) END IF RETURN END C C SUBROUTINE SYSTEM ( MXN, NBW, NNODE, A , B ) IMPLICIT REAL*8 ( A-H , O-Z ) DIMENSION A(MXN,NBW) , B(MXN) B(1) = B(1) / A(1,1) A(1,1) = A(1,2) / A(1,1) DO I = 2 , NNODE P = A(I,1) - A(I-1,2) * A(I-1,1) A(I,1) = A(I,2) / P B(I) = ( B(I) - A(I-1,2)*B(I-1) ) / P END DO C------ BACK SUBSTITUTION ---- DO I = NNODE-1, 1,-1 B(I) = B(I) - A(I,1) * B(I+1) END DO RETURN END