C PROGRAM PLFNL2 BY MIKE FOLSE LAST MODIFIED 12/19/01 C ADDED DOUBLE PRECISION DUE TO NUMERICAL STABILITY PROBLEMS C STATIC ANALYSIS OF PLANE FRAMES AND TRUSSES, C INCLUDING P-DELTA EFFECTS C (NOT VERY LARGE DISPLACEMENTS) C USE CONSISTENT UNITS C C REVISED 1989 TO INCLUDE P-DELTA EFFECTS ON FIXED END C ACTIONS C C INPUT FROM FILE PLFNL2.DAT OR FROM KEYBOARD C YOU MAY BUILD PLFNL2.DAT IN AN EDITOR OR WITHIN PROGRAM C PROGRAM AUTOMATICALLY STORES KEYBOARD INPUT IN PLFNL2.DAT C OUTPUT IS TO FILE PLFNL2.OUT C C INPUT LIST FOR FILE PLFNL2.DAT(FREE FORMAT-SEPERATE INPUT C ITEMS WITH ONE COMMA OR ONE OR MORE BLANKS,INPUT C ALL ITEMS EVEN IF ZERO) C C LINE 1: TITLE(15A4) C LINE 2: NUMBER OF JOINTS(NJT),NUMBER OF MEMBERS(NMEM), C NUMBER OF JOINTS WITH APPLIED LOADS(NLDJT), C MODULUS OF ELASTICITY(E), C NUMBER OF PIN END MEMBERS(NPE) C NPE=-1 INDICATES ALL MEMBERS PINNED IE. A TRUSS C NPD = P-DELTA CODE =0 NO P-DELTA C =1 P-DELTA C C LINE 3 (REPEAT NJT TIMES): JOINT NUMBER, X COORDINATE, C Y COORDINATE, SUPPORT CODE(NS) C (=1 IF AT LEAST 1 DOF AT JOINT IS SUPPORT,0 OTHERWISE) C NOTE:INPUT JOINTS IN ASCENDING ORDER,MISSING JOINTS WILL C BE GENERATED BASED ON EQUAL SPACING C LINE 3A IF NS=1,INPUT SPRING SUPPORT CONSTANTS IN X,Y, AND C Z ROTATION. C IF NEGATIVE NUMBER INPUT,DEFAULT = 1.E30 ASSIGNED C LINE 4 (REPEAT NMEM TIMES): MEMBER NUMBER, JOINT NUMBER OF C JOINT TO BE CONSIDERED J NODE THIS MEMBER, C JOINT NUMBER OF K NODE THIS MEMBER, C CROSS-SECTION AREA THIS MEMBER, C MOMENT OF INERTIA ABOUT BENDING AXIS, C MEMBER LOAD CODE(MLC) C (=1 IF MEMBER HAS LOAD, =0 OTHERWISE) C LINE 4A IF MLC=1 INPUT VALUE OF DISTRIBUTED LOAD AT J NODE C THEN VALUE AT K NODE, C THEN NUMBER OF CONCENTRATED LOADS (3 MAX) C INPUT POSITIVE LOAD IN DOWN DIRECTION C LINE 4B IF THERE ARE CONC LOADS,INPUT C ( LOAD VALUE , DISTANCE FROM J) FOR EACH LOAD C LINE 5 (REPEAT NLDJT TIMES): JOINT NUMBER OF LOADED JOINT, C APPLIED LOAD IN X DIRECTION, C APPLIED LOAD IN Y DIRECTION, C APPLIED COUNTERCLOCKWISE MOMENT ABOUT Z AXIS C LINE 6 (IF NPE .NE. 0) INPUT MEMBER NUMBERS OF C PIN ENDED MEMBERS IMPLICIT REAL*8 (A-H, O-Z) DIMENSION X(98),Y(98),JEND(130),KEND(130),RI(130),D(6), 1 SM(6,6),RT(6,6),RTRANS(6,6),SX(98),SY(98),SRZ(98), 1 B(294),TITLE(15),AX(130),ALEN(130),SMR(6,6), 1IN(6),SMRT(6,6),CXX(130),CYY(130),EACT(6),NS(98),MLC(130), 1 WL(130),WR(130),FEM(6),P(6),NP(130),NPP(130),NCLD(130), 1 PV(130,3),PD(130,3),AF(130),AFP(130),FEMS(130,6) COMMON/AA/A(294,63) DIMENSION AJ(294) OPEN(6,FILE='PLFNL2.OUT',STATUS='UNKNOWN') C PROGRAM DIMENSIONED FOR MAXIMUM 98 JOINTS AND 130 MEMBERS IDIM=294 JDIM=63 30 CONTINUE NCD=0 IB=0 C MAXIMUM DIFFERENCE IN END NODE NUMBERS=20 INP=5 IOUT=6 IP=7 NGMP=0 NGM=0 WRITE(*,5000) 5000 FORMAT(5X,'PROGRAM PLFNL2 BY MIKE FOLSE, P-DELTA ANALYSIS', 1' PLANE FRAMES',/,5X,'SELECT INPUT MODE', 1 //,15X,'(1) KEYBOARD', 1//,15X,'(0) FILE PLFNL2.DAT',//,15X,'(2) TERMINATE',/) READ(*,*)NCD IF(NCD.NE.0.AND.NCD.NE.1)STOP WRITE(*,5005) 5005 FORMAT(/,5X,'OUTPUT IS TO FILE PLFNL2.OUT',//,5X, 1 ' "EDIT PLFNL2.OUT" TO VIEW CONTENTS',//,5X, 1 ' "PRINT PLFNL2.OUT" TO GET PRINTED OUTPUT') IF(NCD.EQ.1)GO TO 50 OPEN(5,FILE='PLFNL2.DAT',STATUS='OLD',ACCESS='SEQUENTIAL') READ(INP,1000)TITLE GO TO 60 50 CONTINUE OPEN(7,FILE='PLFNL2.DAT',STATUS='UNKNOWN') WRITE(*,5010) 5010 FORMAT(/,5X,'ENTER TITLE UP TO 60 CHARACTERS',/) READ(*,1000)TITLE WRITE(IP,1001)TITLE 1001 FORMAT(1X,15A4) 1000 FORMAT(15A4) 60 CONTINUE WRITE(IOUT,2000)TITLE 2000 FORMAT(//,10X,15A4) IF(NCD.EQ.1)GO TO 80 READ(INP,*)NJT,NMEM,NLDJT,E,NPE,NPD GO TO 85 80 CONTINUE WRITE(*,5020) 5020 FORMAT(/,5X,'ENTER NUMBER OF JOINTS(98 MAXIMUM)',/) READ(*,*)NJT WRITE(*,5021) 5021 FORMAT(/,5X,'ENTER NUMBER OF MEMBERS(130 MAXIMUM)',/) READ(*,*)NMEM WRITE(*,5022) 5022 FORMAT(/,5X,'ENTER NUMBER OF JOINTS WITH APPLIED LOADING',/) READ(*,*)NLDJT WRITE(*,5023) 5023 FORMAT(/,5X,'ENTER MODULUS OF ELASTICITY,E (USE CONSISTENT', 1 ' UNITS)',/) READ(*,*)E WRITE(*,5025) 5025 FORMAT(/,5X,'ENTER NUMBER OF PIN ENDED MEMBERS', 1 /,10X,'VALUE OF -1 INDICATES ALL PINNED(A TRUSS)',/) READ(*,*)NPE WRITE(*,5078) 5078 FORMAT(/,5X,'DO YOU WANT TO CALCULATE P-DELTA EFFECTS?', 1' (0=NO,1=YES)',/) READ(*,*)NPD WRITE(IP,1010)NJT,NMEM,NLDJT,E,NPE,NPD 1010 FORMAT(3I5,E12.3,2I5) WRITE(*,5024) 5024 FORMAT(///,5X,'INPUT JOINT COORDINATES AND SPRING ', 1 'SUPPORT CONSTANTS',//) 85 CONTINUE WRITE(IOUT,2010)NJT,NMEM,NLDJT,E,NPE,NPD 2010 FORMAT(//,10X,'NUMBER OF JOINTS =',I5, 1 /,10X,'NUMBER OF MEMBERS =',I5, 1 /,10X,'NUMBER OF JOINTS WITH APPLIED LOAD(S) =',I5, 1/,10X,'MODULUS OF ELASTICITY =',E12.3,/,10X, 1'NUMBER OF PIN END MEMBERS =',I5,' ( -1 INDICATES TRUSS)',/ 110X,'NPD=',I5,'(1 MEANS CALCULATE P-DELTA EFFECTS)' 1 //,40X,'SPRING CONSTANTS',/,1X,'JOINT X Y', 1 ' SX SY SRZ',/) NUMEQ=3*NJT DO 87 I=1,NJT SX(I)=0. SY(I)=0. SRZ(I)=0. 87 CONTINUE MM=0 88 CONTINUE IF(NCD.EQ.1)GO TO 90 READ(INP,*)M,X(M),Y(M),NS(M) IF(NS(M).EQ.1)READ(INP,*)SX(M),SY(M),SRZ(M) GO TO 100 90 CONTINUE WRITE(*,5030) 5030 FORMAT(/,2X,'ENTER JOINT NUMBER,X COORD,', 1'Y COORD AND SUPPORT CODE (=1 IF ANY SUPPORTS)',/) READ(*,*)M,X(M),Y(M),NS(M) IF(NS(M).EQ.1)WRITE(*,5031) 5031 FORMAT(/,5X,'ENTER SPRING SUPPORT CONSTANTS X,Y,Z',/, 19X,'(ENTER -1 FOR RIGID SUPPORT,O IF DOF NOT SUPPORTED)',/) IF(NS(M).EQ.1)READ(*,*)SX(M),SY(M),SRZ(M) WRITE(IP,1020)M,X(M),Y(M),NS(M) 1020 FORMAT(I5,2E13.4,I5) IF(NS(M).EQ.1)WRITE(IP,1021)SX(M),SY(M),SRZ(M) 1021 FORMAT(6E12.3) 100 CONTINUE IF(SX(M).LT.0.)SX(M)=1.E30 IF(SY(M).LT.0.)SY(M)=1.E30 IF(SRZ(M).LT.0.)SRZ(M)=1.E30 IF(M.EQ.(MM+1))GO TO 116 C GENERATE MISSING NODES M1=MM+1 M2=M-1 DIFF=FLOAT(M-MM) XINC=(X(M)-X(MM))/DIFF YINC=(Y(M)-Y(MM))/DIFF AINC=0. DO 115 K=M1,M2 AINC=AINC+1. X(K)=X(MM)+AINC*XINC Y(K)=Y(MM)+AINC*YINC C NOTE:GENERATED NODES CANNOT BE SUPPORTS NS(K)=0 115 CONTINUE 116 IF(M.EQ.NJT)GO TO 119 MM=M GO TO 88 119 CONTINUE DO 120 I=1,NJT 120 WRITE(IOUT,2020)I,X(I),Y(I),SX(I),SY(I),SRZ(I) 2020 FORMAT(1X,I5,2E13.4,3E12.3) DO 140 I=1,NMEM NP(I)=0 IF(NCD.EQ.1)GO TO 130 READ(INP,*)M,JEND(M),KEND(M),AX(M),RI(M),MLC(M) IF(MLC(M).EQ.0)GO TO 135 READ(INP,*)WL(M),WR(M),NCLD(M) NCL=NCLD(M) IF(NCLD(M).GT.0)READ(INP,*)(PV(M,J),PD(M,J),J=1,NCL) GO TO 135 130 CONTINUE WRITE(*,5040) 5040 FORMAT(/,5X,'ENTER MEM NUM,J NODE,K NODE,AREA,INERTIA', 1 ' AND LOAD CODE', 1 /,6X,'(LOAD CODE =1 IF ANY MEMBER LOAD, 0 OTHERWISE)',/) READ(*,*)M,JEND(M),KEND(M),AX(M),RI(M),MLC(M) WRITE(IP,1030)M,JEND(M),KEND(M),AX(M),RI(M),MLC(M) 1030 FORMAT(3I5,2E12.3,I5) IF(MLC(M).EQ.0)GO TO 135 WRITE(*,5042) 5042 FORMAT(/,2X,'ENTER DISTRIBUTED LOAD VALUE J END AND K END', 1 /,5X,'(-Y LOCAL AXIS IS POSITIVE)',/,5X, 1 'AND NUMBER OF CONCENTRATED LOADS(3 MAX)',/) READ(*,*)WL(M),WR(M),NCLD(M) WRITE(IP,6021)WL(M),WR(M),NCLD(M) 6021 FORMAT(2E12.3,I5) NCL=NCLD(M) IF(NCL.EQ.0)GO TO 135 WRITE(*,5043) 5043 FORMAT(/,2X,'ENTER (LOAD VALUE,DISTANCE FROM J NODE)FOR', 1 ' UP TO 3 CONC LOADS') READ(*,*)(PV(M,J),PD(M,J),J=1,NCL) WRITE(IP,1021)(PV(M,J),PD(M,J),J=1,NCL) 135 CONTINUE C C CALCULATE BAND WIDTH OF STIFFNESS MATRIX NBW=3*(IABS(JEND(M)-KEND(M))+1) IF(NBW.GT.IB)IB=NBW 140 CONTINUE IF(IB.LE.JDIM)GO TO 150 WRITE(IOUT,2025)IB 2025 FORMAT(/,10X,'BAND WIDTH EXCEEDS DIMENSION ALLOWED', 1 'IB=',I5) 150 CONTINUE WRITE(IOUT,2030) 2030 FORMAT(//,5X,'MEMBER J NODE K NODE AREA', 17X, 'INERTIA') DO 160 I=1,NMEM 160 WRITE(IOUT,2040)I,JEND(I),KEND(I),AX(I),RI(I) 2040 FORMAT(3(5X,I5),4E12.3) WRITE(IOUT,2035) 2035 FORMAT(//,5X,'MEMBER LOADS',/,2X,'MEMBER W AT J', 1' W AT K P1 D1 P2 D2 P3', 1' D3') DO 170 I=1,NMEM IF(MLC(I).EQ.0)GO TO 170 NCL=NCLD(I) WRITE(IOUT,2041)I,WL(I),WR(I),(PV(I,J),PD(I,J),J=1,NCL) 170 CONTINUE 2041 FORMAT(1X,I5,2X,8E9.2) IF(IB.GT.JDIM)STOP C INITIALIZE LOAD MATRIX TO ZERO DO 180 I=1,NUMEQ 180 B(I)=0. IF(NLDJT.EQ.0)GO TO 192 C INPUT JOINT LOADS WRITE(IOUT,2045) 2045 FORMAT(//,5X,'JOINT XFORCE YFORCE ZMOM',/) DO 190 I=1,NLDJT IF(NCD.EQ.1)GO TO 185 READ(INP,*)JT,XFORCE,YFORCE,ZMOM GO TO 188 185 CONTINUE WRITE(*,5050) 5050 FORMAT(/,5X,'ENTER JOINT,FORCE X DIRECTION,', 1 'FORCE Y DIRECTION, AND CC MOMENT',/) READ(*,*)JT,XFORCE,YFORCE,ZMOM WRITE(IP,1040)JT,XFORCE,YFORCE,ZMOM 1040 FORMAT(I5,3E12.3) 188 CONTINUE WRITE(IOUT,2050)JT,XFORCE,YFORCE,ZMOM 2050 FORMAT(5X,I5,3E12.3) C ADD FORCES TO APPLIED LOAD MATRIX NXDF=3*JT-2 NYDF=NXDF+1 NZDF=NXDF+2 B(NXDF)=XFORCE B(NYDF)=YFORCE B(NZDF)=ZMOM 190 CONTINUE DO 191 I =1,NUMEQ 191 AJ(I)=DBLE(B(I)) 192 CONTINUE IF(NPE.LE.0)GO TO 209 IF(NCD.EQ.1)GO TO 200 READ(INP,*)(NPP(I),I=1,NPE) GO TO 205 200 CONTINUE WRITE(*,5093) 5093 FORMAT(/,5X,'ENTER MEMBER NUMBERS OF PIN ENDED MEMBERS',/) READ(*,*)(NPP(I),I=1,NPE) WRITE(IP,5055)(NPP(I),I=1,NPE) 5055 FORMAT(10I5) 205 CONTINUE WRITE(IOUT,5056)(NPP(I),I=1,NPE) 5056 FORMAT(//,5X,'THE FOLLOWING MEMBERS ARE PIN ENDED', 1 (/,3X,16I5)) DO 207 I=1,NPE NP(NPP(I))=1 207 CONTINUE 209 CONTINUE DO 208 I=1,NMEM IF(NPE.EQ.-1)NP(I)=1 JJ=JEND(I) KK=KEND(I) XC=X(KK)-X(JJ) YC=Y(KK)-Y(JJ) AL=SQRT(XC*XC+YC*YC) CX=XC/AL CY=YC/AL CXX(I)=CX CYY(I)=CY C CX,CY, ARE DIRECTION COSINES FOR MEMBER I, I.E. COSINES OF C ANGLES BETWEEN MEMBER AXIS AND GLOBAL X,Y AXES ALEN(I)=AL C ADD MEMBER END ACTIONS TO LOAD MATRIX (FOR 0 AXIAL LOAD) IF(MLC(I).EQ.0)GO TO 208 CALL FORMRT(RT,RTRANS,CX,CY) CALL INDX(JJ,KK,IN) CALL FM(WL(I),WR(I),AL,FEM,NP(I)) NCL=NCLD(I) IF(NCL.EQ.0)GO TO 265 DO 263 J=1,NCL 263 CALL FM2(PV(I,J),PD(I,J),AL,FEM,NP(I)) 265 CALL MULT(RTRANS,FEM,1,P) C STORE FIXED END ACTIONS DO 267 J = 1,6 267 FEMS(I,J)=FEM(J) DO 270 J=1,6 AJ(IN(J))=AJ(IN(J))-P(J) 270 CONTINUE 208 CONTINUE 210 CONTINUE DO 220 I=1,NUMEQ C REINITIALIZE GLOBAL STIFFNESS MATRIX TO ZERO DO 220 J=1,IB A(I,J)=0. 220 CONTINUE C FORM MEMBER STIFFNESS MATRICES, ROTATE TO GLOBAL COORDINATES, C INDEX TO GLOBAL DEGREE OF FREEDOM NUMBERS, AND ADD TO C GLOBAL STIFFNESS MATRIX DO 300 I=1,NMEM CX=CXX(I) CY=CYY(I) AL=ALEN(I) CALL FORMRT(RT,RTRANS,CX,CY) CALL FORMSM(SM,AL,AX(I),E,RI(I),NP(I)) IF(NGM.GT.0.AND.AX(I).GT.0.)CALL GSM(SM,AL,NP(I),AF(I)) CALL MULT(SM,RT,6,SMRT) CALL MULT(RTRANS,SMRT,6,SMR) CALL INDX(JEND(I),KEND(I),IN) C FORM UPPER HALF OF STIFFNESS MATRIX,STORE IN RECTANGULAR C ARRAY WITH DIAGONAL ELEMENT IN COLUMN 1 DO 260 J=1,6 JG=IN(J) DO 260 K=1,6 KG=IN(K) IF(KG.LT.JG)GO TO 260 LG=KG-JG+1 A(JG,LG)=A(JG,LG)+SMR(J,K) 260 CONTINUE 300 CONTINUE C CORRECT GLOBAL STIFFNESS MATRIX FOR BOUNDARY CONDITIONS C I.E. SUPPORT SPRINGS WHERE INPUT DO 360 I=1,NJT NXDF=3*I-2 NYDF=NXDF+1 NZDF=NXDF+2 IF(NS(I).EQ.0)GO TO 350 A(NXDF,1)=A(NXDF,1)+SX(I) A(NYDF,1)=A(NYDF,1)+SY(I) A(NZDF,1)=A(NZDF,1)+SRZ(I) 350 IF(NPE.EQ.-1)A(NZDF,1)=1.E30 360 CONTINUE C SOLVE AX=B FOR DISPLACEMENT MATRIX X CALL UDU(AJ,NUMEQ,IB,IDIM,JDIM,NER) IF(NER.EQ.0)GO TO 370 WRITE(IOUT,4200) 4200 FORMAT(/,10X,'STIFFNESS MATRIX IS SINGULAR') GO TO 550 370 CONTINUE C NOTE: DEFLECTIONS ARE RETURNED IN B MATRIX WHICH WENT INTO C UDU AS APPLIED LOADS IF(NGM.GT.0)GO TO 420 WRITE(6,4203) 4203 FORMAT(///,5X,'OUTPUT NOT INCLUDING P-DELTA EFFECTS',//) 375 WRITE(IOUT,4100) 4100 FORMAT(//,9X,'JOINT DISPLACEMENTS AND REACTIONS',//,5X,'JT ' 1,' X DISPL Y DISPL Z ROT X REACT Y REACT', 1' ZM REACT') DO 400 I=1,NJT NXDF=3*I-2 NYDF=NXDF+1 NZDF=NXDF+2 XREACT=-SX(I)*AJ(NXDF) YREACT=-SY(I)*AJ(NYDF) ZREACT=-SRZ(I)*AJ(NZDF) WRITE(IOUT,4110)I,AJ(NXDF),AJ(NYDF),AJ(NZDF), 1 XREACT,YREACT,ZREACT 4110 FORMAT(1X,I5,3D11.3,3E11.3) 400 CONTINUE WRITE(IOUT,4115) 4115 FORMAT(//,13X,'MEMBER END ACTIONS',//,1X,'MEMBER J ', 1'K AX J SHR J MOM J AX K SHR K', 1 6X,'MOM K') 420 CONTINUE DO 500 I=1,NMEM AL=ALEN(I) JJ=JEND(I) KK=KEND(I) CALL INDX(JJ,KK,IN) DO 450 J=1,6 450 D(J)=AJ(IN(J)) CALL FORMRT(RT,RTRANS,CXX(I),CYY(I)) CALL FORMSM(SM,AL,AX(I),E,RI(I),NP(I)) IF(AX(I).LE.0.)GO TO 455 IF(NGM.GT.0.OR.NGMP.EQ.1)CALL GSM(SM,AL,NP(I),AF(I)) 455 CONTINUE CALL MULT(SM,RT,6,SMRT) CALL MULT(SMRT,D,1,EACT) AFP(I)=EACT(4) IF(NGM.GT.0)GO TO 500 C ADD FIXED END ACTIONS IF(MLC(I).EQ.0)GO TO 480 DO 460 K=1,6 460 EACT(K)=EACT(K)+FEMS(I,K) 480 CONTINUE WRITE(IOUT,4120)I,JJ,KK,(EACT(J),J=1,6) 4120 FORMAT(1X,I3,2I4,6E11.3) 500 CONTINUE 550 CONTINUE IF(NPD.EQ.0)STOP NGM=NGM+1 IF(NGM.LT.5)GO TO 555 WRITE(IOUT,*)' DOES NOT CONVERGE' STOP 555 CONTINUE IF(NGM.EQ.1)GO TO 570 C COMPARE NEW AND OLD VALUES OF MEMBER AXIAL FORCE DIFMAX=0. DO 560 I=1,NMEM PCR=9.87*E*RI(I)/ALEN(I)**2 RAT=AFP(I)/PCR IF(RAT.GT.-.02)GO TO 560 DIF=ABS((AFP(I)-AF(I))/AFP(I)) IF(DIF.GT.DIFMAX)DIFMAX=DIF 560 CONTINUE IF(DIFMAX.LT.0.01)GO TO 590 570 CONTINUE DO 580 I=1,NMEM 580 AF(I)=AFP(I) C C RECREATE LOAD MATRIX INCLUDING EFFECTS OF AXIAL LOAD ON C FIXED END ACTIONS DO 582 I=1,NUMEQ C INITIALIZE WITH JOINT LOADS AJ(I)=B(I) 582 CONTINUE C ADD MEMBER FIXED END ACTIONS C DO 585 I = 1,NMEM IF(MLC(I).EQ.0)GO TO 585 C IF AXIAL STRESS IS LESS THAN .1 PCR COMPRESSION, IGNORE C EFFECTS OF AXIAL FORCE ON FIXED END ACTIONS AL=ALEN(I) EI=E*RI(I) PCR=9.87*EI/AL**2 RAT=AFP(I)/PCR CX=CXX(I) CY=CYY(I) JJ=JEND(I) KK=KEND(I) CALL FORMRT(RT,RTRANS,CX,CY) CALL INDX(JJ,KK,IN) NCL=NCLD(I) IF(RAT.LT.-0.1.AND.AX(I).GT.0.)GO TO 588 CALL FM(WL(I),WR(I),AL,FEM,NP(I)) IF(NCL.EQ.0)GO TO 665 DO 663 J=1,NCL 663 CALL FM2(PV(I,J),PD(I,J),AL,FEM,NP(I)) 665 CONTINUE GO TO 680 588 CONTINUE C INCLUDE AXIAL EFFECTS OF FIXED END ACTIONS CALL FMN(WL(I),WR(I),AL,FEM,NP(I),AFP(I),EI) IF(NCL.EQ.0)GO TO 675 DO 670 J=1,NCL 670 CALL FM2N(PV(I,J),PD(I,J),AL,FEM,NP(I),AFP(I),EI) 675 CONTINUE C 680 CALL MULT(RTRANS,FEM,1,P) C SAVE FIXED END ACTIONS DO 690 J=1,6 690 FEMS(I,J)=FEM(J) C DO 700 J=1,6 700 AJ(IN(J))=AJ(IN(J))-P(J) 585 CONTINUE 586 CONTINUE GO TO 210 590 CONTINUE NPD=0 NGMP=1 WRITE(6,5001)NGM 5001 FORMAT(///,10X,'OUTPUT INCLUDING P-DELTA EFFECTS', 1 /,5X,'NUMBER OF CYCLES=',I5,//) NGM=0 GO TO 375 END C C C SUBROUTINE UDU SUBROUTINE UDU(X,NUMEQ,IB,IDIM,JDIM,NER) C UDU FOR BANDED SYMMETRIC OR MATRICES C BY ERIK THOMPSON OF CSU IMPLICIT REAL*8 (A-H, O-Z) COMMON/AA/A(294,63) DIMENSION X(IDIM) NER=0 NEM1=NUMEQ-1 DO 450 I=1,NEM1 JEND=NUMEQ-I+1 IF(JEND.GT.IB)JEND=IB DO 440 J=2,JEND J1=J+I-1 IF(A(I,1).EQ.0.)GO TO 480 FAC=A(I,J)/A(I,1) K1=0 DO 430 K=J,JEND K1=K1+1 A(J1,K1)=A(J1,K1)-FAC*A(I,K) 430 CONTINUE A(I,J)=FAC 435 CONTINUE X(J1)=X(J1)-A(I,J)*X(I) 440 CONTINUE 450 CONTINUE X(NUMEQ)=X(NUMEQ)/A(NUMEQ,1) DO 470 I=1,NEM1 I1=NUMEQ-I X(I1)=X(I1)/A(I1,1) JEND=NUMEQ-I1+1 IF(JEND.GT.IB)JEND=IB DO 460 J=2,JEND J1=I1+J-1 X(I1)=X(I1)-A(I1,J)*X(J1) 460 CONTINUE 470 CONTINUE RETURN 480 NER=1 RETURN END SUBROUTINE FORMSM(SM,AL,AX,E,RI,NP) C STIFFNESS MATRIX PLANE FRAME OR TRUSS IMPLICIT REAL*8 (A-H, O-Z) DIMENSION SM(6,6) DO 10 I = 1,6 DO 10 N = I,6 SM(I,N)=0. 10 SM(N,I)=0.0 SM(1,1)=ABS(AX)*E/AL SM(4,4)=SM(1,1) SM(1,4)=-SM(1,1) SM(4,1)=SM(1,4) SM(3,3)=0.000001*E*RI/AL SM(6,6)=SM(3,3) IF(NP.EQ.1)RETURN SM(2,2)=12.*E*RI/AL**3 SM(5,5)=SM(2,2) SM(2,5)=-SM(2,2) SM(2,3)=6.*E*RI/AL**2 SM(2,6)=SM(2,3) SM(3,5)=-SM(2,3) SM(5,6)=SM(3,5) SM(3,3)=4.*E*RI/AL SM(6,6)=SM(3,3) SM(3,6)=.5*SM(3,3) DO 20 I=1,5 J=I+1 DO 20 M=J,6 20 SM(M,I)=SM(I,M) RETURN END SUBROUTINE FORMRT(RT,RTRANS,CX,CY) C FORM ROTATION MATRICES IMPLICIT REAL*8 (A-H,O-Z) DIMENSION RT(6,6),RTRANS(6,6) DO 10 I = 1,6 DO 10 N = I,6 RT(I,N)=0. 10 RT(N,I)=0.0 RT(1,1)=CX RT(2,2)=CX RT(4,4)=CX RT(5,5)=CX RT(1,2)=CY RT(4,5)=CY RT(2,1)=-CY RT(5,4)=-CY RT(3,3)=1. RT(6,6)=1. C FORM RTRANS DO 30 I = 1,6 DO 30 N = 1,6 30 RTRANS(N,I)=RT(I,N) RETURN END SUBROUTINE MULT(A,B,JB,C) IMPLICIT REAL*8 (A-H,O-Z) C MATRIX MULTIPLY DIMENSION A(6,6),B(6,JB),C(6,JB) DO 10 I =1,6 DO 10 N =1,JB SUM=0.0 DO 20 L=1,6 20 SUM=SUM+A(I,L)*B(L,N) 10 C(I,N)=SUM RETURN END SUBROUTINE INDX(JJ,KK,IN) DIMENSION IN(6) C INDX GIVES GLOBAL DEGREE OF FREEDOM NUMBERS CORRESPONDING C TO THE SIX MEMBER LOCAL DEGREES OF FREEDOM J3=3*JJ K3=3*KK IN(3)=J3 IN(2)=J3-1 IN(1)=J3-2 IN(6)=K3 IN(5)=K3-1 IN(4)=K3-2 RETURN END SUBROUTINE FM(WL,WR,AL,FEM,NP) C SUBROUTINE TO CALC FIXED END ACTIONS FOR LINEAR LOAD IMPLICIT REAL*8(A-H, O-Z) DIMENSION FEM(6) FEM(1)=0. FEM(4)=0. IF(NP.EQ.1)GO TO 100 FEM(3)=1./12.*WL*AL**2+1./30.*(WR-WL)*AL**2 FEM(6)=-1./12.*WL*AL**2-1./20.*(WR-WL)*AL**2 GO TO 110 100 CONTINUE FEM(3)=0. FEM(6)=0. 110 CONTINUE FEM(2)=WL*AL*.5+(WR-WL)*AL*.5/3.+(FEM(3)+FEM(6))/AL FEM(5)=WL*AL*.5+(WR-WL)*AL*.333333-(FEM(3)+FEM(6))/AL RETURN END SUBROUTINE FM2(PV,PD,AL,FEM,NP) C SUBROUTINE TO CALC FEM DUE TO CONC LOAD IMPLICIT REAL*8(A-H,O-Z) DIMENSION FEM(6) PDP=AL-PD IF(NP.EQ.1)GO TO 100 FEM3=PV*PD*PDP**2/AL**2 FEM6=-PV*PD**2*PDP/AL**2 GO TO 110 100 CONTINUE FEM3=0. FEM6=0. 110 CONTINUE FEM(2)=FEM(2)+PV*PDP/AL+(FEM3+FEM6)/AL FEM(5)=FEM(5)+PV*PD/AL-(FEM3+FEM6)/AL FEM(3)=FEM(3)+FEM3 FEM(6)=FEM(6)+FEM6 RETURN END SUBROUTINE FMN(WL,WR,AL,FEM,NP,PS,EI) IMPLICIT REAL*8 (A-H,O-Z) C SUBROUTINE TO CALC FIXED END ACTIONS FOR LINEAR LOAD C INCLUDING P-DELTA EFFECTS DIMENSION FEM(6) FEM(1)=0. FEM(4)=0. IF(NP.EQ.1)GO TO 100 P=-PS RK=SQRT(P/EI) RKL=RK*AL RL2=AL*AL CO=COS(RKL) SI=SIN(RKL) DENOM=(SI-RKL)*(RK*SI) + RK*(1.-CO)**2 Q=WR-WL A=1./DENOM*(-Q*RL2/6./P*RK*SI +Q*AL/2./P*(1.-CO)) B=1./DENOM*((SI-RKL)*Q*AL/2./P + RK*(1.-CO)*Q*RL2/6./P) U=RK*AL*.5 F1=3.*(TAN(U)-U)/(U**3) C END MOMENTS DUE TO UNIFORM WL FEM(3)=1./12.*WL*(AL**2)*F1*U/TAN(U) FEM(6)=-FEM(3) C END MOMENTS DUE TO TRIANGULAR LOAD 0 AT LT, WR-WL AT RT C Q0=WL-WR C IF(ABS(Q0).LT.0.001)GO TO 80 C CON=Q0/(P*AL*SI)*(1./RK**2*SI-AL/RK*CO) C TH0A=CON- Q0*AL/3./P C TH0B=-Q0*CO/P/RK/SI + Q0/P/RK/SI - CON - Q0*AL/6./P C F2=1.5/U*(.5/U - 1./TAN(2.*U)) C F3 = 3./U*(1./SIN(2.*U) - .5/U) C AA=AL*F2/3./EI C BB=AL*F3/6./EI C CC=BB C DD=AA C DENOM=AA*DD-BB*CC C RML=(-TH0A*DD+TH0B*BB)/DENOM C RMR=(-TH0B*AA+TH0A*CC)/DENOM RML=B*P RMR=-A*P*SI-B*P*CO+Q/P*EI FEM(3)=FEM(3)-RML FEM(6)=FEM(6)-RMR 80 CONTINUE GO TO 110 100 CONTINUE FEM(3)=0. FEM(6)=0. 110 CONTINUE FEM(2)=WL*AL*.5+(WR-WL)*AL*.5/3.+(FEM(3)+FEM(6))/AL FEM(5)=WL*AL*.5+(WR-WL)*AL*.333333-(FEM(3)+FEM(6))/AL RETURN END SUBROUTINE FM2N(PV,PD,AL,FEM,NP,PS,EI) IMPLICIT REAL*8 (A-H,O-Z) C SUBROUTINE TO CALC FEM DUE TO CONC LOAD DIMENSION FEM(6),A(4,4),B(4) IF(NP.EQ.1)GO TO 100 Q=PV P=-PS AA=PD RL=AL B(1)=0. B(2)=Q/P B(3)=Q/P*(AA-RL) B(4)=B(2) RK=SQRT(P/EI) AK=AA*RK RKL=RL*RK SI=SIN(RKL) CO=COS(RKL) SIA=SIN(AK) COA=COS(AK) A(1,1)=-COA/P A(1,2)=SIA/P/RK A(1,3)=-COA A(1,4)=-SIA A(2,1)=RK*SIA/P A(2,2)=COA/P A(2,3)=RK*SIA A(2,4)= - RK*COA A(3,1)= 1./P A(3,2)= -RL/P A(3,3)=CO A(3,4)=SI A(4,1)=0. A(4,2)=1./P A(4,3)=RK*SI A(4,4)=-RK*CO CALL LDU2(A,B,1,4,4) RML=B(1) VL=B(2) C SUM MOMENTS ABOUT RIGHT END RMR=VL*RL - RML - Q*(RL - AA) C C PREVIOUSLY USED PROCEDURE COMMENTED OUT - BREAKS DOWN IF C AXIAL LOAD GREATER THAN SIMPLE BEAM BUCKLING LOAD PDP=AL-PD C IF(NP.EQ.1)GO TO 100 C P=-PS C RK=SQRT(P/EI) C RL=RK*AL C CO=COS(RL) C SI=SIN(RL) C U=RK*AL*.5 C END MOMENTS DUE TO CONC LD C TH0A=PV*SIN(RK*PDP)/P/SI - PV*PDP/P/AL C TH0B= PV*SIN(RK*PD)/P/SI - PV*PD/P/AL C F2=1.5/U*(.5/U - 1./TAN(2.*U)) C F3 = 3./U*(1./SIN(2.*U) - .5/U) C AA=AL*F2/3./EI C BB=AL*F3/6./EI C CC=BB C DD=AA C DENOM=AA*DD-BB*CC C RML=(-TH0A*DD+TH0B*BB)/DENOM C RMR=(-TH0B*AA+TH0A*CC)/DENOM C FEM3=-RML FEM3=RML FEM6=RMR GO TO 110 100 CONTINUE FEM3=0. FEM6=0. 110 CONTINUE FEM(2)=FEM(2)+PV*PDP/AL+(FEM3+FEM6)/AL FEM(5)=FEM(5)+PV*PD/AL-(FEM3+FEM6)/AL FEM(3)=FEM(3)+FEM3 FEM(6)=FEM(6)+FEM6 RETURN END C SUBROUTINE GSM IS TO CALCULATE P-DELTA EFFECTS SUBROUTINE GSM(SM,AL,NP,AF) IMPLICIT REAL*8(A-H,O-Z) DIMENSION SM(6,6) IF(NP.EQ.0)GO TO 50 C NP = 1 PIN ENDED MEMBER SM(2,2) = AF/AL SM(2,5) = - SM(2,2) SM(5,5) = SM(2,2) SM(5,2) = SM(2,5) RETURN 50 CONTINUE A=1.2*AF/AL B=0.1*AF C=2.*AL/15.*AF D=AL/30.*AF SM(2,2)=SM(2,2)+A IF(SM(2,2).LE.0.)WRITE(6,100) IF(SM(2,2).LT.0.)SM(2,2)=0. SM(5,5)=SM(2,2) SM(3,3)=SM(3,3)+C IF(SM(3,3).LE.0.)WRITE(6,100) 100 FORMAT(5X,' MEMBER STIFFNESS TERM.LT.0.') IF(SM(3,3).LT.0.)SM(3,3)=0. SM(6,6)=SM(3,3) SM(2,5)=-SM(2,2) SM(2,3)=SM(2,3)+B SM(2,6)=SM(2,3) SM(3,5)=-SM(2,3) SM(5,6)=SM(3,5) SM(3,6)=SM(3,6)-D DO 20 I=1,5 J=I+1 DO 20 M=J,6 20 SM(M,I)=SM(I,M) RETURN END C SUBROUTINE LDU2 SUBROUTINE LDU2(A,X,LU,NUMEQ,IDIM) IMPLICIT REAL*8 (A-H,O-Z) C LDU2 BEST FOR UNBANDED SYMMETRIC OR NON-SYMMETRIC MATRICES DIMENSION A(IDIM,IDIM),X(IDIM) IOUT=6 NEM1=NUMEQ-1 DO 450 I=1,NEM1 C WRITE(IOUT,5000)I,A(I,I) C 5000 FORMAT(2X,I5,F15.3) JEND=NUMEQ-I DO 440 J=1,JEND JP1=J+I IF(LU.EQ.0)GO TO 435 IF(A(I,I).EQ.0.)GO TO 480 FAC=A(JP1,I)/A(I,I) A(JP1,I)=FAC DO 430 K=1,JEND KP1=K+I A(JP1,KP1)=A(JP1,KP1)-FAC*A(I,KP1) 430 CONTINUE 435 CONTINUE X(JP1)=X(JP1)-A(JP1,I)*X(I) 440 CONTINUE 450 CONTINUE X(NUMEQ)=X(NUMEQ)/A(NUMEQ,NUMEQ) DO 470 I=1,NEM1 IROW=NUMEQ-I DO 460 J=1,I JROW=IROW+J X(IROW)=X(IROW)-A(IROW,JROW)*X(JROW) 460 CONTINUE X(IROW)=X(IROW)/A(IROW,IROW) 470 CONTINUE GO TO 490 480 WRITE(IOUT,4000) 4000 FORMAT(10X,35HDENOMINATOR EQUAL ZERO-END OF JOB ) STOP 490 CONTINUE RETURN END