C_____________________________________________________________________ C | C ***** **** **** ****** ****** ** ** | C * * * * * * * * * * * * | C * * * **** **** ****** * ** * | C * * * * * * * * * * | C ***** **** **** ****** * * * * | C | C Program Name: COBEAM | C | C Purpose: For linear stress, linearized buckling, natural | C vibration, nonlinear and post-buckling analysis | C of laminated composite curved beams using the | C the Layerwise Curved Beam Theory. | C | C Coded by: Samuel Kinde Kassegne | C | C Virginia Polytechnic Institute and State Univ. | C Blacksburg, VA 24061 | C February 1992. | C_____________________________________________________________________| C C C i) parameter variables C C NUX -- Size of global stiffness matrix C NBX -- Half-band width of global stiffness matrix C NGX -- Size of element stiffness matrix C NEX -- Size of global mass matrix, used only for eigenvalue C analysis, (i.e. NEX=NUX) otherwise, NMX = 10 C C ___________________ C C C GSTIF(NUX,NBX) -- Global stiffness matrix C GMASS(NEX,NEX) -- Global mass matrix C STIF(NGX,NGX) -- Element stiffness matrix C SM(NGX,NGX) -- Element mass matrix C ELP(NGX) -- Element force vector C ELFC(NGX) -- Element unbalanced force C GF(NGX) -- Global force vector C W(NGX) -- Displacement vector C EGNVAL(NEX) -- Array of eigen values C ___________________ C C ii) parameter variables C C NEIGEN -> = 1 for eigenvalue analysis C IPROB -> = 0 for natural vibration, C = 1 for buckling C ISINE -> = 0 for uniform load. C = 1 for sinusoidal load, C IR -> = 0 for Full Newton-Raphson, C = 1 for Modified Newton-Raphson algorithm C ICODE -> = 0 for calculating stresses at the nodes, C = 1 for calculating stresses at the Gauss points C ILOD -> = 0 for external pressure, C = 1 internal pressure C IAXIAL -> = 1 for critical axial load C ILATER -> = 1 for critical lateral load C ISHEAR -> = 1 for critical shear load C NRED -> = 1 for full integration C = 2 reduced integ. of transverse normal terms C = 3 reduced integ. of transverse shear terms C = 4 uniform reduced integration of all terms C ISTRES -> = 0 for no printouts of stress C = 1 for printing stresses and force resultants C of the beam elements C C ___________________ C C C iii) data variables C C NELEM -- Number of beam elements in the mesh C NNODE -- Number of nods in the mesh C NPE -- Number of nodes/element C MAXIT -- Maximum iterations numbers for nonlinear analysis C NLOAD -- Number of load cases C PX -- Uniform axial load C PW -- Uniform pressure load C EPSLON -- Tolerance limit for nonlinear analysis C NPRT -- Code for printing, = 2 for printing stiffness. C C_______________________________________________________________________ IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMPLEX*16 ALPHA(175) COMMON/ELAST/ THETA(30),ELX(3) COMMON/CCC/ E1,E2,E3,G12,G13,G23,ANU12,ANU13,ANU23 COMMON/MSH/ NOD(20,3),S(30),DX(15),DY(15) COMMON/DEPTH/ Z(15),THICK(15),SM(50,50),THT COMMON/SEL/ STIF(75,75),ELP(75),WO(175),ELFC(175) COMMON/GLOBAL/ GSTIF(175,175),GF(175),GFT(175),W(175) COMMON/BCS/ VBDY(75),IBDY(75) COMMON/TIT/ TITLE(20),IBSF(20),VBSF(20) COMMON/SMAS1/ GMASS(175,175),EGNVAL(175) COMMON/SITEM/ SPAN,RADIUS,PX,PW,RHO COMMON/WORKSP/ RWKSP DIMENSION T(6) C REAL RWKSP(6498) C CALL IWKIN(4498) CHARACTER*14 FLINP,FLOUT DATA IN,IT / 5,6/ DATA NRMAX,NCMAX /175,175/ C-------------------------------------------------------- C.....READ IN THE FINITE ELEMENT MESH DATA C--------------------------------------------------------- WRITE(6,222) WRITE(6,325) WRITE(6,333) 333 FORMAT(//,' Enter COBEAM Input File --->') READ (*,334) FLINP 334 FORMAT(\A) WRITE(*,335) 335 FORMAT(/,' Enter COBEAM Output File -->') READ (*,336) FLOUT 336 FORMAT(\A) OPEN(UNIT=5,FILE = FLINP,STATUS='OLD') OPEN(UNIT=6,FILE = FLOUT,STATUS='UNKNOWN') WRITE(6,222) READ (IN,300) TITLE READ(IN,*) NELEM, NPE,NNODE,ITYPE,MAXIT,NLOAD,NEIGEN,IPROB READ(IN,*) PX,PW,EPSLON,NPRNT,NRED,STEP,ICODE,ILOD,ISINE DO 20 N = 1,NELEM 20 READ(IN,*) (NOD(N,I), I=1,NPE) READ(IN,*) (S(I),I=1,NNODE) READ(IN,*) NBDY,NBSF IF(NBDY.NE.0) THEN READ(IN,*) (IBDY(I), I=1, NBDY) READ(IN,*) (VBDY(I), I=1, NBDY) END IF IF(NBSF. EQ. 0) GO TO 5 READ(IN,*) (IBSF(I), I=1, NBSF) READ(IN,*) (VBSF(I), I=1, NBSF) 5 CONTINUE C--------------------------------------------------------- C.... READ IN SKIN LAMINATE DATA C--------------------------------------------------------- READ (IN,*) NPLY READ (IN,*) (THETA(I), I=1,NPLY) READ (IN,*) (THICK(I), I=1,NPLY) 123 READ (IN,*) E1, E2, E3, G12, G13, G23, ANU12, ANU13,ANU23,RHO READ (IN,*) RADIUS, SPAN C------------------------------------------------------------- C.... COMPUTE THE HALF BAND-WIDTH C-------------------------------------------------------------- 535 ALOAD = PW NIF=NPLY+1 NDF=2*NIF NEQ = NNODE*NDF NN = 2*NPE*NIF NHBW = 0 DO 105 N=1,NELEM DO 105 I=1,NPE DO 105 J=1,NPE NW = (IABS(NOD(N,I)-NOD(N,J))+1)*NDF 105 IF(NHBW.LT.NW) NHBW = NW C----------------------------------------------------------- C.... WRITE OUT FINITE ELEMENT MESH DATA C------------------------------------------------------------ WRITE(IT,325) WRITE(IT,340) TITLE WRITE(IT,325) WRITE(IT,310) WRITE(IT,345) NELEM,NNODE,NPE,NPLY,NDF,NEQ,NHBW,SPAN,RADIUS WRITE(IT,362) E1,E2,E3,G12, G13, G23, ANU12,ANU13,ANU23 WRITE(IT,348) NRED,ICODE,MAXIT,EPSLON IF(NPRNT.NE.1) GO TO 344 WRITE(IT,315) WRITE(IT,605) (I,S(I),I=1,NNODE) WRITE(IT,615) DO 80 I = 1,NELEM 80 WRITE(IT,620) I, (NOD(I,J), J=1,NPE) WRITE(IT,354) DO 356 I =1,NBDY 356 WRITE(IT,360) IBDY(I), VBDY(I) 344 IF(NBSF.NE.0) THEN WRITE(IT,357) DO 358 I=1,NBSF 358 WRITE(IT,360) IBSF(I), VBSF(I) END IF C------------------------------------------------------ C.... WRITE OUT LAMINATE DATA C------------------------------------------------------ WRITE(IT,361) WRITE(IT,371) DO 391 I = 1,NPLY 391 WRITE(IT,381) I, THETA(I), THICK(I) WRITE(IT,373) C-------------------------------------------------------- C.... PROCESS INPUT DATA C-------------------------------------------------------- C C C..... CALCULATE THE THROUGH-THICKNESS STIFFNESSES C.... (VALID FOR ALL SECTIONS AND ITERATIONS) C NPE1 = 2 CALL SECOND(T(1)) C---------------------------------------------------------------- 635 CALL STFPLY(NPLY,RADIUS,RHO) C----------------------------------------------------------------- CALL SECOND(T(2)) DO 95 IJ =1,NEQ 95 GFT(IJ) = 0.0 DLOAD = 0.0 DO 67 IK = 1,NBSF 67 DLOAD = DLOAD + VBSF(IK) C------------------------------------------------------------------- C.... ITERATIONS FOR NON LINEAR ANALYSIS C-------------------------------------------------------------------- DO 55 LL =1,NLOAD NITER = 0 IF(NBSF.EQ.0) DLOAD = PW IF(NEIGEN.NE.0) GO TO 158 IF(PW.NE.0) WRITE(6,394) PW IF(NBSF.NE.0) WRITE(6,394) DLOAD 158 NITER = NITER+1 IF(NITER.LT.MAXIT) GO TO 159 CALL SOVER GO TO 94 159 CONTINUE C------------------------------------------------------------------ C.... INITIALIZE GLOBAL MATRICES C------------------------------------------------------------------ DO 72 I=1,NEQ GF(I)= 0.0 DO 72 J=1,NHBW GMASS(I,J)= 0.0 72 GSTIF(I,J)= 0.0 C------------------------------------------------------------------ C.... DO LOOP ON NUMBER OF ELEMENTS BEGINS C------------------------------------------------------------------- DO 150 N=1,NELEM L = 0 DO 100 I=1,NPE NI=NOD(N,I) ELX(I) = S(NI) LI = (NI-1)*NDF DO 100 J=1,NDF LI = LI+1 L = L+1 ELFC(L) = GFT(LI) 100 CONTINUE C------------------------------------------------------------------- C C.... THE TANGENT STIFFNESS MATRIX IS CALCULATED C C-------------------------------------------------------------------- C------------------------------------------------------------------- CALL STFSTR(NRED,NPLY,NPE,N,ILOD,IPROB,ISINE,HF) C--------------------------------------------------------------------- IF(NPRNT.EQ.2.AND.NITER.LE.3) THEN DO 777 I=1,2*NDF 777 WRITE(IT,878)(STIF(I,JJ),JJ=1,8) WRITE(IT,177) WRITE(IT,*) END IF 177 FORMAT(//,72('_'),///) 878 FORMAT(8(E9.2,1X)) C------------------------------------------------------------------- C.....ASSEMBLE INTO BANDED GLOBAL MATRIX C------------------------------------------------------------------- DO 140 I=1,NPE NR = (NOD(N,I)-1)*NDF DO 140 II=1,NDF NR = NR+1 L = (I-1)*NDF +II GF(NR) = GF(NR) +ELP(L) DO 140 J=1,NPE NCL = (NOD(N,J)-1)*NDF IF(NEIGEN.EQ.1) NC = NCL DO 140 JJ=1,NDF M=(J-1)*NDF+JJ IF(NEIGEN.EQ.0) THEN NC = NCL-NR+JJ+1 IF(NC.GT.0) THEN GSTIF(NR,NC) = GSTIF(NR,NC)+STIF(L,M) END IF ELSE NC = NC+1 GSTIF(NR,NC) = GSTIF(NR,NC)+STIF(L,M) GMASS(NR,NC) = GMASS(NR,NC) +SM(L,M) END IF 140 CONTINUE 150 CONTINUE CALL SECOND(T(3)) C--------------------------------------------------------------------- C FOR NEIGEN.GE.1,CARRY OUT EIGENVALUE ANALYSIS C NPROB = 0 NATURAL VIBRATION C NPROB = 1 LINEARIZED BUCKLING ANALYSIS C---------------------------------------------------------------------- 167 IF(NEIGEN.EQ.0) GO TO 157 C CALL EGNBOU (GSTIF,GMASS,IBDY,NEQ,NEQR,NBDY,NRMAX) C CALL DGVLRG (NEQR,GSTIF,NRMAX,GMASS,NRMAX,ALPHA,EGNVAL) T11 = REAL(ALPHA(1))/EGNVAL(1) IF(IPROB.EQ.0) THEN WRITE (IT,755) WRITE(IT,840) DO 276 I = 1,NEQR FRQNOR = DSQRT(DABS(E2*EGNVAL(I)/(1-ANU23**2)/SPAN**4/RHO)) FREQ = DSQRT(DABS(E2*T11/(1-ANU23**2)/SPAN**4/RHO)) WRITE(IT,845) I,T11,FREQ,FRQNOR 276 CONTINUE ELSE WRITE (IT,750) WRITE(IT,850) DO 376 I = 1,NEQR PCR1 = DSQRT(E2*T11/(1-ANU23**2)/SPAN**2) PCR2 = DSQRT(E2*EGNVAL(I)/(1-ANU23**2)/SPAN**2) IF(EGNVAL(I).NE.0) THEN TGG = REAL(ALPHA(I))/EGNVAL(I) ELSE TGG = 1.0 END IF WRITE(IT,855) I,T11,PCR1,PCR2,TGG 376 CONTINUE END IF STOP C------------------------------------------------------------------- C.....IMPOSE THE B.C'S C-------------------------------------------------------------------- 157 IF(NBSF.EQ.0) GO TO 200 DO 170 NF =1,NBSF NB = IBSF(NF) 170 GF(NB) = GF(NB)+VBSF(NF) 200 CONTINUE IF(NBDY.EQ.0) GO TO 395 CALL BNDY(NRMAX,NCMAX,NEQ,NHBW,GSTIF,GF,NBDY,IBDY,VBDY) IF(NITER.GT.1.OR.LL.GT.1) GO TO 156 DO 148 IH=1,NBDY 148 VBDY(IH) = 0.0 156 CONTINUE C------------------------------------------------------------------- C.... SOLVE THE BANDED SYMMETRIC EQUATIONS C-------------------------------------------------------------------- IRES=0 395 CONTINUE IF(NPRNT.EQ.2.AND.NITER.LE.3) THEN DO 771 I=1,NEQ 771 WRITE(IT,878)(GSTIF(I,JJ),JJ=1,NHBW) WRITE(IT,177) WRITE(IT,*) WRITE(IT,838)(GF(JJ),JJ=1,NEQ) END IF 838 FORMAT(E10.3,1X) CALL SOLVE(NRMAX,NCMAX,NEQ,NHBW,GSTIF,GF,IRES) CALL SECOND(T(4)) DO 630 IS=1,NEQ 630 GFT(IS)=GFT(IS)+GF(IS) C------------------------------------------------------------------- C.... CHECK FOR CONVERGENCE C------------------------------------------------------------------- SUM1 = 0.0 SUM2 = 0.0 DO 99 I=1,NEQ SUM1 = SUM1+GF(I)**2 99 SUM2 = SUM2+GFT(I)**2 CUT = DSQRT(SUM1/SUM2) IF(CUT.GT.EPSLON) GO TO 158 C-------------------------------------------------------------------- C..... PRINT OUT DISPLACEMENTS C--------------------------------------------------------------------- 94 WRITE(6,505) WRITE (6,515) DO 537 I = 1,NEQ 537 IF(DABS(GFT(I)).LT.1E-08) GFT(I) = 0.0 DO 518 J=1,NNODE J1=(J-1)*NDF+1 DO 528 JJ=1,NIF J2 = J1+NIF IF(PW.EQ.0) PW1 = 1.0 IF(PW.NE.0) PW1 = PW WN =(GFT(J2)*100*E1*2.5**3)/(PW1*SPAN**4) UN = E2*GFT(J1) IF(JJ.EQ.1) THEN WRITE(IT,538) J,JJ,GFT(J1),GFT(J2) ELSE WRITE(IT,548) JJ,GFT(J1),GFT(J2) END IF 528 J1=J1+1 518 CONTINUE WRITE(6,98) NITER,CUT 98 FORMAT(/,5X,'Iteration number =',I2,/,5X,'Measure of error =', + 1X,F7.5/) 54 CONTINUE C_______________________________________________________________________ C C.... POST PROCESSING OF SOLUTION C C_______________________________________________________________________ C_______________________________________________________________________ C.... COMPUTE THE DERIVATIVES OF THE SOLUTION C_______________________________________________________________________ 230 CONTINUE IF(LL.EQ.NLOAD) THEN DO 240 N=1,NELEM,4 L=0 DO 242 I=1,NPE NI=NOD(N,I) ELX(I) = S(NI) LI=(NI-1)*NDF DO 242 J=1,NDF LI=LI+1 L=L+1 GF(L) = GFT(LI) 242 CONTINUE C_______________________________________________________________________ 240 CALL STRESS(NPE,NDF,NIF,N,ICODE,NRED) C_______________________________________________________________________ END IF IF(NLOAD.GT.1) THEN READ(IN,*) PX,PW,(VBSF(I),I=1,NBSF) END IF 55 CONTINUE CALL SECOND(T(5)) TT = T(5) - T(1) DO 500 I=1,4 500 T(I)=T(I+1)-T(I) T(5) = TT WRITE (6,1000) (T(i),i=1,5) C------------------------------------------------------------------- C.....F O R M A T S C------------------------------------------------------------------- 222 FORMAT(//, 4X,67('_'),//, + 10X,5('C'),4X,4('O'),5X,4('B'),5X,6('E'),2X,6('A'),3X, + 2('M'),4X,2('M'),/, + 9X,1('C'),8X,1('O'),4X,1('O'),4X,1('B'),3X,'B',4X, + 1('E'),7X,1('A'),4X,'A',3X, + 1('M'),1X,1('M'),2X,'M',1X,'M',/, + 9X,1('C'),8X,1('O'),4X,1('O'),4X,4('B'),5X, + 4('E'),4X,6('A'),3X,'M',2X,'MM',2X,'M',/, + 9X,1('C'),8X,1('O'),4X,1('O'),4X,1('B'),3X,'B',4X, + 'E',7X,1('A'),4X,'A',3X,'M',6X,'M',/, + 10X,5('C'),4X,4('O'),5X,4('B'),5X,6('E'),2X,1('A'),4X, + 'A',3X,'M',6X,'M',//,24X,'By Samuel Kinde Kassegne',//) 300 FORMAT(20A4) 310 FORMAT(/,5X,'Mesh information :',/,5X,19('-')) 315 FORMAT(//,5X,'________ Global Coordinates Follow ______',//, + /,5X,'Node',34X,'X',/,5X,42('-'),/) 394 FORMAT(/,5X,22('_')/,5X,'Load level =',1PE10.3,/,5X,22('_')/) 325 FORMAT(4X,67('_')) 340 FORMAT(4X,17A4) 345 FORMAT(/,5X,'Num of Elements in the mesh.......=',I3,/ + ,5X,'Num of nodes in the mesh..........=',I3,/ + ,5X,'Num of nodes per element..........=',I3,/ + ,5X,'Num of layers.....................=',I3,/ + ,5X,'Num of dof/node...................=',I3,/ + ,5X,'Num of equations .............=',I3,/ + ,5X,'Half-band width .............=',I3,/ +,5X,'Length of the beam ...............=',1PE10.2,/ +,5X,'Radius of the beam................=',1PE10.2,/) 348 FORMAT(/,5X,'Code for order of integration.....=',I3,/ + ,5X,'(2-Reduced;1-Mixed;0-Full) ',// + ,5X,'Code for locus of calc. stresses =',I3,/ + ,5X,'(1-Gauss points;0-Corners) ',// + ,5X,'Maximum number of iterations .....=',I3,/ + ,5X,'Error tolerance ................=',1X,F7.5,/) 354 FORMAT(/,5X,'Primary dof Value',/, + 5X,48('-'),/) 357 FORMAT(/,5X,'Secondary dof Value',/, + 5X,48('-'),/) 360 FORMAT(5X,I3,31X,E12.5) 361 FORMAT(/,5X,'Beam Laminate Information :',/,5X,28('-'),/) 362 FORMAT(/,5X,'Material Property Information :',/,5X,31('-'),//, + 5X,'E1............... = ',1PE10.3,/, + 5X,'E2............... = ',1PE10.3,/, + 5X,'E3............... = ',1PE10.3,/, + 5X,'G12.............. = ',1PE10.3,/, + 5X,'G13.............. = ',1PE10.3,/, + 5X,'G23.............. = ',1PE10.3,/, + 5X,'Nu12............. = ',1PE10.3,/, + 5X,'Nu13............. = ',1PE10.3,/, + 5X,'Nu23............. = ',1PE10.3,//) 371 FORMAT(5X,'Ply number Ply angle Ply thickness', + /,5X,49('_'),/) 373 FORMAT(5X,49('_'),/) 381 FORMAT(5X,I1,12X,F6.2,17X,F5.2) 380 FORMAT(4X,6E12.5) 390 FORMAT(/,5X,6E12.5) 415 FORMAT(3X,I3,3X,3(1PE10.3,3X)) 505 FORMAT(/,26X,'Interface Displacements',/,26X,23('^'),/) 515 FORMAT(5X,61('_'),/,5X,'Node',3X,'Interface',8X,'U-j',12X, + 'W-j',12X,'Wn',/,5X,61('_'),/) 538 FORMAT(3X,I3,7X,I3,8X,3(1PE10.3,6X)) 548 FORMAT(13X,I3,8X,3(1PE10.3,6X)) 570 FORMAT(/5X,I3,3X,E12.5,5X,4E12.5,/,3(28X,4E12.5,/)) 605 FORMAT(3X,I5,31X,1PE10.3) 615 FORMAT(/,5X,'Boolean (connectivity) Matrix Nod(I,J)',/) 620 FORMAT(1X,20I5) 718 FORMAT(5X,4(E10.3,2X)) 750 FORMAT(//,22X,'Linear Buckling Analysis',/,22x,24('-'),//) 755 FORMAT(//,22X,'Eigen-frequency Analysis',/,22X,24('-'),//) 840 FORMAT(/,15X,51('_'),/,15X, + 'I Frequency Normalized Frequency Egnval(I)',/, + 15x,51('-'),/) 845 FORMAT(13X,I3,4X,2(E10.3,2X),8X,E10.3) 850 FORMAT(/,5X,58('_'),/,5X, + 'I',11X,'Pcr',9x,'Nor. Pcr ',3x,'Nor. Egnv',3x,'Egnval(I)', + /,5X,58('-'),/) 855 FORMAT(3X,I3,8X,E10.3,4X,3(E10.3,2X)) 1000 FORMAT (///19H Time Log (Minutes) ,// + 41H Form Ply Stiffness ................... = ,E10.3/ + 41H Form System Stiffness ................ = ,E10.3/ + 41H Solve Static Load Case ............... = ,E10.3/ + 41H Compute and Print Stresses and Disp. . = ,E10.3/ + 41H Total Time Required ..................= ,E10.3/) STOP END C------------------------------------------------------------------ C C C... THIS SUBROUTINE CALCULATES THE ELEMENT DIRECT AND C TANGENT STIFFNESS MATRIX OF A LAMINATED BEAM ELEMENT C C C------------------------------------------------------------------- SUBROUTINE STFSTR(NRED,NPLY,NPE,N,ILOD,IPROB,ISINE,HF) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /ABDF/A(9,9,3,3), ABAR(9,9,3,3), DBAR(9,9,3,3),G(9,9) COMMON /FBAR/D(9,9,9,3,3), DB(9,9,9,3,3), F(9,9,9,9,3,3) COMMON /MSH/NOD(20,3),S(30),DX(15),DY(15) COMMON /SEL/ST(75,75),ELP(75),WO(175),ELFC(175) COMMON /DOUBLE/ DVDS(20),DWDS(20),TANK(50,50) COMMON /DEPTH/ Z(15),THICK(15),SM(75,75),THT COMMON /SHP/ SF(4),GSF(4),GDDSF(4),DET COMMON /SITEM/ SPAN,RADIUS,PX,PW,RHO COMMON /ELAST/ THETA(30),ELX(3) DIMENSION WT(4,4),GAUSS(4,4) DATA GAUSS /4*0.0D0,-.57735027D0,.57735027D0,2*0.0D0, + -.77459667D0,0.0D0,.77459667D0,0.0D0,-.86113631D0, + -.33998104D0,.33998104D0,.86113631D0 / DATA WT /2.0D0,3*0.0D0,2*1.0D0,2*0.0D0,.55555555D0,.88888888D0, + .55555555D0,0.0D0,.34785485D0,2*.65214515D0,.34785485D0/ IT=6 NIF = NPLY+1 R=RADIUS NDF =2*NIF NN = NPE*NDF PI = DACOS(-1.0D0) IF(NPE.EQ.2) THEN NGP=2 IELEM = 1 NET = 2 ELSE NGP=3 IELEM = 2 NET = 3 END IF IF(NRED.EQ.1) NGP = NGP-1 KGP = 4 DO 5 I = 1,50 ELP(I) = 0.0D0 DO 5 J = 1,50 TANK(I,J) = 0.0D0 SM(I,J) = 0.0D0 5 ST(I,J) = 0.0D0 IF(N.EQ.1) THEN END IF H1 = ELX(NPE)-ELX(1) C C.... Load vector for sinusoidal loads is calculated here C IF(ISINE.EQ.1) THEN DO 50 NI = 1,KGP XI = GAUSS(NI,KGP) CALL SHAPE(XI,H1,NPE,NET,IELEM) CNST = DET*WT(NI,KGP) XX = 0.0 DO 52 I=1,NPE 52 XX = XX + SF(I)*ELX(I) DO 62 I=1,NPE I2 = (I-1)*NDF+2*NIF PWT = PW*DSIN(PI*XX/SPAN)*SF(I)*CNST ELP(I2) = ELP(I2) + PWT 62 CONTINUE 50 CONTINUE END IF C------------------------------------------------------------------- C C......FULL INTEGRATION OF TANGENT STIFFNESS MATRIX FOLLOWS C C--------------------------------------------------------------------- DO 100 NI = 1,NGP XI = GAUSS(NI,NGP) CALL SHAPE(XI,H1,NPE,NET,IELEM) CNST = DET*WT(NI,NGP) C C.... LOAD VECTOR FOLLOWS C IF(ISINE.EQ.0) THEN DO 110 I=1,NPE L1 = (I-1)*NDF +1 IF(ILOD.EQ.2) GO TO 115 IF(ILOD.EQ.0) THEN I1 = L1+NIF-1 I2 = L1 + 2*NIF-1 ELSE I1 = L1 I2 = L1 + NIF END IF ELP(I1) = ELP(I1)+PX*SF(I)*CNST ELP(I2) = ELP(I2) + PW*SF(I)*CNST GO TO 110 115 DO 118 I1=L1,L1+NIF-1 I2=I1+NIF ELP(I1) = ELP(I1)+PX*SF(I)*CNST 118 ELP(I2) = ELP(I2) + PW*SF(I)*CNST 110 CONTINUE END IF DO 15 I = 1,20 DVDS(I) = 0.0 15 DWDS(I) = 0.0 DO 130 I=1,NPE L1 =(I-1)*NDF+1+NIF L2 =(I-1)*NDF+1 DO 140 J= 1,NIF DVDS(J) = DVDS(J)+GSF(I)*ELFC(L2) DWDS(J) = DWDS(J)+GSF(I)*ELFC(L1) L1 = L1+1 L2 = L2+1 140 CONTINUE 130 CONTINUE II=1 DO 70 I=1,NPE JJ=1 DO 60 J=1,NPE I1=II JN=JJ+NPLY NE = 0 DO 55 IN=JJ,JN NE = NE+1 IF(IPROB.EQ.0) THEN CNST1 = RHO*SF(I)*SF(J)*CNST CNST2 = CNST1 ELSE CNST2 = GSF(I)*GSF(J)*CNST/THT CNST1 = 0.0 END IF CON11 = 0.0 CON22 = 0.0 CON22I = 0.0 CON33 = 0.0 CON311 = 0.0 CON312 = 0.0 CON3 = 0.0 CON4 = 0.0 CON44 = 0.0 CON55 = 0.0 CON66 = 0.0 CON555 = 0.0 CON666 = 0.0 C C------------------------------------------------------------------- C C..... CALCULATE K11(TANGENT AND DIRECT) C C------------------------------------------------------------------- C ST(I1,IN) = ST(I1,IN) +(A(NE,NE,1,1)*GSF(I)*GSF(J) + + DBAR(NE,NE,2,2)*SF(I)*SF(J) + - ABAR(NE,NE,2,2)*SF(I)*SF(J)/R + - ABAR(NE,NE,2,2)*SF(I)*SF(J)/R + + A(NE,NE,2,2)*SF(I)*SF(J)/R/R)*CNST TANK(I1,IN) = ST(I1,IN) SM(I1,IN) = SM(I1,IN) +G(NE,NE)*CNST1 IF(IN.LT.JN) THEN ST(I1,IN+1) = ST(I1,IN+1) +(A(NE,NE+1,1,1)*GSF(I)*GSF(J) + +DBAR(NE,NE+1,2,2)*SF(I)*SF(J) + -ABAR(NE,NE+1,2,2)*SF(I)*SF(J)/R + -ABAR(NE+1,NE,2,2)*SF(I)*SF(J)/R + +A(NE,NE+1,2,2)*SF(I)*SF(J)/R/R)*CNST ST(I1+1,IN) = ST(I1+1,IN) +(A(NE+1,NE,1,1)*GSF(I)*GSF(J) + +DBAR(NE+1,NE,2,2)*SF(I)*SF(J) + -ABAR(NE,NE+1,2,2)*SF(I)*SF(J)/R + -ABAR(NE+1,NE,2,2)*SF(I)*SF(J)/R + +A(NE+1,NE,2,2)*SF(I)*SF(J)/R/R)*CNST TANK(I1,IN+1) = ST(I1,IN+1) TANK(I1+1,IN) = ST(I1+1,IN) SM(I1,IN+1) = SM(I1,IN+1) +G(NE,NE+1)*CNST1 SM(I1+1,IN) = SM(I1,IN+1) END IF C C------------------------------------------------------------------- C C..... CALCULATE K12 (TANGENT AND DIRECT) C C------------------------------------------------------------------- C C C ....LINEAR TERMS FOLLOW C IP = IN+NIF CON1 = A(NE,NE,1,1)*GSF(I)*SF(J)/R + + ABAR(NE,NE,2,2)*SF(I)*GSF(J) - + A(NE,NE,2,2)*SF(I)*GSF(J)/R + + ABAR(NE,NE,1,3)*SF(J)*GSF(I) IF(IN.LT.JN) THEN CON2 = A(NE,NE+1,1,1)*GSF(I)*SF(J)/R + + ABAR(NE+1,NE,2,2)*SF(I)*GSF(J) - + A(NE,NE+1,2,2)*SF(I)*GSF(J)/R + + ABAR(NE,NE+1,1,3)*SF(J)*GSF(I) CON2I= A(NE+1,NE,1,1)*GSF(I)*SF(J)/R + + ABAR(NE,NE+1,2,2)*SF(I)*GSF(J) - + A(NE+1,NE,2,2)*SF(I)*GSF(J)/R + + ABAR(NE+1,NE,1,3)*SF(J)*GSF(I) END IF C C.... NON-LINEAR TERMS OF K12 FOLLOW C ME = NE-1 DO 150 IC=1,3 IF(ME.EQ.0) ME = NE CON11 = CON11+D(NE,NE,ME,1,1)*GSF(I)*GSF(J)*(DWDS(NE)+DWDS(ME))/4 144 ME = ME+1 150 CONTINUE IF(IN.LT.JN) THEN IZ = NE ME = NE+1 DO 160 IC=1,2 CON22=CON22+D(NE,NE+1,IZ,1,1)*GSF(I)*GSF(J)*(DWDS(ME)+DWDS(IZ))/4 IZ = IZ+1 160 CONTINUE END IF ST(I1,IP) = ST(I1,IP)+(CON1+0*CON11)*CNST TANK(I1,IP) = TANK(I1,IP)+(CON1+0*2*CON11)*CNST IF(IN.LT.JN) THEN ST(I1,IP+1) = ST(I1,IP+1)+(CON2+0*CON22)*CNST ST(I1+1,IP) = ST(I1+1,IP)+(CON2I+0*CON22)*CNST TANK(I1,IP+1)=TANK(I1,IP+1)+(CON2+0*2*CON22)*CNST TANK(I1+1,IP)=TANK(I1+1,IP)+(CON2I+0*2*CON22)*CNST END IF C------------------------------------------------------------------- C C CALCULATE K21 (TANGENT AND DIRECT) C C------------------------------------------------------------------- C IP = IN IQ = I1+NIF CON1 = A(NE,NE,1,1)*GSF(J)*SF(I)/R + + ABAR(NE,NE,2,2)*SF(J)*GSF(I) - + A(NE,NE,2,2)*SF(J)*GSF(I)/R + + ABAR(NE,NE,1,3)*SF(I)*GSF(J) IF(IN.LT.JN) THEN CON2 = A(NE,NE+1,1,1)*GSF(J)*SF(I)/R + + ABAR(NE,NE+1,2,2)*SF(J)*GSF(I) - + A(NE,NE+1,2,2)*SF(J)*GSF(I)/R + + ABAR(NE+1,NE,1,3)*SF(I)*GSF(J) CON2I= A(NE+1,NE,1,1)*GSF(J)*SF(I)/R + + ABAR(NE+1,NE,2,2)*SF(J)*GSF(I) - + A(NE,NE+1,2,2)*SF(J)*GSF(I)/R + + ABAR(NE,NE+1,1,3)*SF(I)*GSF(J) END IF C C.... NON-LINEAR TERMS OF K21 FOLLOW C ME = NE-1 DO 195 IC=1,3 IF(ME.EQ.0) ME = NE CON33 =CON33+D(NE,NE,ME,1,1)*GSF(J)*GSF(I)*(DWDS(NE)+DWDS(ME))/4 194 ME = ME+1 195 CONTINUE IF(IN.LT.JN) THEN IZ = NE ME = NE+1 DO 196 IC=1,2 CON44=CON44+D(NE,NE+1,IZ,1,1)*GSF(J)*GSF(I)*(DWDS(ME)+DWDS(NE))/4 IZ = IZ+1 196 CONTINUE END IF ST(IQ,IP) = ST(IQ,IP)+(CON1+0*2*CON33)*CNST TANK(IQ,IP) = TANK(IQ,IP)+(CON1+0*2*CON33)*CNST IF(IN.LT.JN) THEN ST(IQ,IP+1)=ST(IQ,IP+1)+(CON2+0*2*CON44)*CNST ST(IQ+1,IP)=ST(IQ+1,IP)+(CON2I+0*2*CON44)*CNST TANK(IQ,IP+1)=TANK(IQ,IP+1)+(CON2+0*2*CON44)*CNST TANK(IQ+1,IP)=TANK(IQ+1,IP)+(CON2I+0*2*CON44)*CNST END IF C C------------------------------------------------------------------- C C CALCULATE K22 (TANGENT AND DIRECT) C C------------------------------------------------------------------- C C C ....LINEAR TERMS FOLLOW C IQ = I1+NIF IP = IN+NIF CON5 = A(NE,NE,1,1)*SF(I)*SF(J)/R/R + A(NE,NE,2,2)*GSF(I)*GSF(J) + +ABAR(NE,NE,1,3)*SF(I)*SF(J)/R + +DBAR(NE,NE,3,3)*SF(I)*SF(J) IF(IN.LT.JN) THEN CON6 = A(NE,NE+1,1,1)*SF(I)*SF(J)/R/R + + A(NE,NE+1,2,2)*GSF(I)*GSF(J) + + ABAR(NE+1,NE,1,3)*SF(I)*SF(J)/R + + DBAR(NE,NE+1,3,3)*SF(I)*SF(J) CON6I= A(NE+1,NE,1,1)*SF(I)*SF(J)/R/R + + A(NE+1,NE,2,2)*GSF(I)*GSF(J)+ + ABAR(NE,NE+1,1,3)*SF(I)*SF(J)/R + + DBAR(NE+1,NE,3,3)*SF(I)*SF(J) END IF C C.... SINGLE NON-LINEAR TERMS OF K22 FOLLOW C ME = NE-1 DO 190 IC=1,3 IF(ME.EQ.0) ME = NE CON55 = CON55+D(NE,NE,ME,1,1)*GSF(J)*SF(I)*(DWDS(NE)+DWDS(ME))/4/R + +D(NE,NE,ME,1,1)*GSF(I)*SF(J)*(DWDS(NE)+DWDS(NE))/2/R + +DB(NE,ME,NE,1,3)*SF(I)*GSF(J)*(DWDS(NE)+DWDS(ME))/4 185 ME = ME+1 190 CONTINUE IF(IN.LT.JN) THEN IZ = NE DO 200 IC=1,2 CON66 = CON66+D(NE,NE+1,IZ,1,1)*GSF(J)*SF(I) + *(DWDS(ME)+DWDS(IZ))/4/R + +D(NE,NE+1,IZ,1,1)*GSF(I)*SF(J)*(DWDS(ME)+DWDS(ME))/2/R + +DB(NE+1,IZ,NE,1,3)*SF(I)*GSF(J)*(DWDS(IZ)+DWDS(ME))/4 IZ = IZ+1 200 CONTINUE END IF C C.... DOUBLE NON-LINEAR TERMS OF K22 FOLLOW C DO 220 IC=1,4 IK=NE IF(IC.GT.2) IK=NE+1 IL=NE IF(IC.EQ.2.OR.IC.EQ.4) IL=NE+1 TERM1 = 0.125*(DWDS(NE)+DWDS(NE))*GSF(I)* + F(NE,NE,IK,IL,1,1)*GSF(J)*(DWDS(IK)+DWDS(IL)) CON555 = CON555+TERM1 220 CONTINUE IF(IN.LT.JN) THEN DO 230 IC=1,4 IK=NE IF(IC.GT.2) IK=NE+1 IL=NE IF(IC.EQ.2.OR.IC.EQ.4) IL=NE+1 TERM2 = 0.125*(DWDS(NE)+DWDS(NE+1))*GSF(I)* + F(NE,NE+1,IK,IL,1,1)*GSF(J)*(DWDS(IK)+DWDS(IL)) CON666= CON666+TERM2 230 CONTINUE END IF C C.... NON-LINEAR CONTRIBUTIONS FROM (DK21/DW)*UJ FOLLOW C ME = NE-1 DO 550 IC=1,3 IF(ME.EQ.0) ME = NE CON311=CON311+D(NE,NE,ME,1,1)*GSF(J)*GSF(I)*(DVDS(NE)+DVDS(NE))/2 544 ME = ME+1 550 CONTINUE IF(IN.LT.JN) THEN IZ = NE ME = NE+1 DO 560 IC=1,2 CON312= CON312+D(NE,NE+1,IZ,1,1)* + GSF(J)*GSF(I)*(DVDS(ME)+DVDS(NE))/2 IZ = IZ+1 560 CONTINUE END IF ST(IQ,IP) = ST(IQ,IP) +(CON5+(CON55+CON555))*CNST TANK(IQ,IP) = TANK(IQ,IP)+(CON5+CON311 + +2*CON55+3*CON555)*CNST SM(IQ,IP) = SM(IQ,IP) +G(NE,NE)*CNST2 IF(IN.LT.JN) THEN ST(IQ,IP+1) = ST(IQ,IP+1)+(CON6+(CON66+CON666))*CNST ST(IQ+1,IP) = ST(IQ+1,IP)+(CON6I+(CON66+CON666))*CNST TANK(IQ,IP+1)=TANK(IQ,IP+1)+(CON6+CON312 + +2*CON66+3*CON666)*CNST TANK(IQ+1,IP)=TANK(IQ+1,IP)+(CON6I+CON312+2*CON66+3*CON666)*CNST SM(IQ,IP+1) = SM(IQ,IP+1) +G(NE,NE+1)*CNST2 SM(IQ+1,IP) = SM(IQ,IP+1) END IF I1= I1+1 55 CONTINUE 60 JJ= NDF*J+1 70 II=NDF*I+1 100 CONTINUE DO 101 NI = 1,NGP XI = GAUSS(NI,NGP) CALL SHAPE(XI,H1,NPE,NET,IELEM) CNST = DET*WT(NI,NGP) DO 11 I = 1,20 DVDS(I) = 0.0 11 DWDS(I) = 0.0 DO 131 I=1,NPE L1 =(I-1)*NDF+1+NIF L2 =(I-1)*NDF+1 DO 141 J= 1,NIF DVDS(J) = DVDS(J)+GSF(I)*ELFC(L2) DWDS(J) = DWDS(J)+GSF(I)*ELFC(L1) L1 = L1+1 L2 = L2+1 141 CONTINUE 131 CONTINUE II=1 DO 71 I=1,NPE JJ=1 DO 61 J=1,NPE I1=II JN=JJ+NPLY NE = 0 DO 51 IN=JJ,JN NE = NE+1 IF(IPROB.EQ.0) THEN CNST1 = RHO*SF(I)*SF(J)*CNST CNST2 = CNST1 ELSE CNST2 = GSF(I)*GSF(J)*CNST/THT CNST1 = 0.0 END IF CON11 = 0.0 CON22 = 0.0 CON22I = 0.0 CON33 = 0.0 CON311 = 0.0 CON312 = 0.0 CON3= 0.0 CON4= 0.0 CON44 = 0.0 CON55 = 0.0 CON66 = 0.0 CON555 = 0.0 CON666 = 0.0 C C------------------------------------------------------------------- C C..... CALCULATE K12 (TANGENT AND DIRECT) C C------------------------------------------------------------------- C IP = IN+NIF C C.... NON-LINEAR TERMS OF K12 FOLLOW C ME = NE-1 DO 151 IC=1,3 IF(ME.EQ.0) ME = NE CON11 = CON11+D(NE,NE,ME,1,1)*GSF(I)*GSF(J)*(DWDS(NE)+DWDS(ME))/4 ME = ME+1 151 CONTINUE IF(IN.LT.JN) THEN IZ = NE ME = NE+1 DO 161 IC=1,2 CON2 =CON22+D(NE,NE+1,IZ,1,1)*GSF(I)*GSF(J)*(DWDS(ME)+DWDS(IZ))/4 IZ = IZ+1 161 CONTINUE END IF M1 = 1 M2 = 1 M3 = 1 IF(TANK(I1,IP).LT.0) M1 = 1 IF(TANK(I1,IP+1).LT.0) M2 = 1 IF(TANK(I1+1,IP).LT.0) M3 = 1 ST(I1,IP) = ST(I1,IP)+M1*CON11*CNST TANK(I1,IP) = TANK(I1,IP)+M1*2*CON11*CNST IF(IN.LT.JN) THEN ST(I1,IP+1) = ST(I1,IP+1)+M2*CON22*CNST ST(I1+1,IP) = ST(I1+1,IP)+M3*CON22*CNST TANK(I1,IP+1)=TANK(I1,IP+1)+M2*2*CON22*CNST TANK(I1+1,IP)=TANK(I1+1,IP)+M3*2*CON22*CNST END IF C------------------------------------------------------------------- C C CALCULATE K21 (TANGENT AND DIRECT) C C------------------------------------------------------------------- C IP = IN IQ = I1+NIF C C.... NON-LINEAR TERMS OF K21 FOLLOW C ME = NE-1 DO 191 IC=1,3 IF(ME.EQ.0) ME = NE CON33 =CON33+D(NE,NE,ME,1,1)*GSF(J)*GSF(I)*(DWDS(NE)+DWDS(ME))/4 ME = ME+1 191 CONTINUE IF(IN.LT.JN) THEN IZ = NE ME = NE+1 DO 193 IC=1,2 CON44=CON44+D(NE,NE+1,IZ,1,1)*GSF(J)*GSF(I)*(DWDS(ME)+DWDS(NE))/4 IZ = IZ+1 193 CONTINUE END IF M1 = 1 M2 = 1 M3 = 1 IF(TANK(IQ,IP).LT.0) M1 = 1 IF(TANK(IQ,IP+1).LT.0) M2 = 1 IF(TANK(IQ+1,IP).LT.0) M3 = 1 ST(IQ,IP) = ST(IQ,IP)+M1*2*CON33*CNST TANK(IQ,IP) = TANK(IQ,IP)+M1*2*CON33*CNST IF(IN.LT.JN) THEN ST(IQ,IP+1)=ST(IQ,IP+1)+M2*2*CON44*CNST ST(IQ+1,IP)=ST(IQ+1,IP)+M3*2*CON44*CNST TANK(IQ,IP+1)=TANK(IQ,IP+1)+M2*2*CON44*CNST TANK(IQ+1,IP) =TANK(IQ+1,IP)+M3*2*CON44*CNST END IF I1= I1+1 51 CONTINUE 61 JJ= NDF*J+1 71 II=NDF*I+1 101 CONTINUE SUM =0.0 DO 260 I=1,NN DO 270 J=1,NN 270 SUM=SUM+ST(I,J)*ELFC(J) ELP(I)=ELP(I)-SUM SUM = 0.0 260 CONTINUE DO 280 I=1,NN DO 280 J=1,NN 280 ST(I,J)= TANK(I,J) RETURN END C-------------------------------------------------------------------- C C THIS SUBROUTINE CALCULATES THE PLY STIFNESS MATRICES C [A], [B], [C], [D], AND [G] C C--------------------------------------------------------------------- SUBROUTINE STFPLY(NPLY,RADIUS,RHO) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /ABDF/ A(9,9,3,3), ABAR(9,9,3,3), DBAR(9,9,3,3),G(9,9) COMMON /FBAR/ D(9,9,9,3,3), DB(9,9,9,3,3), F(9,9,9,9,3,3) COMMON /SS/ SQ(15,3,3) COMMON/SHP/ SF(4),GSF(4),DDSF(4),GJ COMMON /ELAST/ THETA(30),ELX(3) COMMON /DEPTH/ Z(15),THICK(15),SM(50,50),TH NIF=NPLY+1 RAD = RADIUS * RHO NGP=2 PI= DACOS(-1.0D0) TH = 0.0 DO 5 I=1,NPLY 5 TH = TH+THICK(I) DEP = TH Z(1)= -0.5*TH DO 10 IK=1,NPLY 10 Z(IK+1) = Z(IK) +THICK(IK) C C---- INITIALIZE THE STIFFNESS MATRICES ----------------------- C DO 40 I=1,6 DO 40 J=1,6 G(I,J) = 0.0 DO 40 K=1,3 DO 40 L=1,3 A(I,J,K,L) = 0.0D0 ABAR(I,J,K,L) = 0.0D0 DBAR(I,J,K,L) = 0.0D0 DO 40 M = 1,NIF+1 D(I,J,M,K,L) = 0.0D0 DB(I,J,M,K,L) = 0.0D0 DO 40 N = 1,NIF+1 F(I,J,M,N,K,L)= 0.0D0 40 CONTINUE DO 43 L=1,15 DO 43 I=1,3 DO 43 J=1,3 43 SQ(L,I,J)=0.0 CALL CIJ(NPLY) THK= THICK(1) G(1,1) = THK/3.0D0 DO 50 I =1,3 DO 50 J =1,3 A(1,1,I,J) = SQ(1,I,J)*THK/3.0D0 ABAR(1,1,I,J) =-SQ(1,I,J)/2.0D0 DBAR(1,1,I,J) = SQ(1,I,J)/THK D(1,1,1,I,J) = SQ(1,I,J)*THK/4.0D0 DB(1,1,1,I,J) =-SQ(1,I,J)/3.0D0 F(1,1,1,1,I,J)= SQ(1,I,J)*THK/5.0D0 50 CONTINUE THK = THICK(NPLY) N = NPLY M = NIF G(N,M) = THK/6.0D0 G(M,N) = THK/6.0D0 G(M,M) = THK/3.0D0 DO 70 I = 1,3 DO 70 J = 1,3 A(N,M,I,J) = SQ(N,I,J)*THK/6.0D0 A(M,N,I,J) = SQ(N,I,J)*THK/6.0D0 A(M,M,I,J) = SQ(N,I,J)*THK/3.0D0 ABAR(N,M,I,J) = SQ(N,I,J)/2.0D0 ABAR(M,N,I,J) =-SQ(N,I,J)/2.0D0 ABAR(M,M,I,J) = SQ(N,I,J)/2.0D0 DBAR(N,M,I,J) =-SQ(N,I,J)/THK DBAR(M,N,I,J) =-SQ(N,I,J)/THK DBAR(M,M,I,J) = SQ(N,I,J)/THK D(N,N,M,I,J) = SQ(N,I,J)*THK/12.0D0 D(N,M,N,I,J) = SQ(N,I,J)*THK/12.0D0 D(N,M,M,I,J) = SQ(N,I,J)*THK/12.0D0 D(M,M,N,I,J) = SQ(N,I,J)*THK/12.0D0 D(M,N,M,I,J) = SQ(N,I,J)*THK/12.0D0 D(M,N,N,I,J) = SQ(N,I,J)*THK/12.0D0 D(M,M,M,I,J) = SQ(N,I,J)*THK/4.0D0 DB(N,N,M,I,J) = SQ(N,I,J)/3.0D0 DB(N,M,N,I,J) =-SQ(N,I,J)/6.0D0 DB(N,M,M,I,J) = SQ(N,I,J)/6.0D0 DB(M,M,N,I,J) =-SQ(N,I,J)/3.0D0 DB(M,N,M,I,J) = SQ(N,I,J)/6.0D0 DB(M,N,N,I,J) =-SQ(N,I,J)/6.0D0 DB(M,M,M,I,J) = SQ(N,I,J)/3.0D0 F(N,N,N,M,I,J) = SQ(N,I,J)*THK/20.0D0 F(N,M,N,N,I,J) = SQ(N,I,J)*THK/20.0D0 F(N,N,M,N,I,J) = SQ(N,I,J)*THK/20.0D0 F(N,N,M,M,I,J) = SQ(N,I,J)*THK/30.0D0 F(N,M,N,M,I,J) = SQ(N,I,J)*THK/30.0D0 F(N,M,M,N,I,J) = SQ(N,I,J)*THK/30.0D0 F(N,M,M,M,I,J) = SQ(N,I,J)*THK/20.0D0 F(M,M,M,M,I,J) = SQ(N,I,J)*THK/5.0D0 F(M,M,M,N,I,J) = SQ(N,I,J)*THK/20.0D0 F(M,N,M,M,I,J) = SQ(N,I,J)*THK/20.0D0 F(M,M,N,M,I,J) = SQ(N,I,J)*THK/20.0D0 F(M,M,N,N,I,J) = SQ(N,I,J)*THK/30.0D0 F(M,N,M,N,I,J) = SQ(N,I,J)*THK/30.0D0 F(M,N,N,M,I,J) = SQ(N,I,J)*THK/30.0D0 F(M,N,N,N,I,J) = SQ(N,I,J)*THK/20.0D0 70 CONTINUE IY = 1 IF(NPLY.EQ.1) GO TO 120 DO 80 K = 2,NPLY IX = IY IY = IY+1 THK1 = THICK(K-1) THK2 = THICK(K) K1 = K-1 G(IY,IY) = (THK1+THK2)/3.0D0 G(IX,IY) = THK1/6.0D0 G(IY,IX) = THK1/6.0D0 DO 90 I = 1,3 DO 90 J = 1,3 A(IX,IY,I,J) = SQ(K1,I,J)*THK1/6.0D0 A(IY,IX,I,J) = A(IX,IY,I,J) A(IY,IY,I,J) = (SQ(K1,I,J)*THK1 + SQ(K,I,J)*THK2)/3.0D0 ABAR(IY,IY,I,J) = (SQ(K1,I,J) - SQ(K,I,J))/2.0D0 ABAR(IX,IY,I,J) = SQ(K1,I,J)/2.0D0 ABAR(IY,IX,I,J) =-SQ(K1,I,J)/2.0D0 DBAR(IY,IY,I,J) = SQ(K1,I,J)/THK1 + SQ(K,I,J)/THK2 DBAR(IX,IY,I,J) =-SQ(K1,I,J)/THK1 DBAR(IY,IX,I,J) =-SQ(K1,I,J)/THK1 D(IX,IX,IY,I,J) = (SQ(K1,I,J)*THK1+SQ(K,I,J)*THK2)/12.0D0 D(IX,IY,IX,I,J) = (SQ(K1,I,J)*THK1+SQ(K,I,J)*THK2)/12.0D0 D(IX,IY,IY,I,J) = (SQ(K1,I,J)*THK1+SQ(K,I,J)*THK2)/12.0D0 D(IY,IY,IX,I,J) = (SQ(K1,I,J)*THK1+SQ(K,I,J)*THK2)/12.0D0 D(IY,IX,IY,I,J) = (SQ(K1,I,J)*THK1+SQ(K,I,J)*THK2)/12.0D0 D(IY,IX,IX,I,J) = (SQ(K1,I,J)*THK1+SQ(K,I,J)*THK2)/12.0D0 D(IY,IY,IY,I,J) = (SQ(K1,I,J)*THK1+SQ(K,I,J)*THK2)/4.0D0 DB(IX,IX,IY,I,J) = SQ(K1,I,J)/3.0D0 DB(IX,IY,IX,I,J) = -SQ(K1,I,J)/6.0D0 DB(IX,IY,IY,I,J) = SQ(K1,I,J)/6.0D0 DB(IY,IY,IX,I,J) = -SQ(K1,I,J)/3.0D0 DB(IY,IX,IY,I,J) = SQ(K1,I,J)/6.0D0 DB(IY,IX,IX,I,J) = -SQ(K1,I,J)/6.0D0 DB(IY,IY,IY,I,J) = (SQ(K1,I,J)- SQ(K,I,J))/3.0D0 F(IY,IY,IY,IY,I,J) = (SQ(K1,I,J)*THK1+SQ(K,I,J)*THK2)/5.0D0 F(IX,IX,IX,IY,I,J) = (SQ(K1,I,J)*THK1+SQ(K,I,J)*THK2)/20.0D0 F(IX,IY,IX,IX,I,J) = (SQ(K1,I,J)*THK1+SQ(K,I,J)*THK2)/20.0D0 F(IX,IX,IY,IX,I,J) = (SQ(K1,I,J)*THK1+SQ(K,I,J)*THK2)/20.0D0 F(IX,IX,IY,IY,I,J) = (SQ(K1,I,J)*THK1+SQ(K,I,J)*THK2)/30.0D0 F(IX,IY,IX,IY,I,J) = (SQ(K1,I,J)*THK1+SQ(K,I,J)*THK2)/30.0D0 F(IX,IY,IY,IX,I,J) = (SQ(K1,I,J)*THK1+SQ(K,I,J)*THK2)/30.0D0 F(IX,IY,IY,IY,I,J) = (SQ(K1,I,J)*THK1+SQ(K,I,J)*THK2)/20.0D0 F(IY,IY,IY,IX,I,J) = (SQ(K1,I,J)*THK1+SQ(K,I,J)*THK2)/20.0D0 F(IY,IX,IY,IY,I,J) = (SQ(K1,I,J)*THK1+SQ(K,I,J)*THK2)/20.0D0 F(IY,IY,IX,IY,I,J) = (SQ(K1,I,J)*THK1+SQ(K,I,J)*THK2)/20.0D0 F(IY,IY,IX,IX,I,J) = (SQ(K1,I,J)*THK1+SQ(K,I,J)*THK2)/30.0D0 F(IY,IX,IY,IX,I,J) = (SQ(K1,I,J)*THK1+SQ(K,I,J)*THK2)/30.0D0 F(IY,IX,IX,IY,I,J) = (SQ(K1,I,J)*THK1+SQ(K,I,J)*THK2)/30.0D0 F(IY,IX,IX,IX,I,J) = (SQ(K1,I,J)*THK1+SQ(K,I,J)*THK2)/20.0D0 90 CONTINUE 80 CONTINUE 120 CONTINUE RETURN END C---------------------------------------------------------------------- C C THIS SUBROUTINE CALCULATES THE PLY STIFFNESS COEFFICIENTS C IN THE GLOBAL DIRECTION C----------------------------------------------------------------------- SUBROUTINE CIJ (NPLY) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON / SS / SQ(15,3,3) COMMON /CCC/ E1,E2,E3,G12,G13,G23,ANU12,ANU13,ANU23 COMMON /ELAST/ THETA(30),ELX(3) PI = 3.141592654 ANU21 = ANU12*E2/E1 ANU31 = ANU13*E3/E1 ANU23 = ANU13 ANU32 = ANU23*E3/E2 DELTA = 1-ANU12*ANU21-ANU13*ANU31-ANU23*ANU32-ANU21*ANU32*ANU13 + -ANU12*ANU23*ANU31 C11 = E1*(1-ANU23*ANU32)/DELTA C12 = E2*(ANU12+ANU13*ANU31)/DELTA C13 = E3*(ANU13+ANU12*ANU23)/DELTA C22 = E2*(1-ANU13*ANU31)/DELTA C23 = E3*(ANU23+ANU21*ANU13)/DELTA C33 = E3*(1-ANU12*ANU21)/DELTA C44 = G23 C55 = G13 C66 = G12 C21 = C12 C31 = C13 C32 = C23 DO 30 L=1,NPLY ANGLE = PI*THETA(L)/180.0 CN = DCOS(ANGLE) SN = DSIN(ANGLE) SN2 = SN**2 SN3 = SN**3 SN4 = SN**4 CN2 = CN**2 CN3 = CN**3 CN4 = CN**4 S11 = C11*CN4 + 2*(C12+2*C66)*SN2*CN2 +C22*SN4 S12 = (C11-C22 - 4*C66)* SN2*CN2+C12*(SN4+CN4) S16 = (C11- C12-2*C66) * SN*CN3 + (C12-C22+2*C66)*SN3*CN S22 = C11*SN4 + 2*(C12+2*C66)*SN2*CN2 + C22*CN4 S26 = (C11-C12-2*C66)*SN3*CN + (C12-C22+2*C66)*SN*CN3 S66 = (C11+C22 -2*C12-2*C66)*SN2*CN2 +C66*(SN4+CN4) S21 = S12 S61 = S16 S62 = S26 S13 = C13*CN2 + C23*SN2 S23 = C13*SN2 + C23*CN2 S33 = C33 S36 = (C13-C23)*CN*SN S44 = G23*CN2 + G13*SN2 S45 = (G13-G23) *CN*SN S55 = G13*CN2 +G23*SN2 S31 = S13 S32 = S23 S63 = S36 S54 = S45 CON1 = S16*S26-S12*S66 CON2 = S22*S66-S26*S26 CON3 = S12*S26-S16*S22 CON4 = S33*S26-S23*S66 CON5 = S26*S12-S16*S22 CON6 = S36*S26-S23*S66 CON7 = S23*S26-S36*S22 S11 = S11+S12*CON1/CON2+S16*CON3/CON2 S13 = S13+S23*CON4/CON2+S36*CON5/CON2 S33 = S33+S23*CON6/CON2+S36*CON7/CON2 S16 = S16+S12*CON4/CON2+S13*CON5/CON2 S66 = S66+S26*CON4/CON2+S36*CON5/CON2 SQ(L,1,1) = S11 SQ(L,2,2) = (S44-S45*S45/S55) SQ(L,1,3) = S13 SQ(L,3,1) = S13 SQ(L,3,3) = S33 30 CONTINUE RETURN END C C .................................................................. C SUBROUTINE USED TO IMPOSE BOUNDARY CONDITIONS ON BANDED EQUATIONS C .................................................................. C SUBROUTINE BNDY(NRMAX,NCMAX,NEQ,NHBW,S,SL,NBDY,IBDY,VBDY) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION S(NRMAX,NCMAX),SL(NRMAX) DIMENSION IBDY(NBDY),VBDY(NBDY) DO 300 NB = 1, NBDY IE = IBDY(NB) SVAL = VBDY(NB) IT=NHBW-1 I=IE-NHBW DO 100 II=1,IT I=I+1 IF (I .LT. 1) GO TO 100 J=IE-I+1 SL(I)=SL(I)-S(I,J)*SVAL S(I,J)=0.0 100 CONTINUE S(IE,1)=1.0 SL(IE)=SVAL I=IE DO 200 II=2,NHBW I=I+1 IF (I .GT. NEQ) GO TO 200 SL(I)=SL(I)-S(IE,II)*SVAL S(IE,II)=0.0 200 CONTINUE 300 CONTINUE RETURN END C ........................................................... C SOLVING A BANDED SYMMETRIC SYSTEM OF EQUATIONS C IN RESOLVING, IRES .GT. 0, LHS ELIMINATION IS SKIPPED C ........................................................... C SUBROUTINE SOLVE(NRM,NCM,NEQNS,NBW,BAND,RHS,IRES) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION BAND(NRM,NCM),RHS(NRM) MEQNS=NEQNS-1 IF (IRES .GT. 0) GO TO 90 DO 500 NPIV=1,MEQNS NPIVOT=NPIV+1 LSTSUB=NPIV+NBW-1 IF(LSTSUB.GT.NEQNS) LSTSUB=NEQNS DO 400 NROW = NPIVOT,LSTSUB C.... INVERT ROWS AND COLUMNS FOR ROW FACTOR NCOL=NROW-NPIV+1 FACTOR=BAND(NPIV,NCOL)/BAND(NPIV,1) DO 200 NCOL=NROW,LSTSUB ICOL=NCOL-NROW+1 JCOL=NCOL-NPIV+1 200 BAND(NROW,ICOL)=BAND(NROW,ICOL)-FACTOR*BAND(NPIV,JCOL) 400 RHS(NROW)=RHS(NROW)-FACTOR*RHS(NPIV) 500 CONTINUE GO TO 101 90 DO 100 NPIV=1,MEQNS NPIVOT=NPIV+1 LSTSUB=NPIV+NBW-1 IF(LSTSUB.GT.NEQNS) LSTSUB=NEQNS DO 110 NROW=NPIVOT,LSTSUB NCOL=NROW-NPIV+1 FACTOR=BAND(NPIV,NCOL)/BAND(NPIV,1) 110 RHS(NROW)=RHS(NROW)-FACTOR*RHS(NPIV) 100 CONTINUE C.... BACK SUBSTITUTION 101 DO 800 IJK=2,NEQNS NPIV=NEQNS-IJK+2 RHS(NPIV)=RHS(NPIV)/BAND(NPIV,1) LSTSUB=NPIV-NBW+1 IF(LSTSUB.LT.1) LSTSUB=1 NPIVOT=NPIV-1 DO 700 JKI=LSTSUB,NPIVOT NROW=NPIVOT-JKI+LSTSUB NCOL=NPIV-NROW+1 FACTOR=BAND(NROW,NCOL) 700 RHS(NROW)=RHS(NROW)-FACTOR*RHS(NPIV) 800 CONTINUE RHS(1)=RHS(1)/BAND(1,1) RETURN END SUBROUTINE SHAPE (XI,H,NPE,NET,IELEM) C C ................................................................. C EVALUATION OF THE INTERPOLATION FUNCTIONS AND THEIR GLOBAL C DERIVATIVES AT THE GAUSS POINTS FOR LINEAR, QUADRATIC AND CUBIC C LAGRANGE ELEMENTS C .................................................................. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/SHP/SF(4),GDSF(4),GDDSF(4),GJ DIMENSION DSF(4),DDSF(4) IF(IELEM-2)30,40,40 C C INTERPOLATION FUNCTIONS FOR TWO, THREE, AND FOUR NODE ELEMENTS C C LINEAR (LAGRANGE) INTERPOLATION FUNCTIONS C 30 SF(1) = 0.5*(1.0-XI) SF(2) = 0.5*(1.0+XI) DSF(1) = -0.5 DSF(2) = 0.5 GO TO 65 C C... QUADRATIC (LAGRANGE) INTERPOLATION FUNCTIONS C 40 SF(1) = -0.5*XI*(1.0-XI) SF(2) = 1.0-XI*XI SF(3) = 0.5*XI*(1.0+XI) DSF(1) = -0.5*(1.0-2.0*XI) DSF(2) = -2.0*XI DSF(3) = 0.5*(1.0+2.0*XI) GOTO 65 C C.... COMPUTE THE DERIVATIVES OF SF(I) WITH RESPECT TO X C 65 DO 66 I=1,NPE 66 DDSF(I)=0.0 60 GJ = H*0.5 DO 80 I = 1,NPE GDSF(I) = DSF(I)/GJ 80 GDDSF(I)=DDSF(I)/GJ/GJ RETURN END C-------------------------------------------------------------------- C STOPS THE PROGRAM WHEN THE ITERATION LIMIT IS EXCEEDED C-------------------------------------------------------------------- SUBROUTINE SOVER WRITE(6,10) 10 FORMAT(///,5X,'Max number of iterations exceeded!') RETURN END C C---------------------------------------------------------------------- C C.....THIS SUBROUTINE IMPOSES BOUNDARY CONDITION ON THE FULLY POPULATED C STIFFNESS AND MASS MATRICES C----------------------------------------------------------------------- SUBROUTINE EGNBOU (A,D,IBDY,NEQ,NEQW,NBDY,NRMAX) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(NRMAX,NRMAX),D(NRMAX,NRMAX),IBDY(125) DO 20 I=1,NBDY IMAX=IBDY(I) DO 10 J=I,NBDY IF (IBDY(J).LT.IMAX) GO TO 10 IMAX=IBDY(J) IKEPT=J 10 CONTINUE IBDY(IKEPT)=IBDY(I) IBDY(I)=IMAX 20 CONTINUE NEQW = NEQ DO 80 I=1,NBDY IB=IBDY(I) IF (IB .GE. NEQW)GO TO 70 NEQW1=NEQW-1 DO 60 II=IB,NEQW1 DO 40 JJ=1,NEQW D(II,JJ)=D(II+1,JJ) 40 A(II,JJ)=A(II+1,JJ) DO 50 JJ=1,NEQW D(JJ,II)=D(JJ,II+1) 50 A(JJ,II)=A(JJ,II+1) 60 CONTINUE 70 NEQW=NEQW-1 80 CONTINUE RETURN END C---------------------------------------------------------------------- C C----- THIS SUBROUTINE CALCULATES THE NORMAL AND TRANSVERSE C----- STRESSES THROUGH THE THICKNESS OF THE LAMINA. STRESSES ARE C----- CALCULATED AT THE GAUSS POINTS DEPTH WISE AND SURFACE WISE. C C----------------------------------------------------------------------- SUBROUTINE STRESS(NPE,NDF,NIF,N,ICODE,NRED) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NMX=10,NEX=25,NUX=250,NGX=150,NBX=150) COMMON /DOUBLE/ DVDS(20),DWDS(20),TANK(50,50) COMMON /SINGLE/ U(15),V(15),WS(15) COMMON /DWDY/ DVDX(20),DVDY(20) COMMON / SS/ SQ(15,3,3) COMMON /SITEM/ SPAN,RADIUS,PX,PW,RHO COMMON /SHP/ SF(4),GSF(4),GDDSF(4),GJ COMMON /MSH/ NOD(20,3),S(30),DX(15),DY(15) COMMON /DEPTH/ Z(15),THICK(15),SM(50,50),THT COMMON /ELAST/ THETA(30),ELX(3) COMMON /SEL/ STIF(75,75),ELP(75),WO(175),ELFC(175) COMMON /GLOBAL/ GSTIF(175,175),W(175),GFT(175),WG(175) DIMENSION GAUSS(4,4),WT(4,4),SPONT(4,4) DATA GAUSS/ 4*0.0D0,-.57735027D0,.57735027D0,2*0.0D0, + -.77459667D0,0.0D0,.77459667D0,0.0D0,-.86113631D0, + -.33998104D0,.33998104D0,.86113631D0/ DATA SPONT /4*0.0D0,-1.0D0,1.0D0,2*0.0D0,-1.0D0,0.0D0,1.0D0,0.0D0, + -1.0D0,-.5D0,.5D0,1.0D0/ DATA WT /2.0D0,3*0.0D0,2*1.0D0,2*0.0D0,.55555555D0,.88888888D0, + .55555555D0,0.0D0,.34785485D0,2*.65214515D0,.34785485D0/ NPLY = NIF-1 R = RADIUS NET = 2 IF(NPE.EQ.2) THEN NGP = 2 ELSE NGP = 3 END IF IF(NRED.EQ.1) NGP = NGP-1 WRITE(6,90) N DO 10 NI = 1,NGP IF(ICODE.EQ.1) THEN XI = GAUSS(NI,NGP) ELSE XI = SPONT(NI,NGP) END IF IELEM = 1 H1 = ELX(NPE)-ELX(1) CALL SHAPE(XI,H1,NPE,NET,IELEM) X = 0.0 DO 20 J = 1,NIF V(J) = 0.0 WS(J) = 0.0 DVDS(J) = 0.0 20 DWDS(J) = 0.0 DO 30 I = 1,NPE L1 = (I-1)*NDF+1 X = ELX(1) + 0.5*H1*(1+XI) DO 40 J = 1,NIF L2 = L1+NIF V(J) = V(J)+SF(I)*W(L1) WS(J)= WS(J)+SF(I)*W(L2) DVDS(J) = DVDS(J) + GSF(I)*W(L1) DWDS(J) = DWDS(J) + GSF(I)*W(L2) L1 = L1+1 40 CONTINUE 30 CONTINUE WRITE(6,120) NGP1=2 J2 = 1 DO 50 J=1,NPLY J1 = 2*(J-1) DO 60 N1 = 1,NGP1 J1 = J1+1 H = Z(J+1)-Z(J) IF(ICODE.EQ.1) THEN X1 = GAUSS(N1,NGP1) ELSE X1 = SPONT(N1,NGP1) END IF NPE1= 2 CALL SHAPE (X1,H,NPE1,NET,1) IF(ICODE.NE.1) THEN IF(N1.EQ.1) THEN SF(1) = 1.0 SF(2) = 0.0 ELSE SF(2) = 1.0 SF(1) = 0.0 END IF END IF NL = J SIGXX = SQ(NL,1,1)*((DVDS(J)+WS(J)/R)*SF(1)+(DVDS(J+1)+WS(J+1)/R) + *SF(2)+0.5*(DVDS(J)*SF(1))**2+0.5*DVDS(J)*DVDS(J+1)*SF(1) + *SF(2)) + +SQ(NL,1,3)*( WS(J)*GSF(1)+WS(J+1)*GSF(2)) SIGZZ = SQ(NL,1,3)*((DVDS(J)+WS(J)/R)*SF(1)+(DVDS(J+1)+WS(J+1)/R) + *SF(2)+0.5*(DVDS(J)*SF(1))**2+0.5*DVDS(J)*DVDS(J+1)*SF(1) + *SF(2)) + +SQ(NL,3,3)*( WS(J)*GSF(1)+WS(J+1)*GSF(2)) SIGXZ = SQ(NL,2,2)*(((DWDS(J)-V(J)/R)*SF(1)) + V(J)*GSF(1)) + + (((DWDS(J+1)-V(J+1)/R)*SF(2))+V(J+1)*GSF(2)) IF(J1.EQ.1) WRITE(6,45) J1,X,SIGXX,SIGXZ,SIGZZ IF(J1.GT.1) WRITE(6,46) J1,SIGXX,SIGXZ,SIGZZ 60 CONTINUE 50 CONTINUE WRITE(6,68) WRITE(6,85) DO 80 J=1,NIF IF(J.EQ.1) WRITE(6,45)J,X,DVDS(J),DWDS(J) IF(J.GT.1) WRITE(6,46)J,DVDS(J),DWDS(J) 80 CONTINUE WRITE(6,68) 10 CONTINUE 45 FORMAT (3X,I3,1X,1(F8.2,3X),E10.3,2X,2(E10.3,1X)) 46 FORMAT (3X,I3,12X,E10.3,2X,2(E10.3,1X)) 68 FORMAT(5X,49('-')) 85 FORMAT(//,5X,'[INTERFACE STRAINS] :',//,5X,49('_'),/,5X,'J', + 3X,' X ',7X,'Dvds(J)',5X,'Dwds(J)',/,5X,49('_')) 90 FORMAT(/,26X,18('_'),//,26X,'Element Number=',I2,/,26X,18('_'),/) 100 FORMAT(/,5X,'x=',F6.2,2X,'y=',F6.2,/) 120 FORMAT(/,5X,'[INTERAFCE STRESSES] :',//,5X,49('_'),/,5X,'J',3X, + ' X ',7X,'Sigxx', + 7X,'Sigxz', 7X,'Sigzz',/,5X,49('_')) RETURN END C GETTIM IS A FORTRAN-callable SUBROUTINE IN MICROSOFT FORTRAN C VERSION 5.0 RETURNING THE SYSTEM TIME WITH ARGUMENTS IN INT*2 SUBROUTINE SECOND(T) REAL*8 T INTEGER*2 HOUR, MINUTE, SECND, HUNDREDTH CALL GETTIM (HOUR, MINUTE, SECND, HUNDREDTH) T = ((DBLE(HOUR)*3600.0)+(DBLE(MINUTE)*60.0)+DBLE(SECND)+ + (DBLE(HUNDREDTH)/100.0)) RETURN END